Transcriptome analysis of tertiary lymphoid structures (TLSs)-related genes reveals prognostic value and immunotherapeutic potential in cancer

Data source and preprocessing

This study utilized data from The Cancer Genome Atlas Project (TCGA), which included transcriptome expression data from 32 tumors (n=10,942) and corresponding clinical information (except for acute myeloid leukemia (LAML)). Additionally, we collected the METABRIC cohort from Bioportal, which included expression profiles and clinical information from 1,904 breast invasive carcinoma (BRCA) patients [12], as well as data from Chin et al. for 112 breast cancer patients [13]. Furthermore, we obtained data from the GEO database for breast cancers GSE20685 [14] and GSE21653 [15], as well as GSE8894 [16], GSE31210 [17], and GSE30219 [18] in the lung adenocarcinoma (LUAD) dataset, GSE19188 [19] in the lung squamous cell carcinoma (LUSC) dataset, and GSE41613 [20] data from the oral squamous cell carcinoma (OSCC) dataset. Lastly, we obtained data from the “IMvigor210” cohort [21], GSE91061 cohort [22] and GSE135222 cohort [23] for cancer immunotherapy. See Table S2.

Calculation of enrichment score

To analyze the enrichment score (ES), we performed a gene set variation analysis using the “GSVA” software package (version 1.40.1; Universitat Pompeu Fabra; Barcelona, Catalonia, Spain) [24], where the list of TLSs, EMT and cell cycle-related genes is provided in Table S1.

Construction of prognostic TLS-related signature

LASSO algorithm is a linear regression method aimed at feature selection and model sparsification. We utilized LASSO model to find relatively linearly independent sets of genes to construct our prognostic signature. We divided the TCGA dataset into training and testing sets in a 5:5 ratio, and used the LASSO algorithm to process the training set data for TLS-related genes to ensure the simplicity of the model and reduce the risk of overfitting.

Then, we used a multivariate Cox proportional hazards regression model to link the expression of each gene in the TLS-related signature and the OS of patients in the training set. As shown in the following formula:

where 35 represents the number of TLS-related signature genes, Expr (g _i ) represents the normalized expression value of gene i, and coef (g _i ) represents the corresponding coefficient estimated by the multiple Cox regression model (see Table S3). Patients were divided into high- and low-risk stratifications using the optimal cutoff inferred by the surv_cutpoint function in the “survminer” package (version 0.4.9; Datanovia; Montpellier, France) [25].

Functional enrichment analysis

The biological processes associated with differentially expressed genes (DEGs) in the high- and low-risk groups were investigated. The bitr function from “clusterProfiler” R package (version 4.0.5; Jinan University; Guangzhou, China) [26] was applied to convert gene symbols. Significant enrichment was defined as gene ontology (GO) terms with adjusted p-values less than 0.01.

Pathological sections of HCC

A total of 41 pathological sections from patients with HCC were collected in this study, and were obtained from Nanjing Drum Tower Hospital, Nanjing, China (Ethical approval number: 2021-566-02). The section of each sample contained the tumor region as well as the peritumoral region. All patients had pathologically confirmed HCC and underwent curative resection.

Immunofluorescence

Immunofluorescence (IF) staining was performed on pathological sections of 41 HCC samples. Primary antibodies were divided into two groups: CD3, CD20, CD21, and CD23 for TLS; HLA-ABC, HLA-DR+DP+DQ, HLA-E, and HLA-G for HLA (see Table S4). HLA-ABC represents the classical HLA-I molecule; HLA-DR+DP+DQ represents the HLA-II molecule; HLA-E and HLA-G represent non-classical HLA-I molecules. Staining protocols were based on the antibody manufacturer’s protocol. Slides were scanned and imaged by the Vectra3.0 System (version 3.0.5; PerkinElmer; Hopkinton, MA, USA).

Quantification of TLSs and calculation of TLS density

TLSs were detected by the combination of different marker signals on the whole slide images (WSIs), as described previously [27]. The surface areas of the different regions were measured using ImageJ software (version 1.0; National Institutes of Health; Bethesda, MD, USA) [28], and the density of TLSs was calculated as the number of TLSs per square mm² in the peritumor and intratumor regions.

Calculation of HLA positive rate

Five areas were randomly selected in peritumoral or tumoral regions of each WSI, which were imaged as multispectral HPFs (200×). To calculate the positive rate of HLA in different regions for each sample, the cell segmentation algorithm of Inform software (version 2.6.0; Akoya Biosciences; Marlborough, MA, USA) was trained to identify cells expressing HLA-ABC, HLA-D, HLA-E and HLA-G. The number of HLA positive cells and the total number of cells were calculated for each HPF image. The number of HLA-positive cells and the total number of cells were summed for the five HPF images in the peritumoral or tumoral region, and the ratio of the two represents the positivity rate of the corresponding HLA molecule in the peritumoral or tumoral region.

Statistical analysis

Pearson and Wilcoxon rank tests were utilized to measure the similarity or difference. Additionally, the survival probabilities of patients were assessed using Kaplan–Meier analysis, and the log-rank test was used to determine any differences between the groups. A nomogram was created and validated by combining multiple predictive indicators to present the results of the multivariable Cox regression analysis in a clear and understandable way. ns: not significant; * p≤0.05; ** p≤0.01; *** p≤0.001. All statistical analyses were conducted using R software (version 4.2.2; R Foundation for Statistical Computing, Vienna, Austria).

Heterogeneity and prognostic value of TLSs across cancers

Based on previous studies in TLS, 79 TLS-related genes [3], [4], [5, 29], [30], [31], [32], [33], [34], [35], [36], [37] were collected and incorporated into our study, which were composed of chemokine genes and immune cell genes. They are mainly associated with the neogenesis of TLSs and the function of immune cells in TLSs. These TLS-related genes were analyzed in a total of 10,942 samples collected from 32 tumors in the TCGA database (Figure 1A; Table S1), and the results indicated that the distribution of TLSs exhibited significant heterogeneity among different cancer types (Figure 1B). Moreover, the TLS enrichment score (ES) was calculated for the 32 tumors, revealing significant differences among different cancer types (Figure 1C; see Materials and Methods). Some tumors, such as lymphoid, diffuse large B-cell lymphoma (DLBC), LUSC, and kidney renal clear-cell carcinoma (KIRC), presented high ESs, hinting toward a frequent presence of TLS, while other tumors, such as glioblastoma multiforme (GBM) and uveal melanoma (UVM), exhibited low ESs, which was consistent with the characteristics of immune privilege (Figure 1C). Further analysis showed that the high-ES group had a longer OS than the low-ES group in pan-cancer (p≤0.0001, t-test) (Figure 1D).

Figure 1:

Heterogeneous expression pattern of TLSs in 32 cancers that are strongly associated with prognosis. (A) Circle plot shows an overview of samples from 32 cancers collected from the TCGA database, with LAML excluded. The outer circle represents the total number of samples for each cancer, while the inner circle displays the number of tumor samples. (B) The t-SNE plot shows the projection of pan-cancer samples based on TLS-related genes (Table S1). Each dot represents one sample, with cancer types marked by color codes. (C) Violin plot shows the distribution of TLS enrichment scores among pan-cancer, ordered by the median scores. (D) Boxplot shows the distribution of OS between the high- and low-ES groups, as determined by the median scores of pan-cancer. The p-value was obtained by t-test.

Developing a TLS-related signature for pan-cancer clinical prognosis

To obtain a reliable and effective TLS-related signature for predicting the prognosis of pan-cancer, we conducted LASSO regression analysis (Figure S1A and B), and only 35 genes remained effective. We first investigated the expression landscape of these 35 genes in pan-cancer and found significant differences among tumors, indicating the heterogeneity of expression patterns, which was consistent with our previous results (Figure 1B and 2A). To investigate the expression effects on prognosis, univariate Cox regression model based on training pan-cancer cohorts was conducted, demonstrating that all 35 genes were significantly associated with the prognosis (Figure 2B; p≤0.0001, log-rank test) and the per-cancer Cox analyses are shown in Figure S2. Based on these selected TLS-related genes, we constructed a prognostic signature for pan-cancer, namely TLSig. As expected, in the training cohort and testing cohort, TLSig conveniently stratified patients into high- and low-risk groups with significantly different OS (p≤0.0001, log-rank test) utilizing the optimal cutoff point, indicating a worse prognosis in the high-risk group (Figure 2C and D).

Figure 2:

Construction of a TLS-related signature for pan-cancer. (A) Heatmap shows the expression levels of 35 selected TLS genes in a wide range of tumors, where cancers without normal samples were excluded from consideration. Red indicates significantly upregulated expression in tumors compared to normal samples, while green indicates significantly downregulated expression. Gray indicates no significant difference in expression. (B) Univariate Cox regression analysis for 35 selected TLS genes based on the pan-cancer training cohort. (C and D) TCGA cohorts were divided into high and low-risk groups based on the optimal cutoff point for the TLS score. The difference in patient survival between the different risk groups in the training (C) and testing (D) sets was assessed separately, with p-values calculated using a Cox proportional hazards model. (E) Scatter plot showing the association between TLS enrichment scores (ES) and risk scores in pan-cancer. “R” indicates the Pearson correlation coefficient, and p-value was obtained by t-test. (F) Boxplot showing the distribution of enrichment scores between the high- and low-risk groups depicted by TLS risk scores. The p-value was obtained by the Wilcoxon rank test. (G) In TCGA, 32 cancer types were divided into eight groups, and then the impact of the TLSig on the prognosis of individual cancer patients was investigated in the training and testing sets of the pan-cancer cohort. TCNS, Tumors of Central Nervous System; TT, Thoracic Tumors; TDST, Tumors of Digestive System; BT, Breast Tumor; TUS, Tumors of Urinary System and Male Genital Organs; TFRO, Tumors of Female Reproductive Organs; TEO, Tumors of Endocrine Organs; TO, Tumors of others.

Next, a negative relationship was observed between the TLS risk score and the TLS enrichment score. This suggested that low-risk groups had better prognoses and higher TLS enrichment scores, which aligns with previous findings (Figures 1D, 2E and F). The impact of TLSig on prognosis of each cancer type was also examined in both the training and testing groups, and the results indicated that the TLS risk score was mostly associated with poor prognosis across various cancer types (Figures S3 and S4). Specifically, it was found to be significant for LUSC, KIRC, HNSC, liver hepatocellular carcinoma (LIHC), and BRCA (Figure 2G and Figures S3 and S4). Notably, there was no significant impact on the prognosis of GBM, testicular germ cell tumor (TGCT), or UVM, suggesting that TLSs face difficulty in generating an independent microenvironment in immune-privilege cancers.

Evaluation of TLSig’s prognostic performance in independent publicly available cancer cohorts

To assess the prognostic performance of TLSig, multiple independent cancer cohorts were analyzed, which presented similar trends. Specifically, validation results from four external breast cancer cohorts indicated that TLSig was significantly linked to patient prognosis, which was consistent with the results from the TCGA cohort. This suggests that TLSig is a reliable predictor of prognosis for breast cancer patients (Figure 3A–D; Table S2). Furthermore, external datasets for lung adenocarcinoma (LUAD), LUSC, and OSCC were also analyzed, and the results were as expected, except for the GSE31210 cohort (LUAD). In the GSE31210 cohort, no significant difference was observed, although presented with the same trend (Figure 3E–I).

Figure 3:

Validate the performance of TLSig on other independent publicly available cohorts. (A–L) KM curves show the differences in prognosis between the high and low risk groups as portrayed by the TLS risk scores in 12 independent cohorts. The p-values were evaluated by the log-rank test. (M–O) Stack bar plot shows the comparisons of the disease treatment effect between the high-risk and low-risk groups in the cohorts for immunotherapy. Red indicates complete or partial remission of the disease (CR/PR), while blue indicates disease stability or progression (SD/PD).

Given the widespread interest in the role of TLSs in immunotherapy [5], [6], [7], the predictive value of TLSig was tested in the IMvigor210 cohort, GSE91061 cohort and GSE135222 cohort. As expected, there was a potentially better prognosis in the low-risk group, as demonstrated by the KM curve (Figure 3J–L; Table S2). Additionally, the distribution of CR/PR and SD/PD in both groups revealed that a higher percentage of the low-risk group (IMvigor210 cohort: 31.5 vs. 12.5 %; GSE91061 cohort: 52.6 vs. 0 % and GSE135222 cohort: 85.7 vs. 0 %) benefited from immunotherapy (Figure 3M–O). Overall, these results suggest that TLSig has the potential for clinical application in predicting prognosis for pan-cancer and its immune response.

Construction of a prognostic nomogram model based on the TLS risk score for pan-cancer

To provide an intuitive and quantifiable approach to facilitate the clinical prediction of patient survival probabilities, we developed a comprehensive nomogram prediction model that integrates various clinicopathological characteristics, including TLS risk scores, patient age, and cancer type, as shown in Figure 4A. To validate the accuracy of the nomogram, we plotted the calibration curve for 3-year DFS (disease-free survival) and compared it with the standard curve. The results indicated that the prediction of survival probability was closely related to the actual survival (Figure 4B). To further evaluate the effectiveness of our model, we calculated the AUC of TLSig and nomogram prediction in both the training and testing cohorts. As shown in Figure 4C and D, the AUC of the nomogram prediction was higher than that of the TLS risk score alone in both cohorts (training set: 0.710 vs. 0.698, testing set: 0.730 vs. 0.699). This suggested that the predictive efficiency of the nomogram would be improved if clinicopathological factors were included.

Figure 4:

Predictive efficiency of the TLS risk score and nomogram in pan-cancer. (A and B) The nomogram shows the performance of a multivariate Cox regression analysis by combining multiple predictors, including TLS risk score, age, and cancer type. (B) The correlation between predicted and actual values is revealed by a calibration curve that demonstrates the concordance between the nomogram-predicted probability of 3-year DFS and the actual 3-year DFS. (C and D) ROC curves show the clinical benefits of the TLS risk score at 3, and 5 years in TCGA training (C) and TCGA testing (D) cohorts and the clinical benefits of the nomogram at 5 years in TCGA training (C) and TCGA testing (D) cohorts.

Association of TLSig with biological processes of cancers

To investigate the relationship between TLSig and tumor malignant characteristics, we carried out a comprehensive analysis of the correlation between TLS risk scores, epithelial to mesenchymal transition (EMT), and the cell cycle. Our analysis revealed that there were statistically significant positive correlations between the TLS risk score and EMT (R=0.21, p<0.0001) as well as the cell cycle (R=0.28, p<0.0001) in the entire TCGA pan-cancer cohort (Figure 5A and B). Furthermore, we examined individual tumors and discovered extremely significant positive correlations between TLS risk scores and EMT, or the cell cycle, in LIHC and breast cancer (BRCA) (Figure 5C). The higher the abundance of TLSs in LIHC and BRCA suggests that there is a more active immune response compared to other cancers. In summary, our findings indicate that tumors with a higher TLS risk score (i.e., a lower abundance of TLSs) exhibit more active cell proliferation and aggressive tumor cells within the microenvironment.

Figure 5:

Associations between TLS risk scores and biological processes of cancers. (A) Two plots show the correlation between the cell-cycling enrichment score and EMT enrichment score with the TLS risk score for the total sample of 32 cancer types. R represents the correlation coefficient, and p-value was determined by t-test. (B and C) Scatter plots show the correlations between the cell-cycling score (B) and EMT score (C) with the TLS risk score for each cancer type in the TCGA dataset, except for LAML. (D) GSEA shows the top three enriched GO terms in the high-risk (upper panel) and low-risk (lower panel) groups.

In addition, we identified DEGs in the high and low risk groups. The GO analysis revealed that high-risk DEGs were mainly abundant in cellular modified amino acid and phenol-containing compound metabolic processes, while low-risk DEGs were more associated with humoral immune response and defense response to bacterium (Figure 5D).

Differential distribution of TLSs in tumor and nontumor tissues of HCC

Further calculation of TLS risk scores for different types of tumor samples compared with nontumor samples revealed that there were significant differences between the different types of tumors. We also observed that tumor samples from LIHC presented higher TLS risk scores than nontumor samples, implying the presence of fewer TLSs in tumor samples (Figure 6A).

Figure 6:

The heterogeneous distribution of TLSs in tumor and nontumor tissues of HCC. (A) Comparison of the differences in TLS risk scores between tumor and nontumor samples across cancers. (B) Differences in the expression of 12 genes between tumor and nontumor samples of LIHC. (C and D) Statistical analysis of TLS density in HCC sections. Blue dots indicate the density of TLSs in the tumor tissues, and orange dots indicate the density of TLSs in the peritumoral regions. (E and F) Representative images of HCC that present TLSs in the tumoral and peritumoral regions by immunofluorescence staining. Each boxed area represents one TLS.

To validate the results of the analysis of LIHC, we first utilized the 12-Chemokine signature [38] to analyze tumor and nontumor samples of LIHC in TCGA dataset. Our results demonstrated that the expression level of 12 genes in tumor samples of LIHC was lower overall than that in nontumor samples, which was in line with the results of the previous analysis using the TLS risk score and suggested that the TLS risk score could accurately reflect the expression of TLSs (Figure 6B).

Subsequently, immunofluorescence (IF) staining (CD3, CD20, CD21, CD23 on the same sections for TLSs) was conducted on 41 pathological sections of HCC samples (Figure S5A), and the statistical results showed that the density of TLSs was higher in the peritumoral region (nontumor tissue) than in tumor tissue in 39 samples (Figure 6C–F), which was consistent with the results of the previous analysis. Together, these findings indicate that TLSs are less susceptible to form in tumor tissues compared to nontumor tissues of HCC.

HCC may affect the formation of TLSs by regulating HLA expression

To investigate the potential factors for the difference in TLS distribution, we further analyzed the differential expression of HLA. Multiplex immunostaining (HLA-ABC, HLA-D, HLA-E, HLA-G) was performed on 41 identical pathological sections of HCC samples (Figure 7A and Figure S5B). Statistical analysis revealed that tumor tissues presented higher expression of HLA-E and HLA-G compared to peritumoral regions (nontumor tissues), and HLA-ABC and HLA-D were downregulated in tumor tissues, without differences in the coexpression of HLA-E and HLA-G between them (Figure 7B–F).

Figure 7:

Differences in the expression of HLA between tumor and nontumor tissues of HCC. (A) Flow chart of HLA expression analysis. (B–F) Paired statistical analysis of HLA expression in HCC sections. Blue dots indicate the positive rate of HLA in the tumor tissues, and orange dots indicate the positive rate of HLA in the peritumoral regions. (G) Bubble plot of TLS density and HLA expression data. (H) Representative images of HCC showing the expression of different HLA and TLS densities in the tumoral and peritumoral regions.

The TLS density and HLA positive rate for each sample are summarized in Figure 7G. After comparing the TLS density with the expression of each HLA molecule in tumor tissues and peritumoral regions for each sample, we found that the trend of TLS density in paired tumor and nontumor tissues tended to be the same as the trend of HLA-ABC expression, and the opposite trend of HLA-E and HLA-G (i.e., when the TLS density in tumor tissues was lower than that in nontumor tissues, its corresponding HLA-ABC expression in tumor tissues was lower than that in nontumor tissues, and HLA-E and/or HLA-G expression were higher than that in nontumor tissues), and the results of the chi-square test were statistically significant (Figure 7H).

In summary, we speculated that HCC might inhibit the formation of TLSs and evade immune surveillance by downregulating HLA-ABC and upregulating HLA-E and/or HLA-G.