Co-existence of BRAF and NRAS driver mutations in the same melanoma cells results in heterogeneity of targeted therapy resistance

INTRODUCTION

The MAPK pathway, consisting of RAS-RAF-MEK-ERK, is a highly conserved signaling cascade in eukaryotic cells conserved from yeast to humans with many vital cellular functions, such as proliferation, differentiation, migration, and apoptosis [1]. About one-third of cancers have a deregulated MAPK pathway, either due to overexpression of receptor tyrosine kinases (RTKs), increased production of activating ligands, activating mutations in RTKs, RAS or RAF or to failure of pathway control mechanisms [1]. In cutaneous melanoma, deregulation of the MAPK pathway is mainly caused by a hyperactive mutation in BRAF (50% of cases) or NRAS (15% of cases), highlighting the important role of controlled MAPK signaling for melanocyte homeostasis [1, 2].

Targeting a hyperactivated MAPK pathway driven by mutated NRAS or BRAF with specific BRAF- and MEK inhibitors, increases the median overall survival from metastasized melanoma patients from 9 months with no therapy to approximately 14 months with successful inhibitor treatment [3]. Unfortunately, resistance to MAPK inhibition almost invariably develops [4]. Several resistance mechanisms have been described so far, which can roughly be divided into those that reactivate the MAPK pathway by circumventing the inhibitory effects of the MAPK inhibitors, or ones that activate alternative signaling pathways [5].

In the case of BRAF inhibitors, Shi et al identified reactivation of the MAPK pathway (70% of cases), mostly in the form of additional NRAS or KRAS mutation (18% and 7% of cases, respectively), CDKN2A loss (7% of cases), mutant BRAF amplification (19% of cases) or BRAF alternative splicing (13% of cases) as the most common resistance mechanisms. They also identified the PI3K-PTEN-AKT pathway as the second important resistance pathway (22% of their post-treatment samples contained mutations in PI3K-AKT regulatory genes) [5].

One of the more prevalent mechanisms is an additional mutation in NRAS, leading to reactivation of the MAPK pathway [6]. However, it has been published that NRAS and BRAF mutations are mutually exclusive in single cells due to self-induced apoptosis by sustained hyper-activation of the MAPK pathway [7, 8]. Consequently, resistant tumors of patients that contain both mutations concurrently may be comprised of several mutually-exclusive subclones with either the activating BRAF or NRAS mutations [7]. A recent paper showed that both mutations can co-occur in a small area (of approximately 10,000 cells) selected by laser microdissection [9], although this does not prove that the mutations can co-occur within single cells. Likewise, although double-mutated NRAS/BRAF melanoma cultures have been previously reported, these may still represent heterogeneous mixtures of singly-mutated melanoma cells [10, 11, 12], or may have arisen artificially through in vitro drug treatment and selection experiments [13].

Within a patient, various small populations of tumor cells (i.e., subclones) evolve during disease progression, which exhibit different genotypes and/or phenotypes ([14, 15, 16]). Due to these different tumor subclones within a patient (intra-patient heterogeneity), it is believed that diverse resistance mechanisms can co-exist within one patient [17, 5]. However, where these resistance mechanisms originate from and how they evolve under treatment remains poorly understood [6, 11].

To better characterize the evolution of intra-patient heterogeneity under different treatment regimens, we performed exome sequencing on multiple samples from 3 stage IV melanoma patients (cohort 1) who each received a different therapy (BRAF inhibitor (patient 1), MEK inhibitor (patient 3) or multi-receptor tyrosine kinase (patient 2)) but progressed quickly under treatment. We used formal phylogenetic methods on tumor DNA to model the evolution of intra-patient heterogeneity from primary tumors to each individual metastasis for the targeted and non-targeted therapies. In addition, we could detect the presence of an NRAS mutated subclone in 1 of 5 treatment-resistant tumors from the BRAF inhibitor resistant patient (patient 1). Single cell clone sequencing from the cell culture generated from this treatment resistant tumor revealed the co-occurrence of a BRAF and NRAS mutation in a single cell. This was confirmed in an additional group of 5 patients (patient 4, 5, 6, 7 and 8, cohort 2) from whom tumors after BRAF inhibitor treatment were collected and where cell cultures were generated from these tumors and showed secondary NRAS mutations. Sequencing of 65 clonal populations derived from 4 of these cell cultures showed the presence of both activating MAPK mutations in all but one subclone. Further in vitro work with these double-mutated cell cultures demonstrated sensitivity to BRAF inhibition, but heterogeneous responses to downstream MAPK inhibition, as well as to PI3K pathway inhibitors and the multi-receptor kinase inhibitor Pazopanib.

RESULTS

Tumor-type dependent, intra-patient heterogeneity

We sequenced whole exomes of 27 samples from three metastatic melanoma patients (cohort 1) with different mutational statuses and different treatments (Table 1 and Figure 1A–1F). For all patients we performed exome sequencing on all of their samples and confirmed their mutational status (BRAF^V600E mutated, BRAF^WT/NRAS^WT or NRAS^Q61K mutated for patient 1, 2 and 3, respectively) (Table 1, Figure 1, Supplementary. Table S1). In addition to these driver mutations, we looked for other mutated onco- and tumorsuppressor genes: in patient 2 we identified a nonsynonymous germline mutation in the Melanocortin receptor MC1R^V92M and in patient 3 a nonsynonymous germline mutation in the Microphthalmia-associated transcription factor MITF^E318K (data not shown). By using EXCAVATOR and CONTRA algorithms we detected a high number of copy number variations (CNVs) in many chromosomes, with some samples exhibiting large losses throughout the genome (Figure 1G–1I).

Whole-exome phylogenetic analysis identifies monophyletic evolution of therapeutic resistance

In order to investigate the evolutionary relationship between tumor sites within individual patients, we used phylogenetic algorithms to group tumor samples based on their total single nucleotide variants (SNVs) and copy number variations (CNVs) (Figure 2, Supplementary Table S2). As opposed to the phylogenetic tree of patient 2 biopsies, the evolutionary trees from patient 1 and 3 showed monophyletic clades for the post-resistance samples (i.e. late metastases), meaning that these originated from only one node (Figure 2). However, no known and shared mechanism of resistance to BRAF-inhibitor or MEK-inhibitor treatment could be identified (Supplementary Table S2). Confidence is shown by bootstrap supports (arrow) which reflects the percentage of bootstrap trees also placing the clade at the endpoints of that branch. The tree of patient 2 (BRAF^WT and NRAS^WT, patient received non-targeted therapy) did not show this strong monophyletic clade of late tumor metastases, instead the post-resistant samples originated from multiple nodes (Figure 2, arrows). Since there might be multiple different resistance mechanisms present in one patient, we also sought to identify explanatory protein-coding changes in any of the post-treatment samples. However, no known mechanisms of resistance were identified in the exome data of any tumor in the three patients except for the NRAS^Q61K mutation that was present in a cell line (MM121224) derived from a resistant metastasis from patient 1 (Figure 2).

Resistant subclones are present heterogeneously and at low frequencies

In order to determine the origin of this NRAS^Q61K mutation from patient 1 we performed digital PCR as well as ultra-deep sequencing on all tumor samples (data not shown, and Figure 2D). Only one cell culture (MM121224) and one biopsy where the cell culture was derived from had an activating mutation in exon 2 of the NRAS gene (NRAS^Q61K), with a range from 5,473× coverage of the NRAS exon to 26,416× coverage by next generation sequencing (Figure 2D, and Supplementary Table S3). No other sample from this patient had this activating mutation, suggesting other resistance mechanisms to be involved in the metastases of the same patient. Using the program deepSNV, a diverse series of other NRAS coding mutations in exon 2 was present significantly in all samples at low subclonal frequencies compared to the germline blood reference from the same patient (Figure 2D and 2E, Supplementary Table S3). Interestingly, the allele carrying the NRAS^Q61K mutation was only present at a frequency of about 6% in early passage cultures of the MM121224 cell line from this patient (Figure 2D–2F).

Two activating MAPK mutations are present in single melanoma cells

To determine the frequency of double mutations in early passage cultures in general in a larger patient cohort, we Sanger sequenced the NRAS locus in all BRAF-mutated cell cultures generated from 2013—2015 archived in our melanoma biobank [18]. In this way, we identified an additional 8 double-mutated cell cultures, derived from 7 different patients, bringing the total to 9 out of 122 cell cultures (7.4%). Of these 9 cell cultures, 7 had an activating BRAF mutation at position 600 and a co-occurring activating NRAS mutation at position 61 or 12 (Table 1). Two of these cell cultures had a BRAF mutation at another position, namely one culture contained the double mutation BRAF^D594H/NRAS^Q61R and one culture had the double mutation BRAF^L597R/NRAS^Q61R (Data not shown). The BRAF^D594H is an inactivating mutation, and a double mutation of this sort is already described elsewhere [19]. The BRAF^L597R is a less prevalent activating mutation for which less information is available. We therefore decided to focus on the six patients (patients in cohort 2, and patient 1 from cohort 1, see Table 1) with cell cultures that had a BRAF mutation at position 600 and a co-occurring NRAS mutation at position 61 or 12.

We asked if the double-mutated cell cultures consisted of two exclusive populations of cells (one with BRAF and another with NRAS mutations), or if both mutations were present in single melanoma cells. To distinguish between these possibilities, we generated single-cell clones of 4 double mutated cell cultures (from patients 1, 5, 7 and 8) and confirmed by Sanger sequencing the presence of both BRAF^V600E/K and NRAS^{Q61K/H/R or G12A} mutations in 65 independently derived colonies (Table 1).

Both alleles (mutated and WT) from BRAF and NRAS were found to be expressed with Sanger sequencing of cDNA and RNA-seq (data not shown).

Double mutations occur heterogeneously within patients after targeted therapy

In order to investigate the evolution of the NRAS/BRAF double mutated cancer cells within a patient, we analyzed tumor DNA for the presence of double mutations in all available histological and frozen material from patients 4, 5, 6, 7 and 8 before and after BRAF inhibitor treatment (Table 1). For patients 4, 5 and 6 we could confirm the presence of the additional NRAS mutation in the post-treatment samples from which the cell culture was derived. In patient 8, we had generated cell cultures from Late met 1 and Late met 2. In Late met 1 we could detect the BRAF^V600E NRAS^Q61R double mutation, however in Late met 2 we could not detect the additional NRAS^Q61K mutation. For patient 7, the histology block corresponding to the tumor from which the cell culture was derived was no longer available. In patient 4 and patient 8, we could also detect the NRAS/BRAF double mutation in additional post-treatment metastasese. Interestingly, the NRAS mutation could not be detected in any of the pre-treatment biopsy samples.

Double-mutant cells have heterogeneous MAPK pathway inhibitor and alternative pathway responses

In order to gain more insight into the biology of NRAS/BRAF double mutated cells, we performed viability assays with (control) single and double mutated cell cultures under MAPK inhibitor treatment. Double-mutated cell cultures from all 6 identified patients were resistant to three different BRAF inhibitors (Figure 3A, Table 2). The response to MEK and ERK inhibitors was heterogeneous. The cell cultures MM140906 and MM150850 were resistant to MEK and ERK inhibitors, whereas other cell cultures were sensitive or partially resistant for one or for both of the inhibitors (Figure 3A, Table 2). We also tested the response for different inhibitors of the PI3K-AKT pathway, as this pathway is often involved in MAPK pathway inhibitor resistance [5], and for the multi-RTK inhibitor Pazopanib. Here the cells also showed a heterogeneous response, except for the multi-RTK inhibitor, for which all cell cultures were resistant (Figure 3A, Table 2). Therapies combining a MEK inhibitor with different PI3K pathway inhibitors worked synergistically in all cell cultures, albeit with different strengths (Table 2). As it was previously thought that BRAF and NRAS mutations are mutually exclusive due to the growth disadvantage of double mutated cells [8], we analyzed the proliferation rate of double-mutated cells in vitro, compared to single mutated control cell cultures (Figure 3B). Although the cell cultures MM121224 and MM140307 showed a higher proliferation rate compared to the control cell cultures, the other double-mutated cells had reduced proliferation rates.

As a downstream read-out for BRAF and/or MEK activation, western blot analysis for total ERK and phosphorylated ERK was performed in the presence of BRAF, MEK, or ERK inhibition. This showed that basal pERK levels were mostly higher in the double mutated cells compared to the single NRAS mutated control cell cultures (but not in MM150423) (Figure 3C). Upon treatment with the BRAF inhibitor LGX818 (encorafenib), the pERK levels of MM121224, MM150423, MM150849 and MM140906 stayed the same compared to the untreated control, whereas the levels of MM140307, MM150850 and M130903 decreased (Figure 3C). Levels of pERK expression upon MEK162 (binimetinib) treatment decreased in all cell cultures except for MM140906. Upon ERK inhibition, pERK levels were reduced in all cell cultures.

When we compared pERK levels and pAKT levels between the different double mutated cell cultures without treatment, we found that MM140307, MM140906, MM150849 and MM150850 expressed relatively high levels of pERK, whereas M130903, MM150423 and MM121224 express relatively low levels (Figure 3D). pAKT expression is relatively constant between the different cell cultures, except MM150849, which has a relatively high expression (Figure 3D)

DISCUSSION

Genetic or transcriptional heterogeneity in tumors is a major obstacle to obtaining durable responses to targeted therapy for metastatic melanoma. In order to better understand how individual cancer patients respond to standard therapies, we conducted multiple-sample, whole-exome sequencing from multiple time-points in 3 patients receiving different therapeutic regimens.

The sequencing results were used to infer the evolutionary relationships between the tumors within each patient, and to determine how each therapeutic regimen affected the evolution of genetic heterogeneity. Unlike previous studies that showed a branching evolution of clones subsequent to targeted therapy, we could see a strong, well-supported monophyletic evolution of metastases following both BRAF and MEK inhibitor treatment and relapse with phylogenetic analysis [5]. In contrast, patient 2, who received a multi-kinase inhibitor (i.e. pazopanib), did not have a monophyletic topology of late tumor metastases, which is suggestive of genetic drift between the late metastases. Despite the monophyletic segregation of late metastases in the patient who received the BRAF inhibitor, no known genetic mechanism of resistance was shared between all sequenced biopsies accounting for the inter-patient heterogeneity and subsequent treatment difficulties. In fact, the additional activating mutation in NRAS^Q61K found in patient 1 with a BRAF mutation background was only present in a single metastasis of patient 1 and absent in all other resistant tumor samples from that patient. This is consistent with previously published data showing heterogeneity in resistance mechanisms within individual patients [5].

We went on to check if the double mutation could also be found in the cell culture obtained from this resistant metastasis, and if these mutations occurred in the same cell. By isolating and sequencing colonies derived from 23 single-cell clones of the resistant late metastasis 6 from patient 1, we could show for the first time that both activating MAPK mutations (NRAS and BRAF) were present in a single tumor cell.

In order to confirm our finding, we identified a second cohort of BRAF inhibitor resistant patients from whom we had obtained a double mutated cell culture. In screening 121 cells from our live-cell biobank, we could identify an additional 8 cell cultures with double BRAF/NRAS activating mutations, bringing the total to 9 out of 122 cell cultures. By sequencing colonies derived from single-cell clones, we confirmed the presence of the double mutation in 64 out of 65 colonies. However, in our standard biobank protocol we establish cell cultures without additional treatment, also when they are derived from a therapy resistant patient. Therefore, due to selection, this could be an underrepresentation of the actual number of double mutated subclones in a typical MAPK-inhibitor resistant tumor. Sequencing of colonies derived from single cell clones of three of these double-mutated cultures confirmed that all except for one colony contained both BRAF and NRAS mutations, thus confirming the presence of both activation mutations in the same cells. Sequencing of the additional immunohistochemistry blocks from these patients only identified double mutations in the post-treatment samples, confirming the finding in patient 1 from cohort 1.

To understand the general resistance mechanisms of the double-mutated cells, we conducted viability assays with different MAPK inhibitors. The double-mutated cells grew in normal culturing conditions, notably without any MAPK inhibition, were all resistant to BRAF-inhibitors, but showed heterogeneity in their response to MEK or ERK inhibition, possibly because of co-existing mutations in other pathways. Combination treatment with MEK and BRAF inhibitors, as it is now clinical practice, showed synergism in MM121224 and M130903, but no synergistic or additional effect in the other cell cultures,Suggesting that simultaneous or second-line treatment with other MAPK-pathway inhibitors might still be effective in controlling progression in selected patients, but not in all. However, a MEK inhibitor combined with a PI3K, AKT or mTOR inhibitor was synergistic in all of the cell cultures, albeit with different strength. It has to be kept in mind however, as the double-mutated genotype was only present in one or two metastases from each patient, it is likely not the most important resistance mechanism in these patients and the efficacy of these second-line or combination treatments in controlling overall tumor burden is questionable. Since no common mechanism of resistance was found in any patient, it is possible that the other resistant tumors activated different pathways.

Except for 1 cell line (MM150423), all double mutated cell lines showed higher expression of pERK at the basal level compared to the single NRAS mutated cell line M010817. However, the basal level of pERK among the double mutated cells varied, with relatively high expression in MM140307 and MM140906 and relatively low expression in MM150423. It has been argued that BRAF and NRAS mutations are mutually exclusive due to a growth deficit of double mutated cells, because of senescence-inducing high levels of pERK [8]. In our experiment, MM121224 and MM140307 grew faster than the single mutated control cell lines and MM140906, MM150423 and M130903 grew slower, not supporting the view the cells with high pERK level grow slower. However, in vitro growth behavior might not represent in vivo growth, for instance depending on how well the cells have adapted to a 2D culture system, what growth factors are present or missing in the cell culture medium compared to the in vivo situation and how well the immune system can control metastasis formation.

The pERK levels in the cell lines under treatment of various MAPK pathway inhibitors showed some discrepancy with the proliferation assay. Most profound was the strong reduction in pERK upon treatment with the BRAF inhibitor and MEK inhibitor in MM140307, although the cell line was resistant to these two inhibitors. We hypothesize that this is due to the very high pERK levels in this cell line at baseline, in such a way that even a strong reduction compared to baseline does not suffice to block the pathway. This would also be true for MM140906 under ERK inhibitor treatment, although the cell line is resistant to ERK inhibition, levels of pERK show a decrease compaired to the baseline, but the levels are still high. MM150423 pERK levels are relatively low in MEK and ERK treated cells, although the cell line is partially resistant to these inhibitors. However, the pERK levels in the MEK and ERK treated cells do not differ considerably from the untreated cells.

In this study, we show that known-resistance mechanisms are present at low frequencies and heterogeneously within individual patients. Furthermore, we show that finding an additional NRAS mutation in a tumor sample following BRAF inhibitor treatment could indicate the presence of a double mutated subpopulation that is not necessarily sensitive to MEK inhibition or ERK inhibition, rendering the MEK inhibitor therapy in all such cases suboptimal. This study indicates that genetic analysis of one tumor biopsy does not fully define the resistance mechanism for the whole patient, which has important implications for secondary therapy strategies in case of primary resistance.

MATERIALS AND METHODS

Patients and sample preparation

Patients were selected after written consent from the patient, given through the university biobank program according to ethical approval numbers 647 and 800. We collected surplus material before and after therapy at autopsy. Samples were processed immediately after collection to ensure best possible DNA and RNA quality. Primary cell cultures were established as described in [18]. Notably, upon generation of cell cultures, all cultures were kept under standard conditions without additional treatment, even if they were derived from a therapy resistant patient.

DNA was isolated from paraffin embedded tissue, fresh frozen tissue, cultured cells and PMCs stored in the biobank of the institute of Dermatology of the University Hospital of Zürich. Germline DNA from PBMCs was sequenced for all patients if available as a reference [20].

DNA from paraffin blocks was isolated using the FFPE DNA isolation kit from Qiagen (QIAamp DNA FFPE Tissue Kit #56404) and optimized protocols developed by Ultan McDermott at the Sanger institute. Prior to DNA isolation, each block was evaluated by a trained dermato-histopathologist, and punches were made in tumor regions to ensure reduced contamination with stromal tissue.

For DNA isolation from non-paraffin embedded samples we followed standard DNA isolation protocols published earlier.

Library preparation and sequencing

DNA quality was measured by an Agilent 2100 Bioanalyzer or Agilent 2200 Tapestation. One to three ug of high quality DNA was used to prepare the whole exome library using the Agilent SureSelect V4 or V5 kit. Sequencing was performed on an Illumina Hiseq 2000 machine in the Functional Genomics Center at University of Zürich. For the whole exome sequencing we sequenced 0.25 lanes per sample, paired-end, with 100 bp reads.

Whole exome sequencing analysis

Bioinformatics analysis was conducted with a modified GATK pipeline [21–23]. Quality control was done with „FASTQC” [24]. Alignment of the FASTQ file to the reference genome “hg19” [25] Lander, Linton et al. 2001) and transformation from SAM to BAM was done with “BWA” [26]. PCR duplicates were marked by MarkDuplicates from “Picard” [22], Local realignment around indels with RealignerTargetCreator (GATK), realigning with IndelRealigner (GATK), fix mate information with FixMateInformation (Picard), base quality score recalibration with Baserecalibrator (GATK) and PrintReads (GATK). Variant calling was done with UnifiedGenotyper (GATK). For annotation of the VCF files we used Annovar [27], Samtools [28] and Bedtools [29]. For data interpretation we used Microsoft Access, Microsoft Excel, Venny [30], ConSet [31] and IGV [32, 33].

For copy number analysis we used Excavator [34] and Contra [35], results were visualized with Circos [36].

SNVs were filtered according to the following read count criteria: A base must have at least four mutant reads and at least 10 total reads, if less than 10 total reads, at least half of them must be mutated. Also all SNVs with a phred-scaled quality score of < 50 were excluded from further analysis. A SNV was called somatic if the unfiltered blood sample from the same patient did not show any mutant read for this position.

Mutant allele ratios (MAR) were calculated by dividing mutant read counts by total read counts for each called SNV. Frequencies for these ratios were calculated and trendlines were plotted in Excel with the Moving Average method (period: 3). To reduce the number of false positive SNVs we applied more strict filtering on the private SNVs. Quality threshold was raised to a phred score of 100, and the SNV needed to have at least 10 total reads. Genes that had more than 8 SNVs were excluded.

Deep sequencing of PCR amplicons containing NRAS exon 2

DNA of 7 tumor samples (EMG P5 cell culture, M121224, 401/II, 404/II, 403, H12.684, H12.12640/1/B) were amplified with primers containing a NRAS specific sequence (see chapter sanger sequencing), adaptor sequences and a unique multiplex-identifier (MID) sequence (according to eurofins protocol). Each tumor sample analyzed is carrying therefore the adaptor sequence and a unique MID sequence. The PCR product was gel purified and 200 ng of each amplicon was sending for deep sequencing. Amplicons were subjected to Roche 454 sequencing using emulsions-PCR. Data were analyzed using DeepSNV [37].

Sanger sequencing

After DNA amplification of NRAS and BRAF with the following primers: BRAF forward: 5′CTAAGAGGAAAGATG AAGTACTATG reverse: 5′CTAGTAACTCAGCAGCATCTCAG NRAS forward: 5′GATAGGCAGAAATGGGCTTGA reverse: 5′ATCATCCTTTCAGAGAAAATAATGC using a touchdown program going from 60 C to 55 C in 10 cycles, followed by 40 cycles at 55 C, the PCR product was diluted 100× and send to Microsynth for sequencing.

Generation of single cell clones and single cell clone sequencing

Cells were distributed over 96-well plate, containing 1 cell per well, via FACS cell sorting or serial dilutions. Cells were grown for several weeks under standard conditions [18] until visible colonies had formed. Then, medium was removed and wells were washed with PBS. Colonies were directly lysed in the well, with 10 ul lysis buffer (2.5% 1 M Tris pH 8.0 (Ambion), 0.1% 0.5 M EDTA (Sigma-Alderich), 0.25% Tween 20 (Sigma-Alderich), 1% proteinase K (Qiagen), Aqua dest.), and incubated at 55°C for 1 hour, followed by 5 min at 95°C. Afterwards, 10 ul 25 mM MgCl2 was added and the total volume was devided over 2 PCR reactions for NRAS and BRAF Sanger sequencing.

Cell sorting

Around 1 × 10⁷ melanoma cells were resuspended in 100 μl FACS buffer (1% FBS, 5 mM EDTA pH 8, 0.01% NaN₃/ddH₂O in PBS). Cells were incubated for 20 minutes at 4°C with the following photosensitive antibodies: Anti-human MCSP-FITC (1:20 dilution) (Miltenyi Biotec 130-098-794, Bergisch Gladbach Germany) and Anti-human Fibroblasts/Epithelial-PE (1:200 dilution) (ABIN319868, Aachen Germany). After washing, cells were resuspended in 200 μl FACS buffer and sorted using the Aria IIb (BD Biosciences, Franklin Lakes, New Jersey, USA).

Phylogenetic analysis

Maximum Parsimony, Bayesian and Maximum likelihood (ML) phylogenies were constructed with the POSIX-threads version of RAxML v8.0.19 (7). To correct for among-site rate heterogeneity using the Γ distribution, we used an ascertainment bias correction and a general time reversible (GTR) substitution model. Four rate categories (ASC_GTRGAMMA model) were used to calculate the optimal tree. Node support was evaluated with 100 nonparametric bootstrap pseudoreplicates, they therefore indicate the percentage of bootstrap trees that contained a given internode branch.

Variants diagnostic for a given clade are defined as existing solely in that clade and nowhere else for that position. All leaves emanating from the node in question must share a variant and all other leaves must contain a different character for a variant to be diagnostic. Diagnostic variants can therefore also be termed an apomorphy.

Cell viability assay

1 × 10⁴ cells were seeded and treated for 72 hours with different concentrations of either a BRAF inhibitor (PLX4032, LGX818 or GSK2118436), a MEK inhibitor (MEK162), an ERK inhibitor (SCH772984). DMSO treatment was used as a control. After 72 hours, the medium was removed and fresh RPMI1640 supplemented with 10% FCS and 8% MTT reagent (Sigma, 5 mg/ml in PBS) was added, and the cells were incubated at 37°C. After 1 hour, the RPMI1640 with MTT reagent was removed and 10% SDS (Sigma) and 95% isopropanol/ 5% Formic Acid (Sigma) (ratio 1:1) were added. After 5 min of incubation at 37°C, absorbance was measured at 595 nm (reference 620 nm) using a microplate reader.

Synergy calculations

Combination treatments and subsequent calculation of synergy were carried out according to the method from Chou and Talalay, with the compusyn software, available at http://www.combosyn.com/index.html. We have taken the mean of the raw CI values for the different concentration combinations, in order to determine it as overall synergistic (CI value < 0.9) or overall not synergistic (CI value > 0.9).

Proliferation assay

5 × 10⁴ cells/ml were seeded per T75 flask. After 24 hours, 72 hours, 144 hours and 240 hours the cells were counted. From the linear growth fase, the doubling time was calculated with http://www.doublingtime.com/compute.php.

Westernblot

Total protein was collected by washing cells twice with ice cold PBS and subsequent lysis in RIPA buffer (20 mM Tris-HCl (pH 7.5), 1% Triton X-100 (Sigma), 137 mM NaCl, 10% glycerol and protease inhibitors (Roche)). Concentration of the protein was measured with the Bio-Rad Dc Protein Assay (Bio-Rad, Hercules, CA, USA) according to the manufacturer’s protocol. SDS-Page was used to separate the proteins, after which they were transferred onto a nitrocellulose membrane. Membranes were probed with a rabbit anti-pERK antibody (Cell Signaling, product nr #4376S), a rabbit anti-ERK antibody (Cell Signalling, product nr#9102) and a rabbit anti-GAPDH antibody (Abcam, Cambridge, UK, product nr ab9385), followed by horseradish peroxidase-conjugated goat anti-rabbit IgG (Santa Cruz, product nr sc-2030)Bound antibodies were detected using chemiluminescence (ECL, GE Healthcare, Chalfont St. Giles, UK). Afterwards, band intensity was measured using ImageJ software (imagej.nih.gov/ij/) and pERK band intensity was corrected for corresponding GAPDH band intensity.

CONFLICTS OF INTEREST

The authors have no conflicts of interest to disclose.

GRANT SUPPORT

DW was partially funded by a fellowship from the Roche Hub Project (F-85807-09-01) (http://www.roche.com/index.htm), ML was partially funded by the Society for Skin Cancer Research (http://www.skincancer.ch/), MIGR was partially funded by the Swiss Cancer League (KLS-3151-02-2013), and the biobank was partially funded by the University of Zurich Research Priority Program (URPP) in translational cancer research (U-402-02-01) (http://www.cancer.uzh.ch/aboutus.html). The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.

REFERENCES

1. Dhillon AS, Hagan S, Rath O, Kolch W. MAP kinase signalling pathways in cancer. Oncogene. 2007; 2:3279–3290.

2. Hodis E, Watson IR, Kryukov GV, Arold ST, Imielinski M, Theurillat JP, Nickerson E, Auclair D, Li L, Place C, Dicara D, Ramos AH, Lawrence MS, et al. A landscape of driver mutations in melanoma. Cell. 2012; 15:251–263.

3. McArthur GA, Chapman PB, Robert C, Larkin J, Haanen JB, Dummer R, Ribas A, Hogg D, Hamid O, Ascierto PA, Garbe C, Testori A, Maio M, et al. Safety and efficacy of vemurafenib in BRAF(V600E) and BRAF(V600K) mutation-positive melanoma (BRIM-3): extended follow-up of a phase 3, randomised, open-label study. Lancet Oncol. 2014; 1:323–332.

4. Chapman PB, Hauschild A, Robert C, Haanen JB, Ascierto P, Larkin J, Dummer R, Garbe C, Testori A, Maio M, Hogg D, Lorigan P, Lebbe C, et al. Improved survival with vemurafenib in melanoma with BRAF V600E mutation. N Engl J Med. 2011; 36:2507–2516.

5. Shi H, Hugo W, Kong X, Hong A, Koya RC, Moriceau G, Chodon T, Guo R, Johnson DB, Dahlman KB, Kelley MC, Kefford RF, Chmielowski B, et al. Acquired resistance and clonal evolution in melanoma during BRAF inhibitor therapy. Cancer Discov. 2014; 80–93.

6. Van Allen EM, Wagle N, Sucker A, Treacy DJ, Johannessen CM, Goetz EM, Place CS, Taylor-Weiner A, Whittaker S, Kryukov GV, Hodis E, Rosenberg M, McKenna A, et al. The genetic landscape of clinical resistance to RAF inhibition in metastatic melanoma. Cancer Discov. 2014; 4:94–109.

7. Sensi M, Nicolini G, Petti C, Bersani I, Lozupone F, Molla A, Vegetti C, Nonaka D, Mortarini R, Parmiani G, Fais S, Anichini A. Mutually exclusive NRASQ61R and BRAFV600E mutations at the single-cell level in the same human melanoma. Oncogene. 2006; 2:3357–3364.

8. Petti C, Molla A, Vegetti C, Ferrone S, Anichini A, Sensi M. Coexpression of NRASQ61R and BRAFV600E in human melanoma cells activates senescence and increases susceptibility to cell-mediated cytotoxicity. Cancer research. 2006; 6:6503–6511.

9. Chiappetta C, Proietti I, Soccodato V, Puggioni C, Zaralli R, Pacini L, Porta N, Skroza N, Petrozza V, Potenza C, Della Rocca C, Di Cristofano C. BRAF and NRAS mutations are heterogeneous and not mutually exclusive in nodular melanoma. Appl Immunohistochem Mol Morphol. 2015; 2:172–177.

10. Romano E, Pradervand S, Paillusson A, Weber J, Harshman K, Muehlethaler K, Speiser D, Peters S, Rimoldi D, Michielin O. Identification of multiple mechanisms of resistance to vemurafenib in a patient with BRAFV600E-mutated cutaneous melanoma successfully rechallenged after progression. Clin Cancer Res. 2013; 1:5749–5757.

11. Nazarian R, Shi H, Wang Q, Kong X, Koya RC, Lee H, Chen Z, Lee MK, Attar N, Sazegar H, Chodon T, Nelson SF, McArthur G, et al. Melanomas acquire resistance to B-RAF(V600E) inhibition by RTK or N-RAS upregulation. Nature. 2010; 46:973–977.

12. Wong DJ, Robert L, Atefi MS, Lassen A, Avarappatt G, Cerniglia M, Avramis E, Tsoi J, Foulad D, Graeber TG, Comin-Anduix B, Samatar A, Lo RS, et al. Antitumor activity of the ERK inhibitor SCH772984 [corrected] against BRAF mutant, NRAS mutant and wild-type melanoma. Mol Cancer. 2014; 13:194.

13. Greger JG, Eastman SD, Zhang V, Bleam MR, Hughes AM, Smitheman KN, Dickerson SH, Laquerre SG, Liu L, Gilmer TM. Combinations of BRAF, MEK, and PI3K/mTOR inhibitors overcome acquired resistance to the BRAF inhibitor GSK2118436 dabrafenib, mediated by NRAS or MEK mutations. Mol Cancer Ther. 2012; 1:909–920.

14. Tirosh I, Izar B, Prakadan SM, Wadsworth MH, 2nd, Treacy D, Trombetta JJ, Rotem A, Rodman C, Lian C, Murphy G, Fallahi-Sichani M, Dutton-Regester K, et al. Dissecting the multicellular ecosystem of metastatic melanoma by single-cell RNA-seq. Science. 2016; 35:189–196.

15. Widmer DS, Eichhoff OM, Dummer R, Levesque MP. Melanoma’s next top model, it is in the air. Exp Dermatol. 2015; 2:659–660.

16. Sanborn JZ, Chung J, Purdom E, Wang NJ, Kakavand H, Wilmott JS, Butler T, Thompson JF, Mann GJ, Haydu LE, Saw RP, Busam KJ, Lo RS, et al. Phylogenetic analyses of melanoma reveal complex patterns of metastatic dissemination. PNAS. 2015; 11:10995–11000.

17. Takata M, Morita R, Takehara K. Clonal heterogeneity in sporadic melanomas as revealed by loss-of-heterozygosity analysis. Int J Cancer. 2000; 8:492–497.

18. Raaijmakers MI, Widmer DS, Maudrich M, Koch T, Langer A, Flace A, Schnyder C, Dummer R, Levesque MP. A new live-cell biobank workflow efficiently recovers heterogeneous melanoma cells from native biopsies. Exp Dermatol. 2015; 2:377–380.

19. Heidorn SJ, Milagre C, Whittaker S, Nourry A, Niculescu-Duvas I, Dhomen N, Hussain J, Reis-Filho JS, Springer CJ, Pritchard C, Marais R. Kinase-dead BRAF and oncogenic RAS cooperate to drive tumor progression through CRAF. Cell. 2010; 14:209–221.

20. Böyum A. Isolation of mononuclear cells and granulocytes from human blood. Isolation of monuclear cells by one centrifugation, and of granulocytes by combining centrifugation and sedimentation at 1 g. Scand J Clin Lab Invest Suppl. 1968; 97:77–89.

21. McKenna A, Hanna M, Banks E, Sivachenko A, Cibulskis K, Kernytsky A, Garimella K, Altshuler D, Gabriel S, Daly M, DePristo MA. The Genome Analysis Toolkit: a MapReduce framework for analyzing next-generation DNA sequencing data. Genome Res. 2010; 2:1297–1303.

22. Van der Auwera GA. From FastQ Data to High-Confidence Variant Calls: The Genome Analysis Toolkit Best Practices Pipeline. In: G. A. C, ed. 2013.

23. DePristo MA, Banks E, Poplin R, Garimella KV, Maguire JR, Hartl C, Philippakis AA, del Angel G, Rivas MA, Hanna M, McKenna A, Fennell TJ, Kernytsky AM, et al. A framework for variation discovery and genotyping using next-generation DNA sequencing data. Nat Genet. 2011; 4:491–498.

24. Andrews S. FastQC A Quality Control tool for High Throughput Sequence Data.

25. Lander ES, Linton LM, Birren B, Nusbaum C, Zody MC, Baldwin J, Devon K, Dewar K, Doyle M, FitzHugh W, Funke R, Gage D, Harris K, et al. Initial sequencing and analysis of the human genome. Nature. 2001; 40:860–921.

26. Li H, Durbin R. Fast and accurate short read alignment with Burrows-Wheeler transform. Bioinformatics. 2009; 2:1754–1760.

27. Wang K, Li M, Hakonarson H. ANNOVAR: functional annotation of genetic variants from high-throughput sequencing data. Nucleic Acids Res. 2010; 3:e164.

28. Li H, Handsaker B, Wysoker A, Fennell T, Ruan J, Homer N, Marth G, Abecasis G, Durbin R, Subgroup GPDP. The Sequence Alignment/Map format and SAMtools. Bioinformatics. 2009; 2:2078–2079.

29. Quinlan AR, Hall IM. BEDTools: a flexible suite of utilities for comparing genomic features. Bioinformatics. 2010; 2:841–842.

30. Oliveros JC. VENNY. An interactive tool for comparing lists with Venn diagrams. 2007.

31. Kim B, Lee B, Seo J. Visualizing Set Concordance with Permutation Matrices and Fan Diagrams. Interact Comput. 2007; 1:630–643.

32. Thorvaldsdóttir H, Robinson JT, Mesirov JP. Integrative Genomics Viewer (IGV): high-performance genomics data visualization and exploration. Brief Bioinform. 2013; 1:178–192.

33. Robinson JT, Thorvaldsdóttir H, Winckler W, Guttman M, Lander ES, Getz G, Mesirov JP. Integrative genomics viewer. Nat Biotechnol. 2011; 2:24–26.

34. Magi A, Tattini L, Cifola I, D’Aurizio R, Benelli M, Mangano E, Battaglia C, Bonora E, Kurg A, Seri M, Magini P, Giusti B, Romeo G, et al. EXCAVATOR: detecting copy number variants from whole-exome sequencing data. Genome Biol. 2013; 1:R120.

35. Li J, Lupat R, Amarasinghe KC, Thompson ER, Doyle MA, Ryland GL, Tothill RW, Halgamuge SK, Campbell IG, Gorringe KL. CONTRA: copy number analysis for targeted resequencing. Bioinformatics. 2012; 2:1307–1313.

36. Krzywinski M, Schein J, Birol I, Connors J, Gascoyne R, Horsman D, Jones SJ, Marra MA. Circos: an information aesthetic for comparative genomics. Genome Res. 2009; 1:1639–1645.

37. Gerstung M, Beisel C, Rechsteiner M, Wild P, Schraml P, Moch H, Beerenwinkel N. Reliable detection of subclonal single-nucleotide variants in tumour cell populations. Nat Commun. 2012; 3:811.