Spatial transcriptomics reveals tumor niche heterogeneity in advanced gastric cancer and suggests a tumor-promoting role for KRT6B

Abstract

BACKGROUND

Gastric cancer (GC) is highly heterogeneous, which complicates its diagnosis and treatment. Although single-cell transcriptomic approaches can characterize cellular heterogeneity, they do not preserve spatial information. They therefore cannot adequately resolve spatial heterogeneity within tumors. We hypothesized that full-thickness spatial transcriptomic profiling would reveal intratumoral spatial heterogeneity and identify distinct tumor niches in advanced GC (AGC).

AIM

To characterize intratumoral spatial heterogeneity and identify distinct tumor niches in AGC.

METHODS

Full-thickness AGC tissues were profiled using Stereo-seq. Spatial niches were identified using BANKSY, expression trends were analyzed using Mfuzz, and niche-specific signaling was assessed through ligand-receptor analysis. Functional programs were evaluated using gene set enrichment analysis, gene set variation analysis, and enrichment analysis of differentially expressed genes. KRT6B was assessed by immunofluorescence, immunohistochemistry, analysis of public datasets, and survival analysis, and its function was examined in KRT6B-overexpressing GC cells using proliferation, migration, and invasion assays.

RESULTS

We identified 13 spatial niches, demonstrating marked intratumoral spatial heterogeneity. Three tumor-dominant niches (N5, N4, and N2) formed a superficial-to-deep axis and exhibited graded transcriptional programs, with stronger stemness-related and motility-related features in superficial regions and increased immune-related and metabolism-related programs in deeper regions, along with niche-specific intercellular signaling patterns. We also identified a spatially isolated tumor nest (N9) with relatively stronger transcriptional features for proliferation, metabolism, and immune/inflammation, and preferential enrichment for KRT6B. KRT6B was enriched in tumor nests and at the invasive front, upregulated in GC tissues, and associated with poor prognosis. In vitro, KRT6B overexpression promoted GC cell proliferation, migration, and invasion.

CONCLUSION

Full-thickness spatial transcriptomics revealed marked heterogeneity in AGC and implicated KRT6B as a candidate tumor-promoting molecule for further investigation.

Key Words: Gastric cancer; Spatial transcriptomics; Stereo-seq; Spatial heterogeneity; KRT6B

Core Tip: The spatial organization of advanced gastric cancer remains insufficiently characterized. Using full-thickness Stereo-seq profiling, this study identified a superficial-to-deep tumor niche axis, niche-specific transcriptional gradients, and a spatially isolated tumor nest with preferential enrichment for KRT6B. External validation and in vitro assays further supported KRT6B as a candidate molecule associated with aggressive tumor phenotypes. Although the study was based on a limited number of samples, these findings provide an initial spatial framework for investigating intratumoral heterogeneity in advanced gastric cancer.

INTRODUCTION

Gastric cancer (GC) remains a major global health burden, ranking fifth in incidence and fourth in cancer-related mortality worldwide[1]. Despite advances in multimodal therapies, including surgery, chemotherapy, radiotherapy, and immunotherapy, the prognosis remains poor, with approximately 800000 deaths annually[2].

GC exhibits marked intertumoral and intratumoral heterogeneity at both the morphological and molecular levels, which complicates diagnosis and treatment. Although single-cell RNA sequencing (scRNA-seq) has substantially improved our understanding of tumor heterogeneity, tissue dissociation disrupts spatial organization and limits the analysis of in situ tumor architecture. Spatial transcriptomics (ST) addresses this limitation by enabling gene expression mapping in intact tissues and interrogation of spatial structure-function relationships. Among current ST platforms, Stereo-seq (SpaTial Enhanced REsolution Omics-sequencing) provides nanoscale spatial resolution (down to 500 nm) and a large field of view, enabling both fine structural delineation and global tissue mapping[3]. However, the application of ST in GC remains at an early stage.

In this study, we applied Stereo-seq to advanced GC (AGC) tissues and performed full-thickness ST profiling from the mucosa to the serosa, enabling high-resolution spatial mapping across the gastric wall. We identified marked spatial heterogeneity and gradients in tumor cell states and tumor-microenvironment interactions and, in one specimen, detected a spatially isolated tumor nest with distinct functional features and marked KRT6B upregulation. Further clinical and in vitro validation showed that KRT6B was associated with GC clinicopathological features and promoted GC cell proliferation and invasion. Collectively, these findings extend the application of Stereo-seq in GC, provide new insights into the spatial heterogeneity of GC, and suggest that KRT6B may serve as a potential molecular target involved in GC progression and therapeutic intervention.

MATERIALS AND METHODS

Sample collection

Tumor tissues were collected from five patients (P1-P5) with primary advanced gastric antral carcinoma. None had received neoadjuvant therapy, and all were clinically staged as T3 before surgery. Within 30 minutes after resection, full-thickness gastric wall specimens, extending from the mucosa to the serosa, were obtained from tumor regions, trimmed to ≤ 10 mm × 10 mm × 5 mm, embedded in optimal cutting temperature compound, and stored at -80 °C. Frozen sections were stained with hematoxylin and eosin (HE), and tumor regions were annotated by a pathologist. Based on tissue quality, preservation of full-thickness gastric wall architecture, and the suitability of the tumor regions for ST analysis, P2 and P4 were selected for ST profiling (Supplementary Table 1).

Tumor tissues and matched adjacent non-tumor tissues (≥ 5 cm from the tumor margin) were also collected from 54 patients with GC (Supplementary Table 2). These specimens were fixed in 40 g/L paraformaldehyde, paraffin-embedded, and used for immunohistochemistry (IHC) and immunofluorescence (IF).

Stereo-seq data processing and integration

ST data were generated using Stereo-seq as previously described[3]. Tissue sections were captured on Stereo-chips pre-coated with DNA nanoballs (DNBs), and Bin 50 (50 × 50 DNBs; 25 μm square) was used as the analysis unit, hereafter referred to as a spot. Raw FASTQ reads were aligned to the human reference genome (hg38). Regions of interest were manually delineated on HE images to exclude background and were analyzed using Seurat (v4.4.0)[4]. Spots with fewer than 300 detected genes or more than 20% mitochondrial transcripts, as well as genes detected in fewer than five spots, were filtered out. The data were log1p-normalized, the top 2000 variable genes were selected, and principal component analysis was performed. Batch effects across samples were corrected using Harmony[5]. The top 30 principal components were then used for spot-level Louvain clustering (resolution = 0.3).

Niche clustering

Niche clustering of ST data was performed using BANKSY[6], which integrates spot-intrinsic transcriptional profiles with the spatial context of neighboring spots. Specifically, each spot was embedded into a composite space that combined its own transcriptome with those of neighboring spots. The RunBanksy function was run with k_geom = 30 and lambda = 0.3. The resulting BANKSY assay was used for downstream dimensionality reduction and spatial domain identification.

Differentially expressed genes and pathway enrichment analysis

Differential expression analysis was performed using the Presto package, which implements a Wilcoxon rank-sum test with area under the receiver operating characteristic curve. Genes with a log fold change ≥ 0.5 and an adjusted P < 0.05 were defined as differentially expressed. Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes enrichment analyses were performed using KOBAS 2.0[7], with significance assessed by hypergeometric testing and Benjamini-Hochberg false discovery rate correction.

Integration with scRNA-seq data (deconvolution)

Cell-type deconvolution of Bin 50 spots was performed using the RCTD algorithm in the spacexr R package (v2.2.1), with scRNA-seq data from GSE163558[8] used as the reference. For each Stereo-seq slice, raw gene-count matrices and spatial coordinates were used to construct a SpatialRNA object, and RCTD was run in doublet mode to account for mixed cell populations within spots. The estimated cell-type proportions were then pruned and renormalized, and the dominant cell type, defined as that with the highest posterior probability, was assigned for downstream visualization.

Cell-cell communication

Intercellular communication was inferred using CellChat (v2.1.2)[9]. A CellChat object was created from the Seurat object using the built-in human protein–protein interaction network as prior knowledge. Communication probabilities between cell types were estimated with computeCommunProb, restricting spatial interactions to within a 250 μm radius to define a local neighborhood and focus the analysis on spatially proximal interactions, while reducing the influence of distant regions. Signaling networks were inferred and aggregated using default parameters, and the visualizations focused on the total number of interactions between cell types and the outgoing signaling pathways from each cell cluster.

Gene set enrichment analysis and gene set variation analysis

Gene set enrichment analysis was performed using clusterProfiler (v4.8.3) with default parameters. Gene set variation analysis (GSVA) was performed on differentially expressed genes using GSVA (v1.46.0), and differential analysis of GSVA scores was conducted using limma (v3.54.2). All analyses were performed in R (v4.3.1) using curated gene sets from the Molecular Signatures Database (MSigDB).

Public dataset-based expression validation

RNA-seq data for TCGA-STAD tumors and GTEx normal gastric tissues were obtained from the UCSC Xena platform and processed using the Toil pipeline[10]. TPM values were log2(TPM + 1)-transformed. Differences in KRT6B expression were analyzed in R (v4.2.1), visualized with ggplot2, and assessed using the Wilcoxon rank-sum test.

Survival analysis

Overall survival and recurrence-free survival were analyzed using gene expression and clinical data from the GEO datasets GSE62254[11], GSE29272[12], GSE15459[13], and GSE14210[14] through the KMplot platform. Patients were stratified according to KRT6B expression level. Survival differences were assessed using the log-rank test, and hazard ratios with 95%CI were obtained.

IHC and IF staining

Paraffin sections were deparaffinized, rehydrated, and subjected to citrate-based antigen retrieval. For both IF and IHC, sections were treated with 30 mL/L H₂O₂ for 15 minutes; for IHC, goat serum blocking was additionally performed. For IF, sections were incubated overnight at 4 °C with anti-KRT6B (Proteintech, United States), followed by repeated antigen retrieval and overnight incubation with anti-EPCAM (Cohesion, China). Alexa Fluor 594-conjugated goat anti-rabbit immunoglobulin G (IgG) (H + L) and Alexa Fluor 488-conjugated goat anti-mouse IgG (H + L) secondary antibodies (both from Sparkjade, China) were applied at 37 °C for 1.5 hours. Sections were counterstained with DAPI (Servicebio, China), mounted with antifade medium, and imaged using a fluorescence scanner. For IHC, sections were incubated overnight at 4 °C with anti-KRT6B (Proteintech, United States), followed by an enhanced enzyme-labeled IgG polymer (ZSGB-Bio, China); signals were developed with DAB, and sections were counterstained and mounted. IHC scores were independently evaluated by two pathologists using the H-score method: H-score = (% weak × 1) + (% moderate × 2) + (% strong × 3).

Cell culture and lentiviral transduction

AGS and SNU-601 human GC cells were purchased from Procell Life Science & Technology Co., Ltd. (Wuhan, Hubei Province, China) and cultured in RPMI-1640 medium (Gibco, Thermo Fisher Scientific, United States) supplemented with 100 mL/L fetal bovine serum and 10 mL/L penicillin-streptomycin at 37 °C in 50 mL/L CO₂.

To establish stable KRT6B-overexpressing cells and empty-vector control cells, AGS and SNU-601 cells were transduced with lentiviruses carrying KRT6B or the corresponding empty vector (HANBIO, China; sequence information is provided in Supplementary material) in the presence of polybrene (5 μg/mL). The MOI was 100 for AGS cells and 10 for SNU-601 cells. The medium was replaced 24 hours after transduction, and transduction efficiency was assessed at 72 hours by green fluorescent protein fluorescence microscopy. The cells were then selected with puromycin (1 μg/mL) for 48 hours, and KRT6B overexpression was confirmed by Western blotting and reverse transcription quantitative polymerase chain reaction.

Western blotting

Total protein was extracted using a Total Protein Extraction Kit (KeyGEN, China) and quantified with a BCA Protein Assay Kit (KeyGEN, China). Equal amounts of protein were separated by SDS-PAGE on 100 g/L gels, transferred to PVDF membranes (Millipore, United States), and incubated overnight at 4 °C with primary antibodies against KRT6B (1:1000, Proteintech, United States) or GAPDH (1:10000, Proteintech, United States), followed by HRP-conjugated secondary antibodies (1:10000, Proteintech, United States). Signals were detected by chemiluminescence and imaged using a Tanon system (Tanon, China).

Cell proliferation assays

Cell proliferation was assessed using CCK-8 and EdU incorporation assays. For the CCK-8 assay, cells were seeded into 96-well plates, and wells containing cell-free medium were used as blanks. At 0 hour, 24 hours, 48 hours, and 72 hours, 10 μL of CCK-8 solution (MCE, United States) was added to each well and incubated at 37 °C for 3 hours. Absorbance was measured at 450 nm, with 620 nm as the reference wavelength, using a microplate reader (BioTek, United States).

For the EdU assay, cell proliferation was further evaluated using an EdU Flow Cytometry Proliferation Assay Kit (Elabscience, China). Cells were incubated with EdU (10 μM) for 1.5 hours, collected, fixed with 40 g/L paraformaldehyde, and permeabilized with 1 mL/L Triton X-100. Fluorescent labeling was performed using the click reaction according to the manufacturer’s instructions, followed by flow cytometric analysis to determine the proportion of EdU-positive cells. A no-EdU sample was used as a negative control for gate setting.

CCK-8 and EdU assays were performed in three independent biological replicates, each with at least three technical replicates. Technical replicate values were averaged within each biological replicate before statistical analysis.

Cell migration and invasion assays

For the wound-healing assay, cells were seeded in 6-well plates and scratched with a 200 μL pipette tip after reaching confluence. After washing with PBS, images were captured at 0 hour and 24 hours, and the wound area was measured using ImageJ software.

For the Transwell invasion assay, Matrigel (Corning, United States) diluted 1:8 was used to coat the upper surface of Transwell inserts (Corning, United States). Cells were resuspended in serum-free medium and seeded into the upper chamber, while medium containing 100 mL/L FBS was added to the lower chamber. After 48 hours, non-invading cells were removed from the upper surface with a cotton swab. Cells on the lower surface were fixed, stained with crystal violet, and imaged under a microscope (× 200).

Wound-healing and Transwell invasion assays were performed in three independent biological replicates, each with at least three technical replicates. Technical replicate values were averaged within each biological replicate before statistical analysis.

Statistical analysis

Statistical analyses were performed using GraphPad Prism 10.6 and R (v4.2.1). Data from cell-based assays are presented as the mean ± SD from three independent biological replicates, each with at least three technical replicates. Technical replicate values were averaged within each biological replicate before statistical analysis. Sample sizes for the in vitro assays were based on common practice for exploratory cell-based experiments and on reproducibility across independent biological replicates. EdU flow cytometry, wound-healing, and Transwell invasion assay data were analyzed by one-way analysis of variance followed by Tukey’s multiple-comparisons test, whereas CCK-8 data were analyzed by repeated-measures two-way analysis of variance followed by Tukey’s multiple-comparisons test at each time point; adjusted P values were reported. Gene expression trend analysis was performed by soft clustering using the Mfuzz R package (v2.66.0)[15]. Comparisons of gene expression levels and IHC scores were performed using Wilcoxon tests (rank-sum for unpaired data and signed-rank for paired data, as appropriate). All tests were two-sided, and P < 0.05 was considered statistically significant.

RESULTS

Cellular composition and spatial architecture of AGC

To characterize the ST features of AGC, tumor tissues were collected from five treatment-naïve patients (P1-P5), and two full-thickness specimens (P2 and P4) were selected for Stereo-seq profiling based on tissue quality and pathological evaluation. Using Bin 50 (50 × 50 DNBs; 25 μm square) as the spot unit, we detected a mean of 1085 transcripts per spot in P2 (median, 887) and 1842 in P4 (median, 1254). After quality control, 94156 and 83588 spots were retained for P2 and P4, respectively.

HE staining showed that tumor cells were mainly located in the mucosa, submucosa, and muscularis propria (yellow boundaries in Figure 1A), with higher gene counts in tumor areas than in adjacent non-tumor areas. Deconvolution using a public GC scRNA-seq reference (GSE163558)[8] identified ten major cell types – epithelial cells, fibroblasts, endothelial cells, T cells, B cells, plasma cells, macrophages, mast cells, neutrophils, and smooth-muscle cells – with representative markers shown in Supplementary Figure 1A and B. Their spatial distributions are shown in Figure 1B and C, and Supplementary Figure 1C.

Figure 1 Cellular composition and spatial architecture in advanced gastric cancer. A: Hematoxylin and eosin-stained full-thickness advanced gastric cancer tissue sections from samples P2 and P4 (left). Spatial maps of detected gene counts (right); B: Spatial mapping of predicted major cell types in P2 and P4; C: Spatial distribution maps of major cell types in P2 and P4. Color intensity indicates the deconvolution-predicted proportion of each cell type in each spot.

In sample P2, tumor epithelial cells (TECs) were widely distributed across the mucosa, submucosa, and muscularis propria. In P4, TECs were concentrated in the mucosa and submucosa, with a distinct oval tumor nest in the muscularis propria. Fibroblasts were relatively more abundant in P4 and were enriched at the tumor periphery, whereas in P2, they formed linear strands between tumor compartments.

Immune composition also differed markedly between the two samples. P2 contained more immune cells, with T and B cells forming aggregates resembling tertiary lymphoid structure (TLS)-like clusters, consistent with follicle-like patterns on HE staining. Notably, these T/B-cell clusters expressed several naïve immune-associated genes at elevated levels (e.g., IL7R, TCF7, PAX5, and IGHM), suggesting that a subset of lymphocytes may be in naïve or early differentiation states. In contrast, P4 showed lower immune-cell abundance, with sparsely distributed plasma cells and relatively few other immune subsets.

Spatial niches identification and characterization using BANKSY

To further investigate functional properties and spatial interactions across regions, we jointly clustered spots from P2 and P4 using BANKSY[6], which integrates the transcriptome of each spot with those of its neighbors to delineate “spatial niches” that combine spatial continuity with transcriptomic similarity.

We identified 13 niches (Figure 2A), and their distributions across tissue sections are shown in Figure 2B; marker genes and GO enrichment are shown in Supplementary Figure 2A and B. In P2, niches N0, N2, N5, and N9 were predominantly assigned to TECs (each > 95%), with N9 containing relatively few spots, and TECs were also highly represented in N1, N3, and N4. In P4, TECs were enriched in N2, N4, N5, and N9 (Figure 2C). Supplementary Table 3 summarizes the dominant cell type and spot counts for each niche. Notably, in P4, the spatial location of N9 closely overlapped with the well-demarcated tumor nest in the muscularis propria, indicating a spatially and functionally distinct tumor subregion.

Figure 2 Spatial niches identification and characterization using BANKSY. A: UMAP of integrated P2 and P4 spatial transcriptomic data. Left: Spots colored by sample (P2, red; P4, blue). Right: Spots colored by BANKSY-defined niche identities (N0-N12); B: Spatial maps of BANKSY-defined niches (N0-N12) in P2 (left) and P4 (right); C: Bar plots showing the number of spots assigned to each first-weight cell type within each niche in P2 (top) and P4 (bottom). Colors denote different cell types; D: Spatial localization of each niche in P2 (top) and P4 (bottom), with red indicating spatial regions assigned to each niche.

Mapping these niches back onto the tissue sections (Figure 2D) showed that in P2, N5 was mainly located in the mucosa, N4 in the muscularis mucosae and submucosa, N0 in the submucosa, and N2 spanned from the submucosa to the muscularis propria, forming a superficial-to-deep spatial arrangement. Although N0 contained the largest number of TEC-assigned spots, tumor-associated marker expression and tumor-cell weights per spot were lower than those in N5 and N2 (Supplementary Figure 2C), suggesting weaker tumor-associated transcriptional signals. A similar layer-wise pattern was observed in P4, with N5 in the mucosa, N4 closer to the submucosa, and N2 at the submucosal boundary. Together, these data reveal multiple spatially defined niches and highlight the complex organization of the AGC tumor microenvironment (TME). We therefore focused subsequent analyses on three tumor-associated niches – N5, N4, and N2 – which showed relatively consistent spatial positions and transcriptional features in both P2 and P4.

Spatial functional profiling of distinct niches in AGC

Spatial mapping showed a continuous arrangement of niches N5, N4, and N2 from the mucosa to the submucosa and muscularis propria (Figure 3A). All three niches showed high expression of canonical tumor epithelial markers (Figure 3B), suggesting that they represent distinct tumor-associated niche states along the depth axis. To explore their transcriptional dynamics, we applied Mfuzz. This fuzzy clustering algorithm allows genes to have graded membership in multiple clusters, and identified four major expression trends across N5, N4, and N2 (Figure 3C): (1) Trend 1 (increase then decrease); (2) Trend 2 (decrease then increase); (3) Trend 3 (continuous increase); and (4) Trend 4 (continuous decrease). Go enrichment results for Trends 1 and 2 are shown in Supplementary Figure 3. Subsequent analyses focused on Trends 3 and 4.

Figure 3 Spatial functional profiling of distinct niches in advanced gastric cancer. A: Spatial distribution of niches N5, N4, and N2 in P2 (red shading with outlined boundaries), showing a continuous distribution from the mucosa to the submucosa and muscularis propria; B: Violin plots showing expression of canonical tumor epithelial markers (EPCAM, KRT18, KRT8, and KRT19) across niches in P2 (red) and P4 (blue); C: Four gene-expression trends identified by Mfuzz based on expression profiles across N5, N4, and N2; D and E: Top enriched Gene Ontology biological process terms for Trend 3 and Trend 4. Dot size represents the number of input genes, and color indicates the adjusted P value; F: Violin plots showing expression of immune-related genes (CD8A, CD4, and CXCL9) in N2, N4, and N5; G: Violin plots showing expression of stemness-associated genes (CD44, ALDH1A1, PROM1, LGR5, and SOX9) in N2, N4, and N5. Statistical significance: ^aP < 0.05, ^bP < 0.01, ^cP < 0.001.

GO enrichment analysis showed that Trend 3 genes were mainly involved in protein synthesis and translation, antigen processing and presentation, particularly major histocompatibility complex class II-related pathways, and T-cell activation and response (Figure 3D). In contrast, Trend 4 genes were enriched in pathways related to cell motility and migration, anti-apoptotic processes, and membrane protein localization (Figure 3E). Trend 3 programs increased progressively from N5 to N2, indicating increased biosynthetic and antigen-presentation-related transcriptional activity with depth, together with a marked increase in adaptive immune-related signatures. Consistently, key immune-related genes (CD8A, CD4, and CXCL9) displayed a spatial gradient, with significantly lower expression in N5 than in N4 and N2 (Figure 3F).

Trend 4 programs decreased from N5 to N2 along the mucosa-to-deep axis, suggesting reduced motility-associated and invasion-associated transcriptional signatures in deeper regions. Transcriptomic profiling further indicated that TECs within N5 exhibited stemness-associated transcriptional features (Figure 3G). Together, these findings indicate depth-associated differences in transcriptional programs across tumor niches, with superficial regions, such as N5, showing stronger stemness-related and motility-related features, and deeper niches, such as N2, showing increased biosynthetic and immune-related transcriptional activity.

Niche-specific intercellular communication in sample P2

To explore mechanisms that might underlie niche-specific differences in stemness-related and invasion-related features across N5, N4, and N2, we analyzed ligand-receptor (L-R) interactions in P2, which contained more spots in these niches and showed clearer spatial stratification, providing a more robust and interpretable context for niche-level interaction inference. We focused on interactions in which TECs acted as receivers and restricted spatial interactions to within a 250 μm radius.

In N5, TEC-to-TEC communication was strongest (Figure 4A). The inferred network included multiple L-R pairs previously implicated in proliferation, invasion, immune evasion, and treatment resistance, including MDK-related axes (MDK-NCL, MDK-LRP1, MDK-SDC4, and MDK-ITGA6/ITGB4) and MIF-CD74/CD44 (Figure 4B and C). Fibroblast-derived NRG1 also signaled to TECs through NRG1-ERBB3, NRG1-ERBB2/ERBB3, and NRG1-ITGA6/ITGB4 axes, which have been reported to promote tumor stemness, migration, and therapy resistance. Together, these interactions suggest that TECs in N5 display communication features consistent with a more invasive, proliferative, treatment-resistance-associated state, with stemness-related characteristics.

Figure 4 Niche-specific intercellular communication in sample P2. A: Heatmaps showing the overall interaction strengths among major cell types within N5; B: Bubble plots showing significant inferred ligand-receptor interactions targeting tumor epithelial cells (TECs) in N5. Dot size indicates the significance level, and color indicates communication probability; C: Spatial interaction networks of the MDK signaling pathway in N5. The background indicates the spatial distribution of major cell types within each niche; nodes denote cell types, and edges indicate inferred ligand-receptor interactions; D: Heatmaps showing the overall interaction strengths among major cell types within N4; E: Bubble plots showing significant inferred ligand-receptor interactions targeting TECs in N4. Dot size indicates the significance level, and color indicates communication probability; F: Spatial interaction networks of the MDK signaling pathway in N4. The background indicates the spatial distribution of major cell types within each niche; nodes denote cell types, and edges indicate inferred ligand-receptor interactions; G: Heatmaps showing the overall interaction strengths among major cell types within N2; H: Bubble plots showing significant inferred ligand-receptor interactions targeting TECs in N2. Dot size indicates the significance level, and color indicates communication probability; I: Spatial interaction networks of the MDK signaling pathway in N2. The background indicates the spatial distribution of major cell types within each niche; nodes denote cell types, and edges indicate inferred ligand-receptor interactions.

In N4, fibroblast-TEC communication was more prominent (Figure 4D), which may be related to the higher fibroblast proportion in this niche. MDK-related signaling in TECs persisted but was reduced compared with N5 (Figure 4E and F). Fibroblast-TEC interactions included pro-migratory axes such as IGF1-IGF1R, IGF1-ITGA6/ITGB4, and FGF7-FGFR2, although their inferred activity was relatively weak. In contrast, the IGFBP3-TMEM219 axis showed relatively high inferred activity in N4 (Figure 4E and Supplementary Figure 4), a pathway reported to have tumor-suppressive and anti-metastatic roles in multiple cancer contexts, consistent with weaker pro-invasive signaling than that in N5.

N2 showed a similar pattern, with reduced activity of MDK-related pairs and increased inferred activity of IGFBP3-TMEM219 (Figure 4G-I and Supplementary Figure 4F). Overall, N5 was characterized by signaling features associated with stemness, invasion, proliferation, and treatment resistance, whereas N4 and N2 showed weaker malignant-associated signaling and relatively increased signaling axes with reported suppressive roles. These communication patterns may help explain the spatial differences in stemness-related and invasion-related features across the three niches.

Molecular and functional profiling of a KRT6B-enriched tumor nest (N9) in P4

In sample P4, we identified a well-demarcated, spatially isolated tumor nest, clearly separated from other tumor regions. BANKSY clustering assigned this structure to a distinct niche, N9 (Figure 5A). To functionally characterize N9, we performed GSVA across tumor-associated niches (N9, N5, N4, and N2). N9 showed a distinct transcriptional profile, with enrichment of cell-cycle and proliferation programs (e.g., E2F_TARGETS and G2M_CHECKPOINT), metabolic pathways (e.g., GLYCOLYSIS and FATTY_ACID_METABOLISM), growth-related signaling (e.g., MTORC1 and PI3K_AKT_MTOR), and immune-response and inflammatory-response signatures (e.g., INTERFERON_GAMMA_RESPONSE and COMPLEMENT) (Figure 5B). Gene set enrichment analysis comparing N9 with the other niches showed concordant enrichment patterns (Figure 5C), indicating that N9 represents a niche with transcriptional features of high proliferation, metabolic activity, and immune/inflammatory enrichment.

Figure 5 Molecular and functional profiling of a KRT6B-enriched tumor nest (N9) in P4. A: Spatial location of N9 (red with outlined boundary) in P4; B: Gene set variation analysis heatmap of tumor-associated niches (N4, N9, N2, and N5), showing enrichment patterns in N9, including cell-cycle/proliferation, metabolic, and immune/inflammatory programs. Red indicates relative upregulation and blue indicates relative downregulation, and color intensity reflects row-scaled enrichment scores; C: Gene set enrichment analysis comparing N9 with N5/N4/N2, showing significantly upregulated (left) and downregulated (right) gene sets in N9, including proliferation-related, metabolism-related, immune-response-related, and hypoxia-related pathways. The X-axis represents GeneRatio, the dot size indicates the number of enriched genes, and the color denotes adjusted P values; D: Volcano plot of differentially expressed genes between N9 and N5, with representative upregulated genes including MMP7, LCN2, UBD, and KRT6B. Red, significantly upregulated genes; blue, significantly downregulated genes; gray, non-significant genes; E: Gene Ontology biological process enrichment of genes upregulated in N9. Dot size represents the number of enriched genes, and color indicates adjusted P values; F: Spatial map of dominant cell types in the region surrounding the N9 niche, showing that fibroblasts are enriched along the left and upper margins of N9; G: Bubble plots of significant inferred ligand-receptor interactions between N9 tumor epithelial cells and fibroblasts. Left: N9 to fibroblasts; right: Fibroblasts to N9. Dot size indicates statistical significance (small: 0.01 < P < 0.05; large: P < 0.01), and color denotes communication probability; H: Spatial heatmap of KRT6B expression in P4, showing high-expression regions overlapping the N9 tumor nest. Color scale indicates relative expression intensity (low to high); I: Immunofluorescence staining of gastric cancer tissue sections: Nuclei (DAPI, blue), epithelial cells (EPCAM, green), and KRT6B (red). In the panoramic view (top), KRT6B signal overlaps with EPCAM in tumor nests and at the invasive front (scale bar, 500 μm). In the magnified view (bottom), a representative tumor nest shows pronounced KRT6B signal overlapping EPCAM, with little to no detectable KRT6B in surrounding stroma (scale bar, 100 μm).

Given that N5 exhibited relatively stronger stemness-related and motility-related transcriptional features, we used it as the reference niche for differential expression analysis between N9 and N5. The results showed that genes such as MMP7, LCN2, and UBD were significantly upregulated in N9 (Figure 5D), and these genes have previously been implicated in tumor invasion and metastasis. GO enrichment analysis revealed that the upregulated genes were mainly involved in immune and inflammatory regulation, TME remodeling, and cytoskeletal organization (Figure 5E). Spatial analysis further showed that N9 was partially surrounded by a semicircular fibroblast-rich region, particularly along its left and upper boundaries (Figure 5F). Consistently, fibroblast-associated marker genes were enriched around the N9 periphery (Supplementary Figure 5A). To explore potential crosstalk between N9 and its surrounding stroma, we extracted spots from N9 and the adjacent fibroblast compartment and constructed a bidirectional L-R interaction network. Predicted signaling from N9 to fibroblasts included L-R axes associated with fibroblast activation and transition toward a cancer-associated fibroblast-like phenotype, whereas signaling from fibroblasts to N9 involved pathways linked to tumor progression, migration, and epithelial plasticity (Figure 5G). In addition, the fibroblasts surrounding N9 showed high expression of extracellular matrix-related genes together with relatively low ACTA2 and MYL9 levels (Supplementary Figure 5B), consistent with matrix cancer-associated fibroblast (mCAF)-like characteristics. Together, these findings suggest that N9 is embedded in a distinct fibroblast-rich stromal context.

Notably, KRT6B was preferentially enriched in N9, and spatial heatmaps showed relative enrichment along the periphery of the tumor nest (Figure 5H). IF staining of paraffin-embedded GC sections showed strong co-localization of KRT6B with the epithelial marker EPCAM, with markedly enhanced signals in tumor nests and at the invasive front, and little to no detectable KRT6B signal in the surrounding stroma (Figure 5I). The preferential enrichment of KRT6B in N9 prompted further evaluation of its clinical relevance and functional role.

In summary, niche N9 is spatially and structurally isolated and characterized by transcriptional features related to proliferation, metabolism, invasion, immunity, and inflammation. The preferential enrichment of KRT6B may be associated with the distinct stromal context, morphology, and molecular features of this tumor microregion.

KRT6B is upregulated in GC, correlates with poor prognosis, and promotes malignant phenotypes in vitro

To assess KRT6B expression in GC, we performed immunohistochemical (IHC) staining on tumor tissues and matched adjacent non-tumor tissues from 54 patients with GC. KRT6B protein expression was significantly higher in GC tissues than in matched adjacent tissues (Figure 6A and B). Consistently, analysis of UCSC Xena-integrated TCGA-STAD and GTEx transcriptomic data showed that KRT6B mRNA was significantly upregulated in GC tissues (Figure 6C). Kaplan-Meier analysis further showed that high KRT6B expression was associated with shorter overall survival (hazard ratio = 1.38, 95%CI: 1.16-1.64; log-rank P = 0.00027; Figure 6D). Consistently, high KRT6B expression was also associated with shorter recurrence-free survival (Supplementary Figure 6).

Figure 6 KRT6B is upregulated in gastric cancer, correlates with poor prognosis, and promotes malignant phenotypes in vitro. A: Representative immunohistochemical staining images showing KRT6B expression in gastric cancer tissues and matched adjacent non-tumor tissues (scale bar, 200 μm); B: Quantification of KRT6B immunohistochemistry scores (H-score) in paired tumor and adjacent non-tumor tissues; C: Analysis of UCSC Xena-integrated TCGA-STAD and GTEx transcriptomic data showing significantly increased KRT6B mRNA expression in gastric cancer tissues compared with normal gastric tissues; D: Kaplan-Meier overall survival analysis of patients stratified by KRT6B expression level; E: Western blot validation of stable KRT6B overexpression in AGS and SNU-601 cells (GAPDH as a loading control); F and G: CCK-8 assays showing increased proliferative capacity in KRT6B-overexpressing AGS (F) and SNU-601 (G) cells; H and I: Representative flow cytometry histograms of EdU incorporation in AGS (H) and SNU-601 (I) cells. A no-EdU control was used to define the APC-positive gate; J and K: Quantification of EdU-positive cell proportions in AGS (J) and SNU-601 (K) cells; L and M: Representative Transwell invasion assay images of AGS (L) and SNU-601 (M) cells; N and O: Quantification of invading cell numbers in AGS (N) and SNU-601 (O) cells; P and Q: Representative wound-healing assay images (0 hour and 24 hours) of AGS (P) and SNU-601 (Q) cells; R and S: Quantification of wound closure rates in AGS (R) and SNU-601 (S) cells. Data are presented as mean ± SD from independent experiments. NS: Not significant; ^aP < 0.05, ^bP < 0.01, ^cP < 0.001.

To investigate the biological role of KRT6B in GC cells, we selected the AGS and SNU-601 GC cell lines, which exhibit very low endogenous KRT6B expression, and established stable KRT6B-overexpressing cell lines along with corresponding empty-vector controls. Western blotting confirmed successful KRT6B overexpression (Figure 6E). Functional assays showed that KRT6B overexpression enhanced malignant phenotypes in GC cells. Specifically, CCK-8 assays and EdU incorporation combined with flow cytometric analysis consistently showed that KRT6B overexpression significantly increased the proliferative capacity of AGS and SNU-601 cells (Figure 6F-K). In addition, Transwell invasion assays demonstrated that KRT6B overexpression markedly increased the number of invading cells in both cell lines (Figure 6L-O). Consistent with these findings, wound-healing assays showed that KRT6B overexpression increased the wound-closure rate of AGS and SNU-601 cells, suggesting enhanced migratory capacity (Figure 6P-S). Taken together, KRT6B is highly expressed in GC tissues and is associated with poor prognosis; in vitro, KRT6B overexpression promotes proliferation-related, invasion-related, and migration-related phenotypes in GC cells.

DISCUSSION

In this study, we used Stereo-seq to generate full-thickness ST profiles from two AGC samples and to characterize intratumoral spatial heterogeneity across anatomical layers. Although both tumors were of an advanced stage, their immune microenvironments differed markedly: P4 showed limited immune infiltration, whereas P2 contained abundant immune cells, including multiple TLS-like aggregates. TLSs are organized lymphoid structures associated with chronic inflammation and have been linked to enhanced anti-tumor immune responses in several cancers[16]. We also observed elevated expression of genes associated with early lymphocyte states in subsets of immune cells infiltrating P2, suggesting spatial heterogeneity in the differentiation states of lymphocyte-rich aggregates, including TLS-like aggregates. Prior work in GC has shown that TLSs at different locations (tumor core, invasive front, and peritumoral region) can display distinct functional states that influence prognosis and response to immunotherapy[17], supporting the idea that TLS functionality may be influenced by spatial context. Collectively, these findings suggest that, even among AGCs of similar stage, immune spatial organization may vary substantially. These differences in immune context also informed our subsequent communication analysis. Together with the greater number of spots in N5, N4, and N2 and the clearer spatial stratification in P2, P2 was the most suitable sample for niche-specific CellChat analysis in the present study.

Spatial functional heterogeneity has been described in multiple solid tumors, including oral squamous cell carcinoma, kidney cancer, and glioblastoma, in which region-dependent or depth-dependent differences in transcriptional programs, stemness, metabolism, and immune activity have been reported[18-20]. To dissect such patterns in AGC, we used BANKSY, an ST-specific clustering algorithm that integrates the transcriptome of each spot with that of its neighborhood to identify spatially coherent niches with similar transcriptional programs[6]. This analysis resolved 13 spatial niches, including multiple TEC-dominant niches with distinct transcriptional profiles. Consistent with scRNA-seq studies showing that GC TECs comprise multiple molecularly and functionally distinct subpopulations[21], our study places this heterogeneity in spatial context and further suggests that intratumoral diversity is associated with spatial organization and local niche context.

We further focused on three tumor-dominated niches (N5, N4, and N2) arranged along a superficial-to-deep axis, using this framework to examine depth-associated spatial heterogeneity in GC. Mfuzz-based analysis revealed a spatial gradient in tumor-associated transcriptional states across these niches. Genes that were progressively upregulated from N5 to N2 were enriched in antigen presentation, immune activation, and metabolic/biosynthetic pathways, consistent with deeper local immune and metabolic transcriptional programs. In contrast, genes that decreased from N5 to N2 were associated with motility-related and adhesion-related pathways, suggesting weaker motility-associated and invasion-associated signatures in deeper niches. TECs in N5 also displayed stemness-associated transcriptional features, suggesting a relatively more plastic transcriptional state than TECs in N4 and N2. Together, these findings support an association between tumor-cell transcriptional states and spatial location. We speculate that these depth-associated shifts may be shaped by local gradients in oxygen, nutrients, pH, and immune-cell abundance within the TME.

Cell-cell communication analysis further suggested a potential signaling framework underlying this spatial gradient in stemness-related and invasion-related features. In N5, MDK-related signaling showed high inferred activity in TECs. MDK is frequently upregulated in cancers and has been reported to promote proliferation, invasion, immune evasion, and resistance[22-24]. MIF-CD74/CD44 signaling, which has been reported to activate MAPK and AKT pathways and to support invasion, survival, and angiogenesis[25], was also prominent. Fibroblast-derived NRG1 signaling to TECs via ERBB3 and integrin α6β4 has been linked to stemness, proliferation, invasion, and drug resistance[26,27]. By contrast, the inferred activities of these pro-tumor signaling axes were reduced in N4 and N2, where IGFBP3-TMEM219 signaling became more evident. This axis has been reported to induce apoptosis and inhibit metastasis in several cancer contexts[28]. Together, these observations suggest that intercellular signaling networks vary with spatial context and may contribute to the spatial differences in stemness-related and invasion-related features observed in AGC.

In sample P4, we identified a spatially isolated tumor nest that BANKSY classified as a distinct niche (N9). Compared with other tumor-associated niches, N9 showed relatively stronger enrichment of proliferation-related, metabolism-related, and immune/inflammatory response-related programs, together with upregulation of genes such as MMP7, LCN2, and UBD, which have been widely implicated in tumor invasion and metastasis[29-31]. Taken together, these findings suggest that N9 represents a spatially distinct tumor microregion with stronger transcriptional features related to proliferation, metabolism, and invasion. N9 was embedded in a fibroblast-rich stromal compartment and showed predicted bidirectional communication with adjacent fibroblasts. We further analyzed the transcriptomic features of these surrounding fibroblasts and found high expression of multiple extracellular matrix-related components, including COL1A1, COL1A2, and COL3A1, along with relatively low expression of ACTA2 and MYL9. This molecular profile is consistent with mCAF-like features. Previous studies in skin cancer have also shown that mCAFs can localize around tumor nests, form a physical barrier, and may help shield tumor nests while limiting T-cell infiltration[32]. In this context, the biological significance of N9 may need to be interpreted not only in terms of tumor-intrinsic transcriptional programs but also in relation to local tumor-stroma interactions within this niche.

KRT6B was highly expressed and preferentially enriched in N9. KRT6B encodes keratin 6B, a type II intermediate filament protein involved in maintaining cytoskeletal integrity[33]. Studies of KRT6B in cancer remain limited overall. A small number of reports suggest that KRT6B may participate in bladder cancer cell migration and invasion and may also be associated with chemoresistance in GC[34,35]. In addition, several bioinformatics-based studies have linked KRT6B expression to poor prognosis[36,37]. In our study, KRT6B was highly expressed in GC tissues and was associated with poor prognosis, and in vitro functional assays further suggested that KRT6B overexpression promotes GC cell proliferation, migration, and invasion. However, the precise molecular mechanisms by which KRT6B contributes to malignant phenotypes in GC remain unclear. Its upstream regulatory mechanisms, downstream effector molecules, and related signaling networks warrant further investigation.

This study has several limitations. First, the limited sample size and relative homogeneity of the Stereo-seq cohort restrict the generalizability of our findings, and larger cohorts and/or integration with public spatial datasets are needed to better characterize spatial heterogeneity in AGC. In addition, larger cohorts will be needed to evaluate the associations of KRT6B expression and key spatial niche signatures with clinical outcomes. Second, the depth-associated differences across N5, N4, and N2 were inferred primarily from ST analyses and require further orthogonal validation and functional investigation. Third, although Stereo-seq provides high spatial resolution, its transcript-capture efficiency is lower than that of scRNA-seq, which may limit the detection of low-abundance genes and rare cell populations. Fourth, ST computational methods are still evolving and may not fully capture multimodal spatial information, highlighting the need for improved integration with single-cell data. Finally, further in vitro and in vivo studies are required to clarify the biological roles of KRT6B and to evaluate its potential as a biomarker or therapeutic target in GC.

CONCLUSION

In this study, we applied Stereo-seq to perform full-thickness ST analysis of the gastric wall in AGC, delineating tumor niche features across distinct spatial locations. We found spatial differences in tumor-cell states, immune-related transcriptional features, and tumor-stroma interactions, including depth-associated gradients in stemness-related, migration-related, immune-related, and metabolism-related programs, as well as niche-specific intercellular communication patterns. We also identified a spatially isolated, KRT6B-enriched tumor nest (N9) with relatively stronger transcriptional features for proliferation, metabolism, and immune/inflammatory processes. Further validation showed that KRT6B is upregulated in GC and associated with poor prognosis. In vitro assays suggested that KRT6B promotes proliferation-related, migration-related, and invasion-related phenotypes in GC cells. Overall, this study provides a spatial framework for understanding AGC heterogeneity and suggests that KRT6B warrants further mechanistic and translational investigation. Future studies should evaluate the role of KRT6B in tumor growth and metastasis in vivo and further define its upstream regulatory mechanisms and downstream effector pathways in AGC.