RNA sequencing of long-term label-retaining colon cancer stem cells identifies novel regulators of quiescence.

Recent data suggest that therapy-resistant quiescent cancer stem cells (qCSCs) are the source of relapse in colon cancer. Here, using colon cancer patient-derived organoids and xenografts, we identify rare long-term label-retaining qCSCs that can re-enter the cell cycle to generate new tumors. RNA sequencing analyses demonstrated that these cells display the molecular hallmarks of quiescent tissue stem cells, including expression of p53 signaling genes, and are enriched for transcripts common to damage-induced quiescent revival stem cells of the regenerating intestine. In addition, we identify negative regulators of cell cycle, downstream of p53, that we show are indicators of poor prognosis and may be targeted for qCSC abolition in both p53 wild-type and mutant tumors. These data support the temporal inhibition of downstream targets of p53 signaling, in combination with standard-of-care treatments, for the elimination of qCSCs and prevention of relapse in colon cancer.

Introduction

Molecular and functional intratumoral heterogeneity contributes to differences in treatment outcomes between patients with colon cancer of similar mutational profiles (Kreso et al., 2013). Studies of functional heterogeneity, as defined by phenotypic differences between cells, suggest that cancer stem cells (CSCs) are responsible for tumor growth, metastasis, and therapy resistance (De Angelis et al., 2019; Barker et al., 2008; Brock and Huang, 2017; Moore et al., 2011; O'Brien et al., 2007; Shackleton et al., 2009). CSCs share many of the characteristics of normal tissue stem cells, including unlimited self-renewal, the ability to generate differentiated daughter cells, and chemoresistance (Reya et al., 2001; Sell, 2004).

The normal intestine is maintained by highly clonogenic crypt base LGR5Positive stem cells and also contains a population of rare quiescent (G0 phase) stem cells that act as a clonogenic reserve capable of re-entering the cell cycle upon perturbation of tissue homeostasis, e.g., after injury, leading to loss of the cycling crypt base stem cells (Ayyaz et al., 2019; Barker et al., 2007; Buczacki et al., 2013; Cheng and Leblond, 1974; Clevers, 2013; Potten, 1977). Cancer often recapitulates the cellular hierarchy of the tissue in which it arises and recent evidence suggests that many tumor types contain rare slow cycling/quiescent cancer stem cells (qCSCs) (Chen et al., 2012; Dembinski and Krauss, 2009; Ebinger et al., 2016; Gao et al., 2010; Kabraji et al., 2017; Lagadinou et al., 2013; Lin et al., 2013; Pece et al., 2010; Roesch et al., 2010; Saito et al., 2010; Vanner et al., 2014; Zeuner et al., 2014). Conventional chemotherapies and radiotherapies target proliferating cells and require active cycling for induction of apoptosis (Kreso et al., 2013). In addition, cellular quiescence has been shown to facilitate immune evasion (Malladi et al., 2016). Thus, non-dividing qCSCs may escape conventional therapeutic strategies and represent the source of disease relapse after treatment (De Angelis et al., 2019; Chen et al., 2016; Moore and Lyle, 2011; Shen et al., 2020).

Cell cycle activation in qCSCs has been proposed as a therapeutic strategy to sensitize qCSCs to treatment and lead to long-term disease-free survival without relapse (Chen et al., 2016; Moore and Lyle, 2011). However, the molecular profiling of qCSCs for the identification of novel cell cycle regulators that do not also perturb cellular homeostasis in healthy tissues has been limited by both the rarity of qCSCs and the small number of suitable experimental assays available for their detection. Patient-derived organoids (PDOs) echo the morphological, differentiation, intratumor mutational, and drug sensitivity status of the original tumor (Fujii and Sato, 2021; Roerink et al., 2018; Vlachogiannis et al., 2018) and thus provide an excellent model for the prospective isolation and profiling of qCSCs.

Strategies for the identification of quiescent cells employ pulse-chase approaches, including label retention (e.g., bromodeoxyuridine, PKH26, carboxyfluorescein succinimidyl ester [CFSE]), wherein dividing cells lose the label and quiescent or slow cycling cells retain the label for an extended period of time, or the dilution of histone 2B-GFP (H2B-GFP) (Blanpain and Simons, 2013). In contrast to the H2B-GFP approach (Puig et al., 2018), which can identify transient quiescent cells, label retention allows for the identification of cells that remain quiescent from the early stages of tumorigenesis. This is important because cells selectively surviving chemotherapy have been shown to be the same cells that are quiescent/slow cycling in untreated tumors and not cells that became quiescent upon drug treatment (Francescangeli et al., 2020). Such label-retaining cells (LRCs) have previously been reported in colon cancer cell lines, xenografts, and, more recently, in PDOs (Buczacki et al., 2018; Francescangeli et al., 2012, 2020; Moore et al., 2011).

However, to date, the transcriptomic profiling of qCSCs in patients with colon cancer has been limited to microarray analyses of transiently slow-cycling H2B-GFPPositive cells from a single patient with colorectal cancer by Puig et al., (2018) and of PKH26Positive LRCs from two colon cancer patient-derived (via spheroid culture) xenograft models by Francescangeli et al., (2020). In addition, the LRCs reported in the latter study were not directly tested for proliferative or self-renewal capacity before molecular profiling and instead relied on expression of CD133 (AC133) as evidence of a stem cell phenotype.

Here, we report the identification and first whole-transcriptome RNA sequencing analyses of label-retaining qCSCs in a panel of PDOs encompassing primary colon tumors and metastases. These cells maintain a large proliferative capacity, persist long term in vivo, and display the molecular hallmarks of quiescent tissue stem cells (Cheung and Rando, 2013), including enrichment for p53 pathway and developmental gene sets alongside downregulation of cell cycle, transcription, biosynthesis, and metabolism genes. In addition, we show that qCSCs are enriched for p53-interacting negative regulators of cell cycle that we propose may be targeted for cell cycle activation and the elimination of qCSCs in both wild-type and p53 mutant cancers. These data provide a valuable resource for the development of novel therapeutic strategies geared toward the elimination of minimal residual disease and the prevention of relapse.

Results

Colon cancer PDOs contain rare label-retaining qCSCs that persist long term in vivo

PDOs faithfully recapitulate the cellular heterogeneity, architecture, somatic copy number, and mutation spectra of the original tumor and predict treatment response (Amaral et al., 2020; Choi et al., 2019; Dijkstra et al., 2021; Ganesh et al., 2019; Maru et al., 2019; Nagle et al., 2018; Ooft et al., 2019; Sachs et al., 2018; Shannon and Pitelka, 1981; Tiriac et al., 2018; Vlachogiannis et al., 2018; Weeber et al., 2015; Yao et al., 2020). They are thus widely and increasingly used as preclinical models and in translational studies, including personalized therapy design (Wensink et al., 2021). To determine whether PDOs contain non-cycling LRCs, we performed an initial 72-h pulse-chase experiment using CM-DiL dye. PDOs were established as previously described (Regan et al., 2017; Schütte et al., 2017), processed to single cells, uniformly labeled with CM-DiL dye, and seeded in Matrigel culture. CM-DiL is diluted with each cell division, halving its fluorescence between each daughter cell until it becomes undetectable. Non-cycling cells can thus be identified by their label retention. After 72 h the majority of PDO cells had lost the CM-DiL dye but some PDOs contained non-cycling LRCs (Figure 1A). To determine the frequency of these non-cycling (G0) cells we performed 5-ethynyl-2′-deoxyuridine (EdU) cell cycle analysis on a panel of colon cancer PDOs (Table S1). This analysis demonstrated that PDOs contain non-cycling cells that do not proliferate and remain in G0 within a 72 h period (Figures 1B and 1C).

Figure 1

Colon cancer PDOs contain a subpopulation of non-cycling cells

(A) Phase contrast image of colon cancer PDOs labeled with cell-tracker dye CM-DiL after 72 h (scale bar, 75 μm) (see also Table S1).

(B) Representative FACS plots of EdU cell cycle analysis of 151-ML-M PDO cells at 2 h (left hand side) and 72 h (right hand side) after labeling.

(C) Percentage of cells in G0/1, G2/M, and S phase at 2 and 72 h post EdU labeling in PDO models 151-ML-M, 162-MW-P, 195-CB-P, 249-CB-P, 278-ML-P, and 302-CB-M (±SD; data from three independent experiments).

To determine the long-term proliferative capacity of these non-cycling cells, we labeled cells with the lipophilic fluorescent dye PKH26. Unlike CM-DiL, which is suitable for short-term label retention studies, PKH26 labeling can be used to identify non-cycling cells for up to 6 months (in vitro and in vivo) (Cicalese et al., 2009; Horan et al., 1990). PDOs were dissociated to single cells, labeled with PKH26, and replated in Matrigel culture. After 12 days, PDOs were re-processed to single cells and analyzed by fluorescence-assisted cell sorting (FACS). These data demonstrated that PDOs contain rare, non-cycling, long-term LRCs (Figures 2A and 2B). Crucially, FACS isolation and replating of PKH26Positive DAPINegative (live) cells from 12-day cultures demonstrated that they are not label retaining due to terminal differentiation or senescence but can re-enter the cell cycle to generate organoids and have a large proliferative capacity (Figures 2C–2F). In addition, non-adherent spheroid formation assays, the gold standard assay for testing stem cell function in vitro (Ricci-Vitiani et al., 2007; Weiswald et al., 2015), showed that PKH26Positive cells are enriched for self-renewing CSCs (Figure 2G).

Figure 2

Non-cycling PDO cells are quiescent CSCs that can re-enter cell cycle and persist long term in vivo

(A) Representative FACS plots of PKH26-labeled 278-ML-P PDO cells after 12 h (middle panel) and 12 days (right side panel) compared with non-labeled control (left side panel).

(B) Frequency of PKH26Positive LRCs in PDO models after 12 days (±SD; data from 5 independent experiments).

(C) FACS histograms demonstrating frequency of PKH26Positive cells in 151-ML-M PDOs at 12 h (left side panel) and 12 days (middle panel) after staining and 24 days (right side panel) after FACS isolation and serial replating of PKH26Positive cells from 12-day cultures.

(D and E) (D) Phase contrast of unlabeled PDOs (negative control) (scale bar, 100 μm) and (E) immunofluorescence images of PKH26-labeled PDOs at 12 h and 12 days (middle panels) and 24 days after FACS isolation and serial re-plating of PKH26Positive LRCs from 12-day cultures (right side panel). Cells are stained for F-ACTIN (green), and nuclei are counterstained with DAPI (blue) (scale bars, 20 μm).

(F) Mean colony size of PKH26Negative and PKH26Positive cell-derived PDOs in Matrigel culture (±SD; data from three independent experiments). ∗∗p value: < 0.01 (t test).

(G) Limiting dilution spheroid formation assay of PKH26Negative and PKH26Positive cells (data from three independent experiments). The p values for pairwise tests of differences in CSC frequencies between PKH26Negative and PKH26Positive cells in 151-ML-M, 162-MW-P, 195-CB-P, 249-CB-P, 278-ML-P, and 302-CB-M tumors are 1.27 × 10−13, 1.87 × 10−5, 6.42 × 10−11, 1.12 × 10−10, 3.5 × 10−14, 6.14 × 10−12, respectively. ∗∗∗p < 0.001 (t test).

(H) Immunofluorescence image of a frozen PDX section derived from 1,000 PKH26-labeled 195-CB-P PDO cells 80 days post transplantation. Magnified region indicates a long-term label-retaining PKH26Positive cell. Cells are stained for α-tubulin (green), and nuclei are counterstained with DAPI (blue) (scale bar, 100 μm).

To test whether these cells also persisted long term in vivo we generated xenografts by transplanting PKH26-labeled cells. Long-term tracking of LRCs in xenografts requires the slow growth of the tumor. Cells were therefore transplanted at a low cell number based on knowledge of tumor growth rates from previous limiting dilution xenotransplantation assays, in which xenografts were generated from 1,000 PDO cells (Regan et al., 2017).

Unlabeled cells, lacking the burden of carrying a fluorescent dye, may be at a competitive advantage over labeled cells. Therefore, immediately before transplantation, PKH26-labeled cells were processed by FACS to exclude unlabeled cells and thus ensure that only live (DAPINegative) PKH26-labeled cells would give rise to tumors. Significantly, analysis of xenograft tissue demonstrated the presence of PKH26Positive LRCs for up to 80 days after transplantation (Figure 2H). Previous studies have observed quiescence to be a transient state (Puig et al., 2018). However, these data demonstrate that quiescence can be stable and persist long-term from the initial stages of tumor development.

RNA sequencing of PKH26Positive cells reveals the molecular signature of qCSCs

To generate a molecular profile of qCSCs we carried out RNA sequencing analyses of PKH26Negative (cycling) cells and PKH26Positive (non-cycling) qCSCs isolated from a panel of six different PDO models (Table S1) after 12 days in Matrigel culture. These data demonstrated that PKH26Positive qCSCs are enriched for stem cell-associated gene sets, such as embryonic development, organ development, placenta, nervous system development, epithelial-mesenchymal transition, Wnt, and hedgehog signaling (Figure 3A).

Figure 3

qCSCs display the molecular hallmarks of quiescent tissue stem cells, including enrichment for p53 pathway and genes common to damage-induced quiescent revSCs of the regenerating intestine

(A) RNA sequencing-generated gene set enrichment analysis for organ development (nominal p value = < 0.0005), cell development (nominal p value = < 0.0005), nervous system development (nominal p value = < 0.0005), embryonic development (nominal p value = 0.03), placenta (nominal p value = < 0.0005), epithelial-mesenchymal transition (nominal p value = < 0.0005), p53 pathway (nominal p value = < 0.0005), TNFα signaling via NF-κB (nominal p value = < 0.0005), Wnt signaling pathway (nominal p value = 0.002), and hedgehog signaling pathway (nominal p value = 0.002) in 12-day PKH26Positive LRCs (compared with PKH26Negative cells) from PDO models 151-ML-M, 162-MW-P, 195-CB-P, 249-CB-P, 278-ML-P, and 302-CB-M (n = 4 separate cell preparations).

(B) Gene ontology (GO) groups downregulated in PKH26Positive LRCs.

(C) Cell cycle, transcription, and protein synthesis GO terms downregulated in PKH26Positive LRCs.

(D) Venn diagram shows the number of upregulated RNA sequencing-generated transcripts identified in intestinal revSCs (50 genes; log fold change >0.25, p value < 0.05) by Ayyaz et al. (2019) and in PKH26Positive qCSCs (255 genes; log2 fold change >0.586, p value < 0.05) (see also Data S1) and upregulated in both revSCs and PKH26Positive qCSCs (14 genes; representation factor 21.8, p value < 1.452 × 10−15). The representation factor is the number of overlapping genes divided by the expected number of overlapping genes drawn from two independent groups. A representation factor >1 indicates more overlap than expected of the two independent groups.

(E) Table shows the 14 genes upregulated in both revSCs and PKH26Positive qCSCs. ∗ITM2C is a paralog of revSC-enriched Itm2b (see also Figure S1 and Data S1).

At the same time as showing enrichment for genes associated with growth and development, PKH26Positive qCSCs have downregulated cell cycle, transcription, protein synthesis, metabolism, and biosynthesis genes (Figures 3B and 3C). These data are in agreement with the transcriptional profiles of slow cycling/qCSCs reported in previous studies (Buczacki et al., 2018; Francescangeli et al., 2020; Puig et al., 2018) and demonstrate a common molecular signature of qCSCs.

The normal intestine contains quiescent stem cells that can regenerate the damaged intestine upon loss of crypt base stem cells following injury, although whether cellular plasticity or distinct cell types are responsible for this remains unclear. Bmi1 (Yan et al., 2012), Hopx (Takeda et al., 2011), Lrig1 (Powell et al., 2012), and Tert (Montgomery et al., 2011) have previously been reported as markers of quiescent “+4” stem cells, although subsequent studies have shown that actively cycling crypt base stem cells also express these markers at equivalent levels (Muñoz et al., 2012). Similarly, we did not detect enhanced expression of these markers in qCSCs (Data S1). This is also in agreement with a recent single-cell RNA sequencing (scRNA-seq) analyses of the regenerating mouse intestine that identified a damage-induced quiescent cell type, termed revival stem cells (revSCs) (Ayyaz et al., 2019). These cells, required for the regeneration of a functional intestine, are extremely rare during normal homeostasis and are characterized by enhanced expression of the pro-survival stress response gene Clu (Zhang et al., 2014). Interestingly, we find that many of the genes that make up the molecular signature of these quiescent revSCs are also enriched in qCSCs and have been found to regulate therapy resistance in various types of cancer. These common genes include CLU (Koltai, 2014); CTSD (Mahajan et al., 2020; Oliveira et al., 2015); CDKN1A (Koster et al., 2010; Liu et al., 2013; Maiuthed et al., 2018; Morris-Hanon et al., 2017; Xia et al., 2011); EMP1 (Ariës et al., 2014; Jain et al., 2005); MUC3 (Lesuffleur et al., 1993); LAMC2 (Huang et al., 2017); KRT19 (Asfaha et al., 2015); LGALS3 (Wang et al., 2019); F3, ITM2B, and ITGB4 (Folgiero et al., 2008; Stewart and O'Connor, 2015); CDH17 (Atukorala and Mathivanan, 2018; Qiu et al., 2013); and GSN (Chung et al., 2015; Ilmer et al., 2016) (Figures 3D, 3E, and S1 and Data S1). Considering that colon cancer is a heterogeneous tumor that recapitulates the cellular hierarchy of the intestine, these data suggest that the qCSCs identified here may be the tumor equivalent of revSCs. However, in contrast to revSCs and previous studies on qCSCs, our data demonstrate that qCSCs are enriched for p53 signaling (Figure 3, Figure 4A and 4A).

Figure 4

qCSCs are enriched for p53-interacting negative regulators of cell cycle that are indicators of poor prognosis and may be targeted for qCSC elimination

(A) RNA sequencing-generated normalized counts for negative cell cycle regulator and p53 target genes AKAP12, CD82, CDKN1A, FHL2, GPX3, KIAA0247, LCN2, TFF2, UNC5B, and ZMAT3 in PKH26Negative and PKH26Positive cells (±SD; n = 4). ns = not significant; ∗p < 0.05; ∗∗p < 0.01; ∗∗∗p < 0.001; ∗∗∗∗p < 0.0001 (t test). See also Data S1.

(B) Kaplan-Meier survival curves for AKAP12, CD82, CDKN1A, FHL2, GPX3, KIAA0247, LCN2, TFF2, UNC5B, and ZMAT3 in patients with colorectal cancer comparing lower quartile to upper quartile (log rank p values = 2.2 × 10−6, 0.004, 8.1 × 10−5, 0.00012, 0.0003, 0.049, 0.00018, 5.7 × 10−5 and 0.3, respectively). Of these, higher AKAP12, CD82, CDKN1A, FHL2, GPX3, KIAA0247, LCN2, TFF2, and UNC5B are significant at FDR <10%. Results based on data generated by the Kaplan-Meier Plotter (kmplot.com) (Nagy et al., 2018).

(C) Representative FACS plot of PKH26-labeled 151-ML-M control-GFP cells (top row) and shRNA CDKN1A-GFP cells (bottom row) after 12 h and 12 days.

(D) Frequency of PKH26Positive LRCs in shRNA CDKN1A PDO models after 12 days compared with control virus-transduced cells (±SD; n = 3) (see also Figure S2).

Negative cell cycle regulators downstream of p53 are indicators of poor prognosis and required for qCSC maintenance

Loss of p53 signaling in hematopoietic stem cells (HSCs) and neural stem cells (NSCs) causes these cells to exit quiescence and enter the cell cycle (Itahana et al., 2002; Liu et al., 2009; Meletis et al., 2006). Targeting p53 signaling may have the same effect in qCSCs but is complicated by the role of p53 as a tumor suppressor and guardian of homeostasis (Lane, 1992). However, targeting negative cell cycle regulators downstream of p53 may provide novel strategies for qCSC elimination without affecting the role of p53 in healthy cells. In addition, such strategies may also enable cell cycle activation of qCSCs in p53 mutant tumors. Differential gene expression analysis, comparing PKH26Negative and PKH26Positive cells, identified the negative cell cycle regulators AKAP12 (Gelman, 2010; Lin et al., 2000; Liu et al., 2011; Reggi and Diviani, 2017), CD82 (Hur et al., 2016), CDKN1A (El-Deiry, 1993; Vogelstein et al., 2000; Wade Harper, 1993; Xiong et al., 1993), FHL2 (Hellerbrand, 2010; Labalette et al., 2008; Lee et al., 2006; Martin et al., 2007), GPX3 (An et al., 2018; Barrett et al., 2013; Wang et al., 2012), KIAA0247 (Huang et al., 2011; Polato et al., 2014), LCN2 (Chakraborty et al., 2012; Chiang et al., 2016; Kim et al., 2017), TFF2 (Bossenmeyer-Pourié et al., 2002; Dubeykovskaya et al., 2016, 2019; Thim, 1997; Tu et al., 2009), UNC5B (Huang et al., 2020; Kong et al., 2016; Okazaki, 2011), and ZMAT3 (Bersani et al., 2014; Hellborg et al., 2001) to be enriched in qCSCs (Figure 4A and Data S1). Significantly, each of these genes is a target of p53 (Bersani et al., 2014; Fischer, 2017; Gelman, 2010; Lane, 1992; Lee et al., 2006; Liu et al., 2011; Marreiros et al., 2005; Miyamoto et al., 2016; Polato et al., 2014; Rouillard et al., 2016; Soutto et al., 2014), and with the exceptions of LCN2 and ZMAT3, is associated with reduced survival in colorectal cancer (CRC) (Figure 4B). Interestingly, CD82, KIAA0247, and UNC5B proteins localize to the cell surface and may therefore have potential as new markers for the prospective isolation of qCSCs in CRC. Indeed, CD82 has previously been identified as a marker for prospectively isolating stem cells from human fetal and adult skeletal muscle and is a functional surface marker of long-term HSCs (Alexander et al., 2016; Hur et al., 2016).

Deletion of CDKN1A (p21), which is the downstream mediator of p53 induced cell-cycle arrest (El-Deiry, 1993; Georgakilas et al., 2017), leads to cell cycle activation and exhaustion of quiescent HSCs and NSCs (Cheng, 2000; Kippin et al., 2005). In addition, CDKN1A is highly expressed in noncycling intestinal crypt base stem cells (Powell et al., 2012) and revSCs (Ayyaz et al., 2019). We therefore selected CDKN1A as a candidate gene to determine whether targeting the p53 pathway would eliminate qCSCs. Significantly, short hairpin RNA (shRNA)-mediated knockdown of CDKN1A (Figure S2) in PKH26-labeled qCSCs resulted in the elimination of PKH26Positive label-retaining qCSCs (Figures 4C and 4D).

Discussion

Colon cancer is a heterogeneous tumor entity containing a subpopulation of qCSCs that may promote tumor cell heterogeneity, plasticity, and resistance to various types of stress, including resistance to conventional treatments (Moore and Lyle, 2011). However, the rarity and plasticity of qCSCs has made them an elusive and challenging cell state to define and target. Here, we provide the first whole-transcriptome analyses of a population of colon cancer patient-derived long-term label-retaining qCSCs and identify genes that may provide novel targets for their elimination.

Label retention has previously been used as a strategy for the isolation of both healthy quiescent tissue stem cells and qCSCs from a variety of cancer types (Buczacki et al., 2013; Cotsarelis et al., 1990; Dembinski and Krauss, 2009; Moore and Lyle, 2011; Pece et al., 2010; Roesch et al., 2010; Smith, 2005; Wilson et al., 2008). In agreement with these studies, we show that PKH26Positive LRCs isolated from colon cancer PDOs are qCSCs capable of entering the cell cycle and self-renewing after replating in adherent and non-adherent cell culture conditions and maintain long-term quiescence in xenograft models. Interestingly, in vivo these cells were located at the tumor border, suggesting that quiescence may be induced at the invasive tumor front where such cells may be primed for metastatic dissemination. This is in agreement with previous studies showing cell-cycle arrest/decreased proliferation and increased levels of Wnt signaling at the invasive front of colorectal tumors (Dawson et al., 2014; Harbaum et al., 2011; Jung et al., 2001; Rubio, 2008; De Smedt et al., 2017).

RNA sequencing of qCSCs demonstrated that they display the molecular hallmarks of quiescence (Cho et al., 2019) while also being enriched for the same developmental and stem cell-associated gene sets previously described for actively cycling ALDHPositive CSCs (Regan et al., 2017), which unlike PKH26Positive LRCs are enriched in PDOs. We previously reported that hedgehog signaling in active colon CSCs is non-canonical (SHH-dependent, PTCH-dependent, SMO-independent, GLI-independent) and acts as a positive regulator of Wnt signaling for CSC survival (Regan et al., 2017). In agreement with our work, a subsequent study from (Buczacki et al., 2018) demonstrated that qCSC survival in CRC is also dependent on non-canonical hedgehog signaling-mediated regulation of Wnt signaling (Buczacki et al., 2018). In addition, several of the genes common to both the revSCs reported by Ayyaz et al. (2019) and qCSCs, namely, CLU (Schepeler et al., 2007), CTSD (Basu et al., 2019), CDKN1A (Xu et al., 2016), EMP1 (Yao et al., 2011), MUC3 (Pai et al., 2016), LAMC2 (Sánchez-Tilló et al., 2011), KRT19 (Saha et al., 2017), LGALS3 (Korkmaz et al., 2016), F3 (Camps et al., 2020; Kinchen et al., 2018), ITGB4 (Avvisato et al., 2007), CDH17 (Wang et al., 2013), and GSN (Shimura et al., 2004), are targets and/or regulators of Wnt signaling. Overall, these data demonstrate that both cycling and non-cycling CSCs share overlapping molecular profiles and further support the targeting of non-canonical hedgehog signaling to prevent disease relapse (Buczacki et al., 2018; Regan, 2018; Regan et al., 2017).

However, the molecular mechanisms that distinguish non-cycling qCSCs from cycling CSCs required further elucidation. p53 signaling plays a crucial role in regulating cellular stress responses such as DNA damage repair, senescence, apoptosis, and cell-cycle arrest in virtually all cell types (Vogelstein et al., 2000; Vousden and Lane, 2007). In addition, it is an important regulator of stem cell self-renewal and differentiation in embryonic and adult tissue stem cells (Jain and Barton, 2018; Meletis et al., 2006; Tosoni et al., 2015) and CSCs (Freed-Pastor et al., 2012; Tschaharganeh et al., 2014; Zhao et al., 2010). Significantly, it has also been demonstrated to be essential for the maintenance of quiescence in HSCs, NSCs, muscle stem cells, and lung progenitor cells (Flamini et al., 2018; Itahana et al., 2002; Liu et al., 2009; McConnell et al., 2016; Zheng et al., 2008).

Here we show that qCSCs, in contrast to cycling ALDHPositive CSCs (Regan et al., 2017), are enriched for p53 signaling genes (Figures 3A and 3D). p53 is mutated in 40%– 50% of CRCs. Reflecting this, half the tumors included in our study contain a p53 mutation (Table S1). However, regardless of mutation status, all PDO models expressed CDKN1A, the downstream mediator of p53-induced cell-cycle arrest (El-Deiry, 1993; Georgakilas et al., 2017), which we demonstrate is enriched in and required for the maintenance of qCSCs. This is in agreement with other studies showing CDKN1A expression to be independent of p53 mutation status (van Oijen et al., 1998; Perez et al., 2016; Phalke et al., 2012; Sato et al., 1997; Zhang et al., 2015). For example, in head and neck squamous cell carcinoma CDKN1A is also overexpressed in non-proliferating cells in both p53 wild-type and p53 mutant tumors (van Oijen et al., 1998). The mechanisms regulating the enrichment of p53 signaling genes in p53 mutant qCSCs are unclear but undoubtedly reflect the complexity of the p53 pathway, in which cellular context (cell state, epigenetic status, microenvironment, activating signal) may be the predominant determinant of signaling activity and biological outcome (Kastenhuber and Lowe, 2017). These data suggest that targeting negative regulators of cell cycle downstream of p53, in both wild-type and mutant p53 tumors, may provide novel therapeutic “lock-out” strategies to induce the proliferation of qCSCs and thereby sensitize them to chemotherapeutics and prevent relapse (Cho et al., 2019; Kobayashi et al., 2011; Takeishi et al., 2013).

Considering the role of p53 as a tumor suppressor and guardian of homeostasis in healthy tissues, as well as its inactivation in many cancers, most strategies to date have focused on the development of p53 activators (Levine, 2019). However, our data, and others, suggest that strategies that activate p53 may lead to therapy resistance (Kim et al., 2009). For example, in glioma there is evidence that p53 DNA repair activities contribute to overall survival potential and drug resistance (Batista et al., 2007) and in breast cancer p53 induces senescence, drives resistance to therapy, and is associated with poor therapeutic response and overall survival (Jackson et al., 2012; Ungerleider et al., 2018).

Inhibiting p53 could interfere with its role in normal tissue homeostasis or lead to the activation of senescent cancer cells in other tissues. However, healthy cells have lower p53 expression levels than cancer cells (Rogel et al., 1985) and single-dose treatments, that avoid the unwanted consequences of sustained p53 inhibition, may be sufficient to eliminate qCSCs. This was recently demonstrated by Webster et al. (2020) in melanoma, where a single dose of p53 inhibitor during the early stage of BRAF/MEK inhibitor treatment resulted in improved response to therapy (Webster et al., 2020).

In addition, targeting negative cell cycle regulators downstream of p53, such as those identified here (AKAP12, CD82, CDKN1A, FHL2, GPX3, KIAA0247, LCN2, TFF2, UNC5B, and ZMAT3), may provide novel strategies for activating cell cycle in qCSCs without affecting the role of p53 in healthy cells. For example, p53-dependent activation of p21 (CDKN1A) is an important axis in senescence-dependent tumor suppression. However, despite p21 playing an important role in mediating the p53-dependent cellular response to stress, lack of p21 does not promote tumor development (Choudhury et al., 2007). Furthermore, p21 maintains CSC self-renewal, limits proliferation, and confers therapy resistance in numerous cancer types in which its temporal inhibition has been proposed as a strategy to overcome resistance to DNA-damaging chemotherapy and radiation (El-Deiry, 2016; Tian et al., 2000; Kokunai et al., 2001; Koster et al., 2010; Maiuthed et al., 2018; Morris-Hanon et al., 2017; Viale et al., 2009; Weiss, 2003; Xia et al., 2011). Indeed, several small molecule inhibitors of p21 have been reported, including butyrolactone I (Joanna et al., 2002), LLW10 (Park et al., 2008), sorafenib (Inoue et al., 2011), and UC2288 (Wettersten et al., 2013), that could serve as novel drugs for the elimination of therapy-resistant qCSCs. Significantly, these p21 inhibitors have been shown to be p53 independent and so may be effective in both wild-type and p53 mutant cancers without perturbing p53 function in healthy cells (Liu et al., 2013).

These data demonstrate the existence of long-term p53 signaling-dependent qCSCs in colon cancer and provide evidence supporting the temporal inhibition of p53 signaling, in combination with standard-of-care treatments, for the elimination of qCSCs and prevention of disease relapse. The p53 target genes identified here, along with the publication of our qCSC whole-transcriptome data, will provide a valuable resource for the development of such therapeutic strategies in the future.

Limitations of the study

Our data suggest that qCSCs may be the tumor equivalent of revSCs, which were recently identified by scRNA-seq of the normal mouse intestine during homeostasis and injury (Ayyaz et al., 2019). However, revSCs have not yet been identified in the human intestine and their potential role in carcinogenesis and tumor heterogeneity remain to be determined. In addition, xenotransplantation of PKH26Positive cells demonstrated that they maintain quiescence in vivo for up to 80 days. Unfortunately, due to the rarity of these cells in vivo, we were unable to isolate sufficient numbers for RNA sequencing. Future scRNA-seq studies can be used to overcome the obstacle of qCSC rarity in vivo and will be required to determine lineage relationships, e.g., between revSCs and qCSCs, and cell-state heterogeneity within colon cancer (Giraddi et al., 2018; Wahl and Spike, 2017). Furthermore, the molecular mechanisms that maintain high levels of p53 signaling genes in p53 mutant qCSCs remain unclear, but are likely influenced by cellular context, including cell state and epigenetic factors. The elucidation of such mechanisms may be facilitated through future single-nucleus transposase-accessible chromatin sequencing (snATAC-seq), which examines global chromatin accessibility and has been shown to be a better indicator of cell state than scRNA-seq (Chung et al., 2019; Regan and Smalley, 2020). The transcriptome data of functionally tested qCSCs presented here will help to inform cell state designations in such future large-scale single-cell profiling endeavors.

STAR★methods

Key resources table

Resource availability

Lead contact

Further information and requests for resources should be directed to and will be fulfilled by the lead contact Joseph L. Regan (joseph.regan@charite.de).

Materials availability

This study did not generate new unique reagents.

Data and code availability

The accession number for the RNA-sequencing data reported in this paper is ArrayExpress: E-MTAB-8924.

Experimental model and subject details

Human tissue samples

Tumor material was obtained with informed consent from colon cancer patients (3 male and 3 female; detailed in Table S1) under approval from the local Institutional Review Board of Charité University Medicine (Charité Ethics Cie: Charitéplatz 1, 10117 Berlin, Germany) (EA 1/069/11) and the ethics committee of the Medical University of Graz and the ethics committee of the St John of God Hospital Graz (23-015 ex 10/11). Tumor staging was carried out by experienced and board-certified pathologists (Table S1). Age of subjects can be inquired through info@oncotrack.eu.

Patient-derived organoids

Cancer organoid cultures were established and propagated as previously described(Sato et al., 2011; Schütte et al., 2017). Briefly, resected tumor samples were enzymatically digested with Collagenase IV (C9407, Sigma-Aldrich), DNaseI (A3778,0050, AppliChem) and Dispase (07913, Stem Cell Technologies) at 37°C for 60 min. Suspensions were washed, filtered, and depleted of red blood cells using Red Blood Cell Lysis Solution (00-4333-57, Invitrogen). Cells were mixed with phenol-red free growth factor-reduced Matrigel (356231, Corning) and seeded into 24-well plates. Solidified droplets were overlaid with culture medium consisting of Advanced DMEM/F12 (12634-010, Gibco) supplemented with 1% penicillin/streptomycin, 1% HEPES buffer (1064859, Fisher Scientific), 1% Glutamax, 1x N2 (#17502-048, Invitrogen), 1x B27 (17504-044, Invitrogen), 50 ng/ml EGF (E9644, Sigma), and 1mM N-acetylcysteine (A9165-5G, Sigma) and maintained at 37°C. Organoids were released from Matrigel and passaged by adding 5 ml Advanced DMEM/F12 followed by centrifugation and digestion of pellets with TrypLE Express (12604-013, Gibco).

Xenotransplantation

Housing and handling of animals followed European and German Guidelines for Laboratory Animal Welfare. Animal experiments were conducted in accordance with animal welfare law, approved by local authorities, and in accordance with the ethical guidelines of Bayer AG. PDOs were processed to single cells and labelled with PKH26 (PKH26GL, Sigma-Aldrich) following manufacturer’s instructions and DAPI (to exclude dead cells). PKH26Positive DAPINegative cells were collected by FACS and immediately transplanted by injecting subcutaneously in PBS and Matrigel (1:1 ratio) into female 8 – 10-week-old nude-/- mice at 1000 cells per animal. The purity of the sorted PKH26Positive cell population was confirmed by post-sort FACS analysis.

Method details

Cell cycle analysis and colony forming assays

Cell cycle analysis was carried using the Click-iT EdU assay (Invitrogen, #C10337) and assessed by FACS on a BD LSR II analyzer. For colony forming assays, PDOs were processed to single cells and labelled with CellTracker™ CM-DiI fluorescent dye (C7000, Thermo Fisher) or PKH26 following manufacturer’s instructions and DAPI (to exclude dead cells). PKH26Positive DAPINegative (live) cells were sorted by FACS (BD FACS Aria II) into Matrigel culture. After 12 days, PDOs were once again processed to single cells and sorted by FACS, seeding PKH26Positive DAPINegative cells and PKH26Negative DAPINegative cells separately at limiting dilution into 96-well adherent Matrigel and 384-well non-adherent ultra-low attachment plates at a frequency of 100 and 1 cell per well, respectively. The purity of the sorted PKH26Positive cell population was confirmed by post-sort FACS analysis. PDO sizes were determined by ImageJ software analysis. Ultra-low attachment wells containing spheroids were counted and used to calculate the CSC frequency using ELDA software (http://bioinf.wehi.edu.au/software/elda/index.html)(Hu and Smith, 2009).

Immunofluorescence staining

Tumors were fixed in 4% paraformaldehyde overnight and cryopreserved in OCT compound. Immunohistochemistry of frozen sections was carried out via standard techniques with α-Tubulin (T5168, mouse monoclonal, Sigma; diluted 1:1000) and a secondary conjugated antibody at room temperature for 2 hours. For immunofluorescence imaging of PDOs, cultures were fixed in 4% paraformaldehyde for 30 min at room temperature and permeabilized with 0.1% Triton X-100 for 30 min and blocked in phosphate-buffered saline (PBS) with 10% bovine serum albumin (BSA). F-actin was stained with Alexa Fluor® 647 Phalloidin (#A22287, Thermo Fisher; diluted 1:20) for 30 min at room temperature. Nuclei were counterstained with DAPI. Negative controls were performed using the same protocol with substitution of the primary antibody with IgG-matched controls. Cancer organoids were then transferred to microscope slides for examination using a Zeiss LSM 700 Laser Scanning Microscope.

RNA sequencing

PKH26Negative DAPINegative and PKH26Positive DAPINegative cells from four different biological replicates were isolated by FACS. Cells were lysed in RLT buffer and processed for RNA using the RNeasy Mini Plus RNA extraction kit (Qiagen). Samples were processed using NuGEN’s Ovation RNA-Seq System V2 and Ultralow V2 Library System and sequenced on an Illumina HiSeq 2500 machine as 2x125nt paired-end reads. The raw data in Fastq format were checked for sample quality using our internal NGS QC pipeline. Reads were mapped to the human reference genome (assembly hg19) using the STAR aligner (version 2.4.2a) (Dobin et al., 2013). Total read counts per gene were computed using the program featureCounts (version 1.4.6-p2) in the subread package, with the gene annotation taken from Gencode (version 19). Variance-stabilizing transformation from the Bioconductor package DESeq2 (Love et al., 2014) was used for normalization and differential-expression analysis.

Viral transduction

Cells were seeded in 100 μl volumes of antibiotic free culture media at 1.0 x105 cells per well in ultra-low attachment 96-well plates. Control and shRNA lentiviruses were purchased from Sigma-Aldrich (Table S2). Viral particles were added at a multiplicity of infection of 1. Cells were transduced for up to 96 h or until GFP positive cells were observed before being embedded in Matrigel for the establishment of lentiviral transduced cancer organoid cultures. Puromycin (2 μg/ml) was used to keep the cells under selection.

Gene expression analysis

For quantitative real-time RT-PCR analysis RNA was isolated using the RNeasy Mini Plus RNA extraction kit (Qiagen). cDNA synthesis was carried out using a Sensiscript RT kit (Qiagen). RNA was transcribed into cDNA using an oligo dTn primer (Promega) per reaction. Gene expression analysis was performed using TaqMan® Gene Expression Assays (Applied Biosystems) (Table S3) on an ABI Prism 7900HT sequence detection system (Applied Biosystems). GAPDH was used as an endogenous control and results were calculated using the Δ-ΔCt method. Data were expressed as the mean fold gene expression difference in three independently isolated cell preparations over a comparator sample with 95% confidence intervals. Survival curves were generated using the Kaplan-Meier Plotter (www.kmplot.com/analysis) (Nagy et al., 2018). Gene ontology enrichment analysis was carried out using the Gene Ontology Resource (www.geneontology.org)(Ashburner et al., 2000; The Gene Ontology Consortium, 2018).

Quantification and statistical analysis

GraphPad Prism 8.0 was used for data analysis and imaging. All data are presented as the mean ± SD, followed by determining significant differences using the two-tailed t test. Significance of RT-PCR data was determined by inspection of error bars as described by Cumming et al. (2007) (Cumming et al., 2007). Gene set enrichment analysis was carried out using pre-ranked feature of the Broad Institute GSEA software version 2 using MSigDB v5.1 gene sets (Liberzon et al., 2015; Subramanian et al., 2005). The ranking list was derived from the fold changes calculated from the differential gene expression calculation and nominal p-values. P-values <0.05 were considered as statistically significant. The representation factor and the associated probability of finding an overlap were calculated using http://nemates.org/MA/progs/representation.stats.html. Survival curves were generated using the Kaplan-Meier Plotter (www.kmplot.com/analysis)(Nagy et al., 2018). For the final list of significant genes, False Discovery Rate (FDR) was computed using the Benjamini-Hochberg method (Benjamini and Hochberg, 1995).

REAGENT OR RESOURCE	SOURCE	IDENTIFIER
Antibodies
Mouse monoclonal anti-α-Tubulin	Sigma	T5168; RRID: AB_477579
Alexa Fluor® 647 Phalloidin	Thermo Fisher	#A22287; RRID: AB_2620155
Biological samples
Patient-derived organoids (PDOs)	Charité Universitätsmedizin Berlin, Germany	OncoTrack
Chemicals, peptides, and recombinant proteins
CellTracker™ CM-DiI fluorescent dye	Thermo Fisher	C7000
PKH26 Red Fluorescent Cell Linker	Sigma-Aldrich	PKH26GL
Critical commercial assays
Click-iT EdU assay	Invitrogen	#C10337
RNeasy Mini Plus RNA extraction kit	Qiagen	ID:74136
Sensiscript RT kit	Qiagen	ID:205213
Ovation RNA-Seq System V2	NuGEN	7102-A01
Ultralow V2 Library System	NuGEN	0344NB-A01
Deposited data
Array data	This paper; www.ebi.ac.uk/arrayexpress	E-MTAB-8924
Oligonucleotides
CDKN1A MISSION shRNA Lentiviral Transduction Particles: TRCN0000040123, TRCN0000287021, TRCN0000287091	Sigma-Aldrich	SHCLNV-NM_000389
CDKN1A (Hs00355782_m1)	Thermo Fisher	#4331182
GAPDH (Hs02758991_g1)	Thermo Fisher	#4331182
Software and algorithms
ELDA software	Hu and Smith, 2009	http://bioinf.wehi.edu.au/software/elda/index.html
ImageJ	Schneider et al., 2012	https://imagej.nih.gov/ij/
DESeq2	Love et al., 2014	https://bioconductor.org/packages/release/bioc/html/DESeq2.html
STAR aligner (version 2.4.2a).	Dobin et al., 2013	https://github.com/alexdobin/STAR
Kaplan-Meier Plotter	Nagy et al., 2018	www.kmplot.com/analysis
Gene Ontology Resource	Ashburner et al., 2000	(www.geneontology.org)
GSEA software	Liberzon et al., 2015; Subramanian et al., 2005	https://www.gsea-msigdb.org/gsea/index.jsp
Prism 8.0	GraphPad	https://www.graphpad.com/scientific-software/prism/