Clinical Proteomics of Metastatic Melanoma Reveals Profiles of Organ Specificity and Treatment Resistance

Melanoma comprises only 5% of skin cancers but causes the majority of skin cancer–related deaths (1, 2), primarily due to high rates of growth and metastasis (3). The most common sites for initial melanoma metastases are the skin, subcutaneous tissue, and lymph nodes, while distant and visceral metastases are most common in the lung, brain, liver, gastrointestinal tract, and bone (4). Organ-specific colonization suggests that the tumor cells adapt their cellular phenotype to suit the demands of the organ colonized (5, 6). The diversity of phenotypes of different metastatic locations is reflected in the patient's clinical outcome; for example, patients with lung metastases present better overall survival (OS) than patients with other visceral metastases, as reflected in the melanoma tumor–node–metastasis staging (7–9). However, the underlying molecular differences between metastatic locations and the association with treatment response and prognosis have not been studied on a global scale.

In the past decade, there has been a significant improvement in the treatment and OS of patients with advanced-stage melanoma, with the development of two main therapeutic strategies: MAPK pathway–targeted therapy (e.g., vemurafenib) for patients with BRAF V600E mutation and immunotherapy. Immunotherapy includes two main approaches: immune checkpoint inhibitors, primarily targeting the cytotoxic T lymphocyte-associated protein 4 (CTLA-4) or the programmed death 1 (PD-1) and adoptive cell transfer (ACT) of tumor-infiltrating lymphocytes (TILs). Although targeted therapies achieve rapid clinical responses, most patients develop acquired resistance, leading to tumor relapse. In contrast, immunotherapy achieves long-term tumor regression in approximately 50% of patients (10–14). Interestingly, the treatment of BRAF-mutant tumors with MAPK inhibitors (MAPKis) has been shown to elevate the immune activity in the tumor microenvironment (15), and it has been suggested that long-term benefit from MAPKis is strongly associated with immune-related adverse events (16). In agreement, a combination of MAPKis and checkpoint inhibitors led to a synergistic antitumor effect in mouse models (17–20). In contrast, prolonged treatment and developed resistance to MAPKis reduced immune activity in the tumor microenvironment (21–24). Therefore, sequential treatment combinations require deeper molecular understanding of short- and long-term treatment responses.

In a previous study of melanoma response to immunotherapy, we showed a connection between mitochondrial activity and lipid metabolism, and immune activity and response to TILs and anti–PD-1 treatment (25). Association between immune and metabolic processes has also been shown to be affected by metastatic location, and specifically brain metastases (26), by the BRAF mutation status (27, 28), and resistance to MAPKis (21). In this work, we assembled an integrated cohort from both published (25) and unpublished resources. It includes 185 formalin-fixed, paraffin-embedded (FFPE) samples originating from patients with advanced-stage melanoma treated with TILs, anti–PD-1, or anti–CTLA-4 therapy. We performed clinical proteomic analysis of melanoma to associate the clinical parameters of metastatic location, BRAF status, survival, and immunotherapy response with the molecular tumor profiles. Together, using these analyses, we found protein features associated with clinical parameters that may influence response and resistance to immune and targeted therapies.

Clinical samples were prepared as described previously (25). Briefly, 185 stage IV samples were macrodissected from FFPE tumor sections to enrich for tumor cells. Of the 185 samples, 116 were included in Harel and colleagues (25), and are referred to as the Harel cohort, while 69 samples are newly added and are referred to as the Beck cohort. Samples were taken mainly shortly before (n = 173) or after (n = 12) the indicated immunotherapy regimen (Supplementary Table S1A). The study was conducted in accordance with the Declaration of Helsinki, upon approval of the Institutional Review Board (IRB) Committees of the Sheba Medical Center (Tel Hashomer, Israel) and Tel Aviv University (Tel Aviv, Israel), and upon the patients' written informed consent. To obtain accurate proteome quantification, we designed a super-SILAC mix, which was composed of five SILAC-labeled melanoma cell lines that served as a reference for normalization (29, 30). For the Beck cohort, extracted tissues were lysed with 50% 2-2-2 trifluoroethanol and 25 mmol/L ABC. Samples were incubated for 1 hour at 99°C, sonicated, and centrifuged at 20,000 × g. Following protein determination using the fluorescamine protein quantification assay, the heavy-labeled super-SILAC standard lysate and the tumor lysates were mixed at a 1:1 ratio. The super-SILAC mix was prepared as described previously (25). Proteins were reduced using 5 mmol/L DTT, alkylated with 15 mmol/L iodoacetamide, diluted with 50 mmol/L ammonium bicarbonate, followed by overnight digestion with Lys-C-Trypsin mix (1:100 enzyme-to-protein ratio) and Trypsin (Promega, 1:50 enzyme-to-protein ratio). Peptides were then acidified with 1% trifluoroacetic acid, separated into five fractions using strong cation exchange (SCX) chromatography in a stage tip format, and purified on C₁₈ stage tips.

Mass spectrometry (MS)–based proteomic analyses were performed as described in Harel and colleagues (25). Briefly, peptides were separated by reverse-phase chromatography using the nano-ultra high-performance liquid chromatography system (Easy-nLC1000, Thermo Fisher Scientific), coupled to the Q-Exactive HF or Q-Exactive Plus Mass Spectrometers (Thermo Fisher Scientific; ref. 31). In the Beck cohort, fractionated peptides were separated with 140 minutes linear gradients of water/acetonitrile; one fractionated sample was analyzed with 220 minutes gradient and seven samples with low amounts were analyzed with 240 minutes gradients without prefractionation. The resolutions of the MS and MS-MS spectra were 70,000 and 17,500 for Q-Exactive Plus, respectively. The resolutions of the MS and MS-MS spectra were 60,000 and 30,000 for Q-Exactive HF, respectively. The m/z range was set to 300–1,700 or 380–1,800 Th. MS data were acquired in a data-dependent mode, with target values of 3E+06 and 1E+05 or 5E+04 for MS and MS-MS scans, respectively, and a top-10 method.

Raw MS data are available via ProteomeXchange with identifiers PXD006003 and PXD020618.

MS raw files of all samples were jointly analyzed by MaxQuant (32) version 1.6.2.6 and the integrated Andromeda search engine (33). MS-MS spectra were referenced to the UniProt human proteome (http://www.uniprot.org/) released in April 2018. FDR of 0.01 was used on both the peptide and protein levels. Peptides were allowed to have methionine oxidation and N-terminal acetylation as variable modifications and cysteine carbamidomethyl as a fixed modification. The “match between runs” option was enabled to transfer identification between separate LC/MS-MS runs based on their accurate mass and retention time after retention time alignment. The settings for the SILAC-labeled tumor sample runs included Lys-8 and Arg-10 as heavy labels.

Analyses were performed using the Perseus software version 1.6.2.3 (34), R environment, GraphPad Prism (www.graphpad.com), and IBM SPSS software. For all proteomic analyses, the proteinGroups output table was used. Reverse proteins, proteins that were only identified by site, and potential contaminants (excluding keratins) were filtered out, and normalized ratio tumor/SILAC data were log₂ transformed. We identified 10,178 proteins in total, but to retain high-quality quantifications, protein groups were filtered to have valid values in at least 60% of the samples, reaching a list of 4,883 protein groups, which were further used for all downstream analyses. To minimize technical variability and eliminate the risk of global differences between samples, data were normalized by subtracting most frequent value (modal value) in the distribution of each sample. The initial dataset consisted of 185 samples and we excluded samples with less than 3,500 proteins (11 samples). The remaining 174 samples were taken for further analysis. Low expression proteins can lead to missing values, therefore, missing value imputation was performed sample wise by drawing values from a normal distribution, with a downshift of 1.6 SDs and a width of 0.4 of the original ratio distribution width of each sample (34). Batch effect originating from dataset integration was corrected by performing linear modeling using the R limma package (ref. 35; Supplementary Fig. S1A and S1B; Supplementary Table S1A). For the organ tissue contamination test (Supplementary Fig. S1C and S1D), we used organ-specific transcriptomic-based gene signatures from the Human Proteome Atlas database (http://www.proteinatlas.org). Signatures were based on the “enriched genes” signature, a list of genes that have at least 4-fold higher mRNA level in a selected tissue, when compared with any other tissue. For the skin comparison, we omitted subcutaneous samples from the “others” group due to their inherent similarity. For the lung comparison, we used the elevated genes signature because the enriched genes signature contains only 13 genes that were mostly absent in our data. The median expression level of the proteins in each signature was calculated separately for each sample and compared using Student t test. Student t tests and ANOVA tests were performed using the R limma package and Prism software. One-dimensional (1D) annotation enrichment analyses, two-dimensional (2D) annotation enrichment analyses (36), and Fisher exact test enrichment analyses were performed in the Perseus software with Benjamini–Hochberg FDR < 0.02. The 1D annotation enrichment analyses were performed on the median fold changes between groups and tested for significant differences between the distribution of proteins of any category [gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG)] and the overall protein ratio distribution (Wilcoxon–Mann–Whitney test). The 2D annotations enrichment analysis is similar to the 1D test, but examines the differences between two-sample ratios simultaneously (36). Kaplan–Meier and log-rank tests were performed using Prism and the R's survival and survminer packages. Logistic regression was performed in SPSS (Supplementary Table S2E). For the upregulated protein network analysis, we used the STRING database to derive the protein–protein interaction information within each group. The STRING database determines protein–protein interactions based on publicly available sources and data mining (37). Next, to combine two networks, we used the Cytoscape software (38) and the DyNet plugin (39). Biological annotations were taken from GO (40) and KEGG (41).

Weighted gene correlation network analysis (WGCNA) was performed in the Perseus software (42, 43). We used a soft-threshold beta power = 12 to create a robust signed network (Supplementary Fig. S2A). We calculated the Pearson correlations between the module eigengene and clinical annotations. For the unsupervised clustering approach of the lung metastases, we applied the consensus clustering algorithm (44) using R package ConsensusClusterPlus (45). Parameters used were: maximal number of clusters, six; number of iterations, 1,000; subsampling fraction, 0.8; clustering algorithm, hierarchical; and distance matrix, Pearson correlation. Before applying the algorithm ratios across proteins, samples were Z-score normalized.

FFPE tumor blocks were obtained from the Institute of Pathology at the Sheba Medical Center (Tel Hashomer, Israel). For the HLA-DR staining, five slides were stained per group (lung and lymph node) with anti–HLA-DR antibody (ab166777, Abcam). Staining was performed using the BOND-III Automated Staining Platform (Leica Biosystems), following selected protocol for the Bond Polymer Refine Red Detection Kit (Leica Biosystems). For the CD8 staining, three slides were stained per group (lung and lymph node) with anti-CD8 antibody (CRM311, BioCare Medica), one slide from the lung group was already published in Harel and colleagues (25). Staining was performed by using Benchmark XT Staining Module (Ventana Medical Systems). All slides were scanned using the Leica Aperio VERSA Digital Pathology Scanner (Aperio Technologies). Staining quantification was performed by using the Aperio eSlide Manager software. Percentages of positively stained cells were used for downstream statistical analysis (Student t test).

Patient cohort was derived from the National Cancer Database (NCDB), a hospital-based cancer registry, with data assessed from 2004 to 2015 (46). The NCDB captures more than 70% of cancer diagnoses in the United States from >1,500 hospitals with cancer programs accredited by the American College of Surgeons and the American Cancer Society. The cohort included all individuals with stage IV melanoma who received immunotherapy as first-line treatment. The study was approved by the IRB at the Hospital of the University of Pennsylvania (Philadelphia, PA).

The primary exposure of interest was the site of metastasis. All four metastatic sites defined in the NCDB were used for the analysis (lung, liver, bone, and brain). Additional covariates included age, sex, race (categorized as White, African American, or other/unknown), and patient comorbidities (Charlson–Deyo comorbidity condition).

The primary outcome of interest was median OS, measured from the time of cancer diagnosis until death from any cause or last follow-up.

Patients were grouped, according to time of diagnosis, into those who were diagnosed between 2004 and 2013, when anti–PD-1 was not standard of care in the United States, and 2014–2015, when anti–PD-1 became FDA approved. Baseline characteristics were compared using χ² test for categorical variables and Wilcoxon rank-sum test for nonnormal continuous variables. Differences in OS were compared between patients with either site of metastasis. HR and 95% confidence intervals were calculated using the Cox proportional hazards model. Survival curves for time-to-event variables were estimated using the Kaplan–Meier method. All statistical analyses were performed using STATA/IC software 15.0 (StataCorp). A two-sided P ≤ 0.05 was used to define significance.

For the BRAFi-resistant (BRAFi-R) model, we chose the A375 cell line, which is BRAF mutated and BRAFi sensitive. The cells were cultured in DMEM (Biological Industries) supplemented with 10% FBS (Biological Industries) and 1% antibiotics (Biological Industries). The cells were routinely verified to be Mycoplasma free by using PCR Detection Kit (Hy-Mycoplasma Detection Kit, Biolabs). Cell line authentication was performed at the Genomics Core Facility of BioRap Technologies and the Rappaport Research Institute (Technion, Israel), using the Promega Power Plex 16 HS kit to determine short tandem repeat profiles. To obtain resistant cell line, we challenged the cells with increasing concentrations of vemurafenib (PLX4032, Selleckchem). Starting with a concentration matching the IC₅₀ values of the naïve cells, (0.14 μmol/L), vemurafenib concentrations were increased every 7 days, for 9 weeks. We ended the selection protocol after 9 weeks when we observed a 100-fold increase in the IC₅₀ value, finally treating A375 cell line with 5 μmol/L of the drug (Supplementary Fig. S7A and S7B). In parallel, we maintained the parental cell lines for 9 weeks as a control.

A375 vemurafenib–sensitive cell line and A375 vemurafenib–resistant cell lines were cultured in DMEM supplemented with 10% FBS, 1% antibiotics, and with or without the indicated vemurafenib treatment (0.5 or 20 μmol/L) for 72 hours. Cells were gently detached from the plates with 1 mmol/L EDTA, washed with PBS, and incubated with HLA antibodies (BD Biosciences, catalog Nos. 560896 and 560965) at a 1:10 ratio in FACS buffer containing 1% dialyzed FBS, 1 mmol/L EDTA, and 25 mmol/L HEPES in PBS. Measurements were performed using the CytoFLEX-4L Flow Cytometer (Beckman Coulter). Only live cells were analyzed by negative selection using DAPI. Three biological replicates were analyzed.

Aiming to associate molecular profiles of melanoma with clinical parameters, we assembled an integrated melanoma proteomic dataset, consisting of 185 metastatic melanoma tumors. We combined published and unpublished datasets analyzed in our laboratory. The published data were obtained from Harel and colleagues (25), with a total of 116 samples (42 TIL-treated and 74 anti-PD1–treated patients). Our final cohort was composed of 185 advanced-stage melanoma samples and included 86 patients treated with TIL-based therapy, 79 patients treated with anti–PD-1, and 20 patients treated with anti–CTLA-4. Eleven of the samples were filtered because of poor quality, and 174 samples were used for downstream analyses. Samples were taken mainly shortly before (n = 162) or after (n = 12) the indicated immunotherapy regimen (Fig. 1A; Supplementary Table S1A). Each treatment group was separated into immunotherapy responders (complete and partial responders) and nonresponders (stable and progressive disease; Fig. 1B). Each of the immunotherapy treatment groups presented high variability of prior treatments and metastatic locations, which provided another layer of analysis of the tumor proteome (Fig. 1B). First, we characterized the tumors according to the metastatic site (a minimum of five tumors/group) and focused on metastases excised from lymph nodes, subcutaneous, skin, lung, liver, brain, bone, and small bowel (Fig. 1C; Supplementary Table S1A). Next, we characterized the cohort according to the BRAF mutation status and prior treatment with MAPKis (Fig. 1D; Supplementary Table S1A). From the 174 samples, 64 were BRAF mutant and 90 were BRAF wild-type (WT; 20 patients were with unknown BRAF status). Among the BRAF-mutant samples, 19 were treated with MAPKis, and were considered treatment resistant at the time of tumor resection (4 patients had primary resistance and 15 developed acquired resistance), 42 were treatment naïve, and three were with unknown treatment status (Fig. 1D; Supplementary Table S1A). First, aiming to validate the reliability of the combined cohort, we examined the influence of known biomarkers on prognosis, and found longer survival in patients with higher expression levels of proteins involved in IFNγ signaling, MHC class I, and MHC class II (Fig. 1E).

Initial bioinformatic analysis examined the functional differences between different metastatic locations. A comparison between the median fold change of the protein levels in each site with all other sites showed location-specific enriched processes (1D annotation enrichment analysis, FDR q-value < 0.02; Fig. 2A; Supplementary Table S1D). Remarkably, we found a general association between normal tissue function and the processes enriched in the metastases (e.g., fatty acid metabolism in the liver). Given the precise dissection of the tumor regions, we speculated that it reflects the effects of the microenvironment on the tumor cells. To verify that enriched processes do not result from tissue contamination, we used organ-specific signatures and compared the median expression levels of the signature proteins in every organ with all other tissues (Supplementary Fig. S1C; Supplementary Table S1E). Analysis of the enriched genes signature, based on the Human Protein Atlas gene expression datasets, showed insignificant differences in signature expression between the organ-specific location and the rest of the dataset, except for brain metastases that were found to have low expression levels of the organ signature. Significance, in this case, was obtained because of one outlier brain metastasis sample (Supplementary Fig. S1C). The liver signature showed higher expression levels in liver metastases, although the difference was insignificant (Supplementary Fig. S1C). To better understand whether this phenomenon is derived from liver contamination or tumor cells' adaptation, we divided the liver signature into two groups: proteins involved in coagulation and inflammation, and proteins involved in metabolic processes (Supplementary Table S1E). Interestingly, the metabolic signature was significantly higher in liver metastases, while coagulation and inflammation signatures were similar to the metastases from other tissues, supporting tumor cell adaptation and negligible tissue contamination (Supplementary Fig. S1D).

Examination of different metastatic locations found the highest enrichment of immune-related processes in lung metastases, represented by MHC and TAP protein complexes and IFNγ signaling (Fig. 2A; Supplementary Fig. S1E). These results may provide a functional explanation for the favorable clinical outcome of melanoma lung metastases (7–9). High MHC levels and immune activity were also found in skin and subcutaneous metastases, while low immune activity was observed in lymph node, brain, and small bowel metastases (Fig. 2A; Supplementary Fig. S1E). In addition to low immune-related protein expression, and in agreement with Fischer and colleagues (26), brain metastases also showed high oxidative phosphorylation (OXPHOS) and NADH dehydrogenase complex (Fig. 2A; Supplementary Fig. S1F). Higher mitochondrial activity was also enriched in lung metastases, suggesting a positive association between immune and mitochondrial activities in the lung, while having a negative association in brain metastases. Liver metastases showed high fatty acid metabolism (Fig. 2A; Supplementary Fig. S1G), while the highest protein translation activity and extracellular matrix (ECM)-related processes were shown in lymph node metastases (Fig. 2A). The highest enrichment of proliferation-related annotations was observed in bone and small bowel metastases (Fig. 2A; Supplementary Fig. S1H).

Next, we took a complementary approach and performed WGCNA. This method separates the proteomic data into modules of highly correlated proteins. Every module has a module eigengene, which enables the calculation of correlation to the clinical features of the samples (42). We found nine protein modules that had significant correlations between the module's eigengene and the metastatic location (positive or negative correlations, Pearson correlation P < 0.05; Fig. 2B; Supplementary Fig. S2A and S2B; Supplementary Table S1F). Enrichment analysis of the proteins in each module reinforced the previous analysis (Fisher exact test enrichment, FDR q-value < 0.02). For example, higher immune activity was found in lung and subcutaneous metastases and lower immune activity in lymph node, brain, and small bowel metastases (Fig. 2B; Supplementary Table S1F).

Supervised comparison among all metastatic locations identified 44 significantly changing proteins (ANOVA FDR q-value < 0.1; Fig. 2C; Supplementary Table S1G). Among them, we found higher expression of the SLC2A1/GLUT1 glucose transporter in the brain metastases, higher expression of translation initiation factors (EIF4E2 and EIF2A) and RNA binding protein (RBM3) in subcutaneous, skin, and lymph node metastases, and higher expression of HLA class II (HLA-DRB1) in lung metastases. In addition, we detected higher expression of proteins involved in lipid metabolism in liver metastases. Surprisingly, among these proteins, we found ACAT1, which we have shown previously to be involved in improved immunotherapy treatment response (25). A closer examination of ACAT1 expression across all metastatic sites showed that it is significantly associated with better survival in the entire dataset, and especially in lung metastases (Fig. 2C; Supplementary Fig. S2C). These results show that biomarkers may largely depend on the metastatic location.

To validate the finding of increased immune activity in lung metastasis, we performed IHC staining to compare the expression of HLA-DR and CD8 in lung and lymph node metastases (Fig. 3A and B). Immunostaining showed significantly higher levels of both markers in lung metastases, in agreement with the proteomics data, and emphasized the variability of cancer–immune interactions in distinct metastatic locations.

Given the known immune and metabolic effects of BRAF mutation status and targeted therapy, we examined whether these clinical parameters may have confounded our conclusions. To that end, we examined the percentage of BRAF mutation carriers and prior MAPKi-resistant (primary or acquired) individuals in each metastatic location group (Supplementary Fig. S2D). Although most of the metastatic locations had similar percentages of BRAF-mutant and prior MAPKi-treated patients, the small bowel group had a high proportion of patients with MAPKi-acquired resistance (three of five patients), suggesting a potential confounding factor to the small bowel results. Taken together, these results show the importance of the metastatic microenvironment in determining the tumor phenotype and the potential influence on treatment response and patient prognosis.

To better understand the clinical effect of the unique processes identified for each metastatic location, we analyzed clinical data from the NCDB, a hospital-based cancer registry (46). We analyzed the clinical data from 526,167 individuals that were diagnosed with melanoma during the years 2004–2015. In this cohort, 19,141 individuals had a metastatic disease upon diagnosis, of whom 3,407 received first-line immunotherapy and were used for further analysis (Supplementary Fig. S3A). Because NCDB does not provide the detailed treatment information for each patient, we grouped the patients according to the time of diagnosis, in the context of the approved treatments at each timepoint. The two patient groups included (i) patients diagnosed between 2004 and 2013, when anti–PD-1 immunotherapy was not standard of care in the United States and (ii) patients diagnosed between 2014 and 2015, when anti–PD-1 was FDA approved. A total of 55% of the study cohort (1,871) were treated during the years 2004–2013 and 45% (1,536) during the years 2014–2015 (Supplementary Fig. S3A). Patient characteristics are presented in Supplementary Fig. S3B according to the year of diagnosis. The median follow-up time was 12.7 months [interquartile range (IQR), 6.2–26.5 months]. OS for the entire study population was 15.8 months (IQR, 7–50.4 months). In agreement with Eggermont and colleagues (10), survival improved among individuals diagnosed from 2014 to 2015 compared with 2004–2013 (median OS, 18.2 months; IQR, 6.6–46.8 vs. 14.8 months; IQR, 7.1–40.9 months, respectively; P < 0.001), however, beyond these known differences, we examined the association of this improvement with the site of metastasis. We, therefore, compared patients with first-line immunotherapy and single-site metastases diagnosed before and after anti–PD-1 became standard of care in the United States. In a subgroup analysis according to the metastatic site, only individuals who had single-site metastasis to the lungs at the time of diagnosis had an improved survival upon anti–PD-1 routine implementation (median OS, 29.3 vs. 22.1 months; P = 0.01), while individuals with metastasis to the liver, brain, or bone did not have a similar benefit (median OS, 16.7 vs. 13.3; 20.2 vs. 21; and 18.3 vs. 16.5, respectively; Fig. 3C). These results show the benefit for survival of anti–PD-1 treatment for lung metastases compared with other visceral metastases. In addition, results stand in line with the proteomic results and highlight the contribution of metastatic site microenvironment to disease prognosis and response to immunotherapy.

Focusing on single-metastatic sites, we examined whether we can further cluster the data in an unsupervised manner and associate clusters with clinical parameters. Because of the lung's unique immunophenotype, we focused our analysis on this site. Unsupervised clustering of lung metastases resulted in three clusters that best separate these tumors. The optimal number of clusters was selected using quantitative evaluation of the area under the cumulative distribution curve (Fig. 4A; Supplementary Fig. S4A). Kaplan–Meier analysis depicted a significant difference between the patient clusters (P = 0.0382); with cluster 1 having the highest survival and cluster 3 having the lowest survival (Fig. 4B). Next, we examined whether the clusters differ in their immune infiltration using a modified ImmuneScore (26, 47), which originally quantified the presence of immune cells in the tumor on the basis of transcriptomic expression values of 141 immune-related genes. Of these 141 genes, 36 were identified in the proteomic dataset and showed significantly higher expression in clusters 1 and 2 compared with 3 (ANOVA P = 0.013; Fig. 4C). The 36 proteins identified in our data included proteins from both adaptive and innate immune activity (Supplementary Table S2A).

Examination of the significantly changing proteins between the three lung metastasis clusters identified 277 significantly different proteins (ANOVA FDR q-value < 0.1). Hierarchical clustering showed that the significantly changing proteins divided into three protein clusters in accordance with the sample clusters, and can indicate cellular processes that differ between clusters (Fig. 4D; Supplementary Table S2B). ANOVA results showed higher expression level of proteins involved in adaptive immune activity in cluster 1 (e.g., TAP2 and HLA-DRB1), while higher expression of proteins involved in innate immunity (e.g., ORM2, MPO, and GSDMD), proteasome ubiquitin, regulation of nitric oxide synthases (e.g., DDAH1 and DDAH2), and ribosome was found in cluster 2 (Supplementary Fig. S4B; Supplementary Table S2B). 2D annotation enrichment analysis (36) to compare the three clusters showed similar results, for example, relatively higher immune activity (e.g., antigen presentation and IFNγ signaling) and mitochondrial activity in cluster 1 and higher proteasome activity in cluster 2 (Fig. 4E; Supplementary Fig. S4C; Supplementary Table S2C). Cell growth–promoting proteins and cell growth–related pathways (e.g., chromosome organization and DNA binding) were higher in cluster 3 (Fig. 4D and E; Supplementary Table S2B–S2C), in agreement with the more aggressive phenotype and lower patient survival.

Next, we examined whether the lung metastasis clusters can highlight proteins that correlate with survival. To that end, we conducted a low stringency Student t test with a nominal P value cutoff (P < 0.01) to look for differentially expressed proteins between clusters 1 and 3 (high and low survival, respectively; Supplementary Table S2D). Kaplan–Meier test of the 133 differentially expressed proteins identified eight proteins with a significant positive correlation with survival and 17 proteins with a significant negative correlation with survival (P < 0.05; Fig. 4F; Supplementary Fig. S5A). Next, we examined the specificity of these 25 proteins with lung metastasis, and found only six of them to be associated with survival in the whole dataset; two were also significant in subcutaneous metastases and one was significant in lymph node metastases (Fig. 4F; Supplementary Fig. S5B–S5D). Similar unsupervised clustering of lymph node and subcutaneous metastases, which also form sufficiently large groups of tumors, did not identify robust clusters and showed no significant difference in survival (Supplementary Fig. S6A–S6D). Together, the focused analysis of lung metastasis unraveled the clinical importance of proteomic heterogeneity even within a single metastatic site and the importance of tissue specificity predictive markers.

Our results showed the effect of the metastatic location on the proteomic profiles and response to immunotherapy. We next examined the proteomic effects of additional clinical parameters, primarily BRAF status. Given that 38% of the patients (n = 64) in our cohort were BRAF mutation carriers, and of those, 23% (n = 15) had acquired MAPKi resistance and 6% (n = 4) had primary MAPKi resistance, we postulated that examination of the clinical and proteomic data may highlight novel functional associations with immunotherapy response. Given the small number of samples that had primary resistance to MAPKi, we focused on the effect of acquired resistance in the analyses described below. To examine the independent influence of BRAF mutation and acquired MAPKi resistance on response to immunotherapy, we performed logistic regression analysis and controlled for age, gender, and immunotherapy treatment group (TIL, anti–PD-1, anti–CTLA-4; logistic regression P < 0.05; Fig. 5A; Supplementary Table S2E). In agreement with Ackerman and colleagues (23), statistical analysis of the clinical data showed a significantly higher response to immunotherapy in the BRAF mutation carriers that were MAPKi naïve, compared with MAPKi-resistant BRAF mutation carriers and compared with the WT BRAF patients.

To associate the BRAF groups with their molecular functionalities, we performed a 1D annotation enrichment analysis (FDR q-value < 0.02) on the median fold change of the proteomic data. A comparison of the BRAF WT with treatment-naïve patients with BRAF-mutant tumors showed higher levels of immune processes, mitochondrial respiration, and ECM in the BRAF-mutant treatment-naïve group (Fig. 5B; Supplementary Table S2F). To assess the influence of developed resistance on MAPKi, we compared the MAPKi-resistant group with the naïve group. Similar to the previous results, we found higher immune processes, movement, and ECM in the naïve group, however, we found mitochondrial translation to be enriched in the resistant group (Fig. 5C; Supplementary Table S2G). These results suggest that developed resistance to MAPKi reduces the tumor immunogenicity, but increases mitochondrial activity. To further examine this phenomenon, and to identify the key proteins contributing to these processes, we plotted a unified protein–protein interaction network of the proteins that were higher in the MAPKi-naïve group compared with the two other groups: BRAF WT and MAPKi-resistant patients (P < 0.05; Fig. 5D; Supplementary Table S2H–S2I). In agreement with the enrichment analysis results, we found an elevation in mitochondrial OXPHOS proteins only in comparison with the BRAF WT group and elevation in HLA proteins and connection to the ECM and cell movement proteins in comparison with both groups.

To further examine our hypothesis that developed resistance to MAPKi reduces immune activity, we established a BRAFi-R model by challenging A375 melanoma cell line with increasing vemurafenib concentrations (Supplementary Fig. S7A and S7B). Following this selection protocol, we examined HLA class I and HLA class II expression in the vemurafenib-naïve and -resistant cells with and without acute vemurafenib treatment (naïve and resistance cells IC₅₀ values). In agreement with the clinical proteomics data, we found a significant decrease in HLA expression upon developed resistance to vemurafenib (Fig. 5E and F; Supplementary Fig. S7C and S7D). However, upon acute treatment with vemurafenib, HLA-ABC and HLA-DR were elevated primarily in the control cells, and resistant cells required 40-fold higher concentrations of the drug to increase HLA-ABC (Fig. 5E and F; Supplementary Fig. S7C and S7D). Altogether, these results provide the molecular associates with treatment resistance and suggest a clinical advantage to coordinate administration of MAPKis and immunotherapy and emphasize the need to improve the understanding of treatment combination timing.

Despite the dramatic improvements in metastatic melanoma treatment and outcome, a large proportion of patients remains without a suitable treatment solution. Therefore, advances in the understanding of resistance mechanisms to those treatments is an urgent need to provide suitable treatment in an individualized manner and to develop new treatments for resistant patients. In our previous work, we showed that the metabolic state of melanoma cells influences the immune activity within the tumor, suggesting a mechanism of improved response to immunotherapy (25). In this work, we expanded our cohort to 185 patients with advanced-stage melanoma, aiming to associate clinical subgroups with immunotherapy response, and identify the molecular determinants of each group. The proteomic data provided novel molecular explanations to critical clinical phenomena, including the effects of sequential treatments and metastatic locations.

Examination of the proteomic differences between tumors from distinct metastatic locations identified a unique combination of cellular processes for each organ. In the lung, we found higher immune activity and higher mitochondrial activity (Fig. 2A–C). Orthogonal IHC analysis strengthened this finding and found higher levels of HLA-DR and CD8 in lung metastases when compared with lymph node metastases (Fig. 3A and B). These results provide a potential explanation for the favorable clinical outcome of patients with lung metastases in comparison with other metastases (Fig. 3C; refs. 8, 9). Low expression of immune-related proteins, together with high mitochondrial activity, was observed in brain metastases, in agreement with Fischer and colleagues (26). Furthermore, we found differences in lipid metabolism, translation, proliferation, and adhesion processes in the different metastatic sites (Fig. 2A–C). Because our proteomic cohort included only a small number of samples originating from brain, liver, bone, and small bowel metastases, follow-up studies could extend the functional analysis of metastases to these organs. In addition, further research of different metastatic locations from the same patients can improve the understanding of tumor adaptation process to different environments.

Unsupervised clustering analysis identified clinically relevant heterogeneity in lung metastases based on the proteomic data (Fig. 4A–E). In agreement with our previous work (25), high levels of adaptive immune and mitochondrial proteins were associated with the group with the highest survival. Interestingly, high levels of proteins related to innate immunity were associated with the middle survival cluster, suggesting that different immune-related pathways affect the clinical outcome (Fig. 4D; Supplementary Fig. S4B). Furthermore, one of the main characteristics of the middle survival cluster is high levels of proteins related to the ubiquitin proteasome system, which emphasizes the protein-level regulation of this cluster (Fig. 4D and E). The lower survival cluster was mainly enriched with nuclear proteins and DNA replication processes, thus suggesting a more proliferative phenotype (Fig. 4D and E).

The high proteomic variability among metastatic sites suggests that these have to be considered in the pathological examination, treatment decisions, and the search for predictive biomarkers. Our analyses identified high liver expression levels of ACAT1, which we have shown previously to increase immunotherapy response in the analysis of the entire cohort (25). Here, the tissue-specific analysis showed higher ACAT1 expression in the liver, which is an organ commonly associated with worse survival and shows ACAT1 correlation with better survival mainly in lung metastases (Fig. 2C; Supplementary Fig. S2C). Analysis of proteins that correlate with survival in the lungs showed that only a minority of the proteins are associated with survival in other tissues (Fig. 4F). This conundrum emphasizes the importance of the tissue context in the identification of predictive markers. In total, these analyses highlight the potential to develop tissue-specific therapies (or combination therapies), specifically for those presenting poorer responses. This heterogeneity may be specific to metastatic melanoma, which presents with diverse metastatic locations. However, similar heterogeneity may be relevant in other highly metastatic tumors as well, and should be accounted for in clinical diagnostics and treatment decision-making.

Cohort segregation according to BRAF treatment status unraveled the functional processes associated with resistance to MAPKi treatment. A similar decrease in immune activity and an increase in mitochondrial activity upon acquired resistance to MAPKis have been shown previously using in vitro and immunofluorescence methods (21, 22, 48). In addition, preclinical models demonstrated that combination of MAPKis and immunotherapy agents can have an synergistic effect (17, 23, 49). Our results provide molecular foundation to the potentially additive effect of immunotherapy and MAPKis, when given together (11, 50), but reinforces the importance of treatment regimen and timing because developed resistance to MAPKis potentially reduces tumor immunity. Since the standard of care changed in the past few years from a single-BRAFi agent to a combination of BRAFis and MEKis, this needs to be taken into account in future laboratory and clinical experiments.

In conclusion, the work presented here shows the critical significance of the metastatic site, mutation status, and prior treatment resistance on the pathogenesis and the phenotype of advanced-stage melanoma, and sheds light on the complex heterogeneity of those tumors. The results reveal phenotypic and molecular foundations for urgent clinical phenomena. The clinical proteomic analysis also shows the importance of protein-level analysis, which integrates the genetic information of the tumors with the effect of the microenvironment. Importantly, these findings may contribute to the understanding of resistance mechanisms to immunotherapy treatments and the prognosis of those patients.

T. Geiger reported grants from Melanoma Research Alliance, Israel Science Foundation, Israel Innovation Authority, Samueli Foundation, and Emerson Foundation during the conduct of the study, and grants from European Research Council outside the submitted work. G. Markel reported grants and personal fees from BMS and Novartis, and personal fees from MSD, 4c Biomed, Biond Biologics, and Ella Therapeutics outside the submitted work. No disclosures were reported by the other authors.

L. Beck: Conceptualization, formal analysis, validation, investigation, writing–original draft. M. Harel: Investigation, methodology. S. Yu: Data curation. E. Markovits: Data curation. B. Boursi: Data curation, formal analysis, validation. G. Markel: Resources, data curation, funding acquisition. T. Geiger: Conceptualization, resources, supervision, funding acquisition, writing–original draft.

We thank the Geiger laboratory for fruitful discussions, and specifically Gali Yanovich-Arad for computational assistance. We thank Prof. Michal Lotem from the Hadassah Medical Center for consultation. We thank the Lemelbaum family for the generous continuous support. T. Geiger and G. Markel acknowledge funding from the Melanoma Research Alliance (Saban Family Team Science Award), from the Israel Innovation Authority, from the Emerson Collective Cancer Research Fund, and from the Israel Science Foundation- Israel Personalized Medicine Program (grant #no., 3495/19). G. Markel also acknowledges the Samueli Foundation grant for Integrative Immuno-Oncology. This work was supported by grants from Kamin grant of the Israel Innovation Authority (to T. Geiger and G. Markel) and Israel Science Foundation- Israel Personalized Medicine Program, IPMP (to T. Geiger, G. Markel, and B. Boursi).

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.