The identification of genes maintaining cancer growth is critical to our understanding of tumorigenesis. We report the first in vivo genetic screen of patient-derived tumors, using metastatic melanomas and targeting 236 chromatin genes by expression of specific shRNA libraries. Our screens revealed unprecedented numerosity of genes indispensable for tumor growth (∼50% of tested genes) and unexpected functional heterogeneity among patients (<15% in common). Notably, these genes were not activated by somatic mutations in the same patients and are therefore distinguished from mutated cancer driver genes. We analyzed underlying molecular mechanisms of one of the identified genes, the Histone–lysine N-methyltransferase KMT2D, and showed that it promotes tumorigenesis by dysregulating a subset of transcriptional enhancers and target genes involved in cell migration. The assembly of enhancer genomic patterns by activated KMT2D was highly patient-specific, regardless of the identity of transcriptional targets, suggesting that KMT2D might be activated by distinct upstream signaling pathways.
Significance: Drug targeting of biologically relevant cancer-associated mutations is considered a critical strategy to control cancer growth. Our functional in vivo genetic screens of patient-derived tumors showed unprecedented numerosity and interpatient heterogeneity of genes that are essential for tumor growth, but not mutated, suggesting that multiple, patient-specific signaling pathways are activated in tumors. Cancer Discov; 6(6); 650–63. ©2016 AACR.
This article is highlighted in the In This Issue feature, p. 561
Overwhelming evidence suggests that tumor growth is sustained by gene mutations that confer a selective growth advantage to target cells (cancer-driver mutations; refs. 1–3). The advent of next-generation sequencing (NGS) has provided an initial view of the landscape of cancer-associated mutations (∼400,000 nonsynonymous mutations with ∼18,000 genes involved) and has uncovered great interpatient heterogeneity (10–200 mutations per tumor, with very few recurrent mutations; refs. 1–3).
Not all mutations, however, are drivers. The vast majority (>99.5%) have no effect on tumorigenesis and accumulate passively during tumor progression (passenger mutations). Putative driver mutations are currently identified by statistical methods (based on patterns and/or frequency of mutations), which allowed prioritization of a few hundred (1–3). Their mutational frequency is very low (∼60% mutated in <1% of patients), and, as a consequence, only <5% of patients can be treated with targeted therapies directed against mutated drivers (4).
The identification of mutated drivers requires formal demonstration of their tumorigenic role under conditions of in vivo cancer growth. Unfortunately, however, this information is available for only a fraction of the mutated cancer genes. Loss-of-function experiments targeting multiple genes were extensively used to investigate tumor vulnerabilities in vivo, using either cancer cell lines or genetically engineered tumors which, however, do not reflect the genetic diversity of human malignancies (5–14).
Here, we report the first genetic screen performed on human tumors grown in vivo. We used metastatic malignant melanomas, rapidly growing tumors with a 10-year survival rate of <10%. Around 50% of patients carry mutations in BRAF and can be treated with BRAF inhibitors. Strikingly, these drugs induce objective clinical responses in ∼50% of patients. However, virtually every patient eventually experiences disease progression during treatment, due to the emergence of resistant clones carrying secondary mutations, or activation of compensatory signaling pathways (15, 16). Thus, there is an urgent need to identify new cancer drivers and related signaling pathways, to exploit novel approaches for the treatment of melanoma.
Generation of Patient-Derived Xenografts of Metastatic Melanomas
For the genetic screens, we used patient-derived xenografts (PDX) of metastatic melanoma, obtained by injecting patient-derived bioptic samples (Supplementary Table S1) in NSG mice (NOD.Cg-Prkdcscid Il2rgtm1Wjl/SzJ; Fig. 1A). For each sample (PDX1), we generated secondary PDXs (PDX2). A total of 6 patients were included in this study, three carrying BRAF mutations at codon 600 and three carrying NRAS mutations at codon 61 (Supplementary Table S1). PDX1 and PDX2 samples were phenotypically indistinguishable from the patients' original tumors, as evaluated by histopathology (hematoxylin and eosin staining; not shown) and melanoma marker (S100, HMB-45, Melan-A, CD271, and MITF) analyses (Supplementary Fig. S1).
To ensure that our PDXs also retained genomic features of the original tumors, we compared whole-exome sequencing data of patients' tissues and corresponding PDX1s/PDX2s of three patients (two NRAS mutated: MM13 and MM16; one BRAF mutated: MM27). The vast majority of mutations [single nucleotide variants (SNV)/InDels] found in the patients' melanomas were also present in the corresponding PDXs (>98.8% in PDX1s and >97.3% in PDX2s; Supplementary Fig. S2A), including relevant melanoma driver mutations (e.g., NRAS, BRAF, RAC1, CDKN2A, NF1; not shown). Allele frequency of individual SNVs/InDels in the patients' samples was highly variable (5%–60%) and, most notably, was maintained in the PDX1 and PDX2 tumors (Supplementary Fig. S2B), suggesting that the growth pattern of the various cell subclones composing each tumor is retained after transplantation in NSGs. Thus, melanoma PDXs fully recapitulate both genomic and biological complexity of the patient tumors.
Unbiased In Vivo Pooled shRNA Screens in Patients with Metastatic Melanoma
One critical feature for the feasibility of in vivo genetic screens is the number of transduced tumor cells that grow after transplantation in NSG mice [tumor-initiating cells (TIC)], which must be sufficiently large to support the molecular complexity of the library (library representation). The challenge in achieving this condition is due to the fact that TICs might represent a fraction of the entire tumor population and that TIC frequencies (and their growth potential) can vary within the multiple subclones that characterize each tumor. TIC frequency is relatively high in melanoma samples (ref. 17; >1:86 in our tested samples, not shown). However, relative number and growth potential of TICs within each subclone cannot be directly quantified. Thus, we analyzed the degree of biological complexity indirectly, by genetically marking each tumor and testing its capacity to retain the same molecular complexity after transplantation into NSG mice. Cells from MM13, MM16, and MM27 PDX1 tumors were infected with a nontargeting lentiviral library containing ∼12,500 plasmids with unique barcodes (BC; 13K library) and no associated shRNAs, under experimental conditions that allowed high representation of individual BCs (∼400 transplanted cells/BC) and one retroviral integration per cell (see Methods). Transduced cells were then injected into NSG mice to obtain PDX2s. Genomic DNA (gDNA) from PDX1 cells and PDX2 tumors was then analyzed by NGS to assess absolute and relative representations of each BC (Fig. 1A). We found the entire molecular repertoire of BCs in all the analyzed samples (not shown). The relative representation of individual BCs was highly comparable in the two PDX2 tumor replicates (Fig. 1B, left), and the ratio between individual BC reads in the two PDX2 tumors and in the PDX1 cells [log2(ratio)] followed a symmetric distribution, with a median that was centered at around zero (Fig. 1B, right). Thus, tumors were capable of supporting the molecular complexity of a 13K library in vivo, despite their heterogeneous growth properties, with less than 5% of BCs depleted more than 3-fold.
We then screened the three melanomas under the same experimental conditions using a lentiviral shRNA library targeting 236 epigenetic modulators (10 different shRNAs/gene) and 4 screening control genes (a total complexity of 2,410 shRNAs; Fig. 1A and Supplementary Table S2). For all three patients, NGS analyses revealed full representation of the library complexity (≥99% of the 2,410 BCs; not shown) and high correlation (R = 0.81, 0.55, and 0.70; Fig. 1C, left) in PDX2 tumor replicates. The log2(ratio) distribution of the BCs was markedly shifted toward a negative value (medians of −2.72, −2.61, and −0.76, respectively, in patients MM13, MM16, and MM27; Fig. 1C, right), suggesting that the library inhibited melanoma growth in vivo. Notably, in each of the three screens, the positive (PSMA1, KIF11, and RPL30) and the neutral [Renilla luciferase (Luc)] controls behaved as expected (Supplementary Table S2A).
Genes were scored as candidate hits when >6 different shRNAs were found depleted during melanoma growth. Depleted shRNAs were identified by applying the median of the log2(ratio) distribution as cutoff threshold (Fig. 1C, right panel). These analyses led to the identification of 117 genes that were consistently depleted in at least one of the three melanomas (Fig. 1D; Supplementary Table S2A). The depleted genes were equally distributed among the different functional classes of epigenetic modifiers included in the library (Supplementary Table S2B). Analyses of the depleted genes among different melanomas, however, showed a high degree of patient specificity: 66, 69, and 55 genes were counterselected, respectively, in the MM13, MM16, and MM27 melanomas, with only 17 (∼15%) in common (Fig. 1D and Supplementary Table S2A).
Validation of the Epigenetic shRNA Screens
We chose 14 hits with different biological effects in the three samples (Supplementary Table S2A) and different levels of shRNA depletion: (i) strongly depleted [under the first quartile in the corresponding log2(ratio) distribution curve of the epigenetic screen, Fig. 1C; HIT1Q]; (ii) weakly depleted (between the first and the second quartiles; HIT2Q); and (iii) not depleted (no-HIT; Fig. 2A). For their validation, we generated a shRNA library in the pRSI lentiviral vector (80sh library) containing 28 shRNAs against the 14 hits (2 each), 2 shRNAs against control genes (Luc and PSMA1), and 50 scrambled (SCR) shRNAs (see Methods and Supplementary Table S3).
PDX1 cells from the three patients were transduced with the 80sh library and injected in NSG mice. For all patients, we found high correlation of BC frequency in PDX2 tumor duplicates (R = 0.97, 0.99, and 0.99 for MM13, MM16, and MM27, respectively; not shown), and significant depletion in PDX2 tumors versus PDX1 cells of both HIT1Q and HIT2Q shRNAs (Fig. 2A).
Finally, we validated 4 hits (BAZ1B, SMARCA4, CHD4, and KMT2D genes; Supplementary Table S2A; refs. 18–22) by a conventional “single-shRNA” in vivo approach. Three were positive hits in all three melanomas (BAZ1B, SMARCA4, CHD4), whereas the fourth (KMT2D) was positive only in the two NRAS-mutated melanomas.
For each of the four hits and the three tumors, PDX2 cells were independently transduced with pRSI lentiviruses expressing shRNAs against the hits (two each; Fig. 2B), or against Luc and PSMA1. Silencing efficiency was monitored in transduced cells injected into NSG mice (Supplementary Fig. S3A–S3C) and melanoma growth evaluated (Fig. 2B). As expected, the PSMA1 shRNA inhibited tumor growth in all samples, as compared to the Luc shRNA (90%–99% growth inhibition; Fig. 2B). Likewise, the shRNAs silencing BAZ1B, SMARCA4, and CHD4 markedly reduced tumor growth in all three patients (75%–99%; Fig. 2B). The shRNAs silencing KMT2D, instead, inhibited growth of the two NRAS-mutated melanomas (MM13 and MM16), but not of the BRAF-mutated MM27 melanoma (Fig. 2B). Moreover, to ensure that the KMT2D-silenced tumor cells that are still capable of growing in vivo are genetically identical to the original PDX1 cells and to shLuc PDX2 cells, we performed exome sequencing analysis of residual MM13-shKMT2D and MM27-shKMT2D PDX2 tumors and their paired PDX1 and shLuc PDX2 tumor. It is worth noting that residual KMT2D-transduced MM13 and MM27 PDX2 tumors retain the same allelic distribution of the somatic mutations in the corresponding PDX1 and shLuc PDX2 tumors (Supplementary Fig. S2C and S2D), demonstrating that the in vivo subclones' composition of each tumor is also retained after KMT2D silencing and suggesting that KMT2D shRNAs act on the same population of cells. Together, these results provide full validation of the epigenetic shRNA screening.
To investigate the biological mechanisms of the in vivo tumor growth inhibition, shRNA-transduced melanoma cells were analyzed in vitro for their proliferation and migration properties, two key features of the in vivo growth potential of melanomas (23). The effect on proliferation was variable in the different patients, yet all together modest (from 0% to a max of 50% reduction compared to shLuc; Fig. 2C). In contrast, the effect on migration was very strong (50%–90% reduction for all silenced genes), with the expected exception of KMT2D silencing in the MM27 BRAF-mutated melanoma (Fig. 2D). The result was further strengthened by the overexpression of SMARCA4 and CHD4 in short-term cultures of MM16 PDX2 cells (Supplementary Fig. S4A–S4F), which significantly induced cell migration. Together, these results suggest that tumor growth inhibition by silencing of the four epigenetic targets is closely associated with inhibition of cell migration.
KMT2D Activates a Migratory Transcriptional Program in NRAS-Mutated Melanomas
We then investigated the molecular mechanisms underlying the biological effects of KMT2D shRNAs. First, we showed that full depletion of KMT2D by CRISPR/Cas9-mediated targeted deletion in melanoma cells has identical effects on migration of NRAS-mutated melanoma cells (Supplementary Fig. S5A–S5I). Then, we extended our KMT2D shRNA analyses to three additional melanomas, one NRAS mutated (MM23) and two BRAF mutated (MM2 and MM25). KMT2D silencing in PDX2 cells (Fig. 3A) markedly reduced tumor growth in vivo (Fig. 3B) and migration in vitro (Fig. 3C) in MM23, whereas it exerted no effect in MM2 and MM25. As expected, no significant effect on cell proliferation was assessed (Fig. 3D). Thus, KMT2D expression is critical for in vivo growth and cell migration of three NRAS-mutated melanomas (MM13, MM16, and MM23), as compared to the three BRAF-mutated melanomas (MM2, MM25, and MM27). Considering that KMT2D expression was comparable in the six patients (not shown), these data suggest that vulnerability to KMT2D silencing is specific of NRAS-mutated melanomas. Notably, expression of a BRAF-mutated allele in melanoma cells harboring a mutated RAS allele reverted the effect of KMT2D silencing, suggesting that mutant BRAF signals to migration through KMT2D-independent pathways (Supplementary Fig. S6A–S6D).
To characterize the transcriptional program activated by KMT2D in NRAS-mutated melanomas, we performed RNA sequencing (RNA-seq) analyses of parental versus KMT2D-silenced PDX2 cells in the two NRAS-mutated (MM13 and MM16) and the BRAF-mutated (MM27) patients. Comparison of KMT2D-dependent transcription in the three patients revealed 103 genes specifically deregulated in the two NRAS-mutated patients (KMT2D signature; 64 downregulated and 39 upregulated; Fig. 4A and Supplementary Table S4). RNA-seq data were validated analyzing the expression levels of 52 genes (42 deregulated in the two NRAS melanomas; 10 in all three samples). Validation was 100% in the two NRAS melanomas, with a very high correlation between RNA-seq and qPCR quantitative data (R = 0.97 and 0.96 in MM13 and MM16, respectively; Fig. 4B). As expected, the BRAF-mutated MM27 melanoma showed a lower validation rate (88%) and, accordingly, a slightly lower Pearson correlation (R = 0.82).
Ingenuity Pathway Analyses of the 103 genes (IPA, QIAGEN; ref. 24) revealed significant enrichment of genes involved in the regulation of “cell movement” and “migration” (n = 29, P < 0.00001; Supplementary Table S4). A manually curated literature mining enlarged the number of genes involved in migration to 43 (∼42%). A coherent link between deregulation imposed by the KMT2D interference in our melanomas and their effect on migration (i.e., downregulation of genes promoting migration and upregulation of genes counteracting migration) was found in 26 of the 43 genes (7 upregulations and 19 downregulations; not shown). Many of the downregulated genes are overexpressed in different cancers (13/19) or have been demonstrated to be critical for tumor growth in vivo (6/19, e.g., SHC4, AQP1, or RASGRP3 in melanoma; refs. 25–27).
To investigate whether the 103-gene KMT2D signature is specific to NRAS-mutated patients, we analyzed the expression of 52 genes (including 29 migratory genes) in NRAS-mutated and BRAF-mutated melanomas upon KMT2D silencing (Supplementary Table S5). The 52 genes followed the same pattern of regulation in the three NRAS-mutated melanomas, whereas their regulation was nonuniform in the threeBRAF-mutated samples (Fig. 4C). These data demonstrate that KMT2D supports a transcriptional program in the NRAS melanomas, which mainly involves genes responsible for cell migration.
KMT2D Regulates Enhancer Activity in the Patients with Melanoma
KMT2D is a member of the Histone–lysine N-methyltransferase 2 (KMT2) family of proteins that methylate lysine 4 on the histone H3 tail (H3K4) and induce genome accessibility and transcription. KMT2D predominantly promotes H3K4 monomethylation (H3K4me1) at adipocyte-, myocyte-, and macrophage-specific enhancers, and transcriptome changes during adipogenesis and transdifferentiation of preadipocytes into myocytes (20, 22, 28).
Thus, we investigated if the effect of KMT2D on transcription is associated with its monomethyltransferase activity on enhancers, using the MM16 NRAS-mutated melanoma. To map KMT2D genomic sites, shKMT2D- or control shLuc-expressing cells were analyzed by chromatin immunoprecipitation sequencing (ChIP-seq) using anti-KMT2D antibodies. ChIP-seq data revealed 14,258 KMT2D sites in control cells, 6,545 of which were absent (no peak call) in the KMT2D-silenced cells (Fig. 5A). Strikingly, the vast majority of the shKMT2D-sensitive peaks (89%; n = 5,794) mapped to either intergenic regions (∼50%) or gene bodies (∼50%), whereas most of the KMT2D sites present in both control and KMTD2-silenced cells mapped to the transcriptional start site (TSS) of known genes (82%; Fig. 5A). Thus, our experimental conditions of KMT2D depletion (∼50% reduction; Supplementary Fig. S3C) induced a specific loss of KMT2D genomic sites at regions outside gene promoters.
To map KMT2D-bound active enhancers, we first selected genomic regions with H3K4me1-positivity, H3K27Ac-positivity and a high H3K4me1 to H3K4me3 ratio (28). ChIP-seq analyses revealed 20,714 active enhancers in MM16 PDX2 cells (Fig. 5B). Finally, we intersected the KMT2D sites and active enhancers and identified 2,832 KMT2D-bound active enhancers (Fig. 5B). Thus, in MM16 melanoma, ∼50% of the shKMT2D-sensitive KMT2D peaks that map outside the gene promoters were found at active enhancers.
We then investigated the effect of KMT2D expression on H3K4me1 and H3K27Ac levels in the KMT2D-bound active enhancers. Upon KMT2D silencing, the KMT2D-bound active enhancers showed reduction of H3K4me1 (∼50%) or H3K27Ac (∼40%, >1.5-fold reduction) levels (Fig. 5C and D), with only half showing both H3K4me1 and H3K27Ac reduction. This is probably due to the distinct kinetics followed by H3K4me1 and H3K27Ac during enhancer activation/inactivation (29, 30). Thus, to evaluate the effects of KMT2D binding on enhancer activity, we considered, as read out, H3K27Ac levels. Together, these data demonstrate that KMT2D silencing in MM16 cells leads to inactivation of a subset of the KMT2D-bound active enhancers that we named “KMT2D-dependent enhancers” (H3K27Ac >1.5 fold reduction; n = 1,041; Figs. 5D and 6A), and suggest that this might be due to a KMT2D-dependent reduction in H3K4me1 deposition at the same sites. The lack of any effect of KMT2D on the remaining KMT2D-bound active enhancers (those showing no reduction of H3K27Ac, named “KMT2D-independent enhancers,” n = 897; Fig. 6A) might be due to the presence on the same enhancers of other mono-methyltransferases, such as KMT2C, which is also expressed in MM16 cells (data not shown).
Notably, the number of KMT2D-bound and KMT2D-dependent enhancers in MM16 PDX cells did not increase significantly when all the KMT2D peaks were considered and not only the shKMT2D-sensitive ones (3,274 vs. 2,832 and 1,133 vs. 1,041, respectively; Figs. 5B and 6A). Using the same approach, we mapped 6,510 and 4,117 KMT2D-bound active enhancers in MM13 and MM27 PDX cells, respectively (Fig. 6B) and found that the number of KMT2D-bound active enhancers or KMT2D-dependent enhancers in common among MM16, MM13, and MM27 was very low, suggesting a specific KMT2D-enhancer-regulation pattern among different tumors.
KMT2D Deregulates Specific Enhancers and Target Genes in NRAS Melanomas
We then investigated whether the effect of KMT2D on enhancer activity is mechanistically linked to its effects on gene-specific transcription. Enhancers influence expression of their targets over large distances (tens to several hundreds of kilobases), and different mechanisms—DNA looping, tracking/scanning of intervening sequences—have been proposed to explain how they associate with target TSSs (31). As an initial assessment of the relationship of TSSs of KMT2D downregulated genes with KMT2D-dependent enhancers, we measured their physical distance in linear genomic sequence, and compared this value to the distance between TSSs of KMT2D downregulated genes with KMT2D-independent enhancers in all three melanomas. We found that the KMT2D-dependent enhancers are in greater proximity to the closest downregulated genes (Fig. 6A and C), suggesting that KMT2D-dependent modifications at chromatin of distal enhancers are mechanistically linked to variations in expression of target genes.
To test this hypothesis, we investigated whether patterns of enhancer activation by KMT2D correlate with selectivity of its transcriptional effect. To this end, we first investigated whether the TSSs of 64 NRAS-specific downregulated genes (Fig. 4A) were closer to the KMT2D-dependent enhancers in the NRAS (MM13 and MM16)- versus BRAF (MM27)-mutated melanomas. Twenty-nine and 45 genes scored as the closest to KMT2D-dependent enhancers in MM16 and MM13, respectively (median distance of 303.6 and 72.7 Kb; Fig. 6D). Strikingly, the distances between the same genes and the closest KMT2D-dependent enhancers in MM27 (2.9 and 1.4 Mb; Fig. 6D) were 10- and 20-fold greater than those in MM16 and MM13, respectively (Fig. 6D), suggesting that these enhancers might not be involved in the regulation of the NRAS-downmodulated genes. We then investigated the effect of the binding of KMT2D to the identified KMT2D-dependent enhancer sites on the expression of the closest putative target genes. For this purpose, we selected two KMT2D-dependent enhancers showing strong H3K27Ac drop (2.2-fold reduction; E1: 18 kb upstream of the MFGE8 TSS, and E2: 113 kb upstream of the RPL39L TSS; Fig. 7A and B) among the closest to the 29 MM16 downregulated genes (Fig. 6D), and one control enhancer (E3; Fig. 7C) mapping 24 kb upstream of ITPKB TSS. This enhancer, in fact, is KMT2D bound but shows no significant H3K27Ac reduction upon shKMT2D (Fig. 7D). To inhibit KMT2D activity at the selected enhancer regions, we used the dCas9–KRAB fusion protein with two different single guide RNAs: sgRNA #1a or #1b, #2a or #2b, #3a or #3b, targeting E1, E2, and E3, respectively (Fig. 7D). Krüppel-associated box (KRAB) has been reported to efficiently silence transcription by recruiting KAP1 and HP1 proteins (32) and allows genome-specific targeting when fused to catalytically inactive Cas9 (dCas9). Short-term cultures of MM16 PDX2 cells were independently infected with sgRNAs and then transduced with dCas9-KRAB construct. Strikingly, four days after dCas9–KRAB transduction, targeting of the E1 or E2 enhancers significantly reduced the expression of MFGE8 and RPL39L, respectively, whereas no ITPKB downregulation was observed upon targeting of E3 (Fig. 7D). MFGE8 and RPL39L downregulation was also observed upon targeting of the KMT2D-dependent enhancers E4 or E5 (mapping farther away from MFGE8 and RPL39L TSS) with the corresponding sgRNAs (#4a or #4b and #5a or #5b, respectively; Fig. 7D). Together, these data show that specificity of the transcriptional effect of KMT2D correlates with physical proximity of KMT2D-dependent enhancers to the target genes.
We demonstrated that in vivo genetic screens of patients' tumors are feasible, at least using shRNA libraries of relatively high complexity (2,410 shRNAs) and metastatic melanomas. Melanomas are among the tumors with the highest frequency of TICs, which is a critical limiting factor of the in vivo screens. Our screening protocol, however, can sustain the high subclonal complexity of melanomas, suggesting that under the appropriate experimental conditions, biological complexity may not limit in vivo genetic screens of primary tumors.
In vivo genetic screens allow the identification of genes indispensable for cancer maintenance in patient-derived tumors under in vivo conditions. This is significantly different from the currently used model systems, mainly cancer cell lines (5–14), and might provide novel insights into mechanisms of tumor maintenance.
Our screens revealed unprecedented numerosity and unexpected degree of functional heterogeneity among individual tumors. We screened 236 genes in three patients with melanoma and identified around 60 critical genes per tumor, with only 17 (∼15%) in common and a total of 117 involved (e.g., ∼50% of all tested genes). We also compared the results obtained in our in vivo shRNA screens with in vitro genome-wide RNAi screens performed on A2058 (BRAFV600E), COLO783 (BRAFV600E), SK-MEL-5 (BRAFV600E), and HS944T (NRASQ61K) metastatic melanoma cell lines (see details in Supplementary Methods; ref. 33). We reanalyzed the published in vitro dataset, designating as hits those genes whose shRNA depletion has an ATARiS gene score <-1 corresponding to a gene depletion of 2-fold (see details in Supplementary Methods). We focused our analysis on the epigenetic genes present in the in vivo and in vitro screens (a total of 109 genes) and we found: (i) 39 of 109 (∼36%) positive hits in the genome-wide in vitro screens (Supplementary Fig. S7A), and 58 of 109 (∼53%) in our in vivo screen of PDX melanomas (Supplementary Fig. S7B); (ii) no common hits among the four cell lines, and a low frequency of common hits among the three PDX melanomas (8 of 58); (iii) a total of 22 (20%) common hits among cell lines and PDX tumors, with few common hits when comparing single PDX tumors with the melanoma cell lines (Supplementary Fig. S7C–S7E). The comparison of our in vivo genetic screen in melanoma PDXs with the available in vitro screen done in melanoma cell lines showed the same high heterogeneity of epigenetic genes critical for tumor growth or cell viability. We think that this is unlikely due to a peculiar role of epigenetic genes in the regulation of the transformed phenotype. Indeed, the same melanomas were screened with shRNA libraries of metabolic genes and deubiquitinating enzymes/helicases (335 and 287 genes, respectively), which gave similar results (data not shown). Alternatively, the observed functional complexity and patient heterogeneity might reflect multiple layers of adaptation of individual tumors to the continuously changing tumor environment or genomic context, with consequent activation of multiple and nonredundant signaling pathways and high frailty of the transformed phenotype. Notably, the identified melanoma hits do not appear to be activated by somatic mutations in our patients, despite the fact that epigenetic genes are frequently mutated in cancer, including melanomas (34, 35). The frequency of SNVs in the identified hits, in fact, was slightly lower than in the nonhits (Supplementary Figs. S8A and S8B).
Genes carrying biologically relevant mutations (cancer driver genes; ref. 2) are considered the best candidates for the development of targeting drugs. Our data, however, suggest that somatically mutated genes are not necessarily the most critical genes for tumor maintenance, because virtually all the identified hits in our melanomas were not somatically mutated. This might have important clinical implications, because it would significantly expand the pool of druggable genes, and, for each patient, the number of critical genes for which a drug is available (actionable genes). For example, 10 of the melanoma essential genes that we have identified are indeed actionable (Supplementary Fig. S8C). This is particularly relevant for the NRAS-mutated melanomas, which currently have few treatment options (36).
We investigated the biological and molecular mechanisms in tumorigenesis of KMT2D, one of the NRAS melanoma essential genes that we have identified. KMT2D is frequently mutated in a variety of cancers (20). Most cancer-associated KMT2D mutations are frameshift and nonsense alterations, suggesting that KMT2D functions as a tumor suppressor. Several lines of evidence, however, suggest that KMT2D might also exhibit oncogenic properties for various tumors. It is overexpressed in breast and colon carcinomas, where it is associated with poor prognosis, whereas its silencing significantly reduces migration in colorectal and breast cancer cell lines and growth in a mouse xenograft of bladder cancer (20). We showed here that KMT2D depletion reduces cell migration and inhibits in vivo growth of NRAS-mutated melanomas. Notably, KMT2D was found in its germline configuration in the same patients, suggesting alternative mechanisms of KMT2D activation in these tumors.
We have initially investigated the consequences of KMT2D activation in NRAS-mutated melanomas. KMT2D is the major candidate enhancer H3K4me1 methyltransferase in mammals (22). We showed that KMT2D silencing leads to inactivation of a subset of KMT2D-bound enhancers (reduced H3K4me1 and H3K27Ac) and downregulation of a subset of genes that are critical for cell migration, including MFGE8 and RPL39L. Downregulation of either gene in several cancer types, and in melanoma and breast carcinomas specifically, reduces cell migration, and their expression correlates with a more aggressive phenotype and unfavorable outcome in the patients (37–39). Notably, KMT2D target genes were the most proximal to the KMT2D-dependent enhancers, suggesting that KMT2D promotes tumorigenesis by deregulating enhancer activity.
Patterns of gene expression and enhancer deregulations, however, showed great heterogeneity in the 3 analyzed patients (2 NRAS mutated and 1 BRAF mutated), even in the 2 NRAS mutated where KMT2D functions as melanoma essential gene. We identified ∼3,000 KMT2D-regulated genes, with very few genes common to the 3 patients (∼5%). The overlap in the 2 NRAS-mutated patients was ∼11%, even lower than that observed between each of the 2 NRAS melanomas and the BRAF melanoma (∼21%; Fig. 4A). Likewise for the KMTD2-dependent enhancers: 2,803 in total, 0.07% common to the 3 patients (2/2,803), 1.5% common to the 2 NRAS-mutated patients (35/2,418), 4.7% (71/1,526, MM16/MM27) and 0.6% (10/1,774, MM13/MM27) in common between the NRAS-mutated and BRAF-mutated melanomas (Fig. 6B). Even when analyses of enhancer gene patterns were restricted to the 69 genes upregulated by KMT2D in both NRAS patients, the closest KMT2D enhancers differ in the 2 patients (not shown), suggesting that assembly of enhancer genomic patterns by activated KMT2D is highly tumor specific, regardless of the identity of the selected transcriptional targets. A central feature of enhancers is their ability to function as platforms for the recruitment of multiple transcription factors, thus ensuring integration of multiple intrinsic and environmental signaling pathways. We speculate that activation of KMT2D in the two NRAS melanomas follows completely distinct upstream signaling pathways, which might reflect cellular responses to tumor-specific environmental or genetic context.
NSG male mice (6–12 weeks age; 15–25 g weight) were purchased from The Charles River Laboratories. In vivo studies were performed after approval from our fully authorized animal facility and notification of the experiments to the Ministry of Health (as required by the Italian Law; IACUCs N° 02/2012 and N° 758/2015-PR), and in accordance with EU directive 2010/63.
The in vivo mouse model generated in the present study was similar to that previously described (17), with few modifications. Briefly, tissue biopsies of metastatic melanomas were collected from patients whose informed consent was obtained in writing according to the policies of the Ethics Committee of the European Institute of Oncology and regulations of the Italian Ministry of Health. The studies were conducted in full compliance with the Declaration of Helsinki. Tumors were mechanically dissociated and subsequently digested with an enzymatic combination of Collagenase Type III (1 mg/mL, Worthington Biochem) and Dispase (0.5 U/mL, STEMCELL-Technologies) for 45 minutes at 37°C. After incubation, cells were treated with red blood cell lysis buffer (155 mmol/L NH4Cl, 12 mmol/L NaHCO3, 0.1 mmol/L EDTA) to remove erythrocytes, and then filtered (40-μm cell strainer) to obtain a single-cell suspension. Melanoma cells obtained from dissociation of lymph node metastases were separated from leukocytes by Magnetic Activated Cell Sorting, using CD45 microbeads (Miltenyi Biotec 130-045-801) and LD columns (Miltenyi Biotec 130-042-901), reaching an enrichment of more than 90% (as assessed by FACS). To generate primary PDX1s, 100,000–500,000 dissociated cells were resuspended in a 3:1 mix of L15 medium and Matrigel Matrix (Corning 354248) and subcutaneously injected into the flank of NSG mice. Tumor formation was monitored weekly, and tumor diameters measured with calipers. Mice were sacrificed when tumors reached the volume of ∼0.5 cm3. PDX1 tumors were dissociated as patient biopsies, and either serially transplanted to obtain secondary and tertiary PDXs (PDX2s and PDX3s) or snap-frozen. Tissue biopsies and PDXs were characterized by immunohistochemistry analysis (see below). Purity of PDX-dissociated human melanoma cells (≥95%) was evaluated, analyzing HLA (BD555555) expression by FACS Aria.
Libraries and Plasmids
All libraries and plasmids were purchased by Cellecta (Cellecta Inc.) and assembled into the pRSI vector backbone, containing the puromycin-resistance marker and the GFP fluorescent reporter. shRNAs were under the control of a constitutive U6 promoter and univocally associated with a barcode cassette of 18 degenerated, nonoverlapping nucleotides. The nontargeting library (13K library) consists of 12,415 vectors, each carrying a unique BC and no associated shRNA. The epigenetic shRNA library (see Supplementary Table S2) contains 2,410 vectors, targeting 236 genes (10 different shRNAs per targeted gene), 3 positive controls (PSMA1, KIF11 and RPL30), and 1 neutral control (Luciferase). Eighty pRSI vectors, each one expressing an shRNA associated with a unique BC (see Supplementary Table S3), were purchased by Cellecta and either pooled together to obtain a small-scale shRNA library (80sh library) or used individually. Fifty SCR shRNAs were synthesized starting from sequences of the following genes: E. coli lactose operon (GenBankJ01636.1), trp operon (GenBankV00372.1), and kanamicin resistance (GenBankAJ002684.1); S. aureus chloramphenicol resistance (GenBankAB481130.1); vector pHV1249 ampicillin resistance (GenBankAF307748.1); Anabaena nitrogen fixation (GenBankJ05111.1); P. abies rubisco large subunit (GenBankX75478.1); C. arabica rubisco large subunit (GenBankAJ419827.1); S. cerevisiae genes ACH1 (Gene ID 852266), TPS1 (Gene ID 852423), PNC1 (Gene ID 852846), ILV1 (Gene ID 856819), and YMC2 (Gene ID 852401). Each gene was divided in a pool of 21-bp oligonucleotides using the sliding window scan through the sequence. All sequences were aligned against human (hg18) and mouse (mm9) genomes using two aligner programs, Bowtie (40) and blastn. SCR shRNAs were selected among those sequences that did not align against human and mouse genomes using both aligner programs.
Cell Culture and Infection
The A375 and SK-MEL-30 melanoma cell lines were obtained from IZSBS in 2001 and from DSMZ in 2002, respectively. Cells were routinely tested (every two months) and authenticated by Gene Print 10 System Promega in 2015 (A375) and 2011 (SK-MEL-30). A375 were cultured in DMEM 10% FBS and SK-MEL-30 in RPMI 10% FBS. Media were supplemented with 200 mmol/L glutamine, 100 U/mL penicillin, and 100 μg/mL streptomycin. PDX cells were maintained in Iscove's modified Dulbecco's medium (IMDM) supplemented with 200 mmol/L l-glutamine and 10% FBS. Concentrated lentiviral particles (TU, transducing units) from libraries or single plasmids were either purchased by Cellecta or produced by transfecting 293T cells, as described in the Cellecta User Manual. Lentiviral particles, together with 4 μg/mL polybrene, were added to the PDX cell standard medium for 16 hours. Forty-eight hours after infection, the medium was replaced and puromycin (2 μg/mL) added for 3 days. For library infection, PDX1 cells were infected using a multiplicity of infection (MOI) =∼0.2 TU/cell. Conversely, in the in vivo validation, in vitro studies, RNA-seq and ChIP-seq experiments, PDX2 cell were infected at an MOI of ∼3, with single or pooled shRNAs silencing specific target genes (see the main text).
In Vivo shRNA Screening
Five million 13K library or 1 million epigenetic library transduced cells were injected subcutaneously in duplicate or triplicate, respectively, in NSG mice. 80sh library transduced cells (400,000) were transplanted in duplicate. Tumors were harvested when they reached ∼0.5 cm3 in volume. BC representation was determined by NGS using the Illumina/Solexa platform (HiSeq 2000). NGS libraries were obtained according to the Illumina manual. A BC-specific, sequencing primer was utilized. BCs were identified by aligning each sequencing read to the 13K library and to the barcoded libraries using the Bowtie aligner (40), and by considering only those BCs having, at most, three mismatches for each alignment. BC frequency was calculated in PDX1 cells (fc) and PDX2 tumors (ft) by dividing each BC count by the total number of aligned reads. For each BC, log2(ratio) was calculated as the base 2 log of the ft to fc ratio. Details of gDNA extraction, nested PCR condition, and NGS are reported in the Cellecta User Manual.
In Vivo Validation
For the individual in vivo validation, 100,000 PDX2-infected cells were transplanted in quadruplicate in NSG mice. Each gene was silenced by two different shRNAs. Tumor volume was calculated using the modified ellipsoid formula: 1/2 (length × width2).
In Vitro Studies
PDX2-infected cells (5,000) were incubated for 72 hours, and proliferation was measured by CellTiter-Glo assays (Promega). The migration assay was performed using 8.0-μm pore size, fibronectin precoated inserts in 24-well plates. Duplicates of 100,000 cells in serum-free medium were plated in the upper chamber and complete medium was added to the lower chamber. After 36 hours of incubation, cells that had migrated to the lower surface of the inserts were stained with 0.5% crystal violet. Five images of each insert were acquired and analyzed with ImageJ software as described (41).
gDNA of patients' samples was extracted from frozen tissues, containing at least 85% melanoma cells (MM13 and MM27), or formalin-fixed, paraffin-embedded (FFPE) tissues, containing around 60% of neoplastic cells (MM16). gDNA was also obtained from matched patients' nontumor tissues (normal counterpart) and paired xenograft (PDX1 and PDX2) tissues. gDNA was prepared and whole-exome sequencing performed according to standard protocols (see Supplementary Methods). Sequencing alignment to human genome and subsequent bioinformatic analysis is fully described in Supplementary Methods.
Total RNA extracted from PDX2 infected cells was purified and libraries sequenced on an Illumina HiSeq2000. Bioinformatic analysis is fully described in Supplementary Methods.
ChIP Sequencing and Bioinformatics Analysis
PDX2-infected cells were subjected to ChIP. Briefly, after cross-linking, cells were lysed, sonicated, and incubated with specific antibodies (see Supplementary Methods for details). Samples were sequenced and aligned to human genome. Bioinformatic analysis is fully described in Supplementary Methods.
Tissue fragments from patient biopsies or PDXs were formalin fixed and paraffin embedded. After deparaffinization, sections were treated according to standard protocols, as detailed in Supplementary Methods.
Enhancer Repression by CRISPRi Technology in PDX Cells
The experiments were performed using the CRISPRi technology (31). Detailed experimental protocol is indicated in Supplementary Methods.
Total RNA was extracted with the Qiagen RNeasy Mini Kit, and reverse transcribed using random hexamers (Improm-II, Promega). RNA expression of KMT2C, KMT2D, and BAZ1B was determined by qRT-PCR using the Fast SYBR Green Master Mix (Applied Biosystems) and the iCycler iQ real-time detection system software (Bio-Rad). RPLP0 was used as a housekeeping gene. RNA expression of KMT2D-regulated genes was determined by TaqMan Gene Expression assays (Applied Biosystems), using single tube assays or 384-well microfluidic cards. Up to 8 different housekeeping genes were utilized as normalizers. Three biological replicates were included for each experiment. Primer and TaqMan assay details are available upon request.
Sequencing data have been submitted to the National Center for Biotechnology Information (NCBI) Gene Expression Omnibus (GEO) under accession N° GSE71854.
Disclosure of Potential Conflicts of Interest
No potential conflicts of interest were disclosed.
Conception and design: D. Bossi, A. Cicalese, A. Testori, G.F. Draetta, S. Minucci, P.G. Pelicci, L. Lanfrancone
Development of methodology: D. Bossi, A. Cicalese, G.I. Dellino, L. Luzi, C. D'Alesio, G.R. Diaferia, A. Carugo, E. Cavallaro, M. Barberis, S. Punzi, I. Pallavicini, T.P. Heffernan, L. Lanfrancone
Acquisition of data (provided animals, acquired and managed patients, provided facilities, etc.): D. Bossi, A. Cicalese, G.I. Dellino, A. Carugo, R. Piccioni, M. Barberis, A. Testori, G. Tosti, L. Lanfrancone
Analysis and interpretation of data (e.g., statistical analysis, biostatistics, computational analysis): D. Bossi, A. Cicalese, G.I. Dellino, L. Luzi, L. Riva, A. Testori, L. Giacó, G. Melloni, S. Minucci, P.G. Pelicci, L. Lanfrancone
Writing, review, and/or revision of the manuscript: D. Bossi, A. Cicalese, G.I. Dellino, L. Luzi, G. Tosti, G. Natoli, S. Minucci, P.G. Pelicci, L. Lanfrancone
Administrative, technical, or material support (i.e., reporting or organizing data, constructing databases): D. Bossi, A. Cicalese, L. Luzi, L. Riva, L. Lanfrancone
Study supervision: D. Bossi, A. Cicalese, P.G. Pelicci, L. Lanfrancone
Other (immunohistology and microscopy analysis): G. Mazzarol
The authors thank M. Varasi for insightful discussion, A. Gobbi and M. Capillo for excellent support in animal work, L. Rotta for excellent sequencing support, and G. Giardina, C. Spinelli, A. Papait, D. Di Gesto, F. Cataldo, I. Davidson, and E. Pasqualucci for providing reagents and protocols. They also thank all members of the Department of Experimental Oncology for discussion and reagents. They thank the Genomic Unit (IEO), the Mouse Facility (Cogentech), the DNA service (Cogentech), and the Cell Biology Unit (IEO). PierGiuseppe Pelicci and Luisa Lanfrancone are members of the EurOPDX Consortium.
This work was supported by the European Research Council Advanced Grant 341131 and the Italian Association for Cancer Research Investigator Grant 14216.