Quan Jiang ,Chenyu Tian ,Hao Wu ,Lingqiang Min ,Hao Chen ,Lingli Chen ,Fenglin Liu,Yihong Sun

1Department of General Surgery,Zhongshan Hospital,Fudan University,Shanghai 200032,China;2Cancer Center,Zhongshan Hospital,Fudan University,Shanghai 200032,China;3Gastric Cancer Center,Zhongshan Hospital,Fudan University,Shanghai 200032,China;4Department of Laboratory Medicine,Zhongshan Hospital,Fudan University,Shanghai 200032,China;5Department of Pathology,Zhongshan Hospital,Fudan University,Shanghai 200032,China

Abstract Objective:Recent studies have highlighted the distinct value of tertiary lymphoid structure (TLS) for immunotherapeutic response prediction.However,it remains unclear whether TLS could play such roles in gastric cancer (GC).Methods:In this study,tumor tissue slices from 292 GC patients from Zhongshan Hospital were firstly reviewed to explore the correlation between TLS and clinical characteristics.Subsequently,we curated 38 reported genes that may function as triggers of TLS and performed consensus molecular subtyping in public RNA-seq datasets to determine TLS patterns in GC.Based on the differentially expressed genes acquired from two TLS patterns,we quantified TLS-related genes on the principal component analysis (PCA) algorithm to develop TLS score.A Zhongshan immunotherapy cohort including 13 patients who received programmed cell death 1 (PD1) blockade therapy was established to conduct RNA sequencing analysis and multiplex immunohistochemistry (mIHC) tests using formalin-fixed and paraffin-embedded (FFPE) tissues.The corresponding TLS score and immune cell counts were further compared based on therapeutic response variations.Results:Mature TLS was revealed as an independent prognostic factor in 292 GC patients.Patients with higher TLS score was characterized by prolonged survival time and superior response to immunotherapy.TLS score was correlated with immunotherapy-related characters,such as microsatellite instability (MSI) and tumor mutation burden (TMB).In addition,RNA-seq data analysis in the Zhongshan immunotherapy cohort indicated that a higher TLS score was correlated with a superior response to PD1 blockade therapy.mIHC tests also revealed that PD1+CD8+ T cell counts were significantly increased in the high-TLS score group.Conclusions:This study highlighted that TLS was significantly associated with immune landscape diversity and complexity.Quantitatively evaluating TLS patterns of individual tumor will strengthen our understanding of TME characteristics and promote more effective immunotherapy strategies.

Keywords: Gastric cancer;tertiary lymphoid structures;immune landscape;immunotherapy

Introduction

Gastric cancer (GC) is the fourth leading cause of cancerrelated deaths worldwide (1).As the major type of infiltrating immune cells,tumor-infiltrating lymphocytes(TILs) have been reported to be related to favorable prognosis in various tumors such as GC,melanoma and nasopharyngeal carcinomas (2-4). Both the cellular composition and organization of tumor infiltrating lymphocytes are crucial for inhibiting cancer progression(5,6).As an important source of TILs,tertiary lymphoid structure (TLS) was characterized by ectopic aggregated lymphocytes with high endothelial venules (7).Considering that TLS was not encapsulated and embedded within the tumor microenvironment,lymphocytes in TLS might have easier access to encounter tumor antigens (8).Recent studies suggested that B cells and CD8+T cells in TLS were associated with better prognosis in GC patients (7,9).Surgery is the most important and effective method for the treatment of gastric cancer (10).Immune checkpoint blockade (ICB) therapy has provided a new approach for cancer therapy in recent years.However,not all patients could benefit from ICB therapy according to the results of relevant clinical trials (11).It is urgently needed for us to explore adequate biomarkers for proper selection of patients who may benefit from ICB therapy.

In this study,we initially uncovered the relationship between TLS and the prognosis of GC from pathological examination.Subsequently,we explored genes associated with the development of TLS in GC and integrated the genomic information to classify patients into two distinct TLS patterns.Immune landscape analysis further suggested that TLS played a crucial role in shaping individual tumor microenvironment characterization.Therefore,we established a set of scoring systems to quantify the TLS pattern.Ultimately,the TLS score was validated and used for guiding the application of immunotherapy based on RNA-seq analysis and multiplex immunohistochemistry(mIHC) tests in the Zhongshan immunotherapy cohort.

Materials and methods

Tumor hematoxylin &eosin (H&E) slices collection and preprocessing of public expression datasets

Tumor slices from 292 GC patients who underwent surgery from December 2008 to June 2019 were collected from Pathological Department of Zhongshan Hospital.Clinical data and follow-up of the patients were also acquired (last follow-up time: January 2020).The studies involving human participants were reviewed and approved by Ethics Committee of Zhongshan Hospital Affiliated to Fudan University.The participants provided their written informed consent to participate in this study.Written informed consent was obtained from the individual(s) for the publication of any potentially identifiable images or data included in this article.Gene expression data,somatic mutation variations (SNVs) and clinical features of GC samples were retrospectively collected from public datasets of The Cancer Genomic Atlas (TCGA) database(https://cancergenome.nih.gov/).GSE84437,GSE26899,GSE78220 and GSE176307 were also gathered from the NCBI GEO database (https://www.ncbi.nlm.nih.gov/geo/)for further analysis. To render the TCGA cohort comparable with the other GEO cohorts,the expression profile (FPKM values) of TCGA-STAD datasets was transformed into transcripts per million (TPMs) (12-14).TCGA (n=375) and GSE84437 (n=433) were integrated to establish a merged cohort for signature establishment.Batch effects from non-biological technical biases were examined and corrected using the “ComBat” algorithm of“sva” package in R software (Version 3.6.1;R Foundation for Statistical Computing,Vienna,Austria) as previously reported (15-17).

Pathological evaluation and RNA sequencing of Zhongshan cohort

Formalin-fixed and paraffin-embedded (FFPE) tissues derived from primary surgeries of 13 patients who received programmed cell death 1 (PD1) blockade therapies after cancer recurrence were acquired from the Pathology Department of Zhongshan Hospital.Response Evaluation Criteria in Solid Tumors version 1.1 (RECIST 1.1) was applied as the criteria for calculating the immunotherapy response of patients (18).The programmed cell death 1 ligand (PD-L1) expression level was defined as the percentage of PD-L1 positive cells based on an immunochemistry (IHC) test.Three slices for each patient were reviewed by two pathologists.The areas for each slice were acquired using Aperio ImageScope software (Version 12.3.3;Leica Biosystems,Nussloch,Germany) and TLS counts were calculated artificially by pathologists.The qualitative classification divided TLS into mature and immature TLS,as previously reported.TLS with follicles was defined as mature TLS and TLS without follicles was defined as immature TLS (19).

Total RNA was extracted from recently cut 10 mm FFPE sections using the miRNeasy FFPE kit (Qiagen,Valencia,CA) according to the manufacturer’s protocol,using 1-4 sections (10-40 mm) per case depending on the assay.RNA yield and quality were determined by UV absorption on a NanoDrop 1000 spectrophotometer and fragment size was analyzed using the RNA 6000 Nano assay (Agilent Technologies,Santa Clara,CA,USA) run on the 2100 Bioanalyzer.Base call files were converted to fastq format using Bcl2Fastq (Version 1.8.4,Illumina,San Diego,CA,USA).All RNA seq reads were aligned to the human reference genome (GRCh38,release 84) using STAR (Version 2.5.2b).

Consensus molecular clustering of TLS-related genes

A total of 38 genes (CCL2,CCL3,CCL4,CCL5,CCL8,CCL18,CCL19,CCL21,CXCL9,CXCL10,CXCL11,CXCL13,CD200,FBLN7,ICOS,SGPP2,SH2D1A,PDCD1,CD4,CCR5,CXCR3,CSF2,IGSF6,IL2RA,CD38,CD40,CD5,MS4A1,SDC1,GFI1,IL1R1,IL1R2,IL10,CCL20,IRF4,TRAF6,STAT5AandTNFRSF17) related to TLS formation extracted from previous publications were collected to identify different TLS patterns in the merged cohort (20-28). Unsupervised clustering analysis was applied to identify distinct TLS patterns based on the expression of 38 regulators of TLS.The number of clusters and their stability were determined by the consensus clustering algorithm (29). We used the“ConsensuClusterPlus” package to perform the above steps and 1,000 times repetitions were conducted to guarantee the stability of classification (30).

Gene set variation analysis (GSVA) and gene ontology(GO) annotation

To investigate the difference in biological processes between TLS patterns,we performed GSVA enrichment analysis using the “GSVA” R package.As a non-parametric and unsupervised method,GSVA is commonly employed for estimating the variation in pathway and biological process activity in the samples of an expression dataset (31).The gene sets of “c2.cp.kegg.v6.2.symbols” were downloaded from the MSigDB database for running GSVA analysis.An adjusted P＜0.05 was considered as statistically significant.

Immune cell infiltration estimation by single sample gene set enrichment analysis (ssGSEA)

ssGSEA was introduced to quantify the relative abundance of 28 immune cell types in the tumor microenvironment.Special feature gene panels for marking each immune cell type were curated from a recent study (32).The relative abundance of each immune cell type was represented by an enrichment score in ssGSEA analysis and normalized to unity distribution from 0 to 1.The biosimilarity of infiltrating immune cells was estimated by multidimensional scaling (MDS) and Gaussian fitting models.

Quantifying immune response predictor: immunophenoscore (IPS),neoantigens and tumor immune dysfunction and exclusion (TIDE) algorithms

IPS is a superior predictor of the response to anti-CTLA-4 and anti-PD1 regimens,which quantifies the determinants of tumor immunogenicity and characterizes the intratumoral immune landscapes (32).The scoring scheme was developed from a panel of immune-related genes belonging to one of the following four classes: MHCrelated molecules (MHC),checkpoints or immunomodulators (CP),effector cells (EC) and suppressor cells (SC).The IPS results and neoantigens calculations of STAD patients were downloaded from The Cancer Immunome Atlas (TCIA) (https://tcia.at/home).The TIDE algorithm proposed by Jianget al.was utilized to model distinct tumor immune evasion mechanisms (33),including dysfunction of tumor infiltration cytotoxic T lymphocytes (CTLs) and exclusion of CTLs by immunosuppressive factors.A higher TIDE score indicated that tumor cells are more likely to induce immune escape,thus indicating a lower response rate to ICB treatment.

Identification of differentially expressed genes (DEGs)between two TLS patterns

To identify TLS-related genes,we classified patients into two distinct TLS patterns based on the expression of 38 regulators of TLS.The empirical Bayesian approach of the“l(fā)imma” package was applied to determine DEGs between different modification patterns (34).The significance criteria for determining DEGs was set as an adjusted P＜0.001.

Calculations of TLS score

We developed a TLS scoring scheme to quantify the TLS patterns of individual patients by using principal component analysis (PCA).Specifically,the overlapping DEGs identified from different TLS Clusters were selected and employed to perform prognostic analysis for each gene using a univariate Cox regression model.The genes with a significant prognostic impact were extracted for further feature selection by using recursive feature elimination(RFE) with random forest and the 10-fold cross-validation method in the “caret” package.Then we curated the expression profile of the final determined genes to perform PCA analysis.Principal components 1 and 2 were extracted and served as the signature score.We adopted a formula similar to previous studies to define the TLS score: TLS score=∑(PC1i+PC2i) (15).

mIHC staining

mIHC staining was performed on some of the tumor samples.The samples were fixed in 4% paraformaldehyde solution and embedded in paraffin.Slides were made using 4-μm sections of the tumor samples.Deparaffinization and rehydration were performed with xylene and ethanol,respectively,followed by microwave antigen retrieval using heated citric acid buffer (pH=6.0) for 10 min and endogenous peroxidase blocking in 3% H2O2for 20 min.Goat serum (Vector,MP-7451) was used to block nonspecific binding sites. Then,relevant primary antibodies were incubated for 1 h at room temperature,followed by the corresponding secondary antibodies(Vector,MP-7451;MP-7452) for 20 min.Slides were then incubated with fluorescein TSA plus for 10 min,after which microwave antigen retrieval was repeated as described in the above steps until the last antibody was added.After multiplexing,DAPI (Sigma,D9542) was used to stain the nuclei.The antibodies and fluorescent dyes used for multiplexing are listed in.Slices of 12 patients were selected from the Zhongshan immunotherapy cohort(Supplementary Table S1). Responders (n=6) included patients who were diagnosed as complete remission (CR),partial remission (PR) or stable disease (SD) after PD1 blockade therapies while non-responders (n=6) were diagnosed as progression disease (PD).One of the FFPE tissues in the Zhongshan immunotherapy cohort was not available due to procedural reasons.The slides were scanned by a Vectra 3 automated high-throughput multiplexed biomarker imaging system (Perkin Elmer,Waltham,MA,USA) and analyzed using the inform image analysis software (Perkin Elmer,Waltham,MA,USA).

Statistical analysis

The statistical analyses in this study were generated by R software and GraphPad Prism 8.0.1 (GraphPad Software,Inc.,San Diego,CA,USA).GraphPad Prism software was applied in the multiplexed immunohistochemistry staining result interpretation,while the remaining data analysis was conducted with R software.For quantitative data,statistical significance for normally distributed variables was estimated by Student’sttests,and non-normally distributed variables were analyzed by the Wilcoxon rank-sum test.For comparisons of more than two groups,Kruskal-Wallis tests and one-way analysis of variance were used as nonparametric and parametric methods,respectively (35).Two-sided Fisher’s exact tests were used to analyze contingency tables.Kaplan-Meier survival analysis and the Cox proportional hazards model were used to analyze the association between the TLS pattern and prognosis using the R package “Survminer” (0.4.6).The surv-cutpoint function from the “survival” package was applied to stratify samples into high-and low-TLS score subgroups.Patients with detailed clinical information were included and adjusted for confounding factors in the univariate regression model.All comparisons were two-sided with an alpha level of 0.05,and the Benjamini-Hochberg method was applied to control the false-discovery rate (FDR) for multiple hypothesis testing (36).

Results

Correlation between TLS classification and GC prognosis

To investigate the correlation between TLS and clinicopathological features of GC,tumor H&E slices of tumor tissues from 292 patients were reviewed for quantitative and qualitative TLS assessments.The baseline data of patients are shown inTable 1.To reduce the influence of tumor heterogeneity on the reading results,three H&E slices of each patient’ tumor tissues from different sampling sites were randomly selected for retrospective reviewing.The main research goals included counting the number of TLS and classifying TLS according to its morphological variations.Representative pictures are shown as lymphoid aggregates (Supplementary Figure S1A),primary follicles (Supplementary Figure S1B)and secondary follicles (Supplementary Figure S1C).As was mentioned above,primary follicles and secondary follicles were defined as immature TLS while lymphoid aggregates were defined as mature TLS (19,37,38).In addition,the number of different TLS was numerated,including total TLS counts (integral TLS),immature TLS counts(immature TLS) and mature TLS counts (mature TLS).Finally,H&E sections were scanned and the areas of tumortissues were calculated by Aperio ImageScope software(Aperio ImageScope,Leica Biosystems,Nussloch,Germany,version 12.3.3).The number of TLS was divided by the tumor area to obtain the density value.

Table S1 Information about antibodies used in mIHC

Table 1 Baseline data of patients with gastric cancer (N=292)

Setting the medians of integral TLS as the cut-off value(0.0293/mm2),we firstly classified the 292 patients into TLS-high and TLS-low groups.Survival analysis revealed no significant difference between the two groups (P=0.934)(Figure 1A).Then,we used the medians of immature TLS as the cut-off value (0.0121/mm2).Survival analysis also presented with no significant difference (P=0.287) (Figure 1B).Finally,we set the medians of mature TLS as the cutoff value (0.0083/mm2).Survival analysis revealed a significantly improved prognosis in TLS-high patients compared with TLS-low patients (P=0.004) (Figure 1C).To further explore its prognostic predictive value,we further performed Cox regression analysis of various clinicopathological characteristics with TLS.The results suggested that mature TLS was an independent protective factor for GC.The above findings suggested that TLS with different maturation levels might exert different effects on GC.

Mature TLS and clinicopathological features of GC

Since mature TLS was revealed as indicative of superior prognosis in GC,we further explored the correlation between mature TLS and clinicopathological characteristics.A total of 292 patients were grouped according to different clinicopathological characteristics and mature TLS density was compared among the groups.The results showed that the mature TLS density was significantly increased in patients with deeper tumor infiltration (T),tumor in advanced stages (TNM),lymphovascular invasion (LVI),higher tumor grade(Grade),nonelderly patients (≤65 years) and female patients. In addition,mature TLS was significantly increased in mixed and intestinal types compared with diffuse types of GC in Lauren classification.Mature TLS was also significantly increased in proximal cases compared with the middle and distal ones (Figure 2).The above results suggested that mature TLS was closely related to the development of GC.

Figure 1 Mature TLS suggested favorable prognosis for GC.Survival comparison analysis of high-and low-TLS GC grouped by integral TLS (P=0.934) (A),immature TLS (P=0.287) (B) and mature TLS (P=0.003) (C);(D) Univariate and multivariate Cox regression analysis of TLS and various clinicopathological indices.TLS,tertiary lymphoid structure;GC,gastric cancer;AJCC,American Joint Committee on Cancer;LVI,lymphovascular invasion.

Figure 2 Mature TLS comparisons in different clinicopathological subgroups.Mature TLS comparisons in different groups of clinicopathological features,including depth of tumor infiltration (T),lymph node metastasis (N),distant metastasis (M),AJCC TNM stage,LVI,tumor grade,HP infection,Lauren type,age,gender and tumor location.TLS,tertiary lymphoid structure;AJCC,American Joint Committee on Cancer;LVI,lymphovascular invasion;HP,Helicobacter pylori;ns,no significance;＊,P＜0.05;＊＊＊,P＜0.001.

Genetic variation landscape of TLS-related genes in GC

Considering that mature TLS could serve as an indicator of superior prognosis in GC,we hypothesized that genes correlated with TLS development and maturation might play crucial roles in GC development.Therefore,we firstly explored TLS-related genes.By reviewing the previous reports,we identified 38 TLS-related genes and investigated their roles in GC (20-28).Seventy-five of 433(17.32%) samples contained genetic alterations of 38 TLS regulators,primarily including missense mutations and frameshift mutations in TCGA STAD datasets(Supplementary Figure S2A).However,most of the gene mutation rates remained low.The locations of TLS regulators on chromosomes are presented inSupplementaryFigure S2Band copy number variations (CNVs) of TLS regulators are shown inSupplementary Figure S2C.CNVs of TLS regulators were also not prevalent.However,further analysis demonstrated that CCL3,CCL4,CCL18,CCL20,CXCL9,CXCL10,CXCL11,CXCL13,CD200,ICOS,CD4,CCR5,IGSF6,IL2RA and STAT5A were significantly upregulated in tumor samples,whereas CCL21,FBLN7 and TNFRSF17 were significantly downregulated in tumor samples (Supplementary Figure S2D).GO enrichment analyses of 38 TLS-related genes were conducted and the most enriched pathways were all correlated with TLS development (Supplementary Figure S2E).Most of the regulators were also revealed as positively correlated based on Spearman analysis (Supplementary Figure S3). In the absence of significant genomic alterations,the aberrant expression of TLS-related genes in GC suggested that TLS development was more likely to be affected by tumor progression process.

Identification of TLS-related genes and their distinct immune landscapes in GC

To better elucidate unrevealed TLS-related genes,the GSE84437 and TCGA cohorts were integrated into merged cohort (n=808).Based on the expression matrix data of 38 TLS-related genes,unsupervised clustering analysis was performed on the merged cohort and patients were divided into two clusters,namely Cluster A and Cluster B (Consensus clustering algorithm,Supplementary Figure S4A).As shown inSupplementary Figure S4B,the expression level of most TLS-related genes in Cluster A was higher than that in Cluster B.Meanwhile,Cluster A subtype presented with relatively better prognosis (Figure 3A).Two clusters could be distinctly separated based on PCA algorithm,confirming the efficacy of clustering(Figure 3B).To explore the biological molecular variations underlying the two distinct TLS-related subtypes,we performed GSVA enrichment analysis against the Kyoto Encyclopedia of Genes and Genomes (KEGG) gene set(Figure 3C).Compared with TLS Cluster B,TLS Cluster A was significantly enriched in immune activation-related processes,including antigen processing and presentation,the T cell receptor signaling pathway and B cell receptor signaling pathway.The enriched pathways were also correlated with TLS maturation.

To strengthen the research results,we further conducted GO and KEGG functional enrichment analysis.In the KEGG analysis,the cytokine-cytokine receptor pathway,chemokine signaling pathway and other immune-related pathways were enriched in Cluster A (Figure 3D).Meanwhile,the T cell activation,leukocyte cell-cell adhesion,lymphocyte differentiation and other immunerelated pathways were also enriched in the GO analysis(Figure 3E).Infiltrative immune cell analyses also indicated that Cluster A was remarkably rich in immune cell infiltrates except for CD56 dim NK cells and neutrophil(Figure 3F).The results above confirmed that the clusters defined by the expression levels of TLS-related genes were characterized by tumor microenvironment variations.

Role of TLS score in predicting immunotherapeutic benefits

To quantify the TLS cluster variations,the PCA algorithm was applied to develop the TLS score using merged cohort as the training cohort.To identify other TLS-related genes,differential expression analysis was conducted between the two clusters.The differentially expressed genes were further screened based on Cox regression analysis to finally collect TLS-related genes in GC.TLSrelated genes were integrally calculated by the PCA algorithm to obtain the TLS score for each case.The flowchart is presented inSupplementary Figure S5.Inspired by the results above,we explored the practical relevance of the TLS score in independent GC cohorts to verify its prognostic value.In merged cohort,higher TLS score suggested favorable survival (Figure 4A). Then,we conducted survival analysis based on clinical character variations.The results confirmed that higher TLS score presented with favorable prognosis except for T1-2 cases(Supplementary Figure S6A). Besides,no significant difference for pathological characters (T and N stages) was found between high-and low-TLS score subgroups(Supplementary Figure S6B). Cox regression analysis confirmed TLS score as an independent prognostic factor in merged cohort (Supplementary Figure S6C).Furthermore,the TCGA and GSE84437 cohorts were used for internal validations while the GSE26899 cohort was set as external validation.Higher TLS scores were consistently associated with better survival in all GCrelated cohorts (Figure 4B-D).The survival prediction accuracy of this model in training and validation sets was also texted using calibration curve and time-dependent C index curve.Despite that high consistency was observed in calibration curves of all the examined cohorts,the C-index of 3 cohorts (merged,TCGA and GSE84437) were relatively low (Supplementary Figure S7).Therefore,we may come to the conclusion that TLS score was an important prognostic factor for GC.As one of the characteristics acknowledged as indicators for immunotherapy response,microsatellite instability (MSI) states were compared between high- and low-TLS score subgroups.The proportion of MSI-high was higher in GC patients with higher TLS scores than in patients with lower TLS scores in TCGA cohort,while the proportion of MSS was higher in GC patients with lower TLS scores (Figure 4E).The TLS score was also found to increase along with the MSI state (Figure 4F).Tumor mutation burden (TMB)is another well-known indicator for immunotherapy response.As shown inFigure 4G,TMB was positively correlated with TLS score in the TCGA cohort.TMB was also revealed to be significantly higher in the high-TLS score subgroup (Figure 4H).TMB alone has been reported to predict better prognosis in GC (13) while the combination of TMB and TLS score may enhance the predictive efficacy (Figure 4I).To further explore the potential value of the TLS score as an immunotherapy indicator,two algorithms (IPS and TIDE) were respectively applied in TCGA cohort.Based on the cut-off value in merged cohort,TCGA GC patients were divided into high-and low-TLS score subgroups.The high-TLS score subgroup was found to have a higher IPS score in all estimated immunotherapy-responding groups (Figure 4J).As for TIDE algorithm,the TIDE values were negatively correlated with the TLS scores (Figure 4K).The results all indicated that higher TLS scores might predict better therapy responses to immunotherapy.Meanwhile,we also compare the TLS score between the two TLS clusters.The score value of Cluster A was higher as expected(Supplementary Figure S8A).As one of the indicators of immunotherapy,neoantigens were also positively correlated with TLS scores (Supplementary Figure S8B).The expression levels of immune checkpoints were also correlated with the immunotherapy results.Most of the immune checkpoints were significantly higher in high-TLS score subgroup (Supplementary Figure S8C).Regarding immune spectrum,immune cells with anti-cancer function,such as CD8+activated memory T cells,were positively correlated with TLS score while cells with pro-cancer function,such as macrophage M2 cells,were negatively correlated with TLS scores (Supplementary Figure S8D).The above results also suggested that the TLS score might be a potential indicator for immunotherapy in GC.

Figure 3 GC clustering based on 38 TLS regulators in RNA-seq data.(A) Survival analyses for two TLS clusters (P=0.020);(B) Principal component analysis for two TLS clusters;(C) Heatmap shows GSVA scores of representative pathways curated from KEGG datasets in two TLS clusters.GC cohort composition (GSE84437,TCGA cohort) was used as sample annotations;(D) KEGG enrichment analyses for TLS subtype-related genes;(E) GO enrichment analyses for TLS-related genes;(F) Abundance of each TME infiltrating cell in two TLS patterns based on ssGSEA analysis.GC,gastric cancer;TLS,tertiary lymphoid structure;GSVA,gene set variation analysis;KEGG,Kyoto Encyclopedia of Genes and Genomes;TCGA,The Cancer Genome Atlas;GO,gene ontology;ssGSEA,single sample gene set enrichment analysis;＊,P＜0.05;＊＊,P＜0.01;＊＊＊,P＜0.001.

Figure 4 TLS patterns in the role of prognosis and immunotherapy.Kaplan-Meier curves for high-and low-TLS score patient groups in merged cohort (training cohort) (P＜0.001) (A),TCGA cohort (internal validation) (P＜0.001) (B),GSE84437 cohort (internal validation)(P＜0.001) (C) and GSE26899 cohort (external validation) (P=0.028) (D);(E) Proportion distribution of three MSI status for high-and low-TLS score subgroups in TCGA dataset (Chi-square test);(F) TLS score value comparisons among three MSI status (Kruskal-Wallis test);(G) Linear regression of TMB and TLS score in TCGA dataset (Spearman test,R=0.5;P＜0.001);(H) TMB comparison between high-or low-TLS score subgroups (Log-rank test);(I) Prognostic value for TMB and TLS score combination analysis in TCGA cohort (Log-rank test,P＜0.001);(J) IPS was compared between high-and low-TLS score cases in CTLA4-PD1-,CTLA4-PD1+,CTLA4+PD1-,CTLA4+PD1+subgroups (Wilcoxon rank-sum test);(K) Linear regression of TIDE and TLS score (Spearman test,R=-0.54;P＜0.001).TLS,tertiary lymphoid structure;MSI,microsatellite instability;MSS,microsatellite stable;TMB,tumor mutation burden;IPS,immunophenoscore;TCGA,The Cancer Genome Atlas;TIDE,tumor immune dysfunction and exclusion;＊,P＜0.05;＊＊,P＜0.01;＊＊＊,P＜0.001.

To further validate the results,two GEO datasets concerning immunotherapy were used to test the TLS score based on real-world results.In both immunotherapy cohorts,patients with high-TLS scores exhibited clinical benefits and better survival rates,although the P value of GSE176307 was not statistically significant (Supplementary Figure S9).The proportion of patients diagnosed as CR or PR after immunotherapy was also higher in patients with higher TLS scores,which suggested that higher TLS scores predicted better responses to immunotherapy.However,both immunotherapy cohorts were not GCrelated.To make our conclusion more convincing and explain why the TLS score could serve as an immunotherapy indicator,we needed to further validate the results in a GC-related immunotherapy cohort.

Comparisons of PD-L1 expression,pathological TLS tests and TLS score as PD1 blockade therapy indicators in GC

For PD1 blockade therapy of GC,scoring system established based on the percentage of PD-L1-positive cells,such as the Combined Positive Score (CPS) system(39-41),is currently widely applied in clinical practice.However,considering the relatively small sample size,we calculated the PD-L1 expression level based on the percentage instead of the CPS score to render comparability.As mentioned above,the existence of TLS revealed by pathological tests has been reported to be a solid predictor of superior immunotherapy responses(24,26-28).In addition,there were also studies mentioning that TLS with follicles presented with better prognosis(42).Patients were grouped according to their treatment response and compared in five metrics.We found that integral TLS,mature TLS and TLS score were significantly higher in responders (Figure 5A).To further understand the correlation among the five indices,we performed correlation analysis and found that the TLS score showed a significant positive correlation with both the percentage of PD-L1 positive cells and mature TLS.Among them,the TLS score also showed the highest correlation coefficient value with mature TLS (R=0.623,P＜0.05) (Figure 5B).Finally,we plotted ROC curves to compare the predictive efficacy of the 5 indicators for immunotherapy in GC.We found that the TLS score had the highest area under the curve (AUC=0.976),indicating that the TLS score had the best diagnostic efficacy in Zhongshan immunotherapy cohort (Figure 5C).The results initially validated that the TLS score might be a new predictor for PD1 inhibitor therapy response in GC.Considering that the TLS score was also found to show the most significant correlation with the mature TLS count obtained from pathological examination,the TLS score may be a new method to predict the mature TLS in GC.

Immune microenvironment variations of GC in response to high-and low-TLS score subgroups

In our previous work,we found that there were significant differences in immune cell fractions predicted by the ssGSEA-based algorithm in two TLS Clusters.To explain why patients with high TLS score had better treatment responses to PD1 inhibitors,we further investigated the differences in immune cell fractions in response to high and low TLS scores in GC.Six patients with high TLS scores and six patients with low TLS scores were enrolled in the study and paraffin sections were co-stained with mIHC technique for seven indicators (including CD4,CD8,CD20,PD-L1,PD1,CK and DAPI).One of the FFPE tissues in the Zhongshan immunotherapy cohort was not available due to procedural reasons.We found that the number of CD8+T cells was significantly increased in both tumor and normal gastric mucosa tissues in high-TLS score cases,but CD4+T cells and CD20+B cells were not significantly different (Figure 6A,B).We also compared the cell fractions in mature TLS and immature TLS separately.We found that the numbers of CD8+T cells and CD20+B cells were significantly increased in mature TLS for high-TLS score cases (Figure 6C),but no difference was observed in immature TLS (Figure 6D).Considering that we included patients treated with PD1 inhibitors,we further investigated PD1 expression in immune cells as well as PD-L1 expression in tumor cells.The results revealed that PD1-positive immune cells were mostly located around mature TLS.The number of PD1+CD8+T cells was significantly increased in high-TLS score cases (Figure 7A).PD-L1 expression on tumor cells was also significantly increased in high-TLS score cases (Figure 7B).The above results not only suggested that a high TLS score could reflect the difference in tumor-infiltrating immune cell counts within GC,but also revealed a significant increase in PD1+CD8+T cells and PD-L1-positive tumor cells in high-TLS score cases.This partially explained why high TLS scores could predict the effectiveness of PD1 inhibitors for GC,and it could also provide a theoretical basis for our further study of the TLS scoring system.

Figure S1 Representative images of TLS based on pathological examination.(A) Lymphoid aggregates;(B) Primary follicles;(C) Secondary follicles.TLS,tertiary lymphoid structure.

Figure S2 Landscape of regulators of TLS in GC.(A) Seventy-five of 433 GC patients experienced genetic alterations of 38 TLS regulators,with a frequency of 17.32%.The number on the right indicated the mutation frequency for each regulator.Each column represented individual patients;(B) Location of TLS regulators on chromosomes;(C) CNV mutation frequency of 38 TLS regulators.The column represented the alteration frequency.Green dot,deletion frequency;Red dot,amplification frequency;(D) mRNA expression levels comparison of 38 TLS regulators between normal and tumor samples (Student’s t test,＊,P＜0.05;＊＊,P＜0.01;＊＊＊,P＜0.001);(E) GO enrichment analyses of 38 TLS-related genes.TLS,tertiary lymphoid structure;GC,gastric cancer;CNV,copy number variation;GO,gene ontology.

Figure S3 Spearman correlation analysis of 38 TLS regulators.TLS,tertiary lymphoid structure.

Figure S4 Unsupervised clustering analysis of merging cohorts and expression of TLS-related genes.(A) Unsupervised clustering analysis was performed on the merged cohort and divided the patients into two clusters,Cluster A and Cluster B (Consensus clustering algorithm);(B) Heatmap presented with the expression variations of 38 TLS-related genes in Cluster A and Cluster B.TLS,tertiary lymphoid structure.

Figure S5 Flowchart of TLS score establishment.TLS,tertiary lymphoid structure;GC,gastric cancer.

Figure 5 Correlation of TLS-related characters and PD1 inhibitor therapeutic response in Zhongshan immunotherapy cohort.(A)Differences of PD-L1-positive cell proportion,total TLS count,immature TLS count,mature TLS count and TLS score between responders and non-responders (Mann-Whitney U text);(B) Correlation between PD-L1-positive cell proportion,total TLS count,immature TLS count,mature TLS count and TLS score (Spearman test);(C) Predictive efficacy comparisons of five indicators for PD1 inhibitor therapeutic response based on ROC curve.TLS,tertiary lymphoid structure;PD1,programmed cell death 1;PD-L1,programmed cell death 1 ligand;CR,complete remission;PR,partial remission;SD,stable disease;PD,progression disease;AUC,areas under the curve;ROC,receiver operator characteristic curve.

Figure S6 Prognostic value of TLS score.(A) Survival analysis of different TLS score subgroups based on various clinical characters;(B)Proportion comparisons of clinical characters between different TLS score subgroups (Chi-square test,＊＊＊,P＜0.001);(C) Cox regression analysis of TLS score and clinical characters.TLS,tertiary lymphoid structure.

Figure S7 Prediction accuracy examinations of TLS score.Calibration curve and time-dependent C index curve drawn for merged cohort(A),TCGA cohort (B) and GSE84437 cohort (C).TLS,tertiary lymphoid structure;TCGA,The Cancer Genome Atlas;OS,overall survival.

Figure 7 Comparisons of PD1/PD-L1 expression level in response to TLS score variations.Difference in the number of PD1-positive immune cells (A) and PD-L1-positive tumor cells (B) between high-and low-TLS score subgroups (Mann-Whitney U test).PD1,programmed cell death 1;PD-L1,programmed cell death 1 ligand;TLS,tertiary lymphoid structure.

Figure S8 TLS score and immune spectrum of GC.(A) TLS score value comparison of two TLS clusters (Student’s t test);(B) Correlation analysis of neoantigens and TLS score (Spearman text;R=0.23,P=0.038);(C) Expression level comparisons of immune checkpoints between high-and low-TLS score subgroups (Student’s t test);and (D) Correlation analysis of TLS score and immune cell proportions based on various estimated software using TCGA dataset (Spearman text).TLS,tertiary lymphoid structure;GC,gastric cancer;TCGA,The Cancer Genome Atlas;＊＊＊,P＜0.001.

Figure S9 TLS score for immunotherapeutic response prediction in non-GC cohorts.Survival analysis and proportion comparisons of therapeutic results between high-and low-TLS score in GSE78220 (A) and GSE176307 (B).TLS,tertiary lymphoid structure;GC,gastric cancer;CR,complete remission;PR,partial remission;PD,progressive disease;SD,stable disease.

Discussion

TLSs are lymphocyte aggregates that occur in nonhematopoietic organs.Reviewing TLS-related studies of GC in recent years,we found that most of the studies were exploring the correlation between TLS and prognosis.TLS might serve as an indicator of superior prognosis of GC in most previous studies (7,9,43,44).As for the classification of TLS,spatial localization difference was most frequently used.However,there are still many deficiencies in terms of its practicality. Tumor heterogeneity is obvious in GC.In diffuse type GC,the tumor boundary was not easy to identify.Therefore,we classified TLS according to its degree of maturity,and found that mature TLS indicated a favorable prognosis in GC.Considering that most of the previous studies believed that TLS could play an anti-tumor immune function,we speculated that the increase of mature TLS might be a manifestation of the body’s specific immune response against tumors (26).Therefore,genes related to the occurrence and maturation of TLS may exert significant impact on GC progression.

TLS is increasingly considered as a predictive biomarker of responses to other anti-cancer therapies.In gastrointestinal stromal tumors,a high density of TLS has been associated with lower imatinib resistance,recurrence,and a more favorable survival (42).In triple-negative breast cancer (TNBC),the high density of high endothelial venules (HEV,as a surrogate marker of TLS) correlates with the pathologic complete response after neoadjuvant chemotherapy (45).However,the correlation between TLS and tumor immunotherapy attracted the most attention in TLS-related researches.

Studies in recent years have highlighted the obvious relevance between immunotherapy and TLS.Voabilet al.established anex vivotumor fragment platform to dissect the response to PD1 blockade in cancer and the baseline presence of TLS and their components correlated with the capacity of cancers to undergo intertumoral immune cell reactivation (46).Helminket al.performed bulk RNA sequencing and found that B cell markers were the most differentially expressed genes in the tumors (melanoma and renal cell carcinoma) of immune checkpoint blockade respondersvs.non-responders and the accumulated B cells were most frequently located in TLSs (27).Rodriguezet al.hypothesized that tumor-associated TLS could serve as novel targets for antitumor immunity and cancer immunotherapy (47).However,it remained unclear whether TLS was induced by immunotherapy since most of the pathological examinations were conducted after treatment.Inspired by the pathological results,we hypothesized that genes related to the occurrence and maturation of TLS may exert significant impact on GC progression.To better elucidate why TLS was associated with immunotherapy,RNA-seq data of GC were processed to explore genes which might play crucial roles in tumor progression based on prognosis and therapeutic result analysis.Mounting evidences have demonstrated that TLS plays an indispensable role in inflammation as well as antitumor effects through interactions with various TLS regulators (45,46).As most studies have focused merely on a single TLS count or a single relevant regulator,the overall immune cell infiltration characterizations mediated by the integrated roles of multiple TLS regulators have still not been comprehensively recognized (47,48).Identifying the role of distinct TLS patterns in TME cell infiltration will contribute to enhancing our understanding of the TME antitumor immune response.In this study,we identified two distinct TLS patterns characterized by distinct immune phenotypes.The Cluster A was characterized by immune activation and tumor-infiltrating lymphocytes,corresponding to an immune-inflamed phenotype.While the Cluster B was characterized by an immunosuppressive TME.Baseline levels of tumorinfiltrating CD4+/CD8+T cells,NK cells and inflammatory cytokine secretion have been shown to be correlated with the likelihood of an immune response.We also identified that the Cluster A pattern was significantly associated with elevated tumor-infiltrating lymphocytes.Furthermore,the mRNA transcriptome differences among distinct TLS patterns have been proven to be significantly associated with immune-related biological pathways.These differentially expressed genes were considered as TLSrelated signature genes.Similar to the clustering results of the TLS phenotypes,two genomic subtypes were identified based onTLSgenes.These findings again demonstrated that TLS was of great significance in shaping different TME landscapes.The TLS pattern characterized by the immunosuppressive phenotype exhibited a lower TLS score,while the pattern characterized by the immuneinflamed phenotype showed a higher TLS score.

Gastric cancer is a multifactorial,multi-stage complex disease with marked heterogeneity (49).Our study also showed that TLS patterns played a non-negligible role in shaping different immune landscapes,implying that TLS could affect the therapeutic efficacy of ICB.The TLS gene signature with various integrated various biomarkers including mutation load,stromal and immune TME status,could be a more effective predictive strategy for immunotherapy.The clinical application value of the TLS score we established was further examined based on prognosis and therapeutic result analysis.The prognostic analysis of three gastric cancer related cohorts (TCGA,GSE84437 and GSE26899) all came down to the conclusion that patients with higher TLS scores had better survival rates.Besides,the survival analysis of the melanoma cohort (GSE78220) also presented with similar results.

To further validate its practicality,13 patients who received PD1 blockade therapy were enrolled to establish the Zhongshan immunotherapy cohort.Compared with previous studies,we collected tumor samples sourced from primary surgery of cancer patients receiving palliative immunotherapy after recurrence.The corresponding analytical results of their primary surgical tumor samples indicated that TLS counts and TLS scores were both significantly higher in therapy responders.Compared with TLS counts calculated by pathological examinations,the TLS score presented with superior efficacy for therapeutic response prediction.The analytical results not only suggested that TLS could be an indicator for PD1 blockade application in GC,but also showed that the TLS scoring system we established might have even better performance.To further elucidate why the TLS score could serve that role,we performed mIHC staining to compare the infiltrative immune cell variations between the responders and non-responders.The results revealed that PD1+CD8+T cells were more enriched in high-TLS score subgroup.This may partially explain why patients with higher TLS scores presented with higher response rates to PD1 blockade.It is worthwhile to mention that a recent study concerning the immunotherapy of GC also came to a similar conclusion (48).

Conclusions

This study demonstrated that the difference in TLS patterns was a factor that could not be ignored causing the heterogeneity and complexity of the individual tumor microenvironment.The TLS scoring system can effectively predict the prognosis and immunotherapy effect of GC patients.Patients with high TLS score had favorable survival and better immunotherapy response than those with low TLS score.

Acknowledgements

This study was supported by grants from the National Natural Science Foundation of China (No.82172803 and No.82072679) and the 2020 Zhongshan Hospital Clinical Research Special Fund (No.2020ZSLC15).

Footnote

Conflicts of Interest: These authors have no conflicts of interest to declare.