Single-cell and spatial transcriptome profiling identifies the immunosuppressive spatial niche in -mutant colorectal cancer.

Introduction
Colorectal cancer (CRC) is one of the most prevalent malignancies globally.1
KRAS mutations, which are among the most common driver gene mutations in CRC, occur in approximately 45% of patients and are closely associated with aggressive tumor growth, metastatic potential, and poor prognosis.2 3
KRAS mutations have been identified as critical molecular markers of the refractory CRC, particularly leading to primary resistance to epidermal growth factor receptor (EGFR)-targeted therapies, such as cetuximab and panitumumab, thereby preventing mutant-type patients from benefiting from these treatments.4 Additionally, KRAS-mutant tumors often develop secondary resistance through clonal evolution or compensatory activation of signaling pathways (eg, reactivation of the EGFR pathway), which further restricts clinical efficacy.5 6 Although research has indicated that KRAS mutations may drive resistance through abnormal activation of signaling axes such as RAS/ERK and PI3K/AKT, the precise molecular mechanisms underlying immune microenvironment suppression and metabolic reprogramming remain incompletely understood.79 These unresolved mechanisms present significant challenges for the development of targeted therapies.
The combined application of single-cell RNA sequencing (scRNA-seq) and spatial transcriptomics (ST) has provided revolutionary tools for elucidating the spatial heterogeneity of the tumor microenvironment (TME) in CRC.1012 Conventional bulk transcriptomics loses the spatial localization information of cells during tissue dissociation. In contrast, scRNA-seq reveals gene expression profile differences among tumor cells, immune cells, and stromal cells at a single-cell resolution.13 14 However, tumor resistance and progression often depend on the intercellular spatial interaction networks. For example, Qi et al revealed the co-localization mechanism of FAP+fibroblasts and SPP1+macrophages in the resistant CRC microenvironment.15 Furthermore, Feng et al demonstrated that dMMR (deficient mismatch repair) patients with CRC respond more favorably to immune checkpoint blockade therapy, primarily due to the interaction between LAMP3+dendritic cells (DCs) and CXCL13+T cells, which facilitates the formation of an effective tumor-stroma boundary. In contrast, pMMR (proficient mismatch repair) tumors show poor therapeutic responses due to barriers formed by CXCL14+cancer-associated fibroblasts.16 In summary, ST offers significant advantages for studying the TME and cellular spatial interaction networks.
However, the ST of KRAS-mutant CRC has not been previously explored. This study is the first to perform ST and analysis on KRAS-mutant CRC to elucidate the spatial heterogeneity. By combining scRNA-seq with ST, we identified an immunosuppressive niche within KRAS-mutant CRC, which may offer new perspectives and potential targets for the treatment of KRAS-mutant tumors. The flowchart is created using BioRender (figure 1A).

Material and methods

Patient specimens
Tumor tissues were collected from four patients who underwent surgical resection for CRC at the Affiliated Suzhou Hospital of Nanjing Medical University. Clinical information was collected at the Affiliated Suzhou Hospital of Nanjing Medical University (online supplemental table S1).

Sample preparation and processing for 10x Visium ST
Fresh tumor tissues were first snap-frozen in isopentane and subsequently embedded in OCT compound. Cryosections were prepared using a cryostat, maintaining the frozen state, and sectioned to fit the dimensions of the 10x Genomics Visium spatial gene expression slides. Sections were placed onto both the Visium Spatial Gene Expression Slide and the Tissue Optimization Slide. H&E staining and imaging were performed to assess tissue morphology.
Sample fixation and imaging were carried out prior to permeabilization. The optimal permeabilization time was determined using the tissue optimization protocol. First-strand complementary DNA (cDNA) synthesis was initiated via reverse transcription, followed by second-strand synthesis through PCR. Subsequent denaturation separated the cDNA from the slide. The spatially barcoded, full-length cDNA was then amplified by PCR to obtain sufficient material for library preparation. Enzymatic fragmentation and size selection were used to optimize the cDNA amplicon length. During library construction, sample indexes (P5, P7, i7, i5) and the TruSeq Read 2 primer sequence were incorporated through end-repair, A-tailing, adaptor ligation, and PCR amplification. Final libraries were compatible with standard Illumina amplification protocols and contained both P5 and P7 primers.

ST sequencing and data processing
Libraries were sequenced on the Illumina NovaSeq 6000 platform following the 10x Genomics Visium protocol.17 Each library contained standard paired-end constructs flanked by P5 and P7 adapters. The 16 bp spatial barcode and 10 bp Unique Molecular Identifier (UMI) were encoded in Read 1, while Read 2 captured the cDNA insert. The i7 index read recorded the sample index.
Demultiplexing, alignment, UMI counting, and gene quantification were performed using the Space Ranger software suite (V.3.1.1). The mkfastq module was used to convert raw BCL files into FASTQ format and filter low-quality reads. The count module aligned reads to the human reference genome (GRCh38), using GENCODE V.32 for gene annotation.18 The resulting expression matrix included UMI counts and gene expression values for each spatially barcoded spot on the slide.
ST data were processed using Seurat (V.4.3.0) with stringent quality control. Spots containing fewer than 200 detected genes were excluded from subsequent analysis. Gene expression normalization was performed using SCTransform, and 2,000 highly variable genes were selected for downstream analyses. To address technical batch effects, Harmony (V.0.1.1) was implemented using sample ID as the batch covariate, with the first 20 principal components serving as input for integration.

Preprocessing of scRNA-seq data
We integrated publicly available scRNA-seq datasets from CRC tissues (online supplemental table S2), including GSE132465, GSE144735, GSE161277, GSE200997, GSE221575, and GSE231559, to construct a comprehensive reference dataset.1923 In total, we analyzed tumor samples from 27 KRAS-wild-type and 27 KRAS-mutant patients.
Data preprocessing was performed using the Scanpy (V.1.9.3) pipeline.24 Scrublet (V.0.2.3) was applied to identify and remove potential doublets.25 Cells were filtered by excluding those with fewer than 200 detected genes or with over 50% mitochondrial or ribosomal gene content.26 The remaining cells were normalized and scaled27 using sc.pp.normalize_total() and sc.pp.scale(), and mitochondrial and ribosomal effects were regressed out using sc.pp.regress_out(). After quality control, 5,000 highly variable genes were selected using sc.pp.highly_variable_genes(). Batch effect correction and data integration were performed using Harmony (V.0.0.9), with sample ID as the batch key.28 Dimensionality reduction was carried out using principal component analysis (PCA).
Cell clustering was performed using the Leiden algorithm (resolution 0.2–1.2), and the results were visualized with Uniform Manifold Approximation and Projection (UMAP). Each cluster was manually annotated based on established marker genes reported in previous literature. Subsequently, we performed subclustering of the initial clusters, ultimately identifying 34 distinct cell subpopulations. Marker genes specific to each subcluster were identified using the FindAllMarkers function and used for subpopulation annotation.29
The single-cell sequencing analysis pipeline for mice was similar to that for humans.

Cell–cell communication analysis
We employed the CellChat package (V.1.6.1) to infer direct intercellular interactions and compare cell–cell communication networks between KRAS wild-type and KRAS-mutant samples.30 Following the official workflow, we created a CellChat object using the gene expression matrix of all 34 identified subclusters. The data were normalized, and potential signaling pathways were identified using the built-in ligand-receptor interaction database.
Subsequently, we constructed the intercellular communication network and used the computeCommunProb function to calculate communication probabilities and signal strengths between subclusters. To visualize the results, we generated ligand-receptor interaction bubble plots using netVisual_bubble().

Single-sample Gene Set Enrichment Analysis
To evaluate the enrichment scores of the Epi_01, Fib_CTHRC1, and Mono_S100A8 cellular signatures in bulk transcriptomic data, we performed single-sample Gene Set Enrichment Analysis (ssGSEA).31 The analysis was implemented using the GSVA R package (V.1.44.4). Signature genes for each cell subpopulation were defined as the top 30 differentially expressed genes identified from our scRNA-seq analysis using the FindAllMarkers function.
For prognostic and clinical correlation analyses, ssGSEA scores were calculated for each tumor sample in the TCGA-COADREAD cohort. Patients were stratified into high and low expression groups based on median ssGSEA scores for each cellular signature. Associations between categorical clinical variables and cellular signature expression groups were analyzed using three-way contingency tables and assessed by χ2 tests. For drug sensitivity analysis, ssGSEA scores were computed using the sample of the GSE200407 dataset. Differences in ssGSEA scores between chemotherapy-sensitive and chemotherapy-insensitive groups were assessed using two-sided Wilcoxon rank-sum tests.

InferCNV analysis
To infer large-scale chromosomal alterations, we used InferCNVpy (V.0.4.2) to compute copy number variation (CNV) scores at the single-spot level. Immune cell-enriched spots served as the reference population, and mitochondrial chromosome interference was excluded from analysis.32 CNV gains were visualized in red and losses in blue, with the intensity of the color representing the magnitude of variation.

Cell-type deconvolution
To resolve the spatial distribution of cell types in the TME, we performed deconvolution of the ST data using the Cell2location algorithm (V.0.1.4). This probabilistic model integrates single-cell reference data with spatial gene expression to estimate the abundance of each cell type per spot.33
Specifically, 27 samples with clear KRAS status and 34 subpopulations were selected from public datasets to construct a reference profile. Cell2location modeled gene expression at each spatial spot as a linear combination of the cell-type-specific signatures, inferring the relative contribution of each cell type. The resulting cell type abundance scores were normalized to sum to one per spot, facilitating direct comparison across tissue regions.

Multimodal intersection analysis
To investigate gene-level correspondence between modalities, we compared differentially expressed gene sets identified from both ST (region-specific niches) and scRNA-seq (cell-type-specific clusters). Differentially expressed genes were identified using the FindAllMarkers() function with thresholds of avg_logFC>0.1 and adjusted p<0.01.34

Identification of cellular ecotypes in spatial data
To delineate spatially organized cellular communities within the TME, we used the ISCHIA algorithm (V.1.0.0.0), which identifies spatial co-occurrence patterns in both healthy and inflamed tissues.18 The optimal number of clusters (k) was selected as six, based on the elbow point derived from the mean variance across clustering results for k values ranging from 2 to 20.
Before clustering, PCA was employed to reduce data dimensionality. The Composition.cluster() function was then applied to the top six principal components, resulting in the classification of spatial spots into 10 distinct compositional groups. These groupings, referred to as cellular ecotypes, represent distinct spatial configurations of co-localized cell types, reflecting the principle that biologically interconnected cells frequently occupy similar niches within the tumor landscape.

Spatial mapping of cell dependencies
To investigate the spatial interdependencies among various cell types within the CRC TME, we applied MISTy, a multiview machine learning framework designed to infer how the abundance of one cell type spatially predicts the presence or behavior of others.35 We evaluated both intrinsic (within-spot) and extrinsic (neighboring spots, or “juxta view”) relationships by modeling co-localization patterns. For the juxta view, we aggregated the cell type abundances from immediate spatial neighbors of each spot to model short-range intercellular interactions. MISTy was executed using default parameters, and the resulting predictive importance scores were visualized as heatmaps.

Homotypic and heterotypic score calculation
To assess spatial aggregation based on cell identity, we calculated homotypic and heterotypic interaction scores. Homotypic scores were computed using the kNN() function from the dbscan package (V.1.2.0), by comparing the mean degree of adjacency networks derived from the actual data to those from randomly permuted networks. Heterotypic scores were calculated using the sq.gr.spatial_neighbors() function in Squidpy (V.1.6.1), which generates Z-scores to quantify intercellular spatial proximity by comparing observed neighbor matrices to randomized controls. The detailed algorithms used for score computation have been previously described.35

Spatial ligand-receptor communication
Spatial ligand-receptor interaction scores were computed using the stLearn package (V.0.4.12). The curated interaction database connectomeDB2020_lit was loaded via the st.tl.cci.load_lrs() function. Interaction scores between all spatial spots were computed using st.tl.cci.run(). The spatial distribution scores of selected ligand-receptor pairs were visualized on a shared scale using st.pl.lr_result_plot().

Generation of KRASG12D/+; Villin-Cre Mice
Mice harboring the conditional KRASG12D allele (KRASLSL-G12D/+) were crossed with Villin-Cre transgenic mice, in which Cre recombinase is expressed specifically in intestinal epithelial cells under the control of the villin-1 promoter. Offspring were genotyped by PCR to identify those carrying both the KRASLSL-G12D allele and the Villin-Cre transgene. Cre-mediated recombination leads to activation of oncogenic KRASG12D specifically in the intestinal epithelium. All mice were maintained on a C57BL/6 background.

Animal models
All animal experiments were approved by the Animal Care and Use Committee of Nanjing Medical University and were performed according to the guidelines of the National Institutes of Health (IACUC-2406103). In this study, mice were randomly allocated into experimental groups, with each group consisting of n=6 animals. To minimize potential confounding variables, the order of treatment administration was randomized. Male KRASLSL-G12D and KRASG12D/+; Villin-Cre mice received a single intraperitoneal injection of azoxymethane (AOM, 12 mg/kg body weight; Sigma, USA) at 8 weeks of age. 1 week following AOM exposure, mice were given 2% dextran sulfate sodium (DSS, MP Biologicals, USA) in their drinking water for 7 days, followed by 2 weeks of regular water. This treatment cycle was repeated three times. At 20 weeks of age, mice were euthanized to assess and document the development of colorectal tumors.36 Tumors were dissociated into single-cell suspensions and processed for subsequent scRNA-seq.

Multiplex immunohistochemistry
Formalin-fixed paraffin-embedded tissue samples were sectioned at 3–5 µm thickness and mounted onto microscope slides. Tyramide Signal Amplification (TSA) fluorescence staining was performed with the Melady Biosciences TSA kit in accordance with the instructions provided by the manufacturer. Sections were baked at 60°C for 1 hour, dewaxed in xylene, and rehydrated through a graded ethanol series (100%, 95%, 80%, and 70%), followed by rinsing in distilled water.
Heat-induced antigen retrieval was conducted in EDTA buffer (pH 9.0) using a microwave oven. Endogenous peroxidase activity was quenched with 3% hydrogen peroxide, followed by blocking in a blocking solution for 20 min. Slides were then incubated overnight at 4°C with primary antibodies. After washing, sections were incubated with HRP-conjugated secondary antibodies for 20 min at room temperature, followed by TSA fluorophore staining. This process was repeated for each target antigen. After all targets were stained, sections were washed in phosphate-buffered saline and counterstained with 4′,6-diamidino-2-phenylindole (DAPI).37 Primary antibodies used in this study were as follows: EPCAM (Abcam, Cat# ab223582, 1:500); CEACAM1 (Abcam, Cat# ab300061, 1:2000); CTHRC1 (Abcam, Cat# ab256458, 1:500); S100A8 (Abcam, Cat# ab92331, 1:1500).

Statistical analysis
Statistical analyses were performed using Seurat (V.4.1.1) in R (V.4.2.0) and Scanpy (V.1.9.3) in Python (V.3.9.0). Unpaired, two-tailed Student’s t-tests were used to assess statistical significance; for multiple comparisons, Bonferroni correction was applied. Differences in cell-type proportions between groups were analyzed using the χ² test. Statistical significance was defined as p<0.05 (p<0.05; p<0.01; p<0.001; p<0.0001).

Data availability
The publicly available scRNA-seq data could be downloaded through the following links: https://www.ncbi.nlm.nih.gov/geo/ (accession numbers GSE132465,19 GSE144735,19 GSE161277,20 GSE200997,21 GSE221575,22 and GSE231559).23 The bulk transcriptomic data of the TCGA-COADREAD cohort and the dataset GSE20040738 used for chemotherapy sensitivity analysis in this study are available through the TCGA portal (https://www.cancer.gov/tcga) and the GEO database (accession number GSE200407), respectively. The scRNA-seq data generated from the mouse models in this study are available on request. The raw sequencing FASTQ files and processed data of the ST reported in this paper have been deposited in the Genome Sequence Archive39 in National Genomics Data Center,40 China National Center for Bioinformation/Beijing Institute of Genomics, Chinese Academy of Sciences (HRA011642) that are publicly accessible at https://ngdc.cncb.ac.cn/gsa-human.

Code availability
All essential code used for data analysis and result generation in this study is publicly available at github (https://github.com/yqcj455/KRAS_Mutant_Spatial_Niche).

Results

ScRNA-seq identifies KRAS mutation-associated epithelial cells
To begin, we integrated scRNA-seq data from 27 KRAS wild-type and KRAS-mutant CRC samples obtained from public scRNA-seq databases. We then used Scrublet to remove potential doublets, followed by filtering out low-quality cells with UMI <200 or with mitochondrial and ribosomal gene content exceeding 50%. This resulted in 136,890 high-quality cells. We proceeded with standardization and normalization, selecting the top 5,000 highly variable genes, and performed linear regression for mitochondrial and ribosomal gene content. Using Harmonypy, we corrected for batch effects in each sample. We visualized the data through UMAP, showcasing both database source (online supplemental figure S1A) and sample source (online supplemental figure S1B) distributions. These visualizations demonstrated uniform distribution across cell subpopulations in each sample, suggesting effective batch correction. We then performed PCA on the first 50 principal components and used Leiden clustering with a resolution of 1.2, which identified 29 subpopulations (online supplemental figure S1C).
Subsequently, we annotated each Leiden subpopulation using classical markers41 42 (online supplemental figure S1D). We identified subpopulations as follows: B cells (CD79A+, MS4A1+), endothelial cells (VWF+, CLDN5+), epithelial cells (EPCAM+, KRT19+), fibroblasts (LUM+, DCN+, ACTA2+), mast cells (TPSAB1+, TPSB2+), myeloid cells (CD14+, CD68+), pericytes and smooth muscle cells (SMCs) (PDGFRB+, RGS5+, MYH11+, TAGLN+), T and NK cells (CD3D+, CD3E+, NCAM1+, KLRD1+). The final UMAP plot for the two groups of cell types is shown in figure 1B, with a dot plot for marker expression (figure 1C).
Given the substantial intratumoral heterogeneity observed, not all epithelial cells exhibited features indicative of KRAS pathway activation. Thus, we subdivided the tumor epithelial cells into 10 subpopulations by Leiden clustering (figure 1D). We then analyzed 50 hallmark pathway features across these 10 subpopulations. Notably, the Epi_01 cell subpopulation had the highest KRAS signaling score, and the score for KRAS-mutant Epi_01 was higher compared with wild-type cells (figure 1E). Based on this, we identified Epi_01 as the KRAS mutation-associated epithelial cells. To identify the specific markers of this subpopulation, we used the FindAllMarker tool and identified the most specific differential gene for Epi_01 as CEACAM1 (carcinoembryonic antigen-related cell adhesion molecule 1) (online supplemental figure S1E).
Further analysis of the cellular composition revealed that the proportion of the Epi_01 subpopulation was similar between KRAS-mutant and KRAS wild-type groups (figure 1F). This observation led us to hypothesize that KRAS mutations may exert their biological effects predominantly through remodeling of the TME rather than through intrinsic expansion of this epithelial subclone.

Fib_CTHRC1 and Mono_S100A8 play key roles in KRAS-mutant CRC
We further investigated the differences in the TME between KRAS wild-type and mutant CRC by analyzing the changes in cell proportions between the two groups (online supplemental figure S2A). We found that in the KRAS-mutant group, the proportions of fibroblasts and myeloid cells were significantly higher, while the proportions of other cell types were similar. Therefore, we focused on fibroblasts and myeloid cells for further subclustering.
We first performed Leiden clustering on all fibroblasts, identifying six subpopulations (online supplemental figure S2B), which we named as Fib_CTHRC1, Fib_CXCL14, Fib_GAL, Fib_IGF1, Fib_MMP1, and Fib_PRKG1 based on the most specific differential genes identified by FindAllMarker. The heatmap of the most specific differential genes for each subpopulation was shown (online supplemental figure S2C). We then generated the UMAP plot for the fibroblasts in both groups (figure 2A) and analyzed the differences in fibroblast proportions (figure 2B). We found that the proportion of Fib_CTHRC1 showed the largest difference and was the most abundant among the fibroblast subpopulations. Therefore, Fib_CTHRC1 may play an important role in the KRAS-mutant microenvironment.
Next, we clustered myeloid cells into five subpopulations (online supplemental figure S2D) using Leiden clustering. Based on the most specific differential genes identified by FindAllMarker, we classified these myeloid subpopulations into macrophages and monocytes. Ultimately, we categorized all myeloid cells into DCs, Macro_C1QC, Macro_TNFRSF17, Mono_MALAT1, and Mono_S100A8. The heatmap of the most specific differential genes for each subpopulation was shown (online supplemental figure S2E). We generated the UMAP plot for myeloid cells in both groups (figure 2C) and analyzed the changes in myeloid cell proportions (figure 2D). We found that the proportion of Mono_S100A8 differed the most. Therefore, Mono_S100A8 may also play a significant role in the KRAS-mutant microenvironment.
Cell–cell communication within the TME is a crucial mechanism through which tumors exert their biological functions. To investigate this further, we first subclustered all major cell types, resulting in 34 subclusters (online supplemental figure S2F). We then used CellChat to analyze the receptor-ligand interactions between the three most important cell subpopulations in the KRAS-mutant microenvironment: Epi_01, Fib_CTHRC1, and Mono_S100A8. We compared the receptor-ligand interactions between the KRAS-mutant and wild-type groups. In the KRAS-mutant group, we observed that more Epi_01 emitted the ligand MDK, which acts on the receptor SDC4 on Mono_S100A8. Additionally, Fib_CTHRC1 secreted more ligands, such as COL1A1 and COL1A2, which act on the receptors ITGA2 and ITGB1 on Epi_01 (figure 2E).
We further evaluated the clinical significance of key cellular subsets—Epi_01, Fib_CTHRC1, and Mono_S100A8—in the KRAS-mutant CRC microenvironment. High expression of Epi_01 was significantly associated with poor overall survival (OS: p=0.014) and disease-free survival (DFS: p=0.035). Similarly, elevated Fib_CTHRC1 abundance correlated with unfavorable OS (p=0.019) and DFS (p=0.0057). Additionally, high Mono_S100A8 expression was also associated with poorer OS (p=0.022) in the TCGA-COADREAD dataset (online supplemental figure S3A and B).
In chemotherapy sensitivity analysis, all three cell types demonstrated significantly elevated abundance in chemotherapy-insensitive tumors (online supplemental figure S3C), with statistically significant differences observed for Epi_01 (p<0.01), Fib_CTHRC1 (p<0.01), and Mono_S100A8 (p<0.05).
To further evaluate the clinical relevance of the three key cellular subsets—Epi_01, Fib_CTHRC1, and Mono_S100A8—we analyzed their association with clinicopathological features in the TCGA-COADREAD cohort. We found that high Epi_01 abundance was significantly associated with advanced T stage (p=0.024) and the presence of lymph node metastasis (p=0.001) (online supplemental table S3). Similarly, elevated Fib_CTHRC1 levels were correlated with higher T stage (p=0.007) and lymph node metastasis (p<0.001) (online supplemental table S4). In contrast, Mono_S100A8 abundance showed no significant associations with standard tumor node metastasis staging parameters (online supplemental table S5). This disconnect suggested that the protumorigenic role of Mono_S100A8 is not primarily in driving local tumor invasion, but rather in establishing an immunosuppressive microenvironment that facilitates metastatic spread and treatment resistance, as reflected in its strong association with poor survival outcomes and chemotherapy insensitivity.
In summary, KRAS-mutant tumors exhibit distinct cell interaction patterns compared with wild-type tumors, with Fib_CTHRC1 and Mono_S100A8 playing crucial roles in KRAS-mutant CRC.

Spatial analysis reveals greater heterogeneity in KRAS-mutant CRC
Tumor heterogeneity is widely recognized as one of the key factors driving therapeutic resistance. The greater the heterogeneity within a tumor, the poorer the efficacy of drug treatments.43 44 To analyze tumor heterogeneity, we employed the inferCNV method to assess the CNV status of each spot. We selected spots expressing high levels of immune cell markers such as PTPRC, CD3D, and CD3E as reference spots for calculating the CNV score of each point.
Using Leiden clustering based on the CNV scores for each spot, we observed that both CRC_MUT_A1 and CRC_MUT_B1 in the KRAS-mutant group exhibited five distinct CNV status subgroups. In contrast, the CRC_WT_C1 and CRC_WT_D1 subgroups in the wild-type group only had four CNV status subgroups (figure 3A). This suggests that KRAS-mutant tumors exhibit a greater number of mutation states, reflecting higher heterogeneity. Such heterogeneity may be associated with the therapeutic resistance observed in KRAS-mutant tumors.
Additionally, we displayed the CNV scores for each spot in the ST analysis (figure 3B). Regions with higher CNV scores were hypothesized to correspond to tumor areas, while those with lower CNV scores were more likely to represent stromal regions.
Next, we categorized each tissue section into different niches based on the gene expression profiles of individual spots (figure 3D). We then presented heatmaps showing the average CNV levels of each niche for both KRAS-mutant (online supplemental figure S4A) and KRAS wild-type (online supplemental figure S4B) CRC.
In conclusion, the CNV scoring reveals that KRAS-mutant tumors exhibit greater heterogeneity, which may imply a higher degree of malignancy and greater therapeutic resistance.

Cell-type deconvolution and multimodal intersection analysis
To achieve higher-resolution analysis at the cellular level, we performed deconvolution of the ST data using the scRNA-seq dataset containing 34 cell subtypes. Using these datasets, we applied the Cell2location deconvolution algorithm to map the scRNA-seq-defined cell states onto our ST data. This process provided us with the cell abundance scores for each spot across the 34 cell subtypes.
We specifically focused on the abundance of Epi_01, Fib_CTHRC1, and Mono_S100A8 (figure 3C). We found that both CRC_MUT_A1 and CRC_MUT_B1 in the KRAS-mutant group had higher abundances of these three cell types compared with the wild-type groups, CRC_WT_C1 and CRC_WT_D1. This further confirmed the differences in the microenvironment between KRAS-mutant and wild-type tumors.
Next, we extracted differentially expressed genes from both scRNA-seq and ST and performed multimodal intersection analysis to calculate the overlap between cell-type-specific and region-specific gene sets. We used a p value threshold of <0.01 and performed hypergeometric testing to assess significant enrichment. Finally, we presented the enrichment levels of 34 cell types across different niches in the form of a heatmap (figure 3E).
In summary, we accurately described the abundance of the 34 cell types at the ST level, providing a foundation for our subsequent analysis of cell type co-localization and cellular network interactions.

Spatial architecture of KRAS-mutant CRC
To investigate how cell types interact and distribute spatially in KRAS-mutant CRC, we employed MISTy, a machine learning framework designed to interpret spatial relationships. This approach aimed to assess whether the composition of key cellular populations within transcriptomic spots could be inferred based on their surrounding microenvironment. We specifically evaluated whether the cellular composition in the local vicinity—defined by relative cell type abundances—could predict the presence of specific tumor-associated compartments across juxta-view spatial scales (~200 µm, five-spot radius).
Notably, within tumor regions, spatially co-occurring cell types within single tissue spots were found to be strong predictors, highlighting recurrent spatial associations between KRAS-mutant epithelial cells and non-malignant components, including fibroblasts and myeloid cells. We generated co-localization heatmaps for two KRAS-mutant tumors and two wild-type tumors (figure 4A). Notably, in the two KRAS-mutant tumor tissues (CRC_MUT_A1 and CRC_MUT_B1), there was a strong co-localization between Epi_01 and Mono_S100A8, whereas such co-localization was weaker in the two wild-type tumors (CRC_WT_C1 and CRC_WT_D1). However, we did not observe co-localization between Fib_CTHRC1 and either Epi_01 or Mono_S100A8. We hypothesize that Fib_CTHRC1 might form a fibrous barrier surrounding tumor cells. To test this hypothesis, we first constructed homotypic cell networks by calculating the cell degree of each cell type, which was scored on a scale of 0–6 (figure 4B). The higher the score, the darker the spatial color representation. We found that in CRC_MUT_A1 and CRC_MUT_B1, the homotypic cell network of Epi_01 and Mono_S100A8 was close in spatial proximity, further confirming their co-localization. Additionally, we found that the homotypic cell network of Fib_CTHRC1 partially overlapped with Epi_01 and Mono_S100A8, suggesting that Fib_CTHRC1 may act to form a fibrous barrier around Epi_01. In contrast, although CRC_WT_C1 and CRC_WT_D1 displayed homotypic cell networks for Epi_01, they lacked networks for Mono_S100A8 and Fib_CTHRC1, suggesting that wild-type CRC tumors do not feature the same cellular network involving Epi_01, Fib_CTHRC1, and Mono_S100A8.
Previous studies have indicated that the disruption of cellular neighborhoods in cancer reflects disturbances in the natural ecological system, which impacts disease progression and treatment response. To identify these multicellular communities and their spatial distribution in CRC, we used ISCHIA, a probabilistic methodology for identifying cellular neighborhoods through co-occurrence analysis, inspired by models of species co-occurrence in ecology. This composition-based clustering approach allowed us to reconstruct the local TME architecture based on spatial co-occurrence and potential interactions between diverse cell types within each spot. Clustering the spots from all samples based on their cell-type compositions revealed six heterogeneous compositional clusters, which we termed “spatial ecotypes” (figure 5A).
To further characterize the spatial architecture and cellular composition of each spatial ecotype, we analyzed the distribution patterns and cell type profiles across all six composition clusters (CC1–CC6). As shown in figure 5A, CC5 exhibited a distinct spatial localization pattern, predominantly concentrated in the tumor core regions of KRAS-mutant tumors (CRC_MUT_A1 and CRC_MUT_B1), where Epi_01 and Mono_S100A8 were co-localized at high abundance. In contrast, CC6 displayed a more peripheral and patchy distribution, often bordering or surrounding CC5-rich areas, suggesting its potential role as a fibrous barrier enclosing the tumor mass. This spatial arrangement was absent in wild-type tumors (CRC_WT_C1 and CRC_WT_D1), which lacked both CC5 and CC6 structures.
We next examined the cellular composition of each ecotype using box plots and stacked bar charts (figure 5B,C). CC1 was primarily composed of SMCs and endothelial cells, forming a vascular-associated stromal region. CC2 consisted of a mixture of various epithelial subtypes (including Epi_02, Epi_05, Epi_08, Epi_09, Epi_10) at low abundances, indicating a heterogeneous epithelial community. CC3 was dominated by Epi_02, Epi_03, Epi_04, Epi_09, and proliferative T cells, potentially reflecting an early tumor expansion zone. CC4 shared a similar cellular composition with CC2 but exhibited significantly higher abundances of these cell types, suggesting it may represent a more active or advanced state of CC2. Notably, compositional cluster CC5 possessed the highest proportions of Epi_01 and Mono_S100A8, defining it as the KRAS-mutant spatial ecotype, primarily characterized by the co-localization and cellular interactions between Epi_01 and Mono_S100A8. We also found that CC6 was mainly composed of Fib_CTHRC1-rich stromal cells and contained numerous CD4+ and CD8+ lymphocytes. This suggests that CC6 may represent a fibrous barrier spatial ecotype that excludes lymphocyte infiltration.
Mono_S100A8 has been reported to inhibit CD8+T cell function by recruiting myeloid-derived suppressor cells, thereby promoting immune suppression in tumors.45 46
CTHRC1 (collagen triple helix repeat containing 1), a marker of myofibroblasts, can secrete collagen to form a fibrous barrier that resists CD8+T cell infiltration.47 Moreover, CTHRC1 can upregulate CCL15 to recruit tumor-associated macrophages, facilitating CRC progression.48
In conclusion, we identified that the spatial ecotype of KRAS-mutant CRC is primarily characterized by the co-localization of Epi_01 and Mono_S100A8, with stromal cells like Fib_CTHRC1 surrounding the periphery to form a fibrous barrier that excludes lymphocytes. Both Mono_S100A8 and Fib_CTHRC1 contribute to immune suppression, together forming an immunosuppressive niche in KRAS-mutant CRC.

Spatial heterotypic cell network in KRAS-mutant CRC
To study the spatial interaction networks between different cell types, we first normalized the 34 cell types at each spot using 0–1 standardization. We then assigned each spot to the cell type with the highest score. Next, we calculated the adjacency matrix for each cell type using Squidpy (figure 6A). To visualize the cell-type adjacency relationships more clearly, we used a chord diagram (figure 6B).
We found that in CRC_MUT_A1 and CRC_MUT_B1, Epi_01, Fib_CTHRC1, and Mono_S100A8 exhibited more frequent adjacency relationships. In contrast, these interactions were significantly absent in CRC_WT_C1 and CRC_WT_D1. This further supports the idea that the adjacency between Fib_CTHRC1 and Epi_01/Mono_S100A8 may allow Fib_CTHRC1 to protect tumor cells from lymphocyte-mediated killing by surrounding the tumor cells at the periphery.
The spatial network of cell types clearly illustrated that in CRC_MUT_A1 and CRC_MUT_B1, Epi_01, Fib_CTHRC1, Fib_MMP1, and Mono_S100A8 were in close proximity and formed a network (online supplemental figure S5A). On the other hand, CRC_WT_C1 was mainly composed of Epi_08 and Epi_10, while CRC_WT_D1 was primarily composed of Epi_02, Epi_03, and Epi_04 (online supplemental figure S5B).
In summary, the spatial heterotypic cell network analysis further validates the immunosuppressive niche in KRAS-mutant CRC.

Spatial ligand-receptor interactions in KRAS-mutant CRC
The preceding analysis had provided preliminary insights into the immunosuppressive niche in KRAS-mutant CRC. However, the specific interactions between Epi_01, Fib_CTHRC1, and Mono_S100A8 remained poorly understood. Previous scRNA-seq analysis using CellChat revealed that in the KRAS-mutant group, a higher proportion of Epi_01 secreted the ligand MDK, which acts on the receptor SDC4 on Mono_S100A8. Furthermore, Fib_CTHRC1 was found to secrete increased levels of the ligands COL1A1 and COL1A2, which interact with the receptors ITGA2 and ITGB1 on Epi_01. To validate the effectiveness of these ligand-receptor interactions in ST, we used the stlearn tool to analyze their interaction scores. Our findings indicated that the ligand-receptor interaction scores for COL1A1_ITGA2 and COL1A1_ITGB1 were significantly higher in the KRAS-mutant group compared with the KRAS wild-type group (online supplemental figure S6A). Similarly, the interaction scores for COL1A2_ITGA2 and COL1A2_ITGB1 were also significantly higher in the KRAS-mutant group (online supplemental figure S6B). The COL1A1 gene encodes the α1 chain, while COL1A2 encodes the α2 chain. The native collagen protein is composed of two α1 chains and one α2 chain, forming a heterotrimeric structure.49 ITGA2 and ITGB1 can form a heterodimer (α2β1), a complex that plays a pivotal role in various biological processes. For instance, ITGA2 and ITGB1, through binding with collagen, participate in extracellular matrix formation and tissue regeneration.50 Moreover, they play crucial roles in the TME, such as promoting tumor cell migration and invasion.51 Therefore, Fib_CTHRC1 in KRAS-mutant tumors may secrete more collagen, interacting with integrin receptors on KRAS-mutant tumor cells, thereby facilitating their migration and invasion.
Additionally, we observed that the ligand-receptor interaction score for MDK_SDC4 was also significantly higher in the KRAS-mutant group compared with the KRAS wild-type group (online supplemental figure S6C). Previous studies have reported that MDK can bind to SDC4, promoting the chemotaxis of regulatory T cells, thus suppressing antitumor immune responses. Consequently, the increased secretion of MDK by Epi_01 in KRAS-mutant tumors, which interacts with the SDC4 receptor on Mono_S100A8, may contribute to the chemotaxis of Mono_S100A8 toward the tumor cells, thereby exerting an immunosuppressive function.
Overall, KRAS-mutant tumor cells may use MDK_SDC4 interactions to attract Mono_S100A8, while surrounding Fib_CTHRC1 may modulate KRAS-mutant tumor cell behavior through collagen-integrin interactions, collectively establishing an immunosuppressive microenvironment.

Fib_CTHRC1 and Mono_S100A8 play key roles in KRASG12D/+; Villin-Cre mice
To simulate the TME of KRAS mutations in vivo, we generated KRASG12D/+; Villin-Cre mice (figure 7A). Subsequently, we established AOM/DSS models for both KRASLSL-G12D and KRASG12D/+; Villin-Cre mice (figure 7B). On gross examination, we observed more extensive tumor formation in the colons of KRASG12D/+; Villin-Cre mice compared with KRASLSL-G12D mice (figure 7C). Quantitative analysis revealed a significant difference in tumor burden between the two groups. The colon length of KRASG12D/+; Villin-Cre mice was notably shorter, indicating more severe disease progression. Furthermore, these mice had a higher number of tumors per colon, particularly those larger than 2 cm in size (figure 7D). These findings suggested that activation of KRASG12D in intestinal epithelial cells promoted the development and progression of colitis-associated cancer.
Next, we performed scRNA-seq on tumor tissues from the mouse intestines. Using an analytical approach similar to that used for human tumors, the mouse tumor cells were categorized into B_cells, endothelial_cells, epithelial_cells, fibroblasts, myeloid_cells, and T and NK_cells (online supplemental figure S7A). Cell-type-specific markers were displayed using a dot plot (online supplemental figure S7B). We further analyzed the TME in mice with KRAS mutations, presenting UMAP plots of the cell types in the two groups (figure 7E) and analyzing their proportional changes (figure 7F). Notably, the proportion of fibroblasts and myeloid cells was significantly higher in tumors from the KRAS-mutant group, closely resembling the TME observed in humans.
To explore whether Fib_CTHRC1 and Mono_S100A8 also played key roles in KRASG12D/+; Villin-Cre mice, we extracted and subclustered the fibroblasts and myeloid cells. All fibroblasts were clustered via Leiden algorithm into three subgroups: Fib_Col12a1, Fib_CTHRC1, and Fib_Mettl7b. These subgroups were named based on their most specific differentially expressed genes identified by FindAllMarkers, and their marker genes were shown in a heatmap (online supplemental figure S7C). We displayed the UMAP of fibroblasts from both groups (figure 7G) and analyzed the proportional changes (figure 7H), revealing that Fib_CTHRC1 was more abundant in KRASG12D/+; Villin-Cre mice.
Similarly, the myeloid cells were clustered into two subgroups: Macro_C1qa and Mono_S100A8, with their specific marker genes visualized in a heatmap (online supplemental figure S7D). UMAP plots for Myeloid cells in both groups were shown (figure 7I), and proportional analysis demonstrated an increase in both the proportion and number of Mono_S100A8 cells in the KRASG12D/+; Villin-Cre mice (figure 7J).
To further investigate the functional crosstalk between these key stromal and immune subsets and mutant epithelial cells, we performed cell–cell communication analysis using CellChat on the mouse scRNA-seq data. Differential interaction analysis revealed significantly enhanced ligand-receptor signaling between Fib_CTHRC1, Mono_S100A8, and KRAS-mutant epithelial cells in KRASG12D/+; Villin-Cre mice compared with controls (online supplemental figure S8A). Consistently, the interaction heatmap also indicated markedly strengthened ligand-receptor signals between Fib_CTHRC1, Mono_S100A8, and KRAS-mutant epithelial cells (online supplemental figure S8B).
In human scRNA-seq data, we previously observed that epithelial cells signal to Mono_S100A8 via the MDK ligand binding to its receptor SDC4, and that Fib_CTHRC1 interacts with epithelial cells through COL1A1/COL1A2 binding to the integrin complex ITGA2/ITGB1. Guided by these findings, we examined the corresponding ligand-receptor pairs in the mouse scRNA-seq dataset and obtained consistent results (online supplemental figure S8C).
Our results indicate that Fib_CTHRC1 and Mono_S100A8 are active contributors, rather than mere biomarkers, to a protumorigenic niche in the KRAS-mutant TME. Furthermore, this functionally validated intercellular crosstalk provides a mechanistic basis for the aggressive tumor phenotype in KRAS-mutant mice, thereby defining a coordinated ecosystem supportive of immune evasion and stromal activation that mirrors the spatial ecotypes in human KRAS-mutant CRC.

Cellular co-localization in KRAS-mutant CRC
To validate the co-localization of the three key cell populations—Epi_01, Fib_CTHRC1, and Mono_S100A8—that define the immunosuppressive microenvironment in KRAS-mutant tumors, we performed multiplex immunohistochemistry (mIHC) on tissue samples from three KRAS-wild-type and three KRAS-mutant CRC tissues. We selected CEACAM1 as the most specific marker for Epi_01, along with the classical epithelial marker EPCAM, and included CTHRC1 and S100A8 (S100 calcium binding protein A8) as targets for staining. In KRAS-wild-type tumors, the fluorescence intensities of CEACAM1, CTHRC1, and S100A8 were markedly lower compared with those observed in KRAS-mutant tumors (figure 8A). Notably, in KRAS-mutant tumors, CEACAM1, CTHRC1, and S100A8 exhibited substantial spatial co-localization. CEACAM1 and S100A8 were largely co-distributed, while CTHRC1 was partly located together with them and also partly found surrounding them at the edges (figure 8B).
Furthermore, ST data visualized via the SpatialFeaturePlot() function confirmed the expression patterns of these four genes (online supplemental figure S9A-D), showing a high degree of consistency with the mIHC results and further supporting the distinct cellular architecture of the KRAS-mutant TME.
Taken together, our findings suggest a novel mechanism by which KRAS-mutant epithelial cells recruit Mono_S100A8 via MDK_SDC4 signaling. In parallel, Fib_CTHRC1 secretes increased levels of collagen that not only promote tumor cell invasion through integrin engagement, but also form a physical barrier that limits lymphocyte infiltration. Through these coordinated actions, these three cell populations establish an immunosuppressive niche characteristic of KRAS-mutant CRC. A schematic illustration of this mechanism is presented in figure 8C.

Discussion
KRAS is one of the most frequently mutated oncogenes in CRC and plays a pivotal role in tumor initiation, progression, immune evasion, and therapeutic resistance.52 53 Given the substantial heterogeneity associated with KRAS mutations, there is an urgent need for more precise and individualized treatment strategies for this patient population.54 Single-cell and ST technologies are powerful tools for dissecting tumor heterogeneity.5557 This study represents the first integrated application of scRNA-seq and ST to systematically characterize the spatial heterogeneity of KRAS-mutant versus wild-type colorectal tumors, uncovering a distinct KRAS mutation-associated immunosuppressive niche. Specifically, this KRAS-mutant spatial niche comprises co-localized KRAS-mutant epithelial cells and Mono_S100A8, surrounded by Fib_CTHRC1. These findings suggest that KRAS mutations may contribute to immune escape and therapy resistance by reprogramming the local immune ecosystem.
Notably, CEACAM1, a transmembrane glycoprotein of the CEA family, is identified as a specific marker of KRAS-mutant epithelial cells. CEACAM1 is widely expressed in tissues and is involved in intercellular adhesion, signal transduction, immune regulation, and angiogenesis.58 59 Previous studies have shown that CEACAM1 regulates TIM3-mediated tolerance and exhaustion, and that its upregulation in liver cancer stem cells confers resistance to natural killer cell-mediated cytotoxicity.60
S100A8, also known as calgranulin A or MRP8, is a calcium-binding protein of the S100 family.61 Li et al demonstrated that S100A8 promotes epithelial-mesenchymal transition and metastasis in CRC via the TGF-β/USF2 axis.62 In esophageal cancer, S100A8 has been implicated in chemoresistance through the reprogramming of cancer-associated fibroblasts.63
CTHRC1, a protein closely associated with fibrosis, has recently been implicated in tumor progression and metastasis.64 It activates fibroblasts, enhances extracellular matrix production, and contributes to the formation of an immunosuppressive microenvironment.65
Together, our findings suggest that these KRAS-mutant spatial niche populations—KRAS-mutant epithelial cells, Mono_S100A8, and Fib_CTHRC1—collectively contribute to immunosuppression. Targeting these factors may provide a therapeutic strategy for patients with KRAS-mutant CRC.
Several previous studies have elucidated KRAS-driven mechanisms in CRC. For instance, Hsu et al reported that oncogenic KRAS promotes adipofibrosis and angiogenesis to drive cancer progression.23 Zhou et al revealed that mutant KRAS induces immunosuppression by impairing DDX60-mediated dsRNA accumulation and viral mimicry.66 Although Chong et al have characterized the features of KRAS mutations through multiomics approaches, no studies to date have explored the ST landscape of KRAS-mutant tumors.67 In this study, we provide a novel spatial perspective by delineating the KRAS mutation-associated spatial niche using ST and identify potential ligand-receptor interactions that may underlie cellular communication within this niche. These findings offer a theoretical foundation and potential molecular targets for the development of innovative immunomodulatory strategies tailored to KRAS-mutant CRC.
The clinical implementation of any novel biomarker strategy must overcome significant practical hurdles. While ST serves as a powerful discovery tool, its direct clinical adoption is currently constrained by high costs, analytical complexity, and lengthy turnaround times. Therefore, our translational strategy does not propose the routine use of ST in clinics. Instead, we focus on translating the core molecular signature we identified—CEACAM1, S100A8, and CTHRC1—into a clinically feasible format. The most immediate and practical path is the development and validation of an mIHC assay targeting these three proteins.68 Such an assay would be cost-effective, rapidly deployable in most pathology departments, and could identify patients with this high-risk immunosuppressive niche for tailored therapeutic strategies.
However, we acknowledge that the current lack of real-world evidence for this biomarker signature is a major limitation for its immediate clinical availability. Its future reimbursement by healthcare systems will be contingent on demonstrating a clear impact on treatment decisions and patient outcomes in large, prospective clinical trials. Further steps for clinical integration include rigorous analytical and clinical validation of the mIHC assay, standardization of scoring protocols, and ultimately, the demonstration of cost-effectiveness compared with conventional treatment.
The growing use of liquid biopsy for monitoring tumor dynamics provides an opportunity to integrate spatial information with temporal molecular data.69 While liquid biopsy excels at detecting molecular alterations over time, it cannot capture the spatial organization of the TME. Our spatial characterization of the KRAS-mutant niche provides important contextual information that complements liquid biopsy findings. Future diagnostic frameworks could combine baseline spatial assessment of the TME via tissue-based assays with longitudinal monitoring through liquid biopsy, offering a more comprehensive understanding of tumor biology.70
Nevertheless, our study has several limitations. First, the scRNA-seq datasets used for spatial deconvolution were derived from public databases rather than matched samples from the same patients as the ST, which may introduce discrepancies due to inter-patient heterogeneity. We attempted to reduce the impact of this limitation by integrating a wide range of publicly available single-cell datasets. Second, the ST analysis was conducted on a limited number of samples—two KRAS-mutant and two KRAS-wild-type CRC tissues—due to constraints in tissue availability and the high cost of ST library preparation and sequencing. Third, the spatial resolution of the ST platform used restricts our analysis to multicellular “neighborhoods” and precludes the precise delineation of direct cell-to-cell contacts. Future studies employing higher-resolution spatial technologies, such as Visium HD or multiplexed imaging platforms like Xenium, will be crucial to resolve these interactions at a single-cell level. Fourth, our study cohort was not stratified by specific KRAS mutation subtypes (eg, G12D, G12V, G13D). Given emerging evidence that different KRAS alleles can drive distinct biological programs and TMEs, subsequent investigations with larger, genotypically detailed cohorts are warranted to dissect these potential subtype-specific spatial architectures. Moreover, while we propose that collagen secreted by Fib_CTHRC1 excludes CD8+T cell infiltration, this hypothesis lacks direct experimental validation. Finally, although mIHC offered supportive spatial evidence for the co-localization of key cell populations and marker expression, we did not perform in vitro or in vivo functional assays to directly validate the immunosuppressive roles of CEACAM1, S100A8, and CTHRC1. Future studies will aim to further delineate the molecular mechanisms underpinning cellular interactions in the KRAS-mutant TME.
To constructively address these limitations, we propose several specific directions for future research. First, prospective studies using patient-matched scRNA-seq and high-resolution ST data would help minimize technical artifacts and biological variance. Additionally, systematic exploration of the spatial-immune landscape across major KRAS allele-specific mutations is needed to delineate potential subtype-specific therapeutic vulnerabilities. Finally, direct functional validation of the proposed immunosuppressive axis remains essential, particularly through assessing the efficacy of combined immunotherapeutic strategies targeting both the MDK-SDC4 and CTHRC1-integrin pathways in preclinical models of KRAS-mutant CRC.
In conclusion, by integrating scRNA-seq and ST, we have identified a KRAS mutation-associated immunosuppressive niche in CRC. KRAS-mutant epithelial cells recruit Mono_S100A8 via MDK_SDC4 signaling, while surrounding Fib_CTHRC1 fibroblasts secrete collagen that engages integrin receptors on KRAS-mutant epithelial cells and forms a physical barrier to lymphocyte infiltration. These interactions collectively establish an immunosuppressive spatial niche in KRAS-mutant CRC.

Supplementary material
10.1136/jitc-2025-013763online supplemental file 110.1136/jitc-2025-013763online supplemental file 210.1136/jitc-2025-013763online supplemental file 310.1136/jitc-2025-013763online supplemental file 410.1136/jitc-2025-013763online supplemental file 510.1136/jitc-2025-013763online supplemental file 610.1136/jitc-2025-013763online supplemental file 7