Colorectal cancer progression is intrinsically linked to stepwise deregulation of the intestinal differentiation trajectory. In this process, sequential mutations of APC, KRAS, TP53, and SMAD4 enable oncogenic signaling and establish the hallmarks of cancer. Here, we use mass cytometry of isogenic human colon organoids and patient-derived cancer organoids to capture oncogenic signaling, cell phenotypes, and differentiation states in a high-dimensional single-cell map. We define a differentiation axis in all tumor progression states from normal to cancer. Our data show that colorectal cancer driver mutations shape the distribution of cells along the differentiation axis. In this regard, subsequent mutations can have stem cell promoting or restricting effects. Individual nodes of the cancer cell signaling network remain coupled to the differentiation state, regardless of the presence of driver mutations. We use single-cell RNA sequencing to link the (phospho-)protein signaling network to transcriptomic states with biological and clinical relevance. Our work highlights how oncogenes gradually shape signaling and transcriptomes during tumor progression.

Colorectal cancer (CRC) is one of the most prevalent neoplastic diseases. It progresses by cumulative acquisition of mutations in oncogenes and tumor suppressors, each modulating one or more intracellular signaling cascades (Fearon and Vogelstein, 1990). Activation of the Wnt/β-Catenin pathway by functional loss of adenomatous polyposis coli (APC) protein is the initiating event in a majority of CRCs and promotes a stem cell–like phenotype (Van de Wetering et al., 2002; Fearon, 2011). Constitutive activation of the EGFR/MAPK pathway components KRAS and BRAF occurs in 40 and 5–10% of CRCs, respectively, and promotes proliferation among other functions (Cancer Genome Atlas Network, 2012; Rajagopalan et al., 2002). Functional loss of TP53 occurs in more than 55% of CRCs, disrupts the DNA damage response, and is closely tied to the adenoma-carcinoma transition (Baker et al., 1990; Fearon, 2011; Mamlouk et al., 2020). TGF-β pathway function is impaired by SMAD4 loss, as a late event in the progression of 10–15% of CRCs (Fearon, 2011). It is not well investigated how these key drivers shape signal transduction network activities and differentiation states on the protein level.

Here we analyze a human colon organoid progression series encompassing driver mutations in APC, KRAS, TP53, and SMAD4 (Drost et al., 2015), using mass cytometry (MC), also known as Cytometry by time-of-flight (CyTOF). High-dimensional (phospho-)protein data shows specific effects of cell-intrinsic mutations on the signaling network which in turn anchor cells on different positions on a differentiation axis. Using single-cell RNA sequencing (scRNA-seq) as an orthogonal technique, we find that well-defined transcriptome signatures are correlated to specific cell signaling network states during CRC progression.

A differentiation trajectory links cell signaling and phenotype in normal colon organoids

The intestinal epithelium is a complex tissue consisting of a hierarchy of cell types. We employed MC, a single-cell proteomics method, to measure cell-type and cell-state markers and simultaneously quantify the activity of various intracellular signaling pathways (Bandura et al., 2009; Ornatsky et al., 2010; Spitzer and Nolan, 2016). We therefore established a custom MC antibody panel (Fig. 1 A and Table 1). Pilot experiments established that the panel is suited for simultaneously assessing cell states from stemness to differentiation, phenotypes such as proliferation or apoptosis, and key cell signaling pathway activity (Fig. 1, B and C).

Figure 1.

A CyTOF panel to analyze cell type, state, and signaling. (A) Visual representation of the antibody panel. Intracellular signaling markers in red, cell-state markers in yellow, cell-type markers in blue. IKK, IκB kinase; BMPR, BMP receptor; EMT, epithelial-mesenchymal transition. (B and C) Functional control MC experiments for selected antibodies, using split violin plots with medians. (B) HCT116 cells were subjected to small molecule stimulants and inhibitors to create positive and negative control conditions. Insulin-like growth factor (IGF; 100 ng/ml, 30 min), EGF (25 ng/ml, 10 min), AZD6244 (MEK inhibitor; 10 µM, 60 min), TNF-α (20 ng/ml, 30 min.), Neocarzinostatin (500 ng/ml, 30 min). (C) Cell lines with previously known differential expression of total proteins in question were used.

Figure 1.

A CyTOF panel to analyze cell type, state, and signaling. (A) Visual representation of the antibody panel. Intracellular signaling markers in red, cell-state markers in yellow, cell-type markers in blue. IKK, IκB kinase; BMPR, BMP receptor; EMT, epithelial-mesenchymal transition. (B and C) Functional control MC experiments for selected antibodies, using split violin plots with medians. (B) HCT116 cells were subjected to small molecule stimulants and inhibitors to create positive and negative control conditions. Insulin-like growth factor (IGF; 100 ng/ml, 30 min), EGF (25 ng/ml, 10 min), AZD6244 (MEK inhibitor; 10 µM, 60 min), TNF-α (20 ng/ml, 30 min.), Neocarzinostatin (500 ng/ml, 30 min). (C) Cell lines with previously known differential expression of total proteins in question were used.

Close modal
Table 1.

Antibody marker panel

TypeTargetVendorCloneHostLabel
Surface CD133 (PROM1) Miltenyi AC133 ms 154Sm 
Surface CD24 Fluidigm ML5 ms 169Tm 
Surface CD326 (EpCAM) BioLegend 9C4 ms 167Er 
Surface CD44 Fluidigm BJ18 ms 166Er 
Surface EphB2 BD 2H9 ms 158Gd 
Surface LGR5 Fluidigm 4D11F8 rt 161Dy 
Surface PTK7 Miltenyi 188B.12.45 ms 146Nd 
Intracellular Cleaved Caspase 3 Fluidigm D3E9 rb 142Nd 
Intracellular Cleaved PARP Fluidigm F21-852 ms 143Nd 
Intracellular IκBα Fluidigm L35A5 ms 164Dy 
Intracellular Ki-67 Fluidigm B56 ms 162Dy 
Intracellular Krt20 CST D9Z1Z rb 176Yb 
Intracellular p-p38 [T180/Y182] Fluidigm D3F9 rb 156Gd 
Intracellular p-p53 [S15] CST 16G8 ms 172Yb 
Intracellular p4e-BP1 [T37/46] CST 236B4 rb 170Er 
Intracellular pAkt [S473] Fluidigm D9E rb 152Sm 
Intracellular pChk1 [S345] CST 133D3 rb 148Nd 
Intracellular pChk2 [T68] CST C13C1 rb 141Pr 
Intracellular pERK1/2 [T202/Y204] Fluidigm D13.14.4E rb 171Yb 
Intracellular pH2A.X [S139] Fluidigm JBW301 ms 147Sm 
Intracellular pMEK1/2 [S217/221] CST 41G9 rb 151Eu 
Intracellular pNF-kB p65 [S536] CST 93H1 rb 155Gd 
Intracellular pS6 [S235/236] Fluidigm N7-548 ms 175Lu 
Intracellular pSmad1 [S463/S465]/pSmad8 [S465/S467] BD N6-1233 rt 149Sm 
Intracellular pSmad2 [S465/467]/pSmad3 [S423/425] CST D27F4 rb 153Eu 
Intracellular Vimentin CST D21H3 rb 165Ho 
Intracellular YAP CST D8H1X rb 150Nd 
TypeTargetVendorCloneHostLabel
Surface CD133 (PROM1) Miltenyi AC133 ms 154Sm 
Surface CD24 Fluidigm ML5 ms 169Tm 
Surface CD326 (EpCAM) BioLegend 9C4 ms 167Er 
Surface CD44 Fluidigm BJ18 ms 166Er 
Surface EphB2 BD 2H9 ms 158Gd 
Surface LGR5 Fluidigm 4D11F8 rt 161Dy 
Surface PTK7 Miltenyi 188B.12.45 ms 146Nd 
Intracellular Cleaved Caspase 3 Fluidigm D3E9 rb 142Nd 
Intracellular Cleaved PARP Fluidigm F21-852 ms 143Nd 
Intracellular IκBα Fluidigm L35A5 ms 164Dy 
Intracellular Ki-67 Fluidigm B56 ms 162Dy 
Intracellular Krt20 CST D9Z1Z rb 176Yb 
Intracellular p-p38 [T180/Y182] Fluidigm D3F9 rb 156Gd 
Intracellular p-p53 [S15] CST 16G8 ms 172Yb 
Intracellular p4e-BP1 [T37/46] CST 236B4 rb 170Er 
Intracellular pAkt [S473] Fluidigm D9E rb 152Sm 
Intracellular pChk1 [S345] CST 133D3 rb 148Nd 
Intracellular pChk2 [T68] CST C13C1 rb 141Pr 
Intracellular pERK1/2 [T202/Y204] Fluidigm D13.14.4E rb 171Yb 
Intracellular pH2A.X [S139] Fluidigm JBW301 ms 147Sm 
Intracellular pMEK1/2 [S217/221] CST 41G9 rb 151Eu 
Intracellular pNF-kB p65 [S536] CST 93H1 rb 155Gd 
Intracellular pS6 [S235/236] Fluidigm N7-548 ms 175Lu 
Intracellular pSmad1 [S463/S465]/pSmad8 [S465/S467] BD N6-1233 rt 149Sm 
Intracellular pSmad2 [S465/467]/pSmad3 [S423/425] CST D27F4 rb 153Eu 
Intracellular Vimentin CST D21H3 rb 165Ho 
Intracellular YAP CST D8H1X rb 150Nd 

All antibodies not purchased from Fluidigm were conjugated in-house. Channels were chosen to minimize potential spillover of high-abundance markers into low-abundance channels.

We assessed cell state heterogeneity in human normal colon organoids (NCOs), after 4 d of culture in either complete medium (+Wnt) or Wnt-deprived conditions (−Wnt). We prepared a multiplexed MC sample of two biological replicate experiments (Fig. 2 A and Figs. S1 and S2).

Figure 2.

MC reveals differentiation trajectory in NCOs . (A) Experimental workflow. Human NCOs were cultured with or without Wnt for 4 d, subsequently dissociated into single cells, barcoded, pooled, stained with the antibody panel, and measured via MC. (B) UMAP embedding of combined media conditions based on the seven cell-type markers. Color scale depicts normalized marker signals in the medium condition containing (19,921 cells) or not containing Wnt (24,501 cells). Stem cell, differentiated, apoptotic regions are marked in red. (C) Diffusion map embedding of combined media conditions based on the cell-type markers used in B. Fitted principal curve of the diffusion map space is depicted as a dashed blue line. Events downsampled to 10,000 per medium condition. (D and E) Heatmaps showing cell distributions and mean protein marker signals within equally sized bins along the principal curve defined in C. Means rescaled per marker and smoothed out using a three-bin-wide running average, 1,000 cell events per bin across conditions.

Figure 2.

MC reveals differentiation trajectory in NCOs . (A) Experimental workflow. Human NCOs were cultured with or without Wnt for 4 d, subsequently dissociated into single cells, barcoded, pooled, stained with the antibody panel, and measured via MC. (B) UMAP embedding of combined media conditions based on the seven cell-type markers. Color scale depicts normalized marker signals in the medium condition containing (19,921 cells) or not containing Wnt (24,501 cells). Stem cell, differentiated, apoptotic regions are marked in red. (C) Diffusion map embedding of combined media conditions based on the cell-type markers used in B. Fitted principal curve of the diffusion map space is depicted as a dashed blue line. Events downsampled to 10,000 per medium condition. (D and E) Heatmaps showing cell distributions and mean protein marker signals within equally sized bins along the principal curve defined in C. Means rescaled per marker and smoothed out using a three-bin-wide running average, 1,000 cell events per bin across conditions.

Close modal
+ Expand view − Collapse view
Figure S1.
Biological replicate ofFig. 2. (A) UMAP embedding of combined media conditions based on the seven cell-type markers shown. Color scale depicts normalized marker signals in ±Wnt conditions. Stem cell, differentiated, apoptotic regions are marked in red. Number of cells: +Wnt = 14,806, −Wnt = 7,632. (B) Diffusion map embedding of combined media conditions based on the cell-type markers also used in A. Fitted principal curve of the diffusion map space is depicted as a dashed blue line. Events downsampled to 10,000 per medium condition. (C and D) Heatmap showing cell distribution and mean protein marker signals within equally sized bins along the principal curve defined in B. Means rescaled per marker and smoothed out using a three-bin-wide running average. 1,000 cell events per bin across conditions. Bins containing less than 30 cells were excluded and are shown in gray.

Biological replicate of Fig. 2 . (A) UMAP embedding of combined media conditions based on the seven cell-type markers shown. Color scale depicts normalized marker signals in ±Wnt conditions. Stem cell, differentiated, apoptotic regions are marked in red. Number of cells: +Wnt = 14,806, −Wnt = 7,632. (B) Diffusion map embedding of combined media conditions based on the cell-type markers also used in A. Fitted principal curve of the diffusion map space is depicted as a dashed blue line. Events downsampled to 10,000 per medium condition. (C and D) Heatmap showing cell distribution and mean protein marker signals within equally sized bins along the principal curve defined in B. Means rescaled per marker and smoothed out using a three-bin-wide running average. 1,000 cell events per bin across conditions. Bins containing less than 30 cells were excluded and are shown in gray.

Figure S1.

Biological replicate of Fig. 2 . (A) UMAP embedding of combined media conditions based on the seven cell-type markers shown. Color scale depicts normalized marker signals in ±Wnt conditions. Stem cell, differentiated, apoptotic regions are marked in red. Number of cells: +Wnt = 14,806, −Wnt = 7,632. (B) Diffusion map embedding of combined media conditions based on the cell-type markers also used in A. Fitted principal curve of the diffusion map space is depicted as a dashed blue line. Events downsampled to 10,000 per medium condition. (C and D) Heatmap showing cell distribution and mean protein marker signals within equally sized bins along the principal curve defined in B. Means rescaled per marker and smoothed out using a three-bin-wide running average. 1,000 cell events per bin across conditions. Bins containing less than 30 cells were excluded and are shown in gray.

Close modal
+ Expand view − Collapse view
Figure S2.
Signals of non-mapping markers forFig. 2andFig. S1. (A) UMAP embedding of replicate 1 data based on seven cell-type markers. Color scale depicts normalized marker signals in ± Wnt conditions. (B) UMAP embedding of replicate 2 data based on seven cell-type markers. Color scale depicts normalized marker signals in ± Wnt conditions. Number of cells: +Wnt = 14,806, −Wnt = 7,632.

Signals of non-mapping markers for Fig. 2,andFig. S1 . (A) UMAP embedding of replicate 1 data based on seven cell-type markers. Color scale depicts normalized marker signals in ± Wnt conditions. (B) UMAP embedding of replicate 2 data based on seven cell-type markers. Color scale depicts normalized marker signals in ± Wnt conditions. Number of cells: +Wnt = 14,806, −Wnt = 7,632.

Figure S2.

Signals of non-mapping markers for Fig. 2,andFig. S1 . (A) UMAP embedding of replicate 1 data based on seven cell-type markers. Color scale depicts normalized marker signals in ± Wnt conditions. (B) UMAP embedding of replicate 2 data based on seven cell-type markers. Color scale depicts normalized marker signals in ± Wnt conditions. Number of cells: +Wnt = 14,806, −Wnt = 7,632.

Close modal

Dimension reduction using uniform manifold approximation and projection (UMAP) of seven cell-type-specific markers yielded a clustered representation of cell phenotypes (McInnes et al., 2018; Fig. 2 B). High signals of intestinal stem cell markers LGR5 and PTK7 (Barker et al., 2007; Jung et al., 2015) marked a distinct area from regions high in the differentiation marker Krt20 (Jiao et al., 2008). EphB2, known to be closely linked to the crypt–villus axis (Batlle et al., 2002), gradually decreased between the two extremes. Apoptotic cells, marked by high abundance of cleaved Caspase 3, formed clusters near the region of differentiated cells. Wnt withdrawal strongly reduced cells within the stem/crypt region, as expected (Van de Wetering et al., 2002).

To define a continuous differentiation trajectory within the dataset, we re-embedded the cell-type marker data into a diffusion map. Unlike UMAPs, which emphasize clustering of high-dimensional data, diffusion maps preserve transitions between neighboring points (Coifman and Lafon, 2006). In the space defined by diffusion components 1 and 2, cells aligned on an almost one-dimensional path (Fig. 2 C). The phenotypical regions of cells described above in Fig. 2 B were also visible in this representation, with stem cells clustering at the beginning and differentiated and apoptotic cells at the end of the inferred trajectory. As this diffusion map faithfully reproduces the intestinal epithelium’s inherent differentiation trajectory, we fitted a principal curve (blue dashed line) to define a pseudo-time axis (Hastie and Stuetzle, 1989) on which intestinal organoid cells move along toward differentiation.

We binned the data along the previously defined differentiation axis to assess the distribution of cells along this trajectory as well as the dynamics of protein markers (Fig. 2 D). In full medium (+Wnt) conditions, most cells were in the left-sided bins, representing early, stem-like phenotypes. In contrast, in (−Wnt) conditions cells shifted toward the right-sided bins that contained more differentiated cells. Protein marker distribution along the axis was rather constant, irrespective of culture conditions: Markers typically associated with stem-like and proliferating cells such as LGR5, PTK7, PROM1, and Ki-67 were tied to the left side of the pseudo-differentiation axis, while differentiation marker Krt20 had its maximum at the two-thirds waypoint and markers of apoptosis and DNA damage response, such as cleaved Caspase 3, pChk1/2, and pH2A.X, had their peak at the very end of the axis.

Multiple intracellular signaling markers were also stronger coupled to the differentiation state than to the presence or absence of Wnt (Fig. 2 E). For instance, YAP, p4e-BP1, and pS6 were mostly active in undifferentiated cells, indicating graded Hippo signaling and high metabolic activity in this region. However, activity of the MEK/ERK pathway was partly uncoupled from the differentiation axis, as Wnt withdrawal led to increased phosphorylation of both markers in undifferentiated cells. Taken together, the MC analyses show that cell signaling and cell phenotypes are graded and linked in normal colon organoids.

Cell-to-cell variance in CRC is controlled by oncogenes and the differentiation trajectory

To assess interplay between cell phenotypes, signaling, and oncogenic mutations in CRC progression, we employed a previously established set of four isogenic human colon organoid lines derived from NCOs using sequential genetic modifications (Drost et al., 2015). These lines contain the canonical CRC driver mutations APC loss, KRAS gain-of-function, TP53 loss, and SMAD4 loss. The organoids of the progression series are thus termed A, AK, AKP, AKPS, respectively. We also perturbed ligand-based signaling by growing all organoids in full medium containing Wnt and EGF, or medium lacking Wnt and/or EGF for 48 h and subsequently performed MC analysis and scRNA-seq for orthogonal validation (Fig. 3 A).

Figure 3.

CRC progression model shows non-continuous path toward stemness. (A) Experimental workflow. Human NCOs, sequentially mutated with key oncogenes of canonical CRC progression (APC loss of function, KRAS gain of function, TP53 loss of function, SMAD4 loss of function; Drost et al., 2015) were subjected to three different media conditions and processed into multiplexed scRNA-seq and MC samples. (B) Gene signature scores of scRNA-seq dataset per line and medium condition. Wnt, MAPK, TP53, and TGF-β by PROGENy (Schubert et al., 2018), BMP by Qi et al. (2017). (C) Log2 fold changes of mean marker signals relative to complete medium NCO condition. Based on multiplexed MC sample, showing all sequentially mutated lines described in A. Data in rows sampled to 10,000 cell events per condition. Full circles indicate presence of Wnt or EGF, empty circles omission of Wnt or EGF. (D) MC data binned along a pseudo-differentiation axis defined by EphB2 marker signal, normalized per medium and line. High EphB2 (stem/transit-amplifying [TA] region) on the left and low EphB2 (differentiated/apoptotic region) on the right. Data in rows sampled to 10,000 cell events per condition. 7,500 cell events per bin across conditions. Leftmost panel color scale is cells per bin, remaining panels show mean LGR5, Krt20, cleaved Caspase 3, pERK1/2, and pAkt marker signals per bin. (E) ANOVA per marker, using mutations (line), supplied growth factors (medium), differentiation (EphB2), an interaction term of mutations and growth factors, as well as the replicate ID as linear model factors. Color scale shows fraction of variance explained per predictor. Data sampled to 10,000 cell events per condition.

Figure 3.

CRC progression model shows non-continuous path toward stemness. (A) Experimental workflow. Human NCOs, sequentially mutated with key oncogenes of canonical CRC progression (APC loss of function, KRAS gain of function, TP53 loss of function, SMAD4 loss of function; Drost et al., 2015) were subjected to three different media conditions and processed into multiplexed scRNA-seq and MC samples. (B) Gene signature scores of scRNA-seq dataset per line and medium condition. Wnt, MAPK, TP53, and TGF-β by PROGENy (Schubert et al., 2018), BMP by Qi et al. (2017). (C) Log2 fold changes of mean marker signals relative to complete medium NCO condition. Based on multiplexed MC sample, showing all sequentially mutated lines described in A. Data in rows sampled to 10,000 cell events per condition. Full circles indicate presence of Wnt or EGF, empty circles omission of Wnt or EGF. (D) MC data binned along a pseudo-differentiation axis defined by EphB2 marker signal, normalized per medium and line. High EphB2 (stem/transit-amplifying [TA] region) on the left and low EphB2 (differentiated/apoptotic region) on the right. Data in rows sampled to 10,000 cell events per condition. 7,500 cell events per bin across conditions. Leftmost panel color scale is cells per bin, remaining panels show mean LGR5, Krt20, cleaved Caspase 3, pERK1/2, and pAkt marker signals per bin. (E) ANOVA per marker, using mutations (line), supplied growth factors (medium), differentiation (EphB2), an interaction term of mutations and growth factors, as well as the replicate ID as linear model factors. Color scale shows fraction of variance explained per predictor. Data sampled to 10,000 cell events per condition.

Close modal

Driver mutations affected organoid cell-signal transducers and transcriptomes. APC, KRAS, and SMAD mutations resulted in activation of Wnt and MAPK targets and downregulation of TGF-β/bone morphogenetic protein (BMP) target signatures, respectively. TP53 and SMAD4 mutations synergized in downregulation of the TP53 response (Fig. 3 B). MC could trace back many of the gene expression changes to alterations in upstream signals (Fig. 3 C). Most prominently, KRAS mutated organoids (AK, AKP, APKS) showed increased abundance of phosphorylated MEK (pMEK), phosphorylated ERK (pERK), and pS6, but notably not of pAkt and p4e-BP1. Paradoxically, EGF starvation led to an increase of pERK although expression of the MAPK target gene signature was reduced, probably owing to feedback within the signaling network. Similarly, SMAD4 loss resulted in increased phosphorylation of related co-factors SMAD1/8, although BMP target transcription was decreased.

Further effects of oncogenic mutations and growth factors were also correlated between MC and scRNA-seq data (Fig. 3 C): In NCOs, EGF starvation resulted in a decrease of the MC marker Krt20 and the enterocyte differentiation transcriptional signature while DNA damage and apoptosis markers cleaved Caspase 3, pChk1/2, and pH2A.X were increased. The effect of EGF withdrawal on DNA damage and apoptosis markers was also seen in the APC mutant line. APC loss led to increased abundances of proliferation, stem cell, and crypt markers while simultaneously reducing Krt20, as expected (Van de Wetering et al., 2002). Wnt withdrawal had no consistent effect on any progression series organoids, compatible with cell-autonomous activation of Wnt/β-Catenin signaling after APC loss (Fig. S3, A and B). Lines harboring further mutations (AK, AKP, AKPS) were rather inert to phenotypic changes in response to Wnt or EGF withdrawal. The fully mutated APKS line had highest activity of the metabolic integrator 4e-BP1, and also highest proliferative scores in the transcriptome analyses (Fig. 3 C and Fig. S3 C).

+ Expand view − Collapse view
Figure S3.
Plots ofFig. 3split by biological replicate. (A) Log2 fold changes of mean marker signals relative to complete medium NCO (replicates 1 and 2) or APCko (replicate 3) condition. Data sampled to 10,000 cell events per condition (row) and replicate. (B) Cell density distributions along a shared pseudo-differentiation axis defined across replicates by Ephrin B2 marker signal. Normalized per replicate, medium, and line. High Ephrin B2 (stem/TA region) on the left, low Ephrin B2 (differentiated/apoptotic region) on the right. Data sampled to 10,000 cell events per condition (row) and replicate. 21,000 cell events per bin across conditions and replicates. (C) Gene signature scores of scRNA-seq dataset per line and medium condition. LGR5 by Merlos-Suárez et al. (2011) and Muñoz et al. (2012), Enterocyte by Smillie et al. (2019), and S phase by Seurat (Tirosh et al., 2016).

Plots of Fig. 3 split by biological replicate. (A) Log2 fold changes of mean marker signals relative to complete medium NCO (replicates 1 and 2) or APCko (replicate 3) condition. Data sampled to 10,000 cell events per condition (row) and replicate. (B) Cell density distributions along a shared pseudo-differentiation axis defined across replicates by Ephrin B2 marker signal. Normalized per replicate, medium, and line. High Ephrin B2 (stem/TA region) on the left, low Ephrin B2 (differentiated/apoptotic region) on the right. Data sampled to 10,000 cell events per condition (row) and replicate. 21,000 cell events per bin across conditions and replicates. (C) Gene signature scores of scRNA-seq dataset per line and medium condition. LGR5 by Merlos-Suárez et al. (2011) and Muñoz et al. (2012), Enterocyte by Smillie et al. (2019), and S phase by Seurat (Tirosh et al., 2016).

Figure S3.

Plots of Fig. 3 split by biological replicate. (A) Log2 fold changes of mean marker signals relative to complete medium NCO (replicates 1 and 2) or APCko (replicate 3) condition. Data sampled to 10,000 cell events per condition (row) and replicate. (B) Cell density distributions along a shared pseudo-differentiation axis defined across replicates by Ephrin B2 marker signal. Normalized per replicate, medium, and line. High Ephrin B2 (stem/TA region) on the left, low Ephrin B2 (differentiated/apoptotic region) on the right. Data sampled to 10,000 cell events per condition (row) and replicate. 21,000 cell events per bin across conditions and replicates. (C) Gene signature scores of scRNA-seq dataset per line and medium condition. LGR5 by Merlos-Suárez et al. (2011) and Muñoz et al. (2012), Enterocyte by Smillie et al. (2019), and S phase by Seurat (Tirosh et al., 2016).

Close modal

Many of the changes observed in cell signaling and phenotypes were graded along a differentiation axis. For this analysis, we divided the dataset into 20 equally sized bins after sorting by EphB2 abundance (Fig. 3 D; first panel, blue to yellow). In line with the average EphB2 levels, NCOs predominantly occupied the bins on the right side, corresponding to more differentiated cells, while cells of the A line inhabited the left side of the axis. Subsequent KRAS mutation in AK partly reversed this effect, while AKPS cells occupied the leftmost bins of the pseudo-differentiation axis. Unexpectedly, EGF promoted differentiation prior to KRAS mutation but stemness thereafter, based on the cell distributions along the EphB2 axis. Overall, this analysis shows that different oncogenic mutations and external growth factors both shape cell differentiation states and that oncogenic mutations establish an intrinsically driven cancer signaling network and cell phenotype in a stepwise manner. These patterns could also be visualized by mean marker abundance per EphB2 bin (Fig. 3 D, panels 2–6, purple to yellow; Fig. S4). For example, LGR5 and PTK7 are graded along the EphB2 axis, while also varying greatly between cell lines. Krt20 and cleaved Caspase 3 were graded in the opposite direction, i.e., higher in EphB2-low cells, and also fluctuated between mutational and medium conditions. MAPK signal transducers, in contrast, showed a strong increase after KRAS activation and further after subsequent TP53 loss, while SMAD4 knockout boosted AKT signaling.

+ Expand view − Collapse view
Figure S4.
MC data binned along a pseudo-differentiation axis defined by Ephrin B2 marker signal. Normalized per medium and line. High Ephrin B2 (stem/TA region) on the left and low Ephrin B2 (differentiated/apoptotic region) on the right. Data sampled to 10,000 cell events per condition (row). 7,500 cell events per bin across conditions.

MC data binned along a pseudo-differentiation axis defined by Ephrin B2 marker signal. Normalized per medium and line. High Ephrin B2 (stem/TA region) on the left and low Ephrin B2 (differentiated/apoptotic region) on the right. Data sampled to 10,000 cell events per condition (row). 7,500 cell events per bin across conditions.

Figure S4.

MC data binned along a pseudo-differentiation axis defined by Ephrin B2 marker signal. Normalized per medium and line. High Ephrin B2 (stem/TA region) on the left and low Ephrin B2 (differentiated/apoptotic region) on the right. Data sampled to 10,000 cell events per condition (row). 7,500 cell events per bin across conditions.

Close modal

Based on these observations, we hypothesized that variation in this dataset could be categorized using three dimensions: introduced mutations (intrinsic signals), supplied growth factors (extrinsic signals), and differentiation phenotype. ANOVA on mean marker abundances showed that markers of stem cell or crypt progenitor phenotypes LGR5, PTK7, PROM1, CD24, CD44, YAP were all strongly linked to the EphB2 pseudo-differentiation axis and were in addition modulated by mutations (Fig. 3 E). In contrast, some signaling molecules such as pERK1/2, pMEK1/2, and pS6 were mainly explained by the mutational state, but others, including pAkt, p4e-BP1, and pSmad1/8, were also modulated along the cell differentiation trajectory. Only few markers were influenced by the presence of growth factors, most notably cleaved Caspase 3, which mainly indicated apoptosis in KRAS wild-type organoids in conditions lacking EGF. This was reflected by the highest percentage of variance explained by the interaction term between mutational state and growth factors across all markers. Across most markers, variance explained by replicates was generally very low, indicating a high reproducibility. Overall, both the ANOVA and inspection of single markers show that most signaling states in CRC cells are determined by mutation and differentiation state, and are largely decoupled from external signals provided by growth factors such as EGF and Wnt.

High-level integration of organoid transcriptome and protein data

We wanted to correlate main features of the transcriptome to the MC protein measurements, using the organoid progression series datasets. In a UMAP, single-cell transcriptomes separated by driver mutations and external growth factors (Fig. 4 A). We used non-negative matrix factorization (NMF) to de-composite the transcriptome data into six modules representing higher order gene expression information and showing different strengths across regions of the UMAP (Fig. 4 B). NMF module 1 was related to transcriptomes of NCO cells, whereas the modules 5 and 6 marked areas within the fully mutated AKPS line transcriptome clusters. We now quantified activities of well-characterized and clinically relevant transcriptional signatures in the NMF modules (Fig. 4 C). We considered the Broad Institute Hallmark signatures (Subramanian et al., 2005), the Progeny signatures defining transcriptional footprints of signaling pathways (Schubert et al., 2018), and signatures defining colon and CRC cell heterogeneity (Schwitalla et al., 2013; Mustata et al., 2013; Cañellas-Socias et al., 2022; Gregorieff et al., 2015; Wang et al., 2018; Serra et al., 2019; Sato et al., 2010; Barriga et al., 2017; Smillie et al., 2019; Muñoz et al., 2012; Joanito et al., 2022; Cortina et al., 2017; Álvarez-Varela et al., 2022; Qi et al., 2017; Ayyaz et al., 2019; Merlos-Suárez et al., 2011). As expected for a module highest in NCO organoids, NMF module 1 scored high for signatures of normal colon (Joanito et al., 2022), differentiated cell types (Smillie et al., 2019; Merlos-Suárez et al., 2011), and differentiation-associated BMP signaling (Qi et al., 2017), among others. NMF modules 2 and 3, related to transcriptomes of the AK and AKP lines, were defined by strong expression of oncogenic signaling pathway target genes, most prominently Wnt, MAPK, P53, and YAP. The non-canonical and KRAS/BRAF-driven CRC progression signature iCMS3 was distributed across the NMF 1–3 modules, clustering closely with a EGFR-driven regenerative cell signature (Gregorieff et al., 2015) and a signature for high-relapse probability epithelial CRC cells (Cañellas-Socias et al., 2022). Finally, the NMF 4–6 modules, highest in transcriptomes of cells with the full AKPS set of canonical drivers, were related to Ki67/E2F proliferative signatures (Cortina et al., 2017), Myc signaling, the iCMS2 canonical progression signature (Joanito et al., 2022), and EPHB2 gene expression. The transcriptome data also could be broken down by NMF module and signature strength in the single-cell transcriptomes along the progression series (Fig. 4 D), further confirming the relationships between oncogenic mutations and NMF modules, but also highlighting cell-to-cell variation.

Figure 4.

NMF analysis integrates single cell transcriptome and protein data. (A) UMAP representation of scRNA-seq dataset. Organoid lines and medium conditions as previously shown for MC data. (B) NMF module scores on scRNA-seq UMAP. (C) Correlation of NMF module with gene signature scores related to colon and CRC cell heterogeneity. (D) Single-cell NMF module scores grouped by organoid line. Color scale is by gene signature scores or EHPB2 expression. (E) Correlation of NMF module scores with mean MC marker signals across all three replicates.

Figure 4.

NMF analysis integrates single cell transcriptome and protein data. (A) UMAP representation of scRNA-seq dataset. Organoid lines and medium conditions as previously shown for MC data. (B) NMF module scores on scRNA-seq UMAP. (C) Correlation of NMF module with gene signature scores related to colon and CRC cell heterogeneity. (D) Single-cell NMF module scores grouped by organoid line. Color scale is by gene signature scores or EHPB2 expression. (E) Correlation of NMF module scores with mean MC marker signals across all three replicates.

Close modal

Ultimately, we correlated strengths of the RNA-seq–related NMF modules and the MC protein measurements across the organoid line and growth factor matrix (Fig. 4 E). We found that NMF1 was related to apoptosis markers cleaved Caspase 3 and cleaved PARP, while signal transducers were mostly distributed across two clusters: One was defined by the correlation of pMEK1/2, pERK1/2, pS6, and YAP with NMF2–4, whereas the other cluster was defined by the strong correlation of pSmad1/8, p4E-BP1, and pAKT with NMF4–6. Taken together, our integrated analysis of gene expression and protein activity shows that the MC signals are correlated to the main axes of transcriptome heterogeneity in CRC progression and well-characterized transcriptome signatures, thus potentially allowing to infer (phospho-)protein network signaling states from transcriptomes.

Human CRC signaling network activity is also defined by differentiation state and oncogenes

Unlike our isogenic CRC progression model lines, which inherit a defined mutational profile, genotypes of individual patient-derived tumor organoids are complex and heterogeneous. To investigate if our finding that signaling is largely determined by oncogenic mutations and cell states extends to human CRC, we employed 11 well-characterized patient-derived organoid (PDO) lines (Schütte et al., 2017; Uhlitz et al., 2021). Cell-resolved MC data of the PDOs showed graded EphB2 signals in most lines. BRAF-mutant PDOs were an exception as cells had very low levels of EphB2 and thus clustered at the low end of the pseudo-differentiation axis. This finding is in line with non-canonical progression of BRAF-mutant CRCs that does not select for LGR5-positive, EphB2-high stem-like cells.

Many RAS/RAF-wild type or RAS-mutant PDO lines were somewhat susceptible to changes in growth-factor composition in the medium, suggesting that they are not fully independent of Wnt/EGF niche factors, as noted before (Álvarez-Varela et al., 2022). As a general trend, omission of growth factors resulted in a shift toward higher EphB2 levels in the susceptible lines, accompanied by higher average expression of stem-cell marker genes such as LGR5 (Fig. 5 B). However, growth factor–dependent and –independent PDO lines scored similarly in ANOVA of signaling network molecules, as activities varied mostly with the EphB2 axis and not with medium (Fig. 5 C). This analysis implies that extrinsic growth factors in the medium may determine ranges of cell differentiation states in a subset of CRC PDO lines, but signaling network activity within the individual cell is correlated most strongly with its differentiation state.

Figure 5.

Patient-derived CRC organoid signaling networks are largely defined by differentiation state. (A) Average MC marker signals normalized to marker mean for 11 patient-derived CRC organoid lines and 3 medium conditions. Full circles indicate presence of Wnt or EGF, empty circles omission of Wnt or EGF. Data in rows sampled to 5,000 cell events per condition. (B) Cell distributions per line and medium along a pseudo-differentiation axis defined by EphB2 marker signal across all 11 lines, normalized per medium and line. Data in rows sampled to 5,000 cell events per condition. (C) ANOVA per marker and line, using medium conditions and differentiation as assessed by EphB2 as linear model factors. Color scale shows fraction of variance explained per predictor. Values of all ANOVA are shown as heatmap per factor. Data sampled to 5,000 cell events per condition.

Figure 5.

Patient-derived CRC organoid signaling networks are largely defined by differentiation state. (A) Average MC marker signals normalized to marker mean for 11 patient-derived CRC organoid lines and 3 medium conditions. Full circles indicate presence of Wnt or EGF, empty circles omission of Wnt or EGF. Data in rows sampled to 5,000 cell events per condition. (B) Cell distributions per line and medium along a pseudo-differentiation axis defined by EphB2 marker signal across all 11 lines, normalized per medium and line. Data in rows sampled to 5,000 cell events per condition. (C) ANOVA per marker and line, using medium conditions and differentiation as assessed by EphB2 as linear model factors. Color scale shows fraction of variance explained per predictor. Values of all ANOVA are shown as heatmap per factor. Data sampled to 5,000 cell events per condition.

Close modal

Conclusions

Our analysis shows that the cell differentiation state was the most important determinator of cell signaling and that cancer drivers restricted permissible cell differentiation states, thus limiting the impact of extrinsic Wnt and MAPK effectors. Each of the successive mutations in APC, KRAS, TP53, and SMAD4 had a profound effect on cell differentiation as defined by the pseudo-differentiation axis. It was somewhat surprising that cells did not gradually shift toward higher stemness, but were subject to sometimes opposing effects. While the Wnt effector APC moved cells toward high expression of stem-cell markers, mutation of the MAPK signal transducer KRAS promoted differentiation. Antagonistic Wnt and MAPK activity has recently been shown to limit anti-MAPK therapy efficiency in CRC (Zhan et al., 2019; Uhlitz et al., 2021). Our data extend this antagonistic model to activities of the APC and KRAS driver mutations during CRC progression. Additional loss of TGF-β and BMP signal transducer SMAD4 reversed this KRAS-specific effect in our data. Loss of SMAD4 could thus contribute to progression of CRC by promoting cancer cell plasticity toward stemness that was restricted by prior MAPK activating mutations. Interestingly, while consecutive mutations also shift mouse CRC organoids toward stem-cell fate in a recently published MC data set (Qin et al., 2020; Fig. S5), the effect of KRAS activation on stemness appears to be different between human and mouse organoid models. These potential species differences highlight the importance of such studies in human model systems. Collectively, our results suggest that APC, KRAS, and SMAD4 mutations play different and partially opposing roles in controlling cell plasticity, which was recently recognized as a key hallmark of cancer (Hanahan, 2022).

+ Expand view − Collapse view
Figure S5.
Analysis ofQin et al. (2020)mouse CRC organoid dataset. (A) Cell density distributions along a pseudo-differentiation axis defined by CD44 marker signal. Normalized by line. High CD44 (stem/TA region) on the left and low CD44 (differentiated/apoptotic region) on the right. Number of cells: NCO = 28,590, A = 41,218, AK = 36,005, AKP = 42,405. (B) ANOVA per marker, using introduced oncogene and cell differentiation as described by CD44 signal as linear model factors. Color scale shows fraction of variance explained per predictor.

Analysis of Qin et al. (2020) mouse CRC organoid dataset. (A) Cell density distributions along a pseudo-differentiation axis defined by CD44 marker signal. Normalized by line. High CD44 (stem/TA region) on the left and low CD44 (differentiated/apoptotic region) on the right. Number of cells: NCO = 28,590, A = 41,218, AK = 36,005, AKP = 42,405. (B) ANOVA per marker, using introduced oncogene and cell differentiation as described by CD44 signal as linear model factors. Color scale shows fraction of variance explained per predictor.

Figure S5.

Analysis of Qin et al. (2020) mouse CRC organoid dataset. (A) Cell density distributions along a pseudo-differentiation axis defined by CD44 marker signal. Normalized by line. High CD44 (stem/TA region) on the left and low CD44 (differentiated/apoptotic region) on the right. Number of cells: NCO = 28,590, A = 41,218, AK = 36,005, AKP = 42,405. (B) ANOVA per marker, using introduced oncogene and cell differentiation as described by CD44 signal as linear model factors. Color scale shows fraction of variance explained per predictor.

Close modal

We show that cancer drivers limited the impact of extrinsic Wnt and EGF, both essential factors of the normal colonic tissue stem cell niche. A subset of patient-derived cell lines was still susceptible to the growth factors, indicating that CRC signaling networks are not always fully regulated by oncogenes, as noted before (Álvarez-Varela et al., 2022). Yet, omission of Wnt and/or EGF in CRC cultures never resulted in collective cell differentiation as seen in the normal colon organoids. Instead, all susceptible colon lines in our cohort upregulated EphB2 and/or other stem cell markers as a response to the loss of stem-cell niche-associated growth factors. This indicates that all CRCs are hard-wired for the maintenance of a pool of stem-like CRC cells under a wide range of extrinsic conditions.

Wnt signaling is known to be essential for maintaining intestinal stem and crypt cell states, and EphB2 is a known graded marker of crypt cells (Van de Wetering et al., 2002; Batlle et al., 2002). Ephrins, including EphB2, are thought to be Wnt targets and repressors of CRC progression (Batlle et al., 2005). However, we found that EphB2 expression remained graded also in organoid lines with engineered loss of the Wnt signal regulator APC, suggesting that signaling pathways beyond Wnt play a role in the regulation of EphB2. Interestingly, the Hippo transducer YAP correlated well with EphB2 expression and could therefore be a candidate regulator compatible with its role to regulate stem-cell fate in the colon (Ayyaz et al., 2019; Fig. 2 E). High EphB2 abundance remained associated with stem-cell markers and stem-cell-specific cell signaling network states throughout pre-cancerous and cancerous progression states in the canonical pathway but was generally low in PDOs with BRAF mutations that probably progressed via non-canonical mechanisms, shedding light on the limitations of EphB2’s ability to predict CRC relapse in patient cohorts (Merlos-Suárez et al., 2011).

We used scRNA-seq to validate the MC (phospho-)protein analysis. Integration of both data sets allowed us to correlate high-dimensional transcriptome states to MC signaling network states via NMF modules. Among further correlations, we found that transcriptional signatures of revival stem cells (Ayyaz et al., 2019), EGFR/YAP regenerative cells, (Gregorieff et al., 2015), and fetal spheroids (Mustata et al., 2013) shared high correlation with NMF module 3. This module was in the MC signaling data correlated to high activities of YAP, S6, MEK, and ERK, but low activities of 4e-BP1 and Smad1/8. Part of these correlations are backed up by functional analyses, while others predict further, yet unexplored, features of the respective cell signaling network states. Thus, integrated analyses of organoids cells encompassing a greater diversity of organoid models, experimental conditions and signaling molecules could be used to infer signaling network states of transcriptomically well-characterized and clinically relevant cancer cell subpopulations in future experiments.

Our analysis does not shed light on control of DNA replication and repair, apoptosis, and absorptive lineage allocation in CRC progression, as fluctuations of markers related to these processes, namely pChk2, pH2A.X, cleaved Caspase 3, and Krt20 remained largely unexplained in our analysis of organoids with well-defined mutations (Fig. 4, B and C). We propose that larger MC panels and more perturbations are required to uncover how these processes shape CRC progression.

Antibody panel

We designed an antibody panel targeting a selection of cell-type-specific epitopes, multiple cell-state markers, and various members of intracellular signaling cascades (Table 1 and Fig. 1 A). Most of our antibodies were conjugated in-house using Maxpar X8 labeling kits, following the manufacturer’s protocol (201300; Fluidigm).

Antibodies were functionally tested by perturbing HCT116 cells with small molecule substances or by using cell lines with known differential expression of certain protein markers (Fig. 1, B and C; Iorio et al., 2016). Signals of antibodies targeting phosphorylated Akt, S6, and 4-EBP1 (pAkt, pS6, p4-EPB1) were increased after stimulation of the IGFR/PI3K axis compared to unstimulated controls (Fig. 1 B). pMEK was increased after stimulation with EGF and simultaneous inhibition of MEK function, strengthening EGF-induced ERK to MEK signaling by eliminating a well-known negative feedback loop (Fritsche-Guenther et al., 2011). Abundance of pERK was decreased under these conditions due to upstream MEK inhibition. TNF stimulation led to degradation of IκB and an increased phosphorylation of its inhibition target NF-κB (pNF-κB). DNA damage by Neocarzinostatin increased phosphorylation of Chk1, Chk2, H2AX, and p53. Hippo pathway factor YAP and epithelial cell marker EpCAM had higher abundances in CRC line HCT116 compared to mesenchymal neuroblastoma line LAN6, as expected (Fig. 1 C). Epithelial-to-mesenchymal transition marker Vimentin was mostly absent in HCT116 cells yet strongly expressed in LAN6. Apoptosis markers cleaved Caspase 3 (Kuida et al., 1996; Woo et al., 1998) and cleaved PARP (Cohausz and Althaus, 2009) were higher in CRC organoid line OT326 than OT227, as expected (Schütte et al., 2017; Brandt et al., 2019). Cell-type markers LGR5, PROM1, CD24, EphB2, Krt20, and PTK7 all correlated with previously known RNA expression data of the cell lines compared (Fig. S6 A; Iorio et al., 2016).

+ Expand view − Collapse view
Figure S6.
MC quality controls and gating. (A) Functional control MC experiments for selected cell-state antibodies. Split violin plots with medians. Cell lines with previously known differential expressions of total proteins in question were used. (B) Comparison of selected marker values after Maxpar Fix I or Smart tube PROT1 buffer treatment in OT326 CRC organoids. (C) Representative illustration of MC gating strategy for quality assurance.

MC quality controls and gating. (A) Functional control MC experiments for selected cell-state antibodies. Split violin plots with medians. Cell lines with previously known differential expressions of total proteins in question were used. (B) Comparison of selected marker values after Maxpar Fix I or Smart tube PROT1 buffer treatment in OT326 CRC organoids. (C) Representative illustration of MC gating strategy for quality assurance.

Figure S6.

MC quality controls and gating. (A) Functional control MC experiments for selected cell-state antibodies. Split violin plots with medians. Cell lines with previously known differential expressions of total proteins in question were used. (B) Comparison of selected marker values after Maxpar Fix I or Smart tube PROT1 buffer treatment in OT326 CRC organoids. (C) Representative illustration of MC gating strategy for quality assurance.

Close modal

Organoid and cell line culture

The following organoid lines were used in this study: normal colon organoids, established from human intraoperative material (#EA4/015/13; Charité ethics committee approval), isogenic normal and CRISPR-modified human colon organoids (Drost et al., 2015), human colon organoid lines OT227, OT302 (Schütte et al., 2017) and P009T, P013T (Uhlitz et al., 2021).

Organoids were cultured in growth-factor reduced Matrigel (356230; Corning), according to previously published protocols (Sato et al., 2009; Sato et al., 2011). Culture medium was Advanced DMEM/F12 (#12634010; Gibco) supplemented with 1× GlutaMAX (#35050061; Gibco), 1× N-2 (#17502048; Gibco), 1× B-27 (#17504044; Gibco), 500 mM N-Acetylcysteine, 1 M Hepes, 100 U/ml Penicillin-Streptomycin, and 100 µg/ml Primocin (#ant-pm-1; Invitrogen). This base medium was further supplemented with a selection of growth factors, inhibitors, and conditioned media which varied by line (Table 2). A, AK, AKP, AKPS media were designed to maintain selection pressure toward the introduced mutations (Drost et al., 2015). R-Spondin conditioned medium was produced using HEK 293T HA-RSpo-1-FC clone 3B, according to previously published protocols (Ootani et al., 2009). Wnt was provided as conditioned medium for routine culture (Wallaschek et al., 2019) or as recombinant protein prior to measurement (100 ng/ml; Time Bioscience).

Table 2.

Custom culture medium conditions for the different lines measured

SubstanceSourceConc.NCOAAKAKPAKPSPDO
Wnt CM in-house 50%      
Rspo1 CM in-house 5%      
Noggin Peprotech #250-38 100 ng/ml   
EGF Peprotech #315-09 50 ng/ml    
Nutlin3 Cayman #10004372 1 µM     
SB 202190 Sigma-Aldrich #S7067 3 µM 
A 83-01 Tocris #2939 500 nM 
SubstanceSourceConc.NCOAAKAKPAKPSPDO
Wnt CM in-house 50%      
Rspo1 CM in-house 5%      
Noggin Peprotech #250-38 100 ng/ml   
EGF Peprotech #315-09 50 ng/ml    
Nutlin3 Cayman #10004372 1 µM     
SB 202190 Sigma-Aldrich #S7067 3 µM 
A 83-01 Tocris #2939 500 nM 

NCOs were sequentially modified with APC loss (A), KRAS gain-of-function (AK), TP53 loss (AKP), and SMAD4 loss (AKPS). x indicates the presence of the factor in the custom medium.

For medium perturbation, organoids were re-plated without disaggregation and kept in base medium supplemented with combinations of Wnt, R-spondin 1, and EGF for 2 or 4 d before measurement.

MC

Organoids were collected from culture plates and dissociated into single cells using a 1:1 mixture of TrypLE Express (12604013; Gibco) and Accutase (A1110501; Gibco), supplemented with 100 U/ml Benzonase (88700; Pierce). Subsequently, the single-cell suspension was kept in 5 µM Cisplatin/PBS solution (201064; Fluidigm) for 5 min at 37°C to stain dead cells and cell debris. After washing twice with culture medium to remove Cisplatin, cells were kept for 30 min at 37°C in their original medium supernatant, which was preserved prior to harvesting the organoids, to restore culture conditions, including any growth factors secreted by the organoids. Following this resting phase, samples were fixed for 10 min at room temperature and subsequently frozen at −80°C. To preserve surface epitopes, we used proteomic stabilizer PROT1 (NC0627333; Smart Tube, Inc.), previously shown to be compatible with MC (Gaudillière et al., 2014). We adapted a recently published protocol (Böttcher et al., 2019), by mixing PROT1 1:1.4 with 10% PBS/BSA. MC signals of several key surface and intracellular markers were considerably higher compared to cells fixed with Maxpar Fix I buffer (201065; Fluidigm; Fig. S6 B). 1 d prior to measurement, samples were thawed and barcoded using the Cell-ID 20-Plex Pd Barcoding Kit (201060; Fluidigm) according to manufacturer’s instructions and subsequently pooled.

For patient-derived organoid samples, we adapted a previously described alternative multiplexing protocol (Sufi et al., 2021). Processing began 2 d prior to measurement with in situ fixation in 2% formaldehyde (28906; Pierce), followed by an overnight barcoding step at 4°C using tellurium-based reagents (Willis et al., 2018). On the following day, organoids were pooled and jointly dissociated in an enzymatic cocktail of Dispase II (17105041; Gibco), Collagenase I (17100017; Gibco), Collagenase IV (17104019; Gibco), DNAse I (A3778; PanReac AppliChem), and Benzonase (88700; Pierce) combined with mechanical disruption via a gentleMACS Octo Dissociator (130-096-427; Miltenyi Biotec) and corresponding C tube (130-093-237; Miltenyi Biotec). Remaining cell clusters were filtered out using a 70 µm cell strainer.

Multiplexed cell pools were incubated with metal-isotope tagged antibodies in a two-step protocol. First, antibodies targeting cell-surface proteins (Table 1) were added to the cell pool for 30 min at room temperature, followed by 10 min fixation in 2% formaldehyde (28906; Pierce) at room temperature and membrane permeabilization in ice-cold methanol. Subsequently, antibodies targeting intracellular and nuclear proteins were added to the cell pool for 30 min at room temperature. In between each step, cells were washed in Maxpar Cell Staining Buffer (201068; Fluidigm). DNA staining was performed using Cell-ID Intercalator-Ir (201192a; Fluidigm) with a final concentration of 62.5 nM in PBS for 20 min at room temperature. After overnight fixation in 2% formaldehyde (28906; Pierce) at 4°C, stained cells were washed once with Cell Staining Buffer, twice with doubly distilled water, and filtered through a 30 µm cell strainer.

EQ Four Element Calibration Beads (201078; Fluidigm) were added 1:10 to the cell suspension prior to analysis in a Helios mass cytometer (Fluidigm). The instrument’s software was used to normalize measured data for machine-related variance based on the added calibration beads. Further data processing and plotting was performed in R (R Core Team, 2020) using various packages of the Tidyverse family (Wickham et al., 2019) and visually adjusted in Inkscape (Inkscape Project). We used the CATALYST package for de-convolution of barcoded samples and spillover compensation (Chevrier et al., 2018). Further pre-processing steps involved filtering out EQ bead events via the bead-exclusive 140Ce channel, gating for single cells via event length and DNA parameters, and excluding dead cells via a platinum-iridium gate (Fig. S6 C). All MC signals were arsinh transformed to keep zero values while benefiting of roughly normal-distributed signal distributions in logarithmic space.

Biological replicates of the CRC progression series experiment were measured in separate MC runs. We computed scaling factors to match between runs the 85th percentile of each antibody signal across all perturbation conditions measured, eliminating technical variation while preserving treatment-related effects. We sampled equal numbers of cells per perturbation condition to ensure that the adjustment was not dominated by any specific perturbation. Run 2 was adjusted using run 1 as reference, run 3 was adjusted using the normalized run 2 as reference.

MC data was normalized to account for unwanted covariance across measured protein markers that may have arisen from factors such as physical cell size or non-specific staining. For this, we fitted a multivariate linear model for each channel using four surrogates of unwanted covariance (namely, Cell-ID palladium barcode, DNA staining, total-ERK, and pan-Akt) and the dependency was regressed out. This normalization process was carried out independently for each patient-derived organoid, while preserving the average intensity of each channel

scRNA-seq

Organoids were collected from culture plates and debris removed through a process of mechanical stimulation and washing in PBS. Dissociation involved a mixture of TrypLE Express (12604013; Gibco) and Accutase (A1110501; Gibco) as well as 100 U/ml of DNAse I (PanReac, A3778). The reaction was quenched in ice-cold culture medium + 0.5% BSA and cells were filtered through a 40 µm cell strainer. Experimental conditions were tagged using the antibody-based Human Single-Cell Multiplexing Kit (633781; Becton Dickinson) according to the manufacturer’s instructions. Subsequently, cells were fixed in an ice-cold mix of 20% PBS and 80% Methanol for 15 min and stored at −80°C. For library preparation, equal amounts of cells per condition were rehydrated in PBS + 1% BSA + 0.5 U/ml RNAse inhibitor, filtered through a 40 µm cell strainer, and processed for single-cell sequencing as two pooled samples using a 10× Chromium Controller and Chromium Next GEM Single Cell 3′ v3.1 kits (1000269) according to manufacturer’s instructions. Sequencing was performed on an Illumina NovaSeq 6000 system. Raw sequencing data was processed using Cell Ranger (version 6.1.1, 10× Genomics). Processing of the filtered count matrices was done in Seurat (version 4.2.0; Hao et al., 2021) in R (version 4.2.1, R Core Team, 2020). Counts of sample tags were used to de-multiplex samples by first normalizing the counts using centered log ratios and subsequently applying HTODemux. Count data of both libraries was then combined and filtered by excluding cells with less than 2,000 detected genes and more than 20% of reads from mitochondrial transcripts. Data was normalized using the default parameters of Seurat. The 2,000 most variable transcripts were scaled and used for principal component analysis. UMAP representation was calculated based on the top 10 PCs with default parameters. Signatures were scored using the AddModuleScore method of Seurat, PROGENy scores using the progeny library (Schubert et al., 2018) with default parameters. NMF components were calculated using the R package NMF (Gaujoux and Seoighe, 2010) based on the 2,000 most variable genes, after scaling and setting negatively scaled values to zero. An initial calculation of 3–15 components was used to determine the optimal number of NMF components using the co-phenetic correlation coefficient as recommended by the package documentation. Ultimately, six components were used.

Color scales

Throughout this manuscript, we use various unique color scales of the viridis and scico packages to depict different analyses types. CyTOF marker signals are depicted on a black-red-yellow scale (“inferno”), cell distributions in pseudo-time on a dark blue to yellow scale (“viridis”), and ANOVA-based plots use a white to blue scale (“oslo”). All scales use a high contrast range for better readability and are color-blind friendly.

Online supplemental material

Figs. S1 and S2 show biological replicate data and additional markers of MC experiments. Figs. S3 and S4 show additional analyses and alternative visualizations of scRNA-seq and MC experiments performed with the human CRC progression organoid series, Fig. S5 shows a re-analysis of mouse organoid MC data, originally published by Qin et al. (2020), and Fig. S6 describes functional controls for the MC panel, cell processing and MC data acquisition.

All data used in this study and data analysis scripts are available on Zenodo: https://doi.org/10.5281/zenodo.6400082.

The authors thank the Berlin Institute of Health Flow & Mass Cytometry Core Facility, Gudrun Kliem (Robert Koch-Institute), and In-Fah M. Lee (Clinical Physiology at Charité) for their excellent help and technical assistance, Mark Nitz and Yong Jia Bu (Department of Chemistry at University of Toronto, Toronto, Canada) for providing TeMal barcoding reagents.

We acknowledge funding by the German Ministry of Education and Research, projects ZiSSTrans (02NUK047E) and MSTARS (161 L0220A, 16LW0239K), as well as Deutsche Forschungsgemeinschaft via the graduate programme CompCancer (RTG2424) and Sachbeihilfe MO 2783/5 (to M. Morkel).

Author contributions: Conceptualization: N. Blüthgen, M. Morkel, T. Sell; Investigation: T. Sell, C. Klotz, S. Krug; Formal Analysis and Methodology: T. Sell, M.M. Fischer, R. Astaburuaga-García, N. Blüthgen; Validation: T. Sell; Visualization: T. Sell; Resources: J. Drost, H. Clevers, C. Klotz; Funding Acquisition: N. Blüthgen, M. Morkel; Supervision: N. Blüthgen, M. Morkel, C. Sers; Writing—original draft: T. Sell, M. Morkel; Writing—editing and review: T. Sell, M. Morkel, N. Blüthgen.

Álvarez-Varela
,
A.
,
L.
Novellasdemunt
,
F.M.
Barriga
,
X.
Hernando-Momblona
,
A.
Cañellas-Socias
,
S.
Cano-Crespo
,
M.
Sevillano
,
C.
Cortina
,
D.
Stork
,
C.
Morral
, et al
.
2022
.
Mex3a marks drug-tolerant persister colorectal cancer cells that mediate relapse after chemotherapy
.
Nat. Cancer
.
3
:
1052
1070
.
Ayyaz
,
A.
,
S.
Kumar
,
B.
Sangiorgi
,
B.
Ghoshal
,
J.
Gosio
,
S.
Ouladan
,
M.
Fink
,
S.
Barutcu
,
D.
Trcka
,
J.
Shen
, et al
.
2019
.
Single-cell transcriptomes of the regenerating intestine reveal a revival stem cell
.
Nature
.
569
:
121
125
.
Baker
,
S.J.
,
A.C.
Preisinger
,
J.M.
Jessup
,
C.
Paraskeva
,
S.
Markowitz
,
J.K.
Willson
,
S.
Hamilton
, and
B.
Vogelstein
.
1990
.
p53 gene mutations occur in combination with 17p allelic deletions as late events in colorectal tumorigenesis
.
Cancer Res.
50
:
7717
7722
.
Bandura
,
D.R.
,
V.I.
Baranov
,
O.I.
Ornatsky
,
A.
Antonov
,
R.
Kinach
,
X.
Lou
,
S.
Pavlov
,
S.
Vorobiev
,
J.E.
Dick
, and
S.D.
Tanner
.
2009
.
Mass cytometry: Technique for real time single cell multitarget immunoassay based on inductively coupled plasma time-of-flight mass spectrometry
.
Anal. Chem.
81
:
6813
6822
.
Barker
,
N.
,
J.H.
van Es
,
J.
Kuipers
,
P.
Kujala
,
M.
van den Born
,
M.
Cozijnsen
,
A.
Haegebarth
,
J.
Korving
,
H.
Begthel
,
P.J.
Peters
, and
H.
Clevers
.
2007
.
Identification of stem cells in small intestine and colon by marker gene Lgr5
.
Nature
.
449
:
1003
1007
.
Barriga
,
F.M.
,
E.
Montagni
,
M.
Mana
,
M.
Mendez-Lago
,
X.
Hernando-Momblona
,
M.
Sevillano
,
A.
Guillaumet-Adkins
,
G.
Rodriguez-Esteban
,
S.J.A.
Buczacki
,
M.
Gut
, et al
.
2017
.
Mex3a marks a slowly dividing subpopulation of Lgr5+ intestinal stem cells
.
Cell Stem Cell
.
20
:
801
816.e7
.
Batlle
,
E.
,
J.T.
Henderson
,
H.
Beghtel
,
M.M.W.
van den Born
,
E.
Sancho
,
G.
Huls
,
J.
Meeldijk
,
J.
Robertson
,
M.
van de Wetering
,
T.
Pawson
, and
H.
Clevers
.
2002
.
Beta-catenin and TCF mediate cell positioning in the intestinal epithelium by controlling the expression of EphB/ephrinB
.
Cell
.
111
:
251
263
.
Batlle
,
E.
,
J.
Bacani
,
H.
Begthel
,
S.
Jonkheer
,
A.
Gregorieff
,
M.
van de Born
,
N.
Malats
,
E.
Sancho
,
E.
Boon
,
T.
Pawson
, et al
.
2005
.
EphB receptor activity suppresses colorectal cancer progression
.
Nature
.
435
:
1126
1130
.
Böttcher
,
C.
,
S.
Schlickeiser
,
M.A.M.
Sneeboer
,
D.
Kunkel
,
A.
Knop
,
E.
Paza
,
P.
Fidzinski
,
L.
Kraus
,
G.J.L.
Snijders
,
R.S.
Kahn
, et al
.
2019
.
Human microglia regional heterogeneity and phenotypes determined by multiplexed single-cell mass cytometry
.
Nat. Neurosci.
22
:
78
90
.
Brandt
,
R.
,
T.
Sell
,
M.
Lüthen
,
F.
Uhlitz
,
B.
Klinger
,
P.
Riemer
,
C.
Giesecke-Thiel
,
S.
Schulze
,
I.A.
El-Shimy
,
D.
Kunkel
, et al
.
2019
.
Cell type-dependent differential activation of ERK by oncogenic KRAS in colon cancer and intestinal epithelium
.
Nat. Commun.
10
:
2919
.
Cañellas-Socias
,
A.
,
C.
Cortina
,
X.
Hernando-Momblona
,
S.
Palomo-Ponce
,
E.J.
Mulholland
,
G.
Turon
,
L.
Mateo
,
S.
Conti
,
O.
Roman
,
M.
Sevillano
, et al
.
2022
.
Metastatic recurrence in colorectal cancer arises from residual EMP1+ cells
.
Nature
.
611
:
603
613
.
Cancer Genome Atlas Network
.
2012
.
Comprehensive molecular characterization of human colon and rectal cancer
.
Nature
.
487
:
330
337
.
Chevrier
,
S.
,
H.L.
Crowell
,
V.R.T.
Zanotelli
,
S.
Engler
,
M.D.
Robinson
, and
B.
Bodenmiller
.
2018
.
Compensation of signal spillover in suspension and imaging mass cytometry
.
Cell Syst.
6
:
612
620.e5
.
Cohausz
,
O.
, and
F.R.
Althaus
.
2009
.
Role of PARP-1 and PARP-2 in the expression of apoptosis-regulating genes in HeLa cells
.
Cell Biol. Toxicol.
25
:
379
391
.
Coifman
,
R.R.
, and
S.
Lafon
.
2006
.
Diffusion maps
.
Appl. Comput. Harmon. Anal.
21
:
5
30
.
Cortina
,
C.
,
G.
Turon
,
D.
Stork
,
X.
Hernando-Momblona
,
M.
Sevillano
,
M.
Aguilera
,
S.
Tosi
,
A.
Merlos-Suárez
,
C.
Stephan-Otto Attolini
,
E.
Sancho
, and
E.
Batlle
.
2017
.
A genome editing approach to study cancer stem cells in human tumors
.
EMBO Mol. Med.
9
:
869
879
.
Drost
,
J.
,
R.H.
van Jaarsveld
,
B.
Ponsioen
,
C.
Zimberlin
,
R.
van Boxtel
,
A.
Buijs
,
N.
Sachs
,
R.M.
Overmeer
,
G.J.
Offerhaus
,
H.
Begthel
, et al
.
2015
.
Sequential cancer mutations in cultured human intestinal stem cells
.
Nature
.
521
:
43
47
.
Fearon
,
E.R.
2011
.
Molecular genetics of colorectal cancer
.
Annu. Rev. Pathol.
6
:
479
507
.
Fearon
,
E.R.
and
B.
Vogelstein
.
1990
.
A genetic model for colorectal tumorigenesis
.
Cell
.
61
:
759
767
.
Fritsche-Guenther
,
R.
,
F.
Witzel
,
A.
Sieber
,
R.
Herr
,
N.
Schmidt
,
S.
Braun
,
T.
Brummer
,
C.
Sers
, and
N.
Blüthgen
.
2011
.
Strong negative feedback from Erk to Raf confers robustness to MAPK signalling
.
Mol. Syst. Biol.
7
:
489
.
Gaudillière
,
B.
,
G.K.
Fragiadakis
,
R.V.
Bruggner
,
M.
Nicolau
,
R.
Finck
,
M.
Tingle
,
J.
Silva
,
E.A.
Ganio
,
C.G.
Yeh
,
W.J.
Maloney
, et al
.
2014
.
Clinical recovery from surgery correlates with single-cell immune signatures
.
Sci. Transl. Med.
6
:
255ra131
.
Gaujoux
,
R.
, and
C.
Seoighe
.
2010
.
A flexible R package for nonnegative matrix factorization
.
BMC Bioinformatics
.
11
:
367
.
Gregorieff
,
A.
,
Y.
Liu
,
M.R.
Inanlou
,
Y.
Khomchuk
, and
J.L.
Wrana
.
2015
.
Yap-dependent reprogramming of Lgr5(+) stem cells drives intestinal regeneration and cancer
.
Nature
.
526
:
715
718
.
Hanahan
,
D.
2022
.
Hallmarks of cancer: New dimensions
.
Cancer Discov.
12
:
31
46
.
Hao
,
Y.
,
S.
Hao
,
E.
Andersen-Nissen
,
W.M.
Mauck
III
,
S.
Zheng
,
A.
Butler
,
M.J.
Lee
,
A.J.
Wilk
,
C.
Darby
,
M.
Zager
, et al
.
2021
.
Integrated analysis of multimodal single-cell data
.
Cell
.
184
:
3573
3587.e29
.
Hastie
,
T.
, and
W.
Stuetzle
.
1989
.
Principal curves
.
J. Am. Stat. Assoc.
84
:
502
516
.
Iorio
,
F.
,
T.A.
Knijnenburg
,
D.J.
Vis
,
G.R.
Bignell
,
M.P.
Menden
,
M.
Schubert
,
N.
Aben
,
E.
Gonçalves
,
S.
Barthorpe
,
H.
Lightfoot
, et al
.
2016
.
A landscape of pharmacogenomic interactions in cancer
.
Cell
.
166
:
740
754
.
Jiao
,
Y.-F.
,
S.
Nakamura
,
T.
Sugai
,
N.
Yamada
, and
W.
Habano
.
2008
.
Serrated adenoma of the colorectum undergoes a proliferation versus differentiation process: New conceptual interpretation of morphogenesis
.
Oncology
.
74
:
127
134
.
Joanito
,
I.
,
P.
Wirapati
,
N.
Zhao
,
Z.
Nawaz
,
G.
Yeo
,
F.
Lee
,
C.L.P.
Eng
,
D.C.
Macalinao
,
M.
Kahraman
,
H.
Srinivasan
, et al
.
2022
.
Single-cell and bulk transcriptome sequencing identifies two epithelial tumor cell states and refines the consensus molecular classification of colorectal cancer
.
Nat. Genet.
54
:
963
975
.
Jung
,
P.
,
C.
Sommer
,
F.M.
Barriga
,
S.J.
Buczacki
,
X.
Hernando-Momblona
,
M.
Sevillano
,
M.
Duran-Frigola
,
P.
Aloy
,
M.
Selbach
,
D.J.
Winton
, and
E.
Batlle
.
2015
.
Isolation of human colon stem cells using surface expression of PTK7
.
Stem Cell Rep.
5
:
979
987
.
Kuida
,
K.
,
T.S.
Zheng
,
S.
Na
,
C.
Kuan
,
D.
Yang
,
H.
Karasuyama
,
P.
Rakic
, and
R.A.
Flavell
.
1996
.
Decreased apoptosis in the brain and premature lethality in CPP32-deficient mice
.
Nature
.
384
:
368
372
.
Mamlouk
,
S.
,
T.
Simon
,
L.
Tomás
,
D.C.
Wedge
,
A.
Arnold
,
A.
Menne
,
D.
Horst
,
D.
Capper
,
M.
Morkel
,
D.
Posada
, et al
.
2020
.
Malignant transformation and genetic alterations are uncoupled in early colorectal cancer progression
.
BMC Biol.
18
:
116
.
McInnes
,
L.
,
J.
Healy
, and
J.
Melville
.
2018
.
UMAP: Uniform manifold approximation and projection for dimension reduction
.
arXiv
.
Merlos-Suárez
,
A.
,
F.M.
Barriga
,
P.
Jung
,
M.
Iglesias
,
M.V.
Céspedes
,
D.
Rossell
,
M.
Sevillano
,
X.
Hernando-Momblona
,
V.
da Silva-Diz
,
P.
Muñoz
, et al
.
2011
.
The intestinal stem cell signature identifies colorectal cancer stem cells and predicts disease relapse
.
Cell Stem Cell
.
8
:
511
524
.
Muñoz
,
J.
,
D.E.
Stange
,
A.G.
Schepers
,
M.
van de Wetering
,
B.-K.
Koo
,
S.
Itzkovitz
,
R.
Volckmann
,
K.S.
Kung
,
J.
Koster
,
S.
Radulescu
, et al
.
2012
.
The Lgr5 intestinal stem cell signature: Robust expression of proposed quiescent “+4” cell markers
.
EMBO J.
31
:
3079
3091
.
Mustata
,
R.C.
,
G.
Vasile
,
V.
Fernandez-Vallone
,
S.
Strollo
,
A.
Lefort
,
F.
Libert
,
D.
Monteyne
,
D.
Pérez-Morga
,
G.
Vassart
, and
M.-I.
Garcia
.
2013
.
Identification of Lgr5-independent spheroid-generating progenitors of the mouse fetal intestinal epithelium
.
Cell Rep.
5
:
421
432
.
Ootani
,
A.
,
X.
Li
,
E.
Sangiorgi
,
Q.T.
Ho
,
H.
Ueno
,
S.
Toda
,
H.
Sugihara
,
K.
Fujimoto
,
I.L.
Weissman
,
M.R.
Capecchi
, and
C.J.
Kuo
.
2009
.
Sustained in vitro intestinal epithelial culture within a Wnt-dependent stem cell niche
.
Nat. Med.
15
:
701
706
.
Ornatsky
,
O.
,
D.
Bandura
,
V.
Baranov
,
M.
Nitz
,
M.A.
Winnik
, and
S.
Tanner
.
2010
.
Highly multiparametric analysis by mass cytometry
.
J. Immunol. Methods
.
361
:
1
20
.
Qi
,
Z.
,
Y.
Li
,
B.
Zhao
,
C.
Xu
,
Y.
Liu
,
H.
Li
,
B.
Zhang
,
X.
Wang
,
X.
Yang
,
W.
Xie
, et al
.
2017
.
BMP restricts stemness of intestinal Lgr5+ stem cells by directly suppressing their signature genes
.
Nat. Commun.
8
:
13824
.
Qin
,
X.
,
J.
Sufi
,
P.
Vlckova
,
P.
Kyriakidou
,
S.E.
Acton
,
V.S.W.
Li
,
M.
Nitz
, and
C.J.
Tape
.
2020
.
Cell-type-specific signaling networks in heterocellular organoids
.
Nat. Methods
.
17
:
335
342
.
R Core Team
.
2020
.
R: A Language and Environment for Statistical Computing
.
R Found. Stat. Comput.
Rajagopalan
,
H.
,
A.
Bardelli
,
C.
Lengauer
,
K.W.
Kinzler
,
B.
Vogelstein
, and
V.E.
Velculescu
.
2002
.
Tumorigenesis: RAF/RAS oncogenes and mismatch-repair status
.
Nature
.
418
:
934
.
Sato
,
T.
,
R.G.
Vries
,
H.J.
Snippert
,
M.
van de Wetering
,
N.
Barker
,
D.E.
Stange
,
J.H.
van Es
,
A.
Abo
,
P.
Kujala
,
P.J.
Peters
, and
H.
Clevers
.
2009
.
Single Lgr5 stem cells build crypt-villus structures in vitro without a mesenchymal niche
.
Nature
.
459
:
262
265
.
Sato
,
T.
,
J.H.
van Es
,
H.J.
Snippert
,
D.E.
Stange
,
R.G.
Vries
,
M.
van den Born
,
N.
Barker
,
N.F.
Shroyer
,
M.
van de Wetering
, and
H.
Clevers
.
2010
.
Paneth cells constitute the niche for Lgr5 stem cells in intestinal crypts
.
Nature
.
469
:
415
418
.
Sato
,
T.
,
D.E.
Stange
,
M.
Ferrante
,
R.G.J.
Vries
,
J.H.
Van Es
,
S.
Van den Brink
,
W.J.
Van Houdt
,
A.
Pronk
,
J.
Van Gorp
,
P.D.
Siersema
, and
H.
Clevers
.
2011
.
Long-term expansion of epithelial organoids from human colon, adenoma, adenocarcinoma, and Barrett’s epithelium
.
Gastroenterology
.
141
:
1762
1772
.
Schütte
,
M.
,
T.
Risch
,
N.
Abdavi-Azar
,
K.
Boehnke
,
D.
Schumacher
,
M.
Keil
,
R.
Yildiriman
,
C.
Jandrasits
,
T.
Borodina
,
V.
Amstislavskiy
, et al
.
2017
.
Molecular dissection of colorectal cancer in pre-clinical models identifies biomarkers predicting sensitivity to EGFR inhibitors
.
Nat. Commun.
8
:
14262
.
Schubert
,
M.
,
B.
Klinger
,
M.
Klünemann
,
A.
Sieber
,
F.
Uhlitz
,
S.
Sauer
,
M.J.
Garnett
,
N.
Blüthgen
, and
J.
Saez-Rodriguez
.
2018
.
Perturbation-response genes reveal signaling footprints in cancer gene expression
.
Nat. Commun.
9
:
20
.
Schwitalla
,
S.
,
P.K.
Ziegler
,
D.
Horst
,
V.
Becker
,
I.
Kerle
,
Y.
Begus-Nahrmann
,
A.
Lechel
,
K.L.
Rudolph
,
R.
Langer
,
J.
Slotta-Huspenina
, et al
.
2013
.
Loss of p53 in enterocytes generates an inflammatory microenvironment enabling invasion and lymph node metastasis of carcinogen-induced colorectal tumors
.
Cancer Cell
.
23
:
93
106
.
Serra
,
D.
,
U.
Mayr
,
A.
Boni
,
I.
Lukonin
,
M.
Rempfler
,
L.
Challet Meylan
,
M.B.
Stadler
,
P.
Strnad
,
P.
Papasaikas
,
D.
Vischi
, et al
.
2019
.
Self-organization and symmetry breaking in intestinal organoid development
.
Nature
.
569
:
66
72
.
Smillie
,
C.S.
,
M.
Biton
,
J.
Ordovas-Montanes
,
K.M.
Sullivan
,
G.
Burgin
,
D.B.
Graham
,
R.H.
Herbst
,
N.
Rogel
,
M.
Slyper
,
J.
Waldman
, et al
.
2019
.
Intra- and inter-cellular rewiring of the human colon during ulcerative colitis
.
Cell
.
178
:
714
730.e22
.
Spitzer
,
M.H.
, and
G.P.
Nolan
.
2016
.
Mass cytometry: Single cells, many features
.
Cell
.
165
:
780
791
.
Subramanian
,
A.
,
P.
Tamayo
,
V.K.
Mootha
,
S.
Mukherjee
,
B.L.
Ebert
,
M.A.
Gillette
,
A.
Paulovich
,
S.L.
Pomeroy
,
T.R.
Golub
,
E.S.
Lander
, and
J.P.
Mesirov
.
2005
.
Gene set enrichment analysis: A knowledge-based approach for interpreting genome-wide expression profiles
.
Proc. Natl. Acad. Sci. USA
.
102
:
15545
15550
.
Sufi
,
J.
,
X.
Qin
,
F.C.
Rodriguez
,
Y.J.
Bu
,
P.
Vlckova
,
M.R.
Zapatero
,
M.
Nitz
, and
C.J.
Tape
.
2021
.
Multiplexed single-cell analysis of organoid signaling networks
.
Nat. Protoc.
16
:
4897
4918
.
Tirosh
,
I.
,
B.
Izar
,
S.M.
Prakadan
,
M.H.
Wadsworth
II
,
D.
Treacy
,
J.J.
Trombetta
,
A.
Rotem
,
C.
Rodman
,
C.
Lian
,
G.
Murphy
, et al
.
2016
.
Dissecting the multicellular ecosystem of metastatic melanoma by single-cell RNA-seq
.
Science
.
352
:
189
196
.
Uhlitz
,
F.
,
P.
Bischoff
,
S.
Peidli
,
A.
Sieber
,
A.
Trinks
,
M.
Lüthen
,
B.
Obermayer
,
E.
Blanc
,
Y.
Ruchiy
,
T.
Sell
, et al
.
2021
.
Mitogen-activated protein kinase activity drives cell trajectories in colorectal cancer
.
EMBO Mol. Med.
13
:e14123.
van de Wetering
,
M.
,
E.
Sancho
,
C.
Verweij
,
W.
de Lau
,
I.
Oving
,
A.
Hurlstone
,
K.
van der Horn
,
E.
Batlle
,
D.
Coudreuse
,
A.P.
Haramis
, et al
.
2002
.
The beta-catenin/TCF-4 complex imposes a crypt progenitor phenotype on colorectal cancer cells
.
Cell
.
111
:
241
250
.
Wallaschek
,
N.
,
C.
Niklas
,
M.
Pompaiah
,
A.
Wiegering
,
C.-T.
Germer
,
S.
Kircher
,
S.
Brändlein
,
K.
Maurus
,
A.
Rosenwald
,
H.H.N.
Yan
, et al
.
2019
.
Establishing pure cancer organoid cultures: Identification, selection and verification of cancer phenotypes and genotypes
.
J. Mol. Biol.
431
:
2884
2893
.
Wang
,
Y.
,
X.
Xu
,
D.
Maglic
,
M.T.
Dill
,
K.
Mojumdar
,
P.K.-S.
Ng
,
K.J.
Jeong
,
Y.H.
Tsang
,
D.
Moreno
,
V.H.
Bhavana
, et al
.
2018
.
Comprehensive molecular characterization of the hippo signaling pathway in cancer
.
Cell Rep.
25
:
1304
1317.e5
.
Wickham
,
H.
,
M.
Averick
,
J.
Bryan
,
W.
Chang
,
L.D.
McGowan
,
R.
François
,
G.
Grolemund
,
A.
Hayes
,
L.
Henry
,
J.
Hester
, et al
.
2019
.
Welcome to the tidyverse
.
J. Open Source Softw.
4
:
1686
.
Willis
,
L.M.
,
H.
Park
,
M.W.L.
Watson
,
D.
Majonis
,
J.L.
Watson
, and
M.
Nitz
.
2018
.
Tellurium-based mass cytometry barcode for live and fixed cells
.
Cytometry A
.
93
:
685
694
.
Woo
,
M.
,
R.
Hakem
,
M.S.
Soengas
,
G.S.
Duncan
,
A.
Shahinian
,
D.
Kägi
,
A.
Hakem
,
M.
McCurrach
,
W.
Khoo
,
S.A.
Kaufman
, et al
.
1998
.
Essential contribution of caspase 3/CPP32 to apoptosis and its associated nuclear changes
.
Genes Dev.
12
:
806
819
.
Zhan
,
T.
,
G.
Ambrosi
,
A.M.
Wandmacher
,
B.
Rauscher
,
J.
Betge
,
N.
Rindtorff
,
R.S.
Häussler
,
I.
Hinsenkamp
,
L.
Bamberg
,
B.
Hessling
, et al
.
2019
.
MEK inhibitors activate Wnt signalling and induce stem cell plasticity in colorectal cancer
.
Nat. Commun
.
10
:
2197
.

Author notes

Disclosures: H. Clevers reported a patent to PCT/NL2010/000017 WO2010/090513 issued “Foundation HUB”; and “Since March 2022, I am a full-time member of the executive board of F. Hoffmann-La Roche Ltd. as head of Pharma, Research and Early Development (pRED) in Basel, Switzerland.” No other disclosures were reported.

H. Clevers’s current affiliation is Pharma, Research and Early Development of F. Hoffmann-La Roche Ltd., Basel, Switzerland.

This article is distributed under the terms of an Attribution–Noncommercial–Share Alike–No Mirror Sites license for the first six months after the publication date (see http://www.rupress.org/terms/). After six months it is available under a Creative Commons License (Attribution–Noncommercial–Share Alike 4.0 International license, as described at https://creativecommons.org/licenses/by-nc-sa/4.0/).