Multiplexed 3D atlas of state transitions and immune interaction in colorectal cancer

Abstract

Advanced solid cancers are complex assemblies of tumor, immune, and stromal cells characterized by high intratumoral variation. We use highly multiplexed tissue imaging, 3D reconstruction, spatial statistics, and machine learning to identify cell types and states underlying morphological features of known diagnostic and prognostic significance in colorectal cancer. Quantitation of these features in high-plex marker space reveals recurrent transitions from one tumor morphology to the next, some of which are coincident with long-range gradients in the expression of oncogenes and epigenetic regulators. At the tumor invasive margin, where tumor, normal, and immune cells compete, T cell suppression involves multiple cell types and 3D imaging shows that seemingly localized 2D features such as tertiary lymphoid structures are commonly interconnected and have graded molecular properties. Thus, while cancer genetics emphasizes the importance of discrete changes in tumor state, whole-specimen imaging reveals large-scale morphological and molecular gradients analogous to those in developing tissues.

Keywords: 3D microscopy; PD1-PDL1 interaction; cellular; colorectal cancer; intermixed molecular; large-scale; morphological features; multiplexed imaging; spatial gradients; spatial proteomics; spatial transcriptomics; tertiary lymphoid structures; tumor atlas; tumor budding.

Figure 1.. Data overview.

(A) Data collection strategy - 93 CRC specimens available as 3D stack, single whole-slides, and TMAs. (B) Histopathologic annotation of six ROIs and three invasive margins (A: budding, B: mucinous, C: pushing) on H&E (left). Representative images of ROIs (center). Schematic diagram of architectural features (right). (C) CyCIF whole-slide image and cell-type assignment. 21 cell types from 3 main categories (tumor, stroma, and immune; Table S4) were defined and locations mapped. (D) Comparison of cell-type percentages by scRNA-seq and CyCIF. (E) t-SNE of single-cell data (CRC1/097) generated using all markers; 50,000 randomly-sampled cells displayed. Cell-type plot (right) color code same as Figure S1C.

Figure 2.. Spatial heterogeneity and estimation errors for regional sampling.

(A) Length scales for select markers across CRC1–17. (B) Spatial correlations of binarized staining intensities for CK⁺ (red), α-SMA⁺ (blue), and FOXP3⁺ (green) cells, and exponential fits. (C) CyCIF image showing CD20⁺ TLS (pink circle) and CD31⁺ blood vessel (yellow circle). (D) Spatial distribution of CD20⁺ cells (magenta dots, contours) and CD31⁺ cells (cyan dots); #1–6: annotated ROIs. (E) Virtual TMA cores from CRC1/097 and real TMA cores from CRC2–93. (F) Cell-type abundance estimates using vTMA cores or random sampling. (G) Estimation error of vTMAs summarized by fold-reduction in effective sample size, N/N_eff, for marker log-intensities and cell-type compositions. (H) Correlation of select cell-type pairs amongst 10 nearest neighbors (I) Correlation functions of CK⁺ cells, estimated from vTMAs or random sampling. Estimates from four cores also shown. (J) Images of cores highlighted in (I). (K) Fraction of marker-positive cells across CRC2–17 whole-slide or TMA data, or TMAs from CRC18–93. Box plot displays data points and 1^st-3^rd quartiles, whiskers extend at most to 1.5x interquartile range, and proportions <0.0001 are denoted as a single data point along dotted line. Outliers labeled as crosses (F&H) or circles (A&K); medians are indicated.

Figure 3.. Correlation and prediction of morphologic and molecular tumor phenotypes.

(A) Example ROIs corresponding to four tumor morphologies used for training and non-adjacent regions predicted with high confidence. kNN classifiers were trained and validated separately for each section to evaluate model reproducibility. (B) Prediction confidence for assignment of kNN classes as measured by Shannon entropy (0 corresponds to perfect certainty; 2 indicates random assignment (equal mixing). (C) Posterior probability that each CK⁺ cell belongs to the given tumor class. Annotation reflects classifier gradients corresponding to morphologic phenotype. (D) Left: Sample tumor region that transitions from normal to abnormal glandular features coinciding with transition from E-cadherin expression to PCNA (CyCIF, bottom). Contours describe averaged local epithelial cell expression of PCNA. Center and right: Additional examples of transition regions. (E) PCA of 31 spatially-resolved GeoMx transcriptomics regions (areas in Figure S1A). (F) Cumulative distribution of single-cell classification entropy of CRC1–17. Patients with only two classes had only normal epithelial and a tumor morphology class. Different CRC1 sections used different markers for classification. (G) Examples of marker gradients; whole tumor sections. White circles denote TMA cored regions.

Figure 4.. Tumor budding is a distributed phenomenon associated with graded molecular and morphologic transitions.

(A) Left: H&E FOV from CRC1/096 IM-A (Figure 1B); budding cells indicated by boxes/arrowheads. Right: Corresponding CyCIF (CRC1/097). Outlines indicate main tumor mass (red) and canonical tumor buds (yellow). (B) Different magnifications of annotated budding region (CRC1/097). (C) CRC1 IM-A 3D overview. Left: Surface renderings of glandular tumor (blue), α-SMA⁺ stroma (purple), normal mucosa (green), CD68⁺PDL1⁺ cells (yellow), budding cells (red). Right: All annotated buds colored by budding cell density showing interconnected fibril-like networks of budding cells. (D) 3D visualization of annotated buds (purple) relative to connected tumor mass (gray) and cells with uncertain connectivity (green). Corresponding regions in 2D images shown in Figure S4B. (E) Delaunay clusters of CK⁺ cells in a local FOV (CRC1/097). CK⁺ cell neighborhoods are denoted by edges, along with CK⁻ cells (blue) and pathology annotated buds (white). (F) Cluster sizes (log2) in CRC1. Left: Histogram across all 25 sections. Right: Mapped onto section 097. (G) Left: t-SNE of cluster size. Color represents log2 cluster size; black outline denotes small clusters (including annotated buds). Center and right: t-SNE of CK⁺ cell expression of indicated marker intensity. (H-I) Marker intensity and cluster-size. Annotated buds in green. Box plots show 1^st-3^rd quartiles; points beyond not shown. Each box represents ~10⁵-10⁶ tumor cells.

Figure 5.. Small, isolated tumor and mucin structures in 2D are large, connected networks in 3D.

(A) Example transition from main tumor mass into fibrils and ‘bud-like’ cells in stroma; CyCIF (top), H&E (bottom). Na-K ATPase and PCNA decrease with cluster size from main tumor mass to fibril tips (arrows, budding cells). Image oversaturated for visualization. (B) Analogous budding structures in mucinous tumor regions, with fibrils and budding cells (arrowheads) extending into mucin pools. (C) GeoMx data heatmap for selected EMT hallmark genes. Columns correspond to analyzed region from one tissue section (Figure S1A); morphology indicated. (D) Two H&E FOVs from different regions of reconstructed mucin structure with apparently isolated pools in 2D sections (arrowheads). (E) Connectivity of mucin pools across serial sections. Largest contiguous mucin network (red) extends to lumen surface (yellow outline). Image mirrored along Z relative to Figure 1B. (F) Schematic depicting serial sectioning through fibrils at invasive margin, illustrating contiguous 3D structures appearing as isolated cells/small clusters in 2D.

Figure 6.. 3D TLS structure and cell compositions.

(A) 2D TLS domains (CRC1/097); numbers indicate individual TLS/SLO domains in this section. (B) 3D rendering of TLS networks (TLSNs); CRC1. 7 largest TLSNs (A-G) - histogram shows number of individual TLS identified in 2D sections from each. (C) 3D TLSNs projected onto XY-surface. (D) TLS domain clustering by kNN (left) and number of domains in each cluster (right). (E) TLS cluster distribution in CRC1; 7 largest TLSNs are outlined/labeled. (F) Example CyCIF images of TLS clusters 1 and 3. (G) Left: 3D view of TLSN-B from CRC1 with each TLS domain colored by cluster. Right: Cross-sectional views of XY (top) and XZ (bottom) show TLS domains in TLSN-B. (H) 3D view of TLSN-B, colored by principal component 1. (I) Example CyCIF and H&E images of TLS clusters 4, 5, 6, 7. (J) TLS domain counts in CRC1–17 (section 097 for CRC1). (K) TLS cluster heatmap from CRC1–17. (L) 2D TLS domains of CRC16, colored by clusters.

Figure 7.. Immune landscape of CRC and its invasive margins.

Abundance and distribution of (A) CD45⁺, (B) CD4⁺FOXP3⁺(T_reg), CD8⁺(Tc), and (C) PDL1⁺ cells; TB (tumor budding); labels correspond to Figure 1B. (D) Co-occurrence of PDL1⁺ and PD1⁺ using a 20 μm distance cutoff. Panels A-C and K depict CRC1/097. (E) LDA topics and relative abundancies along the tumor margin. (F) PDL1 expression in indicated cell types. Top panel represents relative fractions of PDL1⁺ cells over indicated populations, while bottom panel shows absolute fractions of PDL1⁺ or double-marker positive cells. (G) Representative images of PDL1⁺CK⁺ cells in CRC1 (top) and CRC17 (lower). (H) Plot of PDL1⁺CK⁺ (top) or PDL1⁺CD68⁺ cell fractions in MSI-H or MSI-L samples from TMA data (CRC2–93). (I-J) Fraction of PDL1:PD1 interaction (20 μm) within CK⁺ (top) and CD45⁺ (bottom) cells; P-values from pairwise t-test shown (n=25). (I) In CRC1 (all 25 sections) or (J) CRC1–17 (n=17). (K) Co-occurrence maps using 20 μm distance cutoff. (L) High-resolution 3D imaging of PDL1:PD1 interaction among tumor and myeloid cells. Top: maximum intensity projections. Bottom: 3D rendering, Imaris software. (M) Schematic illustrating tumor-immune interactions at different types of invasive margins. Boxplots 25%−75% with whiskers at 5% and 95%; medians indicated. Outliers labeled crosses (F&H), circles (E).