Selective advantage of epigenetically disrupted cancer cells via phenotypic inertia

The evolution of established cancers is driven by selection of cells with enhanced ﬁtness. Subclonal mutations in numerous epigenetic regulator genes are common across cancer types, yet their functional impact has been unclear. Here, we show that disruption of the epigenetic regulatory network increases the tolerance of cancer cells to unfavorable environments experienced within growing tumors by promoting the emergence of stress-resistant subpopulations. Disruption of epigenetic control does not promote selection of genetically deﬁned subclones or favor a phenotypic switch in response to environmental changes. Instead, it prevents cells from mounting an efﬁcient stress response via modulation of global transcriptional activity. This ‘‘transcriptional numbness’’ lowers the probability of cell death at early stages, increasing the chance of long-term adaptation at the population level. Our ﬁndings provide a mechanistic explanation for the widespread selection of subclonal epigenetic-related mutations in cancer and uncover phenotypic inertia as a cellular trait that drives subclone expansion. trypsin Addgene #86677) (Raab et al., 2016), or mCherry-NLS (pTRIP-SFFV-mCherry-NLS). To generate pTRIP-SFFV-mCherry-NLS, mCherry was ampliﬁed from pUAS-mCherry-NLS (Addgene #87695) (Zhang et al., 2017) by PCR (forward primer: TAGCGGATCCATGGTGAGCAAGGGCGAG, reverse primer: AGGTCTCGAGCGCCCCGACTCTAGATTATAC) and cloned in the pTRIP-SFFV expression vector after GFP removal by BamHI-XhoI digestion. Seven days after infection, cells expressing similar levels of the ﬂuorescent proteins were isolated by ﬂow cytometry. mCherry-labelled lines were transduced with sgRNA constructs targeting epigenetic regulators, while GFP-labelled cells with sgRNAs targeting the non-expressed gene TNP2 . After 10 days of blasticidin selection and induction of KO through doxycycline treatment, mCherry-labelled KO populations were mixed with GFP-labelled control cells in equal quantities and seeded in 96-well Twist Human Core Exome kit to the manufacturer’s instructions (#100578; #101174; #100254). 4000 using paired end 100 bp Transcriptomic analysis was performed in 6 different KO populations of MEXF 2090 cells in which EED , EZH2 , HIST1H1B , SMARCD1 , SUZ12 and the non-expressed gene TNP2 were targeted and that were grown under nutrient deprivation for 12 days. Four and six replicates of unperturbed and nutrient-deprived cells, respectively, grown in distinct wells of 96-well plates were pooled and total RNA extraction was performed using a RNeasy Plus Micro kit. Total RNA was quantiﬁed using an RNA QuantiFluor RNA system (Promega, #E3310) on the GloMax Multi Detection system following the manufacturer’s guidelines. RNA quality was assessed via the High Sensitivity RNA ScreenTape using the TapeStation 4200 (Agilent Technologies, #5067-5579). RNA was normalized to 30 ng and used for cDNA synthesis and library preparation using the QuantSeq 3 0 mRNA-Seq FWD kit (Lexogen, #015) according to the manufacturer’s instructions. Libraries were sequenced on the Illumina HiSeq 4000 with single ended reads. of cells 0.67 conﬁrming experimental technical developed R scripts. Enrichment analysis of each cell population in each cell cluster was performed by hypergeometric test: for each cell type, p values were calculated from the cumulative distribution function, using the phyper function for R and corrected for multiple testing using Benjamini-Hochberg method. Only clusters signiﬁcantly enriched for at least one a cell population (adjusted p value < 0.05) were used in the following steps of the analysis. Gene signatures deﬁning each KO- or control-enriched clusters consisted of DE-Gs identiﬁed by MAST (Finak et al., 2015) between each cell cluster and all other cells of the pair merged set. Genes with FDR <0.05 and a non-zero log 2 -fold-change (with a non-zero-crossing 95% conﬁdence intervals) were considered differentially expressed and used to characterize the cells belonging to the cluster. After performing the clustering analysis on cells for all KO-control comparisons, all the identiﬁed cell-cluster-speciﬁc gene signatures were clustered together to identify DE-Gs shared by different cell-clusters. To increase the robustness of the results, signatures consisting of less than 50 DE-Gs or that corresponded to cell clusters with a size smaller than 5% of the pooled set size were removed. The gene signatures were deﬁned as vectors where each element corresponded to the product of the point estimate of the log2-fold-change and the signiﬁcance (equal to (cid:1) log 10 p ) of each DEG, as determined by MAST . All non DEGs were assigned a signature value equal to 0. Gene signatures were clustered using a hierarchical clustering. Spearman’s correlation was used as similarity measure and the optimal linkage corresponded to the maximum cophenetic correlation coefﬁcient (Sokal and Rohlf, 1962). The optimal number of clusters corresponded to the maximum silhouette coefﬁcient (Rousseeuw, 1987) calculated with a varying number of clusters in f 2 ; 3 ; 4 ; . ; N sig g , where N sig is equal to the total number of gene signatures. Genes sets shared by KO-enriched or control-enriched signatures were denoted as meta-signatures to distinguish them from the cell cluster-speciﬁc for H3K27me3 and the –narrow_peak option for H3K27Ac, H3K4me3 and H4K4me1. This produced the BigWig ﬁles used for subsequent analysis. Metaproﬁles at ﬁtness and stress gene promoters were generated using deepTools: the computeMatrix and plotProﬁles functions were used against the BigWig ﬁles for each sample to provide a matrix of the density of ChIP signal and generate metaproﬁles in a 5Kb window centered on gene TSSs. For quantiﬁcation of differential peaks, the DESeq2 output of the nf-core pipeline was further analyzed to eliminate peaks of unlikely biological signiﬁcance: detected peaks were ﬁrst ﬁltered based on magnitude, retaining regions with pseudo-counts R 250 in at least one sample (maximal value: 25,035). The fraction of differential peaks (FDR <0.01, log2 FC R 1 or % 1) was then calculated, assessing either all peaks genome-wide or peaks at promoters of ﬁtness and stress genes. Promoter regions were identiﬁed using the nf-core pipeline default annotation method. Fitness and stress genes were deﬁned as stated in the ATAC-seq section.

Correspondence paola.scaffidi@crick.ac.uk In brief Loukas et al. show that loss of many epigenetic regulators, frequently mutated across cancer types, promote subclone expansion by conferring resistance to stressful environments. Epigenetically disrupted cells fail to efficiently rewire transcription in response to stress, lowering the probability of cell death at early stages and favoring stochastic population adaptation.

INTRODUCTION
Cancers comprise heterogeneous populations of malignant cells characterized by distinct biological properties, which differentially contribute to disease maintenance and progression. The phenotype and properties of each cell are shaped by a combination of genetic, epigenetic, and environmental factors (Marusyk et al., 2012). Genetic intratumor heterogeneity is prevalent across cancers, and many tumors contain spatially separated subclonal populations that evolve following distinct trajectories (Dentro et al., 2021;McGranahan and Swanton, 2017). Evidence of positive selection of fitter subclones exists, with dominant subclonal lineages detected in various cancer types (Vendramin et al., 2021). However, the mechanisms that drive sublclone expansion remain only partially understood. At the pan-cancer level, sublconal alterations in canonical cancer drivers that promote cancer initiation (e.g., KRAS, TP53, CDKN2A) are present in some patients, but they are relatively rare (Dentro et al., 2021), raising the possibility that distinct mechanisms may drive tumorigenesis and subsequent cancer evolution. Furthermore, only 11% of detected subclones harbor mutations in functionally validated cancer genes, and it has been difficult to identify recurrent, subclonally mutated genes (Dentro et al., 2021). Instead, classes of genes, such as those involved in the DNA damage response, epigenetic regulators, and splicing factors, are enriched among subclonally mutated genes (Dentro et al., 2021;Jamal-Hanjani et al., 2017). These observations suggest that phenotypic traits induced by disruption of the corresponding biological processes (i.e., genome maintenance and gene expression regulation), rather than specific genotypes, may be selected during cancer evolution.
Epigenetic regulator genes encode proteins involved in the establishment and maintenance of chromatin and DNA methylation landscapes and, as a group, are among the most mutated genes across cancer types (Shen and Laird, 2013;Lawrence et al., 2014). Genetic alterations include translocations that create oncogenic fusion proteins (e.g., MLL fusions in leukemia), hot-spot gain-of-function (GOF) mutations (e.g., EZH2 Y646 substitutions in B cell lymphoma), and loss-of-function (LOF) alterations that affect a wide range of regulators (Shen and Laird, 2013). In contrast to oncogenic alterations, which are highly tissue specific, LOF mutations are often observed across multiple cancer types, with varying recurrence (Shen and Laird, 2013;Lawrence et al., 2014;Mittal and Roberts, 2020). Moreover, while translocations and GOF mutations are typically clonal, in line with an early role in driving carcinogenesis (Bö dö r  Nikbakht et al., 2016), LOF mutations can be clonal (e.g., PBRM1 in clear cell renal cell carcinoma) (Turajlic et al., 2018), be mainly present in emerging subclones (e.g., SETD2 in various cancers) (Dentro et al., 2021), or exhibit complex patterns (e.g., ARID1A), with early inactivation in some contexts (Wiegand et al., 2010;Martincorena et al., 2018) and later, subclonal mutations in others (Dentro et al., 2021). These patterns suggest that loss of epigenetic regulators may play distinct role in tumorigenesis and cancer evolution. Mechanisms underpinning early oncogenic functions have been characterized for some genes, and they often promote transformation by interfering with cellular differentiation (Shen and Laird, 2013). On the other hand, how epigenetic deregulation induced by late mutations further enhances cancer cell fitness and favors subclone expansion within established tumors is poorly understood. Cancer cell fitness, defined as the number of descendants generated by a malignant cell in a given time frame, is an integrative phenotype that combines long-term proliferative potential and the ability to withstand both cell-intrinsic stress (e.g., genotoxic, proteotoxic, and metabolic stress) and environmental challenges faced during tumor growth. Regional differences in selective pressures such as hypoxia, acidity, limited availability of nutrients, and physical constraints shape the tumor composition by favoring expansion of subpopulations of cells that can withstand those challenges (Junttila and de Sauvage, 2013). The fitness of cancer cells is thus highly dependent on the environment in which they grow and the stress they experience. In this study, we modeled the widespread epigenetic-related mutations observed in patients and examined their impact on cancer cell fitness under stressful conditions typical of the tumor microenvironment.

RESULTS
Selection of a disrupted epigenetic regulatory network during cancer evolution Recurrently mutated epigenetic regulators identified through unbiased analysis of patient samples belong to diverse functional classes, including chromatin remodelers, histones, histone modifiers, regulators of genome topology, and DNA modifiers (Lawrence et al., 2014). To comprehensively examine the distribution of mutations across the epigenetic regulatory network independently of recurrence at the individual-gene level, we probed the status of more than 300 genes (Table S1) in six pan-cancer datasets available through the cBioPortal (Zehir et al., 2017;Robinson et al., 2017;Miao et al., 2018;Hyman et al., 2018; ICGC/ TCGA Pan-Cancer Analysis of Whole Genomes Consortium, 2020;Nguyen et al., 2022;Gao et al., 2013). Overall we identified 60,907 nonsynonymous mutations spread across the whole network, with no functional class being particularly affected (Figures 1A and 1B). As observed in targeted studies, approxi-mately one-third of these alterations across all classes were truncating LOF mutations, a fraction comparable with that of established tumor-suppressor genes ( Figure 1C). At the pan-cancer level, LOF mutations were also prevalent in genes that act as oncogenic drivers in specific cancer types (e.g., EZH2 in B cell lymphoma; KMT2A, CREBBP, and EP300 in leukemia) or that have been linked to context-dependent synthetic lethalities (e.g., EZH2 in ARID1A-deficient cancers [Bitler et al., 2015] and MEN1 in MLL-driven leukemia [Yokoyama et al., 2005]) (Figures S1A-S1C).
While alterations with high variant allele frequency (VAF) were detected (maximal VAF across patients R0.5 for 258 genes, Figure 1A), mutations across the epigenetic regulatory network displayed overall low VAF (median VAF <0.3 for 271 genes, Figure 1A), raising the possibility that many tumors may only harbor these alterations in subsets of cells. To directly examine subclonal mutations selected at later stages of tumor evolution, we analyzed 100 non-small cell lung cancer (NSCLC) patients profiled by multiregion sequencing in the TRACERx study (Jamal-Hanjani et al., 2017). Subclonal nonsynonymous mutations in epigenetic regulator genes were prevalent, with 62% of the samples harboring mutations in at least one gene and 130 genes being affected in at least one patient ( Figure 1D). These alterations showed a high nonsynonymous/synonymous ratio at the individual-gene level, which suggests their positive selection during tumor growth (Martincorena et al., 2017) ( Figure 1E). No gene or functional class was particularly affected, confirming a broad disruption of the network in emerging subclones (Figure 1D). We also interrogated smaller-scale datasets from glioblastoma multiforme (Suzuki et al., 2015) and kidney cancer (Gerlinger et al., 2014) patients, in which detailed phylogenetic trees were reconstructed to identify subclonal mutations selected over time in evolving tumors. Confirming the patterns observed in NSCLC, subclonal mutations affected a range of epigenetic regulators of diverse function, with up to ten alterations independently selected in individual tumors (Figures 1F, S1D, and S1E). These observations suggest that expanded subclones are often characterized by a disrupted epigenetic regulatory network and that mutations in distinct but functionally related proteins may converge into similar favorable phenotypes that are selected during cancer evolution. We therefore searched for cellular traits gained by epigenetically disrupted cells that could explain the mutational patterns observed in patients.
Disruption of the epigenetic regulatory network enhances cell fitness under environmental stress In physiology, epigenetic mechanisms mediate the cellular response to environmental cues. We therefore asked whether disruption of epigenetic control in cancer cells may affect their interactions with the tumor microenvironment (TME). We selected two cancer types originating from distinct lineages: NSCLC lung (legend continued on next page) ll OPEN ACCESS Article carcinoma, of epithelial origin, and melanoma, of melanocytic origin, and obtained cells from respective patient-derived-xenograft (PDX) models (LXFL 1674 andMEXF 2090) (Table S2). The lines express 275 and 278 epigenetic regulators, respectively, and the absence of nonsense mutations in the corresponding genes make them suitable systems to dissect the functional consequences of disrupting the network through experimental gene inactivation (Tables S1 and S2). To minimize genetic heterogeneity, we isolated clonal populations from each model and avoided prolonged culture. A common challenge faced by cancer cells is the reduced availability of nutrients, especially glutamine, in the poorly vascularized TME (Jiang et al., 2019). As expected, both melanoma and lung cancer cells suffered glutamine deprivation, showing a marked decrease in proliferation rate over 3 days ( Figure S2A). A comparable effect was also observed when glucose was depleted, indicating that reduced levels of distinct nutrients elicit similar cellular responses ( Figure S2B). Despite the apparent cytostasis, the populations were not arrested and reached a steady state where approximately 20% of cycling cells counteracted the death of over 60% of the population (Figure S2C). Prolonged culture revealed the presence of rare stress-resistant cells that reconstituted a growing population ( Figure S2D). For both cancer models, replicate populations grown in independent wells showed a range of responses to nutrient deprivation, reflecting varying numbers of resistant colonies ( Figure S2E). The stochastic emergence of stress-resistant subpopulations, together with the clonal nature of the lines, suggests a nongenetic basis for the observed phenotype. In agreement, treatment of cells with a panel of epigenetic inhibitors altered the response to stress and promoted resistance ( Figures S2F and S2G). Multiple compounds targeting histone modifiers with either repressive or activating functions increased cell fitness under stress, with treated populations showing up to 5-fold more cells than control populations after 9 days (Figures S2F and S2G). Transient inhibition of epigenetic regulators prior to stress was not sufficient to confer the stressresistance phenotype, which only emerged upon sustained treatment ( Figure S2G). The effect was only observed under stress and in fact counteracted a decrease in fitness induced by some drugs in unperturbed conditions (Figures S2F and S2G). Thus, cancer cells with a functional network of epigenetic regulators show a heterogeneous response to nutrient deprivation, with most cells succumbing to stress and a small subset surviving and being selected over time. This functional heterogeneity can be modulated experimentally, and interference with epigenetic regulation can benefit the cell population under stress.
To mimic on a large scale the widespread LOF mutations observed in patients, we systematically inactivated epigeneticrelated genes in both cancer models using an arrayed singleguide RNA (sgRNA) library that targets regulators of chromatin and DNA methylation across 16 functional classes (Figures 2A and 2B; Table S1) (Henser-Brownhill et al., 2017). The approach enables efficient editing in polyclonal populations and is compatible with microscopy-based readouts ( Figures S3A and S3B). The aim of this analysis was not to perform a screen but to model the broadly disrupted epigenetic network that is selected in patients and dissect its functional consequences. We generated knockout (KO) cell populations in 96-well format, including on each plate negative control samples in which non-expressed genes were targeted ( Figure S3C). We grew cells in either unperturbed conditions or under glutamine deprivation over 7 days, and assessed the stress-specific effect of gene KO on cell fitness ( Figures 2B and S3C). To avoid inflated ratios, KO populations with severely compromised fitness under unperturbed conditions were excluded from the analysis, leaving 262 and 250 populations for MEXF 2090 and LXFL 1674, respectively (Table S3). Biological replicates were highly concordant, allowing detection of even mild biological effects ( Figure S3D). Assessment of negative controls showed that interwell technical variability was controlled for (Figures 2C,2D,and S3E). While unperturbed MEXF 2090 KO populations showed fitness values largely within the distribution of negative controls, their fitness was broadly enhanced under stress, with 91 KO populations exhibiting a favorable phenotype (Z score > 1.64) (Figures 2C-2F  and Table S3). As expected, targeting of related subunits of protein complexes induced consistent phenotypes, whereas KO of lowly expressed genes had no effect ( Figure S3F). Furthermore, validation experiments confirmed accurate detection of enhanced fitness ( Figure 2G and Table S3). This analysis also revealed favorable phenotypes that only emerged at later time points, suggesting that the short time frame of the large-scale assays may have underestimated the number of KO populations with enhanced fitness ( Figure 2G). We observed both shared and cancer-type-specific effects ( Figure S3G and Table S3). Importantly, the stress-resistant phenotype was not restricted to KO (C and D)  populations defective in a specific molecular function but was shared across all layers of epigenetic regulation in both cancer models ( Figure 2F), with patterns resembling the distribution of LOF mutations observed in patients ( Figures 1A and 1B).
Notably, 89% and 84% of the inactivated genes that conferred enhanced fitness in MEXF 2090 and LXFL 1674 cells, respectively, were affected by LOF mutations in pan-cancer cohorts, and 38% and 32% were subclonally mutated in the TRACERx cohort. Within the subset of epigenetic regulators with an established role in cancer , 52% were associated with the stress-resistance phenotype, including genes often mutated in expanded subclones (e.g., ARID1A and SETD2) ( Table S4). KO of other frequently inactivated genes such as KDM6A, KMT2C/ D, or CTCF did not enhance cell fitness under stress, suggesting that their loss benefits cancer cells via alternative mechanisms (Table S4). The stress-resistant phenotype of epigenetically disrupted cells was not restricted to melanoma and NSCLC, and cells derived from PDX models of colon, pancreatic, and bladder cancer also exhibited enhanced fitness under stress upon pharmacological interference with multiple epigenetic regulators ( Figure 2H). As observed with genetic inactivation and in agreement with the absence of recurrent subclonal mutations in patients (Dentro et al., 2021), we detected both shared and model-specific effects. Thus, independently of context-specific effects of individual genes, the consequences of broadly disrupting the epigenetic regulatory network are conserved across multiple cancer types. To begin to characterize the stress-resistant phenotype, we selected KO populations that showed varying degrees of enhanced fitness ( Figure 2G, dashes) and belonged to various classes of regulators. Increased percentages of proliferating cells and lower percentages of apoptotic cells 3 days after nutrient deprivation suggested that more cells in the population survived and continued to proliferate under stress ( Figure S3H). In agreement, clonogenic assays showed a higher number of independent stress-resistant subpopulations after 12 days ( Figure 2I). Thus, the enhanced fitness of KO populations results from an increased frequency of surviving cells shortly after stress.
To exclude that the observed phenotypes were due to a decreased dependency of epigenetically disrupted cells on glutamine, we challenged MEXF 2090 cells with another type of stress often faced by cells during tumor growth: hypoxiainduced acidification of the TME (Yoo et al., 2020). More than 60 KO populations were resistant to both types of stress (Figure 2J and Table S3), and acidity-resistant phenotypes were also shared across all functional classes ( Figures 2K and 2L). Since nutrient deprivation and low pH have antithetic consequences on cell metabolism (Yoo et al., 2020), this response is likely not due to alterations in specific molecular pathways but rather to a more general ability to survive in a hostile milieu. Of note, KO populations did not show broadly enhanced fitness in response to DNA-damage-related stress such as hydroxyurea-induced replication defects (median Z score: 0.23 versus 0.74 for nutrient deprivation and acidity) ( Table S3), suggesting that disruption of epigenetic control mainly increases the ability of cells to cope with environmental challenges. We conclude that inactivation of numerous and diverse epigenetic regulators converges into a common phenotype that enables cancer cells to survive and proliferate in harsh environments.

Stress-dependent selection of epigenetically disrupted cells in the tumor microenvironment
We then examined the behavior of epigenetically disrupted cells in the presence of competing wild-type cells. Co-culture experiments assessing the previously selected KO populations indicated a broad selection of KO cells under stress ( Figures 3A and 3B). We then chose EZH2 as a representative regulator frequently inactivated in cancer ( Figure S1B) (Laugesen et al., 2016) and performed in vivo competition assays. Although nutrients and oxygen are often limiting in unperturbed tumors, we induced further deprivation by treatment with bevacizumab, an anti-vascular endothelial growth factor A antibody that inhibits angiogenesis (Ferrara, 2005) ( Figure 3C). Reduced levels of phosphorylated S6 ribosomal protein (pS6), a readout of mammalian target of rapamycin (mTOR) activity that correlates with nutrient availability (Sabatini, 2017), indicated that the treatment was effective ( Figure 3D). EZH2-deficient cells consistently outcompeted control cells in all bevacizumab-treated tumors, independently of host cell abundance, with samples containing up to 95% of EZH2-KO cells (Figures 3E,S3I,and S3J). Immunostaining of tumor sections confirmed the selection of EZH2 LOF during tumor growth, as most cells were depleted of the EZH2-deposited H3K27me3 mark ( Figure 3F). EZH2-KO cells were particularly enriched in pS6-negative tumor regions-those most severely affected by bevacizumab treatment-suggesting that intratumor heterogeneity in the TME influences the expansion of epigenetically disrupted cells ( Figure 3G). Vehicle-treated tumors showed a bimodal behavior: tumors exhibiting uniformly high pS6 levels, indicative of a nutrient-rich TME, showed comparable abundance of the two cell populations, whereas EZH2-KO cells outcompeted control cells in tumors showing reduced pS6 levels ( Figures 3E and 3H). Thus, basal levels of nutrient deprivation in unperturbed tumors are sufficient to promote the selection of EZH2-KO cells, and this effect is maximized under conditions of enhanced stress. Collectively, the evidence from cellular and mouse models, combined with the prevalence of epigenetic-related subclonal mutations in cancer patients, suggests that the increased ability of epigenetically disrupted cells to withstand environmental stress promotes their selection during cancer evolution.
Disruption of epigenetic control does not promote genetic heterogeneity or cell state transitions Multiple cellular mechanisms could underpin the stress-dependent selective advantage of epigenetically disrupted cells (Figure 4A). One possibility is that inactivation of chromatin and DNA modifiers may promote genomic instability and generate genetically defined stress-resistant subclones. In this scenario, restoring proper epigenetic control after an initial perturbation would not alter the emerged resistance. To test this possibility, we selected epigenetic inhibitors that confer resistance to stress ( Figure S2G) and treated cells to induce reversible deregulation under nutrient deprivation; we assessed cell fitness after 9 days, washed out the drugs, and examined the response after further growth under stress ( Figure 4B). At each step, all populations were seeded at identical density as the initial seeding to allow comparison of growth rates across conditions. In all cases, stress-resistant phenotypes were reversed upon compound withdrawal, demonstrating that the observed resistance is not genetically encoded ( Figure 4B).  Figure S4 and Table S5. Alternatively, a more permissive epigenetic landscape in KO populations may enhance phenotypic plasticity by favoring transitions between cellular states. This could either entail acquisition of an alternative, favorable cellular state or reversal from a stressed to the initial state ( Figure 4A). Consistent with enhanced plasticity, epigenetically disrupted populations reversibly adjusted their growth rate in response to fluctuating environments maintaining consistently higher fitness than control cells, both when facing multiple rounds of nutrient deprivation and when challenged sequentially with distinct types of stress ( Figures 4C and 4D).
A third mechanism that we considered, which could also explain the above results, is phenotypic inertia ( Figure 4A). If cells are unable to efficiently respond to unfavorable conditions by halting proliferation or activating an apoptotic program, their ability to withstand stress in the short term may be enhanced, increasing the chance of overcoming transient challenges or acquiring secondary adaptive traits. While conceptually antithetic to cellular plasticity, phenotypic inertia would also be consistent with the observed patterns, as the selection of stress-tolerant cells during environmental fluctuations could promote survival and growth of the population under multiple unfavorable conditions ( Figures 4C and 4D).
To begin to assess whether epigenetically disrupted cells are plastic or inert, we searched for stress-responsive pathways. We profiled cell transcriptomes prior to stress (day 0 [d0]) and 12 days after nutrient deprivation (d12), comparing control cells and 5 of the 17 KO populations characterized in previous assays. Again, we chose genes that belong to distinct functional classes and that exhibited varying degrees of stress resistance (Figure S4A): the Polycomb repressive complex 2 (PRC2) subunits EZH2, EED, and SUZ12 (Laugesen et al., 2016), the chromatin remodeler SMARCD1 (Carlson and Laurent, 1994), and the linker histone H1B, an integral component of chromatin (Scaffidi, 2016). Control cells responded to stress by downregulating fitness signatures that support cell metabolism, proliferation, and survival, including c-MYC target genes, cell-cycle regulators, mTOR signaling targets, and genes involved in oxidative phosphorylation (OXPHOS) in the mitochondria, the primary source of cellular energy ( Figure S4B and Table S5). Concomitantly, stress signatures such as nuclear factor kB or p53 targets and genes driving apoptosis were upregulated ( Figure S4B). Similar changes were observed when cells grew in acidic environments, indicating that distinct triggers induce overall similar stressed phenotypes ( Figure S4C). In the KO populations, stress-responsive genes were only mildly affected, correlating with the degree of resistance observed in the fitness assays, which confirms a less-stressed phenotype ( Figures S4D and  S4E). The milder response in KO populations was not due to distinct basal states in unperturbed conditions, as indicated by comparable levels of stress-responsive genes or energy production and unaltered cell-cycle profiles ( Figures S4F-S4H). The identification of affected molecular pathways allowed us to test the plasticity model by examining cellular transitions at the single-cell level. Using an established fluorescence resonance energy transfer (FRET)-based biosensor that allows monitoring of metabolic states in living cells (Kondo et al., 2021), we compared relative OXPHOS levels in epigenetically disrupted and control cells. To avoid the presence of unedited cells that would confound the single-cell analysis in KO populations, we selected EZH2 for this analysis and used the specific inhibitor GSK126 (EZH2i) to homogeneously block its function. In agreement with previous reports (Kondo et al., 2021), a low FRET signal indicated high mitochondrial activity, as confirmed by tetramethylrhodamine ethyl ester perchlorate (TMRE) staining, with EZH2i-treated stressed cells displaying higher OXPHOS levels than control cells, as predicted by the RNA sequencing (RNAseq) analysis ( Figures S4I and S4J). One day after nutrient deprivation, control cells acquired a stressed phenotype with markedly decreased OXPHOS levels ( Figure 4E). In contrast, most EZH2i-treated cells only displayed minor changes ( Figure 4E). To directly examine cellular transitions, we tracked individual cells by time-lapse imaging. While control cells progressively transitioned from a fit OXPHOS high to a stressed OXPHOS low state, EZH2i cells were more static and retained overall higher OXPHOS levels ( Figures 4F and 4G). This cell state transition is the immediate response to changing environments and precedes effects on cell proliferation and viability. These observations provide initial evidence that epigenetically disrupted cells do not display enhanced plasticity upon stress and in fact are less responsive to changing environments.

Phenotypic inertia of epigenetically disrupted cells
We then performed single-cell RNA-seq (scRNA-seq) analysis of nutrient-deprived cells at four time points: unperturbed (d0), early stress response (d1, preceding detectable phenotypic changes, and d2, showing initial reduction in fitness but unaffected cell viability), and after selection of resistant cells (d12) ( Figure 5A). We minimized technical variability by multiplexing populations and by performing a reversed time course that enabled simultaneous sequencing of all samples ( Figure S5A) (see STAR Methods). Quantification of low-quality (LQ) cells, which include highly stressed cells, showed the expected increase at d12, with lower fractions detected in the KO populations ( Figures S5B and S5C). Since LQ cells were removed for downstream analysis, detected differences between control and KO populations are an underestimation of the real effects. Control cells showed immediate downregulation of fitness signatures by d2 and a concomitant upregulation of stress pathways, which peaked at d12 (Figures S5D and S5E; Table S5). The late stress response (differentially expressed genes [DEGs] between d2 and d12) also included a partial reversal of stress-induced changes, which confirms the selection of stress-resistant cells (Figures S5F and S5G). Cells from different time points segregated from one another in control cells, as assessed by uniform manifest approximation and projection (UMAP)-based dimensionality reduction, whereas they were less separated in epigenetically disrupted cells and subsets of d1, d2, and d12 cells clustered with d0 populations (Figures 5B, S5H, and S5I).
(F) Relative changes in fitness pathway scores over time upon stress. The median value is plotted for each cell population. (G) Fitness assays examining control and KO populations of MEXF 2090 cells grown under the indicated types of stress. Values are mean ± SEM from three biological replicates. p values from one-way ANOVA. See also Figure S5 and Table S5. Since all cell populations at d12 contained a mix of stressed and resistant cells, in varying ratios, we performed clustering analysis and looked for clusters enriched for KO cells (mainly resistant subpopulations) or for control cells (mainly stressed subpopulations) ( Figure 5C). Pairwise DEG analysis was then performed to identify signatures defining each subset of cells, and signatures clustered to identify common meta-signatures. The aim of this analysis was to search for possible subsets of adapted cells, i.e., displaying activation of pathways that allowed them to overcome stress. Signatures defining KO-enriched clusters were enriched for fitness signatures and depleted of stress signatures, confirming that they contained stress-resistant cells ( Figure 5D). However, these cells did not show activation of pathways that could drive rewiring to an adapted state, such as reprogramming to alternative differentiation states or upregulation of compensatory metabolic pathways. Furthermore, analysis of fitness and stress genes at earlier time points revealed progressive, although mild, changes in KO cells ( Figures 5E and 5F), ruling out the possibility that they had restored the initial state after reaching a fully stressed state ( Figure 4Aii). In agreement with the live cell imaging analysis, these observations indicate that the stress-resistant phenotype of epigenetically disrupted cells is not due to a more efficient rewiring of cellular programs that favors cell state transitions. In fact, KO cells change less than control cells and fail to reach a fully stressed state.
The phenotypic inertia model predicts that the inability of epigenetically disrupted cells to respond to changing environments should confer tolerance to diverse types of stress as long as they impact cell fitness via transcriptional changes. Indeed, in addition to metabolic stress induced by nutrient deprivation and acidity, the selected KO populations were also resistant to thermal stress, which induces global transcriptional alterations (Aprile- Garcia et al., 2019), and buparlisib, an inhibitor of phosphatidylinositol 3-kinases (PI3Ki) that blocks cancer cell proliferation by interfering with oncogenic signaling and its downstream transcriptional effects (Maira et al., 2012) ( Figure 5G).
We conclude that epigenetically disrupted cells fail to efficiently rewire gene expression programs in response to unfavorable conditions, thereby maintaining high levels of genes that support cell growth and low levels of genes that promote cell death. This blunted stress response results in enhanced tolerance to multiple types of environmental stress.
Chromatin-mediated changes in global transcriptional activity in response to stress To uncover how inactivation of numerous and diverse epigenetic regulators converges into phenotypic inertia, we examined general features of transcriptional activity that could be broadly affected by disruption of epigenetic control. Gene transcription is discontinuous and occurs in pulses of activity interspersed by relatively long intervals of inactivity, resulting in the production of RNA molecules in discrete bursts (Tunnacliffe and Chubb, 2020). The frequency and size of transcriptional bursts can be inferred from scRNA-seq data using established methods (Larsson et al., 2019), and genome-scale profiles of burst parameters provide a readout of global transcriptional activity. We used this approach to assess (1) how environmental stress affects transcriptional activity and (2) how disruption of epigenetic control alters this response. As expected, individual genes showed a wide range of burst frequencies in unperturbed control cells, with burst intervals of few minutes for some genes and several hours for others (Larsson et al., 2019) ( Figure S6A). However, the profile underwent substantial changes upon stress, with burst frequencies progressively decreasing at d1 and d2, and showing signs of recovery at d12 ( Figures S6A and S6B). While the reduction affected most genes, those characterized by higher burst frequency exhibited the strongest changes ( Figure S6B). We therefore examined whether high-frequency genes (HFGs; K on > 4, Table S6) included functionally related gene sets. Strikingly, multiple fitness signatures (OXPHOS and cell-cyclerelated genes, c-MYC targets, and mTOR signaling) were strongly enriched among HFGs, with false discovery rate q values up to 10 À150 ( Figure S6C). Indeed, entire sets of fitness genes showed particularly high frequencies in unperturbed cells, which strongly decreased upon stress ( Figures S6D and S6E). Thus, intrinsic features of fitness genes make them particularly sensitive to the global changes in transcriptional bursting induced by stress and explain their coordinated behavior, independently of their specific regulatory networks. In line with reduced burst frequencies, fitness genes displayed up to 95% reduction in nascent RNA levels at d2 ( Figure S6F), and chromatin immunoprecipitation of RNA polymerase II phosphorylated on Ser2 (p-Ser2 Pol II) confirmed inhibited transcription of fitness genes but not of stress genes ( Figure S6G). Stress genes, which are low-frequency genes, did not show substantial changes in burst frequency ( Figure S6H). However, their burst size progressively increased over time, likely as a result of binding of stress-activated transcription factors to their promoter (Larsson et al., 2019), showing a distinct behavior compared with all other expressed loci ( Figure S6H). Thus, a normal response to stress entails a global reduction in transcriptional activity, which affects primarily high-frequency fitness genes, and upregulation of low-frequency stress genes via increased burst size. Importantly, these effects precede changes in cellular behavior.
(D) Pairwise comparison of chromatin accessibility at d0 and d2 after stress. Each dot is an ATAC-seq peak detected in merged replicates, at gene promoters or active enhancers. RPKM, reads per kilobase per million.  Table S6. KO populations showed distributions of burst frequency similar to those of control cells in unperturbed conditions (Figure 6A and S7A). However, the response to stress was strongly affected, and by d2 cells had only undergone mild changes ( Figures 6A-6C). When compared with control cells, the distribution of burst frequency was significantly different in KO populations, with numerous HFGs, including fitness genes, retaining high frequency upon stress (Figures 6A-6C and S7A). In agreement, KO populations only showed a minor decrease in p-Ser2 Pol II occupancy at fitness genes, indicating that cells were unable to efficiently inhibit transcription upon stress ( Figure S7B). Similarly, stress genes displayed largely unaltered patterns of burst size in the early response and only a mild increase at d12 ( Figure S6H). Analysis of chromatin remodeling by assay for transposase-accessible chromatin sequencing (ATAC-seq) confirmed that KO populations were refractory to stress-induced changes. Nutrient-deprived control cells exhibited global alterations in chromatin accessibility, which disproportionately affected regions with high basal accessibility, as observed for transcriptional burst frequencies ( Figures 6D-6F and S7C). Interestingly, the genome-wide reduction in burst frequencies was associated with a global increase, rather than a decrease, in chromatin accessibility and was not accompanied by substantial changes in active or repressive histone marks ( Figures 6D, 6E, 6G, and S7D-S7F). These patterns resemble those associated with UV-induced transcriptional shutdown (Liakos et al., 2020) and suggest that stress-induced transcription inhibition is mediated by mechanisms distinct from canonical modalities of gene repression acting in steady-state conditions. KO populations exhibited defective chromatin remodeling, globally at promoters and enhancers, and particularly at fitness genes, which are char-acterized by high basal levels of chromatin accessibility ( Figures 6F, 6G, and S7E-S7G). Thus, the defective transcriptional response in epigenetically disrupted cells is likely due to their inability to efficiently reorganize chromatin structure across the genome in response to stress.
We finally asked whether a defective inhibition of global transcription is causal in establishing the enhanced tolerance of epigenetically disrupted cells to stress. To do so, we disrupted the negative elongation factor complex NELF, which mediates transcriptional inhibition by blocking the release of paused RNA Pol II from promoters (Aprile- Garcia et al., 2019;Rawat et al., 2021). NELFA KO phenocopied the effect of epigenetic disruption, conferring a fitness advantage under stress, which demonstrates that sustained transcriptional activity upon stress is sufficient to confer phenotypic inertia ( Figure 6H). Altogether, these observations indicate that the inability to efficiently rewire global transcription in response to stress, a property that we name ''transcriptional numbness'', underpins the enhanced tolerance of epigenetically disrupted cells to unfavorable conditions. A blunted stress response and an increased ability to survive and proliferate at early stages increase the chance for epigenetically disrupted cells to overcome a transient challenge or to find a way to cope with a sustained one. Indeed, expression of glutamine synthetase (encoded by GLUL), the enzyme responsible for de novo synthesis of glutamine (Bott et al., 2015), was detected in subsets of cells at d12 ( Figure 6I). Other genes known to compensate for glutamine deficiency, such as ASNS (Zhang et al., 2014) and SLC1A3 (Tajan et al., 2018), were also upregulated ( Figures 6I and 6J). The expression pattern of the three genes varied within the populations, indicating that individual Values are mean ± SEM from three biological replicates. *p < 0.05, **p < 0.01, ***p < 0.001, and ****p < 0.0001 (Student's t test).
cells employ various strategies for long-term survival ( Figure 6J). The same genes were also detected in the few surviving control cells, suggesting that stress selects common states independently of their final abundance ( Figure 6I). We conclude that phenotypic inertia in the early response to stress increases the frequency of sporadic events that promote long-term adaptation and results in a competitive advantage at the population level (see graphical abstract).

Hypersensitivity of phenotypically inert cells to CDK9 and BET inhibitors
Independently of the selective advantage it may confer, disruption of cellular homeostasis in cancer cells often creates vulnerabilities. Since a defective transcriptional rewiring underpins the observed stress-resistant phenotypes, we hypothesized that further interference with transcription may allow targeting of phenotypically inert cells. Indeed, stress-resistant cells selected after growth in nutrient-deprived conditions were hypersensitive to two CDK9 inhibitors (CDK9i) and one BET inhibitor (BETi) compared with the whole population ( Figures 7A and 7B). Of note, a similar sensitivity was observed in all populations independently of the extent of stress resistance and their basal sensitivity, suggesting again that a common cellular state had been selected. Confirming the effect of the inhibitors, genetic inactivation of BRD2, BRD4, and BRD8 resulted in a rare stress-hypersensitive phenotype in largescale fitness assay (Table S3). Since therapeutic windows exist for both classes of compounds (Bradner et al., 2017), these findings suggest a strategy to target stress-resistant cells that drive cancer evolution.

DISCUSSION
Cancer cells face many environmental challenges during tumor development, and the immediate response to these stresses dictates whether a cell will find a way to cope with it or die. Evolutionarily conserved mechanisms, including sensors of stress and effectors, orchestrate this response (Fulda et al., 2010). We show that cells deficient in a broad range of epigenetic regulators exhibit a defective transcriptional response to environmental stress, which increases their tolerance to harsh conditions and confers a competitive advantage during tumor growth. Our results provide a mechanistic explanation for the widespread selection of subclonal mutations affecting epigenetic regulators in patients ( Figures 1D-1F). Phenotypic inertia of epigenetically disrupted cells has a dual beneficial effect: first, by failing to halt proliferation, stress-tolerant cells can propagate the favorable trait to their progeny and increase the size of the selectable pool; second, by raising the threshold at which individual cells suffer from stress, it delays the activation of an apoptotic program. As a combined result, the chance of long-term adaptation increases. Indeed, we find that stresstolerant cells employ various strategies to compensate for the lack of glutamine. Notably, the same strategies are also used by rare control cells, but a higher number of stress-tolerant cells in epigenetically disrupted populations increases the frequency of these sporadic events, resulting in overall enhanced fitness at the population level. The observed tolerance to environmental stress is reminiscent of the increased tolerance to DNA damage induced by loss of p53 (Vousden and Lu, 2002). In both cases, a faulty response to either cell-extrinsic or cell-intrinsic stress disables safeguard mechanisms that normally protect cell populations and eliminate unfit cells, thereby turning an LOF into a survival advantage. Phenotypic inertia may extend beyond mutated epigenetic regulators and contribute to the selection of cells lacking other classes of genes involved in sensing or responding to stress. As a proof of principle, we observed enhanced fitness in cells lacking two major cellular sensors that were included in the large-scale assay as controls: the nutrientsensing mTOR kinase and the damage-sensing ATM kinase (Table S3). Furthermore, other genes frequently inactivated across cancers, such as PTEN and FBXW7, have been linked to enhanced resistance to stress via transcriptionbased mechanisms (Nowak et al., 2013;Sanchez-Burgos et al., 2022).
In addition to influencing the selection of cells in the harsh microenvironment of primary tumors, phenotypic inertia may also be relevant to other aspects of tumor biology. For instance, metastatic cells in the circulation or in foreign tissues may benefit from being insensitive to surrounding conditions (Senft and Ronai, 2016). Similarly, the inability to respond to differentiation cues may help cells maintain an undifferentiated and self-renewing state (Caré n et al., 2015). A defective transcriptional response may also promote resistance to targeted therapy that relies on gene expression changes, acting through a mechanism antithetic to phenotypic plasticity (Marine et al., 2020). Interestingly, both mechanisms have nongenetic bases, but while transcriptional numbness promotes the selection of an unresponsive phenotype, the ability to transition between cellular states favors adaptation through transcriptional rewiring. Notably, the enrichment of a favorable phenotype within a resistant population after weeks of growth can be explained by both scenarios. We postulate that epigenetically competent cells exploit both phenotypic plasticity and stochastic inertia, depending on the type of stress and their cellular state, while transcriptional numbness becomes dominant in epigenetically disrupted cells.
Underpinning the phenotypic inertia of epigenetically disrupted cells is their inability to efficiently remodel chromatin and decrease global transcriptional activity, which leads to a sustained expression of fitness genes that control cell proliferation and metabolism. The observation that loss of epigenetic regulators with diverse functions leads to a common defect indicates that multiple layers of epigenetic regulation are required to execute this complex transcriptional response. Several parallelisms exist between the observed effect and the transcriptional response to heat shock or UV irradiation, with seemingly unrelated stressors leading to a similar inhibition of global transcription. Interestingly, transcriptional downregulation in nutrient-deprived cells is accompanied by a genome-wide increase in chromatin accessibility, and is independent of histone modification-patterns strikingly similar to those observed upon UV irradiation (Liakos et al., 2020). Synchronous ''opening'' of regulatory regions across the genome may lead to RNA Pol II or other co-factors becoming limiting, affecting particularly highly expressed fitness genes that are transcribed at high frequency. More accessible promoters may also be the result of reduced nucleosome turnover at gene bodies upon stress (Teves and Henikoff, 2011) or of an accumulation of proximally paused Pol II due to an active block of the transition to productive elongation (Aprile- Garcia et al., 2019). Future studies will be required to characterize the molecular underpinning of this noncanonical modality of gene repression and to clarify how different classes of epigenetic regulators contribute to the global chromatin reorganization induced by stress.
In conclusion, our study uncovers a mechanism through which disruption of epigenetic control confers a selective advantage to cancer cells. We propose that phenotypic inertia shapes the evolution of primary tumors harboring mutations in epigenetic regulators and, potentially, other processes that exert a selective pressure on cancer cells through environmental cues.

STAR+METHODS
Detailed methods are provided in the online version of this paper and include the following:

INCLUSION AND DIVERSITY
We support inclusive, diverse, and equitable conduct of research.    Table S2 PDX LXFL 1674 derived cells This paper Table S2 PDX PAXF 1997 derived cells This paper Table S2 PDX BXF 1352 derived cells This paper Table S2 PDX CXF 269 derived cells This paper  Table S7 Primers used in ChIP-qPCR assays This paper Table S7 Primers used for cloning This paper Table S7 Primers used for library preparation and NGS for the in vivo competition experiment This paper Table S7 Synthetic construct carrying an engineered sgRNA scaffold with a 3 0 terminal Capture Sequence 1 and is compatible with the 10x Genomics Feature barcoding technology This paper Table S7 (Continued on next page) ll OPEN ACCESS

Article
Cancer Cell 41, 1-18.e1-e14, January 9, 2023 e3 Data and code availability Z score tables from the large-scale fitness assays are available as supplemental data. Exome-seq, bulk RNA-seq, ChIP-seq, ATACseq and single-cell RNA-seq data have been deposited at NCBI-GEO (Super-Series GSE178643) and are publicly available as of the date of publication. Accession numbers are listed in the key resources table. Any additional information required to reanalyze the data reported in this paper is available from the lead contact upon request.

Patient datasets
Pan-cancer cohorts. Nonsynonymous mutations affecting genes listed in Table S1  LGG cohort. Mutational data displayed in Figure S1D was obtained from Table S1 in Suzuki e al., (Suzuki et al., 2015) focusing on nonsynonymous mutations with VAF <0.4. Phylogenetic trees from LGG and ccRCC (Gerlinger et al., 2014) were drawn based on published trees, visualizing mutations in epigenetic regulators retrieved from the corresponding supplementary data. MEXF 2090, LXFL 1674, PAXF 1997, BXF 1352, CXF 269 derived cells were derived from the corresponding PDXs obtained from the Charles Rivers tumor model compendium (Table S2). All cell lines (and their derivatives) were cultured in RPMI 1640 (ThermoFisher Scientific, #22400089) with 10% fetal bovine serum (FBS), 100 U/ml penicillin, and 100 mg/mL streptomycin at 37 C in 5% C0 2 . To enable efficient CRISPR-induced KO, MEXF 2090 and LXFL 1674 cells were transduced with a lentiviral pCW-Cas9 vector (Monserrat et al., 2021). Following a 7-day selection with 600 mg/mL hygromycin, clones with minimal background levels of Cas9 and responsive to induction with 1 mg/mL doxycycline were isolated.

Mice
NSG mice (NOD.Cg-Prkdc scid Il2rg tm1Wjl /SzJ) used in the in vivo competition experiment were obtained from the common Francis Crick colony. All animal studies were performed in accordance with the Francis Crick Institute Animal Welfare and Ethical Review Body and regulation by the UK Home Office project license PPL 70/8167 and PC2165EA4. Test naive male mice were used at 11-13 weeks of age for the experiment. Male mice were used to avoid possible interference of female menstrual cycle. The mice were housed in well-ventilated cages at a constant temperature and humidity (23 C ± 2 C, 50-60%) in a pathogen-free controlled environment, with a standard 12-12 h light-dark cycle, and food and water were provided ad libitum. The weight and signs of diseases were monitored in regular intervals.

Generation of KO lines
To generate cell lines expressing sgRNAs, cells were transduced with pLenti_BSD_sgRNA constructs (Henser-Brownhill et al., 2017) and selected with blasticidin (4 mg/mL) for 5 days. To induce gene KO, Cas9 expression was induced by the addition of doxycycline (1 mg/mL). sgRNA constructs were sourced from an available arrayed lentiviral sgRNA library (Henser-Brownhill et al., 2017). Lentiviral particles were produced by transfecting HEK293T packaging cells, in 96 well plates, at 80% confluence with 67.5 ng of the sgRNA lentiviral construct, 50.6 ng of psPax2 packaging plasmid (Addgene, #12260), 16.8 ng of pMD2G envelope plasmid (Addgene, #12259) and 15 ng pAdVantage (Promega, #E1711) at a ratio of 3:1 FugeneHD to DNA (Promega, #E2311). Approximately 24 and 48 h after transfection, the supernatant containing viral particles was recovered, filtered through a 0.45 mm filter (Millipore, #MSHVS4510), pooled together, and frozen at À80 C. After selection, expansion and induction of KO for 10 days, the libraries of KO populations were stored at À80 C in multiple aliquots. To freeze cells in 96-well format, cells were detached from the plate using 30 mL of trypsin per well. Following addition of 80 mL of FBS containing 10% DMSO, plates were sealed and stored at À80 C for up to 4 weeks. To thaw cells, plates were placed in a water bath at 37 C for a few seconds and spun for 5 min at 4 C after addition of 50 mL of medium to each well. Fresh medium (100 mL) was finally added to each well after removal of 120 mL of freezing medium. Isolated clonal KO populations were used for ChIP and ATAC experiments, after confirming that stress-resistant phenotypes were similar to the polyclonal populations. Experiments involving knock-out of NELFA were performed by transfecting synthetic guide RNAs (Edit-R crRNA, Horizon) (Table S7) in Cas9-expressing MEXF 2090 cells. EZH2-and TNP2-targeting crRNAs were used as positive and negative control, respectively. Each CRISPR RNA (crRNA) and the trans-activating CRISPR RNA (tracrRNA) were resuspended in 53 siRNA buffer (Horizon, #B-002000-UB-100) at a concentration of 20 mM. Each crRNA was then mixed with an equal amount of tracrRNA and the mix diluted 1:100 in Optimem medium (Gibco, #31985047). A total of 10 mL of the crRNA/tracrRNA mix was added to a well of a 96-well plate (to achieve a final concentration of 20 nM in a 100 mL of final volume) and mixed with 10 mL of transfection reagent (DharmaFECT Duo: 0.3 mL/96-well (Horizon, #T-2010-01)) diluted in Optimem. After 15 min, 4,000 cells resuspended in 80 mL of complete medium containing 1 mg mLÀ1 doxycycline were added to each well. After 2 days of transfection, cells were split 1:30 and plated in 6 replicate plates. Three days after transfection, cells from three plates were stained with SYTOX Green (d0) (ThermoFisher Scientific, #S7020 and #S11380 respectively) and medium was changed to L-glu-deprived medium (ThermoFisher Scientific, #21870076). At day 9, cells were fixed and stained with SYTOX Green (d9). Fitness was measured as fold change between d9 and d0.
Large-scale fitness assays Plate layout: From the available sgRNA library (Henser-Brownhill et al., 2017), constructs targeting 318 genes, encoding stricto sensu epigenetic regulators were selected and arranged in eight 96-well plates (Corning, #3596). Each plate consisted of 40 KO populations in which distinct epigenetic regulators were inactivated, along with 20 negative control populations transduced with sgRNAs targeting 5 non-expressed genes, four replicates each. Multiple negative controls per plate were included to account for technical variability and well-effects, and to enable robust normalization of the observed phenotypes across plates and experiments. Each plate contained an ARID2-targeting sgRNA as an inter-plate standard. External wells were excluded and filled with PBS to avoid edge effects.
Assay pipeline: Multiple identical replicates, sufficient for the different treatments, were generated from each library plate. When the median cell count/well across plates reached $4000 cells, the plates were either grown under stressful conditions (see below) or maintained in unperturbed conditions. Over the course of the experiment one representative plate was monitored to confirm that the growth kinetics were as expected. At the indicated endpoints for each condition (see below) the population fitness was assessed by quantification of cell count: cells were fixed with 4% paraformaldehyde (PFA, Alfa Aesar, #43368) followed by permeabilization with 0.5% Triton X-100 in PBS and nuclei staining with SYTOX Green Nucleic Acid Stain (ThermoFisher Scientific, #S7020). Imaging and quantification were performed using an Incucyte S3 Live-Cell Analysis System.
Stress conditions: KO and control populations were cultured in the following conditions: a) unperturbed b) glutamine deprivation c) acidic environment and d) replication stress, with each stress applied at two distinct strengths. Glutamine deprivation was sustained for 3 or 7 days after glutamine removal. Media acidification was induced by addition of HCl to a final pH of 6.7 or 6.5, which resulted in a 40-60% reduction in cell counts after 3 days compared to untreated cells. To induce replicative stress, cells were cultured in the presence of 200 mM or 250 mM hydroxyurea (Sigma-Aldrich, #H8627) which resulted in 40% or 60% in cell counts after 3 days compared to untreated cells. Cell count for populations grown in unperturbed conditions was quantified at day 2. The endpoint of each treatment was determined based on when confluence was reached in fittest populations and depended on how severely cells were affected by each stress, with glutamine deprivation being the most deleterious stress. A survival benefit in at least one strength condition for each stress was considered as enhanced fitness. Experiments assessing the cellular response to thermal stress were performed by incubating MEXF 2090 cells at 43 C for 3 h, and then at 37 C for 9 days. To assess the response to oncogenic signaling inhibition, cells were cultured in the presence of 1 nM Buparlisib (Generon, #HY-70063-5mg) for 3 days. At the indicated endpoints for each condition, plates were fixed with 4% PFA followed by nuclei staining with SYTOX Green. Imaging and quantification were performed in Incucyte S3 Live-Cell Analysis System.
Filtering: Before quantifying normalized fitness, various parameters were assessed to remove: 1) wells containing cell clumps that would affect the measurement, identified by visual inspection of all plates; 2) outliers among replicates; 3) epigenetic regulators that are lowly expressed in MEXF 2090 (Log 10 pseudo-counts <1) or LXFL 1674 (Log 10 TPM <0.5) cells; 4) KO populations with severely compromised fitness in the unperturbed condition (20% reduction compared to the plate median for MEXF, 2090 and less than 20,000 cells at endpoint for LXFL 1674). The last step was performed to avoid inflated stress/unperturbed ratios. For steps 1 and 2, we excluded from the analysis 0.85% of imaged wells (61 out of 7200 total wells). 53 wells (87% of excluded wells) were excluded because they contained visible clumps of cells that would have confounded the analysis. The remaining 8 wells were excluded because they were highly inconsistent with the other two replicates (typically differing from the two consistent replicates more than 3 times their difference).
Data analysis. For each KO or control population the normalized fitness was derived from the stress/unperturbed ratio in cell count at endpoint. KO populations with enhanced or reduced fitness were defined based on the formula: Z=(c-m)/s, where c is the fitness of individual KO populations, m is the mean fitness of negative controls, s is the SD of the fitness of negative controls. Populations exhibiting enhanced or reduced fitness were defined as those with a Z score > 1.645 or < -1.645 (90% confidence interval), respectively. Since we used a high ratio of sample vs negative controls (2:1) and the Z score includes information on confidence level, correction for multiple comparisons was not required.

Validation fitness assay
Validation fitness assays were performed in a similar way as the large-scale assays, either measuring cell count over time using replicates fixed at various time points, or at endpoint. Population growth was quantified by normalizing the average cell count at each time point to the pre-treatment count (d0). In experiments where the population growth was monitored over several weeks, smoothing was applied to the curves to account for technical noise introduced by media change. Throughout the study, selection of genes for targeted analysis was driven by three factors: (i) the need of including epigenetic regulators representative of multiple functional classes, to be able to assess network-level effects; (ii) the decision to not limit the analysis to ''top performers'' and instead assess populations exhibiting varying degree of stress-resistance, including mild phenotypes (e.g. HIST1H1B-KO cells); (iii) the complexity of the assay performed, which affected the number of populations that we could analyze in the different experiments. When using inhibitors, availability of specific compounds was an additional criterion for selection. MEXF 2090, PAXF 1997, BXF 1352 and CXF 269 cells were pre-treated for 3 days before being exposed to stress, and drugs and media were replaced every 3 days over the course of the experiment. Variable population growth rate and sensitivity to nutrient deprivation contribute to the effect size detected in each model. When assessing the reversibility of the stress-resistant phenotype, MEXF 2090 cells were grown for 9 days in the presence of the inhibitors, and for an additional 9 days in the absence of the compounds ( Figure 4B). When assessing the effect of pre-treatment, MEXF 2090 were pre-treated for 3 days, then drugs were withdrawn and cells were left to recover for 2 days before being exposed to stress ( Figure S2G)

Proliferative and apoptotic fractions
Two hours before the endpoint, populations were pulsed with 10 mM EdU and 5 mM of Ac-DEVD-NucView488. (Cen et al., 2008) Live cell imaging of Caspase 3 activity was performed using an Incucyte S3 Live-Cell Analysis System. Subsequently, cells were fixed with 4% PFA followed by EdU staining by using the Click-iT EdU Cell Proliferation Kit (ThermoFisher Scientific, #C10340) as per manufacturer's protocol, and SYTOX Green staining for nuclei quantification.
Cell cycle analysis MEXF 2090 cells were grown in unperturbed condition for 3 days. One-hour prior collection, cells were pulsed with 25 mM EdU. Cells were detached, fixed in 4% PFA and permeabilized in 0.5% Triton X-100. EdU was detected using a biotin-azide click reaction by incubating cells in 100 mM Tris HCl pH 7.6, 4 mM CuSO4, 5 mM M Sulfo-Cyanine5 Azide (Lumiprobe, A3330), 100mM Sodium Ascorbate for 5 min at room temperature. After washing, cells were resuspended in 1 mM/mL DAPI and analyzed by flow cytometry.
ATP assays MEXF 2090 were plated in two identical replicates in 96well-plate. After 3 days of growth, one plate was lysed using CellTiterGlo reagent (Promega, G9242) following manufacturer instructions. Briefly, culture medium was removed, and cells were incubated with a 1:1 mixture of CellTiterGlo reagent and medium with gentle shaking. After 10 min, the solution was transferred to a white 96 well plate and luminescence was measured using a PHERAstar microplate reader. To infer the ATP levels per cell, the replicate plate was fixed with 4% PFA, permeabilized and stained with SYTOX Green Nucleic Acid Stain (Thermo Fisher Scientific, S7020). Imaging and quantification were performed using an Incucyte S3 Live-Cell analysis system to obtain the number of cells per well. ATP levels per cells were then calculated dividing luminescence intensity by the number of cells per well.

Sensitivity assays
To assess the basal sensitivity of control and epigenetically-disrupted MEXF 2090 cells to CDK9 and BET inhibitors, cells were treated for 48 h with compounds or DMSO as control: CDK9i: 5,6-Dichlorobenzimidazole 1-b-D-ribofuranoside (DRB) 15 mM (Merck, F3055-5MG) and Flavopiridol hydrochloride hydrate 1 nM (Merck, F3055-5MG). BETi: JQ1 1 mM (Merck, SML1524-5MG). At the indicated endpoint, cells were fixed with 4% PFA followed by nuclei staining with SYTOX Red. Imaging and quantification were performed in Incucyte S3 Live-Cell Analysis System. The percentage of killing was calculated by comparing drug-treated cells with DMSO-treated cells. In parallel, control MEXF 2090 cells and KO populations were cultured for 8 days in nutrient-deprived condition to select stress-resistant cells. At day 8, cells were split at equal density and grown overnight in unperturbed conditions. The following day, stress-resistant cells were treated for 48 h with CDK9 and BET inhibitors or DMSO as control and sensitivity assay was calculated as stated above.
In vitro competition assays MEXF 2090 cells were transduced with lentiviral constructs expressing GFP-NLS (pTRIP-SFFV-EGFP-NLS, Addgene #86677) (Raab et al., 2016), or mCherry-NLS (pTRIP-SFFV-mCherry-NLS). To generate pTRIP-SFFV-mCherry-NLS, mCherry was amplified from pUAS-mCherry-NLS (Addgene #87695) (Zhang et al., 2017) by PCR (forward primer: TAGCGGATCCATGGTGAGCAAGGGCGAG, reverse primer: AGGTCTCGAGCGCCCCGACTCTAGATTATAC) and cloned in the pTRIP-SFFV expression vector after GFP removal by BamHI-XhoI digestion. Seven days after infection, cells expressing similar levels of the fluorescent proteins were isolated by flow cytometry. mCherry-labelled lines were transduced with sgRNA constructs targeting epigenetic regulators, while GFP-labelled cells with sgRNAs targeting the non-expressed gene TNP2. After 10 days of blasticidin selection and induction of KO through doxycycline treatment, mCherry-labelled KO populations were mixed with GFP-labelled control cells in equal quantities and seeded in 96-well plates at a density of 3000 cells per well. A mix of GFP-and mCherry labeled control cells was used as a baseline to account for possible differential fitness of the two labeled lines. Twenty-four hours after plating, the medium was refreshed with phenol-free RPMI lacking L-glutamine (ThermoFisher Scientific, #32404014). Live cell imaging was performed at 12-h intervals using an Incucyte S3 Live-Cell Analysis System. During the course of the experiment, the media was refreshed every 3 days. After 12 days in nutrient deprivation, the mCherry to GFP ratio was calculated and normalized to the pre-treatment ratio.
In vivo competition assays MEXF 2090 cells transduced with sgRNAs targeting either EZH2 or TNP2 (non-expressed control gene) and treated with doxycycline for 10 days were mixed at equal ratio and injected in NOD.Cg-Prkdc scid Il2rg tm1Wjl /SzJ (NSG) mice obtained from the common Francis Crick colony. 5 x 10 5 cells from the mix were intradermally injected in both flanks of 11-13-week-old male NSG mice in 50 mL of PBS. Approximately 3 weeks after injection mice were randomly segregated in two groups and treated with either Bevacizumab (Stratech, #A2006-SEL-5mg; twice a week i.p. at the dose of 2 mg/kg or 8 mg/kg) or vehicle. Tumor volume was measured twice a week using electronic calipers until animals were humanely killed approximately 5 weeks after injection. At the endpoint, tumors were harvested and the relative abundance of EZH2-KO and control cells was estimated by next-generation sequencing of the sgRNAs amplified from tumors. To do so, tumors were cut into small pieces of $2-3 mm in length, transferred into gentleMACS M tubes (GentleMACS, #130-096-335) containing two volumes of ATL buffer supplemented with 1:10 proteinase K (Qiagen, #69506) and blended at high speed with a gentleMACS dissociator (RNA_01.01 Program). (Morales Torres et al., 2020) gDNA from tumors and the injected cells was extracted using a DNeasy Blood & Tissue kit as per manufacturer's protocol (Qiagen, #69506). NGS libraries were prepared by performing a two-step nested PCR using Herculase II Fusion Enzyme Kit (Agilent Technologies, #600679) and primers listed in Table S7. To ensure efficient amplification of the sgRNAs, multiple PCR reactions were run for each sample, using a maximum of 1 mg gDNA in 50-mL reactions with 20 cycles of amplification. Following the first round, the PCR product was cleaned using QIAquick PCR Purification Kit (Qiagen, #28104) and 1/50 of the reaction was used as template for the second PCR, run for 28 cycles. Final products were run on a 2% agarose gel and purified using a QIAquick gel extraction kit (Qiagen, #28706). Libraries were sequenced on the Illumina MiSeq using the MiSeq Reagent Nano Kit V2 (Illumina, #MS-102-2001) with 250 bp paired end reads and generated approximately 6000 251bp reads per sample. Raw reads were trimmed with the fastx_trimmer tool available within the FASTX-Toolkit (version 0.0.14) http://hannonlab.cshl.edu/fastx_toolkit using the parameters "-f 122 -L 141 -m 20 00 to extract the sgRNA sequence. These were then mapped to a reference consisting of the 14 guide sequences of interest using BWA (version 0.5.9-r16) (Li and Durbin, 2009) with the parameters ''-l 20 -k 4 -n 4''. sgRNA counts were obtained after filtering the mapped reads for those that had zero mismatches, and mapped to the sense strand of the guide sequence. To quantify the relative abundance of sgRNAs in each condition, raw reads for each sgRNA were normalized to the overall read counts.

Immunofluorescence microscopy
Immunostaining of cultured cells was performed using standard protocols as described in (Monserrat et al., 2021) using anti-lamin A/C (Santa Cruz Biotechnology, #sc-7292, 1:200), anti-H3K27me3 (EMD Millipore, #07-449, 1:400), anti-H4K16ac (Cell Signaling Technology, #2591, 1:500) primary antibodies and relevant fluorescent secondary antibodies. Imaging was performed using either an IncuCyte S3 system or an Axiovert Zeiss confocal microscope. For analysis of tumors, portions of tumors harvested from mice were fixed in 10% formalin and embedded in paraffin. Sections were deparaffinized with xylene and rehydrated in an ethanol gradient. Antigen-retrieval was performed for 20 min at 95 C in citrate buffer. Slides were then blocked, incubated overnight with anti-H3K27me3 antibody (1:200) or phospho-S6 Ribosomal Protein (Ser235/236) (Cell Signaling Technology, #2211) at 4 C, washed, incubated with the secondary antibodies (ThermoFisher Scientific, anti-rabbit Alexa flour 488, #A-21206 or anti-rabbit Alexa flour 568, #A10042) for 1 h at room temperature, washed 3 times with PBS, incubated with DAPI and mounted with Pro-Long Gold Antifade Mountant (ThermoFisher Scientific, #P36934). Slides stained with only the secondary antibody were used as a negative control. Stained slides were imaged using an Olympus VS120 Slide Scanner and images were processed with QuPath-0.2.2 (Bankhead et al., 2017) and quantified with ImageJ 1.45s.

Clonogenic assays
For qualitative assessment, 1500 cells were plated in 6-well plates, while for quantification 10 cells/well were seeded in 96-well plates. After 24 h, the medium was refreshed to RPMI 1640 without glutamine (ThermoFisher Scientific, #21870076) and cells grown for 12 days until visible colonies appeared. Media was refreshed every 3 days. Cells were fixed with 4% PFA, permeabilized with 0.5% Triton X-100 in PBS stained with SYTOX green and imaged using an Incucyte S3 Live-Cell Analysis System.
Live-cell imaging MEXF 2090 were modified to express the FRET-based glucose biosensor, which indicates relative preference for glycolysis and OXPHOS in living cells. (Kondo et al., 2021) The PiggyBac transposon containing the FRET biosensor (Kondo et al., 2021) was co-transected with a plasmid expressing the PB transposase (PBase) (Liang et al., 2009) at a 1:3 ratio using the FuGENE HD reagent as per manufacturer's protocol (Promega, #E2311). After 7 days, citrine-expressing cells were isolated by flow cytometry. For live-cell imaging, cells were plated in optical grade 96-well microplates (GBO, #655090). Two hours before imaging, the media was replaced with phenol-free RPMI with or without L-glutamine (ThermoFisher Scientific, #11835063 and #32404014, respectively). Cells were imaged at 15-min intervals for 24 h using an inverted Zeiss LSM 880 confocal microscope and Zeiss Zen software (v2.3).
A Plan-Apochromat 20x/0.8 NA objective lens was used and the emission signals were detected using the internal 32-channel GaAsP detector. Excitation light from an argon ion laser set to 3.5% was passed through a 458/514/561/633 multiple beam splitter and emission light was detected between 464 and 506 nm for eCFP and 517-571 nm for sensitized emission (FRET). Image acquisition settings were set to 512 x 512 pixels, zoom 0.6, 8.24 ms pixel dwell and line 4 averaging. Master gain for eCFP detection was set to 850 and 750 for citrine detection. Digital gain was set to 1.0 for both and digital offset was set to 0 for both. 3 x 2 tiling with 5% overlap followed by stitching was used to capture a rectangular field of view. The FRET ratio per cell was calculated from perinuclear areas (4x4 pixels) by dividing the total intensity in the FRET channel by the total intensity in the eCFP channel. For time-lapse single cell analysis, individual cells were manually tracked and the FRET ratio at each time point was calculated as described above. The measurements from time points overlapping with active cell divisions were discarded and replaced with the average of 6 time points (three prior and three after) flanking mitosis. The FRET biosensor is very sensitive to even minor changes in metabolic states induced by various factors, including cell density and motility. (Kondo et al., 2021) Thus, small variations in the number or distribution of seeded cells could account for the different basal level detected in Figure 4F. Other technical details, such as position of the well within the 96-well plate, which influences gas exchange and temperature equilibration after change to nutrient deprived medium, could contribute to the small initial difference. Independently of the source of variation in the baseline, 3h after media change, both samples exhibited comparable FRET signal, which increased at distinct rates. To assess mitochondrial activity, cells were pulsed with 100 nM of Tetramethylrhodamine ethyl ester perchlorate (TMRE, ThermoFisher Scientific, #T669). After 30 min, TMRE fluorescence was imaged and quantified using an Incucyte S3 Live-Cell Analysis system.
RT-qPCR and nascent RNA analysis RNA was extracted using a RNeasy Plus Micro Kit (Qiagen, #74034) according to the manufacturer's guidelines. For standard analysis of RNA levels, 500 ng of RNA was reverse transcribed using a High Capacity cDNA Reverse Transcription Kit (ThermoFisher Scientific, #4368814) as per manufacturer's instructions. cDNA was diluted 1/10 and used as input for RT-qPCR using SsoAdvanced Universal SYBR Green Supermix (Bio-Rad, #172-5274) on a CFX96 real-time PCR detection system (Bio-Rad). PPIA was used as reference gene for normalisation. For nascent RNA analysis, genomic DNA was removed first by passing cell lysates into gDNA eliminator spin columns and then by treatment of eluted RNA with TURBO DNase for 20 min (ThermoFisher Scientific, #AM2238). cDNA was prepared from 1 mg or 2 mg of extracted RNA and nascent transcripts detected using intronic primers. C1ORF43 was used as a reference gene for normalisation. Primers used for RT-qPCR in this study are listed in Table S7.
scRNA-seq analysis sgRNA demultiplexing and cell type assignment After removal of cells with empty CRISPR guide genes, a multinomial logistic regression model was applied to assign cell labels. Cells with only one CRISPR gene counts different from zero were assigned to the corresponding cell type, and used to fit logistic regression model where f i;k is the probability that cell i is assigned to cell type k, and g i;j represents the normalised j-th CRISPR counts of the i-th cell. CRISPR gene counts were normalized by dividing by the total counts of the gene across all cells. The model, fitted using caret package (Kuhn, 2008) for R, was used to predict the assignment probability of cells with two or more non-zeros CRISPR gene counts. Cells with a predicted probability greater than 95% were assigned to the corresponding label. The other cells were considered ambiguous and therefore removed from the dataset. Filtering An iterative procedure was applied to identify and remove highly-stressed, low-quality cells and lowly-expressed genes from the count table of the labeled cells at the four time points. At each iteration, the procedure involved: a) removal of genes detected in less than 20 cells across all time points, b) removal of cells with small total counts and a small number of detected genes (lower than 3 median absolute deviations (MADs) from the median) using the scater package (McCarthy et al., 2017) for R, c) removal of cells with a high percentage of mitochondrial genes counts (higher than 3 MADs from the median). (Ilicic et al., 2016) After identification and quantification of unhealthy low-quality cells in each cell population, those cells were removed from the downstream analysis.

Normalization
The filtered dataset was normalized using the scran package (Lun et al., 2016) for R. The scaling factors were firstly estimated on clusters identified by the function quickCluster, and then deconvoluted using the function calculateSumFactors. The same results were obtained by using the function computeSumFactors. The normalization affected the dataset minimally (median scaling factor across all cells: 0.89, with 50% of cells between 0.67 and 1.23), confirming that the experimental design minimized technical variability ( Figure S5A). The log 2 -transformed normalized counts were used as features for the subsequent modeling, except for the burst kinetics analysis, which used the raw data, and the variance analysis where normalization was carried out with the PAGODA algorithm. Dimensionality reduction and differential gene expression analysis Dimensionality reduction was performed using Seurat (version 3.2.2) (Stuart et al., 2019) after merging all samples into one Seurat object to allow projection of all samples onto the same UMAP space, and using stress-responsive genes defined at d1 and d2 (early stress response). For differential gene expression analysis, pre-filtered single cell samples relevant for each comparison (either all cell populations at each time-point, or one cell population across time) were merged and differentially expressed genes were identified using the FindMarkers function from Seurat. For time-point comparisons, all time points were compared to the d0 sample. In addition, pairwise comparisons between consecutive time points were made (d2 vs d1, d12 vs d2) to assess progressive changes or reversal of gene mis-expression. For cell-population comparisons, all KO populations were compared to control cells within each time point. These comparisons were repeated using the glmGamPoi Biocnductor package (1.2.0) (Ahlmann-Eltze and Huber, 2021) correcting for cell cycle phase which was calculated using the CellCycleScoring function and cell cycle phase specific genes provided by Seurat. 100 cells were sampled from each group prior to analysis to speed up the computation. Statistical significance of the differences was assessed by non-parametric Wilcox test using default parameters followed by Bonferroni correction for multiple testing (based on all tested genes). Differentially expressed genes were defined as genes with an FDR <0.01 between conditions using cellcycle corrected data.

Cell clustering and identification of KO-enriched or control-enriched subpopulations
The general procedure to determine the sub-population-specific genes consisted of two steps: firstly, each pair of KO and control cell populations at d12 were clustered based on their gene expressions, and the differentially expressed genes between cell clusters enriched with either population were determined (see Figure 5C). Secondly, the corresponding signatures defining each enriched cluster were merged through hierarchical clustering across all pairwise comparisons to define gene signatures shared by KO-or controlenriched subpopulations. More in detail, for each cell population P, the N ðPÞ 3M gene expression matrices X ðPÞ were defined, where N ðPÞ represents the size of P, and M represents the number of genes. For each KO population, we then performed a pairwise modeling with the control cells using the N ðd12Þ 3M 0 matrix, X ðd12Þ , where N ðd12Þ = N ðKOÞ + N ðcntrÞ , consisting of the expression of the most variable genes (VGs) within the merged cell populations. VGs were defined by applying the method described in Brennecke et al. (2013) to the normalized (not log2-transformed) gene expression values. Principal Component Analysis (PCA) dimensionality reduction was applied to the merged cell cycle-corrected VG data to generate the feature matrix used for each comparison. A number d < M of PCA dimensions were retained, corresponding to the elbow of the variance explained var exp by each component. The elbow point was defined as the farthest point of the var exp curve from the line connecting the points P 0 = ð1; var exp ðPC1ÞÞ and P 1 = ðPC max ; var exp ðPC max ÞÞ, with PC max = 200. The N ðd12Þ 3d PCA scores matrix T was used to determine the cell clusters. For visualization purposes, a 2-dimensional t-SNE mapping was obtained from the scores T, with a perplexity value set to 30. PCA and t-SNE dimensionality reduction were performed using the Seurat (version 3.2.2) package for R. Cell clustering was performed via nearest-neighbors method (num. neighbors = 20, tree pruning = 1/15) (Jarvis and Patrick, 1973) followed by graph partitioning through the functions available in Seurat. Consensus clustering (Monti et al., 2003) and detection of the optimal resolution were performed with in-house developed R scripts. Enrichment analysis of each cell population in each cell cluster was performed by hypergeometric test: for each cell type, p values were calculated from the cumulative distribution function, using the phyper function for R and corrected for multiple testing using Benjamini-Hochberg method. Only clusters significantly enriched for at least one a cell population (adjusted p value < 0.05) were used in the following steps of the analysis. Gene signatures defining each KO-or control-enriched clusters consisted of DE-Gs identified by MAST (Finak et al., 2015) between each cell cluster and all other cells of the pair merged set. Genes with FDR <0.05 and a non-zero log 2 -fold-change (with a non-zero-crossing 95% confidence intervals) were considered differentially expressed and used to characterize the cells belonging to the cluster. After performing the clustering analysis on cells for all KO-control comparisons, all the identified cell-cluster-specific gene signatures were clustered together to identify DE-Gs shared by different cell-clusters. To increase the robustness of the results, signatures consisting of less than 50 DE-Gs or that corresponded to cell clusters with a size smaller than 5% of the pooled set size were removed. The gene signatures were defined as vectors where each element corresponded to the product of the point estimate of the log2-fold-change and the significance (equal to À log 10 p) of each DEG, as determined by MAST. All non DEGs were assigned a signature value equal to 0. Gene signatures were clustered using a hierarchical clustering. Spearman's correlation was used as similarity measure and the optimal linkage corresponded to the maximum cophenetic correlation coefficient (Sokal and Rohlf, 1962). The optimal number of clusters corresponded to the maximum silhouette coefficient (Rousseeuw, 1987) calculated with a varying number of clusters in f2; 3; 4; .;N sig g, where N sig is equal to the total number of gene signatures. Genes sets shared by KO-enriched or control-enriched signatures were denoted as meta-signatures to distinguish them from the cell cluster-specific signatures, and defined as the average of the signatures belonging to the same cluster. Genes defining the KO-enriched or control-enriched meta-signatures were finally analyzed by GSEA to identify affected pathways.
Pathway-score estimation Pathway scores of fitness and stress signatures were estimated as the first principal component scores of the cell cycle corrected counts of the whole dataset, comprising all cell populations and time points. Only detected genes that belong to the relevant Hallmark gene sets in the MSigDB database were used for the estimation of the principal components.