High-throughput identification of genes involved in single and multi-drug resistance in pancreatic cancer with pooled CRISPR screening

SUMMARY Chemoresistance is the norm in pancreatic ductal adenocarcinoma (PDAC) leading to an abysmal five-year survival rate of less than 10%. We performed genome-wide CRISPR activation (CRISPRa) and CRISPR knock out (CRISPRko) screens to identify mechanisms of resistance to four cytotoxic chemotherapies (gemcitabine, 5-fluorouracil, irinotecan, and oxaliplatin) used to treat PDAC patients in two PDAC cell lines. Mirroring patient treatment outcomes, we found that multi-drug resistance genes were more likely to be associated with patient survival. We identified activation of ABCG2, a well-described efflux pump, as the most consistent resistance gene in four drugs and two cell lines. Small molecule inhibitors of ABCG2 restored sensitivity to chemotherapy. Our screen demonstrated that activation of members of the hemidesmisome complex and transcriptional repressor complexes conferred resistance to multiple drugs. Finally, we describe an approach for applying our results to predict drug sensitivity in PDAC tumors and cell lines based on gene expression.


INTRODUCTION
The era of high-throughput genome sequencing has led to an unprecedented understanding of the genetic and epigenetic alterations underlying the most common cancers (Kandoth et al., 2013;Sanchez-Vega et al., 2018). In several cases this has led to the development of successful targeted therapies. Early examples such as the BCR-ABL gene fusion, which can be targeted with the tyrosine kinase inhibitor imatinib (Goldman and Melo, 2003), and EGFR kinase domain mutations, which have been inhibited with gefitinib (Pao et al., 2004), have been followed by a rapidly growing number of promising targeted therapeutics (Bollag et al., 2012;Katzel et al., 2009;Warr and Shore, 2008). Unfortunately, pancreatic ductal adenocarcinoma (PDAC) has remained largely refractory to the improvement seen in several other cancers with five-year survival rates still less than 10% (Siegel et al., 2018). With largely non-specific symptoms and invasive procedures required for diagnosis, only 20% of PDAC patients are eligible for surgical resection, leaving a majority of patients with chemotherapy and radiation as their sole course of treatment (Kamisawa et al., 2016). Even patients eligible for surgical resection tend to have recurrent disease within 2 years that is refractory to adjuvant chemotherapy (Stathis and Moore, 2010). In an effort to improve patient outcomes, multi-drug combinations such as FOLFIRINOX (fluorouracil, folinic acid, irinotecan and oxaliplatin) are increasingly being used despite associated toxicities and modest improvements in outcomes (Conroy et al., 2011).
Large cohort sequencing efforts have established a relatively short list of commonly mutated genes, including KRAS, TP53, SMAD4, and CDKN2A, however these have thus far proven to be poor drug targets (Bailey et al., 2016;Raphael et al., 2017;Witkiewicz et al., 2015).
Moreover, few genetic perturbations have been reproducibly associated with patient prognosis (Shugang et al., 2016;Witkiewicz et al., 2015). Transcriptomic sequencing has been more successful identifying subtypes of PDACs that may predict prognosis (Bailey et al., 2016;Collisson et al., 2011;Moffitt et al., 2015;Raphael et al., 2017). Of particular interest, is the squamous/quasi-mesenchymal subtype, which consistently exhibits significantly poorer prognosis than other subtypes. Unfortunately, patient subtyping is based on the expression levels of thousands of genes in primary tumor tissue, making the delineation of individual targets difficult.
Complete resistance to potent multi-drug cocktails is common in PDAC patients (Garrido-Laguna and Hidalgo, 2015). Several extra-cellular mechanisms for drug resistance have been proposed including a highly immunosuppressive micro-environment (Bazhin et al., 2013), extensive desmoplasia (Chu et al., 2007), and hypovascularization (Koong et al., 2000). Tumor heterogeneity is another notable mechanism of resistance (Yachida and Iacobuzio-Donahue, 2009). PDAC tumors often consist of several unique clonal populations each of which may harbor a genetic background capable of conferring drug resistance and have the potential to become the dominant tumor population after drug treatment. However, these clones are almost impossible to detect in an agnostic manner a priori without extremely deep levels of sequencing that are cost prohibitive in most settings (Griffith et al., 2015). An alternative approach to this problem is to define the landscape of cellular mechanisms of PDAC drug resistance experimentally, then deeply screen tumors in a targeted manner for the presence of previously identified resistance drivers. Previous insertional mutagenesis-and RNA interference-based screens have successfully identified novel genes whose inactivation leads to gemcitabine sensitivity in PDAC cells (Bhattacharjee et al., 2014;Smith et al., 2014;You et al., 2011). More recently, genome-wide CRISPR-Cas9 screening, which provides complimentary information to previous screening methodologies (Evers et al., 2016;Morgens et al., 2016), for essential genes across RNF43-mutant and wild type PDAC cell lines found mutant lines to be specifically dependent upon FZD5 despite non-aberrant expression of FZD5 or other FZD receptors (Steinhart et al., 2017). We performed CRISPR-Cas9 knock-out (CRISPRko, Wang et al., 2014) and endogenous activation (CRISPRa, Konermann et al., 2015) screening in two PDAC cell lines (BxPC3 and Panc-1) to identify genes whose loss or gain of expression were able to modulate sensitivity to four of the most common cytotoxic chemotherapies used in the treatment of PDAC (gemcitabine, oxaliplatin, irinotecan, and 5-fluorouracil). We identified several novel genes and pathways that contribute to drug specific or multi-drug resistance. We validated our top hits in three cell lines using multiple sgRNAs. We determined that many of our hits are likely relevant to PDAC patients as their expression levels are correlated with overall survival. Finally, we used our screen results to develop a simple algorithm for predicting drug sensitivity in cell lines and patients.

CRISPR screening reveals drug resistance genes
We performed genome-wide CRISPRko and CRISPRa screening in the Panc-1 and BxPC3 cell lines to identify both drug-specific and multi-drug mechanisms of resistance to four drugs commonly used to treat PDAC (gemcitabine, oxaliplatin, irinotecan, and 5-fluorouracil) ( Figure 1A, Figure S1). Drug doses were optimized to inhibit growth to a level of 10-20% of the untreated cells although a broader range of selection strength was observed in some cases ( Figure   S2). We found the sgRNA log 2 fold changes (treated/untreated) of drug replicates were significantly more correlated than samples treated with different drugs ( Figure S3). However, there was a much higher correlation between samples treated with different drugs than expected by chance suggesting mechanisms of resistance were sometimes shared between drugs ( Figure   S3). The degree of replicate agreement also exhibited a strong inverse correlation with drug selection strength (Rho = -0.55, P=0.026, Figure S4). Clustering replicates by the most variable sgRNAs highlights that the strongest mechanisms of resistance are drug-specific, however the two nucleoside analogs, 5-fluorouracil and gemcitabine, share a relatively large number of top sgRNAs ( Figure 1B, Table S1). For example, over-expression of SLC29A1, a nucleoside transporter that is capable of both influx and efflux, confers significant resistance to gemcitabine and 5-FU, but little resistance to oxaliplatin or irinotecan. Other genes showing similar patterns include FMO5, a flavin-containing monooxygenase involved in metabolism xenobiotics that has been linked to prognosis in colorectal cancer (Zhang et al., 2018); several genes that function in cellular metabolism: SCD5, HYAC1, COQ6; and other genes involved in transport or its regulation, SLC22A9, SLC16A7, CCZ1, TBC1D3I.
Drug-specific mechanisms clustered into pathways, some of which were previously described. For example, TGFB1-dependent activation of the SMAD regulatory complex has been associated with oxaliplatin resistance (Mao et al., 2017;Sun et al., 2017) and hedgehog signaling has been linked to irinotecan resistance (Meng et al., 2015;Tripathi et al., 2014). Biosynthesis of heme and several phospholipid molecules were pathways uniquely associated with resistance to our nucleoside analogs. Both of these metabolic pathways have been described as important mechanisms of gemcitabine resistance and likely play a role in oxidative stress response (Miyake et al., 2010;Tadros et al., 2017). Overall, we saw that there was a weak positive correlation between enrichment across our screens and patient survival. We determined that this correlation was driven primarily by associations between survival and enrichment in multi-drug resistance genes. We observed that genes enriched in individual drug screens tended to be depleted for significant associations with PDAC patient survival ( Figure S5). Based on this result and the fact that PDAC patients are routinely treated with complex drug cocktails, we focused the remaining analysis on multi-drug resistance mechanisms under the assumption that they may hold greater clinical relevance.
To prioritize resistance genes, we computed the sum of the replicate minimum log 2 fold change of the two most enriched sgRNAs targeting each gene (L2FC sum) in each cell line (Table S2A-B). This approach requires that prioritized genes have a large increase in representation for two independent sgRNAs in each replicate and produces a relatively conservative list of hits. Multi-drug resistance genes were identified by computing the mean L2FC sum across all four drugs in each cell line. This approach was particularly powerful as it leveraged information from 48 experiments (4 drug screens x 3 replicates x 2 cell lines x 2 sgRNAs). We found L2FC sums to be weakly correlated across cell lines (rho=0.014, P=0.034), likely a reflection of strong phenotypic and genetic differences found between our two cell lines (Kim et al., 2014).
ABCG2 over-expression confers drug resistance A subset of genes possessed high L2FC sums in both cell lines. We performed validation experiments on several genes nominated by our screening data in three cell lines, Panc-1, BxPC-3, and MiaPaca-2. These experiments showed evidence for 8 genes conferring resistance with at least one guide in at least one cell line. Overall, 10% of CRISPRko and 22.7% of CRISPRa screens showed statistically significant validation as measured by increased EC50 (Table S3, Figure S6A). In the CRISPRa system, we observed that guides that resulted in a high degree of over-expression were more likely to significantly affect EC50 ( Figure S6B).
The top multi-drug resistance gene identified by our CRISPRa screen and validated using independent guides in multiple cell lines was the ATP-binding cassette (ABC) transporter, ABCG2 ( Figure 2A). This gene has been associated with multi-drug resistance primarily in breast cancer, and to some extent in pancreatic cancer cell lines (Bhagwandin et al., 2016;Mo and Zhang, 2012;Sun et al., 2016). ABCG2 functions as an efflux pump with broad range substrates and our screen data supported these observations with CRISPRa-mediated activation of ABCG2 promoting resistance to all four drugs in our screen ( Figure 2B). In our screen, ABCG2 activation led to particularly strong resistance to oxaliplatin treatment, with a six-fold enrichment observed for the top ABCG2 sgRNA in treated relative to control cells. This sgRNA was capable of increasing ABCG2 expression to approximately 25 times the baseline expression level exhibited in control cells expressing a scrambled sgRNA and recapitulated the effect observed in the pooled screen in an individual sgRNA validation experiment using a different two vector CRISPRa system ( Figure 2C,D). RNA-sequencing data does not suggest significant off-target effects ( Figure S7). Furthermore, treatment with KO143 and sorafenib, two inhibitors of ABCG2, was able to reverse the resistance to irinotecan observed in the context of ABCG2 overexpression (2E, Allen et al., 2002;Wei et al., 2012).
Although ABCG2 activation proves to be a robust mechanism for multi-drug resistance, its expression level is not strongly correlated with patient survival suggesting it may not be a common mechanism of drug resistance in PDAC patients ( Figure S8). Nonetheless, ABCG2 is expressed at levels above background with expression increased in some patients for whom our findings may be relevant.
Novel pathways are associated with multi-drug resistance Examination of top genes in our CRISPRa and CRISPRko screens revealed enrichment of genes involved in a set of Reactome pathways. Many of these pathways also contain genes whose expression was associated with PDAC patient survival ( Figure 3A-B, Table S4). These pathways, which when activated induce drug resistance and when knocked out confer drug sensitivity, represent potentially useful therapeutic targets. One pathway meeting these criteria is the hemidesmosome assembly pathway ( Figures 4A, S8). The hemidesmisome is a multi-protein complex that adheres epithelial cells to the basement membrane (Walko et al., 2015) and has been strongly linked to cell adhesion-mediated drug resistance (CAM-DR, Correia and Bissell, 2012). One of the strongest hits in this pathway, PLEC, acts as a linker protein for cytoskeletal scaffolding and has been previously associated with BRCA2 activity in breast cancer ( Figure S9, Niwa et al., 2009). Several genes in the hemidesmosome pathway show enrichment in our screen ( Figure 3C). Our validation results showed PLEC activation was most associated with oxaliplatin resistance in MiaPaca-2 and BxPC-3 cell lines ( Figure 3D). In patients, genes of this pathway were highly expressed in the squamous subtype of PDAC, which was described as conferring a particularly poor prognosis ( Figure 3E, Bailey et al., 2016). The mean expression level of genes in this pathway was also significantly greater in treatment naïve tumors from PDAC patients who succumbed to disease within 300 days (bottom quartile) relative to those who survived greater than 900 days (top quartile) in two independent cohorts ( Figure 3F, Kirby et al., 2016;Raphael et al., 2017). These data nominate the hemidesmosome as potential driver of poor prognosis of the squamous subtype.
Activation of several genes involved in chromatin remodeling also fostered drug resistance ( Figure 3A,B). The role of chromatin remodeling in pancreatic cancer has been heavily investigated as several genes involved in histone methylation (MLL2 and MLL3) and members of the tumor suppressing SWI/SNF complex (SMARCA1 and ARID1A) are recurrently mutated in PDAC tumors (Ryan et al., 2014). Moreover, global changes in repressive histone marks has been associated with metastatic PDAC tumors (McDonald et al., 2017). We found that activation of several genes involved in the repression of chromatin via histone deacetylation resulted multi-drug resistance ( Figure S10A). Particularly prominent were members of the Nucleosome Remodeling Deacetylase (NuRD) complex (MTA2/3, CHD3/4, HDAC1/2, and GATAD2A), members of the Nuclear receptor CoRepressor (NCoR) complex (NCOR2, TBL1XR1, and TBL1X), and the repressive transcription factor REST and its interacting partner RCOR1, each of which have been demonstrated to regulate gene expression programs important for chemoresistance (Hou et al., 2017;Wong et al., 2014). Histone deacetylase inhibitors have demonstrated promise in PDAC pre-clinical studies and are currently being investigated in multiple clinical trials (Feng et al., 2014). In agreement with these findings, over-expression of HDAC2 demonstrated resistance to oxaliplatin relative to control cells ( Figure S10B), however our validation was likely limited by the level of HDAC2 overexpression achieved for this already highly expressed gene ( Figure S6B).

CRISPR screen data can be used globally to predict chemosensitivity
Our screen yielded a large number of candidate genes and the validation of top hits suggests that our results are rich with new genes relevant to drug resistance. Since it is not feasible to individually validate all interesting genes, we sought to globally validate data from our screen. We reasoned that if our screen represents a gene's ability to confer drug resistance, then expression of genes that were top hits in our screens could be used to predict drug sensitivity in cell lines. To test this, we summed the product of a cell line's scaled expression data and the square of the L2FC sum for a given drug for each gene to generate a putative resistance score ( Figure 4A). This method was tested on two independent previously published data sets characterizing PDAC cell line drug sensitivity. We found our predicted drug response scores significantly associated with observed irinotecan and gemcitabine sensitivity in the 18 pancreatic cancer lines assayed by cancer cell line encyclopedia (Barretina et al., 2012) and 14 lines from our own published work ( Figure 4B-C, Kirby et al., 2017). Our drug response predictions were also able to respond to dynamic changes within a cell type. An increase in predicted drug resistance was observed in each of seven PDAC cell lines after a weeklong treatment with a sub lethal dose of gemcitabine ( Figure 4D). Finally, we assessed the performance our gemcitabine resistance predictions in patients who received adjuvant gemcitabine chemotherapy, but found our predictions were only meaningful for patients at the extreme ends of predicted sensitivity or resistance ( Figure 4E).

DISCUSSION
The ability to both inhibit and endogenously activate the expression of genes in a modular and programmable manner has revolutionized forward genetic screening and is empowering discoveries that were previously impossible (Gersbach and Barrangou, 2018). We present the largest screen for cellular mechanisms of chemoresistance in PDAC cells to date and the first screen to perform both genome-wide knock out and endogenous activation in PDAC cells. The scale of our screen has allowed us to make a few fundamental observations of CRISPR screening that may be generalizable outside of PDAC. First, we found that the selection strength was the methodological variable most significantly correlated with our screen reproducibility.
The implication of this finding is that the false positive rate of a screen is likely to be dose dependent and investigators should make every effort to perform as strong of a selection as possible. This can be difficult to achieve in screens incorporating cytotoxic chemotherapies such as ours because cells often accumulate damage throughout an extended selection period and will exhibit residual death after the drug is removed from the media. Thus, significant effort should be made to ensure appropriate selection strength, which in our hands was roughly 80% selection relative to control cells, while maintaining sufficient survivorship for sequencing library construction post-selection.
The second general observation was a relatively weak correlation between replicate reproducibility and gene significance levels across cell lines. While this observation may be partially explained by methodological covariates like selection strength ( Figure S4), it is also likely that strong genetic and phenotypic differences observed between these cell lines (Kim et al., 2014) permit different mechanisms of chemoresistance. A gene capable of modulating chemosensitivity across every PDAC cell line would be an ideal target, however, our data suggests that this is not a realistic expectation. Previous screens for chemoresistance in PDAC cell lines have been limited in that they only use one cell line (Bhattacharjee et al., 2014;You et al., 2011) and our study is likely undersampling the natural variation in PDAC with only two cell lines. As sequencing and screening costs decrease, future studies examining orders of magnitude greater numbers of cell lines will likely add to the knowledge catalogued by our study.
Our last general observation is that genes and pathways associated drug-specific resistance tended to be less associated with PDAC patient survival than multi-drug resistance pathways. This result mirrors the reality of current PDAC clinical management in which patients are routinely treated with multiple drugs as frontline therapy (Garrido-Laguna and Hidalgo, 2015) and, at least at the cellular level, drug-specific mechanisms of resistance likely account for a smaller proportion of chemoresistance. As the rational design and delivery of multi-drug cocktails becomes more sophisticated and widely used we expect this observation to be made in other cancers as well (Hu et al., 2016).
We found ABCG2 to be the most consistent gene associated with resistance across each drug and both cell lines. ABCG2 is a well-described gene capable of mediating multi-drug resistance and as a efflux pump capable of acting upon a broad range of substrates, the mechanism of resistance is clear and helps to confirm the validity of our findings. Our validation data confirm that two independent guides confer resistance to multiple drugs in multiple cell lines. However, we were unable to find evidence of ABCG2 expression being correlated with patient survival in PDAC patients, thus the clinical relevance is unknown, however some evidence suggests that ABCG2 expression is localized to cancer stem cells, which represent a minority of the total cell population, thus ABCG2's role may be obscured from bulk tumor sequencing (Bhagwandin et al., 2016). It is for that reason we sought to demonstrate that inhibition of ABCG2 has the potential to sensitize cells to cytotoxic chemotherapy. Our findings suggest that further investigation of these agents would be warranted.
We also identified several components of the hemidesmisome complex as capable of mediating multi-drug resistance and possessing a strong correlation with patient survival across multiple cohorts. The hemidesmisome is critical for cells to adhere to the basement membrane and has been described as a core regulator of cell adhesion-mediated drug resistance (CAM-DR) (Correia and Bissell, 2012). CAM-DR is increasingly recognized for its importance in PDAC and the evolutionary trade offs between cells committing to epithelial-to-mesenchymal transition (EMT) in favor of metastasis as opposed to maintaining connections to the basement membrane and engaging in CAM-DR is an active area of investigation(Makohon-Moore and Iacobuzio-Donahue, 2016). As expression of hemidesmisome genes are strongly associated with quasi-mesenchymal/squamous subtype, genes identified by this study represent candidate markers and drivers of the poor patient prognosis described with this subtype across multiple studies (Bailey et al., 2016;Collisson et al., 2011;Moffitt et al., 2015;Raphael et al., 2017).
While the hemidesmisome is not readily targetable, the quasi-mesenchymal/squamous subtype has been previously characterized as sensitive to the EGFR inhibitor, erlotinib. Erlotinib was approved for PDAC despite marginal survival benefits and its widespread use has been limited by lack of a reliable biomarker capable of identifying responsive patients (Garrido-Laguna and Hidalgo, 2015). Performing tumor RNA-sequencing to cluster tumors into corresponding subtypes, which often have fluid definitions across studies, may be difficult clinically, however measuring a few components of the hemidesmisome complex identified by our study, such as PLEC, may be clinically tractable.
We also identified several transcriptional repressors including REST, the NCoR complex, and the NuRD complex as important regulators of chemosensitivity. Previous complementary studies have found recurrent inactivating mutations in several transcriptional activators including the SWI/SNF complex and the histone methyltransferases MLL2 and MLL3 in PDAC patients (Ryan et al., 2014). The transcriptional repressors identified by our study, such as HDAC2, are responsible for broad epigenetic reprogramming, thus future investigations will involve narrowing down the key epigenetic changes linked to resistance.
Lastly, we used our CRISPRa data to generate an algorithm for predicting sensitivity to drugs assayed in our study. We provided evidence of our algorithm's performance for gemcitabine and irinotecan in two independent studies that cataloged drug sensitivity and transcription of several PDAC cell lines. Our algorithm performed less well in predicting the survival PDAC patients treated with gemcitabine, only holding predictive power at the extreme ends of sensitivity prediction. This is likely explained by multiple non-cellular factors contributing to patient survival, such as comorbidities, patient drug metabolism, immune response, and stage at diagnosis. Thus cellular mechanisms of resistance may only dictate patient response in outlier cases. Furthermore, the expression data used to predict patient response was obtained from treatment naïve tumors and tumor sensitivity profiles may change at disease recurrence. We hope this effort inspires further attempts to integrate screening data with transcriptional profiling to predict drug response and facilitate rational combination chemotherapy.

ACKNOWLEDGEMENTS
We thank the Myers and Cooper labs for helpful feedback. We also thank the HudsonAlpha genome services laboratory for assistance and optimizing the sequencing parameters of our CRIPSR screen. We also acknowledge The Cancer Genome Atlas, Cancer Cell Line Encyclopedia project datasets, which were extremely valuable and without which this study would not be possible. This work was supported by the State Cancer Fund of Alabama (to RMM). RCR and AAH were funded by the UAB MSTP (NIH-NIGMS 5T32GM008361-21. SJC and ERG were supported by the HudsonAlpha Tie the Ribbons Fund. SJC received support from the UAB CCTS grant (NIH 1UL1TR001417-01) and the UAB Comprehensive Cancer Center (NIH 5P30CA013148).

CONTRIBUTIONS
RCR, AAH, SJC and SJC designed the experiments. RCR, AAH, and ERG collected data. RCR, AAH, and ERG analyzed the data. RCR, AAH, and SJC wrote the first draft. All authors contributed to writing of the paper and read and approved the final manuscript.

DECLARATION OF INTERESTS
The authors report no conflicts.

Contact for Reagent and Resource Sharing
Further information and requests for resources and reagents should be directed to and will be fulfilled by the Lead Contact, Sara Cooper (sjcooper@hudsonalpha.org).

Library amplification
The GeCKO A and SAM library were amplified as described previously (Joung et al., 2017). Briefly, for each pooled library 8 electroporations were performed using 1uL of library and 25uL of Lucigen Endura Electrocompetent cells (Lucigen #60242). Electroporations were recovered in 1.975mLs of recovery media and rotated at 250 rpm for 1hr at 37 o C. All electroporations were pooled after incubation and 2mLs of transformation were plated on 8-25cm 2 bioassay plates with LB agar containing ampicillin. Plates were grown for 14hrs at 35 o C.
Transformation efficiency was greater than 1x10 8 for both libraries. Bacterial colonies were scraped off of bioassay plates in 20mL of LB broth, pelleted and stored at -20 o C for less than 1 week. Plasmid DNA was extracted from greater than 3g of pellet using the Qiagen EndoFree Plasmid Maxi Kit (Qiagen #12362).

Viral packaging
HEK 293FT cells (ThermoFisher #70007) were cultured in DMEM (ThermoFisher #11965) supplemented with 10% FBS (VWR #16777), 0.5% Penicillin-Streptomycin (ThermoFisher #15140122), 0.5% Non-essential amino acids (ThermoFisher #11140), and 0.5% sodium pyruvate (ThermoFisher #11360070). Viral packaging of pooled sgRNA libraries was performed in three 225cm 2 flasks. One hour prior to transfections, after cells had reached ~60% confluency, DMEM was removed and replaced with 13mL of pre-warmed OptiMEM (ThermoFisher #31985). For each 225cm 2 flask, 15.3ug of pMD2.G, 23.4ug psPAX2, and 30.6ug of pooled library were combined with 200uL of Lipofectamine PLUS reagent (ThermoFisher #15338) in 4mL of OptiMEM and vortexed to mix. Separately, 100uL of Lipofectamine LTX (ThermoFisher #15338) was diluted in 4mL of OptiMEM, vortexed, and incubated at room temperature for 5 minutes. The lipofectamine and plasmid solutions were then combined, vortexed, and incubated for 5 additional minutes at room temperature. The 8mLs of combined solution were added to each flask. After incubating for 6 hours, PptiMEM was replaced with the original DMEM culture media supplemented with 1% bovine serum albumin (ThermoFisher #A9647). After 60 hours, the media was harvested, centrifuged at 3000rpm for 10 minutes at Viral packaging of single sgRNAs used for validation was performed in 6-well tissue culture treated plates. The day prior to transfection, 750,000 HEK293T cells were seeded into a 6 well plate; minimum one well per sgRNA. One hour prior to transfections, after cells had reached ~90% confluency, DMEM was removed and replaced with 3mL of pre-warmed OptiMEM (ThermoFisher #31985). For each well to be transfected, 1ug of pMD2.G and 1.5ug psPAX2 were combined with 12.5uL of Lipofectamine 2000 reagent (ThermoFisher #11668019) in 200uL of OptiMEM and vortexed to mix. Separately, 10uL of Lipofectamine PLUS reagent (ThermoFisher #15338) was diluted in 200 uL of OptiMEM and vortexed to mix. Both master mixes are combined and added to a tube containing 2.5ug of sgRNA plasmid, vortexed, and incubated at room temperature for 5 minutes. The 400uL of combined solution were added to each corresponding well in the 6-well late. After incubating for 48 hours, the media was harvested, centrifuged at 3000rpm for 10 minutes at 4 o C, filtered through a 0.45um filter (VWR #28145). Virus was either immediately stored at -80 o C or used for transfection.

Drug resistance screening
Panc-1 and BxPC3 cell lines were transduced with LentiCas9-Blast or co-transduced with Lenti-dCAS9-VP46-Blast and lenti-MS2-p65-HSF1-Hygro for knock out and activation screening respectively. 500,000 cells were seeded into each well of a six-well flask with 0.8ug/mL polybrene (Sigma #TR-1003-G). Previously packaged plasmid added at volumes that ranged from 0 to 500uL of virus. Cells were then centrifuged at 2000rpm for 2 hours at 37 o C.
After centrifugation the cells were incubated overnight and media was changed the following morning. The following day cells were treated with appropriate antibiotic. Panc-1 cells were treated 10ug/mL Blasticidin and 1500ug/mL Hygromycin and BxPC3 cells were treated with 2.5ug/mL Blasticidin and 1250ug/mL Hygromycin. Cells were kept under antibiotic selection for one week post-transduction until no-virus control cells were dead. The lowest viral titer that had a sufficient number of cells to carry forward were grown up such that pooled libraries could be transduced at 500x representation at a multiplicity of infection (MOI) of 0.4: Min. cells needed = 70,000 sgRNAs x 500 cells/sgRNA x (1/0.4) = 87.5x10 6 cells Pooled library transduction was performed in at least 30 wells of a 12 well flask at 3x10 6 cells per well with 0.8ug/mL polybrene and sufficient volume of concentrated virus to reach an MOI of 0.4 (approximately 30% of cells surviving after antibiotic selection). Cells were transduced as described above. The following day cells were split out into larger 225cm 2 flasks and either puromycin (Panc-1 10ug/mL, BxPC3 2.5ug/uL) or zeocin (Panc-1 2mg/mL, BxPC3 3mg/mL) was added to select for presence of knock out and activation plasmids, respectively.
Library-transduced cells were under selection for one week post-transduction and expanded to 7x10 7 cells per treatment replicate or 1000x representation for each of the 4 drug and control conditions per replicate. A minimum 500x representation was maintained at all times in control cells. Drug treatments included gemcitabine (Tocris #G6423, Panc-1 25nM, BxPC3 25nM), oxaliplatin (Sigma #O9512, Panc-1 2.5uM, BxPC3 2.5uM), irinotecan (Sigma #I1406, Panc-1 500nM, BxPC3 250nM), and 5-fluorouracil (Sigma #F6627, Panc1 7.5uM, BxPC3 5uM). All doses were optimized to yield ~20% the number of viable cells in drug treated cells relative to untreated control cells after 14 days of culture. After 14 days of drug treatment, cells were counted, pelleted, and stored at -80 o C. The first replicate of the knock out screen for each cell line was performed separately from the second two in order to optimize all drug doses and library preparation methods. All three replicates of the activation screen were performed in parallel for each cell line.

DNA extraction and library preparation
Genomic DNA was extracted from cell pellets adapting a previously published protocol (Chen et al., 2015). The protocol described below assumes 5x10 7 input cells, but can be scaled proportionally to accommodate different cell ranges. Briefly, cells were resuspended in 6mLs of a NK lysis buffer (50mM Tris (ThermoFisher #BP1757)), 50mM EDTA (Fluka #03690), 1% SDS (Invitrogen #15525), and 30uL of 20mg/mL Proteinase K (Zymogen #D3001)) and incubated at 55 o C overnight. The following morning, 30uL of 10mg/mL RNAse A (Qiagen 19101) was added to the lysate, inverted several times, and incubated at 37 o C for 30 minutes.
The supernatant was carefully decanted to a fresh conical tube and 6mLs of 100% isoproponal (Sigma #I9516) was added to the tube inverted 50 times and centrifuged at 4100xg for 10 minutes. The supernatant was then discarded, 6mLs of 70% ethanol (Sigma #E7023) was added to the DNA pellet, and centrifuged at 4100xg for 1 minute after inverted 10 times. The supernatant was pipetted off carefully and the DNA pellet was dried at room temperature for 1 hour. 500uL of TE buffer (Sigma #T9285) was added to the dried pellet and the tube was incubated at 65 o C 1-2 hours. Next the conical tubes were briefly spun at 1000xg to pull down any evaporated TE buffer and incubated for 2-3 days at room temperature vortexing periodically.
If necessary, additional TE buffer was added until the DNA pellet was completely dissolved. The gDNA concentration was measured using Qubit fluorometric quantitation (ThermoFisher #Q32850). Two-step PCR reactions were used to amplify and append sequencing tails to both the activation and knock out libraries. Primer sequences are available in Table S5. To maintain 500x representation in control samples we performed 86 50uL PCR1 reactions (or as many reactions as possible for drug selected samples) per replicate each containing 2.5ug gDNA, 2.5 uL of appropriate forward and reverse primer (25uM), 25uL NEBNext High-Fidelity 2X PCR master mix (NEB #M0541), and the remaining volume required to reach 50uL in nuclease-free water.
For the knock out library, we performed an initial denaturation of 98 o C for 60 seconds, 18 cycles of 98 o C denaturation for 15 seconds, 62 o C annealing for 30 seconds, 72 o C extension for 30 seconds, a final extension of 72 o C for 5 minutes, and a 4 o C hold. The activation screen PCR1 was carried out under identical conditions except the annealing temperature was adjusted to 60 o C and only 16 cycles of PCR were performed. PCR1 reactions for each replicate of each treatment condition were pooled and the second PCR reaction (PCR2) was used to append Illumina sequencing tails.
We performed 12 PCR2 reactions for each replicate of each condition consisting of 2.5uL of pooled PCR1, 2.5uL of appropriate forward and reverse primer (10uM), 25uL NEBNext High-Fidelity 2X PCR master mix, and 20uL of nuclease-free water. The PCR conditions were carried out similarly to PCR1 except it was reduced to 10 cycles for each reaction with an annealing temperature of 65 o C for the knock library amplification and 69 o C for the activation library amplification. All PCR2 reactions for each replicate of each condition were pooled and 50uL was run on a 2% gel at 125V for 2 hours. A 240bp or 280bp band for the knock out or activation library respectively was excised and gel extracted (Qiagen #28704). All controls and drug treated samples for each replicate were prepped in the same batch in order to avoid batch affects.

Single gene validation
The top sgRNAs for genes of interest were cloned into either LentiCrispr-v2 or LentiSAMv2 for knock out or activation screen validation as described previously (Konermann et al., 2015). (Sequences are available in Table S6). Plasmids putatively containing sgRNAs of interest were transformed into One Shot Stbl3 competent E. coli (ThermoFisher #C737303).
Colonies that survived ampicillin selection were cultured overnight in 10mL of LB broth containing ampicillin. Plasmids were harvested with an E.Z.N.A Endo Free plasmid mini kit (Omega #D6950) and confirmed to contain the appropriate sgRNA with Sanger sequencing. To generate stable cell lines expressing each sgRNA, single sgRNA plasmids underwent viral packaging and were transduced as described previously, followed by transduction. For transduction 10 6 cells were seeded per well of a 6-well plate for a total volume of 2mL of cells and media. To that we added, 1.6ul of polybrene (Millipore Sigma #TR-1003-G), and 1mL of packaged virus. Then, we used spinfection for 30 minutes at 2000RPM at 37C. All lines were established using antibiotic selection as described above.
Non-targeting and targeting sgRNAs were plated in 96-well plates at 750 cells/well and treated with a range of gemcitabine, oxaliplatin, or irinotecan concentrations. Plated cells were grown for 6 and 8 days, for the knockout and activation screen, respectively, with the media changed and drug applied every 2 days. The number of viable cells surviving drug treatment was assayed with Cell Titer-Glo cell (Promega #G7571). ABCG2 inhibition was done under these conditions as well, with the only change being the addition of 3uM Sorafenib or 3uM KO143 where indicated. Two NTCs were tested and the complete data sets is available in Figure S11.

RNA-sequencing
RNA was extracted from cell pellets containing 3 x 10 5 cells. Cells were suspended in 400 ul of RL Buffer from the Norgen total RNA extraction kit (cat. # 37500, 25720) plus 1% BME.
Samples were homogenized by vortex. We then continued extraction of RNA with the Norgen Total RNA extraction kit including the Norgen DNase kit (cat. # 37500, 25720). RNA integrity numbers (RIN) were measured using the BioAnalyzer (Agilent) and all measured at a RIN of 10 indicating high quality RNA. RNA was quantified using Qubit 2.0 Fluorometer (Thermofisher).
We used 1000 ng of total RNA as input to the NEBNext Poly(A) mRNA Magnetic Isolation Module (NEB cat# E7490S) followed by the NEBNext Ultra RNA Library Prep Kit for Illumina (cat# E7530S). Libraries were barcoded using the Genomic Services Lab at HudsonAlpha's custom dual index barcode plate. We pooled all samples onto a lane of a S4 NovaSeq flow cell, and sequenced an average of 95.8 million 150bp paired-end reads per sample with an average Q30 score of 95%. Samples were processed using our published primary analysis tool, aRNApipe (Alonso et al., 2017).

Quantification and Statistical Analysis
Sequencing and data processing Three sets of replicates (control and 4 drug treated samples) for each cell line were sequenced on one lane of Illumina NextSeq resulting in an average of 40 million reads per sample. 5' and 3' adapters were observed in >99% of reads and were trimmed using cutadapt (Martin, 2011). Adapter trimmed fastq files were then aligned to the gRNA libraries and raw count tables generated using MAGeCK . We had a perfect alignment rate of 72%-74% of raw reads for each sample. Sequencing the control samples revealed sufficient representation of guides with an average of 99.8% and 99.3% of sgRNAs detected and 98.4% and 94.2% of sgRNAs detected at greater than 1 read per million in our knock out and activation control samples, respectively. A log 2 fold change was computed for each sgRNA in each drug treated sample relative to the untreated control sample for each replicate. sgRNAs with less than 10 counts in the untreated control samples and less than 50 counts in a treated control sample were excluded from further analysis (<1% of sgRNAs)(Table S1A-B). At this step we recognized the first replicate of our knock out screen correlated relatively poorly with the second two replicates and was excluded from downstream analysis. We ranked each of the sgRNAs targeting each gene by the minimum log 2 fold change across each replicate. Top genes were subsequently prioritized for follow up by their "L2FC sum" in each cell line, which is the sum of the replicate minimum log 2 fold changes of the top two sgRNAs targeting each gene. Multi-drug hits were prioritized by computing the mean "L2FC sum" across each of the four drug treatments. Pathways enriched for genes conferring drug resistance were identified by comparing the distribution of log 2 fold changes for the top, second, and third sgRNA targeting each gene within a Reactome pathway to that of all other genes by Wilcox ranked sum test. Reactome (Fabregat et al., 2018) pathways with fewer than 10 genes targeted in our libraries were excluded. A consolidated knock out and activation score was computed for each pathway by summing the -log 10 P-values for the top, second, and third sgRNAs of each pathway (Table S4).

Curve fitting single gene validation
Where EC50 values were calculated, we use cell count data as measured by Cell Titer Glo, as input, and fit those data to a rectangular hyperbola which led to the use of we used Michaelis-Menten equation to curve fit the data and derive the estimated drug concentration at which the drug is half its maximum effectiveness and confidence intervals for each cell line, drug, and guide combination. We used those data to determine the change in EC50 and p-value associated with those changes for each guide compared to each of the two non-targeting controls (NTCs) . These calculations were done using Prism GraphPad 7.

Predicting drug response using genome-wide screening results
Sensitivity to each of the four drugs was computed using a cell line's treatment naïve gene expression levels and the minimum L2FC Sum for each gene across both cell lines. Raw expression data for each cell line, obtained from sources described below, was processed into raw count tables using the aRNApipe RNA-seq processing pipeline (Alonso et al., 2017). Cell line gene expression data was normalized using the R package DESeq2's variance stabilizing transformation  and each gene's expression level was z-scored across each cell line. Next, a scalar by which to weight each gene's expression level was computed using the L2FC Sum described above. The replicate minimum L2FC sum from the relevant drug's CRISPRa screen was used for all genes who's L2FC Sum was greater than 0. The replicate maximum L2FC sum was used for all genes whose L2FC Sum was less than 0. All other genes were assigned a scalar value of zero. The L2FC Sum scalar was squared to ensure a positive value and place more weight on genes with high absolute fold change. Next, the product of each gene's z-scored expression data and L2FC Sum scalar was computed to generate a weighted expression value for each gene. Finally, the cumulative sum of the weighted expression value for each gene with a negative L2FC Sum in each cell line was subtracted from the cumulative sum of the weighted expression value for each gene with a positive L2FC Sum to generate a single value representing a cell line's expected level of resistance to the given drug. Irinotecan sensitivity and gene expression data were obtained for 18 PDAC cell lines with permission from