A Meta-analysis of Somatic Mutations from Next Generation Sequencing of 241 Melanomas: A Road Map for the Study of Genes with Potential Clinical Relevance

Melanoma is a malignant cancer of melanocytes. The disease accounts for 80% of deaths from skin cancer (1), with an estimated 76,690 new patients and 9,480 deaths expected in the United States in 2013 (1). The high mortality rate of advanced melanoma is largely because of its aggressive behavior and limited treatment options. As such, the 5-year survival rate for patients with advanced melanoma has historically been below 15% (2).

Melanoma has been traditionally classified based on clinical and histopathologic characteristics, such as growth pattern and anatomic site of origin. More recently, driver mutations in genes encoding signaling proteins that activate the mitogen-activated protein kinase (MAPK) pathway (3) have been implicated in the pathogenesis of a majority of melanomas. Accordingly, we recently developed a melanoma-specific multiplex mutational profiling assay to detect 43 recurrent mutations in 6 genes relevant to targeted therapy for the disease (4). This SNaPshot assay, used routinely in our clinic, detects mutations in BRAF (at position V600), NRAS (G12/13, Q61), KIT (W557, V559, L576, K642, D816), GNAQ (Q209), and GNA11 (Q209) in tumors. We chose these specific mutations because (i) they occur in melanoma at a frequency of greater than 1% and (ii) they have relevance to emerging or targeted therapies such as the BRAF inhibitors, vemurafenib and dabrafenib, the MEK inhibitor, trametinib, and the KIT inhibitor, imatinib. Using this assay, we have found that one third of tumors lack any of these mutations and are herein referred to as “pan-negative” (4). Continued investigation for novel driver mutations in these cases is critical to improve therapeutic outcomes for patients. For example, we recently found that ∼8% of “pan-negative” cases harbor non-V600 exon 15 BRAF mutations (5) and 4% to 8% contain activating BRAF fusions, both of which may be sensitive to MEK inhibition (6).

Rapid advances in next generation sequencing (NGS) technologies have made comprehensive characterization of cancer genomes feasible (7–9). A large number of whole exome sequencing (WES) and whole genome sequencing (WGS) studies have already been performed in melanoma (10–18). Here, we analyzed 241 paired melanoma tumor/normal tissue samples from 6 recently published WES and WGS studies to identify mutation signatures associated with melanomas with the most common driver genes (i.e., BRAF and NRAS) and those without known recurrent BRAF, NRAS, KIT, GNAQ, or GNA11 mutations. This analysis provides important insight into the molecular events associated with melanomas and has identified potential therapeutic targets among “pan-negative” melanomas.

Similar to data collection procedures in our NGS Catalog database (8), we conducted a comprehensive literature search of melanoma NGS studies using the keywords “exome sequencing AND melanoma” and “whole genome sequencing AND melanoma” through PubMed (http://www.ncbi.nlm.nih.gov/pubmed). We performed a careful manual check of the search results. Our query revealed at least 10 melanoma NGS studies published from 2010 to 2012 (as of September 2012, right before we started the analysis; ref. 8). Studies were excluded if only part of the NGS mutation data was available. The mutation data from (11) was not included in our study because only one tumor-normal pair was sequenced and it harbored the known driver mutation, BRAF V600E. Duplicate data were filtered by examining authors' names and affiliations and tumor name/ID. As a result, 6 melanoma WGS or WES studies (12, 14–17, 19) were collected for our meta-analysis (Fig. 1 and Supplementary Table S1). The sequencing quality of these melanoma genomes/exomes was high, with the validation rate estimated to be >95% in most of these studies.

The number of sequenced tumor samples varied among the 6 studies, ranging from 7 to 121 samples. Here, we only used the NGS data from the tumors that had matched normal tissues in the same study. In addition, 23 of the 25 melanoma samples in (14) were sequenced in another study (19), so these 23 duplicated samples in ref. 14 were removed in our study. The mutation rate is high in melanoma tumor genomes compared with other types of tumor genomes (9). Surprisingly, no somatic mutation data were identified in 10 melanoma samples in ref. 15, most of which (6/10) were mucosal or acral. Therefore, those samples were also excluded. In total, we analyzed NGS data from 241 tumor samples with mutation information, along with their matched normal samples (Fig. 1 and Table 1). Among them, 182 originated from cutaneous sites, 17 from acral sites, 7 from mucosal sites, 6 from uveal sites, and 29 from unknown primary sites (Supplementary Table S2).

The somatic mutation data were downloaded from the supplementary section of each of the 6 publications (12, 14–17, 19). For the current analysis, we primarily focused on somatic single nucleotide variants (SNV) because they constitute the largest fraction of oncogenic drivers identified in melanoma. For these SNVs, we collected all nonsilent parts, including missense, nonsense, nonstop, and splice site mutations, as well as accompanying information such as tumor name, tumor type, and gene name. We generally refer to these SNVs as mutations in this article.

NGS has enabled the identification of tens of thousands of somatic mutations in melanoma genomes (10). As a result, sifting through these mutations to pin down the few driver mutations present in melanoma is a challenge; therefore, we generated a list of candidate genes to further analyze in the published datasets. We first compiled a list of genes that were potentially related to melanoma from literature reports (2, 20–23). Because mutations in protein kinase genes have been frequently reported in cancers (24), and they are often the targets in the design of antitumor therapies (25), we also included the 612 genes from Agilent's SureSelect Human Kinome targeted NGS Kit (http://www.genomics.agilent.com/en/home.jsp) in our candidate gene list. The resulting list contained 632 genes for further analysis (see Supplementary File S1). Unless specifically noted, the analysis below was conducted using these candidate genes.

Currently, Vanderbilt University Medical Center utilizes the SNaPshot genotyping assay to screen for 43 specific point mutations in 6 genes (BRAF, NRAS, KIT, GNAQ, GNA11, and CTNNB1) in melanoma (4) that are relevant to current and emerging targeted therapies in melanoma except for CTNNB1. Excluding CTNNB1 (CTNNB1 mutations commonly co-occur with mutations in the other 5 genes), we analyzed the melanoma NGS data against these drivers to determine mutations associated with these 5 driver genes, as well as to uncover potential novel drivers in “pan-negative” samples [i.e., samples which lack all the known, recurrent mutations in BRAF (V600), NRAS (G12, G13, and Q61), KIT (W557, V559, L576, K642, D816), GNAQ (Q209), and GNA11 (Q209)]. In-house Perl scripts were developed to analyze these data and a single-sided Fisher exact test was used to assess the significance of mutation association.

To classify melanoma genomes according to our clinical SNaPshot-based assay, we queried WGS and WES data from 241 melanoma samples for known driver mutations in BRAF (V600), NRAS (G12/13, Q61), KIT (W557, V559, L576, K642, and D816), GNAQ (Q209), and GNA11 (Q209). Supplementary Table S2 summarizes the number of tumors, the tumor subtypes, and known driver mutation(s) that each tumor harbored. Briefly, 50.2% (121/241) tumors were found to harbor BRAF V600 mutations (Fig. 1). Among them, 86.8% (105/121) had V600E missense mutations. Fifteen had V600K mutations (12.4%) and 1 had a V600R mutation (0.8%). Forty-seven samples (19.5%) had NRAS mutations, including Q61 mutations [44/47 (93.6%): Q61R (22/47, 46.8%), Q61K (12/47, 25.5%), Q61L (6/47, 12.8%), and Q61H (4/47, 8.5%)] and G12 mutations [3/47 (6.4%): G12V (2/47, 4.3%) and G12D (1/47, 2.1%)]. No G13 mutations were detected. Three uveal melanoma samples (3/241, 1.2%) had GNA11 Q209L mutations. Only one tumor (1/241, 0.4%) had a KIT mutation (V559A). No mutations were found in GNAQ. The remaining samples were wild type (WT) for these mutation hotspots. These data indicate that the 241 pooled genomes were generally representative of a typical melanoma population (4, 5, 21).

We next sought to identify mutations that co-occurred with mutant BRAF in this dataset. From a biologic standpoint, co-occurring mutations could contribute to disease progression and/or metastasis. From a clinical standpoint, such mutations could modify clinical responses to single-agent–targeted therapies and/or contribute to drug resistance (26, 27). Table 1 lists mutations that occurred in at least 10% of BRAF-mutated melanomas with a P value < 0.05 (Fisher exact test). Genes were ranked according to their P values for different frequencies between samples with and without BRAF mutations.

Mutations in the gene TTN were the most significantly associated with BRAF mutation, occurring in 64.6% of samples (P = 0.009). Encoding a functional component of striated muscle, TTN is the longest gene in the human genome (24). TTN mutations are often seen in cancer NGS studies and are likely passenger mutations because of gene length bias (9). The next most frequently mutated gene was TP53, found in 21.5% of BRAF-mutated tumors (P = 0.011). Although TP53 is a highly recognized tumor suppressor gene harboring numerous mutations, it is rarely mutated in melanoma relative to other cancers (28). The association between BRAF and TP53 mutations supports a previous report (29) based on Sanger resequencing of selected genes in melanoma cell lines. We used MutationAssessor (30) to predict the functional relevance of these TP53 mutations to melanoma. All missense mutations in TP53 were predicted to have a medium impact on the protein function. Because a nonsense mutation may lead to the loss of function of the encoding gene, we further examined the distribution of nonsense mutations in TP53 in the samples with versus without BRAF mutations. BRAF mutant samples harbored more nonsense mutations than BRAF wild-type samples (5.4% vs. 2.7%). The third most commonly mutated gene was COL1A1, which encodes a component of type I collagen (occurring in 13.1% of BRAF-mutated tumors; P = 0.034). Thus far, there is no evidence suggesting that COL1A1 is associated with melanoma tumorigenesis.

Table 2 shows the mutated genes associated with NRAS mutations (occurred in at least 10% of NRAS-mutated melanomas and P < 0.05, Fisher exact test). Aside from TTN, this list of mutated genes was surprisingly very different from the list associated with BRAF mutant tumors. For example, alterations in the cell-cycle effector, PPP6C were frequently observed in NRAS-mutated melanomas, with a frequency of 17.7% and P = 0.011. This gene was recently found to be a novel candidate melanoma gene harboring putative somatic driver mutations (19). The genes encoding the Rho-GEF kinase, KALRN, and the phosphoinositide-3 kinase, PIK3R4, which ranked third and fourth in this association list based on P values, were mutated in 27.5% and 11.8% of NRAS-mutated melanomas, respectively. Neither of these 2 genes has yet been reported to be involved in melanomagenesis. The fifth gene was TRPM6 (27.5%, P = 0.020); its role in melanoma remains unknown, but other members of the Transient Receptor Potential ion channel subfamily M gene family (TRPM1, TRPM2, TRPM7, and TRPM8) have been proposed to play a crucial role in melanoma (31). The enterotoxin receptor, GUCY2C, which is an emerging prognostic molecular biomarker in colorectal cancer (32), was also found to co-occur with NRAS mutations (13.7%, P = 0.021) in melanoma. Finally, we observed that 13.7% of all NRAS mutations identified were coincident with mutations in PRKAA2 (P = 0.043), which encodes the catalytic subunit of AMP-activated protein kinase (AMPK). We further assessed the functional impact of the 8 PRKAA2 missense mutations that were found to be associated with the NRAS mutations. Four missense mutations in PRKAA2 were predicted to have a medium impact on the protein function by MutationAssessor. Shen and colleagues (33) reported that PRKAA2 (AMPK) is the long-sought kinase for BRAF Ser729. BRAF phosphorylation by PRKAA2 could attenuate BRAF signaling and inhibit cell proliferation. The mutation association we found here might suggest a new mechanism of tumorigenesis involving inhibitory mutations of PRKAA2 in BRAF wild-type, NRAS-mutant melanoma; however, future biologic analyses outside the scope of this study are needed to confirm this concept. In addition, a recent study (34) suggested that PRKAA2 plays a potential role in cancer progression and may represent a possible therapeutic target for the treatment of early gastric cancer.

For the above genes that were found to be significantly associated with BRAF or NRAS mutations, we manually reviewed their results in individual studies for reliability. Because most samples used in this study were from Hodis and colleagues (ref. 19; 121 paired tumor-normal samples) and Krauthammer and colleagues (ref. 15; 89 paired tumor-normal samples), we examined the mutations that co-occurred with mutant BRAF or NRAS only in these 2 studies (the details are available in Supplementary Tables S3–S10). We found some results in Table 2 were consistently observed, including PIK3R4, PPP6C, TRPM6, and KALRN using the Hodis and colleagues data (Supplementary Table S8) and TTN using the Krauthammer and colleagues data (Supplementary Table S10).

It is known that some driver mutations are associated with sun exposure whereas others are more commonly found in sun-protected areas. We explored whether the observed mutation association was because of ultraviolet (UV) induction or a specific driver mutation in genes such as BRAF and NRAS. Specifically, we analyzed the mutation association patterns in each melanoma subtype (acral, mucosal, uveal, cutaneous, or unknown type). Three genes (TTN, TP53, and COL1A1) that were associated with BRAF-mutant melanomas in the whole cohort did not co-occur with mutant BRAF in any melanoma subtype (Supplementary Table S11). This lack of association might be because of small sample sizes in each subtype. However, we did not observe a statistically significant association in cutaneous melanoma either, which had 109 samples. This preliminary analysis suggested that the genes we reported are more likely because of the association with the driver gene BRAF, rather than sun exposure.

Interestingly, almost all genes associated with NRAS mutation (except PRKAA2) in melanomas (Table 2) were consistently observed in cutaneous melanoma (Supplementary Table S12). Because cutaneous melanoma is mainly caused by UV damage, we could not reject the hypothesis that the NRAS-associated genes (TTN, PPP6C, KALRN, PIK3R4, TRPM6, and GUCY2C) are due to sun exposure. Considering that UV specifically causes C>T transitions (i.e., C changes to T), we further compared C>T mutation frequency in NRAS mutant samples versus NRAS wild-type samples. The C>T frequency in NRAS mutant samples (80.7%) was slightly larger than that in NRAS wild-type samples (78.9%). For the purpose of comparison, the C>T frequency in BRAF mutant samples (78.3%) was slightly smaller than that in BRAF wild-type samples (80.5%). Further studies are needed to determine the influence of sun exposure on the association of mutations in NRAS-mutant melanomas.

As stated above, our clinical SNaPshot-based assay examines tumors for mutations in BRAF (V600), NRAS (G12/13, Q61), KIT (W557, V559, L576, K642, D816), GNAQ (Q209), and GNA11 (Q209). A total of 69 tumors (28.6%) of the 241 tumor/normal pairs were “pan-negative” according to this mutation screening strategy (Fig. 1). To identify other potential mutations in these genes, we interrogated the 69 tumor/normal pairs for non–hotspot mutations in BRAF, NRAS, KIT, GNAQ, and GNA11 (Table 3, Fig. 2, and Supplementary Figs. S1–S4). Two samples with distinct BRAF L597 mutations were identified. This frequency (2/69 = 2.9%) is consistent with our previous report, in which we identified 2 BRAF L597 mutations in 49 SNaPshot pan-negative samples (4.1%) upon screening of BRAF exon 15 (5). Six additional non-V600 BRAF mutations were found in these 69 pan-negative melanoma samples, including K601E (2 samples; ref. 5), P75L (1 sample), G469R (1 sample; ref. 35), N581T (1 sample), G593S (1 sample), and P731S (1 sample; ref. 36). Some of these are already known to activate the MAPK pathway and confer sensitivity to MEK inhibition (5, 36). Thus, in this dataset, 10.1% (7 of 69) of pan-negative melanomas harbored non-V600 BRAF mutations, some of which are already suspected to benefit from treatment with targeted agents (37). This proportion (10.1%) was significantly higher than that in driver mutation-positive melanomas (3.5%, 6/172; P = 0.039, Fisher exact test). In addition, 3.3% (4/121) and 4.3% (2/47) of melanomas with BRAF V600 and NRAS G12/G13/Q61 mutations, respectively, harbored non-V600 BRAF mutations. These 2 proportions were not significantly higher than that of their respective wild-type samples harboring non-V600 BRAF mutations: P = 0.960, Fisher exact test, BRAFV600 melanomas versus non-BRAFV600 melanomas (7.5%, 9/120); and P = 0.761, Fisher exact test, NRAS G12/G13/Q61 mutant melanomas versus non-NRAS G12/G13/Q61 melanomas (5.7%, 11/194). The rate of non-V600 BRAF mutations in the whole cohort is 5.4% (13/241).

Two somatic nonsynonymous mutations were identified in NRAS (Q22K and H131R), and 3 somatic nonsynonymous mutations were identified in KIT (L572P, K638E, and N818Y). The NRAS mutations were found in disease from sun-damaged skin, and all KIT mutations were found in melanomas from acral sites. Additional GNAQ mutations included R183X, P185S, S211L, and A302G, of which P185S and A302G were present in a single patient. Interestingly, only 1 of these 4 GNAQ mutations originated from a uveal site; the other 3 GNAQ mutations were found in sun-damaged skin. Only 1 GNA11 mutation (R183C) was identified in a uveal melanoma.

The proteins encoded by GNAQ and GNAQ11 belong to the family of guanine nucleotide-binding protein α-subunits (GNA) (38), which play an important role in cell signaling as well as tumor initiation and progression (39). Recently, activating mutations in GNAS (another GNA-encoding gene) in breast cancer, colorectal cancer, and ovarian cancer were reported (40, 41). Further interrogation of the “pan-negative” melanoma samples revealed 19 mutations in 9 GNA-encoding genes other than GNAQ and GNA11 (Supplementary Table S13). Interestingly, 17 of the 19 GNA mutations were found in cutaneous melanomas, unlike the known recurrent mutations in GNAQ and GNA11, which are primarily associated with uveal melanomas (42, 43). Notably, 13 of these 19 GNA mutations are caused by C>T (G>A) mutations that are compatible with UV-mediated damage, providing evidence for UV exposure in melanoma pathogenesis. Among these 19 GNA-encoding genes, GNAI2 had 6 missense changes, of which 2 were present in a single sample. Taken together, these data indicate that GNA mutations potentially have a more significant role in cutaneous melanomas than previously thought. Future biologic studies are necessary to determine whether these GNA mutations are drivers.

To identify other potential novel driver mutations in pan-negative melanomas, we identified gene mutations that were enriched in the pan-negative samples (69) compared with the driver-positive samples (172). Twelve genes had statistically significant mutations in the pan-negative samples (mutated in ≥10% of pan-negative melanomas and P < 0.05, Fisher exact test; Table 4 and Supplementary Table S14), including ALK, STK31, DGKI, RAC1, EPHA4, ADAMTS18, EPHA7, ERBB4, TAF1L, NF1, SYK, and KDR (Supplementary Figs. S5–S16). Here, the order of the genes was based on the P values. Among them, RAC1 (15, 19), ADAMTS18 (44), ERBB4 (45), and NF1 (19) have already been reported as potential contributors to melanoma tumor progression, providing validation of the effectiveness of this analysis strategy. The most significantly mutated gene was the tyrosine kinase, ALK (altered in 17.4% of pan-negative tumors, P = 0.001; Supplementary Fig. S10). Among the 13 ALK missense mutations in 12 samples, 1 (E1197K) was predicted to have a high impact on the protein function by MutationAssessor and 4 were predicted to have a medium impact on the protein function. The remaining mutations had a low or no impact on protein function. All these missense mutations were C>T (G>A) transitions, which may be the result of UV light damage. Considering that activating ALK rearrangements are found in lung cancer (46) and more recently in colorectal cancer (47), and a subset of lung cancer patients with ALK fusion genes are sensitive to the kinase inhibitor crizotinib (48), the role of ALK in melanoma and the clinical efficacy of ALK inhibition in ALK mutation-positive melanoma warrants further investigation.

Clustering of mutations in these 12 genes indicated that 5 genes (RAC1, ADAMTS18, EPHA7, TAF1L, and NF1) harbored at least 1 recurrent mutation present in at least 2 pan-negative tumor samples. Mutations in the GTPase RAC1 occurred in 8 (11.6%) of the 69 pan-negative tumors compared with 5 of the 172 (2.9%) driver-positive tumors (P = 0.011; Supplementary Fig. S5). Seven of these 8 RAC1 mutations were the alteration P29S, which was previously shown to activate RAC1 (19) and increase melanocyte proliferation (15); therefore, RAC1 P29S may be a potential therapeutic target in melanoma. Notably, the strong recurrent mutation pattern in RAC1 implies that it might act as an oncogene instead of tumor-suppressor gene in melanoma. ADAM metallopeptidase with thrombospondin type I motif, 18 (ADAMTS18), which typically acts as a tumor-suppressor, was recently found to be mutated in melanoma, and mutant ADAMTS18 could increase melanoma cell migration and metastasis (44). We found ADAMTS18 mutations in 23.2% of the 69 pan-negative melanomas (Supplementary Table S14), with a hotspot R172Q mutation in 3 pan-negative tumors (Supplementary Fig. S6). The remaining 3 genes (EPHA7, TAF1L, and NF1) each contained 1 recurrent mutation in 2 pan-negative melanomas. EPHA7, which encodes ephrin receptor A7, has been shown to act as a tumor-suppressor in follicular lymphoma (49). We observed 14 EPHA7 mutations in 12 pan-negative tumors (17.4%, P = 0.017). Among them, a G642E mutation, present in the tyrosine kinase domain of EPHA7, occurred in 2 tumors (Supplementary Fig. S7). This mutation may affect kinase activity of EPHA7, making it a potential clinical target. Furthermore, we found another significantly mutated gene in the ephrin family, EPHA4 (10.1%, P = 0.015; Supplementary Fig. S13). Ephrin receptors have been proposed as potential biomarkers and therapeutic targets in prostate cancer (50). In TAF1L [TAF1 RNA polymerase II, TATA box binding protein (TBP)-associated factor, 210 kDa-like], we observed 15 mutations in 11 of the 69 pan-negative samples (15.9%, P = 0.022; Supplementary Fig. S8). TAF1L was the fourth most frequently mutated gene by statistical ranking among the kinase genes in another study (24).

Many loss-of-function mutations were observed in the genes STK31 (serine/threonine kinase 31) and NF1 (neurofibromin 1). We found 22 STK31 mutations in 18 pan-negative tumors (26.1%, P = 0.001; Supplementary Fig. S11), including 3 nonsense, 1 splice-site, and 18 missense mutations. STK31 is a highly conserved member of the serine/threonine kinase (STK) family. It has been shown that the expression of STK genes is frequently altered in human cancers (51). Considering that previous findings suggested a functional role for STK31 in regulating gastrointestinal (52) and colon cancer development (53), and considering that a variety of STK inhibitors are available, STK31 may represent a potent target in the treatment of melanoma. Moreover, 22 NF1 mutations, including 6 nonsense (include one recurrent mutation W784X), 4 splice-site, and 12 missense mutations, were observed in 12 pan-negative tumors (17.4%, P = 0.024; Supplementary Fig. S9). Four of these 12 missense mutations in NF1 were predicted to have a medium impact on the protein function by MutationAssessor. For the splicing mutations, we did not find any of them in the SpliceDisease database, which is a comprehensive database linking splicing mutations with diseases (54). As a negative regulator of RAS signaling, NF1 inactivation was recently reported to promote BRAF inhibitor resistance in the context of BRAF-V600 mutant melanoma (36). In pan-negative melanoma, NF1-inactivating mutations (nonsense mutations) may also play a role in promoting MAPK pathway activation. Therefore, melanomas with NF1-inactivating mutations may respond to MEK1/2 inhibition.

It has been reported that mutation patterns are different among melanoma subtypes (28). Here, we further examined mutant genes in each melanoma subtype in the cohort of “pan-negative” tumor samples. No mutated genes were found in more than 1 melanoma samples originating from mucosal sites, uveal sites, or unknown primary sites. Two genes, KIT (mutated in 3 samples) and PREX2 (2 samples), were mutated in acral melanoma; however, this result did not reach statistical significance (KIT: P = 0.535; PREX2: P = 0.669). Results for sun-exposed melanomas (cutaneous sites) are summarized in Supplementary Table S15. Specifically, all 12 genes in Table 4 were consistently found in cutaneous melanoma. This is not unexpected because the majority of tumors in “pan-negative” samples were derived from cutaneous sites (44/69 = 63.8%). Besides the 12 genes, many other genes were also found to be enriched in pan-negative cutaneous melanoma samples, including MAPK1 (a.k.a. ERK2) and MARK1 (Supplementary Table S15).

With the rapid advances in NGS technologies, the molecular profiles of many tumors have been extensively characterized (9). Knowledge gained from these mutation profiles can potentially be used to identify relevant cellular targets for the development of targeted therapeutics to implement a personalized medicine approach for patients whose tumors harbor driver mutations in these targets (55). In this study, we systematically examined mutations from recently published melanoma NGS data encompassing 241 paired tumor-normal samples, aiming to identify potential, clinically relevant mutations in melanomas that did not present with any of the well-established and well-characterized driver mutations in BRAF, NRAS, KIT, GNAQ, and GNA11, which we describe as “pan-negative” tumors. Notably, and consistent with other reports (2), we confirmed that approximately half of the melanoma tumors (50.2%) harbored BRAF V600 mutations. Through our statistical analyses, we identified several genes with nonsilent mutations that were likely to co-occur with BRAF- and NRAS-driven melanomas. Genes that co-occurred with BRAF-mutant tumors included TP53 and COL1A1, whereas genes that co-occurred with NRAS-mutant tumors included PPP6C, KALRN, PIK3R4, TRPM6, GUCY2C, and PRKAA2. The co-occurring mutations could modify responsiveness to targeted therapies for BRAF- or NRAS-mutant melanoma. Subsequent analysis of the sequence data from 69 pan-negative samples revealed less common mutations in BRAF, NRAS, KIT, GNAQ, and GNA11 along with many mutations in other genes encoding guanine nucleotide binding protein α-subunits. Furthermore, our statistical analysis identified 12 significantly mutated genes in pan-negative samples, including ALK, STK31, DGKI, RAC1, EPHA4, ADAMTS18, EPHA7, ERBB4, TAF1L, NF1, SYK, and KDR. Among them, we discovered that 5 genes (RAC1, ADAMTS18, EPHA7, TAF1L, and NF1) harbored a recurrent mutation that presented in at least 2 pan-negative tumor samples. Interestingly, somatic nonsense mutations were frequently found in 2 genes: 6 nonsense mutations (including one recurrent mutation, W784X) in NF1 and 3 nonsense mutations in STK31 (for details, see Supplementary Table S14). These newly identified genes in pan-negative melanoma samples may serve as potential targets for therapeutic agents in the treatment of melanoma without known BRAF, NRAS, KIT, GNAQ, and GNA11 mutations.

If tumors containing these mutations are activated and inhibition of the mutation(s) leads to cell death or attenuated growth, then clinically targeting them through pharmacologic approaches may be a useful strategy. For example, the kinase inhibitor crizotinib might be tested in tumors with ALK mutations, STK inhibitors could be tested in tumors with STK31 mutations, and MEK1/2 inhibitors may be tested in tumors with NF1 mutations. Moreover, the rare mutations in the 5 known driver genes (BRAF, NRAS, KIT, GNAQ, and GNA11) in this pan-negative melanoma panel clearly demonstrated that the commonly used SNaPshot screen could overlook a subset of patients with melanoma that might also respond to the currently available, effective BRAF, NRAS, KIT, GNAQ, and GNA11 targeting inhibitors, or at least their effector pathways (MEK1/2 inhibitors). For example, non-V600E BRAF mutations, including L597R/S/Q/V and K601E, are actionable in melanoma with MEK1/2 inhibition (5).

The NGS analysis method in Hodis and colleagues (19) requires mutation data in both the exonic regions (both nonsynonymous and synonymous mutations) and intronic regions. Because some studies did not report mutations in the intronic regions or synonymous mutations in the exonic regions, we developed our method using only nonsynonymous mutations from all the available studies. However, Hodis and colleagues applied a filtering strategy to remove those exonic SNVs whose adjacent noncoding regions had high mutation rates. This would likely affect our pool-based association analyses for those mutations that occurred in the hypermutable regions. Because the raw sequence data from Hodis and colleagues is not immediately available, the results reported in this study are not definitive. Other factors, such as sequence coverage and variant calling pipelines, which affect SNV calling, might affect our mutant gene association analyses. Such factors, however, are currently the main challenges for NGS applications (56). Moreover, recent studies suggested that synonymous mutations might alter a protein's expression, three-dimensional conformation, and function (57). A recent WGS study has detected a recurrent functional synonymous mutation in melanoma (58). However, it is still too early to systematically explore the effects of synonymous mutations on melanoma. In additions, we narrowed the candidate gene list based on their association with known melanoma driver genes. Although this method seems effective, it also likely missed other potential candidate genes. In summary, different from each individual NGS study, which attempted to identify mutations and genes involved in melanoma as a whole, our study uniquely provides a roadmap of potentially druggable mutations and genes in specific melanoma subtypes, especially those lacking any known driver mutations (pan-negative samples).

Another limitation of this study is that only SNVs were examined. In future work, we will extend our analysis to other types of mutations, such as insertions and deletions (indels), copy-number variations (CNV), and structural variants, for a more comprehensive map of genetic alterations in melanoma. In addition, more work is needed to link the promising genomic alterations to clinical outcomes or drug activity. We will attempt to identify genetic variants as the predictors of drug sensitivity by mining datasets from the Cancer Cell Line Encyclopedia (CCLE; http://www.broadinstitute.org/ccle/home) project, the Genomics of Drug Sensitivity in Cancer project (GDSC, http://www.cancerrxgene.org), and other large-scale pharmacogenomic datasets in cancer cells.

To our knowledge, this study represents the largest effort to-date that characterizes potentially clinically relevant mutations in melanomas without known BRAF, NRAS, KIT, GNAQ, and GNA11 driver mutations (i.e., pan-negative patients). Our comprehensive analyses of the mutation patterns in 69 pan-negative melanoma genomes, along with 172 melanoma genomes harboring known driver mutations, demonstrated the effectiveness of discovering underrepresented somatic mutations potentially relevant to targeted therapies in patients with melanoma who otherwise could not receive targeted therapies. As clinical NGS rapidly emerges, the high capacity and low cost of these platforms will soon allow routine screening of more comprehensive sets of tumor mutations in CLIA-certified laboratories, and thus, will benefit numerous patients with cancer.

In summary, we used existing published NGS data from 241 tumor-normal pairs to identify mutations that could influence existing and emerging targeted therapies in both driver mutation positive and pan-negative melanomas. Future preclinical studies will need to be conducted to determine whether these mutations are driver or passenger mutations and how they will affect response to targeted therapies specifically in melanoma. The strategies proposed in this study can also be applied to other types of cancer, like lung cancer, that have extensive NGS data available but still have many cases without well-defined, targetable driver mutations.

K.B. Dahlman has received honoraria from the speakers' bureau of Illumina. No potential conflicts of interest were disclosed by the other authors.

The funders had no role in study design, data collection, and analysis, decision to publish, or preparation of the article.

Conception and design: J. Xia, J. Sosman, W. Pao, Z. Zhao

Development of methodology: J. Xia, Z. Zhao

Acquisition of data (provided animals, acquired and managed patients, provided facilities, etc.): J. Xia, P. Jia, K.B. Dahlman, D. Johnson

Analysis and interpretation of data (e.g., statistical analysis, biostatistics, computational analysis): J. Xia, P. Jia, K.E. Hutchinson, K.B. Dahlman, W. Pao, Z. Zhao

Writing, review, and/or revision of the manuscript: J. Xia, P. Jia, K.E. Hutchinson, K.B. Dahlman, J. Sosman, W. Pao, Z. Zhao

Administrative, technical, or material support (i.e., reporting or organizing data, constructing databases): W. Pao

Study supervision: D. Johnson, W. Pao, Z. Zhao

The authors thank the members of the Bioinformatics and Systems Medicine Laboratory for their valuable discussion in this project.

This project was partially supported by NIH grants (R01LM011177 and R03CA167695 to Z. Zhao, K12CA0906525 to D. Johnson, and K24CA097588 to J. Sosman); a Stand Up To Cancer Innovative Research Grant, a Program of the Entertainment Industry Foundation (SU2C-AACR-IRG0109 to W. Pao); and VICC Cancer Center Core Grant (P30CA068485 to J.A. Pietenpol). The authors also acknowledge support from the Joanna M. Nicolay Melanoma Foundation 2013 Research Scholar Award, the TJ Martell Foundation, the Robert J. Kleberg, Jr. and Helen C. Kleberg Foundation, the Ingram Professorship, the American Cancer Society Mary Hendrickson-Johnson Melanoma Professorship, and the Craig Chair Professorship Funds.

The costs of publication of this article were defrayed in part by the payment of page charges. This article must therefore be hereby marked advertisement in accordance with 18 U.S.C. Section 1734 solely to indicate this fact.