Transcriptomic Analysis Reveals Prognostic Molecular Signatures of Stage I Melanoma

Cutaneous melanoma continues to increase in incidence worldwide. Although earlier diagnosis has been documented with correspondingly better outcomes, the rising incidence of thinner tumors means that, counterintuitively, one-fifth of deaths now occur in patients presenting initially with early disease (1). In the United Kingdom, 91% of melanomas are diagnosed at American Joint Committee on Cancer (AJCC) stages I to II (2). Therefore, better prognostic biomarkers are needed to identify early-stage disease requiring adjuvant therapies, as well as predictive biomarkers of response to checkpoint blockade.

Previous transcriptomic analyses of cutaneous melanoma have generated gene signatures with a prognostic value independent of AJCC stage (3–7). The prognostic signature developed by Jonsson and colleagues (3) classifies metastatic melanomas into four classes (Lund 4 classes), later simplified into two classes (Lund 2 grades; ref. 4), and the signature developed by The Cancer Genome Atlas (TCGA) consortium classified melanomas into three classes (TCGA 3 classes; ref. 8). The prognostic significance of the Lund 4 class and TCGA 3 class signatures have been replicated in relatively small datasets, notably with few AJCC stage I patients (5, 9). Another transcriptomic signature based on 27 genes was developed by Gerami and colleagues (6) to classify patients with primary melanoma into tumors that were high or low risk for metastasis.

In this study, the first aim was to test the prognostic value of the Lund and TCGA signatures, as well as the gene list of signature of Gerami and colleagues (6) in a large population-based cohort of primary melanomas with a good proportion of stage I patients and extensive phenotypic annotations (Leeds Melanoma Cohort, LMC). Because the dataset was well powered for discovery of novel tumor subtypes, unsupervised clustering of the tumor transcriptomes of the LMC was performed and the prognostic value of the resultant signature was compared with that of sentinel node biopsy (SNB) in analyses stratified by AJCC stage. Finally, the association between the Leeds signature and outcome was tested in published data from patients receiving immunotherapy (10–12).

As described previously (13), 2,184 patients with primary melanoma were recruited to the LMC in the period of 2000 to 2012. This was a population-ascertained cohort, which, therefore, recruited patients treated at multiple clinical centers (recruitment rate 67%). During this period, SNB was neither offered nor accepted universally. The study was ethically approved (ethical approval MREC 1/3/57, PIAG 3-09(d)/2003) and in accordance with the Declaration of Helsinki. Participants were consented to sampling of their formalin-fixed, paraffin-embedded tumor blocks that were stored in the NHS (UK National Health Service) histopathology departments of the respective hospitals. Hematoxylin and eosin–stained slides were generated and examined to facilitate subsequent sampling of the blocks using a 0.6-mm diameter tissue microarray needle as previously reported (5, 13). Prior to sampling, all the tumor blocks were reviewed, and if there was only a small amount of tumor left in the block, then the block was not sampled, lest a clinically important block be destroyed. Up to two cores were sampled from each block, and, to increase the comparability between tumors, the samples were consistently extracted from the least inflamed, least stromal regions of the invasive front of the tumor. The tumor-infiltrating lymphocytes (TILs) were scored using the classification system of Clark and colleagues (14). As previously reported (13), 703 tumor transcriptomes were profiled and in this study, 16 samples were removed in quality control leaving a cohort of 687 patients, henceforth referred to as the whole LMC dataset. The dataset contained 251 patients who had an SNB test (Supplementary Table S10), and only 16 patients are known so far to have been treated with checkpoint blockade. The LMC patients were assigned an AJCC stage based on the AJCC staging eighth edition (15). Where patients did not have an SNB, the AJCC staging used was derived from clinical staging and pathologic examination of the wide local excision sample.

Both mRNA and DNA were extracted from the tumor samples derived from cores following a previously described protocol (5, 13). The whole genome mRNA expression profiling was carried out using Illumina's DASL-HT12-v4 array. As described previously, for quality control, the mRNA was extracted from up to two cores for a number of patients (117 duplicates in total); gene expression data from only one extraction per patient was used in subsequent analyses (13). The raw transcriptomic data were extracted from the image files using GenomeStudio (Illumina Inc.) and were preprocessed as previously reported (13). Briefly, after background correction and quantile normalization (R package LUMI; ref. 16), singular value decomposition was used to remove the batch effect [R package SWAMP (refs. 13, 17)].

The array included 29,262 probes corresponding to 20,715 unique genes. For genes with multiple probes, the probe detected in the largest number of tumors was retained, and two additional filters were applied: genes had to be detected with P < 0.05 in at least 40% of tumors and had to have a standard deviation (SD) >0.40. This SD threshold was chosen on the basis of the overall distribution across the 20,715 genes on the log2 scale. The median SD was 0.68. The data were standardized to give each gene a mean of 0 and SD of 1.

The LMC tumors were classified into the Lund 4 classes, Lund 2 grades, and TCGA 3 classes using the supervised nearest centroid classification (NCC) as described previously (5). All the 27 genes of the gene signature of Gerami and colleagues (6) were present in LMC dataset and were analyzed using a univariable survival analysis in the whole LMC dataset and stage I tumors. Unsupervised clustering was performed using the consensus Partitioning Around Medoids clustering method in the R-package ConsensusClusterPlus (18, 19) with Euclidean distance as the dissimilarity measure and a resampling fraction of 0.8 for both genes and samples in 1,000 iterations (Supplementary Methods).

Cox proportional hazard models and Kaplan–Meier curves were used to test the association with survival (R-package Survival; ref. 20). The survival time was calculated from time of diagnosis to the time of last follow-up or time of death from melanoma, whichever occurred first, referred to as melanoma-specific survival (MSS). Patients with deaths caused by factors other than melanoma were censored at the time of death. ROC analysis was performed using AJCC stage pre-SNB and AJCC stage post-SNB for patients who had SNB. Clinical staging prior to SNB is described as AJCC pre-SNB. The AJCC stage post-SNB includes additional information on regional lymph node metastasis. The analysis used AJCC staging eighth edition, and MSS up to 3 years was chosen as cutoff based upon the inclusion of the majority of the deaths without loss of data as a result of censoring (Supplementary Table S11). Patients who were censored before 3 years were not included in the analysis. The analysis was performed using R-packages pROC, plotROC, and ggplot2 (21–23).

The differentially expressed genes (DEG) were identified using the Significance Analysis of Microarrays (R-package SAMR) by comparing each class versus all others (24). Pathway enrichment and biological network analysis of DEGs with a q value equal to 0 were performed using ReactomeFiviz in Cytoscape (25). The central nodes of the biological network were identified using a centrality measure (betweenness) in Gephi (ref. 26; Supplementary Methods).

The copy-number alterations (CNA) data were generated in a subset of LMC tumors using Illumina's next-generation sequencing (NGS) platform as reported in the study by Filia and colleagues (ref. 27; Supplementary Methods). Among the 687 transcriptome-profiled patients of LMC, CNA data were available for 272 patients. The CNA were assessed in the regions spanning the genes identified as hubs in network enrichment analysis. The ratio between mean of the window read counts in the region mapping to a gene and the average read count of the 10 flanking regions around that gene was used to estimate the copy-number changes. The windows (5k) corresponding to a gene locus were identified using the R packages biomaRt (28, 29). The cutoff for calling a region amplified was chosen as a value greater than 0.4 whereas a value less than −0.4 was used to identify a deletion. The 272 samples in the CNA dataset were at AJCC stages I (n = 80), II (n = 147), and III (n = 45; similar distribution to the whole LMC dataset).

For replication, a primary melanoma transcriptomic dataset of 223 tumors from a Lund cohort (Sweden) was used (Harbst and colleagues; ref. 4). The samples were classified using the newly generated signature by the supervised NCC approach (5). Of those 223 patients, 200 had recorded information on melanoma relapse in the follow-up time postdiagnosis and were used to test the association between patient subgroups and relapse-free survival (using Cox proportional hazard models, Kaplan–Meier curves, and log-rank test).

Three publicly available transcriptome datasets (Hugo Cohort: GSE78220, Ulloa-Monotoya Cohort: GSE35640, and Riaz Cohort: https://github.com/riazn/bms038_analysis) were downloaded (10–12) and samples were quantile normalized and classified using the NCC method (Pearson correlation coefficient). The Riaz cohort was a mixture of samples from various melanoma types (cutaneous melanoma, mucosal melanoma, acral melanoma, uveal/ocular melanoma, others). In this study, the samples labeled as cutaneous melanoma were analyzed. In all the three cohorts, the association with response to immunotherapy was tested using Fisher exact test.

The structure of datasets used in this study is depicted in Fig. 1. When applied to the whole LMC dataset (n = 687), the three formerly published signatures (Lund 4 class, Lund 2 grade, and TCGA 3 class) replicated previously observed associations with MSS (Fig. 2A, C, and E). However, upon stratifying LMC patients on the basis of AJCC stage, the Lund and TCGA signatures showed no association with prognosis for LMC stage I patients (Fig. 2B, D, and F). The Lund 2-grade signature had the highest statistical power (since based on only two groups) and showed a statistically significant interaction with AJCC stage (P = 0.04, Supplementary Table S1), suggesting that the lack of association in stage I was not solely due to low sample size. Because the full details of the commercial signature of Gerami and colleagues (6) were not published, we were limited in the scope of its replication in the LMC dataset. However, analyzing the 27 Gerami genes identified 23 genes as predictors of prognosis in the whole LMC dataset (Supplementary Table S2). However, in keeping with the Lund and TCGA signatures, none of these genes showed a significant association with prognosis in stage I tumors (Supplementary Table S2).

Consensus clustering of the LMC dataset was performed, and following additional quality control measures (Supplementary Table S3), six distinct, novel molecular classes were identified (Fig. 3A). These classes were associated with clinicopathologic variables known to have prognostic value, including tumor site (P = 0.03), age at diagnosis (P = 0.03), mitotic rate (P = 0.002), ulceration (P = 0.01), AJCC stage (P = 6 × 10⁻¹⁰), TILs (P = 6 × 10⁻⁴), and Breslow thickness (P = 9 × 10⁻¹⁴; Table 1). The LMC class 1 and 5 tumors tended to be thin and nonulcerated, whereas class 2 and 4 tumors were thicker. Class 3 and 6 tumors were the thickest and most frequently ulcerated. The six classes showed strong association with BRAF (P = 6 × 10⁻⁵) and NRAS mutation status (P = 3 × 10⁻⁴): class 1, 5, and 6 tumors were frequently BRAF mutated, whereas class 2, 3, and 4 tumors were frequently NRAS mutated (Table 1).

The LMC classes predicted MSS in the whole LMC dataset and notably, across AJCC stages I, II and III subsets (Fig. 3B and C, Supplementary Fig. S1). In the unadjusted analysis of the whole dataset (Fig. 3B, Supplementary Table S4), class 1 (baseline) had the best prognosis and class 2 [HR = 1.7, 95% confidence interval (CI): 0.8–3.5] and class 5 (HR = 1.5, 95% CI: 0.7–3.1) showed intermediate prognosis, whereas class 3 (HR = 5.0, 95% CI: 2.5–10.1), class 4 (HR = 2.4, 95% CI: 1.2–4.7), and class 6 (HR = 3.1, 95% CI: 1.6–6.1) had the worst prognosis. In multivariable analysis, classes 3, 4, and 6 remained significant predictors of poor prognosis after including AJCC stage, sex, age at diagnosis, mitotic rate (Supplementary Table S4), and when the AJCC stage was replaced by ulceration and Breslow thickness in the model (Supplementary Table S6). In the LMC stage I subset, class 6 (HR = 6.6, 95% CI: 1.4–31.2) significantly predicted poor prognosis in unadjusted analysis (Fig. 3C and Supplementary Table S5) and it remained significant when sex, age at diagnosis, mitotic rate, ulceration, and Breslow thickness were adjusted (HR = 9.8, 95% CI: 1.1–86.2, Supplementary Table S6). Because Gerami signature was not available to us in full, we ran unsupervised clustering of the LMC dataset using the 27 Gerami genes to generate the two tumor groups analyzed by Gerami and colleagues (6), referred to as the Gerami clusters. This analysis showed that the LMC classes and Gerami clusters had independent prognostic effects in the whole LMC dataset (Supplementary Table S7); however, the Gerami clusters showed no prognostic value in stage I tumors whereas LMC class 6 remained a significant predictor in the multivariable model (Supplementary Table S8).

To validate the prognostic value of the LMC classes in an independent dataset, a 150-gene–based signature (LMC signature), generated after refining approximately 13,000 genes (Supplementary Fig. S2), was applied to the Lund dataset (4). In keeping with the observations made in the LMC dataset, class 3, class 4, and class 6 predicted worse prognosis in the Lund dataset, whereas class 1, class 2, and class 5 predicted better prognosis (Fig. 3D, Supplementary Table S9). Because the Lund dataset had only a few stage I cases (n = 58), the prognostic value of LMC signature could not be replicated in stage I disease.

In the dataset derived from individuals who had an SNB, the prognostic value of combined LMC class signature and pre-SNB AJCC stage was similar to that of AJCC stage with SNB (i.e. stage post-SNB; AUC 0.82 vs. 0.79, P = 0.7, Fig. 3E). Combining the LMC signature with AJCC stage post-SNB, patient's sex, age at diagnosis, and site of tumor increased the AUC to 0.88. Similarly, in the subset of patients at stage I pre-SNB, the LMC signature alone had comparable prognostic value to AJCC stage post-SNB (AUC = 0.88 vs. 0.82, P = 0.1, Fig. 3F). In this stage I subset, addition of stage post-SNB, patient's sex, age at diagnosis, and site of tumor to the LMC signature further increased the AUC to 0.99. However, the limited sample size of stage I dataset and including so many variables clearly overfitted the model, giving near-perfect classification and illustrating that independent datasets are needed to better assess performance.

The six classes of LMC signature showed distinct gene expression profiles (Fig. 4A) and showed partial overlap with the existing Lund and TCGA signatures. LMC classes 1, 3, and 5 overlapped substantially with the high-immune, pigmentation, and normal-like classes of the Lund 4 classes (Fig. 4B), and with the immune, MITF low, and keratin classes of the TCGA 3-classes (Fig. 4C). In contrast, LMC classes 2, 4, and 6 represented a mixture of the Lund 4 classes and TCGA 3 classes. Gene expression pathway enrichment analysis revealed distinctive biological features of the 6 LMC classes: notably class 2 was characterized by increased WNT signaling genes and metabolic pathways; class 4 by decreased expression of immune genes, and class 6 by increased expression of cell cycle and consistent downregulation of cell metabolism pathway genes (Supplementary Table S14).

When applied to the LMC 6 classes, the Lund modules (30) revealed discrimination consistent with enriched gene pathways: LMC class 1 tumors showed higher immune module activity, and class 3 tumors showed higher cell-cycle module activity (Fig. 4D). Interestingly, class 6 tumors had relatively higher cell cycle but also immune module activity and, as expected, the immune, stroma, and interferon modules were positively correlated but they negatively correlated with cell-cycle and MITF modules (Fig. 4D). The tumor-infiltrating immune cell populations imputed for each of the LMC classes (31) were consistent with the Lund immune module, as class 1 had the highest immune cell populations and class 3 the lowest, whereas class 6 appeared to maintain an intermediate level of immune cell populations, having the second highest scores on average (Supplementary Fig. S3).

A comparison with the Consensus Immunome Clusters (CICs), previously generated in the same LMC dataset based on 380 immune genes (13), showed that the two most prognostically contrasted LMC classes (class 1 and class 3) had a near-perfect match with CIC 2 (high immune) and CIC 3 (low immune/β-catenin high), respectively (Supplementary Fig. S4), whereas the rest of LMC classes were a mixture of CICs. Cluster 1 had correspondingly a higher proportion of tumors with histologic evidence of brisk TILs (36% compared with 8% in class 3). Analyzing the correlation between the Gerami genes and LMC signature genes showed that the Gerami genes positively correlated with the genes upregulated in LMC class 5 tumors and negatively correlated with genes upregulated in LMC class 3 tumors (Supplementary Fig. S5). Consistent with this, Gerami clusters 1 and 2 highly overlapped with LMC classes 3 and 5, respectively (Supplementary Fig. S6).

LMC class 6 predicted worse prognosis within AJCC stage I tumors. Further biological network analysis identified JUN as a key upregulated nodal gene in this class (Fig. 5A and B). The NGS-based CNA data from a subset of LMC tumors (n = 272) indicated that class 6 tumors were more likely to have DNA amplifications of JUN than other classes (P = 0.003, Fig. 5C, Supplementary Fig. S7). In melanoma, JUN has been reported to activate epithelial-to-mesenchymal transition (EMT), and accordingly a 6-gene–based (32) and 200-gene–based EMT signature (33) consistently scored higher in LMC class 6 than in all other LMC classes (Fig. 5D, Supplementary Fig. S7). A secondary key nodal gene NFKB1 identified to be upregulated in class 6 had no copy number changes. Further examination of immunohistochemically stained sections showed that all four tumors stained from class 6 were positive for NFKB1 protein expression, and this was similar to other LMC classes (P = 0.4, Supplementary Fig. S7).

The value of the LMC signature in predicting outcome in patients treated with immunotherapy was assessed in three disparate clinical trial cohorts of metastatic melanoma (Fig. 5F; refs. 10–12). In the cohort of Hugo and colleagues, tumors classified as class 6 were mainly nonresponders to PD-1 blockade in comparison to the other LMC classes (P = 0.03). Hugo and colleagues reported that expression of AXL predicts poor response to PD-1 blockade; the gene expression data revealed significantly higher AXL expression in class 6 tumors when compared with other classes within their cohort (Fig. 5G). Similarly, for the cohort reported by Ulloa-Montoya and colleagues, class 6 tumors showed a significantly higher proportion of nonresponders to MAGE-A3 immunotherapy in comparison to other classes. The cohort reported by Riaz and colleagues was predominantly composed of nonresponders to anti-CTLA-4 further treated with PD-1 blockade but LMC classes were not convincingly predictive but class 3 predicted poor prognosis, which was consistent with the LMC dataset when compared with good prognosis class 1 (Fig. 5H).

In this study, transcriptome classification was performed utilizing a large population-ascertained cohort of primary melanomas, revealing classes having prognostic value in stage I disease. In stage I tumors, the LMC signature predicted outcome comparably to AJCC staging including SNB. Furthermore, evidence suggests that the signature predicted outcome in patients treated with immunotherapies.

Given the rising incidence of early-stage tumors and the cost of adjuvant therapies to health services and to patients in terms of toxicity, there is an urgent need to identify better prognostic and predictive biomarkers for early-stage disease. When previous gene signatures were applied to the LMC (3, 8), the signatures robustly predicted outcome when the dataset was analyzed as a whole but failed to do so in stage I tumors alone. Although the full Gerami signature was not available, analyzing the prognostic value of genes reported in that study (6) showed that the genes were predictive of prognosis in the whole LMC dataset but not in stage I tumors. In this work, a six-class signature (Supplementary data file) was identified, which was prognostic not only in the whole LMC dataset but also in patients diagnosed at AJCC stage I. The prognostic value of the LMC signature was validated in an independent cohort of primary melanoma built in Lund (4) although the number of stage I cases in this cohort was insufficient to allow replication of the signature's prognostic value in stage I disease.

The LMC signature showed limited overlap with the Lund and TCGA signatures. When comparing it with previously identified immunome clusters by our group (13), two LMC classes strongly overlapped with immune subgroups. The nonoverlapping classes could not be clearly discriminated using the immunome clusters suggesting that these LMC classes are driven by different genomic mechanisms. Comparison of LMC signature genes with Gerami genes indicated a biological pathway overlap as Gerami genes were found to be strongly correlated with LMC classes 3 and 5.

Although SNB is an important melanoma staging tool, the surgery is associated with morbidity (34, 35). In the LMC, SNB was observed to be of prognostic value in the whole dataset and in stage I tumors. However, the LMC signature performed just as well. Given the morbidity of SNB, it may be argued that the LMC signature should be tested in an independent study as a possible alternative to this procedure especially in stage I disease where the likelihood of a positive result is overall low and must be weighed against morbidity.

In melanoma, increased immune gene expression has been consistently shown to predict good prognosis (5, 9, 13, 36). However, a subset of tumors (LMC class 6) was observed which, despite showing immune gene expression, resulted in the patient's early death. Further biological characterization of this class identified copy-number amplifications and increased expression of JUN. Ramsdale and colleagues (37) have shown that JUN promotes an invasive cell phenotype through activation of the EMT pathway, and a higher scoring EMT signature in LMC class 6 confirmed increased activity of the EMT pathway in this class. Riesenberg and colleagues (38) have reported that increased JUN expression leads to proinflammatory and stress signals that promote cytokine expression in coordination with NF-κB. Again, these findings are consistent with the presented transcriptomic observations of JUN and NFKB1 in defining LMC class 6 (Fig. 5B and E). There was insufficient tissue to carry out immunohistochemistry for JUN; therefore, JUN protein expression in the TCGA dataset was examined and confirmed a positive correlation between JUN gene transcription and protein expression (Supplementary Fig. S7). Collectively, these data are indicative of copy-number gains resulting in both increased gene expression and transcriptional activity of JUN in LMC class 6 tumors, although further proteomic studies would be required to confirm this.

The LMC signature was associated with response to immunotherapies; specifically, class 6 associated with poor outcome in two of the three tested datasets. None of these datasets are sufficiently large to make clear inferences. It is of note that the expression of AXL, a known marker for immune evasion, was significantly upregulated in LMC class 6 in metastatic melanoma samples in the Hugo dataset.

The inherent strength of this study is the relatively large size of the population-ascertained cohort. A corresponding limitation is the lack of a well powered AJCC stage I dataset to allow independent replication of the signature in stage I melanoma. Another limitation of this study is that only one-third of LMC patients had an SNB, limiting the power to compare staging tests. The LMC recruitment period preceded the advent of both immunotherapy and targeted therapy, and only a very small number of the study participants have been treated with these drugs. Excluding the samples from these participants showed no modifying effect of such treatments on MSS in the LMC dataset (data not shown).

In conclusion, this study presents a novel signature with demonstrated prognostic value similar in magnitude to that of AJCC staging of melanoma but having added value in stage I melanoma. The data further confirm that AJCC stage largely captures biological variation associated with survival. The LMC class signature prognostic value was similar to that of SNB in the whole dataset (where their effects were additive) and in stage I disease. The signature predicted poor outcome in patients receiving immunotherapies and, in particular, identified high-JUN/high-AXL as a tumor phenotype with poor prognosis in early- and advanced-stage melanoma albeit in very small datasets. This signature has the potential to be trialed as a biomarker in clinical monitoring programs and may help in early identification of patients who may or may not benefit from adjuvant therapies.

No potential conflicts of interest were disclosed.

The accession number for the microarray data reported in this article is European Genome Archive accession number: EGAS00001002922.

Conception and design: R. Thakur, J. N.-Bishop, J.H. Barrett, J. Nsengimana

Development of methodology: R. Thakur, D.T. Bishop, J. N.-Bishop, J.H. Barrett

Acquisition of data (provided animals, acquired and managed patients, provided facilities, etc.): J.P. Laye, S.J. O'Shea, A. Filia, M. Harland, J. Gascoyne, M. Chan, D.T. Bishop, J. N.-Bishop, J. Nsengimana

Analysis and interpretation of data (e.g., statistical analysis, biostatistics, computational analysis): R. Thakur, M. Lauss, J.M.S. Diaz, J. Poźniak, G. Jönsson, D.T. Bishop, J. N.-Bishop, J.H. Barrett

Writing, review, and/or revision of the manuscript: R. Thakur, J.P. Laye, M. Lauss, J.M.S. Diaz, S.J. O'Shea, J. Poźniak, M. Harland, J.A. Randerson-Moor, G. Jönsson, D.T. Bishop, J. N.-Bishop, J.H. Barrett, J. Nsengimana

Administrative, technical, or material support (i.e., reporting or organizing data, constructing databases): J.P. Laye, M. Harland, J. Gascoyne, J.A. Randerson-Moor, M. Chan, T. Mell, J. N.-Bishop

Study supervision: J.P. Laye, M. Lauss, J. N.-Bishop, J.H. Barrett, J. Nsengimana

We thank all the participants of the LMC study and the research nurses who conducted the recruitment.

This work was funded by Cancer Research UK C588/A19167, C8216/A6129, and C588/A10721, and NIH CA83115. R. Thakur, J.M.S. Diaz, and J. Poźniak are supported by Horizon 2020 Research and Innovation Programme no. 641458 (MELGEN). Copy-number data were generated using AICR grant 12-0023.

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.