https://orcid.org/0000-0001-6403-3926
Abstract
Combinatorial CRISPR-Cas screens have advanced the mapping of genetic interactions, but their experimental scale limits the number of targetable gene combinations. Here, we describe 3Cs multiplexing, a rapid and scalable method to generate highly diverse and uniformly distributed combinatorial CRISPR libraries. We demonstrate that the library distribution skew is the critical determinant of its required screening coverage. By circumventing iterative cloning of PCR-amplified oligonucleotides, 3Cs multiplexing facilitates the generation of combinatorial CRISPR libraries with low distribution skews. We show that combinatorial 3Cs libraries can be screened with minimal coverages, reducing associated efforts and costs at least 10-fold. We apply a 3Cs multiplexing library targeting 12,736 autophagy gene combinations with 247,032 paired gRNAs in viability and reporter-based enrichment screens. In the viability screen, we identify, among others, the synthetic lethal WDR45B-PIK3R4 and the proliferation-enhancing ATG7-KEAP1 genetic interactions. In the reporter-based screen, we identify over 1,570 essential genetic interactions for autophagy flux, including interactions among paralogous genes, namely ATG2A-ATG2B, GABARAP-MAP1LC3B and GABARAP-GABARAPL2. However, we only observe few genetic interactions within paralogous gene families of more than two members, indicating functional compensation between them. This work establishes 3Cs multiplexing as a platform for genetic interaction screens at scale.
INTRODUCTION
Cellular homeostasis results from the orchestrated interplay of multiple genes and varies with environmental context. Accordingly, the systematic exploration of genetic interactions (GIs) in human cells could provide substantial insights into the causative effects of synergistic gene function related to complex phenotypes and diseases. The development of advanced gene editing technologies, most notably the clustered regularly interspaced short palindromic repeats-CRISPR associated protein (CRISPR-Cas) system, has enabled large-scale and systematic investigations of gene dependencies in human cell lines (1,2). However, the required experimental scale for GI screens can limit the number of simultaneously investigated gene combinations. Three parameters determine the scale, robustness and reproducibility of CRISPR screens: library diversity, library distribution skew, and experimental coverage (3). The library diversity is defined as the total number of different gRNAs or gRNA combinations and increases with the number of desired target genes. The library distribution is the abundance of all gRNAs or gRNA combinations. In an ideal library, all gRNA sequences are uniformly abundant, a property unachieved due to practical constraints. However, uniformity can be measured with the distribution skew that is defined by the ratio of the highest and lowest abundant gRNAs, elsewhere defined as skew ratio or distribution width (3,4). The experimental coverage describes the average abundance of each gRNA or gRNA combination during the screen. Current guidelines for performing CRISPR screens suggest library coverages between 200- and 1,000-fold, but lack precise indications which coverage should be used with a given library (5–7). Many CRISPR platforms have been developed for generating monogenic and multiplex CRISPR libraries that meet the needs of diverse applications. Still, current methods for CRISPR library generation require PCR amplification and sequential cloning steps, both of which can lead to wide library distributions associated with high skew (3). Thus, given the increasing demand for monogenic and multiplex CRISPR libraries, innovative methods yielding CRISPR libraries with low skews are needed.
Bulk and selective autophagy are tightly regulated processes that target cellular material for lysosomal degradation, and any misregulation can lead to abnormal cell growth or cell death, and have severe implications in human diseases (8,9). Rationally-engineered fluorescent reporter systems in combination with high-throughput CRISPR screens facilitated the systematic identification of genes essential for autophagic activity (10). As such, TMEM41B, UBA6 and BIRC6 were recently identified as contributors to the autophagy network (11–16). Besides mapping core autophagy genes, CRISPR screens, paired with fluorescent reporters, revealed the stress-dependent regulation of autophagy-related genes (17–19), and provided mechanistic insights into gene specificity in bulk and selective autophagy by identifying PARKIN/PARK2 regulators and the ANT complex as essential mitophagy components (20,21). In combination with proximity biotinylation-coupled mass spectrometry, CRISPR screens also contributed to a spatial proteogenomic understanding of PARK2-dependent mitophagy (22). Thus, unbiased CRISPR approaches, coupled to fluorescent reporters and mass spectrometry, are valuable approaches to study protein and gene function in bulk and selective autophagy (23–25). Moreover, the autophagy network contains several paralogous gene (short paralog) families for which selective or compensatory functions are largely elusive. In line with this, recent work revealed a systematic underrepresentation of paralogs in monogenic CRISPR screens, highlighting the need of multiplex CRISPR approaches for their functional characterization (26,27).
Here, we describe covalently-closed circular-synthesized (3Cs) multiplexing, a rapid and scalable method to generate combinatorial CRISPR libraries. We demonstrate that 3Cs multiplexing decouples the library distribution skew from its diversity, resulting in low distribution skews, even for highly diverse combinatorial CRISPR libraries. In pairwise drop-out screens, we identify the library distribution skew to be the major determinant for its required experimental scale. We demonstrate that, due to their uniform distribution, 3Cs libraries can be applied with minimal screening conditions, resulting in a 10-fold reduction of associated efforts and costs. Applying our minimized screening conditions, we investigated 12,736 pairwise gene knockouts (247,032 paired gRNAs) of human autophagy genes. We screened for synergistic gene combinations affecting cell fitness and identified WDR45B-PIK3R4 and ATG7-KEAP1 to result in reduced and enhanced cell proliferation, respectively. Moreover, we applied the multiplex autophagy library in a fluorescence-activated cell sorting (FACS)-based enrichment screen and identified, among others, the GI of the paralog pair ATG2A-ATG2B, confirming the essential role of this interaction for autophagy flux.
MATERIALS AND METHODS
3Cs multiplex template plasmid DNA and cloning
pLentiCRISPRv2 (Addgene: 98290) was enzymatically digested with AleI (New England Biolabs) and BsiWI (New England Biolabs) and gel purified to remove the human U6 (hU6) gRNA- and Streptococcus pyogenes Cas9 (SpCas9)-expressing cassettes. Likewise, the combinatorial gRNA-expressing cassette of pKLV2.2 (Addgene: 72666) was digested with AleI and BsiWI. The 2030 bp fragment that encoded the combinatorial gRNA-expressing cassettes and a PGK promoter were gel purified and cloned into the Cas9-excised, purified backbone of pLentiCRISPRv2. In order to generate unique annealing homology for the 3Cs oligonucleotides and enable template plasmid removal, the human 7SK (h7SK) promoter-associated tracrRNA was replaced by a previously engineered tracrRNA sequence (tracrRNAv2). The h7SK and hU6 promoter-associated gRNA cloning sites were modified to contain placeholder sequences encoding I-CeuI and I-SceI homing endonuclease restriction sites, respectively (28).
3Cs oligonucleotide design rules
To discriminate between h7SK and hU6 and enable exclusive annealing to only one gRNA-expression cassette, the 3Cs oligonucleotides were designed with two distinct homology regions, flanking the intended 20 nt gRNA sequence for either the h7SK or hU6 expression cassettes. The 3Cs h7SK oligonucleotides were 57 nucleotides in length (Tm above 50°C) and matched the 3′ end of the h7SK promoter region and the 5′ start of the tracrRNAv2, while the 3Cs hU6 oligonucleotides were 59 nucleotides in length (Tm above 50°C) and matched the 3′ end of the hU6 promoter region and the 5′ start of the wildtype SpCas9-tracrRNA in the template plasmids.
Generation of artificially skewed 3Cs libraries
For the generation of biased 3Cs multiplex gRNA libraries, two 3Cs oligonucleotide pools were designed, one for each expression cassette of the 3Cs multiplex template plasmid following the 3Cs oligonucleotide design rules. The first pool was composed of tumor suppressor and essential gene-targeting gRNAs (target pool), whereas the second pool only consisted of non-human targeting (NHT) gRNAs (control pool). To generate different libraries with varying sequence distributions, the two oligonucleotide pools were mixed in different ratios. For the first library, the target and control pool were mixed in a 1:1 ratio, to resemble an evenly distributed gRNA library. For the second library, a 1:10 ratio of target to control pool was applied, while the third library was generated with a 1:100 ratio to resemble a library with highly underrepresented gRNA sequences. The mixed oligonucleotide pools were phosphorylated, annealed to purified deoxyuridine-single stranded DNA (dU-ssDNA) of the 3Cs multiplex template plasmid, and the 3Cs synthesis reactions were performed.
Generation of dU-ssDNA
Chemically competent bacteria Escherichia coli strain K12 CJ236 (New England Biolabs, E4141) were transformed with 100 ng 3Cs multiplex template plasmid and grown on LB agar plates with ampicillin (100 μg/ml) and chloramphenicol (34 μg/ml) overnight at 37°C. The next day, single bacterial colonies were grown in 1 ml 2YT medium (Carl Roth GmbH) supplemented with 1×108 pfu/ml M13KO7 helper phage (New England Biolabs) and ampicillin (100 μg/ml) to maintain the host F′ episome and the phagemid, respectively. After 2 h of shaking at 200 rpm and 37°C, 25 μg/ml kanamycin (Carl Roth GmbH) was added to select for bacteria infected with M13KO7 helper phage. Bacteria were kept at 200 rpm and 37°C for 6 to 8 h. Afterwards, the culture was transferred to 30 ml of 2YT media supplemented with 100 μg/ml ampicillin and 25 μg/ml kanamycin. After an additional 20 hrs of shaking at 200 rpm and 37°C, the bacterial culture was centrifuged for 10 min at 10,000 rpm and 4°C in a Beckman JA-12 fixed angle rotor. The supernatant was subsequently transferred to 6 ml PEG/NaCl solution (20% polyethylene glycol 8000, 2.5 M NaCl) and incubated for 1 hour at room temperature to precipitate phage particles. After 10 min of centrifugation at 10,000 rpm and 4°C in a Beckman JA-12 fixed angle rotor, the phage pellet was resuspended in 1.5 ml Dulbecco's phosphate-buffered saline (PBS, Sigma-Aldrich) and centrifuged at 13,000 rpm for 5 min, before the phage-containing supernatant was transferred to a clean 1.5 ml microcentrifuge tube and stored at 4°C. Circular ssDNA was purified from the resuspended phages with the E.Z.N.A. M13 DNA Mini Kit (Omega Bio-Tek, D69001-01) according to the manufacturer's protocol. Purity of the isolated ssDNA was ensured by agarose gel electrophoresis and purified ssDNA was stored at 4°C.
Generation of heteroduplex dU-3Cs-dsDNA
The protocol for 3Cs multiplex-DNA synthesis was adapted from Wegner et al. (2019) and optimized for reactions on the 3Cs multiplex template plasmid with two specific annealing sites (29,30). The oligonucleotides that were used for 3Cs reactions and the suppliers are listed separately (Supplementary Table S1).
Oligonucleotide phosphorylation and annealing
600 ng of oligonucleotides per annealing site (both, 3Cs h7SK- and hU6-oligonucleotides) were phosphorylated in two separate 20 μl reactions by mixing them with 2 μl 10× TM buffer (0.1 M MgCl2, 0.5 M Tris-HCl, pH 7.5), 2 μl 10 mM ATP (New England Biolabs), 1 μl 100 mM DTT (Cell Signaling Technology Europe), 20 units of T4 polynucleotide kinase (New England Biolabs) and water to a total volume of 20 μl. The mixture was incubated for 1 h at 37°C. Phosphorylated oligonucleotides were immediately annealed to purified multiplex dU-ssDNA template by adding both 20 μl phosphorylation products to 25 μl 10× TM buffer, 20 μg of dU-ssDNA template and water to a total volume of 250 μl. The mixture was denatured for 5 min at 95°C, annealed for 5 min at 55°C and cooled down for 10 min at room temperature.
3Cs multiplex-DNA synthesis
3Cs-DNA was generated by adding 10 μl of 10 mM ATP, 10 μl of 100 mM dNTP mix (Carl Roth GmbH), 15 μl of 100 mM DTT (Cell Signaling), 2,000 ligation units of T4 DNA ligase (New England Biolabs), and 30 units of unmodified T7 DNA polymerase (New England Biolabs) to the annealed oligonucleotide/ssDNA mixture. The 3Cs synthesis mix was incubated for 12 h (overnight) at 22°C. The 3Cs synthesis product was purified and desalted using the GeneJET Gel Extraction Kit (Thermo Fisher), according to the following protocol: 600 μl of binding buffer and 5 μl of 3 M sodium acetate (Sigma-Aldrich) were added to the synthesis product, mixed and applied to two purification columns, which were centrifuged for 3 min at 460 g. The flow-through was applied twice to the same purification column to maximize yield. After two wash steps and final 3 min of centrifugation at maximum speed, the DNA was eluted in 50 μl prewarmed water. The 3Cs reaction product was analyzed by gel electrophoresis alongside the dU-ssDNA template on a 0.8% TAE/agarose gel (100 V, 30 min).
Electroporation and determination of bacterial transformation efficiency
To amplify the 3Cs multiplex libraries, 6 μg of purified 3Cs dsDNA synthesis product was electroporated into 400 μl electrocompetent E. coli (10-beta, New England Biolabs, C3020K) by using a Bio-Rad Gene Pulser (resistance 200 Ω, capacity 25 F, voltage 2.5 kV). After electroporation, cells were rescued in 25 ml of pre-warmed SOC media and incubated for 30 min at 37°C and 200 rpm. After 30 min, the culture was transferred into 400 ml of LB media supplemented with 100 μg/ml ampicillin and shaken overnight at 37°C. To ensure library representation of at least 1,000-fold during and after amplification, the number of transformants was determined by preparing serial dilutions of electroporated bacteria in sterile Dulbecco's phosphate-buffered saline (PBS, Sigma-Aldrich) that were plated in triplicates on LB agar plates containing 100 μg/ml ampicillin and incubated overnight at 37°C. The next morning, the obtained colonies were counted. The number of transformants had to be at least 100-fold higher than the library complexity to maintain library diversity.
I-CeuI and I-SceI clean-up and 3Cs library quality control
Plasmid DNA of overnight liquid cultures was purified using a Maxi Plasmid DNA Prep Kit (Qiagen) according to the manufacturer's protocol to obtain the pre-library. To remove residual 3Cs template plasmid from the multiplex pre-library, 3 μg of purified DNA were digested with 10 units of I-SceI and I-CeuI (New England Biolabs) for a total of 6 h at 37°C, with an additional addition of 10 units of I-SceI and I-CeuI to ensure complete removal of wild type sequences. The digestion reaction was then subjected to gel electrophoresis on a 0.8% TAE/agarose gel (125 V, 40 min) to separate undigested 3Cs synthesis product from linearized template plasmid. The band resembling the undigested 3Cs synthesis product was purified using a Thermo Fisher Scientific GeneJET Gel Extraction Kit, according to the manufacturer's protocol. Then, the purified 3Cs synthesis product was electroporated, according to the electroporation protocol described above. The next day, the resulting final 3Cs multiplex library preparation was purified from liquid culture using a Maxi Plasmid DNA Prep Kit (Qiagen), according to the manufacturer's protocol and quality-controlled by analytical restriction enzyme digests.
Sequencing
Unless otherwise specified, all DNA sequencing experiments were performed with Illumina technology. Gel-purified PCR products of plasmid libraries or screening samples were denatured and diluted according to Illumina guidelines and set to a final concentration of 2.6 pM in a total volume of 2.2 ml and 15% PhiX control and loaded onto a MiSeq, NextSeq500 or NovaSeq sequencer (Illumina) depending on required read counts (500- to 1,000-fold sequencing depth), according to the manufacturer's protocol. Sequencing was performed with single- or paired-end reads, 75 or 150 cycles, respectively, plus 8 cycles of index reading.
Sample preparation for sequencing of plasmid libraries
3Cs multiplex plasmid libraries were prepared for sequencing as follows: 250 ng of plasmid DNA were used per PCR reaction in a final volume of 50 μl, containing 25 μl Next High-Fidelity 2× PCR Master Mix (New England Biolabs) and 2.5 μl of each 10 μM primer. Depending on the library complexity, up to four 50 μl reactions were performed. Primer sequences are listed separately (see ‘DNA oligonucleotides’). Thermal cycler parameters were set as follows: initial denaturation at 98°C for 5 min, 15 cycles of denaturation at 98°C for 30 s, annealing at 65°C for 30 s, extension at 72°C for 40 s, and final extension at 72°C for 5 min. PCR products were purified from a 1.5% TAE/agarose gel using a GeneJet Gel Extraction Kit (Thermo Fisher Scientific), according to the manufacturer's protocol.
Sample preparation for sequencing of screen samples
To prepare samples derived from screens for subsequent sequencing, the required amount of genomic DNA for sufficient coverage was calculated first. For the autophagy single and multiplex FACS samples, the amount of required genomic DNA was calculated as the number of FACS sorted cells x screening coverage × 6.6 pg. For the autophagy multiplex proliferation control samples, the amount of required genomic DNA was determined by calculating the library complexity × screening coverage × 6.6 pg. For samples derived from screening with the biased libraries, the amount of required genomic DNA was calculated as the library complexity × 200 (maximum screening coverage) × 6.6 pg DNA. The calculated amount of genomic DNA was used in the first PCR reaction (PCR1) with 2–4 μg of genomic DNA per 50 μl final reaction, using the Next High-Fidelity 2× PCR Master Mix (New England Biolabs) and 2.5 μl of each 10 μM PCR1 primer. Thermal cycler parameters were set as follows: initial denaturation at 98°C for 5 min, 15 cycles of denaturation at 98°C for 55 s, annealing at 65°C for 55 s, extension at 72°C for 110 s, and final extension at 72°C for 7 min. After PCR1, 25 μl of PCR1 product was transferred to a second PCR reaction (PCR2) in a 100 μl reaction with 50 μl High-Fidelity 2× PCR Master Mix and 5 μl of 10 μM PCR2 primers containing Illumina adaptors and barcodes. Primer sequences for PCR1 and PCR2 are listed separately in Supplementary Table S1. Thermal cycler parameters were set as follows: initial denaturation at 98°C for 5 min, 10 cycles of denaturation at 98°C for 30 s, annealing at 65°C for 30 s, extension at 72°C for 40 s, and final extension at 72°C for 5 min. PCR products were purified from a 1.5% TAE/agarose gel and processed for sequencing as described for plasmid libraries.
Sequencing data quality control and read count table generation
Raw sequencing data were processed and demultiplexed with bcl2fastq v2.19.1.403 (Illumina). Read counts of individual gRNAs and gRNA combinations were determined using cutadapt 2.8, Bowtie2 2.3.0, and custom Python 3 scripts (31,32). In brief, reads were trimmed with cutadapt, truncated to 20 nucleotides, and aligned to the respective gRNA library using Bowtie2 with no mismatches allowed. The uniformity of each library distribution was assessed by plotting the cumulative distribution of all sequencing reads as a Lorenz curve and determining the area under the curve (AUC).
The library distribution skew of each library was determined by plotting the density of read counts and dividing the 90th percentile by the 10th percentile. We applied Cohen's d statistics to assess the quality score (QS) of biased library screening performance by measuring the separation of mean LFC values of non-targeting sequences and sequences targeting core essential genes: QS = ((mean LFC of NHT combinations) - (mean LFC of essential combinations)) / (standard deviation (LFC of essential combinations)) (33). Pairwise sample correlations were determined with Pearson's correlation of the normalized read counts and visualized with hierarchically clustered heat maps using the Seaborn library 0.10.1 (34).
Enrichment analyses
Enrichment analyses using MAGeCK were performed with median or total normalization of read counts with gRNAs having zero counts in the control samples being removed (35). Down-sampling of the 1:1 dataset was performed by randomly choosing 1–16 gRNA combinations per gene combination without replacement followed by individual MAGeCK analyses to obtain false discovery rates (FDR) and log2 fold changes (LFC). gRNA combinations with an FDR ≤ 10% and LFC ≤ -0.5 were counted as statistically significant hits.
Genetic interaction models
GIs were computed based on five different models: SUM, MIN, LOG, MULT (adapted from (36)), and MAX. Each model assigned a GI to a gene pair xy, if the double mutant phenotype W(xy), deviated from a predicted double mutant phenotype that is expected for no interaction between x and y, E(xy). The phenotype was measured as the LFC of gRNA abundance between respective samples. The expected double mutant phenotype for joint mutations of the genes x and y was defined by the neutrality function of each definition with MIN = min(W(x), W(y)), MULT = (W(x) × W(y)), LOG = log2[(2W(x) + 1) × (2W(y) – 1) +1], SUM = (W(x) + W(y)), and MAX = max(W(x), W(y)). The single mutant phenotype for a gene x was defined as the median LFC of all gRNA combinations of NHTs and gRNAs targeting x. For each model, the deviation of observed double mutant phenotypes from their expectation was calculated as their difference and termed delta log2 fold change (dLFC): dLFC = observed – expected. To select the best model for our data, we assumed that GIs were rare for randomly selected gene pairs. Density plots of the dLFC for each model were used to identify the model with the highest number of neutral interactions, indicated by a single large peak around 0 on the x-axis. We kept only combinations with P ≤ 0.05 and a dLFC larger than the standard deviation (SD) of all dLFCs. To generate the GI heatmap of core autophagy paralogs, we used MAX model-derived GIs with dLFC > 0 and LFC > 0. Since possible combinations of two genes x and y are x–y or y–x and both were assumed to show the same phenotype, we averaged xy, yx pairs.
Autophagy gene interaction network
To generate a network visualization based on our derived gene interactions in autophagy, we exported MAX model-dependent dLFC per gene interaction and imported them into the open-source software platform for visualizing complex networks, Cytoscape v3.8.0 (37). The style of the derived network was manually adjusted.
Cell culture
Cell culture work was performed as described previously (29). In brief, HEK293T cells (ATCC, CRL-3216) were maintained in Dulbecco's modified Eagle's medium (DMEM, Thermo Fisher Scientific) and puromycin-sensitive human telomerase-immortalized retinal pigmented epithelial (RPE1) cells (provided by Andrew Holland) in DMEM: Nutrient Mixture F-12 (DMEM/F12, Thermo Fisher Scientific), each supplemented with 10% fetal bovine serum (FBS, Thermo Fisher Scientific) and 1% penicillin-streptomycin (Sigma-Aldrich) at 37°C with 5% CO2. In addition, RPE1 cells were supplemented with 0.01 mg/ml hygromycin B (Capricorn Scientific). THP1 cells were cultured in RPMI media supplemented with 10% fetal bovine serum (FBS, Thermo Fisher Scientific) and 1% penicillin-streptomycin (Sigma-Aldrich) at 37°C with 5% CO2. No method to ensure the state of authentication has been applied. Mycoplasma contamination testing was performed immediately after the arrival of the cells and multiple times during the course of the experiments by using Venor®GeM Classic kit from Minerva Biolabs GmbH, according to the manufacturer's instructions. The RPE1 Cas9 and GFP-LC3-RFP reporter cell line was generated by transducing RPE1 cells with lentiviral particles generated with the transfer plasmids pMRX-IP-GFP-LC3-RFP (Addgene: 84573) and lentiCRISPRv2(NHT) (Addgene: 49535). Single cell clones were isolated and reporter functionality was tested by Torin1 and Bafilomycin A1 treatments.
Genomic DNA extraction
Genomic DNA of cells was purified by resuspending PBS-washed pellets from 40–50 million cells in 12 ml of TEX buffer (10 mM Tris-HCl, pH 7.5; 1 mM EDTA, pH 7.9; 0.5% SDS) supplemented with 300 μl proteinase K (10 mg/ml, Carl Roth GmbH) and 300 μl ribonuclease A (90 U/mg, 20 mg/ml, Carl Roth GmbH), and incubated overnight at 37°C at constant shaking. After complete cell lysis, 4 ml of 5 M NaCl were added, the solution was mixed and incubated at 4°C for 40 min and centrifuged at 14,000 × g for 1 h at 4C. The supernatant was transferred to a fresh tube and 24 ml of ice-cold 96% ethanol were added before the mixture was placed at -20°C overnight. The next day, the tubes were centrifuged at 14,000 × g for 1 h. Afterwards, the supernatant was removed and the precipitated DNA was washed with ice-cold 70% ethanol. After further centrifugation at 14,000 × g for 1 h, the supernatant was removed, and the DNA pellet dried and dissolved in 5 ml of sterile water.
Generation of GFP and mCherry 3Cs libraries and knockout cells
To examine the expression of gRNAs from both expression cassettes of the 3Cs multiplex template plasmid, two oligonucleotide pools were designed following the 3Cs oligonucleotide design rules. One oligonucleotide pool was designed for the h7SK cassette with 50 gRNAs targeting GFP, the second pool was designed for the hU6 expression cassette with 50 gRNAs targeting the mCherry gene. The two pools were used to generate three 3Cs libraries to selectively target either GFP (GFP single library), mCherry (mCherry single library) or both simultaneously (GFP-mCherry multiplex library) following the protocol for generation of 3Cs multiplex gRNA libraries as described above. Lentiviral supernatant of the three libraries was generated. Monoclonal RPE1 cells with stable SpCas9, GFP and mCherry expressions were plated at 40% confluency. The next day, the cells were transduced with viral supernatant of one of the three libraries. After 48 h of transduction, cells were selected with puromycin (2.5 μg/ml) for 10 days. Finally, GFP and mCherry ratios were quantified by FACS analysis.
Generation, quantification and transduction of lentiviral particles
Generation, quantification and transduction of lentiviral particles were performed as described previously (29). In brief, the day before transfection, HEK293T cells were seeded to a density of 25,000 cells/ml. To transfect HEK293T cells, transfection media containing 1/10 of culture volume Opti-MEM I (Thermo Fisher Scientific), 10.5 μl Lipofectamine 2,000 (Thermo Fisher Scientific), 1.65 μg transfer vector, 1.35 μg/ml pPAX2 (Addgene: 12260) and 0.5 μg/ml μg pMD2.G (Addgene: 12259) were prepared. The mixture was incubated for 30 min at room temperature and added dropwise to the media. Lentiviral supernatant was harvested 48 hours after transfection and stored at –80°C.
To determine the lentiviral titer, RPE1 cells were plated in a 6-well plate with 50,000 cells per well. The following day, cells were transduced in the presence of 8 μg/ml polybrene (Sigma-Aldrich) and a series of 0.5, 1, 5 and 10 μl of viral supernatant. After 2 days of incubation at 37°C, cells were subjected to 2.5 μg/ml puromycin selection for a total duration of 2 weeks, after which established colonies were counted per viral dilution. The number of colonies in the highest dilution was then volume normalized to obtain the final lentiviral titer.
To transduce RPE1 cells, they were seeded at an appropriate density for each experiment with a maximal confluency of 60–70%. On the day of transduction, polybrene was added to the media to a final concentration of 8 μg/ml. The volume of lentiviral supernatant was calculated on the basis of the diversity of the respective library, and of the desired coverage and multiplicity of infection (MOI) of the experiment. A MOI of 0.5 was applied to all the screens. The number of cells that were transduced at the beginning of an experiment was calculated by multiplying the diversity of the library with the desired coverage and MOI.
3Cs CRISPR screening
Library distribution and experimental coverage screening
To explore the interdependency of multiplex CRISPR library distribution and experimental screening coverage, three distorted 3Cs multiplex libraries were generated that represented libraries of different gRNA distributions (see ‘Generation of distorted libraries’). All three libraries were screened with a 20-fold and 200-fold coverage, each in triplicates. For the 20-fold screening, for each replicate, 1.1 million SpCas9-expressing RPE1 cells were plated (0.37 million cells per flask) and transduced with the respective library with a MOI of 0.5. After 48 hours, cells were selected with 2.5 μg/ml puromycin and kept in growing conditions for 14 days. At day 14, cells were harvested, pooled and stored at –20°C until their genomic DNA was extracted and processed for sequencing, as described above. For the 200-fold screening, a total of 11 million (0.5 million cells per flask) SpCas9-expressing RPE1 cells were plated and transduced with the respective library with a MOI of 0.5. Further screening was performed identically to the 20-fold screen.
Combinatorial 3Cs gRNA autophagy screening
Single and combinatorial autophagy gRNA screens for single or synergistic autophagy inhibition were performed in biological triplicates in the monoclonal RPE1 cell line stably expressing SpCas9 and the autophagic flux probe (GFP-LC3-RFP) (10). For each replicate, 20 million cells (10 million for each, end time point and day 2 control) were transduced with lentiviral supernatant of the autophagy multiplex library with an MOI of 0.5 and a 1,000- or 20-fold library coverage for single or combinatorial autophagy library screening, respectively. The control time points were harvested 2 days post-transduction. The remaining cells were kept in growing conditions until day 7, at which point the cells were passaged, pooled and reseeded at a density maintaining library diversity. After 13, 14 and 15 days, the cells were treated with the mTOR inhibitor Torin1 (250 nM, InvivoGen) for 24 h to induce autophagy. Cells were collected and 150,000 to 300,000 cells for single, or 4.5–6.75 million cells for combinatorial screening were FACS-sorted to enrich cells with blocked autophagy. The sorted cells were reseeded and expanded for seven days before harvesting, pooled and stored at –20°C until their genomic DNA was extracted and processed for sequencing.
FACS
FACS was carried out with a BD FACSAria Fusion. CRISPR screen hit validation analysis was performed on a FACSCanto II flow cytometer (BD Biosciences). Data was processed with FlowJo (FlowJo, LLC). Gating was carried out on the basis of viable and single cells that were identified on the basis of their scatter morphology.
Arrayed autophagy candidate validation
The validation of single and combinatorial autophagy hits was performed in arrayed conditions (one knockout per well). To do so, single and dual gene-targeting CRISPR constructs were designed and generated. For each gene, the top scoring guide sequence was selected with Azimuth 2.0 of the GPP sgRNA Designer (38) and purchased as forward and reverse oligonucleotide with compatible overhangs for restriction enzyme cloning (Supplementary Table S1). The two oligonucleotides containing the gRNA target site were annealed and cloned into a restriction enzyme-digested and gel purified CRISPR vector. In more detail, single gRNA constructs were cloned into the lenti-sgRNA blast vector (Addgene: 104993) by BsmBI restriction enzyme cloning. For combinatorial hit validation, the dual CRISPR gRNA expression cassettes of pKLV2.2-h7SKgRNA5(SapI)-hU6gRNA5(BbsI)-PGKpuroBFP-W (Addgene: 72666) were cloned into the lenti-sgRNA blast plasmid to enable blasticidin selection of dual gRNA constructs. A silent point mutation was introduced to remove the BbsI (New England Biolabs) recognition site within the blasticidin sequence to allow the subsequent insertion of one gRNA by SapI (New England Biolabs) cloning into the h7SK expression cassette and the second gRNA by BbsI cloning into the hU6 expression cassette. After cloning, correct gRNA integration was verified by SANGER sequencing with Microsynth, Switzerland. Next, lentiviral supernatant was generated for each construct as described before. Monoclonal RPE1 cells with stable SpCas9 and GFP-LC3-RFP reporter expressions were plated in six-well plates with 50,000 cells per well. The following day, cells were transduced with lentiviral supernatant in the presence of 8 μg/ml polybrene (Sigma-Aldrich). After 48 hours, the cells were selected with 10 μg/ml blasticidin (InvivoGen) for 7 days, passaged and cultivated at 40–60% confluency under constant blasticidin selection for additional 7 days. At day 14, cells were treated with Torin1 to induce autophagy for 24 h, until they were collected at day 15 and subjected to FACS to measure single or dual gene knockout-induced autophagy blockage. Single and combinatorial autophagy gene deletions in THP1 cells were validated by lentiviral transduction, followed by puromycin selection. Proliferation was estimated by cell-counting at day 14 and basal autophagy flux was quantified by GFP-LC3-RFP reporter-coupled FACS without Torin1 treatment.
gRNA performance and TIDE assay
Guide RNA performance was evaluated by Tracking of Indels by Decomposition (TIDE) assay, as described previously (39). In short, for each gRNA sequence, PCR primers flanking the gRNA annealing site around 400 bp upstream and downstream were designed, resulting in a PCR product of 800 to 1,000 bp in length. The area around the gRNA-locus was PCR amplified with OneTaq DNA polymerase (New England Biolabs), using 1 μg of genomic DNA, 40 μM dNTPs (final concentration), 0.2 μM of each forward and reverse primer, 10x OneTaq standard buffer, and 2.5 units of OneTaq DNA polymerase. PCR cycles were set up as follows: initial denaturation at 94°C for 3 min, 39 cycles of denaturation at 94°C for 20 s, annealing at 55°C for 30 s, strand extension at 68°C for 2 min, and final strand extension at 68°C for 5 min. The PCR products were analyzed on a 0.8% TAE/agarose gel (100 V, 30 min) and purified using a Thermo Fisher Scientific GeneJET Gel Extraction Kit according to the manufacturer's protocol. The purified PCR product was pre-mixed with forward amplification primer and processed by SANGER sequencing with Microsynth, Switzerland, after which wild type and gRNA-treated SANGER chromatograms were analyzed by TIDE and the percentage of unedited DNA extracted (https://tide.nki.nl/) (39).
Lung squamous cell carcinoma patient survival analysis
Kaplan–Meier curves for lung squamous cell carcinoma (LUSC) patient five-year overall survival rates were estimated with the online tool KM plotter (https://kmplot.com) based on GEO, EGA and TCGA datasets (40–43). KEAP1 and ATG7 gene expression data were obtained from GEO, caBIG and TCGA databases. Using the KM plotter online tool, patients were split using the option ‘Auto select best cutoff’ in high or low KEAP1, ATG7 or KEAP1 and ATG7 gene expression groups. The following cutoff values were used: ATG7: 401, KEAP1: 846, KEAP1 and ATG7: 562. P values for log rank tests of the Kaplan–Meier curves were calculated with the KM plotter.
Survival of LUSC patients was analyzed by obtaining ATG7-KEAP1 gene expression data from TCGA in the UCSC Xena online tool (44,45). LUSC patients were classified into four categories depending on ATG7 and KEAP1 gene expression in tumors: (i) high expression of ATG7 and KEAP1; (ii) high expression of ATG7 and low expression of KEAP1; (iii) high expression of KEAP1 and low expression of ATG7; 4) low expression of ATG7 and KEAP1. Expression of ATG7 or KEAP1 was defined as high (low) when the respective expression levels were higher (lower) than the median expression levels in all LUSC tumors. Overall survival of a single patient was normalized to the median overall survival of all patients and the median survival were calculated for the previously described 4 groups. R2 linear regression was calculated using Microsoft Excel.
RNA-seq
Monoclonal RPE1(Cas9) cells were harvested at 90% confluency and total RNA was purified using the QIAGEN RNeasy Plus Mini Kit (Cat No.: 74134), according to the manufacturer's protocol. RNA was stored at -80°C and analyzed by RNA-Seq in quadruplicates. For library preparation, total RNA was quantified using the Qubit 2.0 fluorometric assay (Thermo Fisher Scientific). Sequencing libraries were prepared from 125 ng of total RNA using a 3'DGE mRNA-seq research grade sequencing service (Next Generation Diagnostics srl) which included library preparation, quality assessment and sequencing on a NovaSeq 6000 sequencing system using a single-end, 100 cycle strategy (Illumina Inc.) (46). The bioinformatics workflow included analysis of raw data by Next Generation Diagnostics srl proprietary 3'DGE mRNA-seq pipeline (v2.0) which involves a cleaning step by quality filtering and trimming, alignment to the reference genome and counting by gene (47–49). We filtered out all genes having < 1 cpm in less than n_min samples and Perc MM reads > 20% simultaneously. Differential expression analysis was performed using edgeR (Supplementary Table S7) (50).
DNA oligonucleotides and oligonucleotide pools
Oligonucleotides used for 3Cs reactions, cloning, and sequencing library preparation were obtained from Integrated DNA Technologies (IDT) or Merck KGaA, and are provided in Supplementary Table S1.
RESULTS
3Cs multiplex CRISPR gRNA libraries
To enable gRNA multiplexing based on our previously reported 3Cs technology (29), we generated a dual gRNA-expressing lentiviral plasmid, pLenti-Multiplex, by placing a h7SK promoter upstream of a previously engineered Cas9-tracrRNA, followed by a hU6 promoter upstream of the wild type SpCas9-tracrRNA sequence (51,52). Both gRNA cassettes contained placeholder sequences with I-CeuI and I-SceI recognition sites, respectively (Figure 1A) (29). In addition to a puromycin selection cassette, the plasmid contained an f1 bacteriophage origin of replication sequence in sense orientation, supporting the CJ236 bacteria and M13KO7 bacteriophage-dependent generation of dU-containing single stranded (ss) DNA (Figure 1A and Supplementary Figure S1A). In contrast to single 3Cs reactions, 3Cs multiplexing was performed with two gRNA-encoding oligonucleotide pools, each containing unique 5’ and 3’ homology sequences for gRNA placeholder-specific annealing (Figure 1A). This enabled gRNA-cassette-specific annealing of each oligonucleotide pool and the generation of dU-3Cs-dsDNAcontaining all possible combinations of the two gRNA pools. To obtain the final gRNA-containing dsDNA, remnants of the 3Cs template were enzymatically digested and the product was electroporated into dut/ung-positive bacteria and to test this strategy, we designed a multiplex gRNA library, targeting GFP (pool 1) and mCherry (pool 2), in addition to their respective single gRNA libraries (Figure 1B). Each of the two pools contained 50 gRNAs, yielding 2,500 possible gRNA combinations. Similar to single 3Cs-DNA (29,53), we observed a three-band pattern for 3Cs multiplex reaction-products (Supplementary Figure S1A), and digesting the libraries with I-CeuI and I-SceI enzymes confirmed the exclusive presence of gRNA-containing plasmids, with undetectable rates of 3Cs template plasmid (Figure 1C). Paired-end sequencing confirmed their completeness, and revealed distribution skews of 1.56, 1.17 and 1.53 for GFP, mCherry and GFP+mCherry libraries, respectively (Figure 1D, Supplementary Figure S1B, C and Table S2). Upon transduction of GFP and mCherry co-expressing RPE1(Cas9) cells, the GFP+mCherry multiplex library induced the simultaneous depletion of GFP and mCherry fluorescence, while the single gRNA libraries depleted GFP or mCherry fluorescence selectively (Figure 1E). Thus, 3Cs multiplexing generated combinatorial CRISPR gRNA libraries.
3Cs multiplexing for combinatorial CRISPR gRNA libraries. (A) 3Cs multiplexing workflow. Two gRNA-encoding oligonucleotide pools are annealed to either one of the gRNA-expression cassettes of the dU-containing single-stranded DNA template plasmid for T7 DNA polymerase-dependent generation of heteroduplex dU-dsDNA. 3Cs dU-dsDNA is amplified in dut/ung-positive bacteria. Wildtype SpCas9-tracrRNA (tracrWT). (B) Cas9 GFP/mCherry multiplex library design. Combinatorial gRNA constructs simultaneously target GFP and mCherry genes with 2,601 gRNA combinations; each gene with 50 gRNAs. The single gRNA GFP- or mCherry-targeting libraries contain the wild type gRNA-placeholder in the hU6 or h7SK cassette, respectively. (C) Gel-electrophoresis of final 3Cs single and multiplex libraries after restriction enzyme digestion. (D) Cumulative distributions of GFP and mCherry 3Cs multiplex gRNA libraries. A uniformly distributed library (ideal) is shown in grey. Percentages indicate library representations at 90% of cumulative reads. Area under the curve values are indicated next to each library identifier. (E) FACS analysis of GFP- and mCherry-positive RPE1 cells after transduction with single or combinatorial libraries. Error bars represent standard deviations over three biological replicates (n = 3).
3Cs multiplexing generates highly diverse and uniformly distributed libraries
Screening high numbers of target genes or gene combinations, requires a robust technology to generate highly diverse CRISPR libraries. Based on our previous work (29), we hypothesized that 3Cs multiplexing would facilitate the generation of uniformly distributed libraries irrespective of the diversity. To test this, we designed four oligonucleotide pools per gRNA cassette in which one, two, three, or four nucleotide positions, were randomized to mimic increasing combinatorial gRNA diversities (1N, 4*4 = 16 combinations; 2N, 16*16 = 256 combinations; 3N, 64*64 = 4,096 combinations; 4N, 256*256 = 65,536 combinations) (Supplementary Table S1). In addition, we also generated two non-randomized libraries with ∼247,000 and ∼913,000 gRNA combinations. Sequencing identified > 99% of gRNA pairs and confirmed uniform distributions of all libraries with AUC values between 0.59 and 0.7 (Figure 2A, Supplementary Figures S2A–B, S4A–C, S5, and Supplementary Table S3). Notably, we identified low distribution skews, ranging from 1.1 to 6.02, values most often unmatched even with single gRNA libraries (Supplementary Figures S2B and S4B). Therefore, we concluded that 3Cs multiplexing was robust, scalable, and obtained combinatorial CRISPR libraries with uniform library distributions and low distribution skew.
The library distribution skew determines the experimental scale. (A) Cumulative distributions of randomized and non-randomized 3Cs multiplex gRNA libraries. A uniformly distributed library (ideal) is shown in grey. Percentages indicate library representations at 90% of cumulative reads. Area under the curve (AUC) values are indicated next to each library identifier. (B) Fraction of NHT-reads in artificially skewed 3Cs multiplex libraries (gene to NHT ratios of 1:1, 1:10, 1:100). Distribution median and quartiles are displayed as red straight or black dotted lines, respectively. (C) Cumulative distributions of artificially skewed 3Cs multiplex libraries. A uniformly distributed library (ideal) is shown in grey. Percentages indicate library representations at 90% of cumulative reads. AUC values are indicated next to each library identifier. (D) Comparison of artificially skewed library performance. Density plots showing the separation of LFCs of combinatorial gRNA constructs targeting core essential genes (blue) or non-essential controls (green) in experimental coverages of 20-fold (20x) (dotted) and 200-fold (200x) (straight). QS: quality score. (E) Correlation of Cohen's d-based screen quality score and library distribution skew. Artificially skewed 3Cs multiplex libraries were screened at 20× (green circles) and 200× (blue circles) coverages. High screen performance (green area) is achieved when screen quality score is above 2 (horizontal grey dashed line). (F) Quantification of depleted gene pairs from (D), detected with MAGeCK at the indicated FDR and LFC cutoffs. (G) Subsampling analysis to determine the optimal number of pairwise gRNAs for statistical hit calling. Shown is the number of essential genes that are called in individual MAGeCK analyses on samples with 1 to 16 randomly chosen gRNAs. Analysis is based on the 1.2 skew library with 20x and 200x coverages.
The CRISPR library distribution skew determines the required coverage
A recent in silico study showed that a uniform library distribution increases replicate correlation (3). To verify this prediction in cells, we generated multiplex libraries with artificially skewed gRNA representations and screened them with different coverages. We selected a panel of 20 genes, comprising 10 core essential (CE) and 10 tumor-suppressor genes (TSGs), each targeted by 4 gRNAs plus 80 pre-validated NHT sequences that served as internal controls (38,54,55). Next, we designed two oligonucleotide pools. Pool-1 contained CE and TSG-targeting gRNAs for h7SK and hU6 promoters in equal ratios (2*80 gRNAs), and pool-2 contained exclusively NHT sequences in equal ratios for h7SK and hU6 promoters (2*80 gRNAs). Then, both oligonucleotide pools were combined in equimolar ratios (1:1), and in ratios of increasing NHT-sequence molarity (1:10 and 1:100). The 1:1 had equal ratios of gRNA combinations (CE/TSG compared to NHTs) and an AUC value of 0.65, while the 1:10 and 1:100 libraries were skewed with an increasing fraction of NHTs and AUC values of 0.85 and 0.9 (Figure 2B–C, Supplementary Figure S2C and Table S4). Of note, library skews increased from 1.2 (1:1) to 2.4 (1:10) and 13.46 (1:100) (Supplementary Figure S2D), demonstrating that the sequence distribution of an oligonucleotide pool directly translates to the 3Cs library distribution skew.
Next, we employed the libraries to cell proliferation screens in RPE1(Cas9) cells with coverages of 20- and 200-fold and evaluated a total of 153,600 paired gRNAs in biological replicates. Pairwise gRNA and gene-level counts correlated well for individual samples across biological replicates (gene level, 20× – 1:1, Pearson r = 0.96; 1:10, Pearson r = 1; 1:100, Pearson r = 0.79–1; 200× – 1:1, Pearson r = 0.91–0.98; 1:10, Pearson r = 0.96–0.97; 1:100, Pearson r = 0.98–1; Supplementary Figure S3A–B and Table S4). Compared to a single gRNA, two gRNAs directed against the same gene yield higher knockout rates (56). We therefore compared the LFCs of single versus dual-gRNA gene targeting and observed a consistently stronger LFC with dual gRNAs, an effect that was independent of the experimental coverage (Supplementary Figure S2E, F). Next, we aimed at analyzing the performance of each library at the 20- and 200-fold coverage by applying Cohen's d quality scores (QS). The QS was recently introduced to compare screen performance by measuring the separation of mean LFC values of gRNAs targeting either essential genes or non-essential genes (33).
The quality score for the 1:1 library in both coverages was high, with 2.95 and 3.25 for 20× and 200×, respectively (Figure 2D). However, we observed a decline in screen quality when the library distribution skew increased to > 2, demonstrating the library distribution skew to directly affect screen quality (1:10 library, QS 20×: 1.5, QS 200×: 2.68; 1:100 library, QS 20×: -0.44, QS 200×: -1.17; Figure 2D, E). Next, we assessed whether higher experimental coverage was able to rescue higher distribution skews by identifying essential gene pairs through MAGeCK analyses and cutoff filters at FDR ≤ 10% and LFC < -0.5. Strikingly, no apparent difference in significant gene pairs for the 1:1 and 1:10 libraries could be observed, while for the 1:100 library we observed an increased number of significantly depleted gene pairs (Figure 2F). This demonstrates that CRISPR libraries with low skews distribution skews can be screened at lower experimental coverage without compromising screen performance and hit calling accuracy. Still, we noted that even the 200-fold representation of the 1:100 library was largely inefficient in retrieving essential gene pairs, suggesting that an experimental coverage above 200-fold is required to rescue a library distribution skew > 10. We therefore concluded that the library distribution skew was an essential parameter that directly contributed to hit calling accuracy, with skew values of 2 or below supporting a 10-fold reduction in experimental efforts.
The optimal number of pairwise gRNAs for hit calling
The number of gRNAs per gene has a profound impact on statistical hit calling and experimental scale, and therefore on the robustness of screening results (57,58). To assess the optimal number of pairwise gRNAs for combinatorial CRISPR screens, we compared the concordance of essential gene pairs identified between the two experimental coverages of the 1:1 library screens (Supplementary Table S4). We down-sampled these data sets and generated a series of read count tables containing between 1 and 16 randomly chosen gRNA pairs for each gene combination and performed MAGeCK analyses with cutoff filters set to FDR ≤ 10% and LFC ≤ -0.5 (57,58). A total of 100 essential gene pairs were possible (10*10 essential genes), and as expected, 16 gRNA pairs consistently retrieved the highest number of significantly depleted gene pairs (200× –70%, 20× –72%) (Figure 2G). In contrast, one to four gRNA pairs resulted in low hit calling accuracy for both libraries (200-fold: 13%, 20-fold: 20%) (Figure 2G). However, hit calling was improved by increasing the number of paired gRNAs and plateaued between 13 and 16 pairs (Figure 2G). These results were consistent with previous observations for single gRNA screens in which four to six gRNAs per gene have been considered optimal for statistical hit calling (58). Furthermore, assuming the generation of balanced combinatorial libraries, our analysis identified 16 paired gRNAs (4*4) as the optimal number of gRNA pairs for statistical hit calling.
A 3Cs multiplex CRISPR library to investigate gene pairs in human autophagy
Having established 3Cs multiplex libraries and their minimized screening conditions, we applied our technology to investigate proliferative phenotypes and synergistic gene effects in human autophagy. To do so, we designed a multiplex autophagy library by assembling a literature-curated list of 64 core autophagy genes comprising (Supplementary Table S1) the ULK1 and PIK3C3 complexes, the conjugation machinery, ATG8 homologs, PI3P effectors and members of the autophagosome-lysosome fusion machinery. For each gene, we selected four gRNA sequences for the core autophagy pool, added 10% NHTs to obtain a total number of 282 sequences, and extended each sequence with 3Cs-homology arms, needed to correctly anneal to the 3Cs-multiplex vector (29). To generate a multiplex CRISPR library, we designed a second oligonucleotide pool, termed extended autophagy, that included all core autophagy genes/sequences plus, among others, autophagy related transcription factors, receptors and ubiquitin-specific proteases (USPs) (Supplementary Table S1). The extended autophagy pool contained 876 gRNAs, including four gRNAs per gene, 80 NHTs, and 28 control gRNAs targeting essential genes (Figure 3A). We applied the extended autophagy pool alone and in combination with the core autophagy pool to 3Cs reactions and generated both gRNA libraries, with the multiplex library yielding 247,032 gRNA combinations (Figure 3A). Sequencing confirmed their completeness and uniform distributions with AUC values of 0.65 and 0.7, and skews of 1.22 and 1.69 for single and multiplex libraries, respectively (Supplementary Figures S4A–C, S5 and Table S5).
Single and combinatorial screens uncover autophagy GIs in cell proliferation. (A) Autophagy multiplex library design. Combinatorial gRNA constructs target two autophagy genes simultaneously. The autophagy multiplex library consists of two gRNA-expression cassettes; the extended autophagy side (h7SK) with 876 gRNAs targeting 192 genes (four gRNAs per gene), and the core autophagy side (hU6) with 282 gRNAs targeting 64 genes (four gRNAs per gene), resulting in 247,023 gRNA combinations. The dedicated single autophagy library contains the extended autophagy gRNAs on the h7SK cassette, while the hU6 cassette is wild type (I-Sce1 placeholder). (B) Autophagy single and multiplex screen setup for proliferation and autophagy flux readouts. Single (s) and multiplex (mpx) autophagy libraries were applied in replicates with 1,000-fold (1,000×) or 20-fold (20×) coverages, respectively. RPE1 autophagy reporter (GFP-LC3-RFP) cells were transduced with an MOI of 0.5 and cultivated for 12 days. Proliferation readout samples were harvested on day 12 (pre-FACS). Autophagy flux read out cells were treated with Torin1 at day 12 for 24 hours and subsequently cell-sorted to enrich cells with blocked autophagy (post-FACS). Harvested cells were processed for sequencing. (C) Volcano plots of MAGeCK-derived log2-fold changes (LFCs) and P-values for dedicated and multiplex-inherent single autophagy proliferation screens. False-discovery-rate (FDR) for positive and negative selections are color-coded yellow-red and yellow-blue, respectively. Significant (P< 0.05) data points with LFC > 1 or LFC < -1 have dashed strokes. Data point sizes indicate the number of gRNA used for MAGeCK hit calling. (D) Volcano plots of MAGeCK-derived LFCs and P-values for the multiplex autophagy proliferation screen. FDR for positive and negative selections are color-coded yellow-red and yellow-blue, respectively. Significant (P< 0.05) data points with LFC > 1 or LFC < -1 have dashed strokes. Data point sizes indicate the number of gRNA used for MAGeCK hit calling. (E, F) Scatter plot of SUM-model derived expected LFC and observed LFC from multiplex autophagy proliferation screen. All data points with LFC < -1 (E) or LFC < 1 (F) are shown. Highlighted are data points with delta log2 fold-changes (dLFC) < -2*0.625 (E, red) and dLFC > 2*0.6655 (F, blue). (G) Bar graph of the extreme negative (red, depleted) and positive (blue, enriched) observed LFC for GIs from (E) and (F), compared to their single gene LFCs derived from the multiplex-inherent single screen. Dotted lines indicate strongest positive (blue) and negative (red) single gene deletion phenotypes of GIs, with WASHC1 and TP53 serving as references. (H) Number of interaction partners for genes of depleted (E, red) and enriched (F, blue) gene pairs. (I) Kaplan-Meier curves of five-year overall survival (OS) rates of lung squamous cell carcinoma patients based on ATG7, KEAP1 and KEAP1 and ATG7 gene expression levels. Hazard ratios (HR). (J) TCGA analysis of lung squamous cell carcinoma patient survival dependent on ATG7-KEAP1 gene expression profile using the UCSC Xena online tool. Regression line in red color.
Autophagy gene interactions enhancing/suppressing cell proliferation
To quantify proliferative and autophagy flux phenotypes by single and multiplex CRISPR in the same genetic background, we transduced RPE1(Cas9) cells with a previously described autophagic flux reporter and generated a monoclonal cell line (Supplementary Figure S4D,E) (10). The fluorescent reporter consists of GFP-LC3 fused to the N-terminus of RFP, and its expression leads to the near-equimolar presence of GFP-LC3 and RFP through cleavage by cellular ATG4 proteases. Upon autophagy activation, GFP-LC3 is lipidated and degraded by autophagy, while RFP remains in the cytosol and can serve as an internal control (Supplementary Figure S4D), enabling the quantification of autophagic flux by determining the cellular GFP/RFP signal ratio.
Using the reporter cell line, we initially focused on identifying proliferative effects induced by single or combinatorial autophagy gene depletion. To do so, we transduced the dedicated single (extended autophagy) and multiplex autophagy libraries with representations of 1,000- and 20-fold, respectively, and determined their respective gRNA abundance after 14 cell divisions (Figure 3B). The multiplex autophagy library contained all extended autophagy sequences paired with NHTs (inherent single) (Figure 3A). We therefore correlated the dedicated and inherent single screens on gRNA and gene level. As expected, essential genes for cell proliferation were depleted and NHTs were slightly enriched with an overall correlation of 0.97 (Supplementary Figures S4F and S6A,C), further supporting our conclusion that CRISPR libraries with distribution skews of less than 2 could be screened with reduced coverage. Next, we applied MAGeCK to the dedicated and inherent single screen data sets and identified WASHC1 (also known as FAM39E or WASH1) to be essential for RPE1 cell proliferation (Figure 3C and Supplementary Table S5). In addition to its function as a nucleation-promoting factor at the surface of endosomes (59), WASHC1 is thought to negatively regulate autophagy by inhibiting BECN1 ubiquitination to inactivate PIK3C3/VPS34 activity (60). However, we note that WASHC1 was likely an essential gene for RPE1 cell proliferation that lacked previous recognition due to the absence of WASHC1 gRNAs from previous genome-scale CRISPR libraries (1,2). Moreover, we observed no negative proliferative effect when core essential autophagy genes were depleted, suggesting that RPE1 proliferation was independent of autophagy during the course of the experiment (Figure 3C).
Next, we focused on gene pairs that were unknown to cause proliferative effects. A MAGeCK analysis resulted in a total of 446 and 732 significantly depleted and enriched gene pairs, respectively, of which WDR45B-PIK3R4 and ATG7-KEAP1 showed the strongest effects (Figure 3D and Supplementary Table S5).
A GI occurs when the observed combinatorial knockout phenotype deviates from a predicted phenotype that is expected for non-interacting genes. Importantly, for WDR45B-PIK3R4 and ATG7-KEAP1, the combinatorial phenotype could not be explained by summing up the effects of their individual knockout phenotypes (Figure 3E–G). Therefore, our screen identified WDR45B-PIK3R4 and ATG7-KEAP1 as GIs.
To identify all significantly depleted and enriched gene pairs that constitute GIs, we computed expected combinatorial knockout phenotypes with the SUM model. In total, we observed GIs with more than two standard deviations of the dLFC for 35.4% of the depleted and 52.3% of the enriched gene pairs (Figure 3E-F). While the decreased proliferation of WDR45B-PIK3R4 was milder compared to the essential gene WASHC1, the increased proliferation of ATG7-KEAP1 depletion was stronger than that of TP53 depletion (Figure 3G). Lastly, we analyzed the number of interaction partners for each gene and identified PIK3R4, ATG7, and ATG5 to have 172, 174 and 157 interaction partners, respectively. These three genes had the highest number of connections compared to all other combinations, identifying them as hub genes for proliferation within our tested autophagy gene set (Figure 3H).
Previous meta-analysis demonstrated that lung tumors with increased proliferative rates were associated with worse prognosis and reduced survival (61). Considering that co-deletion of ATG7 and KEAP1 enhanced proliferation in RPE1 cells (Figure 3D,F), we investigated their role as survival biomarkers. Therefore, we analyzed GEO and EGA datasets and generated Kaplan-Meier curves to correlate the five-year overall survival probability of lung squamous cell carcinoma (LUSC) patients based on KEAP1, ATG7 or KEAP1 and ATG7 gene expression data. Notably, the simultaneous downregulation of both genes significantly correlated with reduced overall survival (Figure 3I). To corroborate our finding, we analyzed ATG7 and KEAP1 gene expression and LUSC patient overall survival obtained from the Pan-Cancer Atlas consortium (TCGA). LUSC patients were classified in four different categories depending on ATG7 and KEAP1 gene expression levels. Interestingly, patients with low expression of ATG7 and KEAP1 had the lowest overall survival (Figure 3I-J). These analyses suggest a potential clinical relevance of ATG7 and KEAP1 as a novel LUSC biomarker that, however, requires further investigation.
Dedicated and inherent single screens identify genes essential for autophagy flux
To retrieve GIs in autophagy flux, we first identified all single gene phenotypes. To do so, we applied the extended autophagy library to the LC3 reporter cell line, coupled to fluorescence-activated cell sorting (Figure 3B). Cells with blocked autophagy were identified and sorted due to their unchanged GFP abundance using a stringent gating criterion to reduce false-positives (Figure 4A and Supplementary Figure S4E). Sequencing of the sorted cells revealed high gRNA and gene level correlations between post-FACS samples among biological replicates (Supplementary Figures S5 and S6A-B). We then applied MAGeCK and retrieved 12 significantly enriched genes with well-established functions in autophagy flux and LFCs ranging from 3.76 to 9.12 (Figure 4B and Supplementary Table S5) (11,14,15,20,21,62,63). Next, we applied the multiplex autophagy library and noticed an additional population of cells with increased RFP levels (high-gate) during FACS (Figure 4C). We identified the majority of this population to account for ATG4B-containing gRNA pairs (post-FACS 1, 86.94%; post-FACS 2, 78.03%; post-FACS 3, 75.53%) (Supplementary Figure S6A, D, E). Since the autophagy reporter requires processing by protease ATG4, we concluded that targeting ATG4B interfered with its functionality, leading to the appearance of the high-gate cell population. Experimental correlation among the biological replicates of the multiplex screens were high on gRNA and gene level (Supplementary Figure S6A,C). We therefore continued to perform MAGeCK analyses on the inherent single autophagy flux data and identified 10 significantly enriched genes (Figure 4D and Supplementary Table S5). Moreover, in the dedicated single screen, we identified TMEM41B, an autophagy regulator that was only recently described to be essential for bulk autophagy (14,15). Despite the high concordance among dedicated and inherent single autophagy flux screens with seven shared hit genes, both screens identified additional eight complementary hits (Figure 4B,D). The concordance could be improved by FDR-relaxation, resulting in 11 and 15 shared and complementary hits, respectively (Supplementary Figure S7). Shared and complementary hits were experimentally validated with screen-independent gRNAs in RPE1 cells and in the acute monocytic leukemia cell line THP1, with gRNA performance being quantified by TIDE analysis (Figure 4E, G). Interestingly, while the depletion of these genes had no effect on RPE1 cell fitness, the percentage of autophagy flux inhibition correlated well with negative cell fitness in THP1 cells (Figure 4F, G), suggesting that THP1 cells depend on basal autophagy flux for viability. Together, our dedicated and inherent single autophagy flux screens correctly identified genes essential for autophagy flux.
Dedicated and multiplex-inherent single gRNA enrichment screens identify essential autophagy genes. (A) FACS analysis of dedicated single autophagy flux screen, highlighting untreated (red) and Torin1-treated (blue) RPE1 autophagy reporter cells. Sorting gate (#, post-FACS) is shown as a black ellipse at the very boundary of the untreated cell population. (B) Positively selected hits from a MAGeCK analysis comparing pre- and post-FACS samples of the dedicated single gRNA autophagy screen. Hit genes with a false-discovery rate (FDR) below 10% are highlighted in red. Hits that were not identified in the multiplex-inherent single screen are highlighted in blue. (C) FACS analysis of multiplex-inherent single autophagy flux screen, highlighting untreated (red) and Torin1-treated (blue) RPE1 autophagy reporter cells. Sorting gates (#, post-FACS; *, high-gate) are shown as a black ellipse. (D) Positively selected hits from a MAGeCK analysis comparing pre- and post-FACS samples of the multiplex-inherent single gRNA autophagy screen. Hit genes with an FDR below 10% are highlighted in red. Hits that were not identified in the dedicated single screen are highlighted in blue. (E) Autophagy flux validation of selected single hit genes derived from (B) and (D) are shown in red. Evaluation of gRNA activity by TIDE analysis (grey). Error bars represent standard error of mean (SEM) over three biological replicates per autophagy blockage (n = 3). ND: not determined. (F) LFCs of dedicated single autophagy flux and of multiplex-inherent single autophagy proliferation screens do not correlate (r = -0.05). Highlighted are NHTs (blue), single essential genes for autophagy flux (red) and proliferation (yellow, according to Hart et al., 2017). (G) Autophagy block (GFP/RFP to NHT) and proliferation of single essential autophagy gene depletions in THP1 cells correlate (r = -0.8). Highlighted are NHT (blue) and single essential genes for autophagy flux (red).
MAGeCK identifies genetic interactions essential for autophagy flux
To identify enriched gene pairs in the combinatorial screen after FACS, we applied MAGeCK to the entire dataset of 247,032 gRNA combinations and analyzed the results on gene level. We noticed that the majority of the 1,910 gene combinations with a significant positive LFC (p < 0.05) were actually combinations for which one or both genes were essential for autophagy flux, as determined from our single screen data. We identified 170 gene pairs for which both partners alone were essential for autophagy flux and 1,731 gene pairs for which only one partner was essential for autophagy flux (Figure 5A). Next, we were interested in finding gene combinations of essential and non-essential genes. We removed all gene combinations with a LFC < -0.5 and for which both genes were essential and repeated the MAGeCK analysis on the resulting data set. To our surprise, all significantly positive 639 combinations (p < 0.05) were pairs of essential with non-essential genes (Figure 5B). The FACS-enrichment of gene-pairs did not identify them as GIs, because their combinatorial phenotypes might be driven by the strong phenotype of a single essential autophagy gene. We aimed at identifying GIs that blocked autophagy flux, while their corresponding single gene knockouts did not block autophagy on their own. We hypothesized that we could utilize MAGeCK to identify GIs by limiting our dataset to combinations of non-essential genes with a LFC > -0.5 and performed a third MAGeCK analysis. This analysis identified 74 GIs consisting of two non-essential autophagy genes, of which the most significant were: AMBRA1 paired with either TRIM5, WDR45B or PEX13; ULK1 paired with either STX17 or ATG16L2; and ATG2A-ATG2B (Figure 5C). In arrayed validations, we confirmed combinations of two single essential autophagy genes, combinations with one single essential autophagy gene, and combinations of only non-single essential autophagy genes to block autophagy. Targeting combinations of two single essential genes consistently increased the block of autophagy flux compared to targeting one single essential gene, an effect particularly prominent for ATG5-ATG14 (Figure 5D). Similarly, validations of significantly enriched gene pairs consisting of essential and non-essential genes also increased the fraction of cells with blocked autophagic flux (Figure 5E). Most importantly, the arrayed validation confirmed the ATG2A-ATG2B interaction to be essential for autophagy flux in RPE1 and THP1 cells (Figure 5F). Together, our MAGeCK analyses successfully identified essential GIs for autophagy flux consisting of non-essential autophagy genes.
MAGeCK identifies genetic interactions (GIs) in autophagy flux. (A–C) Iterative MAGeCK analysis of the multiplex autophagy flux screen on the entire dataset (A), the entire dataset minus essential:essential combinations (B), and only on non-essential:non-essential combinations (C). Visualized are MAGeCK-derived log2 fold changes (LFCs) and p-values. (A) The majority of essential:essential (ess-ess) combinations are significantly positive (pos.sign., P < 0.05) (170 of 195, dark red). Combinations of significantly positive non-essential:essential genes (pos. sign. ess) are highlighted in dark blue. (B) All significantly positive hits (639 of 639, dark blue) are essential:non-essential combinations. (C) Highlighted are putative GIs of non-essential:non-essential combinations (74, orange), and ATG2 hit combinations (4, green). (D–F) Arrayed validation of block in autophagy flux of single genes and hit gene pair depletions derived from (A-C). Single essential genes highlighted in bold. Error bars represent standard error of mean (SEM) over three biological replicates (n = 3). (D) Co-depletion of two single essential autophagy genes enhances the block in autophagy flux. Single gene knockouts in gray color, and double gene knockouts in red color. (E) Co-depletion of an essential and a non-essential autophagy gene enhances the block in autophagy flux. Control gene knockouts in grey color, and double gene knockouts in blue color. (F) Co-depletion of two non-essential autophagy genes enhances the block in autophagy flux in RPE1 and THP1 cells. Control gene knockouts in grey color, and double gene knockouts in orange color.
MAX model-derived genetic interactions essential for autophagy flux
Four models have been proposed to identify GIs: Multiplicative (MULT), Additive (SUM), Logarithmic (LOG) and Minimal (MIN) (36). In the MIN model, the expected phenotype for combinatorial depletion of non-interacting genes corresponds to the single mutant with the most severe phenotype, while the SUM model defines the expected phenotype as the sum of the measured phenotypic strength of both individual mutants. In the MULT model, the expected phenotype is predicted as the product of both single mutants and the LOG model determines GIs from measurements on a logarithmic fitness scale. However, because the knockout of a second autophagy gene is unlikely to ameliorate the phenotype of the first autophagy gene, we added the maximum (MAX) model that identified a GI between two genes when the observed combinatorial phenotype is stronger than the phenotypes of both single mutants. Each model assigned a GI to a gene pair xy, if the double mutant phenotype W(xy), deviated from a predicted double mutant phenotype that was expected for no interaction between x and y, E(xy). The phenotype was measured as the LFC of gRNA abundance between respective samples, and the expected double mutant phenotype for joint mutations of the genes x and y was defined by the neutrality function of each definition with MIN = min(W(x), W(y)), MULT = W(x) × W(y), LOG = log2[(2W(x) + 1) × (2W(y) – 1) +1], SUM = W(x) + W(y) and MAX = max(W(x), W(y)). We applied all five models to compute expected phenotypes for non-interacting gene pairs based on the single gene effects from the single screens. Then we computed the dLFC, the deviation of observed from expected phenotypes.
Previous screens in Saccharomyces cerevisiae and human cells revealed GIs to be rare (64,65). In line with this, the MAX model detected the fewest deviations of observed combinatorial phenotypes from the predicted value for non-interacting genes (Figure 6A). First, we focused on non-essential gene combinations that were enriched after FACS (observed LFC > 0) and had a positive dLFC according to the MAX model (dLFC > 0), and identified a total of 255 GIs (Figure 6B, blue: 217, orange: 38; and Supplementary Table S6). Next, we analyzed the complete data set of 1,570 MAX model-derived GIs (LFC > 0.5, dLFC > 0.05), including GIs containing single essential autophagy genes. We computed the number of GIs per core autophagy gene and found 11 core essential autophagy genes to be highly connected with an average of 241 interactions (of possible 255) (Figure 6C). The phenotypic strength of autophagy genes correlated well with the number of interactions (Figure 6C). In addition, we identified AMBRA1, TMEM41B, EPG5 and ULK1 as hub genes to interact with 84 to 25 other core and extended autophagy genes (Figure 6C). One essential GI identified by the MAGeCK analysis was the paralog combination ATG2A-ATG2B. Paralogs have been reported to be less frequently identified as essential than expected in monogenic screens (26,27). Assuming a capacity to functionally compensate for each other's loss, it was proposed that combinatorial screens may be able to identify novel gene functions or essentiality for paralog pairs (26,27). To analyze the compensatory capacities of all core autophagy paralogs, we examined their GI profiles in more detail by evaluating their dLFCs across all core autophagy genes (Figure 6D). As expected from our previous analysis, we identified a GI between the ATG2A-ATG2B paralog pair, indicating that the presence of at least one of the paralogs is required for autophagy flux. In addition, within the WIPI gene family, we identified multiple GIs, namely WIPI1-WIPI2, WIPI2-WDR45 and WDR45-WDR45B (Figure 6D). The presence of these GIs suggests specific roles of the WIPI members that cannot be compensated for by other members of that family, as was also previously reported (66–68). In contrast, we observed no GIs within the ULK and ATG4 paralog families, indicating possible buffering capacities of the corresponding three paralogs for each gene. The GIs between GABARAP, MAP1LC3B and GABARAPL2 and their observed higher transcript levels suggested specific roles of these three ATG8s in RPE1 cells compared to the ATG8 paralogs that showed no GIs (Figure 6D). Additionally, we observed very similar GI profiles for the two essential core autophagy paralogs ATG9A and ATG16L1 and only a few weak interactions for their corresponding paralogs ATG16L2 and ATG9B (Figure 6D). To account for context-dependent functions of ATG9B and ATG16L2, we quantified their transcript levels in RPE1 cells by RNA-seq. As expected, ATG9A transcripts were highly abundant, while ATG9B transcripts were undetectable (Figure 6E). In accordance to the distinct tissue-specific roles that have been reported for ATG16L2 and ATG16L1 (69), transcripts levels of ATG16L2 were much reduced compared to ATG16L1 (Figure 6E and Supplementary Table S7). Moreover, interacting members of the ATG8 family GABARAP, GABARAPL2 and MAP1LC3B showed higher transcript levels than non-interacting members of the ATG8 family, suggesting a context-specific dependency of GIs (Figure 6E).
MAX model analysis identifies genetic interactions (GIs) in autophagy flux. (A) Density plot of delta log2 fold change (dLFC) distributions of all gRNA combinations in the autophagy multiplex post-FACS screen according to the MIN, MAX, MULT, SUM, LOG genetic interaction models. (B) MAX model-derived comparison of observed and expected combinatorial phenotypes. Highlighted are combinations with a LFC > 0 and a dLFC > 0 (blue, 217), and with a standard deviation (SD) of the delta LFC (dLFC) above 3 (orange, 38). Combinations with single core essentials for autophagy, viability (Hart et al., 2017), and NHTs are not shown. (C) Correlation of the number of interactions and the LFC of single autophagy genes (r = 0.87). Highlighted in red are single essential autophagy genes. (D) dLFC heat map of averaged AB–BA GIs of core autophagy paralog genes across all core autophagy genes. Core autophagy genes are grouped according to their function in autophagy. Red squares highlight GIs of paralog family members. dLFC values are color-coded white to blue. (E) Mean of GI dLFCs of core autophagy paralogs across all core autophagy genes (blue squares, left Y-axis) compared to their normalized transcript abundance in RPE1 cells (red dots, right Y-axis). Error bars of mean dLFCs represent standard error of mean (SEM) over all GIs. Error bars of normalized transcript abundance represent SEM over three biological replicates (n = 3). Not determined (nd). (F) Network analysis of high-confidence GIs with dLFC > 3*SD and LFC > 0 (orange data points from Figure 6B). Genes and GIs are grouped based on their role in processes, complexes or paralog families. Names of paralog GIs are highlighted in pink. Edge color is set according to GI dLFCs (white to blue). (G) Arrayed validation of block in autophagy flux of control and single genes (grey) and hit gene pairs (orange) depletions derived from (F). dLFC of hit gene pairs, derived from (B), is shown blue. Error bars represent standard error of mean (SEM) over three biological replicates (n = 3).
Lastly, we selected high-confidence synergistic GIs of non-essential autophagy genes with a SD(dLFC) > 3 (Figure 6B, orange). For visualization purposes, we integrated these GIs in a network and identified several to be connected to selective autophagy (NBR1, PEX13, FUNDC1, FUNDC2, PLEKHM1, FAM134B and CCPG1). Among these were three GIs linking ULK4 to the selective autophagy receptors PEX13, FUNDC1 and SQSTM1, and two GIs connecting IRGM to the ER-phagy pathway genes CCPG1 and FAM134B (Figure 6F). Furthermore, we identified a very strong interaction between the transcription factor TFEB and WAC, which is required for starvation-induced activation of ULK kinase (70). To validate our approach and analysis, we generated single and combinatorial gene depletions of selected high-confidence GIs in the RPE1 autophagy reporter cells, induced autophagy by Torin1 treatment and quantified the block of autophagy by FACS. In cells depleted of a single gene, autophagy was largely unaffected, with the exception of the WIPI2 gene that reduced autophagy to 25.0% (Figure 6G). Interestingly, we previously identified WIPI2 as a complementary hit of the dedicated single screen under conditions of FDR-relaxation. This result, thus, validated the approach of FDR-relaxation and identified WIPI2 as a single essential gene for autophagy. In contrast to the single gene depletions, all tested gene combinations blocked autophagy with efficiencies that ranged from 14.6% (SCOC-SIRT2) to 45.6% (WIPI2-GABARAPL2), irrespectively of their dLFC (Figure 6G), validating them as true GIs in autophagy. Together, our analysis and validation identified hitherto unknown GIs in autophagy.
DISCUSSION
Several technologies are available to generate combinatorial CRISPR plasmids and libraries. However, most rely on iterative and pooled cloning of PCR-amplified oligonucleotide pools, contributing to cloning errors and library distribution bias. Our 3Cs multiplexing technology works with single-stranded template plasmids and non-amplified oligonucleotide pools, thereby circumventing traditional cloning and unintended sequence bias to provide a scalable method for generating multiplex CRISPR libraries of arbitrary size without affecting library distribution. Screening coverage and CRISPR library distribution skew have recently been computationally predicted to be critical factors for data quality (3). Indeed, our experimental analysis of combinatorial libraries with varying distribution skews and screening coverages, supports this prediction. We show that distribution skews of 2 or below enable library representations of 20-fold, downsizing current screen sizes, associated efforts and costs by a factor of at least 10-fold. Furthermore, we provide experimental evidence that combinatorial libraries with distribution skews > 2 should be screened with coverages ≥ 200 to compensate for the library distribution skew. However, rescuing a large distribution skew (≥ 10) requires coverages above 200-fold. The impact of the library distribution skew on the required experimental coverage therefore highlights the importance of robust and scalable methods for generating uniformly distributed gRNA libraries.
We demonstrate the feasibility of screening uniform gRNA libraries with minimized coverages by applying our multiplex autophagy library to identify GIs of autophagy genes in cell proliferation and autophagy flux. We identified PIK3R4, ATG5 and ATG7 as hub genes. The presence of hub genes in gene networks that modify the phenotypic consequences of mutations in other genes, has been reported previously and is of clinical relevance (71–73). In line with this, we observed that PIK3R4, ATG5 and ATG7 modified the proliferative phenotypes of more than 80% of all tested genes. Although we identified ATG5 and ATG7 to be connected best in RPE1 proliferation, their single depletion has no effect on RPE1 proliferation but blocks autophagy flux. Thus, we speculate that ATG5 and ATG7 have autophagy-independent functions that, when co-depleted with other autophagy genes, culminate in cell growth phenotypes, a finding also supported by the literature (74). Interestingly, low expression levels of ATG7-KEAP1 correlate with decreased overall survival in lung squamous cell carcinoma (Figure 3I, J), suggesting a clinical relevance of this synergistic gene interaction.
Apart from proliferation, we identified all single essential core autophagy genes as hub genes in autophagy flux, interacting with more than 60% of all tested genes. In addition, we found AMBRA1 to be highly connected (84 interactions, 21% of all tested genes), as the only gene that our and previous monogenic CRISPR screens did not identify as single essential for autophagy. However, the high connectivity of AMBRA1 with other autophagy genes in our combinatorial CRISPR screen for autophagy flux suggests an important buffering role of AMBRA1 in the context of Torin1-induced autophagy that awaits its mechanistic exploration. Moreover, the clear separation of core essential autophagy genes with regard to their connectivity, including TMEM41B that was only recently described to be essential for bulk autophagy (14,15), points to the ability of multiplex screens to reveal single but mild gene phenotypes, such as for TMEM41B, EPG5 or ULK1. We reason that this was because the large number of gRNA combinations in the autophagy multiplex screen provided more data points for each core autophagy gene, when compared to the dedicated single screens. The dedicated single autophagy library contained only four gRNAs per gene, whereas the autophagy multiplex library contained 4 * 282 gRNAs + 4 * 876 gRNAs = 4,632 gRNA combinations per gene. While dedicated and inherent single autophagy flux screens only jointly identified all core essential autophagy genes, the combinatorial autophagy flux screens successfully identified all core essentials on the base of gene connectivity. Thus, the 3Cs multiplex autophagy flux screen had a higher phenotypic sensitivity than the previous monogenic CRISPR screens, and confirms the ability of combinatorial screening to uncover single gene function on the base of gene connectivity.
Despite the high gRNA and gene correlation of dedicated and inherent single autophagy screens, we identified complementary hits in both. Of note, dedicated and inherent single screens had 8 mutually exclusive hits, of which 6 could be experimentally validated. This suggests that although both screens jointly identified autophagy-blocking genes, the dedicated single screen had a higher phenotypic resolution to extract hit genes with milder phenotypes, e.g. TMEM41B and WIPI2, for which the block of autophagy could be validated. In line with this, experimental coverages of dedicated and inherent screens were 1,000- and 20-fold, respectively, and it cannot be excluded that a higher experimental coverage of the multiplex autophagy library would have resulted in an increased phenotypic resolution. However, an increase in library coverage, particularly for FACS-based reporter enrichment applications requires a drastic increase in cell culture and FACS-time. To overcome this issue, we performed a dedicated single screen with a higher coverage than the multiplex screen, thereby accounting for genes with milder phenotypes, and enabling a more stringent cutoff for combinatorial screen analysis and identification of genetic interactions.
Accordingly, genome-wide CRISPR screens demand large numbers of cells, which is why current efforts aim to minimize cell culture demands by providing minimized combinatorial gRNA libraries (56,75). We identified the LFC of dual gRNA-targeted genes to be larger than their single gRNA-targeted counterparts (Supplementary Figures S2E-F, S8A-C), suggesting dual-gRNA targeting to minimize libraries for monogenic screens by simultaneously improving the phenotypic resolution. Paralogs are less frequently identified as hits in monogenic CRISPR screens (26,27). However, for the majority of autophagy paralogs, redundant or specific functions are currently elusive. Of the 82 possible core autophagy paralog GIs in our screen, we only identified ATG2A-ATG2B, GABARAP-GABARAPL2, and GABARAP-MAP1LC3B, confirming the essential interaction of ATG2A with ATG2B for autophagy flux (76). Moreover, our analysis suggests that the loss of GABARAP-GABARAPL2 and GABARAP-MAP1LC3B can not be compensated for by the other ATG8 family members. Interestingly, we did not observe an interaction between GABARAPL2 and MAP1LC3B, suggesting these genes to have non-overlapping functions, which, however, can both be compensated for by GABARAP. Since many core autophagy paralog families include more than two paralogs (ATG8s, ULKs, ATG4s), we reason that redundancy and buffering functions may only be uncovered by higher-order GI screens (77,78). In line with this, we present 3Cs multiplexing for the generation of combinatorial dual-gRNA CRISPR libraries. However, it may be possible to extend the 3Cs technology to higher-order (>2 sgRNAs) combinatorial CRISPR libraries. The availability of different RNA polymerase III promoters from multiple species (U6, 7SK, or H1), together with the existence of engineered gRNA scaffolds for SpCas9 (79), provides sufficient design space to potentially accommodate up to 5 non-recombining gRNA-expression cassettes with unique 3Cs homology on a single lentiviral plasmid. While we do not foresee limitations in the generation of higher-order heteroduplex 3Cs-DNA, coverage-based amplification and high-throughput sequencing impose practical constraints. Bacterial amplification requires a transformation efficiency that outnumbers library diversity by at least 100-fold, and considering a stable transformation efficiency of 1010, the gRNA diversity per gRNA-expression cassette is limited to 10,000, 500, 100, or 50, for libraries with 2, 3, 4, or 5 cassettes, respectively. However, due to the additive nature of transformation efficiencies, the increase in gRNA-expression cassettes and its associated decline in gRNA diversity can be accounted for by paralleling 3Cs reactions and their amplification. Still, their experimental application and analysis is limited by cell culture demands and high-throughput sequencing strategies. 3Cs libraries enable minimized CRISPR screens, but the simultaneous sequencing of up to 5 sgRNAs is currently not established with classical sequencing-by-synthesis approaches. This, however, generates innovative engineering space to integrate long-read sequencing technologies into the analysis of higher-order combinatorial CRISPR libraries. Lastly, it is important to emphasize that, providing sufficient homology for unique oligonucleotide annealing, 3Cs multiplexing is not limited to SpCas9 gRNAs and is likely equally well suited for the generation of gRNA libraries for other Cas nucleases for related and orthogonal applications, as well as the identification of genetic interactions at scale.
DATA AVAILABILITY
Sequencing data are provided as raw read count tables as Supplementary Tables S2-S5. RNAseq data are accessible from GEO (GSE166851). Plasmids encoding extended- and combinatorial autophagy libraries can be requested through the Goethe University
Depository (https://innovectis.de/en/technologies/goethe-depository/, please select ‘Minimized combinatorial CRISPR screens identify genetic interactions in autophagy’) or directly contact the Corresponding Author. Custom scripts are available from GitHub (https://github.com/ManuelKaulich/3Cs-multiplexing).
SUPPLEMENTARY DATA
Supplementary Data are available at NAR Online.
ACKNOWLEDGEMENTS
We thank Andreas Ernst and Svenja Wiechmann for technical advice, Andrea Ballabio and Davide Cacchiarelli of the Telethon Institute of Genetics and Medicine (Naples, Italy), Sebastian Wagner and Khalil Abou Elardat (Cancer Genomics Core Facility Frankfurt, Germany) and Tobias Schmidt (Institute of Biochemistry I (Pathobiochemistry), Frankfurt, Germany) for valuable support with Illumina high-throughput sequencing, and Andrew Holland for providing puromycin-sensitive hTERT-RPE1 cells. lentiCRISPRv2-puro and lenti-sgRNA-blast were a gift from Brett Stringer (Addgene plasmid # 98290; http://n2t.net/addgene:98290; RRID:Addgene_98290; Addgene plasmid # 104993; http://n2t.net/addgene:104993; RRID:Addgene_104993). pKLV2.2-h7SKgRNA5(SapI)-hU6gRNA5(BbsI)-PGKpuroBFP-W was a gift from Kosuke Yusa (Addgene plasmid # 72666; http://n2t.net/addgene:72666; RRID:Addgene_72666). pMRX-IP-GFP-LC3-RFP was a gift from Noboru Mizushima (Addgene plasmid # 84573; http://n2t.net/addgene:84573; RRID:Addgene_84573). psPAX2 and pMD2.G were a gift from Didier Trono (Addgene plasmid # 12260; http://n2t.net/addgene:12260; RRID:Addgene_12260; Addgene plasmid # 12259; http://n2t.net/addgene:12259; RRID:Addgene_12259). Funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.
FUNDING
Hessisches Ministerium für Wissenschaft und Kunst (HMWK) [LOEWE-CGT IIIL5-518/17.004, LOEWE-FCI IIIL5-519/03/03.001]; Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) [CEF-MC-EXC115/2, ECCPS-EXC147/2, CPI-EXC 2026, 259130777 – SFB 1177]. Funding for open access charge: DFG [SFB 1177 - 259130777]
Conflict of interest statement. The Goethe University Frankfurt has filed a patent application related to this work, on which Valentina Diehl, Martin Wegner, Ivan Dikic and Manuel Kaulich are inventors (WO2017EP84625). The Goethe University provides an exclusive license of the 3Cs technology to Vivlion GmbH, for which Ivan Dikic and Manuel Kaulich are co-founders, shareholders and chief officers.
Notes
Present address: Paolo Grumati, Telethon Institute of Genetics and Medicine, Via Campi Flegrei 34, 80078 Pozzuoli (NA), Italy.
REFERENCES
Author notes
The authors wish it to be known that, in their opinion, the first three authors should be regarded as Joint First Authors.
Comments