Acquired Resistance and Clonal Evolution in Melanoma during BRAF Inhibitor Therapy

Approximately 50% of metastatic melanomas harbor BRAF^V600 mutations, most commonly a V600E substitution (1), which hyperactivate the mitogen-activated protein kinase (MAPK) pathway and result in oncogene addiction. BRAF inhibitors, such as vemurafenib and dabrafenib, can elicit highly reproducible objective antitumor responses in the majority of treated patients, but the durability of response is short (median, 7 months) and limited by acquired drug resistance, resulting in disease progression (DP; refs. 2–5).

Earlier studies have identified specific mechanisms of BRAF inhibitor resistance in melanoma (6–17), which may reflect two different biologic processes: (i) early intrinsic resistance or adaptive tumor response, and (ii) late acquired resistance. An example of the adaptive response is the upregulation of receptor tyrosine kinases. Reported mechanisms of acquired BRAF inhibitor resistance in melanoma support the notion of MAPK pathway reactivation (18). Key mechanisms include emergence of mutant BRAF-concurrent RAS (7) or MEK (13) mutations and mutant BRAF amplification (11) or alternative splicing (8, 11), but the relative contribution of these mechanisms to clinical disease progression is unknown. Moreover, a clinical strategy to mitigate acquired BRAF inhibitor resistance by combined BRAF and MAP–ERK kinase (MEK) inhibition, although prolonging tumor suppression, is still beset by acquired drug resistance (19), suggesting MAPK-alternate escape route(s).

Cancer genomic heterogeneity in individual patients and tumors exists on a time-and-spatial continuum, an understanding of which has important consequences for personalized cancer medicine (20, 21). How melanoma genomic heterogeneity under BRAF inhibitor–selective pressure contributes to acquired resistance is unknown. This clonal evolutionary process in response to targeting of a dominant oncogene addiction pathway has not been characterized at the landscape or single-nucleotide levels. Defining the major molecular lesions (both known and novel), their core pathways, and the extent of melanoma genomic diversification underlying acquired BRAF inhibitor resistance represents a key benchmark for the further clinical development of BRAF inhibitor-based therapeutic strategies.

We analyzed 100 tumor samples from 44 patients (median progression-free survival or PFS = 145 days; range, 84–489; Table 1; Fig. 1A; Supplementary Fig. S1; Supplementary Table S1) whose melanomas developed acquired resistance to either vemurafenib or dabrafenib monotherapy. These specimens consisted of 29 baseline (pretreatment) tumors and 71 progressive tumors (Table 1). We first performed detection of known mechanisms of acquired BRAF inhibitor resistance, with the denominator of samples analyzed for each mechanism being slightly variable due to a small number of disease-progressive tumors lacking patient-matched baseline tumors (for the detection of relative BRAF gene copies) or lacking RNA samples (for the detection of BRAF splicing; Supplementary Table S2). We also generated 87 whole-exome sequence (WES) datasets (21 normal tissues, 22 baseline tumors, and 44 disease-progressive tumors) only from patients who donated at least a triad of these tissues. Among these patients, 16 donated multiple geographically and/or temporally distinct disease-progressive melanoma biopsies for this study, providing an opportunity to investigate tumor heterogeneity. We achieved a mean 107× coverage/base, 93.3% exome coverage at ≥ 15× (Supplementary Table S3), and a mutation-calling specificity of 94.4% as estimated by Sanger resequencing (Supplementary Table S4 and Supplementary Fig. S2). Consistent with our earlier study (7), among 56 of 71 (79%) progressive tumors available for analysis by deep sequencing of all 18 BRAF exons, no (0 of 56 or 0%) BRAF^V600E/K-secondary mutations were detected, but persistence of the same BRAF^V600E/K mutations detected at baseline was noted in all 71 of 71 (100%) progressive tumors; that is, BRAF inhibitor therapy did not select for minor, preexisting “contaminating” wild-type (WT) BRAF clones (Supplementary Table S2).

Among MAPK-reactivating mechanisms associated with acquired BRAF inhibitor resistance (Supplementary Table S2), NRAS mutations (G12D/R, G13R, and Q61K/R/L) were detected in 13 of 71 (18%) progressive tumors and KRAS mutations (G12C, G12R, and Q61H) in 5 of 71 (7%) progressive tumors (Supplementary Fig. S3). Mutant BRAF amplification (2–15-fold or 4–75 copies) was detected in 11 of 57 (19%) progressive tumors with patient-matched baseline tumors, and mutant BRAF alternative splice variants were found when novel exon boundaries (between 1 and 9, 1 and 11, and 3 and 9) were detected in 6 of 48 (13%) progressive tumors with available patient-matched baseline tumor RNA. Also, MAP2K1-activating mutations [K57N (22) and C121S] were detected in 2 of 71 (3%) progressive tumors, with the latter detectable only by WES but not by Sanger sequencing due to a low mutant allelic frequency (Supplementary Tables S2 and S3 and Supplementary Fig. S2). From the WES dataset, we defined recurrent CDKN2A loss (3 of 44 or 7% of progressive tumors) as a MAPK-reactivating mechanism, given that one of its gene products, p16/INK4A, negatively regulates the MAPK pathway effector complex consisting of cyclin D and cyclin-dependent kinase 4, which was implicated in BRAF inhibitor resistance (23). The minimal common region of these deletions harbored CDKN2A but few other genes (Supplementary Fig. S4). For all the disease-progressive samples where MAPK-reactivating molecular alterations were detected, their relative distribution of detection frequencies is estimated in Fig. 1B. More broadly among all the disease-progressive samples, RAS, MAP2K1, and CDKN2A mutations as well as mutant BRAF amplification or alternative splicing were detected at a 70% frequency and thus composed a core pathway of acquired BRAF inhibitor resistance through MAPK reactivation (Fig. 1C and Supplementary Table S5).

WES data enabled nomination of the PI3K–PTEN–AKT melanoma pathway as a second core resistance pathway (Supplementary Fig. S5 and Supplementary Tables S2, S5, and S6). AKT1/3 mutations (Q79K and E17K; Fig. 1D) were discovered in 2 of 44 progressive tumors subjected to WES and 0 of 27 remaining progressive tumors by Sanger sequencing of exons 3 and 4 [containing these plekstrin homology domain (PHD) mutations]. AKT3 copy number increase was not detected in any progressive tumor (Supplementary Table S2). Mutations in additional PI3K–AKT positive-regulatory genes (PIK3CA and PIK3CG) and in negative-regulatory genes (PIK3R2, PTEN, and PHLPP1; Fig. 1D–F) were detected in 10 of 44 progressive tumors. Overall, PI3K–PTEN–AKT pathway mutations, all validated by Sanger sequencing, constituted a second core acquired resistance pathway (at a 22% frequency), which overlapped with the MAPK core pathway (Fig. 1C).

We studied in depth the altered signaling effects of several genetic lesions in the PI3K–PTEN–AKT pathway by analyzing their predicted structural impacts (Fig. 2). AKT1^Q79K has been detected in breast and endometrial carcinoma (24), but its functional impact has not been characterized. On the basis of structure modeling (Fig. 2A), AKT1^Q79K in the PHD moved toward negatively charged phosphates of the Ins(1,3,4,5)P₄ by 4.02 Å when bound, suggesting enhanced lipid binding as K79 moved toward the 4- and 5-phosphates of PIP₃. In apo AKT1^Q79K, the acidic E17 forms an ionic interaction with basic K14, “tightening” the closed conformation. With Ins(1,3,4,5)P₄ binding, R86 moved toward the 4-phosphate and K14 toward the 3- and 4-phosphates (25). Superimposition of the apo- and ligand-bound conformations showed that Q79K moved toward the 4- and 5-phosphates in the latter. Notably, the PHD of PDK1, which contains a positively charged arginine (R521) at its AKT1^Q79-homologous position, showed a higher affinity for PtdIns(3,4,5)P₃ than for PtdIns(3,4)P₂ (26), consistent with the prediction that AKT1^Q79K engages the 5-phosphates of PIP₃ during lipid binding. Indeed, a PHD domain harboring the Q79K substitution displayed enhanced recruitment to the cell surface compared with the WT PHD (27). The somatic heterozygous AKT1^Q79K mutation (Supplementary Fig. S5) detected in melanoma during disease progression on BRAF inhibitor therapy is thus predicted to be a gain-of-function mutation.

PIK3CA^D350G has been detected in breast, endometrial, pancreatic, and colorectal carcinomas (24). Our structure modeling predicted that this somatic heterozygous mutant detected in disease-progressive melanoma is also a gain-of-function mutation (Fig. 2B), in that the D350G substitution disrupted a critical interaction with a conserved serine (S565 in PIK3R2) in its negative regulator such as PIK3R1 (p85α) or PIK3R2 (p85β). On the other hand, we detected a somatic PIK3R2 heterozygous mutation (N561D), which was predicted to disrupt a hydrophobic interaction between PIK3R2^N561 and PIK3CA^N345 (Fig. 2C). PIK3R2 mutations, including N561D, are frequently detected in endometrial cancer (28). In addition to PTEN large-scale deletions and frameshift mutations (resulting in early termination), a PTEN^M134 single amino acid deletion was detected in a melanoma with acquired BRAF inhibitor resistance. This M134 deletion was predicted to disrupt the P loop that is critical for its phosphatase activity (Fig. 2D). Consistent with a critical role for PTEN^M134, this mutation recurs in both endometrial and colorectal carcinomas (24).

We then directly analyzed these predicted functional mutants in the PI3K–PTEN–AKT pathway by introducing them into melanoma cell lines and detecting activated AKT (p-AKT Thr308; Fig. 3A). Stable overexpression of AKT1^E17K or AKT1^Q79K and AKT3^E17K in the M229 BRAF^V600E human melanoma cell line increased p-AKT Thr308 levels, whereas WT AKT1 or AKT3 overexpression did not or did so only minimally. In addition, stable overexpression of disease progression-specific PIK3CA mutants (D350G and E545G) and mutant PIK3R2 (N561D) increased the p-AKT level. Three types of PTEN genetic hits correlated with acquired BRAF inhibitor resistance: PTEN copy number (CN) loss (Fig. 1F), frameshift mutation (fs 40), and single codon deletion (M134del; Fig. 1D). The effects of these alterations were tested by stable knockdown of WT PTEN expression in M229 cells, which increased the p-AKT level. In contrast, stable overexpression of WT PTEN, but not PTEN^M134del, in the BRAF^V600E human melanoma cell lines M249 and WM2664 suppressed the p-AKT levels. Melanoma cultures (VUB164 MEL A, B, and C) derived from patient #43 tumor biopsies (baseline, disease-progressive tumors 1 and 2) displayed the same PTEN WT (baseline) and M134del (DP1 and 2) genotypes (Supplementary Fig. S6). Consistently, VUB164 B and C, relative to VUB A, displayed elevated levels of p-AKT, which was reversed upon the stable overexpression of WT PTEN. In available formalin-fixed paraffin-embedded (FFPE) tissues, the detection of AKT1^Q79K and AKT3^E17K, PTEN CN loss or frameshift, and PIK3R2^N561D was associated with elevated levels of p-AKT (Fig. 3B).

We also tested whether the detected functional molecular alterations would alter BRAF inhibitor sensitivity (Fig. 3C). Although AKT1 WT overexpression had no effect on vemurafenib sensitivity, AKT1^Q79K or AKT1^E17K and AKT3^E17K overexpression conferred vemurafenib resistance. PTEN knockdown in the PTEN WT M229 cell line conferred vemurafenib resistance, and PTEN reintroduction into the PTEN nonexpressing WM2664 cell line conferred vemurafenib sensitivity. Also, overexpression of PIK3CA^D350G and PIK3CA^E545G, as well as the positive control mutant PIK3CA^E545K, conferred vemurafenib resistance when compared with the vector or PIK3CA WT. Overexpression of WT PIK3R2 expectedly suppressed the baseline p-AKT level, whereas overexpression of PIK3R2^N561D not only did not suppress the baseline p-AKT level but also elevated it in WM2664 cells (Fig. 3A). Accordingly, when treated with vemurafenib, WM2664 cells overexpressing PIK3R2^N561D were more resistant to BRAF inhibition compared with WM2664 expressing with the vector or WT PIK3R2 (Fig. 3C).

A molecular lesion in the MAPK and the PI3K–PTEN–AKT pathways could be detected in 70% and 22% of the progressive tumors, respectively, but tumor heterogeneity was notable, suggesting clonal heterogeneity and/or collaborative mechanisms (Fig. 1C). In 9 of 44 (20%) patients, at least two mechanisms of resistance were detected in the same patient (same or distinct progressive tumors; Supplementary Table S2). This is an underestimate of the true extent, as only one progressive tumor was sampled in most patients and some progressive tumors were not amenable to all types of analysis. In 8 of these 9 patients, alterations in both core pathways were detected. In the subset of patients from whom multiple progressive samplings were feasible (n = 16 patients; Fig. 4A), 13 of 16 (81%) of patients harbored multiple mechanisms of resistance.

To investigate temporal and spatial tumor heterogeneity underlying acquired BRAF inhibitor resistance, we performed a whole exome-wide phylogenetic analysis of tissue biopsies from patient #37 (Fig. 4B and C) and other patients (Fig. 4D). Patient #37 had multiple cutaneous metastases and was started on the BRAF inhibitor dabrafenib in November 2010. On day 14, every single metastatic nodule had flattened, and the antitumor response was durable until day 383. Biopsies were performed on two baseline, two residual disease (RD1/2), and nine disease-progressive tumors. DP8 and 9 were the first subcutaneous (vs. cutaneous) metastases and presaged the onset of neurologic symptoms and brain metastasis in December 2012 (Fig. 4B).

Whole-exome, deep sequencing of genomic DNAs from blood, two baseline tumors, and nine progressive tumors that spanned 726 days of BRAF inhibitor selection resulted in a mean coverage of 117.5× (SD = 38.7×; range, 44.1–175.5×). To infer clonal relationships among progressive and baseline tumors, we collated unambiguous somatic mutations that were private or unique in at least one progressive tumor. This resulted in a collection of 360 single-nucleotide variants (SNV) and five INDELs that engendered a parsimony-based phylogenetic tree (Fig. 4C). The root node (origin of radii) represented the last common ancestral (LCA) node harboring 2,393 SNVs and 12 INDELs (with respect to the normal sample) that were shared by all the baseline and progressive tumors. Without exception, all progressive tumors followed branched (vs. linear) evolution in this patient and others studied (Fig. 4D and Supplementary Fig. S7 and S8). The average distance of a baseline tumor from the LCA node (5.3 SNVs) of patient #37 was less than that of a progressive tumor from the LCA node (39.3 SNVs), suggesting clonal diversification associated with an increased mutational burden in certain aberrant clones. Interestingly, the extent of genetic diversification of the progressive tumors was not colinear with their timing of clinical emergence.

At least five drivers of acquired BRAF inhibitor resistance accounted for this patient's clinical relapse, including four characterized mechanisms and at least one additional but unknown mechanism (presumably shared by DP8/9). Thus, a single disease-progressive biopsy would have revealed only 1 of 5 (20%) drivers of acquired resistance. These ancestral relationships were further supported by several observations: (i) three KRAS mutations were identical, supportive of DP2/5/6 diverging from a common, most recent ancestral node, (ii) two mutant BRAF amplification events were similar in amplitude and amplicon size (Supplementary Fig. S9), supportive of DP1/3 evolving from a common, most recent ancestral node, and (iii) two mutant BRAF alternative splicing events were associated with distinct novel exon–exon boundaries, supportive of convergent evolution for DP4/7 (Fig. 5A).

To examine how the genetic and mechanistic heterogeneity (Fig. 5A) correlated with tumor phenotypes, we assessed the relative levels of cell-cycle progression marker Ki-67 (Fig. 5B) and the activated MAPK signaling marker phosphorylated extracellular signal–regulated kinase (p-ERK; Fig. 5C). Interestingly, when compared with the baseline tumor, the residual tumors (RD1 and RD2) displayed little to no Ki-67 staining, suggesting tumor proliferative dormancy. Invariably, the fraction of nuclei with Ki-67 staining was dramatically higher in the progressive tumors compared with that in the baseline tumor in patient #37 and the majority of other tumor pairs (Supplementary Fig. S10). This suggests that BRAF inhibitor-related genomic diversification, altered mutational load (e.g., DP9), mutagenic processes (e.g., DP3/6/9), and/or epigenetically alterations (e.g., DNA methylation-driven differential gene expression; unpublished data) may contribute to this observed increase in tumor fitness. Specifically, the extensive genetic divergence of DP3 (marked by LOHs; data not shown) and DP9 (marked by SNVs) correlated with the highest proliferative indices. The highest levels of nuclear p-ERK were observed in those progressive tumors in which KRAS mutations were detected (Fig. 5A). Interestingly, DP8/9, with unknown drivers of BRAF inhibitor resistance, displayed minimal nuclear or cytoplasmic p-ERK staining.

Metastatic melanoma is marked by one of the highest mutational burdens seen across human malignancies, and UV-induced DNA damage, especially C > T transition mutations occurring within a dipyrimidine context, contributes to a dominant (>90%) fraction of all nsSNVs (29). Hence, we investigated the mutational spectra before and after BRAF inhibitor therapy. We found that both baseline tumors from patient #37 (as well as all 22 baseline tumors with WES data) harbored a predominance of C > T transition mutations (Fig. 6A) with preference for a dipyrimidine motif (Fig. 6B). In contrast, the mutation spectra of disease progression–specific SNVs were significantly altered, with a relative reduction in C > T transitions and increase in other transition and transversion mutations (Fig. 6A). Moreover, C > T transitions found in the progressive tumors displayed an attenuated dipyrimidine motif (Fig. 6B and C), indicating disparate, non–UV-related mutagenic processes.

Our study supports MAPK reactivation as a major pathway of acquired BRAF inhibitor resistance in melanoma, but also uncovers multiple genetic “hits” in the PI3K–PTEN–AKT pathway, nominating it as a second, bona fide core pathway of late drug resistance. Tumor heterogeneity contributed substantially to acquired BRAF inhibitor resistance and occurred at the levels of geographically and temporally distinct escape subclones within the same patient and tumor. Among nine disease-progressive tumor biopsies in patient #37, 277 of 360 (77%) of high-confidence somatic SNVs used to reconstruct phylogeny were unique to a single disease-progressive melanoma.

Genomic diversification culminated in multiple molecular lesions driving drug resistance and likely malignancy. Multiple mechanisms of acquired resistance being detected in the same tumor biopsy imply collaborative roles or tumor subclonal heterogeneity. We expect more heterogeneous and functional genetic variants in the MAPK and PI3K–PTEN–AKT pathway components to be described as more patients and their biopsies are subjected to in-depth sequence and molecular analyses. Certain novel genetic variants, such as PHLPP1^K596E, may affect both core pathways. Although PHLPP is a phosphatase for AKT and S6K, it may also repress phosphorylation of ERK (30). Extensive branched evolution, implying a paucity of common driver mutations, challenges efforts to monitor solid cancer heterogeneity, which may be aided by recent advances in liquid biopsies of circulating tumor cells (31) or DNA (32) coupled with genomic analysis. Evidence of non–MEK-dependent survival and increasing clonal fitness in response to BRAF inhibitor therapy may help explain why (i) the combination of BRAF and MEK inhibitors in BRAF inhibitor progressors can induce only low rates of secondary responses with limited durability (33), and (ii) upfront treatment with the combination of BRAF and MEK inhibitors in drug-naïve patients, albeit superior to either agent alone, is still beset eventually with acquired drug resistance (19).

In acquired BRAF inhibitor resistance, the proportion of melanomas harboring MAPK-reactivating molecular lesions was greater than that harboring PI3K–AKT-upregulating lesions. This relative distribution may stem from the differential early effects of therapy on these core pathways. Vemurafenib is known to suppress p-ERK in the majority of melanomas early during treatment (18, 34); hence, a majority of BRAF-inhibited melanomas would be under duress to reactivate p-ERK. In other work, we have shown that MAPK pathway inhibition can induce p-AKT levels within days of therapy and that gain-of-function AKT1 mutants most robustly conferred BRAF inhibitor resistance where this adaptive response was weak (27). Thus, the selective pressure for PI3K–AKT-upregulating genetic lesions may be tempered by tumor context-dependent, AKT-activating adaptive responses.

BRAF inhibitors are not known to be directly mutagenic to the melanoma genome but may engender secondary causes of DNA mutagenesis, such as reactive oxygen species (ROS) production (35) or alterations in DNA methylation (36) and histone modifications (both of which could influence mutation sites and rates (37, 38). More likely, the alterations in the disease progression–specific mutational spectra and signatures reflect preexisting minor (and lower fitness) subclones and/or the consequences of mutagenic processes (distinct from UV) more dominant throughout late-temporal malignant progression.

In conclusion, BRAF-mutant melanomas acquire BRAF inhibitor resistance via upregulation of both MAPK and PI3K–AKT pathways and genomic diversification resulting in branched evolution and increased tumor fitness. These data strongly support upfront cotargeting of both drug escape pathways in BRAF-mutant melanoma and imply a similar rationale to manage genetically complex human malignancies.

We evaluated 100 tumor biopsies from 44 patients. All 71 progressive tumors were examined for BRAF secondary mutations and known mechanisms of MAPK reactivation including RAS/MEK mutations by Sanger sequencing and mutant BRAF amplification by quantitative genomic DNA PCR or its alternative splicing by Sanger detection of novel exon–exon boundaries in the cDNAs (with exceptions noted in Supplementary Table S2). For analysis of known mechanisms such as RAS hotspot mutations, the denominator for analysis included all progressive tumor samples (inclusive of 24% of disease-progressive tumors without patient-matched baseline tumor samples), as such mechanisms have been well-validated previously to be generally disease progression-specific in the acquired resistance setting. To nominate novel genetic alterations causing acquired resistance, baseline (n = 22) coupled with progressive (n = 44) tumors from 21 of 44 (48%) patients with available normal tissues were WESed (Table 1). Genetic variants analyzed functionally were validated by Sanger sequencing. p-AKT (Cell Signaling Technology; #4060), p-ERK (Cell Signaling Technology; #4376), and Ki-67 (Dako; #M7240) levels were assessed by immunohistochemistry. All patients provided written informed consent.

For the MAPK-reactivating mechanisms (Fig. 1B), the detection frequencies of each mechanism (NRAS 18%, KRAS 7%, BRAF amp. 19%, BRAF alt. spl. 13%, MEK1/2 3%, CDKN2A 7%) were tallied, resulting in 67% (18 + 7 + 19 + 13 + 3 + 7 = 67). This normalized for the variable denominators or samples tested for each mechanism, as determined by the particular assay capability or tissue and WES availabilities (Supplementary Table S2). The detection frequency of each mechanism (e.g., RAS 18% + 7% = 25%) was then divided by the total detection frequency (67%) to derive the pie chart shown in Fig. 1B. Accounting of core pathway frequencies and their overlap incorporated samples with unknown mechanisms as shown in Fig. 1C (Supplementary Tables S2 and S4). Disease-progressive samples were designated as “unknowns” if all assays, including WES, were applied successfully. This determination avoided an overestimation of the unknown frequency.

For exome alignment, all samples were pair-end sequenced with length 2 × 100 bps, except for samples from patients #9 and #11 (2 × 76 bps reads; Supplementary Table S3). The sequences were aligned to the hg19 UCSC (University of California, Santa Cruz, Santa Cruz, CA) reference genome using the NovoAlign software (NovoCraft Tech. version 2.08.01). The aligned BAM were filtered for PCR duplicates, INDEL-realigned, and score calibrated using the GATK tools (39), Samtools and Picard software following the protocol listed on the GATK website (39). SNVs were called using a combination of the Unified Genotyper (UG) tool of GATK with a Fisher exact test postprocessing filter (named GATK-UGF; ref. 39), MuTect (40) and VarScan2 (41). INDELS (small insertion and deletions) were called using the combined calls of GATK-UGF, SomaticIndelDetector of GATK (IndelLocator; ref. 40) and VarScan2. We leveraged the combined strength of these programs by nominating all mutations that were called by at least two of the three SNV/INDEL callers, resulting in high-confidence calls. LOH calls on SNVs and INDELs were computed by intersecting calls by GATK-UGF and VarScan2 (as MuTect and SomaticIndelDetector do not report them), resulting in high-confidence LOH calls.

For GATK-UGF, we processed BAM files (normal of all patients, baseline, disease progression) using the GATK UG tool and, using one-sided Fisher exact test (P value cutoff ≤ 0.01), we retained only SNVs/INDELs whose mutant allele fraction at an SNV/INDEL position in the disease progression was significantly larger than that in the baseline and normal samples. Among these retained positions, disease progression–specific SNVs/INDELs were called when the genotype in the disease-progression sample was heterozygous while the genotypes in the baseline and normal were the same and homozygous. Baseline-specific SNVs/INDELs were defined similarly between the baseline and normal samples. DP-LOH events were called on baseline-specific SNVs/INDELs whose mutant alleles (non-normal and nonreference) were significantly enriched in the disease progression (P value cutoff ≤ 0.01). LOH on germline SNV/INDEL positions in the normal sample were excluded from analysis. To avoid artifactual SNVs/INDELs from low-coverage areas, we report only disease progression-specific and DP-LOH events with coverage ≥ 15 and SNV/INDEL score ≥ 30 (in phred scale) in the disease-progression, baseline, and normal samples.

We also ran VarScan2 on each disease progression BAM file with the baseline (first run) and with normal (second run) as the reference BAM file. We allowed low-frequency variants (frequency ≥ 0.05) to pass as long as they met the default P value cutoff. High-confidence calls were selected by the processSomatic and fpfilter.pl postprocessing in accordance with VarScan2. We collected the disease progression-specific calls by intersecting the high-confidence SNV, INDEL, and LOH calls from both runs. For MuTect-based SNV calling, we used the default parameters for both runs (DP vs. baseline and DP vs. normal) and collected the mutations that were annotated as “KEEP” in both runs as disease progression–specific SNVs. For GATK Somatic INDEL Detector-based INDEL calling, we likewise used the default parameters for both runs. Mutations that were annotated as “SOMATIC” in both runs were defined as disease progression–specific. Baseline-specific INDELs were similarly collected by comparing the baseline versus normal BAM files.

Shortlisted SNVs and INDELs were annotated by Oncotator (Broad Institute). Nucleotide and amino acid conservations were annotated by the phyloP (42) and polyphen2 (43) scores respectively. Mutations occurring in globular domain were listed on the basis of domain boundaries in INTERPRO release 39.0. The nominated disease progression–specific and DP-LOH SNVs/INDELs were reviewed manually and visualized using the Integrative Genomics Viewer. SNVs/INDELs in regions with high conservation (polyphen2 ≥ 0.85 or PhyloP ≥ 1.8) were further studied by functional assays and structure-based homology modeling. A subset consisting of 108 nonsilent, coding mutations was selected for Sanger resequencing to estimate the specificity of our SNV/INDEL-calling procedure.

Copy number variation (CNV) analysis intersected calls made by ExomeDepth (44) and ExomeCNV (45). From ExomeCNV, we used a fold ratio cutoff of 1.75 for amplification and 0.625 for deletion (2-fold amplification and deletion with a 25% normal tissue contamination). The CNV call of ExomeDepth was performed using default options, including the removal of common CNVs. Disease progression–specific CNVs were computed using the disease-progression samples as target and both the baseline and normal as references. The CNV calls were visualized using Circos.

For pathway annotations, we extracted canonical MAPK and PI3K–AKT pathway components (Supplementary Table S2) from the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway database. Our gene set of PI3K–PTEN–AKT pathway included the genes directly interacting with PI3K, PTEN, and AKT in KEGG's PI3K–AKT signaling pathway (KEGG ID: hsa04151). RAF1 (CRAF) was excluded as it was already included as a MAPK signaling pathway's gene. Our gene set of (classical) MAPK pathway includes the genes directly interacting with the RAS, RAF, MEK, and ERK genes in KEGG's MAPK signaling pathway (KEGG ID: hsa04151).

On the basis of crystal structure of apo and PIP₄-bound WT AKT1 (46), homology structures of apo and PIP₃-bound AKT1^Q79K were obtained by SWISS-MODEL (47) using apo (1UNP) and PIP₃-bound AKT1 (1UNQ) as templates. Structure alignment was calculated on the basis of the amino acid residues from Lys14 to Arg86 (PIP₄ binding segment). Homology structure of a heterodimer between PIK3CA (p110α) and PIK3R1/2 (p85α/β) was obtained on the basis of the template (PDB entry 3HIZ; ref. 48), and the alignments were based on full-length proteins. The PTEN^M134del structure was modeled by the I-TASSER online server (49) based on the full-length crystal structure bound to tartrate (PDB entry 1D5R; ref. 50).

We chose high-quality, unambiguous, disease progression–specific SNVs and INDELs and selected SNV/INDEL sites (not necessarily falling in coding region) that satisfied the following in all the samples: (i) coverage ≥ 15, (ii) for each site, at least one of the samples must have a different genotype call from the rest of the samples, (iii) sites with homozygous genotype must have at most three reads and at most 5% of any mutant allele (assumed to be caused by sequencing error), (iv) sites with heterozygous genotype must have mutant allele fraction ≥10% and ≥ four mutant allele reads (following a similar cutoff in VarScan2). These stringent criteria were meant to filter out ambiguous GATK genotype calls to produce an accurate phylogenetic tree. We chose the mutant allele (i.e., nonreference, non-normal) of each site as its representative haplotype and used it to build the most parsimonious phylogenetic tree using the PHYLIP program's dnapars subroutine. Tree drawing was done using the drawtree program in PHYLIP. Shared mutations shown in the phylogenetic trees were computed on the basis of the branching structure. We verified our tree construction method by confirming that the phylogenetic tree's branching structure inferred for Pt #37 stayed the same over different cutoffs of coverage depths, mutant allele fractions, and counts on both heterozygous and homozygous sites.

We categorized the disease progression–specific SNVs and baseline-specific SNVs based on types of nucleotide transitions and transversions. To compute the UVB signature, we collected the ±2 bases with respect to the C → T transition SNV sites using samtools' view command on the BAM alignment files and subsequently built a positional nucleotide frequency motif that was visualized using the seqLogo R package. The center nucleotide in the motif showed the WT allele, which, in the C → T transition case, was the nucleotide C. The character height in the motif logo was drawn proportionally to the information content of the position.

Stable expression and knockdown by lentiviral transduction were performed in the BRAF-mutant melanoma cell lines M229, M249, WM2664, and VUB164MEL (A, B, and C). The M cell lines were established at University of California, Los Angeles (UCLA; Los Angeles, CA) with Institutional Review Board approval and routinely authenticated by mitochondrial DNA sequencing. WM2664 and VUB164MEL cell lines were obtained from the Wistar Institute (Philadelphia, PA) and the Ludwig Institute for Cancer Research (Brussels Branch, Belgium) via Material Transfer Agreements and were not further authenticated except for verification of the BRAF-mutant status. Protein lysates were then analyzed by Western blotting (antibodies for p-AKT Thr308, p-AKT Ser473, total AKT, p-ERK1/2 Thr202/Tyr204, total ERK, and total PTEN were purchased from Cell Signaling Technology #4056, 4060, 4685, 9101, 9102, 9188, respectively; for FLAG, and TUBULIN from Sigma F1804 and T9026). Effect of vemurafenib (vs. dimethyl sulfoxide) on cell survival was measured by MTT (3 days of drug exposure) or clonogenic assays (10 days). The P values of compared survival curves were calculated by Graphpad prism 4.0. The dose values were log-transformed and normalized. Nonlinear regressions were performed with the sigmoidal dose–response model. Log EC₅₀ of each fitted curve was compared with that of vector control, and a P value of less than 0.05 (F test) was considered statistically significant.

K.B. Dahlman has received honoraria from the Speakers Bureau of Illumina. M.C. Kelley has received a commercial research grant from GlaxoSmithKline/National Comprehensive Cancer Network (NCCN). R.F. Kefford has received honoraria from the Speakers Bureaus of Roche and Novartis and is a consultant/advisory board member of Roche, GlaxoSmithKline, and Novartis. B. Chmielowski has received honoraria from the Speakers Bureaus of Genentech and Bristol-Myers Squibb and is a consultant/advisory board member of Merck and Genentech. J.A. Sosman is a consultant/advisory board member of GlaxoSmithKline. G.V. Long is a consultant/advisory board member of GlaxoSmithKline and Roche. A. Ribas has ownership interest (including patents) in Entrogen. R.S. Lo is an inventor on a patent. No potential conflicts of interest were disclosed by the other authors.

Conception and design: H. Shi, R.F. Kefford, B. Chmielowski, R.S. Lo

Development of methodology: H. Shi, W. Hugo, X. Kong, R.C. Koya, R.S. Lo

Acquisition of data (provided animals, acquired and managed patients, provided facilities, etc.): H. Shi, X. Kong, A. Hong, R.C. Koya, G. Moriceau, T. Chodon, D.B. Johnson, K.B. Dahlman, M.C. Kelley, R.F. Kefford, J.A. Sosman, N. van Baren, G.V. Long, A. Ribas, R.S. Lo

Analysis and interpretation of data (e.g., statistical analysis, biostatistics, computational analysis): H. Shi, W. Hugo, X. Kong, A. Hong, R.F. Kefford, B. Chmielowski, J.A. Sosman, G.V. Long, R.S. Lo

Writing, review, and/or revision of the manuscript: H. Shi, W. Hugo, R.C. Koya, K.B. Dahlman, M.C. Kelley, R.F. Kefford, B. Chmielowski, J.A. Glaspy, J.A. Sosman, N. van Baren, G.V. Long, R.S. Lo

Administrative, technical, or material support (i.e., reporting or organizing data, constructing databases): H. Shi, W. Hugo, X. Kong, R. Guo, D.B. Johnson, M.C. Kelley, G.V. Long, R.S. Lo

Study supervision: R.S. Lo

The authors thank the patients for their participation. They also thank G. Bollag and P. Lin (Plexxikon Inc.) for providing PLX4032, Gulietta Pupo (MIA), and Donald Hicks (Vanderbilt) for technical assistance, Pam Lyle (Vanderbilt), Richard Scolyer (MIA), and Sharona Yashar (UCLA) for pathology expertise, and A. Villanueva (UCLA) for clinical data management.

R.S. Lo is supported by a Stand Up To Cancer Innovative Research Grant, a Program of the Entertainment Industry Foundation (SU2C-AACR-IRG0409). Additional funding for R.S. Lo came from the National Cancer Institute (K22CA151638 and 1R01CA176111), Burroughs Wellcome Fund, American Skin Association, Melanoma Research Alliance, Sidney Kimmel Foundation for Cancer Research, Eli & Edythe Broad Center of Regenerative Medicine & Stem Cell Research, Ian Copeland Melanoma Fund, the Harry J. Lloyd Charitable Trust, and the National Center for Advancing Translational Sciences UCLA Clinical and Translational Science Institute Grant UL1TR000124. Funding for R.S. Lo and A. Ribas came from National Cancer Institute (1P01CA168585), The Seaver Institute, the Wesley Coyle Memorial Fund, and the Louis Belley and Richard Schnarr Fund. Funding for A. Ribas came from the Dr. Robert Vigen Memorial Fund, the Garcia-Corsini Family Fund, and the Ruby Family Foundation. Funding for H. Shi came from Association of American Cancer Institutes and American Skin Association. Funding for R.F. Kefford and G.V. Long came from the National Health and Medical Research Council of Australia and Translational Research Program of the Cancer Institute New South Wales. Funding for J.A. Sosman came from the American Cancer Society Melanoma Research Professorship and the NCI 5K24 CA097588 Midcareer patient oriented research award.