Melanoma proteomics suggests functional differences related to mutational status

Melanoma is the most lethal cutaneous cancer. New drugs have recently appeared; however, not all patients obtain a benefit of these new drugs. For this reason, it is still necessary to characterize melanoma at molecular level. The aim of this study was to explore the molecular differences between melanoma tumor subtypes, based on BRAF and NRAS mutational status. Fourteen formalin-fixed, paraffin-embedded melanoma samples were analyzed using a high-throughput proteomics approach, combined with probabilistic graphical models and Flux Balance Analysis, to characterize these differences. Proteomics analyses showed differences in expression of proteins related with fatty acid metabolism, melanogenesis and extracellular space between BRAF mutated and BRAF non-mutated melanoma tumors. Additionally, probabilistic graphical models showed differences between melanoma subgroups at biological processes such as melanogenesis or metabolism. On the other hand, Flux Balance Analysis predicts a higher tumor growth rate in BRAF mutated melanoma samples. In conclusion, differential biological processes between melanomas showing a specific mutational status can be detected using combined proteomics and computational approaches.

Melanoma is the most lethal cutaneous cancer, with over 11,000-15,000 estimated deaths in the United States and Europe every year 1,2 . Better understanding of the molecular biology of this tumor has allowed the development of new effective drugs for the treatment of advanced disease, both in the fields of targeted therapies and immunotherapy 3 . However, as not all patients obtain a benefit from new drugs, further insight into the biology of melanoma is needed.
Gene signatures, genomic hybridization, whole-exome genome sequencing, microRNA analysis and other techniques have widely addressed the genomic landscape of melanoma, contributing to significant advances 4,5 . Given the heterogeneity of melanoma and the complex interaction of this tumor with the immune system, the need for combination of biomarkers assays has been recently proposed to properly analyze the disease 6 .
Proteins determine cell phenotype, so proteomics analyses offer the possibility to measure the biological effects caused by genomic abnormalities 7 . The most adapted technique for massive quantification of proteins is mass spectrometry. The recent technological advances in the field allow the identification and quantification of thousands of proteins per sample. Therefore, proteomics offers complementary information to that provided by standard pathology and genomics. We recently demonstrated the feasibility of high-throughput label-free quantitative proteomics to analyze breast cancer from paraffin-embedded samples 8 . In the present study we sought to determine whether high-throughput proteomics combined with computational approaches, such as probabilistic Differential protein expression patterns between subtypes. A Significance Analysis of Microarrays (SAM) was done to find differences among samples at the protein level. Seventeen proteins were found differentially expressed between BRAF mutated and BRAF wild type tumors, all of them underexpressed in BRAF-mutated tumors (Fig. 1, Supplementary Table S2). These proteins are mainly related with fatty acid metabolism.
In addition, delta values between BRAF-mutated and BRAF-wild type, and NRAS-mutated and NRAS-wild type tumors were calculated (Supplementary Table S3). Delta values higher than 1.5 or lower than −1.5 were used to perform gene ontology analyses as well. Proteins related with keratinization, epidermis development and cytoskeleton were underexpressed, whereas proteins involved in melanogenesis and extracellular space were overexpressed in BRAF-mutant as compared with BRAF-wild type samples. SAM and delta analyses did not find significant differences between NRAS-mutant and NRAS-wild type tumors. probabilistic graphical model and functional node activity measurements. A probabilistic graphical model (PGM) was assembled using proteomics data with any a priori information. The resulting network was handled to build a functional structure, as described in previous works [9][10][11] . The network was split into thirteen branches, and gene ontology analyses were performed to establish functional structure. Finally, twelve principal functions were assigned to different branches (from now on, functional nodes) and there was a branch to which no main function could be assigned (Fig. 2, Supplementary Table S4).
Functional node activity measurements were calculated for each functional node using proteins related with the main assigned function and a comparison between BRAF-mutant, NRAS-mutant and double-negative groups was performed. Although the limited number of samples did not allow seeing significant differences, some trends in functional activities were found. For instance, NRAS-mutant had a lower melanosome functional node activity www.nature.com/scientificreports www.nature.com/scientificreports/ than BRAF-mutant or double negative tumors. On the other hand, BRAF-mutant tumors had a higher metabolism functional node activity than NRAS-mutant or double negative (Fig. 3).

Flux balance analysis.
Flux Balance Analysis is a computational approach to assess biochemical networks through the calculation of the flow of metabolites through this network. FBA can be used to calculate the growth rate of an organism or a tumor or the rate of generation of a metabolite. Our model suggest that BRAF mutated tumors may have a higher tumor growth rate than the two other subtypes (Fig. 4, Supplementary Table S5).

Discussion
In this study, proteomics coupled with probabilistic graphical models and flux balance analysis were used to describe differences between melanoma biomarker subgroups in melanoma samples.
Proteomics analyses from FFPE samples are complicated because fixation process provokes chemical alterations. Actually, the number of works in this field in melanoma is scarce. Rezaul et al. analysed by MS one microdissected FFPE melanoma sample and identify 935 proteins, confirming the presence of some of them by immunohistochemistry 12 . A previous work comparing two different FFPE samples of melanocytic nevus and advanced melanoma identified two proteins differentially expressed between those type of tumors 13 . Moreover, Byrum et al. compared primary melanoma, advanced melanoma and benign nevi and identified 171 varying proteins between these diseases 14 . In this work, we demonstrated that it is also possible, using proteomics and FFPE samples, to see more subtle differences such as differential proteins between mutational-status groups.
Mass-spectrometry workflow allowed the detection of 1,606 proteins with detectable expression in at least 75% of the samples and two unique peptides. Differences in fatty acid metabolism, cytoskeleton or keratinization were observed between BRAF-mutant and BRAF-wild type tumors. Also, it seems that differences in functions such as melanogenesis or metabolism may be existed between subgroups.
SAM and gene ontology analysis found 17 proteins differentially expressed between BRAF-mutant and the two other subgroups (NRAS-mutant and double-negative), all of them were underexpressed in BRAF-mutated tumors. These proteins are mainly involved in fatty acid metabolism: acyl-Co A dehydrogenases P11310 (acyl-CoA dehydrogenase medium chain, ACADM), P42765 (acetyl-CoA acyltransferase 2, ACAA2) and P49748 (acyl-CoA dehydrogenase very long chain, ACADVL). It has been previously described that melanoma cells reduce fatty acids from glutamine through tricarboxylic acid cycle 15 . On the other hand, some of the proteins underexpressed in BRAF-mutated tumors have antiproliferative functions. For instance, Q9Y3Z3 (histidine/aspartate (HD)domain containing protein 1, SAMHD1) is implicated in regulation of DNA replication and damage repair and it is proposed to have antiproliferative and tumor suppressive functions in many cancers 16 . Q96CX2 (potassium channel tetradimerization domain containing 12, KCTD12) inhibits proliferation in uveal melanoma cells 17 . O14745 (SLC9A3R1) is involved in suppressing breast cancer cells proliferation 18 . Other proteins of those 17 were previously related with melanoma or melanogenesis processes. For example, P31040 (succinate dehydrogenase complex flavoprotein subunit A,SDHA), which encodes a major catalytic subunit of succinate-ubiquinone reductase, a complex of the mitochondria chain, it was previously related with melanogenesis process 19 . Another protein differentially expressed is P00488 (coagulation factor XIII, F13A1) which it was previously associated with chemotherapy response in melanoma tumors 20 . P08133 (annexin A6, ANXA6) acts as a tumor suppressor in skin cancer and it is involved in in the conversion of melanocytes to malignant melanomas 21 . Lastly, it was previously described that metastatic melanoma tumors have a decreased expression of signal transducer and activator of www.nature.com/scientificreports www.nature.com/scientificreports/ transcription P42224 (signal transducer and activator of transcription 1, STAT1) and it could be one of the mechanism by which melanoma can evade immune detection 22 . Finally, P04899 (G protein subunit alpha i2, GNAI2) contributes to melanoma cell growth 23 . O60749 (Sorting nexin 2, SNX2) is involved in membrane trafficking of growth factor receptors including epidermal growth factor receptor and c-Met 24 . P05107 (integrin subunit beta 2, ITGB2) participates in cell adhesion as well as cell-surface mediated signalling and it is correlated with survival in other cancers such as renal or colorectal tumors 25,26 . As far we know, P14317 (hematopoietic cell-specific Lyn substrate, HCLS1), P63027 (vesicle-associated membrane protein 2, VAMP2), P16402 (histone cluster 1H1 family member d, HIST1H1D) and Q9H0W9 (chromosome 11 open reading frame 54, C11orf54) were not previously related with melanoma or other cancers.
Differential analyses did not show differences between NRAS-mutant and NRAS-wild type tumors, which are attributable to the small sample size. The present study was limited in this regard because it was designed just as a proof of principle that high-throughput proteomics can be used to study clinical samples of melanoma. Future studies with larger sample size will be needed to establish significant differences among subtypes. Interestingly, it seems that delta analyses and SAM provide complementary information about different protein expression patterns, because differential proteins provided by these two analyses were different and they were also related to different biological processes.  www.nature.com/scientificreports www.nature.com/scientificreports/ Besides that, a PGM was used to generate a network based on protein expression data. It is remarkable that, despite the low number of samples, the PGM clearly showed a functional structure. This type of analysis previously demonstrated its utility to characterize other tumor types such as bladder carcinoma or breast cancer and may complement the information provided by genomics 10,11 . Although there are no significant differences in functional node activities between groups, the PGM provides some biological processes as candidates to be deregulated among subtypes. Despite the fact that there are a node without an assigned function, it is remarkable that this node contains proteins related with lipoprotein metabolism and melanine synthesis, both previously related with melanoma progression 27,28 . What is more, despite the reduced number of samples, the high growth rate in BRAF-mutant tumors suggested by FBA, although it is not statistically significant, agrees with previous observations, as BRAF-mutated tumors are more proliferative 29 .
Our study demonstrates that proteomics and computational methods can be applied to the study of formalin-fixed, paraffin-embedded melanoma samples, suggesting that melanoma subgroups, defined by mutational status, present molecular differences. Despite the reduced number of samples analyzed, the probabilistic graphical model showed a functional structure and allowed characterizing differences at biological processes regarding melanoma mutational status. Additionally, flux balance analysis was capable to predict differences at tumor growth rate between these groups. In conclusion, this proof-of-principle work demonstrate the usefulness of proteomics and computational approaches in the molecular characterization of melanoma and suggest some proteins and biological processes that could be used as therapeutic targets after proper validation in larger cohorts. Mass-spectrometry analysis, protein identification and quantification. Proteins were isolated from FFPE samples as previously described 30 . Peptides were desalted using C18 stage tips, dried and resolubilized with 3% acetonitrile, 0.1% formic acid. Samples were analyzed on a QExactive mass spectrometer coupled to a nano EasyLC 1000 (Thermo Fisher Scientific) as previously described 31 Briefly, 3 μL of each sample were loaded. Samples were acquired in a randomized order. Elution was performed at a flow rate of 300 nL/min. Mass spectrometer was operated in data-dependent mode (DDA), acquiring a full-scan MS spectra followed by HCD (higher-energy collision dissociation) fragmentation on the twelve most intense signals per cycle. The samples were acquired using internal lock mass calibration on m/z 371.1010 and 445.1200.

Methods
The acquired raw MS data were processed by MaxQuant (version 1.5.2.8), protein identification was done using the integrated Andromeda search engine as previously described 31 . Briefly, spectra were searched against a forward Swiss Prot-human database, concatenated to a reversed decoyed fasta database (NCBI taxonomy ID9606, release date 2014-05-06). Carbamidomethylation of cysteine, methionine oxidation and N-terminal protein acetylation were set as modifications. Enzyme specificity was set to trypsin/P. The maximum false discovery rate (FDR) was set to 0.01 for peptides and 0.05 for proteins. Intensity, defined as the sum of the precursor intensities of all identified peptides for the respective protein group, was used for protein abundance calculation. All the mass spectrometry raw data files acquired in this study can be downloaded from Chorus (http://chorusproject.org) under the project name "Melanoma proteomics in FFPE samples".
Following MS workflow, identified protein groups were filtered by the presence of at least two unique peptides. Also, detectable expression in at least 75% of the total of samples was required for a protein to be included in subsequent analyses. Additionally, batch effects were removed using limma package 32 and R v 3.2.5 33 .
Protein differential expression analyses. Significance Analysis of Microarrays (SAM) was performed using MeV to find significant differences in protein expression among samples 34 . SAM consists on a t-test corrected by permutations over the number of samples. The significance is determined using a delta parameter, based on the false discovery rate (FDR) 35 . In this analysis, the threshold was a delta = 1.6, equivalent to a FDR = 0.61 (90 th percentile). Protein expression patterns were also compared calculating delta expression values by subtracting the protein expression value for each biomarker status against the expression value for the rest of the tumor samples. Proteins showing a change in their delta expression higher than 1.5 or lower than −1.5 were selected. probabilistic graphical model and activity measurements. A PGM was constructed using R v 3.2.5 33 and grapHD package 36 . The correlation coefficient was selected as the associative method [9][10][11] . PGM are undirected acyclic graphs based on obtaining the spanning tree that maximizes the likelihood and then favoring the graph with the simplest structure preserving the decomposability of tree and minimizing the Bayesian Information Criterion (BIC) as well 37 . Protein expression data was used to build the network with any a priori information. This type of networks reflects relationships between expression patterns, i.e. proteins associated in the network by an edge have related expression patterns. The network was divided into several branches or functional nodes. Network functional structure was explored by Gene Ontology analysis. A major biological function was assigned to each branch. Activity measurements were then calculated by the mean expression of all the proteins in the functional node related to the assigned function. Gene Ontology Analyses were performed in DAVID webtool. "Homo sapiens" was set as background and only Biocarta, GOTERM-FAT and KEGG categories were selected.