The immune response to melanoma improves the survival in untreated patients and predicts the response to immune checkpoint blockade. Here, we report genetic and environmental predictors of the immune response in a large primary cutaneous melanoma cohort. Bioinformatic analysis of 703 tumor transcriptomes was used to infer immune cell infiltration and to categorize tumors into immune subgroups, which were then investigated for association with biological pathways, clinicopathologic factors, and copy number alterations. Three subgroups, with “low”, “intermediate”, and “high” immune signals, were identified in primary tumors and replicated in metastatic tumors. Genes in the low subgroup were enriched for cell-cycle and metabolic pathways, whereas genes in the high subgroup were enriched for IFN and NF-κB signaling. We identified high MYC expression partially driven by amplification, HLA-B downregulation, and deletion of IFNγ and NF-κB pathway genes as the regulators of immune suppression. Furthermore, we showed that cigarette smoking, a globally detrimental environmental factor, modulates immunity, reducing the survival primarily in patients with a strong immune response. Together, these analyses identify a set of factors that can be easily assessed that may serve as predictors of response to immunotherapy in patients with melanoma.

Significance:

These findings identify novel genetic and environmental modulators of the immune response against primary cutaneous melanoma and predict their impact on patient survival.

See related commentary by Anichini, p. 2457

The presence of tumor-infiltrating lymphocytes (TIL) predicts better outcomes from primary melanoma (1, 2) and therapeutic benefit from checkpoint blockade is more likely if tumors are PD-L1–positive (3) in response to T-cell infiltration. Data have been published suggesting that higher mutational load is predictive of response to immunotherapy, and some studies with small numbers of patients have reported gene expression signatures with some predictive value (4, 5). However, the crucial need remains to identify the biological processes underlying “cold” unresponsive tumors. Bioinformatic analysis of large-scale “omic” datasets such as The Cancer Genome Atlas (TCGA) increasingly contribute to our understanding of tumor immunology (6, 7), but the tumors are highly selected/biased, at advanced stage and with limited clinical metadata. In this report, we have used transcriptomic data generated from 703 of the 2,184 participants in a population-based primary melanoma cohort (the Leeds Melanoma Cohort, LMC; refs. 8, 9) to explore the drivers of immune responses/failure at diagnosis, with the aim, ultimately, of improving adjuvant therapeutic choices.

In a previous report, we applied an approach to infer the tumor immune microenvironment described by Bindea and colleagues (10) and identified six immunologically different tumor subgroups (11). The immunome compendium used in that study contained genes specific to 24 immune cells (11). In this study, we utilized a refined version of the immunome compendium derived from a more extensive literature screening and covering 31 immune cell subtypes as published by Angelova and colleagues (12). We defined transcriptomic scores for these immune cells and used them to classify tumors with unsupervised methods to identify immunologically different subgroups. The classification was based on the immune cell scores generated from the expression of genes attributed to each cell subtype rather than on individual genes as we reported in our previous study (11). We postulated that reducing the number of dimensions prior to classification analysis could identify the tumor groupings with a clearer difference in survival, facilitating subsequent in-depth biological and epidemiologic characterization, directed toward the identification of candidate therapeutic targets.

There is evidence that environmental factors may modify immune responses to tumors (13). We have previously reported that smoking was associated with microscopic tumor ulceration and vitamin D was protective (9), and here we demonstrate the effect of smoking as a modifier of outcome within each immune subgroup.

The LMC transcriptomic data

The transcriptomic data from 703 tumors were generated and preprocessed as previously reported using the Illumina DASL whole-genome array (14, 15, 11). These data are accessible for the purposes of melanoma research from the European Genome–Phenome Archive with the accession number EGAS00001002922. All the survival analyses used melanoma-specific survival (MSS). The median follow-up for 703 samples at the time the dataset was fixed was 7.5 years. Detailed information about the cohort is provided in Supplementary Materials and Methods. The participants gave written informed consent and the study received ethical approval (MREC 1/03/57 and PIAG3-09(d)/2003) and was conducted in accordance with recognized ethical guidelines (Declaration of Helsinki).

Immune cell scoring

Angelova and colleagues generated a list of genes identified as specific to certain immune cells in the blood (the compendium of immune genes; ref. 12). These 1,980 genes were identified from the reports of 36 studies comprised of 813 microarrays generated from 30 purified immune cell subtypes [activated and memory B cells, activated CD4+ and CD8+ cells, central memory CD4+ and CD8+, cytotoxic cells, dendritic cells (DC), effector memory CD4+ and CD8+, eosinophils, immature B cells, macrophages, mast cells, monocytes, natural killer (NK) cells, NK 56 bright, 56 dim, and NK T, neutrophils, T cells, T follicular helper, T gamma delta, Th1, 2, 17, T regulatory cells (Treg), immature, plasmacytoid, and memory DCs], as well as genes (searched from literature) for myeloid-derived suppressor cells, resulting altogether in 31 subtypes.

From the initial list of genes, we removed those also strongly expressed (in the top 25%) in a melanocyte cell line (GSE4570) and in the melanoma cell lines, MEWO and SK-MEL28 (in-house data). In a second step, we removed cell subsets for which more than 90% of genes were eliminated (in the previous step) or if there was insufficient published evidence for the remaining genes to be considered representative of those cell types. We expected that expression of the majority of genes specific to a particular cell type would be positively correlated within the cell type, as this was the basis of gene selection in the Angelova and colleagues study. However, in our dataset, this was not always the case, so in a third step of filtering, we removed genes negatively correlating with the majority within each cell subset to reduce the risk of noise due to technical factors. After applying the filters described above, we devised a score for each immune cell type, calculated as the mean of expression values of all genes attributed to that cell, after z-score normalization of the log2-transformed gene expression data as described before (11). The scores were calculated in the LMC primary melanomas and in the TCGA metastatic melanomas. The reciprocal correlations of the genes within each immune cell score were compared between these two datasets.

Clustering of LMC tumors

We applied consensus cluster analysis (16) within the R package ConsensusClusterPlus (17) to classify primary melanomas of LMC based on their immune cell scores. This approach generates a varying number of clusters (to a fixed maximum number, K) using resampling of the data. It is widely used to find stable sample subgroups in high-dimensional data as a better alternative to the standard one-off clustering, which might be affected by random variation. In addition, consensus clustering offers useful metrics (see below) to indicate the optimal number of clusters, unavailable in standard clustering. K-means was chosen as the clustering algorithm with maximum K = 12, Euclidean distance, 5,000 repetitions, 80% genes, and tumor resampling. Examination of the consensus cluster matrices, the cumulative density function (CDF), and delta CDF (the change in the area under the CDF curve) allowed definition of the optimal number of tumor clusters (16).

Clustering replication in TCGA

We downloaded RNA sequencing gene expression and survival data for TCGA metastatic melanoma data (http://www.cbioportal.org/data_sets.jsp; 339 samples downloaded in 2016). We hypothesized that the immune subgroups observed in the primaries would be recapitulated in metastatic melanomas, and to test this hypothesis, we calculated cluster centroids (vector of cell score means within clusters) in the LMC dataset and utilized them to classify TCGA metastatic melanomas using the nearest centroid method, as described elsewhere (15). Briefly, immune cell scores were calculated in the TCGA data in a similar manner as in the LMC and each TCGA tumor was assigned to one of the new clusters, with which it had the strongest Spearman correlation.

Overrepresentation analysis and networks

To test the whole transcriptome differences between the immune subgroups in the LMC, the Kruskal–Wallis test was used for three groups, the Mann–Whitney U test was used for two groups, and Bonferroni correction was applied for multiple testing correction (0.05/29354 = 1.7 × 10−6, the number of probes tested was 29,354). To visualize the expression of significantly differentially expressed genes (DEG; excluding the compendium of immune genes) among the immune subgroups, these DEGs were hierarchically clustered and a heatmap was plotted. Reactome FiViz (18) and Centiscape (19) in Cytoscape (20) were utilized to analyze the protein–protein interaction (PPI) network and infer pathways enriched in the DEGs characterizing each immune subgroup. The networks were created on the basis of existing PPI networks built in Reactome FiViz, which covers over 50% of the human proteins. From the network, pathway enrichment was calculated at FDR < 0.001. To identify the most influential (hub) genes in the networks, the “betweenness” metric (indicating a key role in communication between proteins) was used as a centrality measure in Centiscape (19). Graphical adjustments for network visualization were made in Gephi software (21). The Spearman rank correlation was used to evaluate the correlation between the expression of a hub gene and the whole transcriptome patient-derived primary melanoma cell line cultures.

Primary melanoma cell lines validation experiment

As previously described, primary melanoma cells were isolated from surplus surgical specimens of the consenting patients [approved by the local Institutional Review Board (EK.647/800)] at the University of Zurich (Zurich, Switzerland), and maintained by the University Research Priority Program in Translational Cancer Research at the University of Zurich Hospital (Zurich, Switzerland; ref. 22). The cell lines were authenticated by Sanger sequencing for the oncogenes seen from the tumor and tested for Mycoplasma infection using the PlasmoTest Kit (InvivoGen). Cells were passaged about five times before RNA extraction. Details are provided in Supplementary Materials and Methods.

IHC validation

The most influential genes identified in our analyses were further examined by IHC staining of sections of available primary tumors from the LMC, to assess the protein-level and gene expression correlations. Primary antibodies were supplied by Abcam: anti-MYC (ab32072), anti-NF-κB p105 (ab32360), anti-HLA-B (ab193415). The details of scoring are described in the Supplementary Materials and Methods. The Mann–Whitney U test was used to compare the mRNA level and the IHC scores. The nuclear staining scores in tumor cells and in TILs were compared using the Fisher exact test. The correlation of continuous scoring was tested using Spearman rank correlation. Images presented are of digitally scanned slides generated by the Digital Pathology Group, Leeds Institute of Medical Research at St. James's (Leeds, United Kingdom), using Aperio technology (Leica Biosystems).

Analysis of copy number alterations among the immune subgroups

We extracted copy number profiles in the genomic regions spanning the hub genes from the network analysis and compared them between immune subgroups and with gene expression and patient survival. Next-generation sequencing–derived copy number alteration (CNA) profiles were available from 276 tumor samples among the 703 transcriptomic-profiled samples. To evaluate the full extent of the role of structural variation, we expanded the gene list to other genes of the same family or the same pathway as the hub genes (plus NF-κB and its regulators, and IFNγ signaling genes). The association between CNAs and gene expression was tested using Fisher exact test. Because we have previously reported β-catenin signaling pathway to be upregulated in 42% of primary melanomas overall and in 73% of those with the worst outcome (11), we tested the overlap between the CTNNB1 expression signature with the immunosuppressive mechanisms identified in this study and their joint effect on survival (Cox proportional hazard regression). For the CNA visualization, the ComplexHeatmap package in R was used (23). The CNA data analysis in detail is included in Supplementary Materials and Methods.

Statistical analyses

The univariable Cox proportional hazards model was used to test the association between tumor immune subgroups and MSS in the LMC and overall survival (OS) in TCGA datasets. To test independence between the tumor immune subgroups and clinicopathologic factors, χ2 and Kruskal–Wallis tests were used. A univariable Cox proportional hazard model was used to test the prognostic value of the immune cell scores and the clinical and environmental factors [American Joint Committee on Cancer (AJCC) staging version 7, age at diagnosis (median: 58.34 years), sex, site of melanoma (limbs vs. rest), smoking never/ever (median duration of smoking in the smoking group was 23 years), vitamin D levels at recruitment (median level in winter: 39.5 nmol/L)] and a social status/deprivation index measured by Townsend score (24) in the whole LMC dataset. Subsequently, the significant clinical and environmental predictors were included in a multivariate model (adjusting for the immune clusters). The predictors with the strongest degree of independence in the whole dataset were jointly tested within each immune subgroup.

Devising a list of genes indicative of specific immune cells infiltrating melanoma

The first filtration step resulted in 458 genes representing 30 distinct immune cell subsets (subset 1, Fig. 1; Supplementary Data S1). The second step resulted in the elimination of scores for effector memory CD4+ T cells, activated CD8+ T cells, and activated CD4+ T cells (subset 2, Fig. 1; Supplementary Data S1). The plasmacytoid dendritic cell score (pDCs) was retained despite having only 1 attributed gene (IL3RA), as in the previous version of the immunome compendium (10), because it is known to be highly expressed in pDCs (25, 26). The final filtration step left 376 genes representing 27 immune cell subsets (subset 3, Fig. 1; Supplementary Data S1). We noted that when applied to TCGA transcriptomes, the correlation matrices between genes within each cell type also demonstrated a number of negatively correlated genes although fewer than in our primary melanoma cohort (data available upon request).

The association of 27 immune cell scores with survival

We tested the association of the immune cell scores with MSS (univariable analysis) and the results revealed that the majority of immune cell scores (17 of 27) in the LMC and (23 of 27) in TCGA were associated with improved survival after Bonferroni correction (27 tests, P < 0.002). For eight of the remaining 10 cell scores in the LMC, a similar protective effect was found, but the effects did not withstand adjustment for multiple testing (Supplementary Table S1). The survival analysis was repeated after removal of the 16 participants known to have received checkpoint therapies and the results for all the immune cell scores were essentially unchanged.

Identification of three prognostic immune subgroups

Consensus clustering analysis of tumor samples using the 27 immune cell scores identified three clusters with distinct immune phenotypes (Supplementary Fig. S1), which we termed low, intermediate, and high immune subgroups (Fig. 2A). Importantly, by classifying the TCGA metastatic melanomas in these three immune subgroups (supervised classification), we were able to replicate the results obtained in LMC with strong similarities in overall immunologic profiles (Fig. 2B). Furthermore, the three immune subgroups were associated with survival in both datasets: in the LMC, a significantly lower hazard of melanoma death was observed for patients assigned to the high compared with low and intermediate immune subgroups [HR = 0.5, P = 0.001 (95% confidence interval (CI) 0.3–0.7); HR = 0.6, P = 0.05 (95% CI 0.4–1.0), respectively; Fig. 2C].

For TCGA, tumors classified as high immune also exhibited a lower overall death hazard with HR = 0.3, P = 1.1 × 10−7 (95% CI 0.2–0.5) compared with those classified in the low immune subgroup. Tumors in the intermediate immune subgroup had a HR = 0.5, P = 4.6 × 10−5 (95% CI 0.4–0.7) when compared with those of the low immune subgroup (Fig. 2D).

The 3-class signature reported here (high, intermediate, and low immune) was concordant with the six consensus immunome clusters (CIC) we published previously (Cramer's V = 0.72; see Supplementary Fig. S2A; ref. 11).

However, there was only moderate concordance with another 3-class signature (immune, keratin, and MITF low) published by TCGA, (Cramer's V = 0.47). In essence, our high immune subgroup overlapped well with the TCGA immune class but the intermediate and low immune subgroups had much less overlap with TCGA groups (Supplementary Fig. S2B). Kaplan–Meier curves for our three immune subgroups were more clearly/significantly separated than the three TCGA classes. Generally, we observed the expected prognostic trend of the immune signature, but there was no difference between the immune and keratin groups of the TCGA signature (Supplementary Fig. S2C). There were 70 genes shared between our and TCGA signature (listed in Supplementary Data S2X).

The immune subgroups were associated with tumor thickness, TILs, and mitotic number

In the LMC data, the high immune, in comparison with the low and intermediate subgroups, featured consistently thinner tumors (Kruskal–Wallis, P = 0.004) and, crucially, more TILs reported by both clinical dermatopathologists (χ2, P = 4.0 × 10−7) and a single observer from our research group who was blinded to the transcriptomic data (χ2, P = 3.6 × 10−8; see Supplementary Table S2). The mitotic number was significantly lower in the high immune subgroup (Kruskal–Wallis, P = 2 × 10−4). The low immune subgroup had the lowest proportion of tumors harboring a BRAF mutation (40%) and the highest proportion with an NRAS mutation (30%), although these observations were only marginally significant (Supplementary Table S2). The recorded site of melanoma was significantly different across the immune subgroups, with primary tumors located at “rare” sites (not exposed to the sun) most frequently classified in the low (19%) compared with the intermediate (9.5%) and the high immune subgroups (8%;χ2P = 0.02). AJCC stage, patient sex, age at diagnosis, smoking status, and the levels of season-adjusted serum vitamin D were not significantly different between the subgroups (Supplementary Table S2).

Identification of MYC as the regulator of immune response in network analyses

We compared gene expression among the three immune subgroups (Supplementary Fig. S3). A total of 5607 DEGs across the genome were identified between the high versus low immune subgroups. The genes upregulated in the low immune subgroup (n = 3324) were predominantly associated with high proliferation and metabolic activity (hypergeometric test–adjusted P =10−14–10−7) with lower levels of expression of the genes coding immune checkpoint molecules (expressed by tumor) such as CD274 coding for PD-L1. The most enriched pathways were the citric acid (tricarboxylic acid) cycle and respiratory electron transport, mitochondrial translation, and mitosis pathways (Fig. 3A; Supplementary Data S2A). Network analysis of the genes enriched in the low immune subgroup revealed that the proto-oncogene MYC had the highest centrality (Fig. 3A). Unsurprisingly, the network analysis indicated that the genes upregulated in the high immune subgroup (n = 2283) were mostly involved in immune pathways (hypergeometric test–adjusted P = 10−14–10−10). The top enriched pathways were: IFNα/β signaling, antigen processing and presentation, IFNγ and NF-κB signaling (Fig. 3B; Supplementary Data S2B), with the nodal gene in this network being NFKB1 encoding the p105/p50 subunit of NF-κB (Fig. 3B).

The identification of MYC as the gene with the highest centrality in the low immune network suggested that it might fulfil a key role in immune evasion. We took an agnostic approach to testing correlations between MYC expression and the rest of the genome in transcriptomes from patient-derived primary melanoma cell lines (lacking immune cells; ref. 22). Genes were ranked according to their correlation with MYC and of 50 genes most significantly negatively correlated with MYC, one tenth were involved in antigen processing and presentation (HLA-B, HLA-C, B2M, TAP1, and ERAP1), with HLA-B representing the strongest results (R = −0.57; P = 1.6 × 10−10; Fig. 3C; Supplementary Data S2C). The correlation of MYC with HLA-B gene expression in the LMC was: R = −0.3 (P = 5.5 × 10−14).

An immunosuppressive effect of oncogenic MYC has previously been demonstrated, although by different mechanisms than suggested in this study: MYC was reported to increase the expression of genes encoding CD47 and PD-L1 on lymphoblastic leukemia cells (27). We tested this observation using Spearman rank correlation, but MYC expression did not significantly correlate or correlated negatively with CD47 or PD-L1 expression in either the primary melanoma cell lines (R = −0.17, P = 0.09; R = −0.16, P = 0.1, respectively) or in the LMC tumors (R = 0.04, P = 0.3; R = −0.17, P = 2.3 × 10−6, respectively).

mRNA gene expression correlation with protein scores, IHC

A subset of tumors was immunohistochemically stained using antibodies for key proteins. We found that the protein expression of MYC was localized to the tumor cell nuclei, whereas HLA-B localized to the tumor cell membrane and the expression of both proteins was positively associated with their mRNA transcripts (P = 0.056 and P = 0.002, respectively; Fig. 4). MYC staining was only detected in tumor, not immune cells. Using the Nuance software for calculation of number pixels of positive staining per analyzed image of MYC and HLA-B we observed a negative correlation (R = −0.6, P = 0.02) only for samples where MYC was detected (N = 15; Fig. 4C). For samples where MYC was barely detected (<1%), the correlation was not seen, which indicated that there are other factors regulating HLA-B expression in the absence of MYC. NF-κB p105 was detectable in the nuclei of both tumor cells and TILs and the levels of expression from these were positively correlated (P = 3 × 10−5). Importantly, mRNA expression of NFKB1 was positively correlated with tumor NF-κB p105 nuclear staining (P = 0.02; Fig. 4).

MYC was more frequently amplified, whereas NF-κB and IFNγ signaling genes were more frequently deleted in the low immune subgroup

Given that we observed the upregulation of MYC and downregulation of NFKB1 expression (the nodal genes) in the low immune subgroup, we hypothesized that MYC amplifications and NFKB1 deletions would be more common in this immune subgroup in the LMC. Using next-generation sequencing–derived copy number data from a subset of the LMC tumors, we observed that 29% had amplifications of MYC and 14% deletions of NFKB1 in the low immune subgroup, more than in the intermediate or the high immune subgroup (P = 0.02 for MYC, P = 0.0003 for NFKB1; Fig. 5A; Supplementary Data S2D). Interestingly, both of these copy number changes were strongly predictive of poor prognosis overall (adjusted for AJCC stage) separately [MYC amplifications: HR = 1.8 (95% CI 1.8–2.6), P = 0.006; NFKB1 deletions: HR = 1.5 (95% CI 1.1–2.1), P = 0.007) and when combined [HR = 3.7 (95% CI 1.6–8.5), P = 0.002, adjusted for AJCC; Fig. 5B; Supplementary Data S2E and S2F].

Because the NF-κB and IFNγ pathways were among the most enriched pathways in the high immune subgroup, we then asked whether other genes within these pathways were deleted in the low immune subgroup. Indeed, we found the evidence of deletion of NFKB2 (26% of whole dataset), CHUK (22%), MYD88 (5%), IRAK2 (5%), MAP3K7 (17%), JAK2 (10%), and STAT1 (4%). These copy number changes were not mutually exclusive (Fig. 5A), but were much more frequent in the low immune than in other subgroups (Fig. 5B; Supplementary Data S2D). Deletion of CHUK, MYD88, IRAK2, or JAK2 each were predictive of death from melanoma after adjusting for AJCC stage (Supplementary Data S2E and S2F). As expected, these copy number changes were highly correlated with mRNA expression of corresponding genes (Supplementary Data S2F).

In our previous study (11), we demonstrated that CTNNB1 expression was overexpressed in 30% of all primaries and 59% of in the CIC4 (immune low/β-catenin high). Here, we observed a similar pattern across the three immune subgroups: CTNNB1 was more commonly overexpressed in the low immune subgroup compared with other groups (Fig. 5A). Comparing the CNAs of genes in the NF-κB pathway and CTNNB1 expression, we found some overlap but also heterogeneity in the low immune subgroup. Specifically, 15% of tumors had evidence of increased CTNNB1 expression alone, 32% had a deletion in at least one gene of the NF-κB pathway in the absence of CTNNB1 overexpression, whereas 31% had both (i.e., increased β-catenin and a deletion in at least one gene). In prognostic terms, in the whole dataset, the HR for melanoma death in the presence of CTNNB1 upregulation was HR = 2.2, P = 5 × 10−5, 95% CI 1.5–3.1; for NF-κB pathway deletions was HR = 2.03, P = 2 × 10−4, 95% CI 1.4–3.0; and for the combination of these two pathways was HR = 3.4; P = 5 × 10−5; 95% CI, 2.2–5.5.

These data demonstrate the involvement of genetic factors in modulating the immunity and shaping the tumor immunophenotype. However, it is commonly postulated that environmental factors also play a role in this process, and we, therefore, tested this hypothesis.

Smoking as a strong independent risk factor for melanoma death in the high immune subgroup

In addition to clinico-pathologic tumor characteristics, the LMC has a record of patient smoking behaviors, vitamin D levels from a blood sample at diagnosis, and a deprivation index measured by the Townsend score (24). In a univariable Cox proportional hazard model, AJCC staging, mitotic number, site of primary melanoma, age at diagnosis, sex, and smoking (never/ever) were significantly predictive of MSS in the whole dataset, whereas season-adjusted vitamin D was not. Among these variables, AJCC stage, mitotic number, site of melanoma, age at diagnosis, and smoking remained significant in multivariable analysis of the whole dataset, but different sets of variables were significant in each of the three immune subgroups (Table 1). Body site of the primary melanoma was a strong predictor of MSS in the low immune subgroup along with AJCC stage, driven by tumors arising in the sun-protected body sites that were predominantly classified within this group and are known to have a particularly bad outcome (Table 1; ref. 28). The prognostic effect of smoking was heterogeneous across the three immune subgroups (P < 0.03 for equal HRs across the subgroups). In the high immune subgroup (HR = 4.6 for “ever smoked”, P = 0.003, N = 156), compared with the intermediate immune subgroup (HR = 1.8, P = 0.05, N = 275) and the low immune subgroup (HR = 0.9, P = 0.7, N = 272; Fig. 6). The deleterious effect of smoking in the high immune subgroup was reproduced when the analysis was repeated using two alternative definitions of smoking habits: duration of smoking (number of years) and the cumulative number of smoked cigarettes (packs per year; Supplementary Data S2G and S2H). The negative prognostic effect of cigarette smoking in the high immune subgroup remained significant after adjusting for the deprivation index (HR = 1.6, P = 0.001). To gain a deeper insight into the interplay between smoking and immune responses, we assessed the expression of GPR15, which has previously been described as a biomarker of exposure to tobacco smoke, with increased expression previously demonstrated in a number of immune cell types measured in the peripheral blood (29, 30). In the tumors, we found no significant association between GPR15 expression with smoking (never/ever) across the whole LMC dataset (fold change = 1.07, P = 0.12). However, the association was stronger in the high immune subgroup (fold change = 1.32, P = 0.02; Supplementary Fig. S4A). Furthermore, GPR15 expression was the highest in the high immune subgroup when testing a subset of data of ever smokers: P = 5 × 10−5, whereas the result was not significant for the second subset, never smokers: P = 0.3 (Supplementary Fig. S4B). GPR15 expression in the blood is reported to decrease after cessation of smoking (30) and therefore, we assessed its expression in “still smokers” compared with “non-smokers.” We observed a markedly stronger differential expression in the high immune subgroup for still smoking (fold change = 1.9, P = 0.002) than in the whole dataset (fold change = 1.32, P = 0.01; Supplementary Fig. S4B). We tested the differences in immune cell scores between never/ever smokers (in the high immune subgroup), but we did not find any statistically significant results (Supplementary Table S3). We also examined the association between smoking and tumor histologic features (Supplementary Table S4) as well as the whole-genome expression, including cytokine genes, but no significant associations were identified after multiple testing.

The dramatic survival benefit of checkpoint blockade in melanoma, in around half of patients with advanced disease (31), has highlighted the need to understand the drivers of low immune (“cold” tumors), which are less likely to respond to treatment. In silico immune cell analysis in cancer has been adopted in recent years to better understand these drivers (32, 12, 33); although these approaches do not allow distinction to be made between weak signals coming from numerous infiltrating immune cells and strong signals from fewer cells.

Previously, we reported the survival analysis of 24 immune cell scores (11), derived from an earlier version of the immunome (10). In the current analysis, we used an updated version of the immunome (12) allowing inference of 31 immune cell scores, which we reduced to 27 after the gene filtration. We noted that in the earlier report, 10 of 24 cells were significantly prognostic and replicated in TCGA (41% of cells), but in this study, the number increased to 17 of 27 (63%).

In evaluating the updated immunome, we showed that almost half of the genes proposed to be specific to particular immune cell subsets were also expressed at significant levels in melanoma cell lines, disqualifying them as immune specific. The fact that not all of the genes postulated to characterize a particular immune cell type were positively correlated may represent in part a technical feature of our dataset, as we observed that the correlations between these genes were slightly higher in the TCGA dataset sourced from fresh-frozen tumors rather than archived formalin-fixed paraffin-embedded specimens. These observations suggest that use of “off the shelf” algorithms to infer immune activity may have limited application.

We have identified three immune subgroups with distinct survival profiles indicating better survival in the presence of stronger immune responses. We also showed that these subgroups were stronger in terms of prognosis prediction than the three-class identified in the metastatic melanoma from TCGA, which we applied to our primary tumors from LMC. We recapitulated three immune subgroups in the TCGA metastatic melanomas, suggesting that similar immune infiltration and exclusion mechanisms span the whole spectrum of disease progression.

All the immune scores were highly correlated with each other (including those known to play an immunosuppressive role), the majority being upregulated in the high immune subgroup. We did not observe increasing representation of immunosuppressive cells (e.g., Tregs) or a relative increase in the expression of checkpoint molecules in the low immune subgroup. Rather, our data suggest a coordination of the immune cell populations as a whole. This is entirely in keeping with previously published observations of increased Treg numbers and accompanying expression of the checkpoint molecules as a result of homeostatic mechanisms driven by melanoma-infiltrating CD8+ T cells (34). We cannot, however, exclude the possibility that the inference of immune cell subgroup infiltration from transcriptomic datasets may be insensitive to subtle variations that nevertheless might have an impact on immune function.

The protein–protein network analysis in Reactome FIViz of genes upregulated in the low immune subgroup revealed enrichment for the genes in cell proliferative pathways with MYC as the major node (the gene with the highest centrality). MYC is a pro-proliferative oncogene, which has in recent years been reported to have various immunosuppressive functions (27, 35, 36, 37), and to have specific involvement in melanoma metastasis and invasiveness (38). However, the relation of MYC and immune response within melanoma is unclear.

In our study, we were able to show evidence for a negative relationship between MYC and antigen processing and presentation machinery especially with HLA-B in tumors and patient-derived melanoma cell lines. An inverse relationship between HLA class I genes and MYC expression has previously been reported (39) in melanoma cell lines. Moreover, it was described that MYC downregulates the expression of HLA-B by directly binding to its proximal promoter (40). Our data, therefore, provide strong evidence that MYC contributes significantly to the immune evasion in primary melanoma making it a therapeutic target, notwithstanding the difficulties of achieving this (41). The requirement for MYC in T cells suggests that a more targeted approach may be required, or that regulators/effectors of MYC activity might prove more appropriate targets (41).

NFKB1 was the network hub gene in the high immune subgroup. A number of important NF-κB family genes (NFKB1, NFKB2, c-REL, RELB) were also upregulated in this group, suggesting activation of the pathway. RELA was stable across the immune subgroups, reflecting its constitutive expression in different tissues. IHC staining showed that tumor and TILs nuclear localization of NF-κB significantly correlated with lymphocytic infiltration, suggesting a reciprocal NF-κB–driven phenotype generated between the tumor and its immediate microenvironment, as described in other cancers (42, 43). Conversely, in the low immune subgroup, we found the loss of genes important in NF-κB and IFNγ signaling, resulting in decreased gene expression. JAK2 mutations have recently been reported to be involved in acquired and primary resistance to anti-PD-1 therapy (44, 45). Our hypothesis, therefore, is that a significant proportion of melanoma tumors in the low immune subgroup may have primary resistance to this therapy even in adjuvant usage.

In our study, we report for the first time the association of smoking with immune responses to primary melanoma. Our results implied that smoking had an adverse effect on outcome by reducing the protective value of immune infiltration. That there were no obvious transcriptomic differences between the melanomas in smokers and nonsmokers may, however, suggest that the immune infiltrate in smokers may simply represent nonspecific systemic inflammation or even that we see similar transcriptomic signals from protumorigenic (cigarette driven) and antitumor immune responses.

We did observe a positive correlation between smoking and the expression of the GPR15 gene, which codes for a chemo-attractant receptor that is regarded as a biomarker of smoking known to be hypomethylated and hence, overexpressed in the circulating immune cells in smokers (46, 29). The GPR15 protein is reported to play a role in the trafficking of T cells (47, 48), but its full biological function and significance with respect to smoking is still unknown. The overall pattern of association between reported smoking and death from melanoma, however, reinforces the view that discontinuation of smoking should be strongly recommended in patients with melanoma. As it is not known whether the adverse effects of smoking in melanoma are mediated by nicotine or other components of cigarettes, the recommendation should probably be to avoid vaping (49), despite the acknowledged difficulty of smoking cessation for many.

In conclusion, we report the use of bioinformatics to define broad prognostic immunophenotypes of primary melanoma, with evidence of a prominent role of NF-κB and IFNγ signaling downregulation (including by deletion) and MYC overexpression (including amplification) in driving immunosuppression. We report evidence that a key mechanism in this process is perturbation of the antigen presentation machinery and that smoking predicted significantly worse MSS in patients with evidence for stronger immune responses.

S.J. O’Shea reports receiving a commercial research grant from Cancer Research UK, has received speakers bureau honoraria from L’Oreal UK Ltd, and has provided expert testimony for Janssen-Cilag Ltd and Pierre Fabre Ltd. A.P. Droop is the director at AD Bioinformatics Limited. M.P. Levesque reports receiving commercial research grants from Roche and BMS. J. Newton-Bishop is the expert witness at Legal Reports. No potential conflicts of interest were disclosed by the other authors.

Conception and design: J. Nsengimana, A.P. Droop, D.T. Bishop, J. Newton-Bishop

Development of methodology: J. Poźniak, A.P. Droop, J.A. Randerson-Moor, D.T. Bishop, J. Newton-Bishop

Acquisition of data (provided animals, acquired and managed patients, provided facilities, etc.): J. Nsengimana, J.P. Laye, S.J. O’Shea, A. Filia, M. Harland, T. Mell, S.A. Hogan, S.N. Freiberger, J. Newton-Bishop

Analysis and interpretation of data (e.g., statistical analysis, biostatistics, computational analysis): J. Poźniak, J. Nsengimana, J.M.S. Diaz, A.P. Droop, J.R. Davies, S. Muralidhar, S.A. Hogan, M.P. Levesque, G.P. Cook, D.T. Bishop, J. Newton-Bishop

Writing, review, and/or revision of the manuscript: J. Poźniak, J. Nsengimana, J.P. Laye, S.J. O’Shea, J.M.S. Diaz, A. Filia, M. Harland, J.R. Davies, S.N. Freiberger, M.P. Levesque, G.P. Cook, D.T. Bishop, J. Newton-Bishop

Administrative, technical, or material support (i.e., reporting or organizing data, constructing databases): T. Mell, J.A. Randerson-Moor, J. Newton-Bishop

Study supervision: J. Nsengimana, J. Newton-Bishop

Other (performed IHC staining and scoring): J. Poźniak

Other (PhD supervision): J. Newton-Bishop

We are thankful to Professor Ulf Klein for providing knowledge and insight into NF-κB signaling. We are grateful to all the participants whose samples have been used in this study. This work was funded by Cancer Research UK C588/A19167, C8216/A6129, and C588/A10721 and NIH CA83115. J. Poźniak, J.M.S. Diaz, and S. Muralidhar were funded by Horizon 2020 Research and Innovation Programme no. 641458 (MELGEN).

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.

1.
Clark
WH
,
Elder
DE
,
Guerry
D
,
Braitman
LE
,
Trock
BJ
,
Schultz
D
, et al
Model predicting survival in stage I melanoma based on tumor progression
.
J Natl Cancer Inst
1989
;
81
:
1893
904
.
2.
Day
CL
,
Sober
AJ
,
Kopf
AW
,
Lew
RA
,
Mihm
MC
,
Golomb
FM
, et al
A prognostic model for clinical stage I melanoma of the trunk
.
Am J Surg
1981
;
142
:
247
51
.
3.
Larkin
J
,
Chiarion-Sileni
V
,
Gonzalez
R
,
Grob
JJ
,
Cowey
CL
,
Lao
CD
, et al
Combined nivolumab and ipilimumab or monotherapy in untreated melanoma
.
N Engl J Med
2015
;
373
:
23
34
.
4.
Hugo
W
,
Zaretsky
JM
,
Sun
L
,
Song
C
,
Moreno
BH
,
Hu-Lieskovan
S
, et al
Genomic and transcriptomic features of response to anti-PD-1 therapy in metastatic melanoma
.
Cell
2017
;
168
:
542
.
5.
Snyder
A
,
Makarov
V
,
Merghoub
T
,
Yuan
J
,
Zaretsky
JM
,
Desrichard
A
, et al
Genetic basis for clinical response to CTLA-4 blockade in melanoma
.
N Engl J Med
2014
;
371
:
2189
99
.
6.
Thorsson
V
,
Gibbs
DL
,
Brown
SD
,
Wolf
D
,
Bortone
DS
,
Ou Yang
TH
, et al
The immune landscape of cancer
.
Immunity
2018
;
48
:
812
30
.
7.
Charoentong
P
,
Finotello
F
,
Angelova
M
,
Mayer
C
. 
Pan-cancer immunogenomic analyses reveal genotype-immunophenotype relationships and predictors of response to checkpoint blockade
.
BioRxiv
2016
Available from:
https://www.biorxiv.org/content/10.1101/056101v3.
8.
Newton-Bishop
JA
,
Beswick
S
,
Randerson-Moor
J
,
Chang
YM
,
Affleck
P
,
Elliott
F
, et al
Serum 25-hydroxyvitamin D3 levels are associated with Breslow thickness at presentation and survival from melanoma
.
J Clin Oncol
2009
;
27
:
5439
44
.
9.
Newton-Bishop
JA
,
Davies
JR
,
Latheef
F
,
Randerson-Moor
J
,
Chan
M
,
Gascoyne
J
, et al
25-Hydroxyvitamin D2/D3 levels and factors associated with systemic inflammation and melanoma survival in the Leeds Melanoma Cohort
.
Int J Cancer
2015
;
136
:
2890
9
.
10.
Bindea
G
,
Mlecnik
B
,
Tosolini
M
,
Kirilovsky
A
,
Waldner
M
,
Obenauf
AC
, et al
Spatiotemporal dynamics of intratumoral immune cells reveal the immune landscape in human cancer
.
Immunity
2013
;
39
:
782
95
.
11.
Nsengimana
J
,
Laye
J
,
Filia
A
,
O’Shea
S
,
Muralidhar
S
,
Poźniak
J
, et al
β-Catenin–mediated immune evasion pathway frequently operates in primary cutaneous melanomas
.
J Clin Invest
2018
;
128
:
2048
63
.
12.
Angelova
M
,
Charoentong
P
,
Hackl
H
,
Fischer
ML
,
Snajder
R
,
Krogsdam
AM
, et al
Characterization of the immunophenotypes and antigenomes of colorectal cancers reveals distinct tumor escape mechanisms and novel targets for immunotherapy
.
Genome Biol
2015
;
16
:
64
.
13.
Chen
DS
,
Mellman
I
. 
Elements of cancer immunity and the cancer-immune set point
.
Nature
2017
;
541
:
321
30
.
14.
Jewell
R
,
Conway
C
,
Mitra
A
,
Randerson-Moor
J
,
Lobo
S
,
Nsengimana
J
, et al
Patterns of expression of DNA repair genes and relapse from melanoma
.
Clin Cancer Res
2010
;
16
:
5211
21
.
15.
Nsengimana
J
,
Laye
J
,
Filia
A
,
Walker
C
,
Jewell
R
,
Van den Oord
JJ
, et al
Independent replication of a melanoma subtype gene signature and evaluation of its prognostic value and biological correlates in a population cohort
.
Oncotarget
2015
;
6
:
11683
93
.
16.
Monti
S
,
Tamayo
P
,
Mesirov
J
,
Golub
T
. 
Consensus clustering: A resampling based method for class discovery and visualization of gene expression microarray data
.
Mach Learn
2003
;
52
:
91
118
.
17.
Wilkerson
MD
,
Hayes
DN
. 
ConsensusClusterPlus: a class discovery tool with confidence assessments and item tracking
.
Bioinformatics
2010
;
26
:
1572
3
.
18.
Wu
G
,
Feng
X
,
Stein
L
,
Hanahan
D
,
Weinberg
R
,
Vogelstein
B
, et al
A human functional protein interaction network and its application to cancer data analysis
.
Genome Biol
2010
;
11
:
R53
.
19.
Scardoni
G
,
Petterlini
M
,
Laudanna
C
. 
Analyzing biological network parameters with CentiScaPe
.
Bioinformatics
2009
;
25
:
2857
9
.
20.
Shannon
P
. 
Cytoscape: a software environment for integrated models of biomolecular interaction networks
.
Genome Res
2003
;
13
:
2498
504
.
21.
Bastian
M
,
Heymann
S
,
Jacomy
M
. 
Gephi: an open source software for exploring and manipulating networks
.
Third Int AAAI Conf Weblogs Soc Media
2009
Available from:
https://www.aaai.org/ocs/index.php/ICWSM/09/paper/view/154.
22.
Raaijmakers
MIG
,
Widmer
DS
,
Maudrich
M
,
Koch
T
,
Langer
A
,
Flace
A
, et al
A new live-cell biobank workflow efficiently recovers heterogeneous melanoma cells from native biopsies
.
Exp Dermatol
2015
;
24
:
377
80
.
23.
Gu
Z
,
Eils
R
,
Schlesner
M
. 
Complex heatmaps reveal patterns and correlations in multidimensional genomic data
.
Bioinformatics
2016
;
32
:
2847
9
.
24.
Townsend
P
. 
Deprivation
.
Jnl Soc Pol
1987
;
16
:
125
46
.
25.
Masten
BJ
,
Olson
GK
,
Tarleton
CA
,
Rund
C
,
Schuyler
M
,
Mehran
R
, et al
Characterization of Myeloid and Plasmacytoid Dendritic Cells in Human Lung
.
J Immunol
2006
;
177
:
7784
93
.
26.
Fening
K
,
Parekh
V
,
McKay
K
. 
CD123 immunohistochemistry for plasmacytoid dendritic cells is useful in the diagnosis of scarring alopecia
.
J Cutan Pathol
2016
;
43
:
643
8
.
27.
Casey
SC
,
Tong
L
,
Li
Y
,
Do
R
,
Walz
S
,
Fitzgerald
KN
, et al
MYC regulates the antitumor immune response through CD47 and PD-L1
.
Science
2016
;
352
:
227
31
.
28.
Asgari
MM
,
Shen
L
,
Sokil
MM
,
Yeh
I
,
Jorgenson
E
. 
Prognostic factors and survival in acral lentiginous melanoma
.
Br J Dermatol
2017
;
177
:
428
35
.
29.
Bauer
M
,
Fink
B
,
Thürmann
L
,
Eszlinger
M
,
Herberth
G
,
Lehmann
I
. 
Tobacco smoking differently influences cell types of the innate and adaptive immune system—indications from CpG site methylation
.
Clin Epigenetics
2016
;
7
:
83
.
30.
Kõks
S
,
Kõks
G
. 
Activation of GPR15 and its involvement in the biological effects of smoking
.
Exp Biol Med
2017
;
242
:
1207
12
.
31.
Wolchok
JD
,
Chiarion-Sileni
V
,
Gonzalez
R
,
Rutkowski
P
,
Grob
J-J
,
Cowey
CL
, et al
Overall survival with combined nivolumab and ipilimumab in advanced melanoma
.
N Engl J Med
2017
;
377
:
1345
56
.
32.
Newman
AM
,
Liu
CL
,
Green
MR
,
Gentles
AJ
,
Feng
W
,
Xu
Y
, et al
Robust enumeration of cell subsets from tissue expression profiles
.
Nat Methods
2015
;
12
:
453
7
.
33.
Becht
E
,
Giraldo
NA
,
Lacroix
L
,
Buttard
B
,
Elarouci
N
,
Petitprez
F
, et al
Estimating the population abundance of tissue-infiltrating immune and stromal cell populations using gene expression
.
Genome Biol
2016
;
17
:
218
.
34.
Spranger
S
,
Spaapen
RM
,
Zha
Y
,
Williams
J
,
Meng
Y
,
Ha
TT
, et al
Up-Regulation of PD-L1, IDO, and tregs in the melanoma tumor microenvironment is driven by CD8+ T cells
.
Sci Transl Med
2013
;
5
:
200ra116
.
35.
Kortlever
RM
,
Sodir
NM
,
Wilson
CH
,
Burkhart
DL
,
Pellegrinet
L
,
Brown Swigart
L
, et al
Myc cooperates with ras by programming inflammation and immune suppression
.
Cell
2017
;
171
:
1301
15
.
36.
Topper
MJ
,
Vaz
M
,
Chiappinelli
KB
,
DeStefano Shields
CE
,
Niknafs
N
,
Yen
R-WC
, et al
Epigenetic Therapy Ties MYC depletion to reversing immune evasion and treating lung cancer
.
Cell
2017
;
171
:
1284
1300
.
37.
Schaub
FX
,
Dhankani
V
,
Berger
AC
,
Trivedi
M
,
Richardson
AB
,
Shaw
R
, et al
Pan-cancer Alterations of the MYC oncogene and its proximal network across the cancer genome atlas
.
Cell Syst
2018
;
6
:
282
300
.
38.
Schlagbauer-Wadl
H
,
Griffioen
M
,
van Elsas
A
,
Schrier
PI
,
Pustelnik
T
,
Eichler
HG
, et al
Influence of increased c-Myc expression on the growth characteristics of human melanoma
.
J Invest Dermatol
1999
;
112
:
332
6
.
39.
Versteeg
R
,
Noordermeer
IA
,
Krüse-Wolters
M
,
Ruiter
DJ
,
Schrier
PI
. 
c-myc down-regulates class I HLA expression in human melanomas
.
EMBO J
1988
;
7
:
1023
9
.
40.
Peltenburg
LTC
,
Schrier
PI
. 
Transcriptional suppression of HLA-B expression by c-Myc is mediated through the core promoter elements
.
Immunogenetics
1994
;
40
:
54
61
.
41.
Chen
H
,
Liu
H
,
Qing
G
. 
Targeting oncogenic Myc as a strategy for cancer treatment
.
Signal Transduct Target Ther
2018
;
3
:
5
.
42.
Hopewell
EL
,
Zhao
W
,
Fulp
WJ
,
Bronk
CC
,
Lopez
AS
,
Massengill
M
, et al
Lung tumor NF-kB signaling promotes T cell–mediated immune surveillance
.
J Clin Invest
2013
;
123
:
2509
22
.
43.
Muthuswamy
R
,
Berk
E
,
Junecko
BF
,
Zeh
HJ
,
Zureikat
AH
,
Normolle
D
, et al
NF- B hyperactivation in tumor tissues allows tumor-selective reprogramming of the chemokine microenvironment to enhance the recruitment of cytolytic T effector cells
.
Cancer Res
2012
;
72
:
3735
43
.
44.
Zaretsky
JM
,
Garcia-Diaz
A
,
Shin
DS
,
Escuin-Ordinas
H
,
Hugo
W
,
Hu-Lieskovan
S
, et al
Mutations associated with acquired resistance to PD-1 blockade in melanoma
.
N Engl J Med
2016
;
375
:
819
29
.
45.
Shin
DS
,
Zaretsky
JM
,
Escuin-Ordinas
H
,
Garcia-Diaz
A
,
Hu-Lieskovan
S
,
Kalbasi
A
, et al
Primary resistance to PD-1 blockade mediated by JAK1/2 mutations
.
Cancer Discov
2017
;
7
:
188
201
.
46.
Bauer
M
,
Linsel
G
,
Fink
B
,
Offenberg
K
,
Hahn
AM
,
Sack
U
, et al
A varying T cell subtype explains apparent tobacco smoking induced single CpG hypomethylation in whole blood
.
Clin Epigenetics
2015
;
7
:
81
.
47.
Bilsborough
J
,
Viney
JL
. 
GPR15: a tale of two species
.
Nat Immunol
2015
;
16
:
137
9
.
48.
Lahl
K
,
Sweere
J
,
Pan
J
,
Butcher
E
. 
Orphan chemoattractant receptor GPR15 mediates dendritic epidermal T-cell recruitment to the skin
.
Eur J Immunol
2014
;
44
:
2577
81
.
49.
Martin
EM
,
Clapp
PW
,
Rebuli
ME
,
Pawlak
EA
,
Glista-Baker
E
,
Benowitz
NL
, et al
E-cigarette use results in suppression of immune and inflammatory-response genes in nasal epithelial cells similar to cigarette smoke
.
Am J Physiol Cell Mol Physiol
2016
;
311
:
L135
44
.

Supplementary data