Single-cell and spatial transcriptomics implicate a prognostic function of tertiary lymphoid structures in gastric cancer.

Introduction
Immune checkpoint inhibitors (ICIs) based immunotherapy has been proven to be an effective anti-cancer therapy in a wide range of human malignancies1–4, as ICIs bind to immune checkpoint proteins to relieve the tumor-mediated inhibition and recover T-cell cytotoxic function3,5. Gastric cancer (GC) is a prevalent malignancy of the digestive tract and ranks as the third leading cause of cancer-related mortality globally. Because the early symptoms of GC are insidious, most GC patients are metastatic in the advanced stage at the time of diagnosis, with extremely short lifespans6. In recent years, immunotherapy has brought new hope to treating advanced-stage GC. Large-scale phase III randomized controlled studies such as ORIENT-167, CheckMate6498, and KEYNOTE-8119 have basically determined the status of immunotherapy in the first-line treatment of advanced-stage GC. However, not all patients with advanced-stage GC can benefit from immunotherapy; approximately 25%–40% of patients experience disease progression despite receiving initial first-line therapy with immunotherapy plus chemotherapy or immunotherapy with targeted chemotherapy. How to further screen patients before treatment to exclude and avoid unnecessary treatment for people who do not respond to first-line immunotherapy is one of the urgent challenges to be solved in immunotherapy for advanced-stage GC in the future.
Numerous biomarkers have been proposed to predict the immunotherapies’ prognosis, thereby facilitating the precision-guided therapeutic program and clinical management. GC patients with high microsatellite instability (MSI-H) or mismatch repair-deficient (dMMR) tumors have been reported to benefit from immunotherapy, and MSI can be used as an index for detecting genetic instability within tumors10. Additionally, the tumor cell positive score and combined positive score (CPS), developed based on PDL1 expression, have shown promising prognostic value in predicting tumor immunotherapy efficacy11. The Epstein-Barr virus (EBV12) and tumor mutational burden (TMB13) status have also shown beneficial effects and have been recommended as necessary supplemental testing to determine whether the patients should receive immunotherapy. However, due to the intra- and inter-tumor heterogeneity, the abundance or spatial distribution of these indicators in advanced-stage GC tissues may be heterogeneous. There remains a lack of effective reference values for these prognostic indicators to accurately predict the benefits of immunotherapy, limiting their applicability in advanced-stage GC.
Recently, the tumor tertiary lymphoid structures (TLSs) have received extensive attention as potential mediators of immunotherapy response and prognostic indicators, owing to their crucial role in facilitating the infiltration of immune cells into the tumor microenvironment (TME)14. Resembling lymph nodes, TLSs exhibit unique characteristics, comprising a B lymphocytes (BLCs)-enriched zone with dendritic cells (DCs) surrounded by T lymphocytes (TLCs) follicle exhibiting germinal center characteristics, which potentially wield profound influence over the anti-tumor immune response. Recently, increasing research has shown that TLS status and its prognostic value in tumors are largely determined by their spatial distribution, density, and maturation stage15–17. Ding et al.18 also reported that the distribution and density of TLS can predict clinical outcomes in intrahepatic cholangiocarcinoma, revealing contrasting prognostic implications depending on whether TLSs are located within or outside the tumor. These findings indicated that the spatial distribution of TLS and the TME characteristics are determinants of the immunotherapy prognosis. Hence, a comprehensive understanding of the spatial distribution, abundance, and cellular composition of TLS in tumors is essential for classifying GC responsive to immunotherapy and for informing clinical management strategies. To achieve this, the integration of single-cell RNA sequencing (scRNA-seq) with spatial transcriptomic19 offers a transformative approach, enabling precise dissection of the cellular composition and molecular intricacies with spatial distribution information within the TME across a wide range of human malignancies20–24.
Here, according to the TLS spatial distributions, we enrolled three groups of GC samples: peritumor TLS (pTLS), intratumoral TLS (iTLS), and desert TLS (dTLS) for integrated single-cell and spatially resolved transcriptomic analysis to profile the molecular and cellular heterogeneity of TME within GC with varying TLS status. This approach provided a comprehensive in situ cellular and molecular atlas of GC with different TLS locations. Our study highlighted different TME programs, including key cellular players, pathways, and spatial interactions in GC with different TLS status. We also defined a single-cell/spatial TLS (ssTLS) signature comprising critical cellular components and features enriched in intratumoral TLS GC (iTLS-GC). The ssTLS signature offered an effective biomarker panel to predict the prognosis of GC patients undergoing immunotherapy.

Results

A single-cell and spatially resolved atlas of human GC with varying TLS status
In this study, we performed large-scale human GC tissue pathology examination and survival follow-up analysis (n = 110). Our analysis revealed that GC demonstrates heterogeneity with different degrees of immune cell infiltration and varying TLS status (Fig. 1a, b, Supplementary Fig. 1a, b). GC can be divided into desert TLS GC (dTLS-GC), peritumoral TLS GC (pTLS-GC) and iTLS-GC, based on the spatial distribution of TLS (Fig. 1a, Supplementary Fig. 1b). Kaplan–Meier survival curves showed that iTLS-GC patients exhibited significantly improved overall survival (OS) and progression-free survival (PFS) (log-rank test, p < 0.05) compared with pTLS-GC and dTLS-GC patients (Fig. 1b). We also performed correlation analysis between GC immunotherapy efficacy and TLS spatial distribution status, and found that the iTLS-GC had a better immunotherapy response compared with pTLS-GC and dTLS-GC (Supplementary Fig. 1c), suggesting that the TLS spatial distribution may be a key factor influencing the prognosis of GC. To investigate the single-cell and spatial architecture corresponding to varying GC TLS status, we performed unbiased scRNA-seq on 25 GC specimens and spatial transcriptomic sequencing (ST-seq) on five GC specimens (Fig. 1c, Supplementary Fig. 2a–c, Supplementary Data 1). All enrolled samples were collected from GC patients without chemotherapy or radiotherapy. For scRNA-seq, a total of 88,100 cells were adopted for further bioinformatics analysis following rigorous quality control, which can be decomposed into 11 major cell types (Fig. 1d) based on the expression of specific cell-type markers (Fig. 1e, Supplementary Data 2) and copy number variation analysis (Supplementary Fig. 2d, e). The cell-type density distribution analysis (Fig. 1f) and relative cell-type abundance analysis (Fig. 1g) demonstrated that iTLS-GC was significantly enriched with endothelial cell (Endo), BLC, and the proliferative BLC and TLC. By comparing the differentially expressed genes (DEGs) between the malignant epithelium (mEPI) and benign epithelium (bEPI), we found mEPI specifically expressed CDLN3/4/7, KRT17, and EPCAM (Fig. 1h). Gene set enrichment analysis (GSEA) showed that the upregulated DEGs in mEPI compared to bEPI were enriched in pathways associated with tumor progression such as MYC targets V1, epithelial-mesenchymal transition, hypoxia, and tumor invasiveness (Fig. 1i). Five GC samples including two iTLS-GC, two pTLS-GC, and one dTLS-GC were enrolled for ST-seq (Fig. 2a), unbiased clustering of all spatial spots revealed 13 major ST clusters (Fig. 2b–d). The spatial distribution of malignant epithelium score was consistent with the tumor region confirmed by hematoxylin and eosin (H&E) staining images (Fig. 2e) and the major cell-type spatial localization of five GC samples were shown in Fig. 2f and Supplementary Fig. 3. ST cluster 4 expressed high level of BLC specific gene MS4A1, CD37, and TLS-related cytokines such as CCL19, CXCL13, and CXCR4, indicative of tumor TLS region, and ST cluster 4 was mainly distributed in iTLS-GC and pTLS-GC samples, while they were almost absent in dTLS-GC (Fig. 2g). For iTLS-GC sample, the spatial cell-type abundance of TLC, BLC, proliferative BLC and TLC showed high correlation coefficient with the spatial TLS score (Fig. 2h).

Spatial and single-cell dissection of TLC and NK cell diversity of GC with varying TLS status
We further performed unsupervised clustering of TLC and NK cells to reveal the intrinsic structure and molecular features of TLC and NK cells in GC with varying TLS status. A total of two NK cell subclusters, five CD8+ TLC subclusters, and four CD4+ TLC subclusters were identified (Fig. 3a–c, Supplementary Data 3). Among them, we found that cluster CXCL13+CD8+, GZMB+CD8+, and CXCL13+CD4+ TLC exhibited high expression level of the immune cell-homing signal gene CXCL13, ITGAE (encoding CD103), and cytotoxic effector gene GZMB (Fig. 3b, c), which have been reported to be associated with a favorable outcome in immunotherapy25–27. The cell composition of TLC and NK cell subsets varied across GC TLS status. Detailed, naïve IL7R+CD4+ TLC was specifically enriched in dTLS-GC, while iTLS-GC and pTLS-GC showed dramatic enrichment of CXCL13+CD8+, and CXCL13+CD4+ TLC (Fig. 3d, e). The T-cell receptor (TCR) clone expansion analysis revealed that GZMB+CD8+ TLC, CXCL13+CD8+ TLC, and ZNF683+CD8+ TLC possessed the highest hyperexpanded TCR clonotypes, and CXCL13+CD8+ TLC shared the most expanded TCR clonotypes with ZNF683+CD8+ TLC, indicating close relationships in development hierarchy and crucial roles in immune infiltration (Supplementary Fig. 4a–c). Moreover, slingshot pseudo-time trajectory analysis supported the hypothesis28 that cytotoxic effector CXCL13+CD8+ TLC originate from the differentiation of naïve IL7R+CD8+ TLC and resident ZNF683+CD8+ TLC (Supplementary Fig. 4d). We further compared the DEGs of TLC in GC samples with varying TLS status and found that CXCL13, ITGAE, GZMA, GZMB, GZMK and the immune checkpoint molecule PDCD1 were highly expressed in iTLS-GC than pTLS-GC and dTLS-GC (Fig. 3f, g, Supplementary Fig. 4e). To validate the scRNA-seq findings, we utilized ST-seq to profile the spatial distribution pattern of TLS and CXCL13+ TLC in five GC samples with varying TLS status. The in situ spatial transcriptomic analysis revealed a pronounced spatial expression of CXCL13+CD4+ and CXCL13+CD8+ TLC in the peritumoral and intratumoral regions of iTLS-GC samples and peritumoral regions of pTLS-GC samples, while these cell types were rarely expressed in dTLS-GC (Fig. 3h). Interestingly, the spatial distribution of CXCL13+CD4+, CXCL13+CD8+ TLC was highly correlated with the TLS spatial distribution (Fig. 3h, Supplementary Fig. 4f). We further performed the survival analysis in GC patients with anti-PD1 immunotherapy and demonstrated that patients with higher abundance of CXCL13+ TLC possess a favorable immunotherapy response and outcome (Fig. 3i) in terms of PFS and OS.

Spatial and single-cell dissection of BLC diversity of GC with varying TLS status
BLCs, the major cellular component of infiltrating immune cells within TLS across various cancers, participated in anti-tumor immunity by antigen presentation and formation of TLS, which are predictive of good prognosis for immunotherapy. Detailed sub-clustering analysis revealed four BLC subclusters: germinal center BLCs (gc_B, marked by MS4A1, CD79A, and CXCR5), proliferative gc_B (gc_B_MKI67, marked by MS4A1, MKI67, and CXCR5), naïve BLCs (naïve_B, marked by FCER2), and memory BLCs (memory_B) with high expression of CD27 and negative expression of CD38. Additionally, we identified four SDC1+ (encoding CD138) plasma BLCs subclusters: IGA_B (IGHA1), IGG_B (IGHG1), IGM_B (IGHM), and IGJ_B (IGJ) (Fig. 4a, b, Supplementary Data 4). The proportion of these BLC subtypes differed notably among GC with varying TLS status; the cluster gc_B_MKI67, gc_B, and IGA_B were predominantly found in iTLS-GC, while IGG_B and IGM_B were more abundant in pTLS-GC (Fig. 4c, d). Functional enrichment analysis revealed that gc_B_MKI67 and gc_B were enriched with high levels of immune checkpoint pathways such as PD1 and CTLA4 (Fig. 4e). Further intercellular crosstalk analysis revealed that CXCL13+ TLC interacted with gc_B_MKI67 and gc_B through the CXCL13-CXCR5 ligand–receptor pair (Fig. 4f). The multiplex immunostaining in the iTLS-GC sample further validated that CXCL13+ TLC co-existed with CXCR5+ gc_B (Fig. 4g), which is associated with BLC recruitment and germinal center formation. Additionally, the ST-seq result also revealed that gc_B had a similar spatial distribution pattern with CXCL13+CD4+ and CXCL13+CD8+ TLC in GC with varying TLS status (Fig. 4h, Supplementary Fig. 5a). Interestingly, we found that the spatial distribution of ligand CXCL13 secreted by CXCL13+ TLC colocalized with the receptor CXCR5 secreted by gc_B, and the gc_B-associated signature genes MS4A1, CD19, and CXCR5 showed similar spatial distribution or expression pattern with TLS (Fig. 4i, Supplementary Fig. 5a, b). We further performed the survival probability analysis in GC with anti-PD1 immunotherapy, and found that the abundance of gc_B_MKI67 and gc_B was related to better immunotherapy outcome, in terms of PFS and OS (Fig. 4j).

Spatial and single-cell dissection of myeloid cell diversity of GC with varying TLS status
We next performed detailed sub-clustering analysis to decipher the molecular heterogeneity and crucial roles of myeloid cells in GC with varying TLS status, which generated eleven subclusters: mast cell (Mast) featured by CPA3, neutrophil (Neu) featured by CXCR2, monocytes (Mono) featured by VCAN, four tumor-associated macrophages (TAM) subclusters: CCL18+ TAM (CCL18), SPP1+ TAM (SPP1), FOLR2+ TAM (FOLR2), and C1QC+ TAM (C1QC); four DC subclusters: CLEC10A+ DC (CLEC10A), CD1C+ DC (CD1C), LAMP3+ DC (LAMP3) and LILRA4+ DC (LILRA4) (Fig. 5a, b, Supplementary Fig. 6a, Supplementary data 5). We found the proportion of Neus, Masts, and C1QC+ TAM were significantly increased in pTLS-GC compared to dTLS-GC and iTLS-GC, whereas LAMP3+ DC was more enriched in iTLS-GC compared to dTLS-GC (Fig. 5c). The gene set variation analysis (GSVA) revealed that LAMP3+ DC exhibited high levels of immune checkpoint pathway such as PD1, 41-BB, CTLA4, NECTIN2, and CD28 co-stimulation pathway (Fig. 5d, Supplementary Fig. 6b). Additionally, the functional enrichment analysis revealed that LAMP3+ DC exhibited high M1 score, but also migration and activation scores (Supplementary Fig. 6c), which were reported to serve as mature DC and involved in TLC activation and maturation29. The intercellular crosstalk analysis demonstrated that LAMP3+ DC interacted with TLC through PDCD1-CD274 (PD1-PDL1), CXCL16-CXCR6, CCL19-CCR7, CD80-CD28, and LGALS9-HAVCR2 ligand–receptor pair (Fig. 5e, Supplementary Fig. 6d), both of which were important predictors of immunotherapy response. The ST-seq further confirmed the enrichment of LAMP3+ DC in intratumoral regions of iTLS-GC and peritumoral regions of pTLS-GC (Fig. 5f). The spatial distribution of ligands CXCL16, CCL9, CD274 and CD80 secreted by LAMP3+ DC colocalized with the receptors CXCR6, CCR7, PDCD1, and CD28 secreted by CXCL13+ TLC (Fig. 5g, Supplementary Fig. 6e). Moreover, the LAMP3+ DC was found to colocalized with CXCL13+ TLC and TLS, indicating a significant role of LAMP3+ DC in crosstalk and function regulation with TLC (Fig. 5h). The multiplex immunostaining also validated the co-existence of CD80+PDL1+ DCs with PD1+CD4+ and PD1+CD8+ TLC (Fig. 5h), indicating the critical roles of LAMP3+ DCs in TLS formation. Additionally, survival probability analysis further demonstrated that the presence of LAMP3+ DC was associated with a favorable immunotherapy response in terms of OS and PFS (Fig. 5i).

Spatial and single-cell dissection of stromal cell diversity of GC with varying TLS status
Like secondary lymphoid organs (SLOs), TLS functions are largely supported by a specialized network of stromal cells that share some of the specialization and plastic features ascribed to the SLO network. To elucidate the diversity and functional characteristics of stromal cells in GC with varying TLS status, we performed sub-clustering analysis on cancer-associated fibroblasts (CAFs) and endothelial cells (Endos). The CAF can be classified into three subsets of matrix CAF (mCAF) characterized by high expression of ECM genes containing PDGFRA, and COL8A1, three vascular CAF subclusters featuring the expression of microvasculature genes NOTCH3, GJA4, and MYH11, and inflammatory CAF (iCAF) identified by high expression of inflammatory chemokine genes CXCL5/8, and IL11/24 (Supplementary Fig. 7a, b, Supplementary Data 6). Among three mCAF subclusters, the mCAF1 was enriched in iTLS-GC (Supplementary Fig. 7c), featured with a high expression level of CXCL12, complement genes (C3 and C7) and HLA family member genes (HLA-DRA, HLA-DRB1, and HLA-DPA1), and the expression level of these genes were strikingly elevated in cells of iTLS-GC (Supplementary Fig. 7d). The gene ontology (GO) enrichment analysis indicated that the mCAF1 enriched significant GO terms related to collagen fibril organization and extracellular matrix remodeling, MHC class II antigen presentation, and regulation of leukocyte adhesion (Supplementary Fig. 7e). The intercellular crosstalk analysis revealed that CXCL12-CXCR4 and HLA-CD4/8 were the most frequent ligand–receptor pairs to guide the interaction between mCAF1 and TLC (Supplementary Fig. 7f).
We next divided Endos into six subclusters (Endo1-Endo6, Fig. 6a). The cluster Endo1 was characterized by a high expression level of enriched ACRK1, P-selectin (SELP) and E-selectin (SELE), and cluster Endo6 expressed high level of LYVE1, PDPN, and CCL21 (Fig. 6b, Supplementary Fig. 8a, Supplementary Data 7), both of which were notably enriched in iTLS-GC (Fig. 6c). The GO enrichment analysis indicated that the Endo1 depicted significant immune cell recruitment related biological process such as positive regulation of leukocyte migration, leukocyte adhesion, migration, and extravasation processes, and the cluster Endo6 was enriched with lymph vessel development and lymph angiogenesis (Fig. 6d, Supplementary Fig. 8b). We further performed transcriptional programs similarities analysis between MECA79-positive tumor high endothelial venules (TU_HEV, GSE19806930) and the cluster Endo1. Then, we annotated Endo1 as HEV since Endo 1 has the strongest transcriptional program correlation with TU_HEV (Supplementary Fig. 8c–e). Multiplex immunostaining revealed that HEVs colocalized with CXCL13+CD8+ TLC in the intratumoral regions of iTLS-GC and at the peritumoral region of pTLS-GC (Fig. 6e). Interestingly, the Endo1/HEV cells engaged in more frequent interactions with immune cells through the ACKR1-CXCL13/CCL5/CCL17, SELP/SELE-SELPG, NECTIN2-TIGIT, and LGALS9-HAVCR2 pathway, while Endo5 and Endo6 engaged in more frequent interactions with immune cells through CXCL12-CXCR4 and CCL21-CCR7 pathway (Fig. 6f, Supplementary Fig. 8f). We further found that VCAM1, ACKR1, and SELP exhibited higher expression level in Endo1 from iTLS-GC and pTLS-GC than dTLS-GC, especially in iTLS-GC (Fig. 6g). The above results revealed that the Endo1/HEV cells likely play a prominent role in facilitating immune cell recruitment and infiltration in iTLS-GC. The ST-seq further demonstrated that the spatial distribution of Endo1/HEV cells well corresponding to the TLS, and CXCL13+CD8+ TLC (Fig. 6h). Specifically, the ligands SELPLG, CXCL13, CCL5, and CCL17 secreted by CXCL13+ TLC colocalized with the receptors ACKR1 and SELE secreted by Endo1/HEV cells (Fig. 6i, Supplementary Fig. 8g). The above results provided evidence that Endo1/HEV may promote the formation of TLS by recruiting CXCL13+ TLC into the TLS region. To investigate the role of the ACKR1-CXCL13 axis in mediating TLC recruitment and infiltration, ACKR1-overexpressing ECs (ACKR1OE vECs) were generated through an in vitro differentiation assay. The ACKR1OE vEC were treated for 24 h with PBS (control), 50 ng/mL CXCL13, 5 ng/mL TNF, or a combination of both CXCL13 and TNF. Treatment with CXCL13, alone or in combination with TNF, significantly upregulated the expression of the adhesion molecules VCAM1 and ICAM1 (Fig. 6j). Moreover, TLC migration assay demonstrated that CXCL13-treated ACKROE vECs recruited more TLC compared to PBS-treated controls, confirming the functional role of CXCL13 in enhancing TLC recruitment (Fig. 6k). Collectively, these findings support that ACKR1-expressing Endo1/HEV cells facilitates the recruitment of CXCL13+ TLC by reciprocally upregulating VCAM1 and ICAM1.

A unique TLS signature is indicative of potential immunotherapy response benefit in GC
Although TLS has been consistently associated with a favorable prognosis in various cancers, conflicting prognostic impacts have been reported in some studies. This discrepancy is potentially attributed to the spatial distribution of TLS among individuals. A comprehensive assessment of the spatial distribution, abundance, and cellular composition of the GC TME based on different TLS states is crucial for exploring precision molecular features linked to the immunotherapeutic response in GC. Our investigation, which identified key cellular players and molecular features in GC with varying TLS spatial distributions, enabled further exploration of candidate cellular characteristics for predicting the immunotherapy response in GC. Pan-cancer analysis across a pan-cancer dataset (https://kmplot.com) of key TLS-related genes demonstrates that high levels of TLS-associated genetic signatures, including chemokines (e.g., CCL5, CXCL13, CCL19), receptors (e.g., CXCR5, CCR7), regulatory molecules (e.g., ACKR1, SELE) are robustly associated with better survival with cancer patients receiving immunotherapy (Fig. 7a, Supplementary Fig. 9a). We observed a higher abundance of CXCL13+ TLC, CXCR5+ gc_B, LAMP3+ DC, and SELP+ACKR1+ Endo1/HEV in anti-PD1 immunotherapy-responsive patients, while SPP1+ TAM and iCAF were more abundant in anti-PD1 non-responsive GC cases (Supplementary Fig. 9b). We also found that the immunotherapy response efficacy was influenced by the patients’ pathological features such as TMB, microsatellite instability (MSI), EBV and PDL1 CPS (Supplementary Fig. 10a). By adjusted the above pathological features, the relative cell-type abundance of CXCL13+ TLC, CXCR5+ gc_B, LAMP3+ DC, and SELP+ACKR1+ Endo1/HEV demonstrated a significant favor survival outcome in GC patients with anti-PD1 immunotherapy (Fig. 7b, Supplementary Fig. 9c, d). Importantly, based on the cell-type characteristics identified from scRNA-seq and ST-seq, we defined a single-cell & spatial TLS (ssTLS) signature, which included not only traditional TLS-specific signature genes (CD3D, CD3E, CD19, MASA1, and CD79A), but also the immune-related signature genes (CXCL13, CXCR5, CXCL16, CCL19, CCL22, ACKR1, and SELE). The ssTLS signature spatial distribution was highly consistent with the TLS region distribution (Fig. 7c). Moreover, the ssTLS signature score was elevated in anti-PD1 immunotherapy-responsive GC patients who exhibit high immune signature, CPS scores, and EBV-positive status (Fig. 7d, Supplementary Fig. 10b), and was positively associated with prolonged PFS and OS in GC patients receiving anti-PD1 therapy (Fig. 7e, f). Importantly, the ssTLS signature effectively distinguished anti-PD1 responsive cases from non-responsive cases, which achieved an impressive area under the curve (AUC) of 0.795 (Fig. 7g). In summary, our integration of single-cell transcriptional results of GC with public bulk RNA-seq data from immunotherapy-treated GC identified key cell players and molecules serving as indicators to predict the immunotherapy response in GC. This approach lays the foundation for individualized GC therapy.

Discussion
Immunotherapy has shown its benefits for cancer patients; several biomarkers, such as PDL1, TMB, MSI/MSS, have been proposed to evaluate tumor immunotherapy benefit and predict treatment prognosis4,7,9,31. However, the prognostic value of these indicators and the treatment outcomes can vary significantly among individuals due to tumor heterogeneity. Various research studies have demonstrated that TME characteristics were strongly associated with the immunotherapy efficacy and patients’ outcomes. As an important player in anti-tumor immune response, TLS in the TME can promote immune infiltration and adaptive immune response, leading to the clearance of tumor cells and improving the prognosis of most cancer patients, which has been proven to be associated with improved survival rates across various cancer types32. However, in liver cancer, TLSs have been correlated with poor prognosis18,33. The dual prognostic roles of TLS in tumor prognosis suggest underlying heterogeneity, potentially attributable to factors such as density, composition, maturity stages, and spatial distribution within intra- and peri-tumor regions of TLS34. A comprehensive and deep understanding of the single-cell and spatial characteristics of TME in GC with different TLS status is essential for identifying key cellular players and molecular pathways, and biomarkers with prognostic values in guiding GC immunotherapy.
In this study, through large-scale GC TLS status and follow-up analysis, we found that iTLS-GC and pTLS-GC exhibited better survival and immunotherapy response than those with dTLS-GC. And iTLS-GC had an even more favorable outcome compared to pTLS-GC, which suggests that both the presence and spatial localization of TLSs contribute to immunotherapy efficacy in GC patients. By integration of scRNA-seq and ST-seq survey on GC with varying TLS status, we uncovered that the iTLS-GC showed significant enrichment of CXCL13+ TLC, gc_B, LAMP3+ DC, and HEV cells compared with pTLS-GC and dTLS-GC. The spatial distribution of these cell types colocalized with the TLS-enriched regions and differed largely in GC with varying status, which indicated that the TLS spatial distribution status had a significant influence on the cellular and molecular diversity of GC. We also uncovered important multi cell–cell interaction and important ligand–receptor pairs involved in TLS formation and tumor TLS spatial distribution status. Detailed, we uncovered the Endo molecular diversity and found the Endo1 termed as HEVs expressed high levels of SELP, SELE, ACKR1, and VCAM1. The Endo1/HEV enriched GO terms with immune cell recruitment related biological processes, such as positive regulation of leukocyte migration, leukocyte adhesion, migration, and extravasation processes. Previous studies have reported that HEVs express key genes involved in cell recruitment and immune defense, serving as important drivers to promote the infiltration of TLCs and BLCs into tumors35–37. Our results also found several key chemokine receptor-ligand interactions, such as ACKR1-CXCL13/CCL5, and SELP/E-SELPG between HEVs and CXCL13+ TLC. These interactions play a significant role in promoting intra-tumor immune cell infiltration in GC. The low abundance of HEV cells within GC tissues posed technical challenges in isolating enough for downstream functional validation. Although our in vitro experiments demonstrated that CXCL13-ACKR1 signaling in ACKR1OE EC led to the reciprocal upregulation of adhesion molecules VCAM1 and ICAM1, thereby promoting the recruitment of CXCL13+ TLC. We were still unable to functionally evaluate the contribution of these ligand–receptor pairs, including the CXCL13-ACKR1 axis, to TLS formation in vivo due to the lack of appropriate mouse models. Future studies utilizing genetic or inducible mouse models will be crucial to elucidate the mechanistic role of these pathways in HEV differentiation and TLS development in the gastric TME.
Previous researchers have identified a range of cellular and molecular biomarkers, including diverse types of tumor-infiltrating immune cells, immune checkpoint-related and inflammation-associated gene signatures, and genomic characteristics, for predicting clinical outcomes38,39. However, due to the complexity of TME and the dynamic nature of therapeutic targets that evolve over time and across different tumor locations, a single prognostic factor may exhibit dual prognostic roles across different cancer types. TLS, which comprises numerous tumor-infiltrating immune cells and inflammation-related molecules, has demonstrated a strong correlation between TLS existence and superior prognosis in most cancers, as evidenced by numerous retrospective cohort studies primarily evaluated by immunohistochemistry and H&E staining40. Recent research has defined TLS signatures based on CXCL13 expression in breast cancer and melanoma, or 12-chemokine-related genes in colorectal cancer, hepatocellular cancer, and melanoma41–43. The application of a 29-gene signature, including globulin genes, BLC marker genes, TLC marker genes, fibroblast marker genes, and complement genes, has increased the accuracy of defining TLS imprint signature43. Having defined the cell players and molecular pathways associated with immune infiltration and infiltrated TLS formation by leveraging single-cell and spatial transcriptomics, we proposed a ssTLS signature. This signature integrates characteristics containing immune checkpoint and inflammatory molecules enriched in iTLS-GC TME cell features, providing insights into TLS spatial architecture and predicting GC immunotherapy prognosis with high accuracy.
Our study still has several limitations. First, the relatively modest sample size and inherent interpatient heterogeneity of GC may limit the generalizability of our findings. Second, the resolution of 10X Visium spatial transcriptomic methods is insufficient to resolve the heterogeneity within TLS at the single-cell level. This limitation may introduce inaccuracies in deconvolution analyses that rely on single-cell transcriptomic data. Additionally, 10X Visium spatial transcriptomic captures the 2D spatial information by relying on the tissue slices to capture mRNA, thereby ignoring the intrinsic heterogeneity inherent in the 3D structure of TLS. Further advanced image-based spatial transcriptomic techniques leveraging tissue optical clearing and cell expansion methods are expected to reconstruct the 3D structure of TLS within GC TME. Finally, although we integrated immunostaining and in vitro experiments to validate our transcriptomic results, screened potential ligand–receptor pairs to validate our findings, further investigation incorporating genetic models and functional assays is warranted to elucidate the mechanistic contribution of identified signaling pathways to TLS formation and their therapeutic implications. In summary, our study deciphered the molecular and cellular heterogeneity of the TME within GC with varying TLS spatial distributions. The SELP+ACKR1+ HEVs enriched in iTLS-GC potentially recruit the CXCL13+ TLC via the CXCL13-ACKR1 pathway. These CXCL13+ TLC, in turn, subsequently lead to the attraction of CXCR5+ gc_B recruitment via the CXCL13-CXCR5 pathway, suggesting a potential cellular mechanism underlying the formation of TLS. Based on the cell-type characteristics identified from scRNA-seq and ST-seq, we defined an ssTLS signature, which effectively distinguished anti-PD1 responsive cases from non-responsive cases, and achieved an impressive AUC of 0.795. Our research contributes to a deeper understanding of the spatial architecture of TLS and aids the assessment of GC immunotherapy prognosis and the development of optimal treatment strategies.

Methods

Ethics approval
This research was carried out in compliance with the principles outlined in the Declaration of Helsinki. It was granted clearance by the Ethics Committee of the Fourth Medical Center of the PLA General Hospital (approval numbers: 2021KY011-HS001 and 2024KY0134-KS001). Before their participation, all individuals involved in the study provided acknowledging and agreeing to the study’s terms and conditions. All participants participate in the study without compensation.

Clinical cohorts, immunotherapy response, and prognosis of GC
This retrospective study was conducted after obtaining approval from the ethics committee. A total of 110 patients were enrolled to assess the role of TLS status on the prognosis of GC. Additionally, 39 GC patients who received immunotherapy and had detailed follow-up information were enrolled to assess the role of TLS status on the drug efficacy. The inclusion criteria for the prognosis-analysis cohort were as follows: (1) patients who underwent gastric cancer surgery at the Fourth Medical Center of the PLA General Hospital between December 2017 and December 2019; (2) patients who had not received any preoperative treatment at the time of surgery; (3) postoperative pathological specimens included at least three H&E-stained slides containing both tumor tissue and adjacent non-tumor tissue for the evaluation of TLS spatial distribution; and (4) patients who could be contacted via telephone to obtain data on PFS and OS. The exclusion criteria for the prognosis-analysis cohort were as follows: (1) patients who had received neoadjuvant therapy before gastric cancer surgery; (2) patients with missing postoperative pathological specimens (fewer than three slides) or without adjacent non-tumor tissue available for evaluation; and (3) patients who could not be followed up. Inclusion criteria for immunotherapy response cohort were as follows: (1) patients who underwent gastric cancer surgery at the Fourth Medical Center of the Chinese PLA General Hospital between 2020 and 2024; (2) postoperative pathological specimens included at least three H&E-stained slides containing both tumor tissue and adjacent non-tumor tissue for evaluating the spatial distribution of TLS; (3) received immunotherapy (including sintilimab, tislelizumab, toripalimab, pembrolizumab, camrelizumab) after surgery, with or without the use of other treatments (including targeted therapy and chemotherapy); (4) received immunotherapy for at least three times and had evaluable efficacy. Exclusion criteria or immunotherapy response cohort were as follows: (1) GC surgical patients who did not receive immunotherapy; (2) patients with missing postoperative pathological specimens (less than three slides), without tumor tissue or adjacent non-tumor tissue available for evaluation; (3) patients whose immunotherapy efficacy could not be evaluated. For TLS classification, three pathologists independently evaluated the spatial distribution of TLS in tumor and adjacent tissues based on H&E staining, with each tumor sample assessed using at least three slides. TLSs were categorized as follows: (1) iTLS: defined as TLS present within the tumor boundary on H&E staining images; (2) peritumoral TLS (pTLS): defined as TLS present outside the tumor boundary but not within it; (3) desert TLS (dTLS): defined as H&E staining images without any detectable TLS. If both pTLS and iTLS co-existed (mixed-TLS), meantime, the overwhelming majority of TLS in mix TLS GC was located in the intra-tumor region, the sample was classified as iTLS. Survival definitions: PFS: the time interval from the date of surgery to the date of tumor progression or recurrence; OS: the time interval from the date of surgery to the date of death from any cause. For patients who remained progression-free after surgery, PFS and OS were calculated as the time interval from the date of surgery to the date of the last follow-up. The Kaplan–Meier survival curves were plotted using the R package “survminer” (v0.4.9). Efficacy assessment: complete response (CR): all measurable and non-measurable lesions disappeared completely, with no new lesions, and normal tumor markers, lasting for at least 4 weeks. Partial response (PR): the sum of the longest diameters of target lesions decreased by ≥ 30% compared to the baseline level, with no new lesions, no progression of any non-measurable lesions, and no deterioration of tumor markers, lasting for at least 4 weeks. Stable disease (SD): does not meet the criteria for CR, PR, or PD. Progressive disease (PD): the sum of the longest diameters of target lesions increased by ≥ 20% compared to the smallest measurement value during treatment, and the absolute value increased by ≥ 5 mm, with no new lesions; the appearance of one or more new lesions, regardless of changes in target lesions; non-target lesions showed clear progression, regardless of the situation of target lesions and new lesions. CR and PR were considered as having an immune therapy response (R), while SD and PD were considered as having no immune therapy response (NR).

Acquisition of fresh tissue specimens and preparation of a single-cell suspension
This study collected fresh tumor specimens from 25 patients with GC patients who underwent gastrectomy at the Fourth Medical Center of the PLA General Hospital from April 2021 to March 2023. Following surgical excision with a scalpel, GC tissues were acquired. Subsequently, the freshly obtained tissue specimens were thoroughly rinsed with phosphate-buffered saline and subsequently divided into two distinct sections. One section was processed to generate a single-cell suspension, while the other section was allocated for additional experiments, including histological analysis through 10X Visium spatial transcriptomics, H&E staining, and immunofluorescence staining techniques. The distribution of TLS within cancerous and adjacent tissues was evaluated by three pathology experts through the analysis of H&E staining and immunofluorescence staining images. Each tumor sample was evaluated with at least three H&E staining images containing cancer tissue and adjacent tissue.

Sample processing
The tissue samples designated for single-cell suspension were immediately immersed in tissue preservation solution (Miltenyi Biotec, Cat#130-100-008, Germany) and transported to the on-site laboratory in an ice bath for immediate processing. The fresh tissue samples were thoroughly rinsed three times with Dulbecco’s Phosphate-Buffered Saline (DPBS, Solarbio Cat# D1040, Beijing) that had been precooled to 4 °C. Subsequently, the samples were meticulously dissected into smaller fragments using surgical scissors and transferred into 1.5 mL centrifuge tubes for further processing. The tissue samples were then placed in a shaker and incubated at 37 °C for 30 to 50 min, utilizing the MACS Human Tumor Dissociation Kit (Miltenyi Biotec, Cat# DS130-095-929, Germany) to facilitate dissociation. The incubation process was halted once the digestive fluid became cloudy, signifying the complete disintegration of the tissue block. The resulting cell suspension was filtered through a 40-µm cell strainer and subsequently centrifuged at 400 × g for 5 min at 4 °C. Following the centrifugation, the supernatant was discarded, and the cells were resuspended in 1 mL of DPBS. Additionally, 3 ml of precooled red blood cell lysis buffer (Solarbio, Cat# R1010, China) was added. The cells were gently aspirated to ensure uniformity and incubated at 4 °C for 5 to 10 min, followed by another round of centrifugation. Subsequently, the cells were stained with a 7-Aminoactinomycin D (7-AAD) (eBiocience, Cat# 00-6993-50, United States) staining solution, prepared by mixing 100 μL of 1% BSA/PBS with 5 μL of 7-AAD. The staining procedure was carried out at 25 °C for 5 min. Following the staining procedure, high-quality individual cells were sorted using fluorescence-activated cell sorting (FACS) with a BD Aira II instrument. To verify the quality of the sorted cells, they were stained with trypan blue and examined microscopically prior to constructing the single-cell transcriptomic library.

Multi-labeled immunofluorescence staining and multispectral imaging
For multiplex immunofluorescence staining, we employed the PANO 7-plex IHC kit (Panovue, Cat# 10004100100, China). The slides were incubated in a 65 °C oven overnight to achieve deparaffinization. Subsequently, the tissues underwent a series of treatments with xylene, ethanol, and distilled water sequentially. Afterward, the slides were microwaved in the presence of an antigen retrieval solution (consisting of citric acid solution with pH 6.0/pH 9.0, Panovue, Cat #10020020500, and 10019020500, respectively, China). This was followed by sequential incubations with various primary antibodies and horseradish peroxidase-labeled secondary antibodies (100 μL per tissue). Subsequently, the slides underwent tyramide signal amplification (TSA) using the Fluorescence Kit (Panovue, Cat# 10021001010, China). After each incubation step, the slides were thoroughly washed with 1× PBST. Following each round of TSA, the slides were microwaved to enhance the staining process. Once all the antigens were labeled, the nuclei were counterstained with 4′,6′-diamidino-2-phenylindole (DAPI, Sigma-Aldrich, Cat# D9542, United States). Finally, the stained slides were scanned using a Mantra system (PerkinElmer, Waltham, Massachusetts, United States) to capture multispectral images. The detailed category of antibodies used in this study is provided in Supplementary Table 1. Cell density = number of positive cells/tissue area (mm2), TLS density = number of TLS/tissue area (mm2).

Single-cell RNA and TCR sequencing
The scRNA-seq and TCR libraries were generated using a 10X Genomics Chromium Controller Instrument and a Chromium Single Cell 5′ library & gel bead kit (10X Genomics, Cat# PN-1000165, United States). Briefly, cells were concentrated to approximately 1000 cells/μL and loaded into each channel to generate single-cell gel bead-in-emulsions (GEMs). Following the reverse transcription step, the GEMs were broken, and barcoded cDNA was purified and amplified. The amplified barcoded cDNA was then fragmented, A-tailed, ligated with adapters, and subjected to index PCR amplification. The final libraries were quantified utilizing the Qubit High Sensitivity DNA Assay (Thermo Fisher Scientific, Cat# Q32851 QUBIT, United States), and their size distribution was determined using a High Sensitivity DNA chip on a Bioanalyzer 2200 (Agilent). All the libraries were sequenced using an Illumina sequencer (Illumina, San Diego, CA) with a 150-bp paired-end run. Single-cell RNA-seq and TCR libraries were prepared on the Chromium platform (10X Genomics) in Berry Genomics (China), using the Chromium Next GEM Single Cell 5′ Kit v2 (10X Genomics, Cat# PN-1000263, United States). The FACS sorting live cells (7-AAD negative) were pooled together and washed trice with RPMI-1640, concentrated to 600–1000 cells per microliter, and immediately loaded on the 10X microfluidic chip (10X Genomics) following the manufacturer’s protocol to generate a cDNA library. The amplified cDNA was subsequently utilized for constructing a 50-gene expression library construction and TCR V(D)J targeted enrichment was performed with the Chromium Single Cell V(D)J Enrichment Kit, Human TLC (10X Genomics, Cat# PN-1000005, United States), followed by V(D)J library construction.

Spatial transcriptomics profiling
Fresh GC tissues were collected intraoperatively and processed into 5 μm-thick formalin-fixed, paraffin-embedded sections. RNA quality was assessed using DV200, with samples meeting a minimum threshold of ≥30% selected for downstream analysis. Sections were mounted on Superfrost™ Plus microscope slides (Fisherbrand™), deparaffinized, and stained with H&E. Imaging and decrosslinking were performed following the manufacturer’s protocol (CG000520; 10x Genomics, United States). The Human Probe V2 panel, containing gene-specific probe pairs, was hybridized to the stained sections at 50 °C overnight. Post-hybridization washes and ligation were carried out to join adjacent probe pairs targeting the same transcript. A probe release mix containing RNase and tissue-dissociation enzymes was used to liberate ligated products, which were subsequently captured on Visium CytAssist Spatial Gene Expression Slides using the Visium CytAssist instrument (CG000495; 10x Genomics, United States). Captured probes were extended with Unique Molecular Identifiers (UMIs), spatial barcodes, and partial read sequences. Released products were eluted with 0.08 M KOH and pre-amplified (10 cycles), followed by purification with 1.2× AMPure XP (BECKMAN COULTER, Cat# A63882) select beads. The optimal PCR cycle number for library amplification was determined, and indexed libraries were generated to prevent barcode overlap during multiplexed sequencing. Final libraries were sequenced on the Illumina NovaSeq 6000 platform using 150 bp paired-end reads in Berry Genomics (China).

Generation of ACKR1- overexpressing Endothelial cells (ACKR1OE EC)
The human ACKR1 cDNA sequence was obtained from the pCMV-HA-FLAG-ACKR1 (Neo) plasmid purchased from Miaolingbio (China). The ACKR1 fragment was amplified by PCR using the following primers (Forward primer: ACKR1-F 5′-ATGGGGAACTGTCTGCACAG-3′, Reverse primer: ACKR1-R 5′-GGATTTGCTTCCAAGGGTGT-3′). The PCR product was cloned into the PiggyBac transposon system vector to generate the pB-CAG-ACKR1-EGFP-BSD plasmid. The correct construction of the plasmid was confirmed by sequencing. Subsequently, the pB-CAG-ACKR1-EGFP-BSD plasmid was transfected into human induced pluripotent stem cells (hiPSC) using Lipofectamine Stem reagent (Thermo Fisher, Cat# STEM00003, United States). Forty-eight hours post-transfection, 2 μg/mL Blasticidin S (Beyotime, Cat# ST018, China) was used to select stable ACKR1-overexpressing cell lines for 7 days. The differentiation protocol for generating ACKR1OE EC from human hiPSCs was based on previously published methods44. Briefly, hPSCs were dissociated into single cells and seeded at a density of 100,000 cells per well in a 12-well plate, followed by recovery and adhesion in mTeSR medium for 24 h. For mid-primitive streak induction, cells were cultured in CDM2 medium supplemented with 30 ng/mL Activin A (PeproTech, Cat# 120-14P-500UG, United States), 40 ng/mL BMP4 (PeproTech, Cat# 120-05-5UG, United States), 6 μM CHIR99021 (Selleck, Cat# S2924, United States), 20 ng/mL FGF2 (PeproTech, Cat# 100-18B-50UG, United States) for 24 h. The resulting primitive streak cells were then cultured in CDM2 medium containing 40 ng/mL BMP4, 2.5 μM GDC-0941 (Cellagen Technology, Cat# 957054-30-7, United States), 10 μM Forskolin (Selleck, Cat#S2449, United States), 2 μM SB-505124 (Selleck, Cat#S2186, United States), 100 ng/mL VEGF (PeproTech, Cat# 100-20-10UG, United States), 1 μM XAV939 (Selleck, Cat# S1180, Selleck,), and 200 ng/mL ascorbic acid-2-phosphate (AA2P) for an additional 24 h to induce lateral mesoderm. These cells were further differentiated in CDM2 medium supplemented with 2 μM SB-505124, 250 nM DMH1 (Selleck, Cat# S7146, United States), 2 μM RO4929097 (Selleck, Cat# S1575, United States), 100 ng/mL VEGF, 1 μM XAV939, and 200 ng/mL AA2P for 24 h. Finally, terminal differentiation was carried out in CDM2 medium containing 2 μM SB-505124, 2 μM RO4929097, 500 nM PD0325901 (Selleck, Cat# S1036, United States), 1 μM CHIR99021, and 200 μg/mL AA2P. The efficiency of ACKR1OE EC generation was assessed by fluorescence-based analysis.

Quantitative polymerase chain reaction (qPCR)
ACKR1OE ECs were treated for 24 h with either PBS (control), 50 ng/mL CXCL13, 5 ng/mL TNF, or a combination of both cytokines (50 ng/mL CXCL13 + 5 ng/mL TNF). Total RNA was extracted using TRIzol reagent (Invitrogen, Cat# 15596026, United States) and reverse-transcribed into cDNA using the HiScript III All-in-One RT SuperMix Perfect for qPCR kit (Vazyme, Cat# R333, China). Quantitative PCR was performed with the miScript SYBR-Green PCR Master Mix on an ABI Prism 7900 system. Gene expression levels were normalized to GAPDH, and all assays were conducted in accordance with the Minimum Information for Publication of Quantitative Real-Time PCR Experiments (MIQE) guidelines. Primer sequences are provided in Supplementary Table 2.

TLC recruitment in the TLC/ ACKR1OE EC transwell system
To assess the effect of ACKR1OE EC on the TLC recruitment via the transwell co-culture system, 70% confluent ACKR1OE ECs were treated with either PBS (control), 50 ng/mL CXCL13, 5 ng/mL TNF, or a combination of both cytokines (50 ng/mL CXCL13 + 5 ng/mL TNF). At the same time, PBMCs were isolated by lymphocyte separation medium (Solarbio, Cat# P8610, China), TLCs were enriched by CD3/28 separation and activation magnetic beads (Seafrom, GMP-TL603-1000, China), and were cultured in the VIVO-15 medium containing 200 IU/ mL IL-2. Then, 1 × 105 TLCs were placed in the 5μm upper inserts of a transwell system (LABSELECT, 02225110E). After 24 h of culture, the inserts were transferred into fresh wells of ultra-low attachment plates supplemented with VIVO-15 medium, in which the chemokine-treated ACKR1OE EC were seeded in the bottom chamber of a 24-well ultra-low attachment plate. TLCs were counted under a stereomicroscope using ImageJ after 24 h.

scRNA-seq data processing
The raw sequencing data were aligned to the GRH38 reference genome using cellranger (10X Genomics, v5) count and vdj functions. The resulting count matrices of gene expression for each sample were imported into Seurat (v4.0.2)45. High-quality cells were selected for further analysis based on three criteria: (1) cells contained more than 2001 UMIs, fewer than 6000 and more than 301 expressed genes, and fewer than 10% of UMIs derived from the mitochondrial genome; (2) genes were expressed in more than ten cells within a sample; (3) cell doublets were identified and removed using the DoubletFinder46 package (v2.0.3). The cells-by-genes expression matrices of the remaining high-quality cells were integrated using the RunFastMNN function provided by the SeuratWrappers R package (v0.3.0) and then normalized to the total cellular UMI count. Then, gene expression matrices were scaled by regressing out the total cellular UMI counts and the percentage of mitochondrial genes. The principal component analysis (PCA) was conducted using highly variable genes, and the top 30 significant principal components were selected for subsequent performance of Uniform Manifold Approximation and Projection (UMAP) to achieve dimensionality reduction and visualization of gene expression. The FindAllMarkers function in Seurat was employed to identify the top DEGs between two groups, with the min.pct parameter set at 0.2, thereby considering only those genes expressed in more than 20% of cells. The nonparametric Wilcoxon rank-sum test was utilized to compute the p values for comparisons, and the adjusted p value based on Bonferroni correction was calculated. Cell subclusters exhibiting similar gene expression patterns were annotated as the same cell type, and the cell types in the resulting two-dimensional representation were further annotated to cell types and subtypes using canonical marker genes and the top 30 DEGs.

TCR analysis
We employed the scRepertoire47 R package (v1.7.2) to examine the TCR clonality of TLC. Briefly, the filtered_contig_annotations output from Cell Ranger was imported into Seurat to generate a list object of TCR genes and CDR3 sequences by cell barcodes using the combineTCR function. Subsequently, the combineExpression function was utilized to integrate the clonotype information with the filtered Seurat object. Clonotype frequencies were categorized by patient using the following parameters: single = 1e-4, small = 0.001, medium = 0.01, large = 0.1, and hyperexpanded = 1. Clonotypes were determined using the VDJC genes comprising the TCR. The distribution of clonotype bins was visualized using the Seurat DimPlot function.

Cell‒cell communication analysis at the single-cell level and spatial level
We utilized the R package CellChat48 (v1.6.1) to analyze cell-to-cell communication between tumor cells and other cell types. Initially, a CellChat object was established by grouping defined clusters. The “CellChatDB.human” ligand‒receptor interaction database was used for analysis without additional supplementation, and all preprocessing steps were performed using default parameters. The functions computeCommunProb and computeCommunProbPathway were employed to infer the network of each ligand‒receptor pair and signaling pathway individually. Visualization was achieved through a hierarchy plot, circle plot, and heatmap, each serving as a distinct form of visualization.

Single-cell copy number analysis
Copy number instability was assessed with the R package infercnv49 (v1.8.1), which is specifically designed to infer copy number alterations from tumor single-cell RNA-seq data. This package compares the expression intensities of genes across malignant cells in advanced and recurrent tumor tissues. Epithelial cells in normal tissues and the immune cells were used as reference21,23.

ST-seq data analysis, spatial cell type, and TLS abundance estimation
Raw ST-seq reads were quality-checked and aligned using Space Ranger (v1.2.0). Normalization across spatial spots was conducted using the LogVMR function. Following quality control by Space Ranger software, dimensionality reduction and clustering were performed via PCA, using the top 30 components at a resolution of 0.8. Spatial feature expression patterns were visualized with Seurat’s spatial feature plot function. To infer spatial cell-type and TLS composition or abundance, we first obtained the cell-type-specific signature, using the top 30 DEGs, which were identified using the FindAllMarkers function provided by Seurat. The traditional TLS-specific signature genes included CD3D, CD3E, CD19, MASA1, and CD79A, the ssTLS signature genes included CD3D, CD3E, CD19, MASA1, CD79A, CXCL13, CXCR5, CXCL16, CCL19, CCL22, ACKR1, and SELE. The GSVA50 was then employed to calculate the cell type or TLS-specific signatures score (cell-type abundance) of each spatial spot, which was visualized by the SpatialFeaturePlot function provided by the Seurat package.

Pathway enrichment
To illustrate the enriched signaling pathways in BLC and myeloid subtypes, the GSVA50 package (v1.38.2) was employed to assess pathway differences using the C2 curated gene set provided by the Molecular Signatures Database (MSigDB), which was calculated with a linear model offered by the limma package (v3.46.0).

Immunotherapy response and survival analysis
GC immunotherapy bulk RNA-seq data, along with curated clinical data from patients, were obtained from a previous study51. To evaluate the association between the signatures of immune and stromal subtypes identified in our study and the anti-PD1 immunotherapy response to GC patients’ survival, the GSVA50 (v1.38.2) was used to calculate the cell-type-specific signatures (top 30 marker genes) score. Patients were classified into high- or low-gene expression groups based on the median expression level of the cell-type signature score. Then, the Kaplan–Meier survival curves were plotted using the R package “survminer” (v0.4.9). Subsequently, a Cox proportional hazards model was conducted, incorporating age and tumor stage as covariates. The results of the Cox regression model among the various cell subtypes were visualized using weighted Z-scores.

Statistics and reproducibility
No statistical methods were used to predetermine the sample size; all available clinical GC samples were included in the study. No data were excluded from the analyses. Detailed statistical approaches are described in the respective figure legends. Sex-based differences were not considered, as they were not relevant to the study objectives. Data are presented as mean ± SEM, with sample numbers indicated. Group comparisons were assessed using the Wilcoxon rank-sum test, the Kruskal–Wallis test, followed by Dunn’s multiple comparison tests were used for three or more group comparisons, and survival analyses were performed using the two-sided log-rank test. All analyses were conducted in R (version 4.0.2; R Foundation for Statistical Computing). The following R packages were used: Seurat (for single-cell data processing), survival and survminer (for survival analysis), pwr (for the sample size determination of patients to assess the role of TLS status on the prognosis of GC), and ggplot2 (for data visualization). A p value < 0.05 was considered statistically significant.

Reporting summary
Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article.

Supplementary information