Genetic and Environmental Determinants of Immune Response to Cutaneous Melanoma

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 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).

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 log₂-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.

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).

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.

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⁻⁶, 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.

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.

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).

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.

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, χ² 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.

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).

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.

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⁻⁷ (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⁻⁵ (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).

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 (χ², P = 4.0 × 10⁻⁷) and a single observer from our research group who was blinded to the transcriptomic data (χ², P = 3.6 × 10⁻⁸; see Supplementary Table S2). The mitotic number was significantly lower in the high immune subgroup (Kruskal–Wallis, P = 2 × 10⁻⁴). 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%;χ² P = 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).

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⁻¹⁴–10⁻⁷) 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⁻¹⁴–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⁻¹⁰; 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⁻¹⁴).

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⁻⁶, respectively).

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⁻⁵). Importantly, mRNA expression of NFKB1 was positively correlated with tumor NF-κB p105 nuclear staining (P = 0.02; Fig. 4).

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⁻⁵, 95% CI 1.5–3.1; for NF-κB pathway deletions was HR = 2.03, P = 2 × 10⁻⁴, 95% CI 1.4–3.0; and for the combination of these two pathways was HR = 3.4; P = 5 × 10⁻⁵; 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.

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⁻⁵, 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.