A spatially resolved atlas of gastric cancer characterises a lymphocyte-aggregated region

Domain: neuroscience
Published: 2026-01-27
Editorial outcome: minor_revision

Abstract
The tumour microenvironment is a focal point in cancer immunotherapy: its cellular composition and spatial organisation can affect the clinical outcomes of cancer patients. By integrating single-cell and spatial transcriptomics, we identify four spatial regions and survey how cellular spatial distribution varies across them in gastric cancer. One region, the Lymphocyte Aggregated Region, consists of lymphocyte aggregates and tertiary lymphoid structures. Within it, we observe associations between naive T cell abundance and T cell activation-associated pathways, and correlations exist between distribution patterns of different lymphocytes and two transcriptomically distinct groups—more activated lymphocytes reside in the adjacent cancerous regions of Group A, while more resting lymphocytes settle in those of Group B. Within Group A, PD1+CD27+CD8 T cells cluster in closer proximity to CD70+LAMP3+ dendritic cells. Our study unveils the gastric cancer tumour microenvironment at a spatial resolution and provides insights into the exploration of immunotherapy biomarkers.

Introduction
Deciphering the tumour microenvironment (TME) in detail is challenging due to its complex cellular makeup, including immune and non-immune cells, and especially the complicated and varied cellular spatial organisation within it. Previous studies have profiled the cellular composition of the immune compartment in multiple cancers by single-cell RNA sequencing (scRNA-seq), but the spatial distribution pattern of lymphocytes, particularly how critical tissue architecture like tertiary lymphoid structures (TLSs) affects their distribution, remains ambiguous1. As a hot hub of lymphocytes, the presence and density of TLSs are related to improved prognostics and survival predictions in various cancers2. Recently, growing evidence has revealed that the maturation status of TLSs ranges from mere lymphocyte aggregates to highly organised structures resembling lymphoid nodes3,4, obscuring the exact roles of TLSs over cancer progression. Hence, it is imperative to systematically depict the influence TLSs exert on the TME or vice versa, the anti-tumorigenic immune responses orchestrated by TLSs, and novel makers specifically expressed within TLSs that could potentially govern the anti-tumorigenic activities.

Furthermore, cell-cell interactions and the locations where they occur often modulate biological activities within the TME. Being mediated by ligands and receptors, these cellular communications can trigger anti-tumorigenic activities ranging from naive T cell priming to tumour cell killing by cytotoxic T cells and pro-tumour behaviours like immune evasion exploited by tumour cells. Notably, critical ligand-receptor pairs utilised by immune checkpoint blockade (ICB) therapy, such as CTLA4-CD86 and PD1-PDL1, can coordinate anti-tumorigenic immunity and thus shape the TME5–8. However, the expression variation of immune checkpoint molecules across different tumours cannot fully explain why the treatment response to these molecules varies considerably among patients. One hypothesis is that, aside from expression variation, spatial locations at which these molecules are expressed could influence the subsequent anti-tumorigenic activities9. However, prior studies employing scRNA-seq data to uncover cellular communications in the TME relied mainly on in silico inference based on gene expression rather than in situ exploration based on spatially resolved co-occurrence of ligands and receptors10. Thus, integrating scRNA-seq with spatial transcriptomics (ST) that can capture gene expression across distinct tissue architectures is a suitable approach to investigating cellular interactions in a spatial context. Moreover, despite completing T cell atlases in multiple cancers separately and jointly11–22, we have yet to thoroughly explore the spatial distribution of T cell subsets in the TME, especially for exhausted T cells that are pivotal indicators in prognostics and clinical responses to immunotherapy23–25. Utilising the ST platform26, we are now positioned to examine the spatial distribution pattern of these lymphocytes in detail.

In this work, we construct a fine-grained single-cell atlas from matched samples and use it as a reference to deduce the cellular proportion in each ST spot. By integrating the scRNA-seq and ST data of gastric cancer (GC), our study may pave the way to discover potential immunotherapeutic markers by unveiling the spatial organisation of cell subtypes, the fundamental roles of TLSs in coordinating anti-tumorigenic activities in the TME, and the spatial distribution of immune checkpoint molecules.

Results
Construction of a fine-grained single-cell atlas for GC
We first built a high-resolution single-cell transcriptomic atlas by sequencing cells derived from four different tissue sources, including surgically resected tumour (denoted as T), histologically normal adjacent tissue (N), lymphoid node (L), and peripheral blood (P), from a cohort of 27 treatment-naive donors with GC, covering two clinical phenotypes — subtype (diffuse, indeterminate, and intestinal) and stage (I, II, and III) (Fig. 1a, b, Supplementary Fig. 1a, and Supplementary Data 1, 2). In addition, 19 tumour samples were assayed by 33 ST sections, covering three categories of sampling locations in the tumour (core, in-between, and edge) (Supplementary Fig. 1b). Overall, we obtained transcriptomes for 470,609 single cells and 61,035 ST spots after stringent quality control.

Eight major cell types spanning from immune cells (T cells, B cells, innate lymphoid cells, mast cells, and myeloid cells) to non-immune cells (endothelial cells, stromal cells, and epithelial cells) were identified (Fig. 1c), each validated by canonical gene markers (Supplementary Fig. 1c). The platform-induced batch effect (10×3’ versus 5’) was efficiently removed as the proportion of major cell types was comparable between these two platforms across all tissues sampled (Supplementary Fig. 1d). Each major cell type was identified across all tumour samples (Supplementary Fig. 1e), while the bona fide biological differences were preserved as each major cell type had a distinct gene expression pattern (Supplementary Fig. 1f and Supplementary Data 3). We further clustered the major cell types into a total of 53 subtypes, non e was exclusively derived from or dominant in a single tumour sample (Supplementary Fig. 1g, h). Each cell subtype had a unique tissue distribution preference (Supplementary Fig. 1i).

Since immune populations are a prime player in most anti-tumorigenic activities in the TME, we next shifted our focus from the global cellular composition of GC to its immune compartment, with transcriptomic heterogeneity within T cells (CD8 and CD4 T cells) and myeloid cells observed (Supplementary Fig. 2a–d and Supplementary Data 4). Critical immune subtypes in the TME, including CD8-Tex-LAYN, CD4-Treg-FOXP3, and cDC-LAMP3, were abundant in the tumour, in line with previous research17,19 (Supplementary Fig. 1i). In addition, by scoring each T cell subset with gene programs derived from the TCR and proliferation signalling pathways separately19,24, we identified CD8-Tex-LAYN as potential tumour-reactive T cells that could participate in local immune responses in GC (Supplementary Fig. 2f). Briefly, we generated a fine-grained transcriptomic atlas composing a diverse group of cell subsets from four different tissues of GC donors, particularly for the immune compartment.

Enrichment preference of immune cells across tissues and clinical phenotypes
Since immune cells can migrate between different types of tissue, we next compared the tissue enrichment preference for each immune cell subset of T cells and myeloid cells in GC. Several cell subsets, including CD4-Tcm-ANXA1-P, CD4-Tn-CCR7-P, CD8-Temra-CX3CR1-P, CD8-MAIT-SLC4A10-P, Monocyte-CD14, and Monocyte-CD16, were enriched in the blood but were hardly detectable in solid tissues, suggesting a distinct immune cellular composition in the blood compared with other tissues (Supplementary Fig. 3a). Indeed, these cell subsets constituted the vast majority of immune cells in the blood. T cell subsets like CD8-Temra-CX3CR1-P and CD8-MAIT-SLC4A10-P, enriched in the blood had counterparts enriched in the tumour, respectively (Supplementary Fig. 1i, 2g, and 3a). Additionally, for both T cells and myeloid cells, the compositional diversity measured by the Shannon index27 was significantly higher in tumours than that in adjacent normal tissues and blood (P<0.05), but no significant difference was detected between tumours and tumour-draining lymph nodes (TDLNs) (Supplementary Fig. 3b). A similar pattern was observed for myeloid subsets, but a significant difference was detected between tumours and TDLNs. Our results indicated that, despite a different cellular composition observed, the diversity of T cells in GC tumours was roughly equivalent to that in TDLNs.

Next, we investigated the cellular distribution preferences within each clinical phenotype, including subtype and stage. Unlike the findings in different types of tissue, few significant differences were identified across different subtypes and stages, except that when accounting for the number of patients for each GC group, the diversity of CD8 T cells in the intestinal subgroup was higher than that in subgroups of diffuse and indeterminate (P<0.05), with the diffuse subgroup exhibiting the lowest diversity score (Supplementary Fig. 3b–d). We then performed the same analysis solely for tumour samples: when accounting for the number of patients for each GC group, the diversity of CD8 T cells in the diffuse subgroup was not only significantly lower than that in the intestinal subgroup but also the indeterminate subgroup; minimal differences were observed between different stages across T cell and myeloid subtypes (Supplementary Fig. 3e–h). Such variation in the diversity of CD8 T cells across GC subtypes could be associated with a lower survival rate in the diffuse subgroup28–30, but more supporting data and analyses should be appended in the future.

Construction of a spatially informed cellular atlas to reveal cell distribution
We noticed that the diversity of histological structures observed on an H&E (hematoxylin and eosin) image of an ST slide varied drastically: some were dominated by the cancerous area, like D07T2; some had TLSs mixed within cancerous regions, such as D21T2. The rest have a much more complex pattern; for example, smooth muscle and stromal region, TLSs, and cancerous region were observed all together in D02T2 (Supplementary Fig. 4a). We suspected that if a single histological pattern dominated a slide, ST spots from that slide could have similar transcriptomic profiles, which can be directly measured by the purity score—a higher purity score indicates more similarities were observed in the transcriptomic profiles of all spots within an ST section31. Indeed, we observed the highest purity score in D07T2, followed by D21T2 and D02T2. Furthermore, likely due to the invasive characteristics of cancer cells, each histological structure was tangled and mixed with others without clear boundaries, complicating the cellular composition of each tissue architecture. Thus, we categorized the 33 ST datasets into three groups based on the purity score generated from the above step and determined the resolution score for each group. After decomposing the TME of GC into various cell subtypes, we mapped them (excluding 11 cell subtypes specifically enriched in the peripheral blood) onto the ST slide to characterise their spatial distribution and explore their association with tissue architectures at the transcriptomic level (Supplementary Fig. 4b).

By integrating scRNA-seq and ST data using RCTD32, we clustered all ST spots into four spatial regions, each possessing a distinctive group of cell subtypes: one exhibited intense signals of CD45-positive cells, with the majority being lymphocytes, such as CD8-Tn-CCR7 and Bn-IGHD; one displayed strong signal of SmoothMuscle-MYH11; the rest were populated with Cancer-CEACAM6, with one exhibiting lower Cancer-CEACAM6 signals but higher immune cell signals, whereas the other displaying the opposite trends (Fig. 1d and Supplementary Data 5). Furthermore, we noticed that the distribution abundance of each cell subtype varied in two contrasting tendencies—mostly lymphocyte subtypes exhibited smaller deviations while non-lymphocyte cell subtypes demonstrated larger ones (Fig. 1d). Since the RCTD enrichment scores, to some degree, can reflect the cellular proportion regardless of the tissue morphology, smaller deviations observed in lymphocyte abundance across samples are probably associated with the migration characteristics of T cells and B cells even in solid tissue like tumours. Next, to validate the clustering results, we overlaid these ST groups onto H&E images and found that these groups mirrored the major histological patterns confirmed by pathological physicians — lymphoid regions that comprise TLSs, smooth muscle and stromal region, and cancer regions (Supplementary Fig. 4c). Such alignment between the ST group and the histological pattern suggested the robust integrative approach adopted in our study. Therefore, we termed these groups derived from ST data “Lymphocyte Aggregated Region” (abbreviated as LAR), “Smooth Muscle and Stromal Region” (SMSR), “Immunogenic Cancer Region” (ICR), and “Negative-immunogenic Cancer Region” (NCR), respectively.

Each ST group was identified either in all samples (ICR and NCR) or most samples (LAR and SMSR), composed of spots distributed evenly across all samples, and possessed a unique transcriptional profile (Fig. 1e–g and Supplementary Data 6)—the expression of genes associated with lymphoid nodes (CXCL13, CCL19, CCL21, and FDCSP) was significantly higher in LAR, which occupied the least proportion of spots (Fig. 1f); smooth muscle associate genes (ACTA2 and ACTG2) were up-regulated in SMSR; and the expression of GC-associated genes, which were identified from bulk samples of GC tumours using GEPIA33, was significantly higher in NCR than ICR (Fig. 1e and Supplementary Fig. 4e).

Survival analysis using the TCGA (The Cancer Genome Atlas) database indicated that higher expression of LAR-enriched genes was associated with a better prognosis in multiple cancers, notably in sarcoma (SARC) and skin cutaneous melanoma (SKCM) (Supplementary Fig. 4f–h). In contrast, higher expression of SMSR-enriched, ICR-enriched, and NCR-enriched genes was associated with worse prognosis in multiple cancers (Supplementary Fig. 4g). By performing a meta-analysis on 13 ICB therapy cohorts, we observed a modest association between the higher expression of LAR-enriched genes and response to anti-PD-1 or anti-PD-L1 treatment (Supplementary Fig. 4i and Supplementary Data 7). Taken together, by combining the scRNA-seq and ST data, we constructed a spatially informed cellular atlas of GC whereby the distribution pattern of each cell subtype across numerous tumours was depicted and pinpointed the potential anti-tumorigenic roles of LAR in the TME.

The complexity and diversity of the TME are revealed in a spatial context
Next, we examined whether the spatially resolved TME constitution, defined as the four ST groups (LAR, SMSR, ICR, and NCR), correlates with clinical phenotypes such as subtypes and stages. We spotted that the share of LAR escalated across subtypes, following the order of diffuse, indeterminate, and intestinal. We further quantified these trends using the ratio of observed to expected spot numbers (Ro/e) and odds ratio (OR) (Supplementary Fig. 5a–c). Particularly, the proportion of LAR increased from early-stage to late-stage tumours. However, when comparing the ratio of ICR versus NCR across different stage groups, we observed insignificant differences among groups (P<0.05, two-sided Wilcoxon test), largely due to the small sample size for each group. In addition, we examined if any enrichment preference for the ST groups exists in sub-categories of the above phenotypes by measuring the diversity index and found that the diversity index was roughly equal at both the sample and spot levels (Supplementary Fig. 5d, e).

By measuring effect size, we prioritised the abundance of all cell subtypes within each ST region. Weighted by stringent cutoff values (>0.8 and <-0.8), the LAR region was strongly enriched with immune cells, including GCB-LRMP, cDC-LAMP3, and CD4-Th17-IL17A but depleted with Cancer-CEACAM6 (Fig. 1h). We also observed a positive correlation between signals of the above immune cell types and that of LAR and a negative association between signals of cancer cells and that of LAR in GC tumours collected by TCGA (Supplementary Fig. 5f, g). Apart from Smooth Muscle-MYH11, SMSR was enriched with two types of fibroblasts — MyoFibroblast-ACTG2 and Fibroblast-PI16, and was depleted with Cancer-CEACAM6 using medium cutoff values (>0.5 and <-0.5) (Supplementary Fig. 5f).

Despite displaying stronger cancer cell signals than LAR and SMSR, ICR was mainly occupied by Mucous-MUC5AC, a cell subtype secreting mucus to protect the gastric mucosa from acid erosion34, indicating that a large proportion of ICR somehow still possessed normal stomach functions. Rather, Cancer-CEACAM was the dominant subtype in NCR, suggesting a terminal cancerous stage within this region. Compared with LAR and SMSR, a collection of cell subtypes was depleted in ICR and NCR when using relaxed thresholds (>0.2 and <-0.2). Quiescent states of immune cells such as Bm-GPR183, CD8-Tn-CCR7, and CD8-Tn-CCR7 were depleted in ICR, whereas functional states of immune cells like cDC-LAMP3, GCB-LRMP, CD4-Treg-FOXP3, and Plasma B cells were depleted in NCR (Supplementary Fig. 5h). In conclusion, as opposed to the cellular composition of LAR and SMSR, a transitional, complex cellular composition existed in ICR and NCR, indicated by a larger number of depleted cell subtypes and the less stringent cutoff values used to prioritize cell types in these spatial regions.

Next, we explored the association between the expression of lymph node-associated genes (CXCL13, CCL19, and FDCSP) and sub-categories within each GC tumour phenotype (subtype, stage, and location) by measuring the effect size (Fig. 1i and Supplementary Fig. 5i, j). These gene signatures are actively involved in lymphocyte recruitment and activation2, and higher expression of these genes was associated with improved immunotherapy responses and better survival in multiple cancer types4,35,36; most importantly, in our study, expression of these genes was higher in LAR than other regions, in line with previous research. Despite the importance of these genes, the correlation between the expression of these genes and the sub-categories for each clinical phenotype of gastric cancer is vague. CCL19 and FDCSP were strongly associated with LAR versus other regions across all spots. When categorised by stage, the association was only significant in stage III (cutoff > 1). Moreover, the association between FDCSP and LAR was substantial at the tumour edge, whereas the association between CCL19 and LAR was significant in the diffuse subtype. Conversely, the association between CXCL13 and LAR was insignificant across all spots, but when evaluated by each phenotype, a significant association was found at the tumour core. Next, we assessed the association between these genes and two other genes associated with lymphoid nodes (CCL21 and CR2) within LAR at different tumour locations. We observed that, from the tumour core to the edge, the association between these genes was increasing (Supplementary Fig. 5k), in concordance with the association of LAR with FDCSP and CCL19 at the tumour edge.

LAR associated with T cell-mediated immune responses
Next, we charted the biological activities in which each ST region participated by measuring the intensity of signalling pathways (archived in Reactome37). Distinctive signalling pathways possibly reflecting the function of each ST region were identified, for instance, “Adaptive and Innate Immune System” in LAR, “Smooth Muscle Contraction” in SMSR, and “Cellular Response to Hypoxia” in NCR (Fig. 1j). Particularly, the “TCR Signalling” pathway, which implicates T cell activation for both naive and effector T cells38, was strongly expressed in both LAR and NCR. Furthermore, two other signalling cascades that were critical in coordinating T cell activation, including “Costimulation By The CD28 Family” and “Antigen Processing Cross Presentation”, exhibited intense signals in LAR (Fig. 2a). We further validated this finding at the level of spatial block: an individual spatial block consists of a varying number of spatially grouped spots labelled with the same ST group (Fig. 2b and Supplementary Fig. 6a). Compared with other spatial regions, LAR possessed significantly higher signals of the above pathways and a higher proportion of blocks exhibiting these signals (Fig. 2c).

In addition, the signal of “Antigen Processing Cross Presentation” was also higher in NCR. In NCR, antigen-presenting cells, such as cDC-LAMP3 and Macrophage-C1QC, were depleted (Supplementary Fig. 5h). Other antigen-presenting cells, including but not limited to cDC1-CLEC9A, cDC2-CD1C, and Macrophage-VCAN, were not depleted in NCR, a region mostly dominated by Cancer-CEACAM6. Thus, we speculated that the higher signal of “Antigen Processing Cross Presentation” in NCR could be due to a combined effect of antigen processing activities of these cells.

Next, we explored the association between the abundance of each T cell subtype and signalling pathways, including TCR and Proliferation signalling pathways24 within LAR. CD8-Tn-CCR7 and CD4-Tem-GZMK were the only T cell subtypes with significant associations between their abundance and TCR and Proliferation signalling pathways within LAR (Fig. 2d–e). The same was true for the association between cellular abundance and Costimulation By The CD28 Family (Supplementary Fig. 6b). In contrast, within the non-LAR region, negative associations were observed between the cellular abundance of CD8-Tn-CCR7 and TCR signalling (correlation coefficient: −0.031, P=0.06), Costimulation By The CD28 Family (−0.029, P=0.07), and Antigen Processing Cross Presentation (−0.059, P=0.0002). An insignificant association was observed between the cellular abundance of CD8-Tn-CCR7 and Proliferation signalling pathways within the non-LAR region (0.022, P=0.17). Since these signalling pathways are prerequisite biological processes for T cell activation38, our results indicated that these cell subtypes could be potentially primed within LAR. In addition, we identified two naive forms of T cells, CD8-Tn-CCR7 and CD4-Tn-CCR7, but only the abundance of CD8-Tn-CCR7 was associated with the above signalling pathways in LAR, whereas significant associations were not identified for CD4-Tn-CCR7.

We next assayed the expression of CD8, TCF1, and KI67 in 11 tumour samples from the main study cohort and five samples from an independent validation cohort using multiplex immunohistochemistry (mIHC). LAR-residing CD8+TCF1+KI67+ cells were observed in various samples (Fig. 2f and Supplementary Fig. 6c, d). Within LAR, we counted the number of CD8+ cells and the number of TCF1+KI67+CD8+ cells and found an average of 12.19% (SD=0.0839) CD8+ cells were TCF1+KI67+ cells from samples in the main study cohort, while a lower average of 0.77% (SD=0.007) CD8+ cells were TCF1+KI67+ cells from samples in the independent validation cohort (Supplementary Data 8). In addition, we assayed 39 replicates from the main study cohort using ST (Supplementary Data 1), collected 280 public ST data across multiple cancers (Supplementary Data 9), and annotated them by SciBet39 using the four GC ST groups as reference. Significant associations between CD8A/CD8B, TCF7, and MKI67 were observed across LAR spatial blocks of the replicate samples (Supplementary Fig. 6e). LAR spatial blocks simultaneously expressing CD8A/CD8B, TCF7, and MKI67 were detected in a group of samples across cancer types (Supplementary Fig. 6f).

Cellular heterogeneity observed in LAR and its association with lymphocyte distribution
We observed that the expression of key gene features (FDCSP and CCL19), which was measured by the coefficient of variation (CV)—a higher CV (>1) indicates greater relative variability, varied drastically among the LAR blocks (a CV value of 3.41 for FDCSP and 1.28 for CCL19), and a significant association was observed between these two trends (Fig. 3a). The abundance of cell subtypes (GCB-LRMP and cDC-LAMP3) also varied drastically among the LAR blocks (a CV value of 2.19 for GCB-LRMP and 2.40 for cDC-LAMP3), but relatively low CV values (0.33 and 0.41) were observed for “TCR Signalling pathway” and “Costimulation By The CD28 Family.” Even for LAR blocks within a single sample, varied signals of GCB-LRMP (CV=0.51) were observed, as indicated by an overlay of cellular abundance onto H&E images (Fig. 3b).

We inferred that such a variation among LAR blocks mirrored distinct groups of LARs. To explore the interplay between varied LARs and diverse types of TME, we aggregated LAR blocks by donors and dichotomised them based on their cellular abundance (Fig. 3c). By applying this approach, even though a disparity among LARs existed on a transection of a tumour (Fig. 3b), we could systematically evaluate the correlation between LARs and the TME across donors. One group of donors, termed Group A, consisted of spatial blocks exhibiting intense signals of cell subtypes enriched in TDLNs, such as GCB-LRMP and CD4-Tfh-CXCR5, with the majority originating from donors of stage III; by contrast, LAR blocks from Group B had a weaker signal of the above cell subtypes. In addition, TDLN-associated signalling pathways, cell subtypes, and hallmark genes enriched in LAR were considerably higher in the LAR region of Group A than in Group B (Supplementary Fig. 7a–c), and such a trend was irrelevant to the size of LAR spatial blocks (Supplementary Fig. 7d, e). More transcriptomic similarities were observed between TDLNs and Group A blocks than between TDLNs and Group B blocks (Supplementary Fig. 7f). A higher expression of LAR-enriched genes, mostly also highly expressed in TDLNs, was observed in bulk samples of Group A than B, further substantiating that the gene expression pattern of Group A LARs was more similar to that of TDLNs than Group B LARs (Supplementary Fig. 7g).

Moreover, compared with other ST regions, ICR was more enriched in Group A than Group B, but NCR was comparable between Groups A and B (Supplementary Fig. 7h). Since ICR has higher immune cell signals but lower cancer cell signals than NCR (Fig. 1d), our results suggest a more deteriorated cancerous state in Group B samples. We then systematically surveyed the cellular abundance differences between Groups A and B within cancerous regions (ICR and NCR) rather than in ICR or NCR separately, as ICR and NCR were spatially tangled with each other. The number of cell subtypes with increased cellular abundance in Group A was nearly three times as many as that in Group B. Notably, among these cell subtypes, most activated lymphocytes — including CD8-Tex-LAYN and CD4-Treg-FOXP3 — were significantly elevated in the cancerous regions of Group A than Group B; on the contrary, most resting lymphocytes, such as CD4-Tn-CCR7 and CD8-Tn-CCR7, were increased dramatically in the cancerous region of Group B than Group A (Fig. 3d, e and Supplementary Fig. 7i, j). To verify the above trends, we implemented an artificial intelligence-guided algorithm, RedeHist, to perform nuclei segmentation on H&E images and impute gene expression for individual cells it detects; this approach allows us to validate scientific findings at the single-cell resolution. More CD8-Tex-LAYN and CD4-Treg-FOXP3 were detected in the cancerous regions of Group A than in Group B (Supplementary Fig. 7k). Taken together, our analyses suggested that distinct groups of LARs could affect the distribution of lymphocytes in the TME, as activated lymphocytes were aggregated in the adjacent cancerous regions of Group A LARs while resting lymphocytes were cumulated in those of Group B LARs.

Distinct spatial distribution pattern of immune checkpoint molecules
Immunotherapies targeting immune checkpoint molecules have advanced and evolved substantially over the last decade40,41, yet the spatial distribution pattern of these molecules in tumours remains obscure. Therefore, by scanning immune checkpoint molecules as ligand-receptor pairs, we delineated their spatial distribution across ST regions. In our analysis, we included immune checkpoint ligand-receptor pairs that have been successfully targeted in immunotherapies for cancer treatment, despite varied response rates and resistance being observed in cancer patients for all treatment strategies targeting these ligand-receptor pairs40. We speculated that the spatial distribution pattern of these ligand-receptor pairs may be associated with varied response rates and resistance to immunotherapies targeting them.

Our results indicated that neither stimulatory nor inhibitory ligand-receptor pairs were enriched in either SMSR or ICR (Fig. 4a and Supplementary Fig. 8a). However, a collection of stimulatory ligand-receptor pairs, notably CD27-CD70, was significantly enriched in LAR. In contrast, most inhibitory axes, including LAG3-LGALS3 and TIGIT-NECTIN2, were enriched in NCR. There were strong signals of TCR signalling in NCR, but no other immune activities were notable (Fig. 1j). However, NCR was depleted with T cells; thus, limited TCR signals within NCR were expected. The higher TCR signalling in NCR could be explained by the co-occurrence of higher inhibitory or negative immune responses within NCR, since the inhibitory regulation of T cells by immune checkpoint molecules requires TCR stimulation of T cells6.

Particularly, the expression of CD27-CD70 was higher in LAR blocks of Group A than in Group B (Fig. 4b, c). The expression of LAG3-LGALS3 and TIGIT-NECTIN was higher in Group B NCR blocks than in Group A (Fig. 4d). Both LGALS3 and NECTIN2 were over-expressed by cancer cells40,41, and Cancer-CEACAM6 was increased in the cancerous region of Group B than Group A (Fig. 3d, e and Supplementary Fig. 7i, j). This phenomenon partially explains why the expression of LAG3-LGALS3 and TIGIT-NECTIN was higher in Group B NCR blocks than in Group A. We next examined the expression of these immune checkpoint axes in the scRNA-seq data from tumour samples and delineated the cell subtypes expressing respective receptors and ligands, respectively (Supplementary Fig. 8b). We noted that CD70 was highly expressed in cDC-LAMP3, while CD27 was highly expressed in CD4-Treg-FOXP3 and CD8-Tex-LAYN. The potential cellular communication of cDC-LAMP3 with CD4-Treg-FOXP3 and CD8-Tex-LAYN through the axis of CD27-CD70 was also implicated solely based on the scRNA-seq data (Fig. 4e). We also observed the physical contact between PD1+CD27+CD8+ cells and CD70+LAMP3+ cells within LARs in samples from both the main study cohort and the independent validation cohort (Fig. 4f and Supplementary Fig. 8c). For the main study cohort, we measured the average distance between the PD1+CD27+CD8+ and CD70+LAMP3+ cells and between DAPI and CD70+LAMP3+ cells within each LAR. We found that the average distance between PD1+CD27+CD8+ and CD70+LAMP3+ cells was significantly lower than the average distance between DAPI and CD70+LAMP3+ cells (P=0.00018, paired t-test) (Fig. 4h and Supplementary Data 10). The proportion of PD1+CD27+CD8+ cells adjacent to LAMP3+ cells dropped over the distance ranges within LAR, as observed in the independence validation cohort (Supplementary Fig. 8d). LAR spatial blocks simultaneously expressing CD8A/CD8B, LAMP3, PDCD1, CD27, and CD70 were detected in a group of samples across cancer types (Supplementary Fig. 8e). Intriguingly, the association between the abundance of CD8-Tex-LAYN and the expression of CD27 was stronger in Group A than in Group B, whereas such a phenomenon was not observed for CD4-Treg-FOXP3 (Supplementary Fig. 8f). In addition, both axes of CD28 and CTLA4 were simultaneously expressed in a considerable amount of LAR blocks, and the majority of these blocks were found in Group A (Fig. 4g and Supplementary Fig. 8g, h), suggesting a potential equilibrium established in Group A LARs for CD28 and CTLA4 competing for binding of CD80 and CD86.

Methods
All relevant ethical regulations were complied with throughout the study, and written informed consent was gathered from all donors. Peking University Institutional Review Board and Peking University Cancer Hospital & Institute Ethical Committee approved the research and sample collection plan (IRB00001052-20037 and 2024KT121).

Sample collection
A total of 27 donors diagnosed with gastric cancer at the Peking University Cancer Hospital & Institute were enroled in this study. All donors were treatment-naive, as surgical resection of the nidus was the first choice of treatment when samples were collected. Before surgery, the peripheral blood (8 ml) was collected using anticoagulant tubes and put on ice before further processing. Tissues (tumour, adjacent normal, and lymph node) were immediately immersed into 5 ml MACS Tissue Storage Solution (Miltenyi Biotec, Germany) on ice after surgical resection. Tumour samples were equally cut into two pieces, each subjected to ST and scRNA seq library preparations, respectively.

Tissue dissociation and single-cell processing
The solid tissues (tumour, adjacent normal, and lymph node) were cut into small pieces (~1 mm3) in the RPMI 1640 Medium (Thermo Fisher Scientific Inc., USA) with 10% foetal bovine serum, loaded in gentleMACS C Tubes, and digested with MACS Tumour Dissociation Kit (130095929, Miltenyi Biotec) on a gentleMACS Dissociator with Heaters (Miltenyi Biotec, Germany). Continuous rotation at 37 °C for 30 minutes was applied. After dissociation, tissues that were still not dissociated were filtered by MACS SmartStrainers (70 μm). The resulting cell suspension was stained with 7AAD (prediluted, 00699350, Thermo Fisher) and CD235a (prediluted, 349114, Biolegend) before being sorted by flow cytometry with a BD FACSAria III sorter (BD Biosciences, USA). Roughly 10,000~20,000 cells were prepared for the 10 × 5’ and 3’ single-cell library construction. The guidance for the library construction is detailed in the manufacturer’s protocol (CG000207 and CG000183).

Construction of a single-cell atlas for GC
The sequencing reads of scRNA-seq data were aligned against the GRCh38 human reference genome by Cell Ranger (version 4.0.0, 10x Genomics Inc) to generate matrixes that record UMI counts of each gene across all cells. After alignment, cells with any one of the following features were filtered out using Seurat42,43 (version 4.4.0): (1) less than 500 or more than 50,000 UMI count detected in a cell; (2) less than 200 or more than 5000 genes detected; (3) the proportion of mitochondrial genes greater than 25%. Genes with expression in less than three cells were also filtered out. After quality control, the UMI counts were normalized and scaled using the function NormalizeData implemented in the R package Seurat (version 4.4.0). The score of the cell cycle was calculated using the function CellCycleScoring, and the score of the tissue dissociation effect was measured based on the expression of Heat Shock Protein (HSP) genes using AddModuleScore.

We next performed clustering using Scanpy44 (version 1.9.1) following steps below: selected 4,000 highly variable genes (HVGs) across platforms (10×3’ and 5’) using sc.pp.highly_variable_genes; regressed out UMI count, the proportion of mitochondrial genes, the score of the cell cycle, and the score of HSP using sc.pp.regress_out; performed principal component analysis (PCA)45; removed platform-induced batch effect using BBKNN46 (version 1.5.1); built the k-nearest neighbour graph47; performed clustering by the Leiden algorithm48; generated the Uniform Manifold Approximation and Projection (UMAP) embedding49. We performed a two-round clustering on the scRNA-seq dataset: eight major cell types were first identified based on the expression of canonical marker genes, and a second round of clustering was then performed to identify subtypes for each major cell type. Differentially expressed genes were identified using sc.tl.rank_genes_groups, with filters of adjusted P value < 0.01 and log2FC > 0.5.

Tissue embedding and optimisation
Fresh tissues were immediately embedded in OCT (optimal cutting temperature, Sakura Finetek, USA) after leaving the Tissue Storage Solution. The embedded tissue was submerged in isopentane in a stainless-steel beaker immersed in liquid nitrogen. The tissue was submerged until fully frozen. The OCT-embedded tissues were transferred on dry ice and stored at -80 °C before cryosection using Leica cryostat at -20 °C. Ten sections for each tissue block were pooled together, and total RNA was extracted using RNeasy Kits (Qiagen, Germany). RNA quality was assessed by evaluating RNA Integrity Number (RIN) on an Agilent 2100 Bioanalyzer (Agilent Technologies, USA). Only tissue blocks with RIN≥7 were subjected to spatial transcriptomics.

Visium Spatial Tissue Optimization Slide (10x Genomics, USA) was used to determine the optimal permeabilisation time. The OCT-embedded tissues were cryosectioned at a ten μm thickness and placed on capture areas of the optimisation slide equilibrated to -20 °C before usage. Seven tissue sections were prepared, one for a negative control (no permeabilisation reagents added) and the rest for a different permeabilisation time length (3, 6, 12, 18, 24, 30 minutes). Mouse reference RNA (RIN≥7) was used as a positive control. The tissue sections were first incubated at 37 °C for 1 minute and then fixed using methanol at -20 °C for 30 minutes. After tissue fixation, tissue sections were stained by sequential processing with isopropanol (1 min incubation), hematoxylin (7 min), bluing Buffer (2 min), and eosin mix (1 min) before 5 minutes incubation at 37 °C. Repeated ultrapure water immersion and air drying were adopted after certain staining stages, as illustrated in the manufacturer’s protocol (CG000238). The H&E (hematoxylin and eosin) stained sections were imaged using Leica Aperio Versa (Leica, Germany) with any of the following objectives: 5X, 10X, and 20X. After imaging, the permeabilisation enzyme was directly added onto each tissue section placed on the optimisation slide sequentially according to the above time length. Fluorescent cDNA was synthesised in situ, and the remaining tissue was removed before fluorescence imaging. The detailed imaging guidelines for both H&E and fluorescent sections can be found at CG000241. We selected 22 minutes as the optimal permeabilisation time for GC tumours, as the maximum fluorescence signal with the lowest signal diffusion was observed.

Spatial transcriptomics library construction
Visium Spatial Gene Expression Slides (10x Genomics, USA) were used to capture the poly-adenylated mRNA from intact tissue sections in situ. The steps of tissue fixation, H&E staining and imaging, and tissue permeabilisation in this section are the same as those of the tissue optimisation section. After 22 minutes of tissue permeabilisation, reverse transcription and second-strand synthesis proceeded before cDNA amplification. The cycle number of qPCR was determined on an Applied Biosystems 7500 Fast Real-Time PCR system (Thermo Fisher Scientific Inc., USA) for each sample. The library was then constructed based on the manufacturer’s user guide (CG000239). The spatial transcriptomic libraries were sequenced (2 × 150bp) on a NovaSeq platform at GENEWIZ (Suzhou, China).

Multiplex immunohistochemistry staining of tumour samples
Fresh tumour samples were first subjected to Formalin-fixed paraffin-embedded (FFPE). The FFPE-embedded tissues were cryosectioned at a four μm thickness and placed on poly-lysine-coated slides. We designed three mIHC experiments, and each assayed a different combination of protein markers using PANO 6-plex IHC kits (Panovue, Beijing). A total of 26 samples from the main study cohort were assayed. For each combination of protein markers, primary antibodies were sequentially added before HRP-conjugated secondary antibody incubation and tyramide signal amplification, and nuclei were stained with DAPI. All slides were scanned using the Olympus VS200. The protein markers used for staining were CD8 (1:400, clone C8/144B, 70306, Cell Signalling Technology), TCF1 (1:200, clone C63D9, 14456, Cell Signalling Technology), Ki67 (1:200, clone SP6, RM-9106-s1, Thermo Scientific), CD27 (1:2000, clone LPFS2/1611, ab268144, Abcam), CD70 (1:400, clone E3Q1A, 69209, Cell Signalling Technology), PD-1 (1:300, clone 29F.1A12, 135210, Biolegend) and LAMP3 (1:200, clone 1010E1.01, Novus Biologicals). According to the same experimental design, fifteen samples from an independent validation cohort were assayed at AlphaX, Beijing.

Preprocessing of ST data
ST is a feasible and applicable option since its UMI-based format is identical to that of scRNA-seq; thus, a smooth transition from scRNA-seq to ST can be applied in our study to delineate the spatial distribution of cell subtypes. First, the sequencing reads were aligned against the GRCh38 human reference genome by Space Ranger (version 1.1.0, 10x Genomics Inc), whereby UMI counts and capture spots under tissue were reported. After mapping, each capture spot was annotated with clinical information, and capture spots with a UMI count of less than 500 and features detected of less than 200 were filtered out. After quality control, the UMI counts were normalized and scaled using the function SCTransform implemented in the R package Seurat (version 4.1.1), and a list of variable genes was generated.

Before clustering, we used the function CalculateRogue in the R package ROGUE31 (version 1.0) to measure the purity score for each ST dataset, and we observed heterogeneity across the 33 ST datasets. Thus, we categorized the 33 ST datasets into three groups based on the purity score generated from the above step. Then, based on the categorization, we determined the resolution score for each group. After excluding genes of a blocklist described by Zheng et al.19 from the list of variable genes, we performed clustering for each dataset using Seurat42,43 by resolution scores determined in the above step. Then, functions of RunPCA, FindNeighbors, and FindClusters were applied sequentially to each ST dataset. As a result, the number of clusters ranges from seven to 11 across all ST datasets, with 521 “ST clusters” reported.

Deconvolution of ST spots
To reduce noise as much as possible, we only adopted the scRNA-seq data generated from the 10× 5’ platform. For each cell subtype, we defined a “central point” based on the principal component of all cells within the subtype and selected the 1,000 nearest cells to the central point based on the Euclidean distance. This cell subset was used as the reference dataset for the subsequent deconvolution. The cellular composition of each capture spot was measured using RCTD32 (version 1.2.0) by deconvolving the transcriptome of each capture spot with the above scRNA-seq dataset. Cell subtypes derived from blood samples were excluded from this analysis since these cell populations were unlikely to reside in solid tumour tissues.

Construction of a spatially informed cellular atlas
After deconvolution using RCTD32 (version 1.2.0), a matrix with rows representing cell subtypes and columns representing capture spots was generated for each dataset, with each score indicating the enrichment degree of a cell subtype within a particular capture spot. For each matrix, the enrichment score for each cell subtype was z-score transformed across all capture spots before all matrices were merged into one. Based on the merged matrix, we calculated the mean value of each cell subtype for each ST cluster identified in the above section and then performed a z-score transformation for each cell subtype across all ST clusters.

With the newly generated z-score transformed matrix, we calculated the Euclidean distance for each pair of ST clusters before hierarchically clustering them by Ward’s minimum variance method. Next, we sequentially validated the number of clusters starting from 10 by mapping the cluster results onto the H&E images for each ST dataset. Clustering ST clusters into four groups revealed the optimal result, as they maximally reflected the major pathological patterns (lymphoid region, smooth muscle and stromal region, and cancerous region) and confirmed by pathological physicians for all ST datasets. Furthermore, the cancerous region was divided into two sub-regions, with one exhibiting lower cancer cell scores but higher lymphocyte cell scores and the other displaying the opposite trends. Accordingly, these groups were termed as “Lymphocyte Aggregated Region” (abbreviated as LAR), “Smooth Muscle and Stromal Region” (SMSR), “Immunogenic Cancer Region” (ICR), and “Negative-immunogenic Cancer Region” (NCR), respectively.

LAR was determined by integrating scRNA-seq and ST data and comprised different sizes (measured by the number of spots) of spatial blocks. It consisted of less organised histological patterns, such as immune cell aggregates (especially for lymphocytes) and more functionally developed structures, such as TLSs. Some TLSs were visible directly on H&E images and exhibited strong TLS-associated signals, such as CXCL13, CCL19, and CCL21; the other spatial blocks in the LAR region were small but also showed strong signals of the above markers when compared with other areas using transcriptomic data.

To reveal the spatial distribution of all cell subtypes across these four regions, we used the aforementioned merged matrix to calculate the mean value of each cell subtype for each region and then performed a z-score transformation for each cell subtype across all regions. Finally, the z-score transformed matrix was used to calculate Euclidean distance for each pair of regions and cell subtypes before hierarchical clustering using Ward’s minimum variance method. The result is visualized using a heatmap implemented in the R package ComplexHeatmap (version 2.10.0).

Identification of hallmark genes for ST groups
Since there is no reliable method to eliminate the batch effect in ST data, we adopted the effect size to evaluate gene expression variation across samples. Specifically, UMI counts of each ST slide were normalized by SCTransform implemented in the R package Seurat (version 4.1.1), and all spots on each ST slide were dichotomized based on the ST cluster (e.g., LAR versus non-LAR). The point-biserial correlation coefficient50 was then calculated to determine whether or not the gene expression was unique to an ST group identified on an ST slide. The Benjamini & Hochberg (BH) approach51 implemented in the R function p.adjust was then used to correct p-values, and genes with an adjusted p-value greater than 0.05 were filtered out. Next, the effect_sizes function implemented in the R package esc (version 0.5.1) was used to convert the point-biserial correlation coefficient to a standardized mean difference effect size. Finally, for each ST group, the function rma implemented in the R package metafor (version 3.8.1) was used to combine gene effect sizes across samples. Again, the BH approach was used to correct p-values, and only genes with an adjusted p-value lower than 0.05 were retained.

Identification of spatial blocks for ST groups
After clustering all spots into ST groups, we identified spatial blocks based on spot coordinates provided by Visium. A spatial block consists of a group of spots from the same ST group, and they are spatially clustered, which means that each spot is adjacent to at least another spot within the block. Specifically, for all spots labelled with the same ST group on a slide, the distance between spots of any pair was calculated using the Euclidean distance. A cutoff value (<=2) was used to determine if the two spots were spatially adjacent. We then customised homemade R scripts to isolate spatial blocks, mainly utilising an R package for network analysis (igraph, version 1.5.1). Neither lower nor upper limit was set to restrict the size of spatial blocks.

Function annotation of ST groups
Gene Set Enrichment Analysis (GSEA)52 was performed to measure the intensity of signalling pathways in each ST group, using the function GSEA implemented in the R package clusterProfiler (version 4.2.2) with the “rank” mode. The annotation database (Reactome37, version 7.5.1) was downloaded from MSigDB53 (https://www.gsea-msigdb.org/gsea/msigdb).

The ratio of observed versus expected
We generated a contingency table of ST groups (LAR, SMSR, ICR, and NCR) by stages (I, II, and III), with each table cell indicating the number of spots for a given combination of ST group and stage. For each cell, the expected number was calculated by the product of the row total and column total divided by the total number of all spots. The same analysis was used to measure the enrichment of ST groups across different subtypes. We also applied the chi-squared test to determine whether the associations between the ST groups and clinical sub-categories are significant.

Cellular interactions coordinated by immune checkpoint axes
We compiled a list of immune checkpoint axes (ligands and receptors) based on previous studies40,41. The raw UMI count of each gene in the list was extracted and aggregated by spatial block. Only spatial blocks expressing both ligand and receptor of an axis were considered positive blocks for that axis. To identify immune checkpoint axes enriched in each ST region, a permutation approach similar to CellPhoneDB54–56 was adopted: the actual mean value of axes for each positive spatial block was calculated; then, ST group labels of all spatial blocks were permutated 1000 times, and each time the mean value of an axis was calculated to generate a null distribution; a P-value was obtained based on the actual mean value of axes and the null distribution. A similar approach was used to infer cellular interactions from the scRNA-seq data.

Survival analysis using TCGA data
We downloaded the UCSC TOIL RNA-seq recompute version of the TCGA data from Xena (https://xenabrowser.net), kept samples from patients with valid information on overall survival, and calculated the enrichment score of the top 30 LAR-associated genes for each sample using GSVA57 (version 1.42.0). Next, we dichotomised the patients into high and low groups based on the GSVA enrichment score. We built a Cox proportional-hazards model to examine if an association exists between the patient survival time and these two groups after adjusting for age and gender, with a hazard ratio (HR) for each variable calculated.

Detecting individual cells on H&E images using RedeHist
We implemented an artificial intelligence-guided algorithm, RedeHist, to detect individual cells on H&E images. The diameter of each ST spot was adjusted to 128 pixels before nuclei segmentation and stain colour normalisation. A scRNA-seq reference was generated by selecting 1000 highly variable genes using ROGUE31. The processed H&E images, the nuclei segmentation data, the scRNA-seq reference, and the spatial transcriptomics data were used as input for the RedeHist. Details about RedeHist can be found at doi: 2024.06.17.599464.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_3": "The authors have provided detailed and thoughtful responses to the extensive comments from Reviewer 1. They have addressed key points such as the limitations of ST analysis, the validity of RCTD, mIHC validation, the confusion between LAR and TLS, issues with LAR/TLS annotation, grouping of spatial blocks, and the cross-section of slides. Each comment was met with specific responses and improvements in the manuscript, except for the following points:\n1. The statement \"provides several markers that could be predictive of immunotherapy response\" in the newly added abstract is an overstatement. As the study does not identify specific markers for immunotherapy, it would be more accurate to conclude that the study provides new perspectives in the exploration of immunotherapy biomarkers.\n2. Regarding the definition of TLS and LAR in relation to the argument 2), the legends of Figures 2 and 3 still refer to TLS. Even though HE staining is used for verification, it is difficult to conclusively demonstrate the presence of follicular structures or HEVs indicative of true TLS. It would be more appropriate to revise these descriptions to LAR.", "reviewer_4": "Revision Notes [editorial note: journal name and manuscript number redacted]\n1) Ext Fig 3A-C-D-E-F lack y axis annotation (is it a proportion over total CD8+ T cells?)\n2) L130 - The whole purity analysis is not very clear in terms of intent and problem that the authors try to tackle, since they end up performing deconvolution at the spot level either way. It appears from the methods that the purity is used to stratify samples in different categories as a way to determine clustering resolution. The mention of these concepts in the main text can help the reader transition from the concept of tumor purity, its relation to sample morphology and expected clustering results.\n3) L131: “all spots from an individual ST sample are transcriptomic similar” should be rephrased to be clearer\n4) Why is the clustering returning 4 regions? Is there any justification? This is a major step in the analysis they decided to adapt clustering results to match the expected tissue morphology identified through H&E. – not really an issue per se but they are essentially reducing the potential resolution of ST data to features visible by pathologists conceptually speaking.\n5) L158 - The referred figures (Ext Fig 4 4b-c) do not show what is described in the text. The proportion of spots occupied by LARs is apparently reported in Fig1 F instead. The expression of GC markers is reported in the heatmap in Fig1E and the same concept (albeit with slightly different markers) is reported in Ext. figure 4d (CXCL13, CCL19, CCL21)\n6) L161 - On the same line, the results reported on the differences between ICR and NCR for GC markers are related to plots shown in Ext. Fig 4E and not D as reported.\n7) Always in relation to survival analysis, I think it could make sense to estimate survival probabilities also for the genes enriched in the other 3 compartments, to see if positive/negative are expected for the other areas (e.g. if the ICR signature predictive of poor survival?)\n8) L186 - Minor: The conclusions of the initial phrase of the paragraph could be placed in the context of the RCTD discussion above, around line 145. Placed here the conclusions do not seem to emphasize the message that follows, which could be standalone in this setting. Moreover, I suggest a bit more explanation as to how the conclusion was derived (i.e. RCTD enrichment scores being consistent relate to stable cell proportions less influenced by original tissue morphology?)\n9) Fig 1G does not seem to be referenced in the main text\n10) L245 - There are still some doubts related to the DCs statements. the authors speculate that the higher signal for “Antigen Processing Cross presentation” in NCR and highest in LAR reflect the trajectory between cDC1-CLEC9A to cDC-LAMP3 so they think that the signal in LAR is promoted by the presence of enriched cDC1-CLEC population and the signal in NCR is given by the presence of not enriched nor depleted cDC-LAMP3, this is not sintomatic of a temporal trajectory and also the signal for the GO term in the two areas (especially LARs) could be promoted by the abundance of other antigen-presenting cells including macrophages. Additionally, cDC1-CLEC9 are not enriched in NCR regions. (this was also a point raised in the previous revision)\n11) L263 - results should be backed up by showing whether the same signatures are also correlated to abundance of these cell types within non-LAR blocks.\n12) L267 - Where does TCF7 come from? Do they show it in any previous results related to CD8?\n13) Figure 2 can benefit from the addition of a quantitative assessment of the co-localization of fluorescent signals instead of just showing 2 ROIs on one sample. Nowhere in the supplementary there is a graphical depiction and a formal analysis on the co-localization patterns of fluorescence signals across the cohort of 11 (+ 5 validation) mIHC samples studied by the authors.\n14) L282 - A brief description of the interpretability of the utilized coefficient of variation (CV) would be beneficial in better understanding the magnitude of the reported results.\n15) L291 - Despite having just shown how LAR blocks can greatly vary within a single slide, the authors decide to group them toghether by donor, a choice which in the author’s view is justified by assessing LAR’s TME disparity across donors. I would argue that retaining a finer LAR identity and clustering them retaining donor information could potentially pinpoint inter-LAR differences conserved and not conserved across donors. Ultimately this could provide a higher levels of granularity to the results.\n16) L303 - “reflecting a more mature LAR status” is a phrase which is unsupported by the data given that the authors do not explain what they expect represent an “immature” LAR.\n17) Line 306 - “ICR was more enriched in Group A than Group B, suggesting a more deteriorated cancerous state in Group B” this would be more informative if the authors could relate this information with the different staging distributions across group A (mostly stage III) and group B (others). Also, it would be informative to elaborate on what the authors consider a “more deteriorated cancerous state”.\n18) L307 - Typo with double open parenthesis\n19) Fig 3D needs annotation of x axis (abundance in terms of proportion over total?). It is also overlapping other panels in the figures which impairs overall interpretability.\n20) L329 - I think a more global picture of the LR interactions detected (also ones in ICR/SMSR) regions could be beneficial. A criteria (or hypothesis) for the choice of the specific LR pairs should be clearly stated\n21) L333 - The sentence “Since the inhibitory regulation of T cells by immune checkpoint molecules accompanied by TCR stimulation of T cells, the co-occurrence of higher inhibitory or negative immune responses within NCR could explain the higher TCR signalling within NCR.” should also discuss the fact that NCR are relatively T cell poor regions hence one should expect the minority of interactions including TCR to happen here, no?\n22) L339 - Do NCR regions across the two patient groups (originally determined by LAR composition) differ based on composition? This could be explanatory of the differences in LR couples activity and maybe aid in interpreting the results, strenghtening the author’s conclusions.\n23) Ext. Fig 8 A,B – Minor: the use of a sequential rather than a divergent colormap in these plots could help interpretation and emphasize the expression values of respective LR pairs, given that the values plotted are not diverging but instead only sequential.\n24) L348 - Fig 4F refers to “TLS” whereas the text refers to LAR regions.\n25) Fig. 4f and Extended Data Fig. 8c should be accompanied by a more thorough analysis quantifying the co-occurrence of these signals to better understand their conservation across LARs and not just a couple of ROIs within specific samples (from both the main and validation cohorts), this can further strengthen the message of the authors. Including a graphical depiction of what is stated in the revised phrase hinting at supp. table 10 in the form of a boxplot (both from the main and validation cohorts) can greatly aid in ensuring the reader about the conclusions drawn by the authors."}, "author_response": "Thank you very much for your thoughtful comments. We have rephrased this sentence based on your suggestion in the abstract: \"Our study unveils the GC TME in a spatially resolved manner and provides new insights into the exploration of immunotherapy biomarkers.\"\n2. Thank you for your suggestion. We have revised the related figure legends in the manuscript.\n1) Thank you for your comment. The y-axis of the aforementioned figures shows the proportion of each cell subtype within its respective major cell type, for example, the proportion of CD4-Treg-FOXP3 within CD4 T cells. We have rephrased the figure legend: \"Bar plots show the subtype composition for CD4 T cells, CD8 T cells, and myeloid cells, respectively.\"\n2) Thank you so much for your insightful suggestions. We have added the rationale of using the purity analysis to stratify samples into different groups for cluster resolution determination in the main text (L131-136).\n3) Thank you for pointing out the confusion. We have rephrased this sentence in the main text: \" We suspected that if a single histological pattern dominated a slide, ST spots from that slide could have similar transcriptomic profiles, which can be directly measured by the purity score — a higher purity score indicates more similarities were observed in the transcriptomic profiles of all spots within an ST section.\"\n4) Thank you for your comments. We agree with your opinion that clustering all ST spots into four groups, rather than more groups, could reduce the potential resolution of ST data. However, we argue that the alignment between ST groups and the well-established histological patterns is a crucial standard for validating the robustness of our clustering approach.\nIn the method section of the original manuscript, we provided the steps and rationale of clustering all ST spots into four groups (L574): \"With the newly generated z-score transformed matrix, we calculated the Euclidean distance for each pair of ST clusters before hierarchically clustering them by Ward's minimum variance method. Next, we sequentially validated the number of clusters starting from 10 by mapping the cluster results onto the H&E images for each ST dataset. Clustering ST clusters into four groups revealed the optimal result, as they maximally reflected the major pathological patterns (lymphoid region, smooth muscle and stromal region, and cancerous region) confirmed by pathological physicians for all ST datasets.\"\n5) Thank you for your comments. We have corrected the reference to Fig. 1f in the revised manuscript. Although the same markers were used for plotting, there are major differences between Fig. 1e and Extended Data Fig. 4d. Fig. 1e showed the average expression of these markers for each spatial region. Extended Data Fig. 4d delineated the expression pattern of these markers across all ST spots for each ST region.\n6) Thank you very much. We have corrected this error in the revised manuscript.\n7) Thank you for your suggestion. We have performed survival analysis for enriched genes of the other three ST regions, respectively. Our results indicated that higher expression of SMSR-enriched, ICR-enriched, and NCR-enriched genes was associated with worse prognosis in multiple cancers, respectively. We have upgraded Extended Data Fig. 4g with the new results and added the below sentence in the revised manuscript (L173-L175): \"In contrast, higher expression of SMSR-enriched, ICR-enriched, and NCR-enriched genes was associated with worse prognosis in multiple cancers (Extended Data Fig. 4g).\"\n8) Thank you. We have inserted the first sentence of the original L186 paragraph into the original L145 paragraph as you suggested. We also added how the conclusion was derived in the main text: \"Furthermore, we noticed that the distribution abundance of each cell subtype varied in two contrasting tendencies — most lymphocyte subtypes exhibited smaller deviations while non-lymphocyte cell subtypes demonstrated larger ones (Fig. 1d). Since the RCTD enrichment scores, to some degree, can reflect the cellular proportion regardless of the tissue morphology, smaller deviations observed in lymphocyte abundance across samples are probably associated with the migration characteristics of T cells and B cells even in solid tissue like tumours.\"\n9) Thank you. Fig. 1g was cited at L162-L164: \"Each ST group was identified either in all samples (ICR and NCR) or most samples (LAR and SMSR), composed of spots distributed evenly across all samples, and possessed a unique transcriptional profile (Fig. 1e-g and Supplementary Table 6).\"\n10) Thank you for your comments. According to your suggestion, we have rephrased this paragraph in the revised manuscript (L252-L257): \"In addition, the signal of \"Antigen Processing Cross Presentation\" was also higher in NCR. Antigen-presenting cells, such as cDC-LAMP3 and Macrophage-C1QC, were depleted in NCR (Extended Data Fig. 5h). Other antigen-presenting cells, including but not limited to cDC1-CLEC9A, cDC2-CD1C and Macrophage-VCAN, were not depleted in NCR, a region mostly dominated by Cancer-CEACAM6. Thus, we speculated that the higher signal of \"Antigen Processing Cross Presentation\" in NCR could be due to a combined effect of antigen processing activities of these cells.\"\n11) Thank you for your suggestion. Within the non-LAR region, we have measured the correlation coefficients between the cellular abundance of CD8-Tn-CCR7 and signals of TCR signaling, Proliferation signaling, Costimulation By The CD28 Family, and Antigen Processing Cross Presentation. Negative associations were observed between the cellular abundance of CD8-Tn-CCR7 and TCR signaling (correlation coefficient: -0.031, P = 0.06), Costimulation By The CD28 Family (-0.029, P = 0.07), and Antigen Processing Cross Presentation (-0.059, P = 0.0002). An insignificant association was observed between the cellular abundance of CD8-Tn-CCR7 and Proliferation signaling pathways within the non-LAR region (0.022, P = 0.17).\nWe have added the following sentence in the revised manuscript (L264-L268): \"In contrast, within the non-LAR region, negative associations were observed between the cellular abundance of CD8-Tn-CCR7 and TCR signalling (correlation coefficient: -0.031, P = 0.06), Costimulation By The CD28 Family (-0.029, P = 0.07), and Antigen Processing Cross Presentation (-0.059, P = 0.0002). An insignificant association was observed between the cellular abundance of CD8-Tn-CCR7 and Proliferation signalling pathways within the non-LAR region (0.022, P = 0.17).\"\n12) Thank you for your questions. TCF7 is a master regulator of T cell development and function, and it serves as a marker to distinguish naive T cells. It is not exclusively expressed on CD8 T cells. In our scRNA-seq data, we found that 47.50% of CD8-Tn-CCR7 cells expressed TCF7. Thus, TCF7 (protein: TCF1) is a reliable marker for tagging naive CD8 T cells when designing mIHC experiments.\n13) Thank you for your suggestion. We have counted the number of CD8+ cells and the number of TCF1+KI67+CD8A+ cells in LAR in 11 tumour samples from the main study cohort and five samples from an independent validation cohort within LAR. An average of 12.19% (SD = 0.0839) CD8+ cells within LAR were TCF1+KI67+ from samples in the main study cohort, while a lower average of 0.77% (SD = 0.007) CD8+ cells within LAR were TCF1+KI67+ positive from samples in the independent validation cohort (Supplementary Table 8).\nWe have updated the results in the revised manuscript (L278-L281): \"Within LAR, we counted the number of CD8+ cells and the number of TCF1+KI67+CD8+ cells and found an average of 12.19% (SD = 0.0839) CD8+ cells were TCF1+KI67+ cells from samples in the main study cohort, while a lower average of 0.77% (SD = 0.007) CD8+ cells were TCF1+KI67+ cells from samples in the independent validation cohort (Supplementary Table 8).\"\n14) Thank you for your suggestion. We have rephrased the sentence to introduce more about the coefficient of variation in the main text (L291): \"We observed that the expression of key gene features (FDCSP and CCL19), which was measured by the coefficient of variation (CV) — a higher CV (> 1) indicates greater relative variability, varied drastically among the LAR blocks (a CV value of 3.41 for FDCSP and 1.28 for CCL19), and a significant association was observed between these two trends (Fig. 3a).\"\n15) Thank you for your suggestion that clustering all LARs, rather than aggregating them by donors and then clustering, could likely generate inter-LAR differences conserved and not conserved across donors. However, our study goal is to depict the influence LARs exert on the TME or vice versa. In addition, distinct LARs were observed on a single ST slide. Therefore, it is rational for us to aggregate LAR blocks by donors and dichotomize them based on their cellular abundance.\nIn addition, our results partially revealed the differences between LARs — the expression of key gene features (FDCSP and CCL19) varied drastically among the LAR blocks, and the abundance of cell subtypes (GCB-LRMP and cDC-LAMP3) also varied drastically among the LAR blocks, suggesting the state of LARs is a spectrum.\n16) Thank you for your suggestion. We agree with your opinion that using the term \"mature\" to describe LARs in Group A was not rigorous. We have rephrased this sentence in the revised manuscript to avoid any confusion (L310-L315): \"More transcriptomic similarities were observed between TDLNs and Group A blocks than between TDLNs and Group B blocks (Extended Data Fig. 7f). A higher expression of LAR-enriched genes, mostly also highly expressed in TDLNs, was observed in bulk samples of Group A than B, further substantiating that the gene expression pattern of Group A LARs was more similar to that of TDLNs than Group B LARs (Extended Data Fig. 7g).\"\n17) Thank you for your suggestions. Fig. 3c showed that most Group A LARs, rather than ICR, were derived from Stage III tumor samples. Thus, we were unable to relate ICR enrichment in Group A with patient stage information.\nWe have elaborated more about ICR enrichment in Group A than Group B in the revised manuscript (L317-L320): \"Moreover, compared with other ST regions, ICR was more enriched in Group A than Group B, but NCR was comparable between Groups A and B (Extended Data Fig. 7h). Since ICR has higher immune cell signals but lower cancer cell signals than NCR (Fig. 1d), our results suggest a more deteriorated cancerous state in Group B samples.\"\n18) Thank you. We have corrected this typo in the revised manuscript.\n19) Thanks for your comments. We have added one sentence in the figure legend to clarify the definition of cellular abundance for the x-axis annotation: \"The cellular abundance is the proportion of each cell subtype in all cell subtypes.\" In addition, we have reorganized the panels in Fig. 3d to enhance its interpretability.\n20) Thank you so much for your suggestions. In our analysis, we included immune checkpoint LR pairs that have been successfully utilized in immunotherapies for cancer treatment. However, varied response rates and resistance were observed in cancer patients for all treatment strategies targeting these LR pairs (Morad et al., Cell, 2021). Thus, the spatial distribution pattern of these LR pairs may be associated with varied response rates and resistance of immunotherapies targeting them.\nWe have added the rationale of selecting LR pairs in the revised manuscript at (L341-L346): \"In our analysis, we included immune checkpoint ligand-receptor pairs that have been successfully targeted in immunotherapies for cancer treatment, despite varied response rates and resistance being observed in cancer patients for all treatment strategies targeting these ligand-receptor pairs. We speculated that the spatial distribution pattern of these ligand-receptor pairs may be associated with varied response rates and resistance of immunotherapies targeting them.\"\nOverall, our results showed that none of these ligand-receptor pairs were enriched in either SMSR or ICR (Fig. 4a and Extended Data Fig. 8a). We did not selectively pick ligand-receptor pairs enriched in LAR or NCR to display in the figures.\n21) Thank you for your suggestion. We have rephrased the sentence in the main text (L351-L356): \"There were strong signals of TCR signalling in NCR, but no other immune activities were notable (Fig. 1j). However, NCR was depleted with T cells; thus, limited TCR signals within NCR were expected. The higher TCR signalling in NCR could be explained by the co-occurrence of higher inhibitory or negative immune responses within NCR, since the inhibitory regulation of T cells by immune checkpoint molecules requires TCR stimulation of T cells.\"\n22) Thank you very much for your thoughtful comments. We fully agree with your opinion that the expression differences of LR pairs within NCR between Groups A and B may be associated with the cellular composition differences within NCR between Groups A and B. However, comparisons for the cellular composition of NCR between Groups A and B are less achievable for multiple reasons: as mentioned at L322, NCR and ICR were spatially tangled with each other; all ST samples, regardless of whether they were grouped as A or B, comprise both NCR and ICR (Fig. 1g); NCR was enriched comparably in Groups A and B (Extended Data Fig. 7h). Therefore, for comparisons between Groups A and B, we combined NCR and ICR as cancerous regions.\nOur results showed that the expression of LAG3-LGALS3 and TIGIT-NECTIN was higher in Group B NCR blocks than in Group A (Fig. 4d). Both LGALS3 and NECTIN2 were over-expressed by cancer cells (Morad et al., Cell, 2021), and Cancer-CEACAM6 was increased in the cancerous region of Group B than Group A (Fig. 3d-e and Extended Data Fig. 7i-j). This phenomenon partially explains why the expression of LAG3-LGALS3 and TIGIT-NECTIN was higher in Group B NCR blocks than in Group A. To strengthen our conclusion, we have added the above explanation in the revised manuscript (L360-L363).\n23) Thank you for your suggestion. We agree with your opinion that the expression values of LR pairs displayed in Figs. 8a-b were sequential rather than divergent. However, these two figures aimed to highlight the expression differences of LR pairs across different ST regions or cell subtypes. Thus, we argue that the red-yellow-blue color system, which uses red to represent high values, blue to represent low values, and yellow to represent intermediate values, is better to highlight LR pairs highly expressed in an ST region than a single-color system.\n24) Thank you for your suggestion. We have changed \"TLS\" to \"LAR\" in the legend of Fig. 4f in the revised manuscript.\n25) Thank you for your suggestion. We have generated a box plot comparing the average distance between PD1+CD27+CD8+ and CD70+LAMP3+ cells and the average distance between DAPI and CD70+LAMP3+ cells within each LAR in the main study cohort. In addition, we have generated a bar plot that showed the proportion of PD1+CD27+CD8+ cells adjacent to LAMP3+ cells at different distance ranges within LAR, as observed in the independence validation cohort. We have updated the results in the revised manuscript (L373-L378): \"We found that the average distance between PD1+CD27+CD8+ and CD70+LAMP3+ cells was significantly lower than the average distance between DAPI and CD70+LAMP3+ cells (P = 0.00018, paired t-test) (Fig. 4h and Supplementary Table 10). The proportion of PD1+CD27+CD8+ cells adjacent to LAMP3+ cells dropped over the distance ranges within LAR, as observed in the independence validation cohort (Extended Data Fig. 8d).\""}, {"round": 2, "reviewer_comments": {"reviewer_3": "The authors have adequately addressed the points raised by reviewers #3 and #4."}, "author_response": "Thank you, and we really appreciate the time you took to evaluate our paper."}]

Decision letter