Genome-wide Screens Implicate Loss of Cullin Ring Ligase 3 in Persistent Proliferation and Genome Instability in TP53-Deficient Cells.

TP53 deficiency is the most common alteration in cancer; however, this alone is typically insufficient to drive tumorigenesis. To identify genes promoting tumorigenesis in combination with TP53 deficiency, we perform genome-wide CRISPR-Cas9 knockout screens coupled with proliferation and transformation assays in isogenic cell lines. Loss of several known tumor suppressors enhances cellular proliferation and transformation. Loss of neddylation pathway genes promotes uncontrolled proliferation exclusively in TP53-deficient cells. Combined loss of CUL3 and TP53 activates an oncogenic transcriptional program governed by the nuclear factor κB (NF-κB), AP-1, and transforming growth factor β (TGF-β) pathways. This program maintains persistent cellular proliferation, induces partial epithelial to mesenchymal transition, and increases DNA damage, genomic instability, and chromosomal rearrangements. Our findings reveal CUL3 loss as a key event stimulating persistent proliferation in TP53-deficient cells. These findings may be clinically relevant, since TP53-CUL3-deficient cells are highly sensitive to ataxia telangiectasia mutated (ATM) inhibition, exposing a vulnerability that could be exploited for cancer treatment.

Introduction

Cancer evolves by accumulation of somatic mutations leading to uncontrolled cell proliferation (Cavenee and White, 1995). Mutations that arise during carcinogenesis disrupt cellular circuits, lead cells to continuous replication and escape from cell death, activate growth-stimulating pathways, and alter the tumor microenvironment (Hanahan and Weinberg, 2011). These mutations may act in isolation but also affect each other via epistatic interactions, the investigation of which has been very challenging in human cancers owing to the genetic complexity of the disease. Prior studies aiming to uncover genetic interactions between cancer-related genes or pathways focused either on established tumors or on cancer-derived cell line models (Chen et al., 2015, Eskiocak et al., 2011, Hart et al., 2015, Rauscher et al., 2018, Sanchez-Vega et al., 2018, Shalem et al., 2014, Wang et al., 2014, Westbrook et al., 2005). In such systems, the complex genetic backgrounds of tumor cells may mask epistatic relationships between genes and pathways, thus preventing their discovery. Highly controlled genetic backgrounds can better reveal genetic dependencies, which has been particularly successful in budding yeast (Costanzo et al., 2016, Kuzmin et al., 2018, Tong et al., 2001). Similar studies in human cells have thus far focused primarily on identifying synthetic lethal interactions (Aguirre and Hahn, 2018, Horlbeck et al., 2018, Shen et al., 2017).

Here, we investigated genetic dependencies of TP53 with respect to cell proliferation and tumorigenic growth in a genetically controlled cell-based system. TP53 is a critical regulator of DNA damage and replication stress, is able to promote cell-cycle arrest and apoptosis, and is mutated in 36% of all cancers, with particular high incidence in certain tumors (e.g., over 90% in small cell lung cancer) (George et al., 2015). Cells with inactive TP53 can survive and propagate even with damaging mutations. Hence, inactivation of TP53 is often among the first hits in tumorigenesis; nonetheless, it is widely accepted that second hits are required to facilitate tumor progression (Kastenhuber and Lowe, 2017). In this study, we conducted genome-wide CRISPR-Cas9-based knockout screens to investigate secondary hits that lead to uncontrolled proliferation and tumorigenic growth in TP53-deficient cells. We identified many components of the ubiquitination/neddylation pathway resulting in enhanced proliferation when inactivated. Roles of the ubiquitination/neddylation pathway in tumorigenesis remain poorly understood, despite the central importance of this pathway in cellular homeostasis and the observation of recurrent mutations affecting this pathway in human cancers (Enchev et al., 2015, Ge et al., 2018). To characterize the interplay of TP53 and neddylation, we focused on the cullin ring ligase 3 (CUL3). CUL3 and the CUL3-associated protein adaptors KEAP1 and SPOP were recently inferred as pan-cancer driver genes by statistical association. Mechanistic details of these associations, however, remain unclear (Ge et al., 2018). Here, we uncovered that CUL3 loss in TP53-deficient cells leads to persistent proliferation due to the activation of an oncogenic transcriptional program driven by nuclear factor κB (NF-κB), activator protein 1 (AP-1), transforming growth factor β (TGF-β), and Notch signaling pathway, resulting in partial epithelial to mesenchymal transition (pEMT). In addition, CUL3 knockout cells exhibit increased levels of DNA damage and genomic instability and show increased sensitivity to ataxia telangiectasia mutated (ATM) inhibition (ATMi). These results provide insights into an epistatic relationship of TP53 and CUL3 and expose a vulnerability that may have therapeutic implications.

Results

Screen Methodology

To identify genes that promote tumorigenic growth, we used an hTERT immortalized retina pigment epithelial cell line (hereafter referred to as RPE). RPE is a non-transformed diploid cell line that carries few genetic abnormalities and is amenable to genetic perturbations. We employed two complementary approaches: (1) anchorage-independent growth screens by soft agar (3D) (Freedman and Shin, 1974, Mori et al., 2009) to discover genes that promote tumorigenic growth and (2) proliferation screens (2D) to identify genes that affect the rate of proliferation. We used wild-type (WT) and TP53−/− RPE cells and applied a whole genome CRISPR-Cas9 knockout library targeting 19,050 genes (Shalem et al., 2014). (Figure 1A).

Figure 1

Whole-Genome Knockout Screens for Proliferation and Tumorigenic Growth

(A) Overview of the genome-wide screens.

(B) Genome-wide screens identify genes whose inactivation promotes growth in 2D and 3D assays. Volcano plots depict the fold change (log2) and the q-value (−log10) for each gene. Fold changes represent sgRNA read counts between treatment and control conditions.

(C) Percentage of genes in each database that overlap with the enriched genes that are identified in our screens are plotted. Blue bars, RPE; red bars, RPE.

(D) Intersection of genes between 2D and 3D screens from curated tumor suppressor genes from the tumor suppressor gene (TSG) database.

(E) Gene ontologies and pathway enrichments of genes identified from 2D and 3D screens in RPE and RPE cells.

See also Figures S1–S3 and Table S1.

We developed an R package termed mixed-effects-model-based analysis of CRISPR screens (MEMcrispR) that enables efficient analysis of genome-wide count-based screens starting from raw sequences to data visualization (Figures S1A–S1C). The underlying data modeling is based on linear mixed-effects regression on the initial and final time points (Bates et al., 2015), which allows analysis of complex experimental designs and accounts for technical effects. We evaluated MEMcrispR with MAGeCK (model-based analysis of genome-wide CRISPR-Cas9 knockout) (Li et al., 2014) with data from two previously published CRISPR screens and found that MEMcrispR had comparable specificity but increased sensitivity to MAGeCK in analyses of negative and positive selection screens (Figures S1D and S1E). We therefore employed MEMcrispR for the analyses of our 3D and 2D screens (Figures 1A and 1B).

The 3D Screens Identified Major Developmental Pathways Involved in Anchorage-Independent Growth

In the 3D screen, we identified six overrepresented genes in TP53 proficient and 33 overrepresented genes in TP53-deficient background (false discovery rate [FDR] = 20%, fold change >2), consistent with the notion that TP53 loss is an important early step during tumorigenic growth (Figures 1B and S2; Tables 1 and S1). In both backgrounds, inactivation of NF2 and PTPN14, which are components of the Hippo tumor-suppressor pathway, led to increased growth. We identified additional genes driving 3D growth exclusively in the TP53 deficient background. In particular, consistent with the existing literature, we found TSC1 and TSC2, components of the growth-promoting mTOR signaling pathway, to act epistatically with TP53 (Akeno et al., 2015, Hasty et al., 2013). Lastly, we identified additional known tumor suppressors in the RPE 3D screen such as PTEN (also involved in the mTOR pathway), SAV1 (a member of the Hippo pathway), and yet-uncharacterized potential tumor suppressor genes such as AHR and FRYL. We verified the anchorage-independent growth potential of all aforementioned genes not only in RPE cells but also in MCF10A cells (Figures S3A–S3D).

Table 1

Genome-wide Screen Hits

Rank	RPEWT 2D	Rank	RPETP53−/− 2D	Rank	RPEWT 3D	Rank	RPETP53−/− 3D
1	NF2	1	NF2	1	NF2	1	NF2
2	TAOK1	2	PTPN14	2	PTPN14	3	PTEN
3	AMOTL2	3	KCTD10	5	AHR	4	PTPN14
4	LRRC24	4	TAOK1	6	EP300	6	TSC2
9	TP53	5	CAND1	N/A	N/A	13	SERPINB2
11	PDCD10	6	CUL3	N/A	N/A	16	TMEM184A
17	FRYL	7	AHR	N/A	N/A	26	THRAP3
19	PTPN14	8	UBE2M	N/A	N/A	27	TCEAL7
25	CDKN2B	9	FRYL	N/A	N/A	28	FRYL
32	OR8J1	13	ARNT	N/A	N/A	32	NAPB
33	KEAP1	14	PTEN	N/A	N/A	34	GOLGA8B
36	C5orf38	18	AMOTL2	N/A	N/A	35	ABCA9
71	AHR	19	SAV1	N/A	N/A	41	C6orf222
N/A	N/A	20	KEAP1	N/A	N/A	45	PRAMEF8
N/A	N/A	21	KCTD5	N/A	N/A	46	RSPH10B2
N/A	N/A	25	ACTG1	N/A	N/A	51	PDCD7
N/A	N/A	28	PRKACA	N/A	N/A	53	AHR
N/A	N/A	29	ARIH1	N/A	N/A	55	DDO
N/A	N/A	32	LATS1	N/A	N/A	56	CAP1
N/A	N/A	38	UBE2F	N/A	N/A	61	PVALB
N/A	N/A	39	RNF7	N/A	N/A	75	SYNGR2
N/A	N/A	44	FRMD6	N/A	N/A	78	TEX19
N/A	N/A	45	TSC1	N/A	N/A	93	RAD52
N/A	N/A	47	RNF19A	N/A	N/A	98	SLC9C2
N/A	N/A	48	CUL5	N/A	N/A	N/A	N/A
N/A	N/A	49	DCUN1D3	N/A	N/A	N/A	N/A
N/A	N/A	58	RNF146	N/A	N/A	N/A	N/A
N/A	N/A	60	TSC2	N/A	N/A	N/A	N/A
N/A	N/A	84	TNFAIP1	N/A	N/A	N/A	N/A
N/A	N/A	85	ARRB2	N/A	N/A	N/A	N/A
N/A	N/A	88	KLHL21	N/A	N/A	N/A	N/A
N/A	N/A	95	SOCS6	N/A	N/A	N/A	N/A

Top hits from 2D and 3D screens of RPE WT and TP53 knockout cell lines. Underlining indicates ubiquitination/neddylation-associated genes.

In the 2D proliferation screens, we identified 21 overrepresented genes in a TP53-proficient background and 132 overrepresented genes in a TP53-deficient background (FDR = 5%, fold change >2; Figure 1B; Tables 1 and S1). Using a curated tumor suppressor gene database (https://bioinfo.uth.edu/TSGene/) and two studies focusing on the identification of tumor suppressors (Davoli et al., 2013, Vogelstein et al., 2013), we identified more predicted tumor suppressors in the TP53−/− background in both 2D and 3D screens (Figure 1C). In addition, the 2D screen identified most of the tumor suppressor genes from the 3D screen, suggesting that the proliferation screen can provide valuable information on genes promoting tumorigenic transformation in vitro (Figure 1D).

Gene Ontology (GO) terms showed underrepresented genes from both TP53 proficient and deficient proliferation screens were enriched in essential pathways, such as ribosome biogenesis, transcription, cell cycle (Figure S3E). GO analysis in the overrepresented genes revealed distinct sets for each screen: the hippo pathway was enriched in all screens; the mTOR pathway in the TP53−/− screens; tissue growth and cell proliferation in the 3D screens; ubiquitination and neddylation in the TP53−/− 2D screen (Figure 1E).

Genetic Interaction between TP53 and Genes of the Neddylation Pathway

To identify knockouts that increase the rate of proliferation specifically in TP53 deficient cells, we developed two models to estimate global and cell line specific effects, respectively (Figure S1). The combination of these models together with a pairwise comparison of two cell lines allowed us to quantify gene knockouts that affect specifically one cell line (Figure 2A; Table S2). Interestingly, ubiquitination/neddylation genes were enriched specifically in the 2D RPE cells (19% of the top 100; Figure 2A; Tables 1 and S1), suggesting a genetic interaction between TP53 and neddylation pathway genes such as E2 conjugating enzymes (UBE2M and UBE2F), E3 cullin ring ligases (CUL3 and CUL5), and CUL3-specific adapters (KEAP1 and KCTD10) (Figure 2B). Interestingly, these genes were not identified in the TP53−/− 3D screen, indicating that although they have a significant impact on proliferation, inactivation of these genes alone is insufficient to promote anchorage-independent growth.

Figure 2

CUL3 Loss in TP53-Deficient Cells Results in Growth Defects and Persistent Proliferation

(A) Volcano plot depicting the enrichment of ubiquitination/neddylation pathway specific genes (highlighted in green) in TP53−/− cells.

(B) Scheme of ubiquitination/neddylation pathway depicting the fold change of each gene knockout in the 2D RPE screen.

(C) Bright-field images of three independent RPE cell lines.

(D) Histograms of Annexin V staining after transfection of independent siRNAs targeting CUL3 as well as one nontargeting (siSCR) siRNA into RPE and RPE cells. Mean intensity values of Annexin V derived from histograms as well as the fold changes compared to scrambled siRNA controls are indicated.

(E) Boxplots representing competitive growth assay of RPE cells tagged with mCherry mixed 1:1 with indicated knockout cells that are each tagged with GFP. Higher ratios indicate faster growing cells as compared to the RPE mCherry-tagged cells.

(F) Boxplot depicting exponential growth rates of RPE and RPE cells are seeded in high (90% confluency) or low (30% confluency) density onto plates and every 3 days cells were counted and reseeded. Growth rates are presented as ratios of growth rate in high versus low cell density.

p value is based on Welch two-sample t test, and the standardized effect size (sES) is depicted. In box-and-whiskers plots, boxes show the upper and lower quartiles (25%–75%) with a line at the median, whiskers extend from the 5 to the 95 percentile, and dots correspond to the outliers.

See also Figure S3 and Table S2.

As the neddylation pathway was highly enriched in the 2D screen in TP53−/− cells, we generated individual knockout clones of different components of the neddylation pathway, namely CAND1, UBE2M, UBE2F, CUL5, and CUL3, as well as the CUL3-specific adaptors KCTD10 and KEAP1. We observed specifically in RPE cells a change in cell morphology from epithelioid to mesenchymal/spindle shaped, a known characteristic of EMT (Figure 2C). Consistent with the literature, which has highlighted CUL3 as an essential gene (Kossatz et al., 2010, Singer et al., 1999, Tateishi et al., 2001, Zhou et al., 2013), we failed to generate knockout clones in RPE cells (data not shown) and observed an increase apoptotic signaling by Annexin V staining specifically in RPE cells that are transfected with small interfering RNAs (siRNAs) targeting CUL3 (Figure 2D). Next, we performed a reanalysis of damaging somatic CUL3 mutations using the The Cancer Genome Atlas (TCGA) database. Reassuringly, damaging CUL3 mutations showed frequent co-occurrence with inactivating TP53 mutations (p = 0.008, OR = 2.39; Fisher’s exact test) in the TCGA database, further corroborating the findings from our isogenic cell line-based screens.

We then asked whether CUL3 loss increased the rate of proliferation or led to persistent proliferation due to loss of contact inhibition in RPE cells and examined how CUL3 affected proliferation in the context of different cell densities. When cells were seeded in low densities, knockout of CUL3 had no additional effect on cell proliferation as opposed to the NF2 knockout (Figure 2E). However, the CUL3 knockout allowed cells to proliferate even when they are seeded in very high densities (p = 0.04, Welch two-sample t test), supporting our earlier observations that the effects of CUL3−/− on cellular proliferation are different than that of typical tumor suppressor genes such as NF2 and that increased proliferation occurs primarily due to loss of contact inhibition (Figures 2F and S3F).

CUL3 Loss in TP53-Deficient Cells Induces an Altered Transcriptional Program Leading to Persistent Proliferation

To gain further insights into how CUL3 loss induces persistent proliferation, we analyzed the proteome of three independent RPE cell lines. We identified 758 upregulated and 747 downregulated proteins with an FDR of 1% (Figure S4A; Table S3). Pathway enrichment analysis of upregulated proteins showed enrichment for cell adhesion, regulation of apoptosis, cell division, and cullin ring ligase complexes, whereas downregulated proteins were enriched for anchoring junction, cell junction, and extracellular matrix regulators (Figure S4B; Table S3). In addition, we analyzed the proteome of RPE cells to find potential overlaps between CUL3 and CAND1 knockouts. CAND1 is a cullin-binding protein that associates with unneddylated cullins (Chua et al., 2011, Min et al., 2003) and is required for the regulation of the cullin ring ubiquitin ligase (CRL) assembly and adaptor recycling (Goldenberg et al., 2004, Pierce et al., 2013, Wu et al., 2013, Zheng et al., 2002). Comparison between the deregulated proteins in the CAND1−/− and CUL3−/− samples revealed a strong correlation (R2 = 0.8), with the adaptor proteins of CUL3 representing a notable exception (Figure 3A). Consistent with the literature suggesting that CRLs are responsible for degradation of their own associated adapters, we identified eight adaptor proteins that are stabilized in CUL3−/− cells, but not in CAND1−/− cells (de Bie and Ciechanover, 2011, Emanuele et al., 2011). Interestingly, an established adaptor of CUL3, KEAP1, and its substrate NRF2, whose protein levels are known to be controlled by CUL3 in response to oxidative stress, were not affected by CUL3 loss (Figure S4C), possibly due to the lack of additional oxidative stress.

Figure 3

CUL3 Loss in TP53-Deficient Cells Leads to an Altered Transcriptional Program Associated with Proliferation and EMT

(A) Proteome of RPE and RPE cell lines compared to RPE cells. Log fold change of protein abundance of each cell line compared to RPE control cell line is plotted. Colors represent specified adjusted p value cutoffs. Increased abundance of CUL3-specific adapters as well as oncogenic factors are highlighted.

(B) Immunoblotting of components of the AP-1 (JUN, FOSL1, and FOSL2) complex as well as PD-L1 in RPE and RPE cells.

(C) Immunoblotting of components of the NF-κB and TGF-β pathways in RPE and RPE cells.

(D) iRegulon analysis predicts that 60% of the overexpressed genes are transcriptional targets of NF-κB and 13% are targets of AP-1 (1,269 genes filtered by adjusted p < 0.05 and log2 fold change of >1).

(E) Gene set enrichment analysis of Gene Ontology (GO) cell proliferation, NF-κΒ, and EMT pathways based on mRNA-seq data.

(F) Network of pathways implemented in EMT. Colors mark the gene expression fold changes observed in RPE cells compared to the RPE control cell line.

See also Figures S4 and S5 and Tables S3 and S4.

Apart from the adaptor proteins that are stabilized upon RPE, we identified the levels of tumorigenic transcription factors JUN, FOSL1, and FOSL2 to be increased in RPE cells (Figure 3B). In addition, we observed PD-L1 (CD274) levels to be elevated in the RPE cell lines (Figure 3B). PD-L1 is a ligand that binds to the co-inhibitory molecule programmed cell death 1 (PD-1) expressed on the surface of activated T cells (Fife and Bluestone, 2008). Interaction between PD-1 and PD-L1 induces suppression of T cell proliferation and is essential in regulating immune tolerance (Buchbinder and Desai, 2016). The connection between CUL3 and PD-L1 was also identified in a haploid genetic screen in HAP1 cells, where PD-L1 was shown to be regulated by CUL3 adaptor protein SPOP that binds to and is regulated by CDK4 (Chen and Chen, 2016, Mezzadra et al., 2017). Our results further support that CUL3 regulates PD-L1 in RPE cells, indicating that PD-L1 regulation by CUL3 may occur in different cellular contexts.

Since PD-L1 expression has been shown to increase in response to NF-κB and TGF-β signaling in cells (Asgarova et al., 2018), we also evaluated whether NF-κB signaling is upregulated in RPE cells. NF-κB has a critical role for several cellular processes such as proliferation, differentiation, immunity, and inflammation. Cytosolic, inactive NFκB can be activated by phosphorylation of its inhibitor, IκB, by the IκB kinase complex (IKK) which then releases it to translocate to the nucleus (Oeckinghaus and Ghosh, 2009). IKKβ, which partakes in the IKK complex, was shown to be targeted by KEAP1 (Lee et al., 2009) and interacts also with KLHL21 (Mei et al., 2016). Both proteins are adapters of CUL3 and were identified in the 2D screen in the RPE cell line, providing further evidence of the importance of the neddylation pathway leading to NF-κB activation. We observed in CUL3−/− cells that the inhibitor of NF-κB, IκBα, is phosphorylated, leading to its degradation and the subsequent activation of the NF-κB signaling pathway (Figure 3C).

Since NF-κΒ signaling activates a variety of cellular functions via transcription, we tested whether these effects can be detected in RPE cells. To this end, we analyzed gene expression changes between RPE and RPE cells by mRNA sequencing (mRNA-seq) and identified 1,630 upregulated and 1,558 downregulated genes with an FDR of 1% (Figure S4D; Table S4). Comparison of the mRNA-seq and proteomics data revealed similarities between the up- and downregulated proteins (Figure S4E) with moderate correlation (R2 = 0.43). To assess whether differentially expressed genes are associated with certain gene ontologies, we conducted GO and pathway enrichment analysis. We found that upon CUL3 loss, upregulated genes were enriched for regulation of cell proliferation, growth factor activity, and NF-κB signaling, whereas downregulated genes were enriched for extracellular matrix, cell adhesion, and cell junction processes (Figure S4F; Table S4). In addition, we inferred that NF-κB controlled 60% of the upregulated genes in the gene expression datasets based on iRegulon transcription factor enrichment analysis (Janky et al., 2014) (Figure 3D). Thus, major transcriptional changes occur in RPE cells, implicating NF-κB.

CUL3 Loss Associated with a Partial EMT Phenotype

One of the target genes of NF-κB signaling is TGFB1, an important factor for differentiation and epithelial-to-mesenchymal transition (EMT) (Xu et al., 2009). We thus sought to examine a possible EMT phenotype initiated by active TGF-β signaling (Bakiri et al., 2015, Huber et al., 2004). Gene set enrichment analysis showed that EMT-related genes are significantly altered in the RPE cell lines (Figure 3E). Analysis of the mRNA-seq data revealed activation of three signaling pathways downstream of NF-κΒ implicated in EMT: TGF-β, WNT, and NOTCH (Figure 3F). Based on mRNA-seq, we found that the TGF-β activator SMAD3 was significantly upregulated, whereas the inhibitory proteins SMAD6 and SMAD7 were downregulated. TGFBR1 receptor was overexpressed, as was the TGFB2 ligand, suggesting activation of the TGF-β pathway. In addition, another activator of the TGF-β pathway, activin A (INHBA), which is implicated in EMT-promoting tumorigenicity, was highly expressed in RPE cells (Bashir et al., 2015) (Figure 3F). We additionally found several members and the targets of the EMT-associated WNT and NOTCH signaling pathways to be significantly overexpressed in the RPE cells (Figures 3F and S5A). Lastly, we found the core transcription factors SNAI2 (SLUG), Twist2, ZEB2, ETS1, and ETS, which are implicated in EMT (Yang and Weinberg, 2008, Zeisberg and Neilson, 2009), to be overexpressed in RPE cells (Figure 3F).

Together, our data present multiple layers of evidence resembling an EMT signature. However, the cells have not fully transitioned to a mesenchymal state, as we did not find altered expression of structural EMT markers such as vimentin, fibronectin, E-cadherin, and VE-cadherin. Therefore, we hypothesize that CUL3 initiates a transcriptional program that leads to the induction of the EMT transcriptional program transitioning the cells from a normal state, without commitment, to a complete transition to mesenchymal cells. These results are in line with recent studies demonstrating that EMT is not a binary process but rather a multistep process with several intermediate states (Pastushenko et al., 2018). We compared our mRNA-seq data to the list of genes that are identified to be responsible for a partial EMT (pEMT) state based on single-cell RNA sequencing of developing head and neck squamous cell carcinoma samples (Puram et al., 2017). Doing so, we observed a strong overlap between the genes that are overexpressed in RPE cells and pEMT genes; 59 out of 91 pEMT genes were found to be overexpressed in our dataset (p = 0.004, two-sided Fisher’s test), indicating that our RPE cells represent tumors that display a pEMT signature.

RPE Have Increased DNA Damage and Genome Instability

Induction of EMT program has been shown to induce persistent proliferation leading to genomic instability via suppression of nuclear lamins (Comaills et al., 2016). Given that RPE cells can proliferate persistently and show characteristics of pEMT, we next tested whether CUL3 loss can also trigger genomic instability. First, we imaged RPE and RPE cells stained with actin and γ-H2AX. In line with previous reports suggesting impaired actin cytoskeleton structure (Chen et al., 2009), we observed overall more intense actin bundles in RPE cells (Figure 4A). Strikingly, we also observed very high levels of DNA fragmentation and increased levels of extensive γ-H2AX staining (p = 3.8 × 10−5, Mann-Whitney U test) (Figures 4A and 4B), hallmarks of genomic instability and high levels of DNA damage, respectively. This phenotype was specific to RPE cells and was not detected in RPE or RPE cells (Figure S5B). These cells consistently had higher numbers of micronuclei (p = 1.3 × 10−6, Mann-Whitney U test) and abnormally shaped/fragmented nuclei (p = 3.7 × 10−12, Mann-Whitney U test) (Figure 4C), which together pointed out a potential role of CUL3 in protecting genome integrity by preventing persistent proliferation.

Figure 4

CUL3 Loss Results in Increased Levels of DNA Damage and Genomic Instability

(A) Immunofluorescence staining of γ-H2aX and actin filaments in RPE and RPE cell lines. Inlets show exemplary images of nuclei observed in these cell lines.

(B) Quantification of cells in (A). Boxplots indicate the percentage of cells with extensive DNA damage assessed by increased γ-H2aX signal.

(C) Quantification of cells in (A). Boxplots indicate the percentage of cells with micronuclei (marked by a white arrow) and abnormal nuclei in RPE and RPE cell lines as exemplified in the inlets.

(D) Flow cytometry analysis of DNA content in RPE and RPE cell lines. Percentages of each cell-cycle stage are depicted above the histograms.

(E) Low-coverage whole-genome sequencing of independent RPE and RPE cell lines. Boxplots indicate the occurrence of arm or whole-chromosome losses/gains in RPE and RPE cell lines.

(F) Colony-formation assay in RPE and RPE cell lines upon inhibition of ATM, ATR, or DNA-PK kinases. Boxplots indicate the ratios of colonies formed after inhibitor treatment (normalized to DMSO control).

(G) Colony-formation assay in NCI-H358 and NCI-H1703 cells. Cas9-expressing NCI-H358 and NCI-H1703 cell lines were transfected either with control sgRNA or sgRNAs targeting CUL3 and treated with the ATM inhibitor AZD0156. Boxplots represent the area of the wells occupied upon each treatment (normalized to DMSO control transfected with a control sgRNA (sgCTRL)).

p value is based on Welch two sample t test and the standardized effect size (sES) are depicted; Box-and-whiskers plots: boxes show the upper and lower quartiles (25%–75%) with a line at the median, whiskers extend from the 5th to the 95th percentile and dots correspond to the outliers.

See also Figure S5.

To identify the specific adaptors that may work with CUL3 in this process, we performed a targeted arrayed CRISPR-Cas9 screen (Serçin et al., 2019). We tested individual effects of 11 adapters identified to be stabilized by CUL3 based on the proteomics screens (Figure 3A). We quantified the level of abnormally shaped/fragmented nuclei upon each adaptor knockout and identified the most significant increase upon sgBTBD2 and sgKCTD10 single guide RNA (sgRNA) transfection along with sgKLHL9 and sgKLHL21 (Figure S5C). KCTD10 and KLHL21 were also identified to increase proliferation in our 2D screens (22- and 2.2-fold, respectively), and led to a significant reduction of nuclear lamins along with CUL3 (Figure S5D). Thus, these adapters may function together with CUL3 in promoting persistent proliferation and genomic instability.

Typically, abnormally shaped nuclei are a result of mitotic errors and lead to aneuploidy. To test for the occurrence of changes in ploidy, we analyzed RPE clones by flow cytometry (Figure 4D). In all cases, we observed increase of ploidy in RPE cells, reaching up to 49% of cells with more than 4N DNA content. Increased ploidy and DNA fragmentation suggested that these cells carry more chromosomal abnormalities. To address this, we generated 20 independent RPE clones and performed low-coverage whole-genome sequencing on these cells to quantify large copy-number alterations that are likely stemmed from this heightened level of genomic instability. Indeed, RPE cells led to a significant increase in large copy-number alterations as compared to 12 independent RPE clones (p = 4 × 10−4, Mann-Whitney U test), suggesting an important role of CUL3 in maintaining genome stability in TP53 deficient cells (Figure 4E).

Finally, we tested whether increased levels of genomic instability can create a vulnerability for RPE cells. To this end, we treated the RPE and RPE isogenic lines with inhibitors of the major kinases of the DNA damage response (ATM, DNA-PK, and ATR) (Blackford and Jackson, 2017). We tested whether CUL3 knockout sensitizes RPE−/− cells to these inhibitors. Measuring the ability of these cells to form colonies in response to inhibitor treatments, we identified a damaging effect of ATM inhibitor (ATMi) (p = 0.027, Welch two sample t test), particularly on RPE cells (Figure 4F). While the DNA-PK inhibitors show similar trends, we did not observe any significant difference on colony-formation capacity of CUL3 cells upon ATR inhibition. To search for additional evidence for this potential synthetic lethal interaction between CUL3 deficiency and ATMi, we constructed two additional Cas9-expressing lung tumor cell lines (NCI-H358 and NCI-H1703) that are deficient in TP53. We observed a dramatic loss of viability in cells transfected with CUL3 sgRNAs that are treated with ATMi. This effect was more than an additive effect and was particularly pronounced in NCI-H358 cells. These results are in line with our observations in RPE cells and together suggest that the CUL3-deficient tumors can be targeted via ATMi (Figures 4G, S5E, and S5F).

Discussion

In this study, we conducted genome-wide knockout screens in non-tumorigenic RPE cells with and without functional TP53 to identify genes that promote oncogenic proliferation. We uncovered several previously studied tumor-suppressor genes, such as NF2 and PTPN14, and also identified factors that promote proliferation in a TP53-dependent manner, such as TSC1 and TSC2. The majority of genes identified in our screens enhanced tumorigenic growth specifically in RPE cells, which highlights the importance of TP53 deficiency as a first hit in tumorigenesis.

We focused our study on the role of a cullin ring ligase, CUL3. CUL3 and its adaptor, KEAP1, are frequently mutated in lung cancers (Kandoth et al., 2013, Singh et al., 2006). CUL3/KEAP1 target both NRF2 (Kobayashi et al., 2004), a master regulator in oxidative stress responses, and IKK-β (Lee et al., 2009), an activator of the NFκB pathway, for degradation. Several prior studies demonstrated that CUL3 loss of function is lethal in normal cells (Kossatz et al., 2010, Singer et al., 1999, Tateishi et al., 2001, Zhou et al., 2013). Our data show that in a TP53-deficient background, CUL3 inactivation not only is viable but also leads to persistent proliferation. Reassuringly, our reanalysis of damaging CUL3 mutations in a pan-cancer setting demonstrated co-occurrence with inactivating TP53 mutations, further corroborating these findings.

How does CUL3 promote persistent proliferation? Profiling isogenic RPE and RPE cells, we found that most of the changes in RPE cells occurred at the gene expression level. In RPE cells, we identified an oncogenic transcriptional program leading to heightened growth signaling and EMT in RPE cells. NF-κB and AP-1 transcription factors were likely responsible for the majority of changes of the transcriptome in response to loss of CUL3, leading to an increase in downstream targets such as TGF-β and pEMT. EMT is considered a gradual process with many intermediates (Pastushenko et al., 2018). Even though cells that have progressed through full EMT have higher metastasis rates (Pastushenko et al., 2018, Thompson et al., 1994, Yang and Weinberg, 2008), EMT intermediates have also metastatic potential (Pastushenko et al., 2018). Here, we propose that loss of CUL3 initiates a transcriptional program that promotes EMT without a complete transition to mesenchymal cells.

One of the outcomes of an EMT in cancers is an increase in genomic instability (Comaills et al., 2016). We also observed increased levels of genomic instability in the RPE clones, as demonstrated by increased rates of DNA fragmentation, aneuploidy, and large structural variations. These events can be caused by mitotic defects resulting in genomic instability and consequently somatic copy-number abnormalities that can dramatically alter the genomic landscape (Yates and Campbell, 2012). Combination of an increase in proliferation and in genomic instability may at first appear paradoxical, since genomic instability is presumably associated with a decrease in mitotic rates required to repair insults to the genome. TP53-deficient cells lose checkpoint control mechanisms and likewise show propensity to be genomically instable (Eischen, 2016, Liu et al., 2017, Yeo et al., 2016). In conjunction with CUL3 loss, perhaps owing to downstream effects of TGF-β signaling, we observed both persistent proliferation phenotype and genome instability, providing a potentially very fertile ground for tumor evolution and growth.

Can tumors with CUL3 deficiency be specifically targeted? In recent years, targeting major DNA damage signaling kinases has been realized as a favorable option, since cancer cells with increased genomic instability tend to show a greater dependency on DNA damage response pathways to deal with endogenous or exogenous DNA damage (O’Connor, 2015, Velic et al., 2015). To investigate this possibility in the context of CUL3 deficiency, which was also implicated in regulating DNA end resection (Ferretti et al., 2016), we employed inhibitors against the three major kinases in the DNA damage response initiating signaling cascades upon DNA damage. We found that RPE cells show an increased susceptibility specifically to ATM inhibitors while they were mostly unaffected by the inhibition of ATR, suggesting an increased sensitivity toward double-strand break signaling. We also observed similar effects in lung tumor cell lines with TP53 and CUL3 deficiency. Although we tested the effect of DNA damage response inhibition specifically on TP53−/−,CUL3−/− cells, it is possible that the clinical significance of these observations may extend to a wider range of tumors with heightened TGF-β signaling. Considering the frequency of upregulated TGF-β signaling in different types of cancer (Haque and Morris, 2017), our findings could offer more widely applicable therapeutic strategies for tumor treatment.

STAR★Methods

Key Resources Table

Lead Contact and Materials Availability

Further information and requests for resources and reagents should be directed to and will be fulfilled by the Lead Contact, Dr. Balca R.Mardin (mardin@bio.mx).

Plasmids and cell lines generated in this study are available upon request by the Lead Contact, Dr. Balca R.Mardin (mardin@bio.mx).

Experimental Model and Subject Details

Cell lines

hTERT RPE-1 cells (46,XX) were purchased from ATCC® (CRL-4000). RPE cells were generated in this study with sgRNAs targeting the first exons of the TP53 gene (sequences available in Key Resources Table). The RPE cell lines were grown in DMEM/F-12 (Thermo Fischer Scientific, 11320074) medium supplemented with 10% FBS (Thermo Fischer Scientific, 10500064) and Antibiotic-Antimycotic (Thermo Fischer Scientific, 15240062). MCF10A cells (46,XX) were purchased from ATCC® (CRL-10317). MCF10A cells were generated in this study. MCF10A’s growth medium contained DMEM/F-12 supplemented with 5% Horse serum (Thermo Fischer Scientific, 26050088), 20ng/ml Epidermal growth factor (Sigma-Aldrich, E9644), 500 μg/ml Hydrocortisone (Sigma-Aldrich, H0888), 100ng/ml Cholera toxin (Sigma-Aldrich, C8052), 10 μg/ml Insulin (Thermo Fischer Scientific, 12585014), and 1X Antibiotic-Antimycotic. 293T cells were grown in DMEM (Thermo Fischer Scientific, 31966021) medium supplemented with 10% FBS and 1X Antibiotic-Antimycotic. NCI-H358 and NCI-H1703 lung cancer cell lines (ATCC CRL-5807 and CRL-5889) were maintained in RPMI-ATCC (Thermo Fischer Scientific, CC) medium supplemented with 10% FBS and 1X Antibiotic-Antimycotic. To generate Cas9 nuclease-expressing cells, the Edit-R inducible lentiviral particles (Horizon Dharmacon) was used according to the manufacturer’s protocol. The transduced cells were selected in the presence of 1 μg/mL blasticidin for seven days. The expression of Cas9 was controlled by a doxycycline-inducible promoter and the expression was induced with 1 μg/mL doxycycline. All cell lines were maintained in appropriate densities and were incubated in a humidity-controlled environment (37°C, 5% CO2).

All cells lines tested negative for mycoplasma contamination.

Method Details

3D and 2D CRISPR screens

Lentiviral libraries were produced by co-transfecting 293FT cells with 100 μg of GeCKOv2 library (Addgene, 1000000048), 50 μg of pMD2.G (Addgene, 12259), and 75 μg of psPAX2 (Addgene, 12260) using Lipofectamine 3000 (Thermo Fisher Scientific, L3000001) following manufacturer’s specifications. After 60 hours, viral supernatant was collected, filtered through a 0.45 μm low protein binding membrane Steriflip HV/PVDF (Millipore), and concentrated 100-300x through ultracentrifugation at 24,000rpm, for two hours, at 4°C (Beckman Coulter). Viral pellet was resuspended in 8 mL of DMEM supplemented with 1% BSA overnight, at 4°C with agitation.

100 million RPE and RPE cells per library, aliquoted to one million cells per well of a 6-well plate, were spinfected with the lentiviral library at 2,000rpm, for two hours at 37°C, with the addition of 8 μg/mL polybrene (Sigma-Aldrich, TR-1003) to increase infection potential. Cells were left to recover for 24 hours incubated at 37°C, 5% CO2. The infected cells were selected with 7.5mg/μl puromycin (Thermo Fisher Scientific, A1113803) and puromycin treatment was maintained throughout the screening process.

Before 2D or 3D selection, one third of the cells were collected for obtaining the baselines representation of the sgRNAs.

For the 3D screen, the cells were seeded into 1:1 mixture of 2x DMEM/F-12 media (Thermo-Fischer Scientific) and 0.7% agarose in CellSTACK flasks of 636 cm2 growth area with low attachment surface (Corning, CLS3303-8EA). Cells were incubated for five weeks in a humidified atmosphere. After 3D selection, cells that formed spheres were collected and incubated for seven days in 150 mm × 25 mm dishes at 37°C, 5% CO2.

For the 2D screen, the cells were transferred to eight Falcon™ 5-layer flasks (Corning, 353144) with media supplemented with 7.5mg/μl puromycin. Cells were incubated in a humidified atmosphere to proliferate for two weeks and subdivided every two to three days to maintain appropriate cell densities.

PCR Amplification of the GeCKO Libraries for Sequencing

Genomic DNA was isolated from the initial and final time point collected cells using QIAamp DNA Blood Maxi Kit (QIAGEN, 51194) following manufacturer’s protocol. DNA was amplified for 18 cycles using plasmid-specific primers for the amplification of the sgRNA coding sequences. An additional round of PCR amplification was run for 8 cycles to ensure the addition of Illumina adaptors. Herculase II Fusion DNA Polymerase (Agilent, 600677) was utilized for the PCR reactions according to manufacturer’s specifications. The primers used for the reactions are depicted in the Key Resources Table. The PCR reaction products were purified using 1.8x Agencourt AMPure XP (Beckman Coulter, A63880) and eluted in nuclease-free water (Thermo Fisher Scientific). Samples were run in a 1% agarose gel and a 341bp fragment was extracted and purified using NucleoSpin Gel and PCR Clean-up Kit (Macherey-Nagel, 740609) following manufacturer’s protocol. All samples were pooled and sequenced using Illumina HiSeq 2000 platform (Illumina), using 50 base pair single reads.

CRISPR (sgRNA coding sequence) design

For the establishment of individual knockout cell lines for each of the candidate genes, the most efficient CRISPR (sgRNA coding sequence) from the GeCKO library and two custom-made CRISPRs were used. The custom made CRISPRs were designed utilizing the web tool DeskGen (Desktop Genetics Ltd, 2017).

CRISPR cloning

The CRISPR sequences (sgRNAs – Key Resources Table) were cloned according to a modified protocol from Shalem et al. (2014) to lentiCRISPRv2-GFP (gift from David Feldser (Addgene plasmid #82416) or pLentiCRISPR-2G-mCherry (based on lentiCRISPRv2 with two sgRNAs in tandem). CRISPR pairs were annealed by adding 10 μM of each CRISPR oligo and 1X T4 Ligation Buffer (NEB, B0202S) at 37°C for 30 minutes, heated up to 100°C and ramped down by 1°C/second to room temperature. In parallel, the vector was digested with the BsmBI/Esp3I enzyme (Thermo Fischer Scientific, ER0451) according to manufacturer’s specifications. The ligation reaction of the CRISPR to the vector was set up using Quick Ligation Kit (NEB, #M2200) and incubated at room temperature for 10 minutes. The plasmid was transformed into Stbl3 chemically competent bacteria (Invitrogen, C737303) following manufacturer’s protocol. The bacteria were grown overnight on antibiotic (ampicillin) containing growth media at 37°C. Bacterial colonies were picked and grown in the appropriate liquid culture overnight at 37°C. Plasmid DNA was extracted using QIAprep Spin Miniprep Kit (QIAGEN, #27104) and the purified plasmids were verified by Sanger Sequencing with the primers CAGGCCCGAAGGAATAGAAG (forward) and CAACTTCTCGGGGACTGTGG (reverse).

Virus Production and Establishment of knockout cell lines

The pLentiCRISPR-2G-GFP and pLentiCRISPR-2G-mCherry vectors containing the selected CRISPR sequences were mixed 1:1:1 with 0.75 μg psPAX2 (gift from Didier Trono, Addgene plasmid #12260) and 0.5 μg pMD2.G (gift from Didier Trono, Addgene plasmid #12259) lentiviral plasmids and transfected to 80% confluent 293FT cells using Lipofectamine 3000 kit (Thermo Fischer Scientific, L3000008) according to manufacturer’s protocol. After 72 hours of incubation at 37°C, 5% CO2, viral supernatant was collected.

Viral supernatant was added to an amount corresponding to 20% multiplicity of infection to RPE or RPE cells grown to 70% confluence and supplemented with 8 μg/mL polybrene (Sigma-Aldrich). The cells were infected with Centrifugation for two hours at 2,000rpm. Cells were left to recover in a humidified atmosphere for 24 hours and growth media was replaced. Seven days later, infected cells presenting the fluorescent marks GFP (for candidate genes) or mCherry (for controls) were sorted in bulk using the BD FACSAriaFusion flow cytometer (BD Biosciences, San Jose)).

Sanger sequencing verification

DNA was extracted DNeasy Blood & Tissue Kit (QIAGEN, 69506) and Sanger sequencing was performed by GATC Biotech (Eurofins Genomics). Analysis of sequences was conducted by Indigo-GEAR (https://gear.embl.de/) to identify knockouts by deconvolution of the mixed signal on the chromogram.

Competitive growth assay

GFP-expressing knockout cells of each of the candidate genes and mCherry-expressing control cells were counted using TC20 Automated Cell Counter (Bio-Rad, manufacturer’s protocol) and equal amounts were mixed together. The 1:1 ratio of GFP:mCherry was verified at seeding using LSR-Fortessa flow cytometry instrument (BD Biosciences). Three different densities were considered for each cell line tested: low, where cells were seeded to 20% confluency, medium, where cells were seeded to 50% confluency, and high, were cells were seeded to 80% confluency.

The cells were grown for seven days at 37°C, 5% CO2. GFP and mCherry intensities within each of the mixed cell populations were then measured using the LSR-Fortessa instrument (BD Biosciences). The post-acquisition analysis was performed with FlowJo v10.2 software (FlowJo, LLC).

Growth characterization of CUL3 wild-type and knockout cell lines

RPE and RPE cells were seeded at 30% and 90% confluency on 12 well plate for each cell line. 72 hours later the cells were counted. The experiment was conducted in triplicates.

Soft agar verification assay

5000 cells were seeded in 96 well plates and resuspended in liquid 0.35% soft agarose in DMEM-F/12. Mixture was allowed to solidify during incubation at 37°C, 5% CO2. After seven days of culture, the cell viability was measured with a CellTiter-Glo® 3D Cell Viability Assay (Promega, G9681) following manufacturer’s instructions. The plates were measured in a fluorescent plate reader (Tecan infinite M1000Pro) and analyzed in R.

Cell cycle analysis

Cells were filtered through a Corning® cell strainer size 70 μm (Sigma-Aldrich, CLS431751) and fixed in 80% ethanol. The fixed cells were incubated at −20°C for 30 minutes, then centrifuged at 2,000rpm for five minutes and resuspended in PBS supplemented with 50 μg/ml Propidium iodide (PI) (Thermo Fisher Scientific, P3566) and 100 μg/ml RNase A (Thermo Fisher Scientific, EN0531). PI fluorescence intensity was measured using an LSR-Fortessa flow cytometer (BD Biosciences). The post-acquisition analysis was performed with FlowJo v10.2 software (FlowJo, LLC).

Immunoblotting

Total protein was extracted from cells using RIPA buffer (CST, 9806) and Benzonase® endonuclease (Merck-Millipore, 71206-3). Protein was quantified using Protein Assay Dye Reagent (Bio-Rad, 5000006). Samples were prepared for blotting using 4x Laemmli Sample Buffer (Bio-Rad, 161-0747). Proteins were separated by 4%–15% Mini-PROTEAN® TGX Precast Protein Gels (Bio-Rad, 4568084) and blotted to Trans-Blot® Turbo Nitrocellulose membranes (Bio-Rad, 1704158) using Trans-Blot® Turbo Transfer System (Bio-Rad, manufacturer’s specifications). Membranes were blocked in 10% low fat milk in TBS-T (50 mM Tris pH 7.5, 150 mM NaCL, 0.1% Tween). Primary antibodies were added to suggested concentrations in 5% low-fat milk in TBS-T and incubated overnight at 4°C. Membranes were washed with TBS-T and incubated with horseradish peroxidase-linked secondary antibodies diluted to suggested concentrations for one hour at room temperature. Membranes were washed with TBS-T, Clarity Western ECL Blotting Substrate (Bio-Rad, 1705060) was added and then they were photographed for chemiluminescence using ChemiDoc MP Imaging System (Bio-Rad).

Immunofluorescence and Microscopy

Cells were grown on coverslips and fixed with 4% paraformaldehyde solution. The coverslips were treated with 0.1% Triton X-100 for ten minutes for membrane permeabilization. Cells were washed with PBS and incubated with 5% BSA in PBS for 30 minutes. The primary antibody was diluted to recommended concentration in 3% BSA, added to the coverslips, and incubated for one hour in a humidified chamber. Coverslips were washed with PBS and secondary antibody diluted to recommended concentration in 3% BSA in PBS together with Hoechst dye were added and incubated for 30 minutes. The cells were washed with PBS and fixed with ProLong (Thermo Fischer Scientific, P36930). Cell imaging was conducted with a Nikon Ti-E using CFI P-Apo DM 60x Lambda oil/ 1.40/ 0,13 phase contrast PH3 objective, emission (single band): DAPI: 447/60, GFP: 525/50, Cy3: 593/40, mCherry: 641/75, Cy5: 692/40.

Mass Spectrometric Identification of Proteins

Snap frozen cell pellets of RPE, RPE, and RPE cell lines were resuspended in 50μl PBS supplemented with 1% SDS in 100 mM HEPES, pH 8.4 and protease inhibitor cocktail (Roche, 11873580001) and heated up to 95°C for five minutes and cooled down on ice. DNA was removal through Benzonase® endonuclease (Merck-Millipore, 71206-3) treatment at 37°C for one hour. Protein concentration was then determined and adjusted to 1μg/μl. 20μg of extracts were subjected to an in-solution trypsin digestion following a modified version of the Single-Pot Solid-Phase-enhanced Sample Preparation (SP3) protocol (Hughes et al., 2014, Moggridge et al., 2018). The lysates were added to Sera-Mag Beads (Thermo Scientific, 4515-2105-050250) in 15% formic acid and ethanol and shaken for 15 minutes at room temperature. SDS was removed through four 70% ethanol washes. Proteins were digested overnight at room temperature using sequencing grade modified trypsin (Promega, V5111) in HEPES/NaOH, pH 8.4 in the presence of 1.25mM TCEP and 5mM chloroacetamide (Sigma-Aldrich, C0267). Sera-Mag Beads were separated and washed with 2% DMSO. Combined eluates were dried down, reconstituted in H2O, and reacted with 80μg of TMT10plex (Thermo Scientific, 90111) (Werner et al., 2014) in acetonitrile for one hour at room temperature. Excess TMT reagent was quenched with the addition of 5% hydroxylamine (Sigma, 438227). The peptides were mixed 1:1 across all TMT-channels and subjected to a reverse phase clean-up step (OASIS HLB 96-well μElution Plate, Waters 186001828BA) followed by an offline fractionation under high pH condition (Hughes et al., 2014). Twelve fractions were obtained and analyzed by LC-MS/MS on a Q Exactive Plus (Thermo Fisher Scientific) following the protocol adapted from Becher et al. (2018).

Analysis and normalization of the acquired data were performed using IsobarQuant (Franken et al., 2015) and Mascot V2.4 (Matrix Science) as previously described (Becher et al., 2018). At least 2 unique peptides with a peptide length of minimum seven amino acids and a false discovery rate lower than 0.01 were required on the peptide and protein level (Savitski et al., 2015). The output file from IsobarQuant was processed with the R programming language. Limma was used for quantification, batch cleaning, variance stabilization normalization, and differential expression. Proteins with a false discovery rate below 5% and a fold change of minimum 50% were classified as “hits.”

mRNA-Seq

Total RNA was extracted from the cells using RNeasy MinElute Cleanup kit (QIAGEN, 74204). RNA quality control was performed using the 2100 Bioanalyzer platform (Agilent Technologies). Library preparation was performed with 200ng starting material using TruSeq Stranded mRNA HT (Illumina, 20020594) with a Beckman Biomek FX automated liquid handling system (Beckman Coulter). Library quantification and quality control were carried out on the Advanced Analytics Technologies Fragment Analyzer. The libraries were sequenced on an Illumina HiSeq 2500 platform (Illumina), using 50 base pair single reads for low-pass sequencing.

Low coverage whole genome sequencing

Genomic DNA was extracted from RPE clones using QIAamp DNA Blood Maxi Kit (QIAGEN, 51194) following manufacturer’s instructions. Library preparation for sequencing was performed with 300ng starting material using NEBNext® Ultra II DNA Library Prep Kit for Illumina (New England Biolabs, E7645) with a Biomek i7 automated liquid handling system (Beckman Coulter). Custom six base pair barcodes were designed to enable pooling. Quality control and quantification of the libraries were performed using Agilent 2100 Bioanalyzer (Agilent Technologies). DNA sequencing was performed on an Illumina HiSeq 2000 platform (Illumina), using 50 base pair single reads for low-pass sequencing.

Clonogenic Assay

Cells were trypsinized and counted. For RPE-1 cell line 500, 1000 and 3000 cells, for NCI-H358 cell line 6000 and 9000 cells and for NCI-H1703 cell line 500 and 2000 cells were seeded on 6 well plates in duplicates. The next day, the cells were subjected to treatment with M3814 (MedKoo 206478), VE-822 (MedKoo 406258), and AZD0156 (Selleckchem, S8375) to inhibit DNA-PK, ATR, and ATM, respectively. Inhibitors were used at the following concentrations: 200 nM/ 500 nM (M3814), 50 nM /500 nM (VE-822) and 200-500-1000 ng/μl (AZD0156). 24 hours later, the treatment containing media was replaced with drug-free DMEM/F12 supplemented with 10% FBS and 1% Antibiotic-Antimycotic. The cells were incubated for ten to thirteen days to form colonies in a humidified atmosphere (37°C, 5% CO2). Subsequently, the colonies were washed with PBS, stained with 0.15% crystal violet (Sigma-Aldrich, V5265) for five minutes at room temperature and counted with GelCount™ (Oxford Optronix). Negative control of 250 untreated cells of each cell line were included in duplicates.

Silencing of the CUL3 gene

RPE and RPE cells were seeded to 40% confluency. 24 hours later, cells were transfected with two pre-designed siRNA targeting CUL3 gene (Thermo Fisher Scientific, s16050 and s16049) using Lipofectamine RNAiMAX (Thermo Fisher Scientific, 13778150) following the manufacturer’s instructions. Cells treated with Scrambled siRNA (Thermo Fisher Scientific, 4390843) served as a negative control. 72 hours post transfection, cells were collected and incubated with Annexin V Alexa Fluor 488 conjugate (Thermo Fisher Scientific, A13201) and Sytox blue dead cell stain (Thermo Fisher Scientific, S34857) according to manufacturer’s indications. The stained cells were analyzed by flow cytometry (FACS Fortessa). Cells treated with 10μM Camptothecin (S1288, Selleckchem) for 4 hours prior the analysis were used as a positive control. Data was analyzed with FlowJo V10.

RNA Extraction, cDNA Synthesis, and qPCR

Total RNA isolation was performed from 106 cells using the RNeasy Plus Mini kit (QIAGEN, Düsseldorf, Germany) as described in the manufacturer’s protocol. The RNA samples were diluted to 250 ng/μl final concentration. All RNA samples within an experiment were reverse transcribed at the same time with the qScript cDNA SuperMix (Quanta Biosciences) using 500 ng of RNA as a template and stored in aliquots at −80 oC. Real-time PCR with Fast SYBR® Green (ThermoFisher) detection was performed using a QuantStudio 5 Real-Time PCR system (Applied Biosystems). The relative quantification of each sample was performed using the comparative Ct method. The acidic ribosomal phosphoprotein P0 gene (36B4) is used as a housekeeping gene. To compare the transcript levels between different samples the 2 -ΔCt method was used.

Flow Cytometry

All experiments involving Flow Cytometry were performed on the same instruments with the following specifications: LSR-Fortessa (Becton Dickinson) 355(20mW)/405(50mW)/488(60mW)/561(75mW)/640(40mW) and FACSAriaFusion (Becton Dickinson) 405(100mW)/488(100mW)/561(75mW)/650(100mW). Instruments were tested daily with CST and resulting settings were modified for the highly autofluorescent RPE cells to ensure optimal resolution of dim and high reporter signals. DAPI was measured with 355nm or 405nm excitation using a 450/50 band pass; GFP was measured with 488nm excitation and 530/30 band pass; mCherry or PI (cell cycles) was measured with 561n excitation and 610/20 band pass. In all cases, singlet cells were gated by a combination of FSC-A/H and SSC-A/H when analyzing live cells; live cells were excluded by DAPI exclusion (1 μg/ml). When measuring DNA cell cycles, singlet gating was based on DNA-fluorescence A/H pulse shape analysis. Cell sorts were operated at 20psi/100um Nozzle using FACSFlow as sheath fluid and single cell precision was chosen when sorting cells for clonal expansion.

CUL3 Adaptor Screen

96-well imagining plates (ZellKontakt) were coated with the sgRNAs targeting CUL3 adapters (Key Supplementary Materials) complexed with Lipofectamine2000 (Thermo Fisher Scientific, 11668027). 2000 cells were added to each well. Following two rounds of transfection, cells were fixed with 4% paraformaldehyde solution and stained with the indicated antibodies.

Overview of MEMcrispR

Pooled CRISPR screens are typically performed with sgRNA and Cas9 transfer via virus infection, followed by selection according to a phenotype of interest. The selection enriches or depletes sgRNAs from the cell population. Multiple sgRNAs (5-10) are used to target the same gene. However, sgRNAs targeting the same gene differ in terms of knockout efficiency and potential off-target effect, creating a challenge in calculating a gene-level statistic from the sgRNA counts. Additional parameters influencing the effect of the knockout are batch effects introduced by the use of different library sets, the stochasticity introduced by the small number of biological replicates (typically 2-3), the cell line in which the experiment is performed, and the condition applied for the screen.

MEMcrispR uses mixed effects models to account for these parameters (Bates et al., 2015). Mixed effects models fit the observed values (here: counts at the gene level) to a linear combination of “fixed” and “random” effects. Fixed effects refer to the parameters of interest to the researcher, such as the effect of a gene’s knockout, whereas random effects represent factors that cannot be controlled by the researcher but influence the observed values independently of the fixed effects (i.e., randomly). Differences in the knockout efficiency of individual sgRNAs, as well as batch effects introduced by the use of different library sets were considered random effects. Treatment (the kind of selection applied in the screen) is the primary parameter of interest to the researcher and was therefore modeled as a fixed effect. In the analysis of screens performed in different cell lines, differences and commonalities in the effect between cell lines represent an additional parameter of interest. Thus, the cell line in which the screen is performed can be considered an additional fixed effect in the model.

The overview of the MEMcrispR tool is presented in Figure S1. The software takes raw sequenced data and, using the Subread aligner (Liao et al., 2013), aligns them to custom or genome-wide CRISPR libraries and normalizes according to sequencing depth, or to the median of each sample. The tool offers a variety of normalization methods, each recommended for particular types of screens. The tool assesses quality controls such as reproducibility, PCA analysis, clustering and counts densities before and after normalization. After normalization, each gene runs through a linear mixed effect model and results are visualized by volcano and clustering plots (Figure 1).

Single screen analysis: before and after treatment

The primary focus in the analysis of individual screens is the effect of the treatment (selection) on the cells with a particular gene knockout. We thus define the following models:Counts correspond to the normalized sgRNA counts for each gene according to the method of choice. Treatment is a binary statement where 0 equals the condition before a selection and 1 the condition after the selection (growth, drug treatment and etc.). Guide is the random effect of the sgRNAs, where the efficiency of the sgRNA functionality is due to each sgRNA. The sgRNAs are grouped into library pools, in the case the multiple different libraries were used for the analysis. We compute the T and beta scores for each gene according to their treatment. To estimate P-values we perform an ANOVA likelihood ratio test between the alternative model (M1) and a null model (M0) containing the same random effects but not the treatment term. If the effect of the gene knockout is significant, the alternative and the null model will be substantially different from each other, resulting to a statistically significant ANOVA statistic. The estimated P-values are then corrected for multiple testing using the Benjamini-Hochberg procedure. The beta scores are converted to fold changes and are used for the directionality of the effect of each gene (positive or negative) as well as for the magnitude of the effect. Fold changes represent sgRNA read counts between treatment and control conditions.

Multiple screen analysis: global treatment effects

The power of the MEMcrispR algorithm is displayed when combining multiple screens, with the same experimental design but performed across multiple cell lines. The combination of data from multiple screens increases the sample size for each gene. This leads to higher sensitivity in detecting genes that exhibit the same effect between the different cell lines (for example core essential genes). The cell line effect can be added as a fixed effect to the model, along with the treatment, since the effect of a gene’s knockout is always measured in the background of the respective cell line. The following models are then considered for each gene:Similarly, to the single screen model, this analysis also aims to identify significant treatment effects. P-values are estimated through an ANOVA likelihood ratio test between the alternative model (M1) and a null model (M0). For effects that are consistent between cell lines due to treatment, pooling the data results in a substantial increase in the sample size. Increased sample size leads to more robust statistics and enhanced sensitivity. Hence, the algorithm can detect less pronounced effects and compensate for inefficient sgRNAs that lower the effect size. It is important to note that, in this model, even effects which are not consistent across all cell lines will give rise to significant statistics if the magnitude of the effect in one or more cell lines is large enough. Therefore, this model alone cannot distinguish between global and cell-specific effects, but instead reports any significant treatment effect.

Multiple screen analysis: Cell specific effects

Besides increasing the sensitivity of the analysis, the integration of screens performed in multiple cell lines allows MEMcrispR to detect cell line specific effects. To this end, MEMcrispR employs an additional model that incorporates an interaction term (cell specific effect) between the treatment and the cell line effect. Similarly, to the treatment effect model, the statistic for the cell specific effect is calculated by comparing an alternative (M1) and a null model (M0) through an ANOVA likelihood ratio test. In this model, the treatment and the cell line effect are present in both models, however their interaction (cell specificity) is only included in the alternative model. As before, if the cell specific effect is important, the two models will differ considerably leading to a significant ANOVA statistic.The cell specific effects are estimated as pairwise comparisons between a baseline cell type, which can be determined by the user, and every other cell line in the dataset. Then the cell specific effects are illustrated by P-value and fold change to the baseline cell line. Thus, the directionality of the cell specific effect is always relevant to the baseline cell type and therefore requires careful assessment when analyzing the results.

Integration of the cell line effect models

The two models mentioned above are valuable on their own, increasing the sensitivity in detecting treatment effects and estimating cell type specificity respectively, but a major feature of MEMcrispR lies in integrating the results of both models to draw more complex conclusions about the dataset. The combination of the treatment and interaction statistic for each gene allows the distinction between global effects, such as core essentiality (significant treatment but no interaction), genes with the same effect but differences in its magnitude (both treatment and interaction are significant) and genes with cell specific effects (significant interaction but insignificant overall treatment). MEMcrispR contains a set of functions that allow the user to identify and visualize cell-type specific effects.

Quantification and Statistical Analysis

Statistical tests were performed using Wilcoxon-Mann-Whitney test unless stated otherwise. The standardized effect sizes were calculated with Cohen’s D method (difference of the means in the two groups divided by the pooled standard deviation). An effect size over 0.8 suggests that the two groups have a large effect size (Sawilowsky, 2009).

The number of the independent replicates are written on the figure legends or in the Methods section. Flow cytometry data were quantified using the median from the cell fluorescent intensities. Microscopy imagens were quantified using Fiji-ImageJ (Schindelin et al., 2012). Gene set enrichment analysis was performed using the GSEA tool from the Broad institute. Gene Ontology was conducted by Toppgene.

sgRNAs were aligned and P-values were calculated from the mixed effect model by MEMcrispR as described in the Methods section.

For the analysis of mutation co-occurrence of human tumor genomes, we obtained data from the TCGA Pancan Atlas project (Hoadley et al., 2018) and analyzed data from 12830 patients where the mutation data was available. We inferred damaging somatic CUL3 mutations, identified as stop-gain, frameshift or splice acceptor site mutations to investigate mutational co-occurrence with TP53 loss of function mutations that also include known missense mutations that are damaging. We employed Fisher’s exact test to investigate mutation co-occurrence.

Data and Code Availability

Sequencing count data are available in Tables S1, S3, and S4.

MEMcrispR is available at: https://github.com/grimbough/MEMcrispR

Raw sequencing data is available at:

https://www.ebi.ac.uk/ena/data/view/PRJEB36947

Raw data from Figures 2, 3, S4, and S5 are available at Mendeley Data: https://doi.org/10.17632/pm5v6pwkkg.1

REAGENT or RESOURCE	SOURCE	IDENTIFIER
Antibodies
Anti-GAPDH (D16H11) XP® Rabbit monoclonal Ab	Cell Signaling Technology	Cat#5174; RRID:AB_10622025
Anti-Cul3 Rabbit polyclonal antibody	Bethyl Laboratories, Inc.	Cat#A301-109A; RRID:AB_873023
Anti-p53 Antibody (DO-1) Mouse monoclonal antibody	Santa Cruz Biotechnology	Cat#sc-126; RRID:AB_628082
Anti-IKKα (3G12) Mouse monoclonal antibody	Cell Signaling Technology	Cat#11930; RRID:AB_2687618
Anti-IKKβ (D30C6) Rabbit monoclonal antibody	Cell Signaling Technology	Cat#8943; RRID:AB_11024092
Anti-Phospho-IKKα/β (Ser176/180) (16A6) Rabbit monoclonal antibody	Cell Signaling Technology	Cat#2697; RRID:AB_2079382
Anti-NF-κB p65 (D14E12) XP® Rabbit monoclonal antibody	Cell Signaling Technology	Cat#8242; RRID:AB_10859369
Anti-Phospho-NF-κB p65 (Ser536) (93H1) Rabbit monoclonal antibody	Cell Signaling Technology	Cat# 3033; RRID:AB_331284
Anti-α-Tubulin Mouse monoclonal antibody, clone DM1A	Sigma-Aldrich, Inc.	Cat#T6199; RRID:AB_477583
Anti-Smad4 (D3M6U) Rabbit monoclonal antibody	Cell Signaling Technology	Cat#38454; RRID:AB_2728776
Anti-SRF (D71A9) XP® Rabbit monoclonal antibody	Cell Signaling Technology	Cat#5147; RRID:AB_10694554
Anti-PD-L1 (E1L3N®) XP® Rabbit monoclonal antibody	Cell Signaling Technology	Cat#13684; RRID:AB_2687655
Anti-c-Jun (60A8) Rabbit monoclonal antibody	Cell Signaling Technology	Cat#9165; RRID:AB_2130165
Anti-FRA1 (D80B4) Rabbit monoclonal antibody	Cell Signaling Technology	Cat#5281; RRID:AB_10557418
Anti-Notch2 (D76A6) XP® Rabbit monoclonal antibody	Cell Signaling Technology	Cat#5732; RRID:AB_10693319
Anti-Notch3 (D11B8) Rabbit monoclonal antibody	Cell Signaling Technology	Cat#5276; RRID:AB_10560515
Anti-ADAM9 (D64B5) Rabbit monoclonal antibody	Cell Signaling Technology	Cat#4151; RRID:AB_1903892
Anti-DLL1 Antibody Rabbit polyclonal antibody	Cell Signaling Technology	Cat#2588; RRID:AB_2292961
Anti-DLL3 (G93) Antibody Rabbit polyclonal antibody	Cell Signaling Technology	Cat#2483; RRID:AB_2092839
Anti-Numb (C29G11) Rabbit monoclonal antibody	Cell Signaling Technology	Cat#2756; RRID:AB_2154298
Anti-RBPSUH (D10A4) XP® Rabbit monoclonal antibody	Cell Signaling Technology	Cat#5313; RRID:AB_2665555
Anti-TACE (D22H4) Rabbit monoclonal antibody	Cell Signaling Technology	Cat#6978; RRID:AB_10828385
Anti-phospho-Histone H2A.X (Ser139) Mouse monoclonal antibody, clone JBW301	Merck KGaA, Darmstadt, Germany	Cat#05-636-I; RRID:AB_2755003
Anti-53BP1 Rabbit polyclonal antibody	Bethyl Laboratories	Cat#A300-272A; RRID:AB_185520
Anti- Phospho-Histone H3 (Ser10) (D2C8) XP® Rabbit monoclonal antibody	Cell Signaling Technology	Cat# 3377; RRID:AB_1549592
Anti- β-Actin (13E5) Rabbit monoclonal antibody	Cell Signaling Technology	Cat# 4970; RRID:AB_2223172
Anti-Lamin A/C mouse monoclonal Antibody	Santa Cruz Biotechnology	Cat#sc-7292; RRID:AB_627875
Anti-mouse IgG, HRP-linked Antibody	Cell Signaling Technology	Cat#7076; RRID:AB_330924
Anti-rabbit IgG, HRP-linked Antibody	Cell Signaling Technology	Cat#7074; RRID:AB_2099233
Annexin V Alexa Fluor 488 conjugate	Thermo Fisher Scientific	Cat# A13201
Bacterial and Virus Strains
One Shot Stbl3 Chemically Competent E. coli	Thermo Fisher Scientific	Cat#C737303
Chemicals, Peptides, and Recombinant Proteins
Herculase II Fusion DNA Polymerase	Agilent	Cat#600677
BsmBI/Esp3I	Thermo Fischer Scientific	Cat#ER0451
Benzonase® endonuclease	Merck-Millipore	Cat#71206-3
Puromycin Dihydrochloride	Thermo Fisher Scientific	Cat#A1113803
M3814	MedKoo	Cat#206478
VE-822	MedKoo	Cat#406258
AZD0156	Selleckchem	Cat#S8375
Critical Commercial Assays
QIAamp DNA Blood Maxi Kit	QIAGEN	Cat#51194
NucleoSpin® Gel and PCR Clean-up Kit	Macherey-Nagel GmbH & Co. KG	Cat#740609
Quick Ligation Kit	New England Biolabs	Cat#M2200
QIAprep Spin Miniprep Kit	QIAGEN	Cat#27104
DNeasy Blood & Tissue Kit	QIAGEN	Cat#69504
RNAiMAX	Thermo Fisher Scientific	Cat#13778030
CellTiter-Glo® 3D Cell Viability Assay	Promega	Cat#G9681
NEBNext® Ultra II DNA Library Prep Kit for Illumina	New England Biolabs	Cat#E7645
Deposited Data
Raw sequencing data	https://www.ebi.ac.uk/ena	Accession number: PRJEB36947
Raw data from Figures 2, 3, S4, and S5	Mendeley Data	https://doi.org/10.17632/pm5v6pwkkg.1
Experimental Models: Cell Lines
Human: hTERT RPE-1	ATCC®	CRL4000, RRID:CVCL_4388
Human: MCF 10A	ATCC®	CRL-10317, RRID:CVCL_0598
Human: 293FT	Thermo Fisher Scientific	Cat# R70007
RPETP53−/−	This manuscript	N/A
MCF10ATP53−/−	This manuscript	N/A
Oligonucleotides (See Table S5 for oligonucleotide information)
Recombinant DNA
Plasmid: psPAX2 (gift from Didier Trono)	Addgene	Cat#12260, RRID:Addgene_12260
Plasmid: pMD2.G (gift from Didier Trono)	Addgene	Cat#12259, RRID:Addgene_12259
Plasmid: lentiCRISPRv2-GFP (gift from David Feldser)	Addgene	Cat#82416, RRID:Addgene_82416
Plasmid: pLentiCRISPR-2G-mCherry	This manuscript	Available upon request
gRNA pooled library in lentiCRISPRv2 (gift from Feng Zhang)	Addgene	Cat#1000000048
Software and Algorithms
FlowJo v10.2	FlowJo, LLC	RRID:SCR_008520
Fiji - ImageJ	https://imagej.net/Fiji	RRID:SCR_002285
MEMcrispR	This manuscript	https://github.com/grimbough/MEMcrispR
R Project for Statistical Computing	http://www.R-project.org	RRID:SCR_001905
RStudio	https://rstudio.com/	RRID:SCR_000432
Gene Set Enrichment Analysis (GSEA)	https://www.gsea-msigdb.org/gsea/index.jsp	RRID:SCR_003199
Other
Polybrene Infection / Transfection Reagent	Sigma-Aldrich	Cat#TR-1003
1.8x Agencourt AMPure XP	Beckman Coulter	Cat#A63880
Sera-Mag Beads	Thermo Scientific	Cat#4515-2105-050250
Sytox blue dead cell stain	Thermo Fisher Scientific	Cat# S34857
Lipofectamine 3000 kit	Thermo Fisher Scientific	Cat#L3000001