CRISPR Screens Identify Essential Cell Growth Mediators in BRAF-inhibitor Resistant Melanoma

BRAF is a serine-threonine kinase that harbors activating mutations in ∼7% of human malignancies and ∼60% of melanomas. Despite initial clinical responses to BRAF inhibitors (BRAFi), patients frequently develop drug resistance. To identify candidate therapeutic targets for BRAFi-resistant melanoma, we conducted CRISPR screens in melanoma cells harboring an activating BRAF mutation that had also acquired resistance to BRAFi. The screens identified pathways and genes critical for BRAFi resistance in melanoma cells. To investigate the mechanisms and pathways enabling resistance to BRAFi in melanomas, we integrated expression data, ATAC-seq, and CRISPR screen results. We identified the JUN family of transcription factors and the ETS family transcription factor ETV5 as key regulators of CDK6 that enabled resistance to BRAFi in melanoma cells. Our findings reveal genes whose loss of function conferred resistance to a selective BRAF inhibitor, providing new insight into signaling pathways that contribute to acquired resistance in melanoma.


41
Melanoma is an aggressive malignancy with a poor prognosis. The median survival for 42 patients with stage IV melanoma ranges from 8 to 18 months after diagnosis, depending on 43 the substage [1]. Somatic mutations in BRAF, most commonly V600E or V600K [2], are the 44 most frequently identified cancer-causing mutations in melanoma, and recurrently appear in 45 colorectal cancer, non-small cell lung carcinoma, and many other cancers [3]. BRAF encodes 46 a protein belonging to the RAF family of serine/threonine protein kinases. This protein plays 47 conditions strongly enriched for roles in fundamental biological processes, such as gene 134 expression, RNA processing, and translation ( Figure S3C and D). These results are consistent 135 with a properly functioning CRISPR screen. 136

Identification of genes essential specifically for growth of cells resistant to PLX-4720 137
To explore which genes might play a role in the BRAFi-resistance, we performed further 138 analysis of CRISPR screen data using MAGeCKFlute [29]. MAGeCKFlute facilities 139 comparison of β score between different conditions. We adopted a "quantile matching" 140 approach to robustly estimate σ , which is the standard deviation of the differential β score 141 ( Figure S4A). We identified genes whose β score decreased in the presence of BRAFi 142 treatment compared to DMSO treatment ( Figure S4B and Table S2). Then, we selected 322 143 candidates whose disruption does not normally affect survival but becomes lethal in the drug 144 treatment condition. We ranked the identified hits by the change of the β score ( Figure 1D). 145 Here ,we labeled the top 10 genes, such as SOS1, PURA, HRAS, SAFB, CRKL, ETV5, CDK6, 146 DYNCH1, H2AFX and MAZ. Among these 322 candidate genes, HRAS, SRC, SOS1, EGFR,147 and RAF1 were previously reported to be involved in BRAFi resistance [31,32] (Figure 1E). 148 To further understand the pathways conferring BRAFi-resistance, we performed 149 GO/GSEA/pathway analyses with the 322 candidate genes ( Figure 1F). Among the network 150 of genes whose β score decreased after drug treatment, we found that the ERBB2 signaling 151 pathway, RAS pathway, ERK pathway, MAPK pathway, and EGFR signaling pathway are 152 highly enriched. These results are consistent with previous studies [13,31,33,34]. Besides 153 these known pathways, cell-cycle genes, and G1/ transition S of mitotic cell cycle were the 154 most enriched newly discovered class ( Figure 1F), represented by CDK6, CCND1, PSMB1, 155 and RRM2. 156

CDK6 confer resistance to BRAF inhibition in melanoma cells 157
We next sought to determine whether any genes whose upregulation confers resistance to 158 BRAF inhibition in melanoma cells. To assess this, we analyzed previously generated gene   Table S3). Our re-analyses confirmed the previously reported overexpression 168 of KIT, MET, EGFR, and PDGFRB in M238R1 relative to the parental line [13]. In addition, 169 we found that the cell cycle genes CDK6, CCND1, and transcription factor (TF) JUN were 170 up-regulated in resistant cells compare to the parental cells ( Figure 2A). 171 We hypothesized that genes with elevated expression in BRAFi resistant cells, as well as 172 the loss of function restored the drug sensitivity, may be responsible for the resistance 173 phenotype. We next integrated the expression results and CRISPR screen results to identify the dysregulated genes related with BRAFi resistance. Within the 322 genes whose depletion 175 sensitize cells to BRAFi, there are 12 genes, including CDK6, specifically over-expressed in 176 BRAFi-resistant cells ( Figure 2B). This suggests that 21 genes might be associated with the 177 resistance to BRAFi and mediate cell proliferation in the resistance line. 178 To explore the potential druggable targets for the BRAFi-resistant cells, we further 179  Figure S6A), suggesting that loss-of-function of CDK6 can cause cells 187 sensitive to PLX-4720. To validate this result from the initial screen, we used five 188 independent sgRNAs to knockout CDK6 in the M238R1 cell line ( Figure 2C). Consistent 189 with our screen data, CDK6 knockout cells showed increased sensitivity to PLX-4720 in 190 long-term colony-formation viability assays ( Figure 2D). Most tumors including melanoma 191 have an abnormal G1-to-S transition, mainly due to dysregulation of CDKs activities [38,39]. 192 We wondered if the increased essentiality we observed for CDK6 was a general property of 193 CDKs or was specific to CDK6. We specifically evaluated the changes in essentiality of the 194 other CDKs ( Figure S6B) The ATAC-seq profiles showed the high-quality features according the criteria defined  209 by Cistrome database, which is a data portal for more than 8,000 ChIP-Seq and chromatin 210 accessibility data in human and mouse [41]. 211 In total, 113,725 high-confidence open chromatin regions (or peaks) were identified in the 212 parental line, and 96,038 peaks were identified in resistant line. Of the distinct peaks, we 213 identified the peaks more accessible in parental cells (M238-specific peaks), and the peaks 214 more accessible in resistant cells (M238R1-specific peaks) ( Figure 3A and Table S5). To assess other genes that might act with JUN to regulate CDK6, we examined the set of 241 genes that physically interact with the JUN protein according to the STRING database and 242 genes whose essentiality increased after BRAFi treatment. We identified ETV5 as being in 243 both of these gene sets ( Figure 4C Figure 4D). Consistent with the hypothesis that ETV5, JUN, and JUNB directly regulate 250 CDK6, these TFs have strong binding around the CDK6 gene ( Figure S8C). We found that 251 ETV5 deletion reduced sensitivity to BRAF inhibition by PLX-4720 in melanoma cells and 252 ETV5 was the top hit of the genes that were more essential in the BRAFi treatment condition 253 ( Figure 1D). Similar to CDK6, the normalized sgRNA read counts of ETV5 continually 254 decrease in the DMSO treatment or PLX-4720 treatment ( Figure S8A Figure S9). These results support the potential of CDK6 and BRAF 275 dual inhibition as a therapeutic strategy to overcome BRAFi resistance in our resistant model.

CDK6 expression is negatively associated with overall survival in BRAF-mutant 277 melanomas treated with BRAFi 278
To determine whether the expression of any validated BRAFi-resistant genes we identified 279 correlated with resistance to BRAF inhibitor therapy in melanomas, we analyzed expression 280 data from two independent cohorts [33,51]. In cohort one [51], 18 patients were treated 281 either with BRAFi alone (12 patients) or dual BRAFi and MEKi therapies (6 patients). RNA-282 seq data on serial tumor biopsies of matched pre-treatment and relapsed tumors were 283 available. In cohort two, 22 patients with advanced melanoma were treated with BRAFi (7 284 patients) or BRAFi plus MEKi (15 patients) [33]. RNA-seq data on pre-treatment, on-285 treatment, or relapsed tumors were available, although they were not paired. These samples 286 were classified into 3 groups: 14 pre-treatment specimens, 12 on-treatment specimens, and 12 287 clinical progression specimens. Of the 21 over-expressed genes also identified in our 288 CRISPR screen, CDK6, CCND1, and ETV5 were more highly expressed in the tumors that 289 have relapsed after BRAFi treatment relative to the on-treatment groups ( Figure 6A). 290 We next investigated whether CDK6 upregulation might be associated with clinical 291 resistance in some cases. To facilitate this, we generated a 10-gene CDK6 expression 292 "signature" (Table S6) BRAFi-treated melanoma patients. Overall, these observations provide initial support for the 308 notion that CDK6 upregulation by transcription factors JUN and ETV5 might be associated 309 with clinical resistance to BRAFi in melanoma patients.

Library design 397
To design a smaller-scale CRISPR/Cas9 knockout screen library focusing on cancer-related 398 genes, we selected 6000 genes based on reported relevancies with cancers using multiple 399 sources, including Cosmic and Oncopanel (Table S1). For each gene, we designed ten 19nt 400 single-guide RNA (sgRNA) against its coding region with optimized cutting efficiency and 401 minimized off-target potentials. We used sequence features of the spacers to calculate the 402 cutting efficiency score for each sgRNA using our predictive model. We used BOWTIE to 403 map all candidate sgRNAs to hg38 reference genome, and chose those with fewest potential 404 off-targets. We selected the 10 best sgRNAs for each gene based on the considerations above. 405 The library also contains both positive controls and two types of negative controls: non-406 targeting controls and non-essential-region targeting sgRNAs. 407 a) Positive controls: we included 1466 sgRNAs targeting 147 positive control genes, which 408 are significantly negatively selected in multiple screen conditions. 409 b) Non-targeting negative controls: 795 sgRNAs with sequences not found in genome. c) Non-essential-region-targeting negative controls: 1891 sgRNAs targeting AAVS1, 411 ROSA26, and CCR5, which have been reported as safe-harbor regions where knock-in 412 leads to few detectable phenotypic and genotypic changes. 413

Cloning of individual sgRNAs and sgRNA libraries 414
For the 6K-cancer library, we used the lentiCRISPR v2 vector (also available at Addgene, 415 plasmid #52961) as backbone [58]. We designed ten sgRNAs per gene to target ~6,000 genes 416 and added non-targeting sgRNAs as controls (Table S1) Table S7. 421 After 8-12 h, the media was changed to 25 ml DMEM + 10% FBS+ 1%BSA. Viral 430 supernatant was collected two and three days after transfection, filtered through 0.45-μm 431 membranes, and added to target cells in the presence of polybrene (8 μ g/ml, Millipore). After 432 48h, puromycin (2 μ g/ml) was used to treat cells for two days for selection, which eliminated 433 all cells in an uninfected control group. 434

Pooled CRISPR screen 435
For the pooled CRISPR screen, a total of 1.2x10 8 cells were infected with the pooled 436 lentiviral library at a MOI of 0.3. After puromycin selection, the surviving cells were divided 437 into three groups (day0 control, vehicle, and drug treatment). For the drug treatment group, 438 the cells were treated with 1uM PLX4720. The cells were cultured in medium for ten 439 doubling times and split every 2-3 days before genomic DNA extraction and library 440 amplification. 441

Amplification and sequencing of sgRNAs from cells 442
After cell harvest, DNA was purified using QIAGEN DNeasy Blood & Tissue Kit according 443 to the manufacturer's instruction. PCR was performed as previously described [58], and the 444 PCR products were sequenced on a HiSeq 2500. Each library was sequenced at 30~40 445 million reads to achieve ~300X average coverage over the CRISPR library. The day 0 sample 446 library of each screen could serve as controls to identify positively or negatively selected 447 genes or pathways. 448

CRISPR screen analysis 449
The CRISPR/Cas9 screening data were analyzed using MAGeCK and MAGeCK-VISPR 450 algorithms [30]. MAGeCK-VISPR uses a metric, "β score", to measure gene selections. The 451 definition of the β score is similar to the term of 'log Fold Change' in differential expression 452 analysis, and β >0 (or <0) means the corresponding gene is positively (or negatively) selected, 453 respectively. We considered a β score of >0.5 or <-0.5 as significant. MAGeCK-VISPR 454 models the gRNA read counts as a negative binomial variable, whose mean value is 455 Tagmented DNA was purified using the MinElute Reaction Cleanup Kit (Qiagen, 28204). 484 The ATAC-seq library preparation was performed as described previously [40]. Then, the 485 concentration of the library was determined using Qubit 3.0 (Life Technologies) and the size 486 distribution was assessed using Agilent 4200 TapeStation system. Libraries were paired-end 487 sequenced (35bp) on an Illumina NextSeq 500. 488

ATAC-seq data analysis 489
Quality control, reads alignment, peak calling were performed by ChiLin  the difference between treatment and control beta scores. Red dots are genes whose beta 722 score increased after treatment. Blue dots are genes whose beta score decreased after 723 treatment. Gray dots are genes whose beta score did not change significantly between 724 different conditions. E. Beta scores of SOS1, RAF1, HRAS, EGFR, and SRC in the 725 PLX4720 condition and DMSO condition. F. Pathway enrichment analysis of the essential 726 322 genes whose β score decreased upon the BRAFi treatment compared to DMSO treatment.    A. Density plot of differential beta scores compared PLX4720 treatment condition with 794 DMSO treatment condition. Delta is used to measure the change of beta score in the two 795 conditions. Delta was calculated by the formula shown in the right panel. If the genes' 796 differential beta scores are bigger than delta (the red line), theses genes' essentiality 797 decreased after PLX4720 treatment. Delta is 0.134 in our screen data. Genes' differential beta 798 scores are smaller than minus delta (the blue line), which indicate theses genes' essentiality 799 increased after PLX4720 treatment. B. The beta score of M238R1 cell line with the treatment 800 of DMSO and PLX4720. The two diagonal lines indicate +/-1 delta of the difference between 801 treatment and control beta scores. Red dots are genes whose beta score increased after 802 treatment. Blue dots are genes whose beta score decreased after treatment. 803