Co-clinical assessment identifies patterns of BRAF inhibitor resistance in melanoma

BRAF inhibitors have contributed to a significant improvement in survival rates for melanoma patients whose tumors have a hotspot V600E/K–activating mutation in the BRAF oncogene (1, 2). In addition to a majority of patients experiencing tumor regression and prolonged survival, many reports have documented major improvements in quality of life, including improved physical activity and emotional state (3–5). Unfortunately, it is also now well documented that BRAF inhibitors, and even the superior combination of BRAF and MEK inhibitors, produce primarily short-term responses that typically last less than 1 year, followed by the emergence of resistance (6). Therefore, an improved understanding of the genetic and epigenetic mechanisms that confer resistance is required to prolong the benefits of BRAF inhibition.

Recent whole-exome and RNA sequencing studies have identified a wide array of acquired mutations that confer resistance, including those that reactivate the MAPK pathway (NRAS, KRAS, and MEK1/2 mutations, NF1 loss, BRAF amplification, and BRAF splice variants) (7–9) and those that activate the PI3K pathway (PIK3CA, PIK3R1, and AKT1/2 mutations and PTEN loss) (10–12). Each of these provides insight into candidate second-line therapies that could potentially bypass the resistance mechanism; these include, for example, pan-RAF (13) and ERK inhibitors (14, 15) or PI3K/AKT/mTOR inhibitors (16–19). However, up to one-quarter to two-fifths (11, 12) of patients’ tumors do not harbor any of the known resistance-conferring mutations, making it challenging to identify genomics-based second-line therapies for these patients.

To address this gap in knowledge, we have undertaken a cross-species analysis of BRAF inhibitor–resistant human and mouse melanomas, the latter derived from a genetically engineered BRAF–driven mouse melanoma model. Our hypothesis is that cross-species comparative analysis of resistance based on a combination of protein-signaling patterns and resistance-conferring mutations could provide clinically actionable information and assist in the stratification of patients into defined resistance classes for downstream therapeutic decisions.

A novel mouse model of BRAF inhibitor resistance. To model BRAF inhibitor resistance, we generated a doxycycline- and tamoxifen-inducible mouse model of BRAF^V600E melanoma. Briefly, the mouse has a Tet-inducible human BRAF^V600E transgene (20), a constitutive Cdkn2a-knockout allele (21), a conditional knockout allele of Pten (22), and inducible Cre expression under melanocyte-specific control (23). Upon the topical application of tamoxifen, Pten was specifically deleted only in the treated melanocytes, and rtTA was activated. Subsequent administration of doxycycline in the diet activated the BRAF^V600E transgene only in the cells in which both the LSL-Stop-rtTA cassette and Pten were codeleted (Supplemental Figure 1; supplemental material available online with this article; doi:10.1172/JCI78954DS1). After topical administration of as little as 1 μl of 10 μM 4-hydroxy-tamoxifen, tumors that were BRAF^V600E positive and CDKN2A- and PTEN null developed with a tightly distributed latency (median = 60 days) and high penetrance (85%) (Supplemental Figure 1).

We first demonstrated that after melanoma formation in these “iBIP” (inducible BRAF INK/ARF PTEN (iBIP) mice, withdrawal of doxycycline resulted in extinction of BRAF^V600E transgene expression, leading to rapid tumor regression (Figure 1, A and B, and Supplemental Figure 1), similar to that seen in an inducible NRAS^Q61K melanoma model (24). Next, administration of 417 parts per million (ppm) of the PLX4720 BRAF inhibitor in the chow, with mice remaining on doxycycline to ensure BRAF^V600E transgene expression, reproducibly led to potent tumor growth inhibition. This manifested as a greater than 30% tumor regression by total volume in 56% (9 of 16) of treated mice as the best response (Figure 1C). After continual administration of PLX4720, we observed the emergence of drug resistance in these autochthonous iBIP tumors at a median of 32 days (Figure 1D).

Figure 1

Generation of BRAF inhibitor–resistant iBIP mouse melanomas. (A) Real-time PCR of the human BRAF^V600E transgene on and off doxycycline in autochthonous iBIP melanomas. (B) Representative tumor growth curves of 11 mice on doxycycline, off doxycycline, or on PLX4720 chow. Differing shades of orange are used only for clarity, and all refer to PLX4720 treatment. (C) Waterfall plot of best response of tumors to PLX4720, shown as the percentage of tumor volume change from baseline of all tumors added together. (D) Kaplan-Meier survival plot of the time to PLX4720 resistance. The mice in C and D were from separate cohorts. dox, doxycycline; Rx, treatment.

Since the BRAF extinction phenotype provides a positive control with which to compare pharmacological BRAF inhibition, we asked to what extent PLX4720 reproduced the effects of genetic extinction of BRAF. First, we determined that in both iBIP nude and syngeneic allografts, PLX4720 and BRAF extinction faithfully produced tumor regressions comparable to those of autochthonous iBIP tumors (Figure 2A). Next, using reverse-phase protein array (RPPA) and expression microarray analysis of allograft samples in nude mice, we noted that pharmacological and genetic BRAF inhibition were globally similar, though not identical (Figure 2, B and C). Specifically, in the RPPA dataset, both pharmacological and genetic BRAF inhibition induced downregulation of phosphorylated ERK (p-ERK), CCND1, p-RSK, p-S6, and p-Rb and upregulation of p-AKT, p-GSK3b, and BIM (Figure 2B), consistent with known perturbed pathways. Next, we examined the molecular phenotypes by IHC analysis of time-matched samples and determined that PLX4720 robustly induced both cell cycle arrest and apoptosis, though to a lesser extent than did genetic extinction (Figure 2, D and E). Given the molecular and phenotypic similarities of PLX4720 treatment to genetic BRAF extinction, albeit at the expected lower potency, we conclude that PLX4720-mediated BRAF inhibition in this model is on target and molecularly relevant. Accordingly, we isolated a number of resistant autochthonous mouse melanomas from iBIP mice for downstream analyses, as presented below.

Figure 2

Comparison of PLX4720 and genetic BRAF extinction. (A) Comparison of PLX4720 treatment and BRAF extinction in iBIP allografts. (B) RPPA data comparing BRAF extinction and PLX4720 daily time course treatments. Two tumors per time point. Gray boxes are log₂ = 0. (C) Log₂ microarray data comparing the fold change versus on-doxycycline controls of: 3 days of continual PLX4720 treatment versus 3 days of BRAF extinction via doxycycline withdrawal. Data are averaged from n = 2 each. (D and E) Proliferation and apoptosis in treated nude allograft tumors, as measured by phospho-histone 3 (p-H3) and TUNEL IHC staining, respectively. **P < 0.005 by Student’s t test. PLX, PLX4720.

Comparison of BRAF inhibitor–resistant human and mouse melanomas. Our goal was to collect and characterize a series of human and mouse BRAF inhibitor–resistant melanoma samples for cross-species comparative analysis. From 2009 to 2012, we collected 14 resistant tumor biopsies from 13 human patients, 5 on BRAFi (vemurafenib) and 9 on BRAFi plus MEKi (dabrafenib plus trametinib) treatment, under appropriate IRB approval. We also harvested 10 independent, resistant autochthonous tumors from 5 iBIP mice. In addition, we collected 13 pretreatment and 14 on-treatment (range: 12 ± 5 days) melanomas from patients, including both matched and nonmatched samples (see Supplemental Table 1 for full patient information). Similarly, we isolated 4 unmatched, but isogenic, pretreatment mouse tumors. Hereafter, for convenience, we refer to these 3 defined time points as A, B, and C for each sample, namely: A samples correspond to pretreatment specimens, B to on-treatment specimens, and C to clinical progression (i.e., resistant) specimens.

For each specimen, we isolated protein and submitted all lysates for analysis by RPPA at the same time. At the time of the analysis, the RPPA platform used 200 antibodies targeting 150 total proteins and 50 phospho-specific antibodies. We first sought to confirm the expected action of the BRAF inhibitor by comparing on-treatment B samples to matched pretreatment A samples. Markers of proliferation, specifically p-Rb (2.1 × 10^–8) and cyclin B1 (1.0 × 10^–4), were downregulated across all B samples compared with those in A samples and were among the most significantly affected proteins (Supplemental Figure 2). The proapoptotic protein BIM was upregulated in 6 of 10 samples, and the apoptosis marker cleaved caspase 7 was upregulated in 5 of 10 samples (Supplemental Figure 2). These results demonstrate that the RPPA platform is able to detect therapy-relevant protein changes.

We next asked whether the drug-resistant human C samples could be stratified into distinct classes. As shown in Figure 3A, unsupervised hierarchical clustering of the RPPA data readily separated the C samples into 3 distinct clusters, with good intracluster pattern similarity (Figure 3A). Notably, the assignment to these 3 categories of resistance patterns did not correlate with known clinical parameters or treatment regimens (Table 1), nor did the pretreatment patterns appear to predict the category assignment (n = 5 matched pairs). Importantly, when the 10 resistant mouse C samples were included in the cluster analysis, the human samples retained their original class assignment, indicating that their initial clustering was stable (Figure 3B). Interestingly, the mouse C samples clustered with human samples in clusters 1 and 3, suggesting that the mouse tumors showed 2 distinct, clinically relevant patterns of resistance, despite their isogenic background.

Figure 3

RPPA data on human and mouse melanoma samples. (A) Human resistant C samples clustered using unsupervised hierarchical clustering. (B) Addition of 9 mouse samples to the human clustering. (C) Overview of human A, B, and C samples and mouse C samples showing only selected antibodies that mechanistically describe the C clusters. The color bar applies to all panels.

Table 1

Summary of human samples analyzed by RPPA, RNA-seq, and/or whole-exome sequencing

We next characterized each cluster of resistant tumors on the basis of the combined human and mouse data. Cluster 1 consisted of 3 human and 6 mouse samples that showed a high degree of similarity, all displaying an increase in multiple activated RTKs including p-EGFR, p-ERBB2, p-ERRB3, and p-MET (Figure 3C), relative to other samples. This led us to call cluster 1 the “RTK” cluster. This cluster showed markedly lower levels of MAPK, PI3K, and mTOR phosphoproteins (e.g., p-ERK, p-AKT, p-4EBP1, and p-S6K) compared with levels in the other clusters, suggesting MAPK-independent resistant pathway(s) (Figure 3C). Cluster 2 consisted of 7 human samples and uniformly showed a marked increase in p-ERK, p-MEK, p-Rb, and cyclin B1 compared with that detected in the other samples (Figure 3C). Notably, p-ERK and p-MEK levels in these samples were also higher than in the pretreatment samples, including matched samples from patients 16 and 22 (Figure 4A). We therefore call cluster 2 the “hyper-MAPK” cluster. To validate these 2 cluster assignments, we performed IHC, which showed that the RTK tumors had substantially lower p-ERK signals than did the hyper-MAPK tumors (Supplemental Figure 3). Cluster 3 had greater heterogeneity in its 4 human samples, but appeared to resemble the pretreatment group. Here, the 4 mouse samples in this cluster proved informative, as they clustered tightly with the isogenic pretreatment mouse samples (Supplemental Figure 2). We therefore called cluster 3 the “MAPK-rebound” cluster, as it indicated a possible return to pretreatment MAPK signaling levels. Taken together, these results indicate that the pattern of clustering likely reflects molecular similarities among the samples and that the mouse model recapitulates clinically relevant patterns of BRAF inhibitor resistance.

Figure 4

Comparison of RPPA and RNA-seq data. (A) Alignment of ERK-related probes. (B) Alignment of proliferation-related probes. (C) Alignment of RTK-related probes. (D) BRAF splice forms identified by RNA-seq.

Comparison of RNA-sequencing and RPPA human data. We also performed RNA sequencing (RNA-seq) of 38 patient samples, 31 of which overlapped with the RPPA set (Table 1). This allowed for direct comparison and correlation of RNA-seq and RPPA data. First, we determined whether an ERK signature in the RNA data matched the RPPA profile. As shown in Figure 4A, there was a general overall correlation across all samples (R² = 0.4), including the expected decrease in on-treatment B samples and enrichment in the hyper-MAPK cluster of resistant tumors (Figure 4A). We also demonstrated a correlation among p-Rb, CCNB1, and FOXM1 antibody levels with a cell cycle gene expression signature (R² = 0.4), again consistent with an antiproliferative effect of BRAF inhibition (Figure 4B). Since the data demonstrated a good correlation between RPPA and RNA-seq, we next assessed the RNA levels of the RTKs. Interestingly, not every activated RTK showed concomitant mRNA upregulation at the RNA level (Figure 4C), suggesting that at least some of the increased RTK activation reflected in the RPPA data is post-translational, consistent with the known phosphorylation-dependent activity of RTKs.

We then showed the veracity of our RNA-seq data by confirming the presence of the known BRAF^V600E/K mutation in each sample by visual inspection of the Binary Alignment/Map (BAM) files (Supplemental Table 2). Next, we used these RNA-seq data to attempt identification of putative resistance-causing mutations by analyzing the coding sequence of 16 genes reported to confer resistance to BRAF inhibitor treatment (Supplemental Table 1). Since matched pretreatment A and on-treatment B pairs were not available in many cases, we used the hg19 reference genome as the comparator in those progressed samples to identify putative resistance-conferring mutations. Consistent with previous reports (11, 12), 50% (6 of 12) of the samples appeared to be WT at these 16 previously reported resistance-conferring genes at the sequence and splice levels (Supplemental Table 1). We found no NRAS^Q61 hotspot mutations or activating MEK1/2 mutations in this cohort of samples. Of the 6 samples with putative resistance-conferring alterations, 15C harbored an acquired missense PTEN^R159S mutation in the phosphatase domain, 25C harbored a known acquired MEK^Q60L mutation (7), 16C had BRAF overexpression consistent with a known acquired amplification (7), and 22C, 28C, and 29C showed resistance-conferring BRAF splice forms (Figure 4D and Table 1). To further verify these findings, we leveraged the whole-exome sequencing data available for 6 of these 12 samples (15C, 24C, and 27C performed in this study and 16C, 22C, and 25C performed in ref. 7) to confirm the presence of putative resistance-conferring RNA-seq mutation calls, i.e., PTEN and MEK2 in 15C and 25C, respectively (Table 1 and Supplemental Table 1), as well as the WT status at the 16 known resistance-conferring genes for samples 24C and 27C (Supplemental Table 2).

As the PTEN^R159S mutation, which has been previously described (25), occurs in the RTK class, we asked whether PTEN inactivation correlates with RTK activation in the TCGA melanoma RPPA dataset. We found no evidence of a correlation (Supplemental Figure 2), suggesting that PTEN loss is not generally associated with pan-RTK activation in melanoma. Since amplification of the BRAF locus without secondary coding changes has been reported as a resistance mechanism (8, 26), we sought to infer this from the expression data for the 6 samples without whole-exome data. The hyper-MAPK samples 17C and 20C had highly elevated levels of BRAF expression and of its neighboring genes, but quantitative PCR (qPCR) analysis did not indicate amplification of the BRAF locus (Supplemental Figure 4). To corroborate this finding, we explored whether high-level BRAF expression occurs in the absence of genomic amplification in other datasets, specifically the TCGA melanoma data. Our analysis of the TCGA melanoma data revealed that, while high-level (>2-fold over the median) BRAF expression was common in samples with high-level focal amplification (4 of 6, 66%) as expected, there were comparably high levels of expression in some nonamplified samples (20 of 346, 6%) as well.

Overall, for human samples with observable or inferable resistance-causing changes, the sequencing data offer molecular support for the cluster assignment by RPPA: the hyper-MAPK class correlated with BRAF gene or mRNA changes in 5 of 5 assayed samples including BRAF overexpression; the RTK class showed 1 missense PTEN mutation in 1 of 3 samples; and the MAPK-rebound class showed an MEK2 mutation and an alternative BRAF splice form in 2 of 4 samples, with no obvious candidates for the remaining 2 samples (Table 1). Collectively, these observations suggest that the RPPA profile can inform the category of the molecular mechanism driving resistance, despite the absence of known resistance-conferring genetic mutations.

Given the prominence of MAPK pathway reactivation in resistant samples, we asked whether ERK activity is sufficient to classify the samples. Accordingly, we clustered the samples by RPPA-derived p-ERK/p-MEK, the RNA-seq ERK signature, or the combined RNA-seq ERK and cell cycle signatures. None of these classified the resistant samples in the same way as did the full RPPA dataset (Supplemental Figure 5). Indeed, both the ERK signature and p-MEK were similar between the RTK and MAPK-rebound groups; instead, Figure 3C indicates that the main MAPK pathway difference between the 2 groups was downstream outputs including p-4EBP1, p-EIF4E, p-S6K, and p-S6, which have recently been shown to be prominent markers and mediators of BRAFi resistance (27, 28). Furthermore, unsupervised clustering of the most variable 2.5% of the RNA-seq genes also did not generate the same 3 clusters (Supplemental Figure 5). Together, these analyses suggest that the multiple, mechanistically relevant pathways queried by RPPA are necessary for robust classification: for the supervised analyses, the lack of mTOR or RTK activation signatures hampers stratification, while for the unsupervised analyses, the relevant signatures are likely statistically drowned out by unrelated gene sets, including immune ones. These data argue that more than just immediate ERK effectors are driving the classification of resistant tumors, though we cannot rule out the possibility that robust mRNA signatures are derivable from larger cohorts.

Patterns of resistance induced by oncogenic RTKs. We next assessed whether the oncogenic RTKs activated in the RTK cluster are capable of inducing resistance in an in vitro system. Physiologically relevant levels of overexpression of EGFR, ERBB2, or MET, but not ERBB3 (Supplemental Figure 6), conferred resistance to the BRAF inhibitor dabrafenib in the BRAF^V600E-mutant A375 melanoma cell line (Figure 5A and Supplemental Figure 6). The degree of resistance conferred by each positive-scoring RTK was comparable to that of the positive control RAF1, which encodes for CRAF. These targeted functional validation results are consistent with broader ORF resistance screens performed by other laboratories using the same cell line (29, 30).

Figure 5

Several RTKs confer resistance to BRAF inhibition. (A) Effects of candidate gene overexpression on dabrafenib sensitivity in the BRAF^V600E melanoma cell line A375. Results represent crystal violet measurements normalized to the GFP control. (B) RPPA analysis of candidate gene overexpression. The rightmost columns represent treated lines normalized to treated GFP. All lines were assayed in technical duplicates. (C) Xenografts from the A375 pHAGE lines, untreated or treated with PLX4720 chow. (D) Staining of A375 pHAGE tumors: red, p-ERK; blue, DAPI.

To assess whether overexpression of these RTKs induced a protein activation pattern similar to that of the RTK cluster, we profiled each overexpressing cell line by RPPA, with and without 100 nM dabrafenib. We found that ERBB2, MET, and RAF1, but not EGFR or ERBB3, mainly counteracted MAPK and mTOR signaling inhibition in the presence of dabrafenib (Figure 5B). Indeed, ERBB2 and MET overexpression largely phenocopied RAF1 overexpression, which primarily activated the MAPK pathway including p-ERK and p-MEK, but also the downstream mTOR effectors p-S6, p-S6K, and p-4EBP1 (Figure 5B). We then noted that the overall RPPA resistance pattern for ERBB2 and MET differed significantly from the RTK cluster pattern in the patient samples, wherein MAPK/mTOR signaling was relatively low (Figure 3C). By contrast, EGFR overexpression induced a degree of resistance similar to that of ERBB2, MET, and CRAF, but without reactivation of MAPK and mTOR (Figure 5, A and B). We confirmed these results in vivo, as A375 xenografts overexpressing EGFR were resistant to PLX4720 (Figure 5C). Although EGFR modulated baseline p-ERK, it did not protect p-ERK upon PLX4720 treatment (Figure 5D). These general EGFR observations were reproduced in vitro in an additional BRAF-mutant cell line, WM88, in which EGFR overexpression induced resistance without high-level p-ERK induction (Supplemental Figure 6). These results are consistent with those of a recent study showing that human samples overexpressing EGFR and other RTKs possess a slow-growing resistance phenotype (31), similar to that of the RTK cluster’s relatively low levels of proliferative markers (Figure 4B). These data suggest that activated EGFR contributes to the RTK cluster resistance phenotype.

Identification of molecular signatures predictive of treatment response. To identify molecular signatures predictive of clinical response, we used the RNA-seq data to rank the expression of all genes according to their Pearson correlations with either progression-free survival (PFS) or degree of tumor regression (Response Evaluation Criteria in Solid Tumors [RECIST]). The ranked gene set enrichment analysis (GSEA) allowed for an unbiased assessment of pathways that correlated with the clinical outcomes. We observed an anticorrelation between proliferation and PFS for both pre- and on-treatment samples (Figure 6A, Supplemental Figure 7, and Supplemental Table 3), which could be visualized with a 22-gene signature (Figure 6, B and C, and Supplemental Table 5). This is supported by the RPPA data, in which the proliferation marker cyclin B1 was one of the proteins with the highest correlation with PFS (Figure 6D), in keeping with similar findings in breast cancer (32, 33). In other words, the more highly proliferative the pretreatment tumor, the faster resistance emerged, perhaps due in part to increased genomic heterogeneity acquired through additional rounds of cell division (Figure 6C). Furthermore, a similar correlation with overall survival was seen when we reanalyzed the results of a study using Ki67 IHC to measure proliferation (Figure 6E) (34). We note that in the present study, it is the basal levels of proliferation that correlated with PFS, not the relative decrease in the proliferative signature induced by BRAF inhibition (Supplemental Tables 3 and 4).

Figure 6

Clinical correlations. (A and F) GSEA FWER q values for top-scoring pathways versus PFS or percentage of tumor regression. (B and C) Correlation of a 22-gene cell cycle signature with PFS. (D) CCNB1 RPPA data. (E) Reanalysis of data from ref. 34. (G and H) Correlations of a 97-gene immune signature with tumor regression. (I) LCK RPPA data.

In contrast, RECIST measurement of tumor regression positively correlated with immune activation signatures in both the pre- and on-treatment groups (Figure 6, F–H, Supplemental Figure 7, and Supplemental Table 3). These included genes involved in T cell activation (e.g., CD28 and CD86) and recruitment (e.g., various interleukins) as well as markers of T cells themselves (e.g. CD3E and CD4) (Supplemental Table 5). Again supporting these RNA expression findings are the RPPA data showing on-treatment LCK antibody levels (an immune marker) to be among the most highly correlated with RECIST responses (Figure 6I). Taken together, these observations suggest that tumors with primed immune systems might be better poised for response to treatment, consistent with a recent report on prognostic immune serum markers in a BRAF inhibitor trial (35) and with the observation that microarray-based immune infiltration signatures positively correlate with survival in non–BRAFi-treated melanoma patient populations (36–39). Additionally, the large-scale TCGA (The Cancer Genome Atlas) study has confirmed that elevated lymphocytic RNA-seq signatures and LCK RPPA signals positively correlate with outcome in over 300 patients (https://tcga-data.nci.nih.gov/docs/publications/skcm_2014/). Other high-scoring pathways include ribosome and chromatin-modifying genes that might serve as additional markers for PFS and RECIST, respectively (Supplemental Tables 3 and 4). Surprisingly, we found that levels of MAPK activation or apoptotic induction signatures did not significantly correlate with PFS or RECIST (P > 0.5) on either RPPA or RNA-seq levels. Finally, it is notable that the RNA signatures and protein markers were defined at the pretreatment and 2-week on-treatment time points, both preceding the standard RECIST measurement of best overall response. Thus, we hypothesize that these candidate genes could potentially serve as predictive response biomarkers to offer guidance on therapeutic benefit up front. Larger cohorts will be required for validation.

In this study, we have demonstrated that targeted proteome profiling can provide a rapid and cost-effective view of BRAF inhibitor resistance patterns in human melanoma samples. Resistant human samples displayed at least 3 distinct RPPA patterns, 2 of which were recapitulated in our molecularly similar genetically engineered mouse (GEM) model, offering an opportunity for functional studies and further target validation. Each RPPA cluster is representative of underlying mutational and gene expression changes and can be used to stratify resistant samples into MAPK-dependent and MAPK-independent mechanisms that may be of clinical value in the selection of or exclusion from potential second-line therapies. These resistance classifications are not dependent on knowledge of the underlying resistance-conferring DNA alterations; thus, they may provide complementary information on patients for whom targeted or genome-wide sequencing is uninformative or unavailable, regardless of whether mutations are undetected for biological or technical reasons. In addition, for the 2 key clinical measurements of PFS and tumor regression, we used combined RNA-seq and RPPA data to identify potential predictive biomarkers of patient response to BRAF inhibitors. Specifically, our findings that an early immune signature in pretreatment samples predicts tumor regression and that the proliferative index inversely correlates with PFS offer the intriguing possibility of using pretreatment biopsies in decisions such as patient selection for targeted agents and the timing and/or use of immunotherapies.

Together with several recent publications describing BRAF inhibitor resistance mechanisms identified by sequencing (7, 8, 11, 12), our data broaden the knowledge base and provide an initial look at the protein landscape. Consistent with the current study, these other whole-exome or targeted sequencing reports similarly identified a portion of samples, ranging from 29% to 42%, without obvious detectable resistance-conferring mutations (11, 12). Furthermore, Rizos and colleagues (12) found that 21% of samples had no MAPK reactivation, as assessed by the presence of an RNA signature, consistent with the non–MAPK-reactivating RTK RPPA class, which makes up 21% (3 of 14) of the current samples. These data reinforce the thesis that RPPA can serve both validating and complementary roles in resistance mechanism assessments.

Key advantages of RPPA and other protein array technologies (40, 41) are their relatively low cost per sample and rapid turnaround. Indeed, RPPA is already being used to direct clinical trials (NCT01074814, NCT01919749, and NCT02008994; https://clinicaltrials.gov/) in a clinically relevant time frame. The identification of a smaller, BRAF inhibitor resistance–specific antibody set from larger follow-up cohorts will further streamline the process and decrease the amount of input protein required. Ultimately, information of this type may inform second-line clinical decisions, although follow-up analyses of larger patient cohorts and additional functional studies are necessary. Such possibilities include the addition of ERK or pan-RAF inhibitors (e.g., TAK-632, ref. 13, or LY3009120) in melanomas displaying the hyper-MAPK pattern. Hyper-MAPK tumors might also specifically benefit from intermittent BRAF inhibitor dosing, on the basis of the observation that very high levels of MAPK signaling, further boosted by treatment stoppage, can be detrimental to melanoma cells (42, 43). The RTK cluster also provides the opportunity to illuminate the role of tumor-stroma interactions and potentially offers novel therapeutic opportunities (44, 45). Although it remains unclear by what mechanism(s) these RTKs are activated in the resistant tumors, many factors including microenvironment, epigenetics (46), and feedback regulations (47) are known to drive broad RTK activation in cancer. The data presented in this study represent an initial effort to establish a robust adjunct system for resistance assessment and may also apply to other contexts, such as biopsy-accessible EGFR or ALK inhibitor–resistant lung cancer.

Additionally, the flexibility and faithfulness of the mouse model further underlines the value of cross-species analyses in clinically oriented studies. The iBIP mouse model of melanoma described here closely resembled human samples in terms of the molecular response to BRAF inhibition and provided statistical power in the absence of additional human samples. Such GEM models offer a rapid and tractable system of longitudinal analyses, perturbation of the system, and, eventually, the testing of novel hypotheses. Thus, we envision co-clinical trials using cohorts of either biopsied resistant mice sorted by mechanism or iBIP variants engineered to model specific resistance mechanisms, with each cohort being tested against appropriate therapies. Insights from such studies in the mouse can potentially facilitate the design of human clinical trials.

Overall, this study brings together a co-clinical analysis through multiplex platforms that highlight the need for multiple dimensions of analysis in making clinical decisions and the utility of early serial tumor biopsies to proactively inform therapeutic interventions. Additionally, it provides a potential accessory for enabling the selection of second-line therapies, even when the specific underlying resistance-conferring mutation(s) are unknown. Finally, it underlines the potential to molecularly predict patient responses even prior to therapy, which may guide the timing and choices of additional targeted and/or immune therapies. Such a personalized view of resistance-specific cancer medicine may be a necessity in the dawning era of targeted therapies.

Patient samples. Twenty-seven patients with metastatic stage IV melanoma and positively genotyped for a BRAF^V600 mutation (Cobas 4800 BRAF V600 Mutation Test; Roche Molecular Systems or Snapshot) were enrolled in clinical trials for treatment with either a BRAF inhibitor (vemurafenib) or a combination of BRAF and MEK inhibitor (dabrafenib plus trametinib). Biopsies were immediately divided for snap-freezing, formalin fixation, and immediate extraction of RNA. One slide of formalin-fixed tissue was stained with H&E, and tumor and stroma percentages were estimated by a pathologist. For assessment of response, patients underwent a CT scan every 3 months. Responses were determined according to RECIST, version 1.1.

Animal husbandry and treatments. Generation of the TetO-BRAF Ink^–/– Arf^–/– mouse was described previously (20) and was bred into a floxed-Pten model (22), along with Tyr-CreERT2 (23) and Rosa26-Lox-Stop-Lox-rtTA-IRES-GFP (48) alleles. All iBIP mice used in this study were 70% FVB and 30% mixed background. All mice used to generate BRAFi-resistant tumors had been intercrossed for 19 to 22 generations, thus creating an isogenic recombinant inbred line. To initiate tumors, iBIP mice were treated with 4-hydroxy-tamoxifen (70% Z-isomer, 30% E-isomer, H6278; Sigma-Aldrich) dissolved in 100% EtOH. One microliter of 100 μM 4-hydroxy-tamoxifen was applied once on the tip of each ear. Mice were then continually administered doxycycline either through the drinking water (2 mg/ml) (D43020; Research Products International) or in the chow (200 mg/kg) (S3888; Bio-Serv). No differences were seen in tumor penetrance or latency between the water and chow vehicles. PLX4720 chow (Research Diets Inc.) and its corresponding control diet replaced the normal mouse chow in trial cohorts. For oral gavage, PLX4720 was dissolved in DMSO, then diluted 1:9 in 1% carboxymethylcellulose and 0.4% Tween-80 (Sigma-Aldrich). Progression was defined as the point at which a tumor reached double the size of its smallest volume. For iBIP allografts, iBIP cell lines were established in RPMI, 10%FBS, 2 mg/ml doxycycline, and 10⁵ cells were injected intradermally into nude (Taconic) or syngeneic immunocompetent mice. A single cell line, YC1474, was selected for downstream microarray and RPPA analyses in order to minimize noise and increase replicability across samples.

RPPA. Mouse tumors were homogenized to extract protein in radioimmunoprecipitation assay (RIPA) buffer plus phosphatase and protease inhibitors (Sigma-Aldrich) and sonicated to shear genomic DNA. Samples were aliquoted and stored at –80°C. Human protein lysates were isolated from snap-frozen patient tumors and denatured by 1% SDS with β-mercaptoethanol. Samples were diluted in five 2-fold serial dilutions in lysis buffer containing 1% SDS. Mouse samples were diluted to 1 μg/μl in RIPA buffer. RPPA was performed at the MD Anderson RPPA Core facility using 30 μg protein per sample. All antibodies were validated by Western blot analysis (49) at the RRPA Core facility, including for species specificity. Consequently, for human-only RPPA heatmaps, all human-compatible antibodies are shown, and for heatmaps containing mouse samples, mouse-incompatible antibodies were filtered out.

Purification of total RNA. Fresh human patient biopsies were homogenized and disrupted using a mortar and pestle, followed by use of a QIAshredder (QIAGEN). A QIAcube (QIAGEN) was used to harvest RNA from both patient biopsies and cell lysates using the RNeasy Mini Kit (QIAGEN).

Hierarchical clustering. All hierarchical clustering was performed using MultiExperiment Viewer software (http://www.tm4.org/mev.html). All samples were normalized against the median of the respective controls for all heatmaps.

Correlation analysis. The Pearson correlation value was calculated for every gene in the RNA-seq data, and the genes were ranked from the highest positive to the highest negative correlation. The data were then run through the GSEA preranked program using the c5 gene ontology (GO) gene set (http://www.broadinstitute.org/gsea/index.jsp). A conservative FWER significance cutoff of less than 0.25 was selected to avoid false positives. To condense each pathway into a point for presentation, the median was calculated across all pathway genes within each median-centered sample. The 22 proliferation and 94 immune pathway genes are listed in Supplemental Table 5. The 22-gene proliferation signature was generated by first intersecting 2 expression datasets for which the proliferation status was known: our previous tumor dataset with matched proliferation IHC data in the NCBI’s Gene Expression Omnibus (GEO) database (GEO GSE39984) and a dataset in which proliferation was specifically inhibited through CDK4 knockdown (GEO GSE8866). This yielded 147 genes. We then further tested and refined this signature using a number of expression datasets with known proliferation statuses (e.g., GEO GSE39646 and GSE35230; TCGA tumor versus normal) to identify a 22-gene set that consistently correlated with more highly proliferative samples by a Student’s t test and with each other in additional datasets (e.g., GEO GSE7553 and GSE46517).

Cell culture and plasmids. We isolated iBIP melanoma cell lines from primary tumors by collagenase digestion, followed by continual culture in RPMI and 10% FBS supplemented with 2 mg/ml doxycycline (Research Products International). The cell line was always used as a low-passage (less than 8) culture for both in vitro and in vivo experiments.

A375 and WM88 were Mycoplasma-free, BRAF-mutant human melanoma cell lines maintained in RPMI and 10% FBS. ORFs encoding the candidate RTKs and controls were obtained in pDONR vectors without tags. All ORFs were engineered into the EF1a promoter–driven pHAGE vector using Gateway Technology (Invitrogen). This vector contains an IRES-GFP. After viral transduction, all sublines were flow sorted for GFP on the same gating to obtain uniform GFP expression levels across the sublines.

Real-time PCR. DNA or cDNA was measured for BRAF copy number or expression, respectively, using real-time PCR. The BRAF-specific primers used were: BRAF-F2, GCTACAGAGAAATCTCGATGGAGT and BRAF-R3, GGTAACTGTCCAGTCATCAATTCA.

Sequencing of mRNA. The principle of transcriptome sequencing and its use for mutation detection and evaluation of gene expression have been described previously (50). Data generation was designed based on Illumina chemistry. Starting with approximately 700 ng total RNA, poly-A–containing mRNA molecules were converted into libraries suitable for next-generation sequencing approaches using reagents and the protocol provided by the Illumina TruSeq RNA Sample Prep Kit, version 2. Unique molecular identifiers were used for each sample. The 12 cycles of PCR were performed to selectively enrich the libraries, followed by size selection, stringent quality control, and qPCR quantification. Then, samples were pooled in equal volumes with 3 per lane, and paired-end 2 × 75 bp sequencing was performed using the Illumina HiSeq 2000. The mean sequencing quality was 37, the library insert size 150 bp, the mapping rate was 99%, and expression profiling efficiency was 0.79. The RNA-seq data are deposited in the European Genome-phenome Archive (EGA S00001000992).

Whole-exome sequencing. In brief, genomic libraries from approximately 500 ng genomic DNA (gDNA) were prepared using the standard KAPA paired-end sample preparation kit (KAPA Biosystems) according to the manufacturer’s instructions. SureSelect Human All Exon Kit, version 4 (Agilent Technologies) was used to enrich sequencing libraries for exomes. Samples were pooled 2 per lane, and paired-end 2 × 75 bp sequencing was performed using the Illumina HiSeq 2000, achieving over 100 times the sequencing depth.

Statistics. All significant differences were performed using a 2-tailed Student’s t test. A standard significance cutoff of 0.05 was used.

Study approval. All samples were obtained from participants who signed an informed consent form. This protocol was reviewed and approved by the Dana-Farber Cancer Institute (DFCI) IRB, in accordance with the applicable Code of Federal Regulations (CFR) set forth in 45 CFR, part 46, and 21 CFR, parts 50 and 56. All relevant clinical trials are registered at ClinicalTrials.gov (https://clinicaltrials.gov/) under the following trial numbers: NCT01006980, NCT01107418, NCT01264380, NCT01248936, NCT00949702, and NCT01072175. All mouse procedures were approved by the IACUC of the University of Texas MD Anderson Cancer Institute.