Integrated multi-omics profiling reveals the role of the DNA methylation landscape in shaping biological heterogeneity and clinical behaviour of metastatic melanoma

Abstract

Background: We developed an integrated multi-omics analysis in metastatic melanoma (MM) cohorts to associate DNA methylation profiles with tumor progression, survival, response to adjuvant immunotherapy, structure of the tumor immune microenvironment and transcriptional programs of immunity and melanoma differentiation.

Methods: Lesions (n = 191) from a fully annotated, retrospective cohort of 165 AJCC 8th Stage III and IV melanoma patients (EPICA cohort) were characterized by reduced representation bisulfite sequencing, RNA sequencing, whole exome sequencing, quantitative immunohistochemistry and multiplex immunofluorescence analysis. The TCGA melanoma datasets were used for validation. Pre-therapy lesions (n = 28) from a cohort of MM patients treated with adjuvant immune checkpoint blockade were characterized for the DNA methylation profile. Impact of a DNMT inhibitor on DNA methylation and transcriptomic profiles of melanoma cell lines was investigated by EPIC arrays and Clariom S arrays.

Results: Four tumor subsets (i.e. DEMethylated, LOW, INTermediate and CIMP) with progressively increasing levels of DNA methylation were identified in EPICA, TCGA MM and TCGA primary melanoma cohorts. EPICA patients with LOW methylation tumors exhibited a significantly longer survival and a lower progression rate to more advanced AJCC stages, compared to patients with CIMP tumors. In an adjuvant immune checkpoint blockade cohort, patients with DEM/LOW pre-therapy lesions showed significantly longer relapse-free survival compared to those with INT/CIMP lesions. RNA-seq data analysis revealed that LOW and CIMP EPICA tumors showed opposite activation of master molecules influencing prognostic target genes, and differential expression of immunotherapy response and melanoma differentiation signatures. Compared to CIMP tumors, LOW lesions showed enrichment for CD8⁺ TCF-1⁺ PD-1⁺ TIM-3^- pre-exhausted and CD8⁺ TCF-1^- PD-1⁺ TIM-3⁺ exhausted T cells, more frequent retention of HLA Class I antigens and a de-differentiated melanoma phenotype. The differentiation and immune-related transcriptional features associated with LOW vs CIMP lesions were tumor-intrinsic programs retained in-vitro by melanoma cell lines. Consistently, treatment of differentiated melanoma cell lines with a DNMT inhibitor induced global DNA de-methylation, promoted de-differentiation and upregulated viral mimicry and IFNG predictive signatures of immunotherapy response.

Conclusions: These results reveal the biological, prognostic and therapeutic relevance of DNA methylation classes in MM and support methylome targeting strategies for precision immunotherapy.

Keywords: DNA methylation; DNMT inhibitor; Immune checkpoint blockade; Immune contexture; Melanoma.

Fig. 1

DNA Methylation profiling of MM EPICA cohort identifies four classes: clinical relevance and transcriptional programs. A Consensus clustering of n = 191 lesions from n = 165 AJCC stage III/IV MM patients (EPICA cohort) based on 4,064 most variable CpG sites. B Alluvial plot showing the initial AJCC stage of n = 165 patients in the EPICA cohort and the final AJCC stage at the end of follow-up. C Stacked bar plot showing the percentage of samples in EPICA cohort according to stage progression occurring (TRUE) or not (FALSE) after the surgical resection of the initial lesion used for methylation profiling. D, E Kaplan–Meier survival curves as a function of methylation cluster in EPICA cohort (D, n = 165) and in TCGA MM cohort (E, n = 368). Statistical analysis by Chi-square (C) and log rank test (D, E). F Bar plot for patients in an adjuvant ICB therapy cohort showing the percentage of pre-adjuvant samples in DEM/LOW vs INT/CIMP methylation groups, stratified by clinical outcome. Statistical significance by the chi-square test. G Kaplan–Meier relapse-free survival curves for patients in the pre-adjuvant cohort (n = 28), stratified by methylation class and grouped into DEM/LOW and INT/CIMP categories. The corresponding risk table is shown on the right. H Dot plot of the normalized enrichment score (NES) of significant Biological Processes (BP) obtained from GSEA (adjusted p-value < 0.05) in DEM, LOW, INT and CIMP methylation-defined classes of n = 187 lesions from n = 165 patients of the EPICA cohort. Dot color represent NES and the size represent p-value resulting from each comparison (one group vs all the others). GO terms selected are the top and bottom 25 significant for each comparison. I Unsupervised clustering of 191 samples of the EPICA cohort based on 72 KEGG pathways significantly enriched in at least the 80% of samples (absolute logit NES > 0.58 and p-value < 0.05)

Fig. 2

Top master molecules (Upstream Regulators) activated in LOW and CIMP MM classes regulate target gene signatures with prognostic significance. A Heatmap of IPA computed z scores for top UR predicted to be activated (red) or inhibited (blue) in each of the four methylation-defined clusters. B, C Stacked bar plots showing for each methylation cluster in the EPICA cohort the percentage of samples with expression of the IFNG (B) or TREX1 (C) target genes above or below the median z score value of each signature. D, E Kaplan–Meier survival curves of patients in the EPICA cohort (top plot) or TCGA MM cohort (bottom plot) according to median z score expression (from RNA-seq profiling) of IFNG (D) or TREX1 (E) target gene signatures. In D and E, patients in both cohorts were grouped according to median UR target gene expression above (“HI”) or below (”LO”) the median z score value of the signature. Statistical analysis B, C by Chi-square; in D, E by log rank test

Fig. 3

Expression of ICB predictive gene signatures in methylation-defined classes of the EPICA cohort. A top panels: dot plot of normalized enrichment scores (NES) from GSEA of seven anti-PD1 response signatures in the four methylation subsets of the EPICA and TCGA metastatic melanoma cohorts. Dot size represents the adjusted p-values and scale colors represent the NES (top). Bottom panels: dot plot of immune-response predictor scores (IRP) MIRACLE and IMPRES. Dot colors represent the mean score in each group and size the adjusted p-values (bottom). B-D Heatmaps of median z score expression of genes in the viral mimicry (B), IFN ICB response (C) and MES resistance (D) signatures in the four methylation subsets. Genes identified by dots in panel C represent a core HLA Class I APM gene set. E–G Stacked bar plots showing for each methylation cluster in the EPICA cohort the percentage of samples with expression of the viral mimicry (E), IFNG ICB response (F) and MES (G) signatures above or below the median z score value of each signature. Statistical analysis in B-D by Kruskal Wallis test followed by Dunn’s multiple comparison test; in E–G by chi square.*: p < 0.05; **: p < 0.01, ***: p < 0.001; ****: p < 0.0001

Fig. 4

Infiltrating T cells, T_PEX and T_EX are enriched in LOW melanomas compared to CIMP melanomas. A Violin plots showing expression, by semi-quantitative IHC, of CD3, CD4, CD8, PD-1, PD-L1, CD68 and CD163 in extra-tumor or intra-tumor compartments of n = 191 EPICA lesions classified according to DEM (n = 39), LOW (n = 50), INT (n = 60), CIMP (n = 42) methylation classes. Data expressed as IHC score (see Supplemental Methods). B,C Multiplex immunofluorescence analysis of a representative LOW (B) and CIMP (C) lesions. In B and C, the H&E image (top) shows the area used for tissue segmentation (bottom) with tumor and stroma identified in dark red and black, respectively. A higher magnification field of the same area shows the density and position of 5 main CD8⁺ T cell phenotypes identified based on differential expression of TCF-1, PD-1 and TIM-3 and color-coded as indicated. Visualization of tumor cells was omitted in B,C. D Density (cells/mm²) in tumor, stroma and whole tissue (tumor + stroma) of the 5 CD8⁺ subsets defined by differential expression of TCF1, PD-1 and TIM-3 in LOW (n = 17) and CIMP (n = 16) lesions. Statistical analysis: in A by Kruskal Wallis test followed by Dunn’s multiple comparison test; in D, by Mann Whitney test for LOW vs CIMP comparisons in each microenvironment compartment and by Friedman multiple comparison test for tumor vs stroma vs tumor + stroma comparisons within each methylation subset. *: p < 0.05; **: p < 0.01; ***: p < 0.001; ****: p < 0.0001

Fig. 5

Downmodulation/loss of expression of HLA class I molecules on tumor cells is more frequent in CIMP compared to LOW lesions. A Ranking of EPICA lesions (n = 175), colored by methylation subset, according to median expression z score (from RNA-seq profiling) of the HLA Class I APM signature and according to expression of HLA Class I molecules on tumor cells by IHC. B Stacked bar plots showing for each methylation cluster in the EPICA cohort the percentage of samples in each of six classes of HLA Class I antigen expression on tumor cells (by IHC and quantitative digital pathology analysis). C Association of HLA class I antigen expression on tumor cells in the EPICA cohort with progression to any subsequent AJCC stages. D Kaplan–Meier survival analysis of patients in the EPICA cohort according to expression of HLA Class I antigens on tumor cells. E Association of methylation classes with immune contexture and with expression of HLA Class I molecules on tumor cells. Expression of HLA Class I and of immune markers, by IHC, was dichotomized into “HI” and “LO” groups corresponding to expression above or below the median value of each parameter. Statistical analysis in A by Spearman correlation; in B,E by Chi-square; in C by Mann–Whitney test, in D by log rank test. *: p < 0.05; **: p < 0.01; ***: p < 0.001; ****: p < 0.0001

Fig. 6

Differential expression of melanoma differentiation signatures and of immunotherapy response signatures in MITF/PMEL^HI vs MITF/PMEL^LO cell lines. A Expression of MITF and PMEL genes, by qPCR, in 46 melanoma cell lines. Subsets of cell lines with high (HI), intermediate (MID) and low (LO) MITF/PMEL expression are highlighted by the indicated color code. Cell line SBL118 was generated from the surgical specimen of the lesion corresponding to sample INT177 of the EPICA cohort. B Heatmaps showing expression of the 7 melanoma differentiation sub-signatures described by Tsoi et al. [59] in four melanoma cell lines. C,D,E Heatmaps showing expression of the TEADS invasive signature (C), of the viral mimicry (D) and of the IFNG-ICB response signatures (E) in three sets of melanoma cell lines (MITF/PMEL^HI, MITF/PMEL^MID, MITF/PMEL.^LO). Statistical analysis in A by Spearman correlation; in B by two-way anova followed by Tukey’s multiple comparison test; in C, D, E by Kruskal Wallis test followed by Dunn’s multiple comparison test. *: p < 0.05; **: p < 0.01; ***: p < 0.001, ****: p < 0.0001

Fig. 7

DNA de-methylation, promoted by a DNMTi, shifts the transcriptional profile of differentiated melanoma cell lines towards the “immune-high” phenotype found in LOW lesions. A Global DNA methylation profiles, from EPIC array data, of six differentiated melanoma cell lines untreated or treated with guadecitabine (GUA) for 7 to 21 days. B Changes in methylation level induced by guadecitabine treatment at methylation sites corresponding to genes in the viral mimicry signature in 6 differentiated melanoma cell lines cultured with guadecitabine for 7 to 21 days. C,D Treatment with guadecitabine for 7 to 21 days of six melanoma cell lines modulates genes in the viral mimicry signature (C) and in the IFNG ICB response signature (D). Statistical analysis in C-E by Mann–Whitney test. *: p < 0.05; **: p < 0.01; ***: p < 0.001; ****: p < 0.0001