Spatially defined multicellular functional units in colorectal cancer revealed from single cell and spatial transcriptomics

Abstract

While advances in single cell genomics have helped to chart the cellular components of tumor ecosystems, it has been more challenging to characterize their specific spatial organization and functional interactions. Here, we combine single cell RNA-seq, spatial transcriptomics by Slide-seq, and in situ multiplex RNA analysis, to create a detailed spatial map of healthy and dysplastic colon cellular ecosystems and their association with disease progression. We profiled inducible genetic CRC mouse models that recapitulate key features of human CRC, assigned cell types and epithelial expression programs to spatial tissue locations in tumors, and computationally used them to identify the regional features spanning different cells in the same spatial niche. We find that tumors were organized in cellular neighborhoods, each with a distinct composition of cell subtypes, expression programs, and local cellular interactions. Comparing to scRNA-seq and Slide-seq data from human CRC, we find that both cell composition and layout features were conserved between the species, with mouse neighborhoods correlating with malignancy and clinical outcome in human patient tumors, highlighting the relevance of our findings to human disease. Our work offers a comprehensive framework that is applicable across various tissues, tumors, and disease conditions, with tools for the extrapolation of findings from experimental mouse models to human diseases.

Figure 1.. A single cell atlas of healthy colon and dysplastic lesions in mouse.

A. Study overview. B. Major cell subsets of healthy colon and dysplastic lesions. 2D embedding of 48,115 single cell profiles colored by cluster (top, legend) or annotated cell type (bottom, legend). C,D. Changes in cell composition in dysplastic tissues. C. 2D embedding of single cell profiles, showing only the cells in each condition state, subsampled to equal numbers of cells per condition state, colored by cluster (same legend as in B (top)). D. Proportion of cells (y axis) of each cell type in each sample (x axis). E. multiplex RNA in situ analysis. Representative images of Cartana analysis of normal colon (left) and AV lesions (right) colored by cell type assignment (same as Supplementary item 2A).

Figure 2.. Composition and cell intrinsic expression program changes in dysplastic epithelial cells.

A-C. Compositional changes in epithelial cells in dysplastic tissue. A. 2D embedding of epithelial cell profiles colored by clusters (legend). Cluster Epi06: potential doublets (Methods). B. Proportion of cells out of all epithelial cells (y axis) of each epithelial cell subset in each sample (x axis). C. Fraction of expressing cells (dot size) and mean expression in expressing cells (dot color) of marker genes (columns) for each cluster (rows). D-E. Use of epithelial cell programs changes in dysplastic tissue. D. Weights (x axis) of each of the 20 top ranked genes (y axis) for each program. E. Proportion of program weights summed over all epithelial cells (y axis) in each sample (x axis). F. Stem cell program 16 is induced in epithelial cells in dysplastic tissue. Scaled log-normalized expression (color bar) of the top 100 genes differentially expressed between cells from normal colon and from dysplastic (AV and AKPV) across the 10,812 cells that accounted for 90% of program 16’s expression across all epithelial cells (columns). Selected program genes are marked.

Figure 3.. Altered cell type neighborship in CRC.

A. Cell type distributions in situ. Slide-seq pucks of dysplastic (top) and normal (bottom) tissue colored by TACCO assignment of cell labels (legend, light grey: low quality beads) (x and y axis: spatial coordinates in μm). B,C. Cell type neighborships in normal and dysplastic colon tissue. Short-range (up to 20μm) neighborship enrichment (Z score, color bar) vs. a background of spatially random annotation assignments for each pair of cell annotations (rows, columns) in normal (A) and dysplastic (B) tissue.

Figure 4.. Three cellular neighborhoods associated with tumor progression.

A. Spatial regions. Slide-seq pucks of AV (top) and normal (bottom) mouse colon colored by TACCO regions (legend, light gray: low quality beads) (x and y axis are spatial coordinates in μm). B,C. Enrichment and depletion of cell subsets and epithelial programs across different regions. Significance (FDR, color bar, two-sided Welch’s t-test on CLR-transformed compositions) of enrichment (red) or depletion (blue) of specific cell subsets (rows, B) or epithelial cell programs (rows, C) in the different regions defined by TACCO (columns) as well as all normal (“N (ref.)”, leftmost column) and AV (“AV (ref.)”, leftmost column) samples. D. TACCO defined regions preferentially relate to normal or AV tissue. Significance (FDR, color bar, two-sided Welch’s t-test on CLR-transformed compositions) of enrichment (red) or depletion (blue) of each TACCO defined region (rows) in normal (“N vs. rest”) and AV (“AV vs. N”) samples (columns). E. TACCO reveals normal colon architecture. Left: Slide-seq puck of normal mouse colon colored by TACCO region annotations (legend) (x and y axis: spatial coordinates (μm)). Right: Main epithelial expression programs enriched in each region (FDR<6.3 10⁻⁴, two-sided Welch’s t test on CLR-transformed compositions) except region 2 (muscularis), which is characterized by non-epithelial (stromal) cell types. F. Expression signatures of cells in normal regions 3,5,10 and 12. Scaled log-normalized expression of the top 20 differentially expressed genes (rows) for each bead (columns) in the region. G,H. Malignant-like regions. G. Slide-seq pucks of two AV lesions colored by TACCO annotations of malignant-like regions 6, 8 and 11. H. Scaled log-normalized expression of the top 20 differentially expressed genes (rows) of each bead (left, columns) in the region; or epithelial (middle left), immune (middle right) or stromal (right) fractions of beads (columns) in regions 6, 8 and 11 in dysplastic lesions. I,J. Epithelial cell subsets and programs associated with “malignant-like”, “normal-like” and normal tissues. Significance (FDR, color bar, two-sided Welch’s t-test on CLR-transformed compositions) of enrichment (red) or depletion (blue) of epithelial cell programs (I, rows) or epithelial, immune and stromal cell subsets (J, rows) in different tissue types (columns) based on Slide-seq or scRNA-seq (“sc”) samples. K. Inferred interaction pathways. Enrichment (FDR) in AV vs. normal tissue of corresponding “sender” (x axis) and “receiver” (y axis) (aggregated over ligand-receptor pairs) pathways (dots). Red/blue: pathways significantly enriched in AV/normal samples and (light red: only for either sender or receiver). The top 5 enriched pathways (in each direction) are labeled.

Figure 5.. Mouse tumor regions associated with tumor progression in human colorectal tumors.

A. Expression profiles characterizing mouse regions are recapitulated in human tumors. Significance (FDR, color bar, two-sided Welch’s t-test on CLR-transformed compositions) of enrichment (red) or depletion (blue) of region-associated epithelial, immune or stromal profiles (rows) compared between normal, MMRp, or MMRd samples (columns). B,C. Mouse regions capture malignant features in human tumors. B. Top left: First (PC1, x axis) and second (PC2, y axis) principal components of mouse region scores of mouse and human epithelial pseudo-bulk samples. Top right: PC1 loadings (x axis) of each mouse region score (y axis). Bottom: PC1 values (box plots show mean, quartiles, and whiskers for the full data distribution except for outliers outside 1.5 times the interquartile range (IQR)) for each type of mouse or human sample (x axis). C. Significance (FDR, color bar, two-sided Welch’s t-test) of enrichment (red) or depletion (blue) of region-associated profile scores (rows) in normal and dysplastic samples (columns) in human or mouse. D,E. Expression of malignant like regions 6 and 11 in tumors are associated with PFI (D) and OS (E) in human patients. Kaplan-Meier PFI (e, n = 662) or OS (f, n = 662) analysis of human bulk RNA-seq cohort stratified by malignant-like region profile scores.