Oncogenic signaling is coupled to colorectal cancer cell differentiation state

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.

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).

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).

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.

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).

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).

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.

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.

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.

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.

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.

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).

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.

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).

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).

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.

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

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.

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.

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.