Molecular cartography uncovers evolutionary and microenvironmental dynamics in sporadic colorectal tumors

Abstract

Colorectal cancer exhibits dynamic cellular and genetic heterogeneity during progression from precursor lesions toward malignancy. Analysis of spatial multi-omic data from 31 human colorectal specimens enabled phylogeographic mapping of tumor evolution that revealed individualized progression trajectories and accompanying microenvironmental and clonal alterations. Phylogeographic mapping ordered genetic events, classified tumors by their evolutionary dynamics, and placed clonal regions along global pseudotemporal progression trajectories encompassing the chromosomal instability (CIN+) and hypermutated (HM) pathways. Integrated single-cell and spatial transcriptomic data revealed recurring epithelial programs and infiltrating immune states along progression pseudotime. We discovered an immune exclusion signature (IEX), consisting of extracellular matrix regulators DDR1, TGFBI, PAK4, and DPEP1, that charts with CIN+ tumor progression, is associated with reduced cytotoxic cell infiltration, and shows prognostic value in independent cohorts. This spatial multi-omic atlas provides insights into colorectal tumor-microenvironment co-evolution, serving as a resource for stratification and targeted treatments.

Keywords: colorectal cancer; immune exclusion; microsatellite instability; multiplex imaging; mutations; spatial transcriptomics; stem cells; tumor evolution; tumor progression.

Figure 1.. Spatial atlas queries layers of molecular heterogeneity in sporadic colorectal tumors.

(A) Diagram detailing colorectal specimens chosen for atlas experiments. (B) Patient-level information from A in table format. Number within ”# Visium” squares indicates the number of 6.5 × 6.5 mm ST tiles covering each specimen. Previously collected or newly generated scRNA-seq data were available for most specimens. All specimens have MxIF imaging and bulk or multiregional WES. (C) Experimental design consisting of layered spatial molecular assays from serial sections of FFPE tissue blocks. (D) Diagram of phylogeographical cartography from multiregional sequencing (LCM-WES) data, layered with ST and MxIF images. Black arrows represent progression pseudotime (PPT) inferred from the phylogenetic relationships between ROIs. (E) Summary of gene signature scores by tumor type (left), ST patient (middle), and matched scRNA-seq patient (right). Patient ID colors represent tumor type. Mean signature expression scaled across groups. In this dotplot and hereafter, size of dots represents expression frequency, while shade represents intensity. (F) Somatic mutations detected in LCM-WES samples, summarized by patient and grouped by biological pathway. Top barplot represents overall TMB breakdown by mutation class per patient. See also Figure S1 and Tables S1-S2.

Figure 2.. CNV inference establishes spatially resolved tumor clones and their phylogenetic relationships.

(A) UMAP embeddings generated from inferred CNV profiles of all ST samples colored by tumor type, CNV score, Patient, and PAT71397 CNV clone to accompany panels E-H. Individual points represent ST microwells, which were subsetted to major clone regions prior to embedding. (B) Boxplots of CNV scores for all ST microwells in major CNV clone regions across atlas, grouped by sample type. CRC n = 48,439; NL n = 1,067; SSL/HP n = 1,735; TA/TVA n = 1,951. Samples are distributed by type according to Figure 1B. (C) Boxplots of CNV scores for epithelial cells from Chen, et al. cohort, grouped by sample type. CRC n = 63,371; NL n = 31,917; SSL/HP n = 11,896; TA/TVA n = 21,275. Samples are distributed by CRC: n=65, NL: n=37, SSL/HP: n=11, TA/TVA: n=14. **** p <1.00 × 10⁻⁴ by ANOVA followed by Tukey post-test, showing significance between CRC and TA/TVA group. (D) Summary of MxIF intensities, cell activity, and immune gene signatures by major tissue domains determined through CNV inference. (E) CNV scores (left) and tumor clone regions (right) for PAT71397. (F) MxIF with inferred progression trajectory for PAT71397. Scale bars 500 μm. (G) Summary of TMB, CNV score, and gene signatures for CNV clone regions of PAT71397. (H) Heatmap of inferred CNVs for PAT71397 ST, corresponding to E-G (top), as well as CNVs measured by WGS and WES for PAT71397 blocks and additional selected premalignant tumors (bottom). Brackets connect WES and WGS from PAT71397 malignant (MSS) and benign (TVA) blocks to dominant CNV clones in respective ST to show similarity of measured and inferred CNV profiles. See also Figure S2 and Table S2.

Figure 3.. Multiregional somatic mutational profiles provide phylogeographical topology.

(A) Oncoplot of detected driver mutations within spatially sampled LCM ROIs of PAT71397. (B) Phylogenetic tree for PAT71397. Length of branches are proportional to the number of shared or private somatic mutations in each LCM ROI. (C) Diagram of LCM ROIs in PAT71397 blocks, overlaid on CNV clone regions identified in ST. Black arrow represents inferred progression trajectory from CNVs and mutational phylogeny in B. (D) Diagram of observed modes of tumor evolution. Example phylogenetic trees from representative atlas samples shown to the right of each diagram. Patient ID colors represent tumor type (MMR status). (E) Tumor regions and their clinical and mutational metadata divided by class and ordered left-to-right by corresponding PPT (CNV score for CIN+, TMB for HM). (F) Summary of gene signatures across all tumor regions grouped by evolutionary mode from D. (G) CIN index versus PPT for tumor regions from E. Points are colored by tumor class, except pre-malignant and normal regions, which are colored according to tumor type as in Figure 1A-B. Point shape corresponds to regions with detected TP53 mutation. See also Figure S3 and Table S2.

Figure 4.. Cell-state deconvolution reveals pseudotemporal tissue dynamics.

(A) Example whole slide MxIF surrounded by refNMF usages for seven cell states as well as MILWRM tissue domain projected onto PAT30884 histology. Scale bar 500 μm. (B) refNMF usages of normal absorptive colonocyte (ABS), normal fibroblast (FIB2), serrated-specific cell (SSC), and goblet cell (GOB) states for HTA11_08622_A. (C) MxIF image of HTA11_08622_A. Scale bar 500 μm. (D) MILWRM tissue domains for HTA11_08622_A, surrounded by top cell-state loadings for SSL (D0), normal epithelium (D5), and submucosa (D6) domains. (E) Proportions of MILWRM domains detected in ST from each patient. Patient ID colors represent tumor class. (F) refNMF states grouped by compartment and summarized across tumor stage, tumor class, MILWRM domain, and patient for all ST samples. Patient ID colors represent tumor class. (G) Heatmap of GAM fits for refNMF states in all ST tumor regions ordered by PPT for HM (left) and CIN+ (right) tumors. Color represents scaled expression within each tumor class. See also Figure S4 and Tables S2-S4.

Figure 5.. Gene expression features of CIN+ CRCs predict immune exclusion.

(A) Heatmap of GAM fits for top genes summarized across all ST tumor regions. Color represents scaled expression within each tumor class. Genes are grouped by biological function. Bracket denotes IEX genes. (B) Pairwise Pearson correlations between progression indicators (CNV score and iCMS2), IEX, cytotoxic T cell refNMF state (TL2), and CD8 T cell gene signature in all CIN+ tumor regions. (C) Genes, gene signatures, and refNMF states grouped into pseudotime indicators (”PPT”), immune exclusion markers (”Excl.”), microenvironmental cells (”uEnv.”), infiltrating immune cells (”Inf.”), tumor activity (”Act.”), and epithelial-specific markers of MSS, MSI-H, and normal mucosa summarized by MILWRM domain, tumor class, and patient for all ST samples. Patient ID colors represent tumor class. (D) PAT71662 ST with annotated tissue domains from MILWRM (left) and CNV clone regions (right). (E) Expression overlay and spatial co-occurrence analysis for IEX, helper T cells (TL1), and cytotoxic T cells (TL2) in PAT71662. Line plots at right indicate the conditional probability of high signature or cell-state expression as a function of distance from CNV clone 1 microwells (STAR Methods: Spatial co-occurrence analysis from ST). (F) PAT71662 MxIF showing collagen, CDX2 (marking MSS epithelium), and lymphocytes (CD3 and CD8). Inset highlights CD3/CD8+ cells sequestered to stroma. Scale bars 500 μm. (G) Centroids of segmented single cells from PAT71662 MxIF plotted in whole-slide space, split into lymphoid and myeloid compartments. (H) Same as F for PAT73458. PCNA and MUC5AC mark tumor epithelium. Inset highlights CD8+ cells invading epithelium. Scale bars 500 μm. (I) Same as in G, for PAT73458. (J-K) Same as in D-E, for PAT73458. TL3 represents γδIELs. (L) Census of infiltrating immune cells in PAT71662 (bottom) and PAT73458 (top) from G and I summarized by CNV clone region (STAR Methods: MxIF immune exclusion analysis). (M) Number of infiltrating CD8+ T cells detected in MxIF plotted against IEX score for all tumor regions. Points are colored by tumor class and sized according to PPT ranking. IEX score threshold of 4.99 represents mean + 1 SD of 36,583 microwells in smooth muscle areas across ST atlas, indicating background IEX score level. n=20 CIN+ and n=6 HM tumor regions. (N) CD8+ T cells as a percentage of total non-epithelial cells detected in 66 scRNA-seq specimens plotted against the average IEX score from each sample’s epithelium. Points are colored by tumor pathway (NL:normal, n=27; CONV:conventional, n=20; SER:serrated, n=19). IEX score threshold of 1.26 represents mean + 1 SD of 22,262 non-epithelial cells from the Chen, et al. and Pelka, et al. cohorts.* p <5.00 × 10⁻², *** p <1.00 × 10⁻³ by Pearson correlation to test for mutual exclusivity. See also Figure S5 and Tables S2-S4.

Figure 6.. IEX trends with tumor progression and predicts poor patient outcomes.

(A-C) UMAP of paired scRNA-seq data generated for CRC sample with ST data in the study, overlaid with (A) individual samples, (B) epithelial versus non-epithelial annotation, and (C) IEX score calculated from normalized gene counts. (D) Boxplot of percentage of CD8+ T cells per sample for immune excluded- (n=8) versus immune excluded+ (n=10) tumors based on IEX score derived from paired scRNA-seq. *** p <1.00 × 10⁻³ by 2-sided Student’s t-test. (E) HCR-FISH for TGFBI, DPEP1, PAK4, DDR1 transcript in an IEX+ MSS CRC. Scale bars 100 μm. (F) Heatmap of GAM fits for genes, gene signatures, and refNMF cell states summarized across all ST tumor regions. Color represents scaled expression within each tumor class. Bracket denotes constituent genes in IEX. (G) Boxplots of IEX scores in TCGA COAD and READ samples, stratified by MMR status (MSS n = 301; MSI-H n = 45). **** p <1.00 × 10⁻⁴ by Student’s t-test with Bonferroni correction. (H) Kaplan-Meier PFS curves for TCGA COAD and READ samples from B with high (+) and low (−) IEX scores. (I-J) MxIF images showing epithelial (left) and immune (right) markers from representative (I) IEX+ and (J) IEX- cores from the CRC TMA. Scale bars 100 μm. (K) Pie chart showing the breakdown of MMR status in the CRC TMA (top) and UpSet plot comparing CRC TMA cores with high-scoring DDR1 and/or TGFBI IHC staining to cores with low-scoring stains for both markers (bottom; total n = 166). (L) Kaplan-Meier PFS curves for CRC TMA cores with high (+) and low (−) IHC staining of both DDR1 and TGFBI. (M) Percentages of each growth rate classification in control and DDR1-reduced (DDR1r) syngeneic tumors. (N) Plots of tumors volumes in control and DDR1-reduced groups from M at day 31 (n=10 per group). * p <5.00 × 10⁻² by 2-sided Student’s t-test. See also Figure S6 and Tables S2-S4.