The immune-related role of BRAF in melanoma.

BACKGROUND: The existence of a dichotomy between immunologically active and quiescent tumor phenotypes has been recently recognized in several types of cancer. The activation of a Th1 type of immune signature has been shown to confer better prognosis and likelihood to respond to immunotherapy. However, whether such dichotomy depends on the genetic make-up of individual cancers is not known yet. BRAF and NRAS mutations are commonly acquired during melanoma progression. Here we explored the role of BRAF and NRAS mutations in influencing the immune phenotype based on a classification previously identified by our group.
METHODS: One-hundred-thirteen melanoma metastases underwent microarray analysis and BRAF and NRAS genotyping. Allele-specific PCR was also performed in order to exclude low-frequency mutations.
RESULTS: Comparison between BRAF and NRAS mutant versus wild type samples identified mostly constituents or regulators of MAPK and related pathways. When testing gene lists discriminative of BRAF, NRAS and MAPK alterations, we found that 112 BRAF-specific transcripts were able to distinguish the two immune-related phenotypes already described in melanoma, with the poor phenotype associated mostly with BRAF mutation. Noteworthy, such association was stronger in samples displaying low BRAF mRNA expression. However, when testing NRAS mutations, we were not able to find the same association.
CONCLUSION: This study suggests that BRAF mutation-related specific transcripts associate with a poor phenotype in melanoma and provide a nest for further investigation.

Introduction

Malignant cutaneous melanoma is an aggressive neoplasm characterized by a complex etiology that challenges targeted therapies (Gray‐Schopfer et al., 2007; Hocker et al., 2008; Nikolaou et al., 2012; Tomei et al., 2014). Several molecular alterations occur during melanoma progression. The most commonly mutated pathway is the mitogen‐activated protein kinases (MAPK)/ERK cascade (Dhomen and Marais, 2009; Palmieri et al., 2009). The critical role of MAPK/ERK activation in melanoma development has prompted attempts to exploit this pathway as a therapeutic target (Ascierto et al., 2012; Frankel et al., 2003; Ibrahim and Haluska, 2009; Karasarides et al., 2004; Zambon et al., 2012). The activation of the MAPK/ERK signaling occurs either through gain‐of‐function mutations in the BRAF and NRAS genes or through autocrine growth factor stimulation (Curtin et al., 2005; Edlundh‐Rose et al., 2006; Satyamoorthy et al., 2003; Smalley, 2003).

Documented mutations have been found in the kinase domain of BRAF, encoded by exons 11 and 15 (Brose et al., 2002; Davies et al., 2002) with a frequency of 50–70% (Davies et al., 2002). The majority of these mutations affect one critical amino acid, resulting in the V600E substitution that leads to constitutive kinase activation of BRAF (Davies et al., 2002). The V600E mutation, accounting for more than 90% of all BRAF mutations, is associated with an over 400‐fold greater basal activity as compared to wild type BRAF (Davies et al., 2002; Kumar et al., 2004; Pollock et al., 2003).

In some melanomas without BRAF mutation, the MAPK/ERK pathway is constitutively activated through mutation of NRAS (van Elsas et al., 1996). Activating mutations in NRAS, mostly at codon 61, have been identified in about 15% of melanomas (Curtin et al., 2005) and result in the reduction of intrinsic GTPase activity and in the constitutive activation of NRAS.

Both functional and genetic studies indicate that BRAF and NRAS act linearly in the signaling pathway, which is demonstrated by the almost mutually exclusiveness of mutations in these genes and the consequent downstream activation.

Although the oncogenic potential of BRAF and NRAS alterations has been attributed to reduced apoptosis, increased invasiveness and increased metastatic behavior (Gray‐Schopfer et al., 2005; Maurer et al., 2011), the role of BRAF and NRAS in the immunological landscape of cutaneous melanoma has been poorly investigated and the effects of BRAF and NRAS mutations on global gene expression remain to be understood. A few studies have attempted to test whether a BRAF and NRAS mutation‐associated gene expression signature exists in melanoma but none of them evaluated such signatures according to the immunological phenotypes (Bittner et al., 2000; Johansson et al., 2007; Kannengiesser et al., 2008; Pavey et al., 2004; Schramm et al., 2012). Relevant in this context, gene expression signatures have been shown to predict response to anticancer immunotherapy (Ulloa‐Montoya et al., 2013; Wang et al., 2013a).

Emerging data is now revealing the existence of at least two different immune phenotypes in melanoma, a Th17 phenotype associated with the over‐expression of WNT5A, enhanced cellular motility, a more undifferentiated status and poor prognosis, and a Th1 immune phenotype associated with the expression of melanocytic lineage specific transcripts, a higher responsiveness to immune cytokines, a more differentiated status and better prognosis (Bittner et al., 2000; Murtas et al., 2013; Spivey et al., 2012). The existence of similar phenotypes has also been shown in colorectal cancer where type, density and location of immune cells have been found to predict patient survival (Galon et al., 2013).

It is not yet, however, clear whether the existence of these two different phenotypes depends upon the genetic background of the host or is dictated by the genetics of the tumor or both (Wang et al., 2013, 2012). The finding that selective BRAF inhibitors induce marked T cell infiltration, enhance melanoma antigen expression and improve the anti‐tumor activity of adoptive immunotherapy in mice (Tompers Frederick et al., 2013; Viros et al., 2008; Wilmott et al., 2012) emphasizes the need to study this link. Moreover BRAF and NRAS activating mutations have been associated with the production of various immunosuppressive factors in melanoma cell lines (Castelli et al., 1994; Sumimoto et al., 2006). Furthermore, a recent study demonstrates that re‐activation of MAPK signaling provides a mechanism of therapeutic resistance via modulation of host immune responses (Devitt et al., 2011). Most relevant to our work, Khalili et al. recently showed that the BRAF V600E mutation promotes stromal cell‐mediated immunosuppression via induction of IL‐1 (Khalili et al., 2012), highlighting the immunosuppressive role of BRAF in melanoma.

We recently proposed a genetic classification of melanoma metastases based on copy number variation and consistency of genes expressed in vivo and in vitro (Spivey et al., 2012). We found that genes consistently expressed by 15 melanoma cell lines (CMs) and their parental tissues (TMs) were critical for oncogenesis and their respective copy number influenced their expression. Most importantly, these genes were able to categorize melanoma metastases into two divergent phenotypes (TARA class: transcriptional adjustments related to amplification/deletions): one with prevalent expression of cancer testis antigens, enhanced cyclin activity, WNT signaling, and a Th17 immune phenotype (Class A) and the other one with prevalent expression of genes associated with melanoma signaling and with a Th1 immune phenotype (Class B). An intermediate third class (C) was further identified. The 3 phenotypes were confirmed by unsupervised principal component analysis. Here, we tested whether these phenotypes might be at least in part explained by BRAF and NRAS mutations in melanoma.

We found that BRAF mutation‐related specific transcripts associate with the poor melanoma phenotype and such association resulted particularly strong in BRAF mutant tumors displaying low expression levels of this gene.

Our findings strongly suggest a role of BRAF mutation and mRNA expression in influencing the immune phenotypes in melanoma.

Materials and methods

Specimens

One‐hundred‐thirteen pre‐treatment snap frozen tumor biopsies were collected from patients treated at the Surgery Branch, National Cancer Institute, Bethesda, Maryland, who underwent adoptive immunotherapy with tumor infiltrating lymphocytes and myeloablative conditioning. Fifteen melanoma cell lines were derived from as many metastatic melanomas. Identity confirmation of cell lines and parental tissues was carried out by HLA phenotyping as previously described (Spivey et al., 2012).

BRAF and NRAS genotyping

Genomic DNA was isolated using QIAamp DNA Mini Kit (Qiagen, Germantown, MD) and phenol/chloroform method from cell lines and tumor tissues, respectively and according to the standard protocols. DNA quality and quantity was estimated using Nanodrop (ThermoScientific, Pittsburgh, PA).

Each sample was screened for mutations in exons 11 and 15 for BRAF gene and exon 1 and 2 for NRAS gene. PCR was performed in 20 μl final volume, containing 50 ng of genomic DNA, 10 μl of Qiagen HotStarTaq Master Mix Kit (Valencia, CA) and 500 nM forward and reverse primers with the following cycling conditions: initial denaturation at 95 °C for 10 min; 35 cycles at 95 °C for 30 s, 56 °C for 30 s and 72 °C for 30 s; final step 72 °C for 10 min. Primers were selected using Primer3 software (http://frodo.wi.mit.edu/):

BRAF_ex11_F: 5′‐TCCCTCTCAGGCATAAGGTAA‐3′

BRAF_ex11_R: 5′‐CGAACAGTGAATATTTCCTTTGAT‐3′

BRAF_ex15_F: 5′‐TCATAATGCTTGCTCTGATAGGA‐3′

BRAF_ex15_R: 5′‐GGCCAAAAATTTAATCAGTGGA‐3′

NRAS_ex1_F: 5′‐CACTAGGGTTTTCATTTCCATTG‐3′

NRAS_ex1_R: 5′‐TCCTTTAATACAGAATATGGGTAAAGA‐3′

NRAS_ex2_F: 5′‐ATAGCATTGCATTCCCTGTG‐3′

NRAS_ex2_R: 5′‐CACAAAGATCATCCTTTCAGAGA‐3′

In each PCR reaction distilled water was used as a negative control.

PCR products were purified with Exosap‐IT (USB Corporation, Cleveland, OH) and labeled using Big Dye terminator kit v3.1 (Applied Biosystems, Foster City, CA). Excess dye terminators were removed using DyeEx 96 Kit columns following the manufacturer's instructions (Qiagen). Sequencing was then performed using Biosystems 3730 Genetic Analyzer (Foster City, CA) and analyzed by Sequencher software (Genecodes, Ann Arbor, MI).

BRAF and NRAS allele‐specific PCR (AS‐PCR)

AS‐PCR for the detection of BRAF V600E and NRAS Q61R, Q61K, Q61L and Q61H mutations was conducted as described above. Primers were designed to possess two bases substitution at 3′‐end compared to wild type sequences, as previously published (Linard et al., 2002; Sensi et al., 2006). These primers were, respectively:

BRAF15_R: 5′‐GGCCAAAAATTTAATCAGTGGA‐3′

BRAF15_WT_F: 5′‐TAGGTGATTTTGGTCTAGCTACAGT‐3′

BRAF15_V600E_F: 5′‐GGTGATTTTGGTCTAGCTACAAA‐3′

NRAS2_R: 5′‐TGACTTGCTATTATTGATGG‐3′

NRAS2_WT_F: 5′‐CATACTGGATACAGCTGGAC‐3′

NRAS2_Q61K_F: 5′‐CATACTGGATACAGCTGGGA‐3′

NRAS2_Q61R_F: 5′‐ATACTGGATACAGCTGGAAG‐3′

NRAS2_Q61L_F: 5′‐ATACTGGATACAGCTGGATT‐3′

NRAS2_Q61H1_F: 5′‐TACTGGATACAGCTGGACTT‐3′

NRAS2_Q61H2_F: 5′‐TACTGGATACAGCTGGACTC‐3′

Cycling conditions were as following: initial denaturation at 95 °C for 10 min; 40 cycles at 95 °C for 30 s, 54 °C for 30 s and 72 °C for 45 s; final step 72 °C for 10 min. PCR products were run on 2% agarose gel.

Microarray analysis

Total RNA from the 15 cell lines and 113 melanoma metastases was isolated using miRNeasy minikit (Qiagen) according to the manufacture's protocol (Bedognetti et al., 2013; Spivey et al., 2012). RNA quality and quantity was estimated using Nanodrop (Thermo Scientific) and Agilent 2100 Bioanalyzer (Agilent Technologies, Santa Clara, CA, USA). First‐ and second‐strand cDNAs were synthesized from 300 ng of total RNA according to manufacturer's instructions (Ambion WT Expression Kit). cDNAs were fragmented, biotinylated, and hybridized to the GeneChip Human Gene 1.0 ST Arrays (Affymetrix, Santa Clara, CA). The arrays were washed and stained on a GeneChip Fluidics Station 450 (Affymetrix). Scanning was carried out with the GeneChip Scanner 3000 and image analysis with the Affymetrix GeneChip Command Console Scan Control. Expression data were normalized, background‐corrected, and summarized using the robust multichip average (RMA) algorithm, http://www.partek.com/. Data were log2 transformed for subsequent statistical analysis. For the external validation, GSE22155 array data were downloaded from GEO (http://www.ncbi.nlm.nih.gov/geo/) and log2 transformed. Hierarchical clustering and PCA analyses were performed using Partek software (Partek Inc., St. Louis, MO, USA).

BRAF siRNA (small interference) in melanoma cell lines

Melanoma cell lines were assigned to “high” and “low” groups based on the BRAF average expression. The “high” cell lines were further selected for BRAF RNA interference experiment.

Control siRNA (siGENOME Non‐Targeting siRNA Pool #1, D‐001206‐13‐05), BRAF siRNA (SMARTpool siGENOME BRAF siRNA, M‐003460‐03‐0005) and DharmaFECT 1 Transfection Reagent (T‐2001‐01) were purchased from Thermo Fisher Scientific‐Dharmacon (Pittsburgh PA, USA).

Melanoma cell lines were seeded in 6‐well plates at 5 × 105 cells/well and the day after transfected with BRAF or Control siRNA using DharmaFECT 1 Transfection Reagent (final siRNA concentration, 50 nM) following the manufacturer's instructions. Forty‐eight hour after transfection, cells were collected and total RNA was extracted using Trizol (Invitrogen – Life Technologies Carlsbad, CA, USA). One μg of total RNA was reverse‐transcribed in cDNA using SuperScript III (Invitrogen – Life Technologies) according to the manufacturers' directions. Four genes associated to the Th1 and Th17 based on a previous report (Spivey et al., 2012) were selected for gene expression analysis, namely STAT1, GBP1 IL17A and IL17B. Quantitative Real Time‐PCR (qPCR) was performed on an AB ViiA 7 real‐time PCR system (Applied Biosystems – Life technologies, Grand Island, NY, USA) using TaqMan gene expression assays (Applied Biosystems – Life Technologies) with the following probes: Hs00269944_m1 for BRAF, Hs01013996_m1 for STAT1, Hs00977005_m1 for GBP1, Hs00174383_m1 for IL17A and Hs00975262_m1 for IL17B.

Biological samples were run in quadruplicate. The endogenous reference gene 18S was used to normalize each gene expression level and as internal reference (Eukaryotic 18S rRNA Endogenous Control, 4310893E, Applied Biosystems – Life Technologies). For each qPCR experiment, a no‐template reaction was included as negative control. The thermal cycling conditions were as follows: initial denaturing at 95 °C for 10 min, 40 cycles of 95 °C for 15 s and 60 °C for 1 min. For each gene, the fold difference between control siRNA‐ and BRAF siRNA‐transfected melanoma cell lines was calculated using the comparative 2–ΔΔCt method.

Statistical analyses

To identify transcripts whose expression was associated with BRAF and NRAS mutational status we performed class comparison of BRAF mutated versus BRAF and NRAS wild type samples and NRAS mutated versus BRAF and NRAS wild type samples. Class comparison was based on analysis of variance (ANOVA). Principal component analysis (PCA) was applied for visualization when relevant. Heat maps are presented based on Partek visualization program. Fisher's exact and χ 2 tests were used as appropriate. All statistical tests were two‐sided. We used relatively low‐stringency class comparison for subsequent enrichment analysis (Bedognetti et al., 2011). Functional gene network analysis was performed using the Ingenuity Pathway Analysis system (IPA) tools 3.0 which transforms large data sets into a group of relevant networks containing direct and indirect relationships between genes based on known interactions in the literature (http://www.ingenuity.com, Ingenuity System Inc., Redwood City, CA, USA).

Results

BRAF and NRAS mutations in cell lines and tumors

Eight out of fifteen (53%) cell lines and sixty‐seven out of one‐hundred‐thirteen (59%) tumors harbored a BRAF mutation. Most of the mutations (8/8 in cell lines and 64/67 in tumors) were the classical V600E substitution. Three samples carried K601E, V600K and V600R mutations, respectively. Among the seven cell lines and forty‐six tumors that did not display a BRAF mutation, six cell lines and thirteen tumors harbored activating mutation in NRAS gene (four cell lines displayed Q61R, two cell lines displayed Q61K and Q61L respectively; among tumors, NRAS mutations were: eight Q61R, two Q61K, and three Q61L, Q61H and G13R respectively). BRAF and NRAS mutations were mutually exclusive both in cell lines and tumors.

Interestingly, when comparing BRAF and NRAS mutational status in cell lines and matched parental tumors, three pairs were discordant.

AS‐PCR (allele‐specific PCR) revealed low frequency BRAF and NRAS mutations

To exclude the possibility that BRAF and NRAS mutations could exist at low frequency in melanoma, we probed the samples with a highly sensitive allele‐specific PCR (Supplementary Figure 1).

When comparing the BRAF and NRAS mutational status of cell lines and parental tumors after applying AS‐PCR, we interestingly found that the two tumors whose corresponding cell line was BRAF V600E actually resulted mutated by AS‐PCR (Table 1). However, the NRAS Q61R tumor whose corresponding cell line was wild type remained wild type after applying AS‐PCR, suggesting that in this specific case the mutation might have been acquired during cell culturing.

Table 1

BRAF and NRAS mutational status of cell lines and matched tumors.

Tumor	BRAF	NRAS	Cell line	BRAF	NRAS
SAR‐32	V600E	wt	2492	V600E	wt
SAR‐38	V600E	wt	2448	V600E	wt
SAR‐39	V600E	wt	3104	V600E	wt
SAR‐58	V600E	wt	2523	V600E	wt
SAR‐59	V600E	wt	2224	V600E	wt
SAR‐89	V600E	wt	2035	V600E	wt
SAR‐52	wt	Q61K	2075	wt	Q61K
SAR‐77	wt	Q61L	3107	wt	Q61L
SAR‐17	wt	Q61R	2744	wt	Q61R
SAR‐33	wt	Q61R	2155	wt	Q61R
SAR‐102	wt	Q61R	1866	wt	Q61R
SAR‐121	wt	wt	2805	wt	wt
SAR‐63	wt	wt	2458	V600E	wt
SAR‐78	wt	wt	3025	V600E	wt
SAR‐21	wt	wt	2427	wt	Q61R

wt: wild type.

We next tested the presence of BRAF and NRAS low frequency mutations in the wild type tumors as well by applying AS‐PCR. Among the 33 BRAF and NRAS wild type samples, 12 turned out to be BRAF V600E and 1 NRAS Q61R‐mutated. We thus considered these samples as low frequency‐mutants and excluded them from further expression data analyses. However, we decided to keep them for the graphical analyses.

BRAF and NRAS mutation‐associated gene expression in tumors

To test whether BRAF and NRAS mutations affect the transcriptome of the 113 metastases, we applied PCA to the complete data set. The assignment of the individual metastases to BRAF, NRAS and wild type groups did not predict their distribution in three‐dimensional space suggesting that BRAF and NRAS mutations do not directly affect gene expression at global level (Figure 1a).

Figure 1

Principal component analysis of BRAF and NRAS mutations based on the complete gene dataset (a). Clustering of melanoma metastases based on genes discriminative of BRAF (b), NRAS (c) and MAPK (d) status at p < 0.01.

Gene expression data were next used to identify genes discriminating the 113 tumors between mutated and wild type for both BRAF and NRAS. At a significance level of p < 0.01, 583 transcripts were differentially expressed between BRAF mutant versus wild type samples (Figure 1b) and 186 between NRAS mutant versus wild type samples (Figure 1c). Notably, when assessing genes concordantly deregulated in both BRAF and NRAS mutant samples, it was found that most encoded constituents or regulators of the MAPK/ERK and related pathways.

In accordance with the findings of other groups (Bloethner et al., 2005; Khalili et al., 2012; Packer et al., 2009), the top ranking genes with a potential role, direct or indirect, in the MAPK/ERK pathway included the dual specificity phosphatase genes, DUSP6 and DUSP4; an inhibitor of MAPK signaling, SPRY2; a serine protease inhibitor, SERPINE2; and genes of the 14‐3‐3 family.

BRAF mutation associates with the poor phenotype in melanoma

We next tested whether genes discriminative of BRAF and NRAS mutational status were able to distinguish different immunological phenotypes. Towards this goal, we used a genetic classification which has been identified by our group in a previous study (Spivey et al., 2012) and which we refer to as TARA (“transcriptional adjustments related to amplification/deletions”). Five TARA classes of melanoma metastases were defined according to the clustering based on genes concordantly expressed in vivo and in vitro and whose expression was also predicted by the corresponding copy number (Supplementary Figure 2, Supplementary Table 1). Class A and Aa referred to an unfavorable immune phenotype associated with a Th17 activation and bad prognosis while class B and Bb referred to a more favorable phenotype associated with a Th1 activation and better prognosis. A third category sitting astride the two polar groups was also identified (Class C). It should be stated that the classification used in this paper differs from the already published one in the number of TARA classes defined. Contrary to our previous study (Spivey et al., 2012), we used 5 instead of 3 TARA classes. As the extreme classes “A” and “B” previously identified showed a gradient in the expression of the genes whose transcriptional efficiency was predicted by the corresponding copy number, we thought that it would have been more appropriate stratifying all the sample set in 5 TARA classes according to the self‐organizing map based on the top 500 genes obtained when comparing TARA A versus TARA B, as shown in the Supplementary Figure 2.

When we used the genes discriminative of BRAF and NRAS status to cluster the melanoma metastases, we observed that BRAF‐discriminative genes clearly segregate samples according to TARA classification. Of note, TARA class A samples laid at the extreme opposite to the non‐TARA class A samples. Although this association did not reach statistical significance (A + Aa versus all: χ 2 2.2, p = 0.1), overall these data suggest a link between BRAF mutational status and TARA classification. No association between NRAS mutational status and TARA classification was observed.

We next defined as “MAPK mutant” those samples carrying either BRAF or NRAS mutation, and performed class comparison between MAPK mutant and wild type samples. Five‐hundred‐ninety‐three genes were differentially expressed between the two groups at a significance level of p < 0.01. When using these genes to cluster the 113 melanoma metastases, we found a tendency of the mutant samples to be associated with TARA class A and Aa, though not absolutely. Interestingly, these genes resulted once again to encode constituents or regulators of the MAPK/ERK pathway.

We then sought to assess whether the association of BRAF mutation with TARA classification was restricted to BRAF‐specific genes or whether it was dependent upon transcripts associated with MAPK activation. Thus, we combined the lists of genes obtained by the 3 class comparisons and we found: 52 genes significantly deregulated in all the 3 gene lists (BRAF, NRAS, MAPK), 107 genes NRAS‐specific, 95 genes MAPK‐specific and 112 genes BRAF‐specific. These sub‐groups of genes were further used to cluster the 113 melanoma metastases (Figure 2).

Figure 2

Venn diagram and cluster of melanoma metastases (a) based on 95 MAPK‐restricted transcripts (green), 52 MAPK‐specific transcripts (red) and 107 NRAS‐specific transcripts (blue). Venn diagram and self‐organizing heat map of 67 BRAF mutant and 20 BRAF wild type samples based on 112 BRAF‐specific transcripts (b, the heat map top legend refers to BRAF mutational status). Self‐organizing heat map of the 113 metastases based on 112 BRAF‐specific transcripts (c, the top legend refers to MAPK mutational status before and after AS‐PCR, and TARA classification). Functional interpretation analysis of the 112 BRAF‐specific transcripts (d, green: down regulated in BRAF wild type samples, red: up regulated in BRAF wild type samples).

Quite interestingly, the self‐organizing heat map showed that the expression of BRAF‐specific genes was highly predictive of the TARA classification (Figure 2b). Functional interpretation of the 112 BRAF‐specific transcripts ranked “Role of JAK1, JAK2 and TYK2 in Interferon Signaling” as the top canonical pathway (up regulated in wild type samples as compared to the mutant samples, Figure 2d) supporting the link between BRAF mutation and immune phenotype. The hierarchical cluster of melanoma metastases using the 112 BRAF specific genes according to the BRAF status (excluding NRAS mutants and samples uncertain after applying AS‐PCR) was able to clearly discriminate BRAF mutants versus WT samples (Figure 2b, top panel).

We could not find any association between NRAS/MAPK mutational status and TARA class when clustering the 113 metastases according to NRAS‐ and MAPK‐specific genes, respectively (Figure 2a), suggesting that the transcriptional consequences resulting from mutation of BRAF and NRAS are different, presumably through their differential capacity to receive and transduce input signals through various effectors. Moreover, some of the genes that discriminated BRAF or NRAS mutant samples were independent of the MAPK activation (112 and 107 genes were BRAF‐ and NRAS‐specific, respectively). This finding implies that some of the genes in the BRAF discriminating gene list may not be necessarily the direct targets of the transcription factors that are ultimately activated by MAPKs. This hypothesis suggests that novel therapeutic targets outside the MAPK pathway could be used to treat melanoma carrying BRAF mutations.

In order to validate our findings, we used GSE22155 as an independent dataset. GSE22155 included data from 57 melanoma metastases (28 BRAF mutated, 12 NRAS mutated and 17 wild type) defined as “high‐immune”, “proliferative”, “pigmentation” and “normal‐like” based on a molecular classification (Jonsson et al., 2010). First, we determined whether this classification was comparable to TARA classification. To do so, we selected genes representative of both TARA A/Th17 (IL17F, IL17B, WNT5A, IL17A) and TARA B/Th1 phenotype (TYR, CCL5, CXCL10, GBP1, STAT1) (Spivey et al., 2012). The unsupervised clustering clearly showed the high immune samples to behave similarly to TARA B/Th1 samples (expressing TYR, CCL5, CXCL10, GBP1, STAT1) and the proliferative samples to behave similarly to TARA A/Th17 samples (expressing IL17F, IL17B, WNT5A, IL17A, Figure 3a). As an additional proof, we also plotted the GSE22155 samples in a PCA graph based on the top 500 genes which were used to define TARA classes (Figure 3b) and we found that such genes divided clearly the high immune versus proliferative groups, allowing us to consider these samples comparable to TARA B and TARA A, respectively. When using the 112 BRAF specific transcripts the high immune samples clearly distributed separately compared to the proliferative ones, suggesting that the BRAF‐specific transcripts are indeed related to an immune phenotype (Figure 3c).

Figure 3

Self‐organizing map of GSE22155 dataset based on genes discriminative of TARA class A/Th1 phenotype and TARA class B/Th17 phenotype (a). PCA of GSE22155 dataset based on the 500 TARA specific transcripts (b). PCA and clustering of GSE22155 dataset based on the 112 BRAF specific genes (c).

The association between BRAF mutation and the poor phenotype is stronger in samples displaying low BRAF mRNA expression

We hypothesized that samples belonging to TARA class A and Aa may have a stronger activation of the BRAF downstream pathway through a higher BRAF mRNA expression. To test this hypothesis, we divided the mutant metastases into “high” and “low” according to the average BRAF mRNA expression among all mutant samples. Surprisingly and in discordance with our hypothesis, we observed that, overall, mutant TARA A and Aa samples displayed low BRAF expression while TARA B and Bb samples displayed high BRAF expression (A + Aa versus B + Bb, Fisher's exact test, p = 0.0009). Such an association was absolute when considering the only TARA A group, in other words all the A samples displayed low BRAF expression (Figure 4a). We also tested whether this association existed in wild type samples (Supplementary Figure 3); even though we found a tendency of low expressing mRNA BRAF samples to be TARA A and Aa and high expressing mRNA BRAF samples to be TARA B and Bb, such association did not reach statistical significance (A + Aa versus B + Bb, Fisher's exact test, p = 0.09), suggesting that the association between TARA classification and BRAF mRNA expression is stronger in BRAF mutant samples.

Figure 4

Cluster of melanoma metastases based on 112 BRAF‐specific transcripts according to BRAF mutational status, TARA classification and BRAF mRNA expression (a). Class comparisons between BRAF mutant HIGH versus wild type (WT, b), BRAF mutant LOW versus WT (c), BRAF mutant HIGH versus LOW (d) at p < 0.01; “HIGH” and “LOW” refer to BRAF mRNA expression. Functional interpretation analysis of 6296 genes discriminative of HIGH and LOW BRAF mRNA expression (e; green: down regulated in high expressing BRAF samples, red: up regulated in high expressing BRAF samples).

In order to obtain insights about this intriguing finding, we compared mutant samples expressing high and low levels of BRAF mRNA versus the wild type samples. Class comparison between BRAF mutant samples displaying high BRAF mRNA expression and wild type samples identified only 418 transcripts significantly deregulated at p < 0.01. When comparing samples displaying low BRAF mRNA expression levels and wild type samples, 2252 genes were significantly differentially expressed at p < 0.01, suggesting that, overall, BRAF high expressing samples are more similar to the wild type than the BRAF low expressing samples. Finally, when comparing low and high BRAF expressing mutant samples, we found 6296 transcripts differentially expressed at p < 0.01 (Figure 4). Functional interpretation analysis of these 6296 transcripts revealed several cancer‐related pathways as well as two immune‐related pathways, namely, IL‐2 and Jak/Stat pathways (upregulated in high expressing BRAF samples; Figure 4e).

BRAF small interference (siBRAF) experiment on 8 melanoma cell lines expressing high BRAF mRNA (higher than the average BRAF mRNA expression) was performed to check the expression of two genes previously associated to TARA class B/Th1 phenotype, namely STAT1 and GBP1 (Spivey et al., 2012). As expected, we could not find a reproducible trend in the 8 “high” cell lines when BRAF was silenced (Supplementary Figure 4), suggesting that the tumor microenvironment plays an important role and cell lines do not account for the complex interplay between tumor cells and tumor microenvironment. Of note, 3 cell lines (3025, 2035 and 3104) had very low or absent GBP1 expression. Similarly, the cell line 2035 had very low STAT1 expression; these findings may be explained by the fact that the expression of STAT1 and GBP1 may be driven by the immune or stromal cells rather than by tumor cells (Galon et al., 2013; Jin et al., 2014). It should be noted that we also assessed IL17A and IL17B which are known to be expressed by Th17 lymphocytes; expectably, we could not find any expression in the cell lines tested (data not shown).

Discussion

In the recent years much has been learned about the biology of BRAF and NRAS mutations in melanoma. However, what is not clear yet is whether and how BRAF and NRAS mutations affect the immune context of tumors. The availability of a large cohort of melanoma metastases allowed us to explore this point.

The sequencing of 113 melanoma metastases showed that 59% and 12% were mutated in BRAF and NRAS genes, respectively. These frequencies are in agreement with the reported mutation rate in melanoma. No sample with BRAF mutation also carried an NRAS mutation, which is in keeping with previous reports that have found BRAF and NRAS mutations to be almost mutually exclusive in a variety of cancer types (Brose et al., 2002; Cohen et al., 2002; Davies et al., 2002; Kimura et al., 2003). This suggests that BRAF and NRAS double mutants not only may not provide an advantage for tumor growth but they may even be selected against during tumorigenesis.

The genotyping of 15 paired cell lines and parental tumors demonstrated discrepancies in BRAF/NRAS mutation patterns in 3 pairs, specifically 2 cell lines were BRAF V600E and one NRAS Q61R while their corresponding tumors were wild type, strongly suggesting that, in a fraction of cases, mutations in these two genes might be acquired during progression. In order to exclude the possibility of low frequency mutations in tumors we applied a highly sensitive allele‐specific PCR. Using this technique, the two tumors, assigned previously as wild type by Sanger sequencing and whose corresponding cell line was V600E, revealed in fact a BRAF mutant phenotype. In contrast, the tumor from the tumor‐cell line pair discrepant for NRAS remained wild type after performing AS‐PCR, suggesting that in this specific case the mutation in the NRAS gene might have been acquired and become prevalent during culture expansion. On the other hand, discrepancy in BRAF and NRAS mutation patterns have also been shown between primary tumors and metastases (Colombino et al., 2012), providing additional evidence that mutations can be acquired during tumor progression.

When assessing genes concordantly deregulated in BRAF and NRAS mutant samples, we found that many of them encoded constituents or regulators of the MAPK/ERK and related pathways. In concordance with the findings of other groups (Bloethner et al., 2005; Kannengiesser et al., 2008), the set of genes prominently up regulated in mutant BRAF and NRAS samples with a potential role, direct or indirect, in the MAPK/ERK pathway included the dual specificity phosphatase genes, DUSP4 and DUSP6; a serine protease inhibitor, SERPINE2; an inhibitor of MAPK signaling, SPRY2 and 14‐3‐3 binding proteins.

Initially, we postulated that there might be a common MAPK activation signature (resulting from either BRAF or NRAS mutation); however, we found no overabundance of discriminatory genes for the combined group of samples displaying either BRAF or NRAS mutations. This suggests that the transcriptional consequences resulting from mutations of BRAF or NRAS might be different, although there was overlapping of some genes, presumably due to their differential capacity to receive input signals and transduce them through different effectors. This is in concordance with recent clinical and preclinical studies which suggested that NRAS mutant melanoma are characterized by patterns of signal transduction and biological behavior different from BRAF mutant melanomas (Devitt et al., 2011; Fedorenko et al., 2012; Jakob et al., 2012; Viros et al., 2008).

We further tested whether genes subsets specific for BRAF, NRAS and MAPK activation were able to explain the different phenotypes of melanoma metastases. Towards this goal, we applied these specific genes to the melanoma metastases and we surprisingly found that the 112 BRAF‐specific genes were highly predictive of TARA classification especially in regards to TARA A samples.

We were not able to find the same association when testing NRAS‐ and MAPK‐specific genes, suggesting that the only BRAF‐specific genes are able to explain the different biological behaviors of the metastases.

The most crucial aspect of this study is, however, the finding that the association between BRAF mutation and TARA class A is particularly strong in cases in which BRAF mutation is present in tumors displaying low expression levels of this gene. It may be likely that the mutation affects the half‐life of the protein, however, we cannot exclude that the lower BRAF mRNA expression in the mutant samples is due to a negative feedback loop established by the presence of BRAF mutation.

Class comparison between BRAF mutant samples displaying high levels of BRAF mRNA expression and wild type samples revealed that there is much more similarity between the two than with the low expressing BRAF mutant samples. To gain further insights concerning this surprising finding we performed functional interpretation analysis of 6296 genes differentially expressed between BRAF mutants with high and low expression of the same gene, deregulated pathways included IL‐2 and JAK/Stat signaling pathways, supporting the immunoregulatory role of BRAF. Other reports have also shown an immune‐related role of BRAF mutation as essential for cancer‐immune evasion and able to promote stromal cell‐mediated immunosuppression via induction of IL‐1 (Khalili et al., 2012). A very recent investigation showed that the treatment of metastatic melanoma patients with either BRAF inhibitor alone or BRAF plus MEK inhibitor was associated to a higher expression of melanoma antigens and to an overall more favorable tumor microenvironment, further supporting the implementation of combined BRAF targeted therapy and immunotherapy (Frederick et al., 2013).

Studies thoroughly investigating the association between BRAF/NRAS status and clinical response to immune‐based therapies are lacking. Joseph and colleagues (Joseph et al., 2012) performed a retrospective analysis assessing the correlation between clinical outcome and BRAF/NRAS mutation in melanoma patients treated with IL‐2. Patients with NRAS mutations were more likely to respond to IL‐2 than wild type patients (p = 0.04). However, the advantage in term of overall survival and progression‐free survival did not reach the statistically significance level (p = 0.30 and 0.13, respectively). Although no significant association between BRAF and outcome parameters was detected, the study was underpowered to detect difference between BRAF and wild type patients. Therefore, additional investigation assessing the predictive role of BRAF and NRAS mutation in the context of immunotherapy (e.g. adoptive therapy, anti‐CTLA4 and anti‐PD1/PD‐1 ligands) are warranted.

Conclusion

In conclusion we provide novel insights into the effect of BRAF and NRAS mutations on gene expression according to the immune classification. However, gene expression analysis alone cannot provide a comprehensive molecular understanding of the mechanisms underlying such an association, thus further deeper analyses are warranted to elucidate the mechanism underlying the association of BRAF status with immunological phenotypes (Tomei et al., 2014).

The following are the supplementary data related to this article:

Supplementary data

Click here for additional data file.

Supplementary data

Click here for additional data file.