Skip to main content
Download PDF

Genome-wide promoter methylation profiling in a cellular model of melanoma progression reveals markers of malignancy and metastasis that predict melanoma survival

Clinical Epigenetics volume 14, Article number: 68 (2022) Cite this article

Abstract

The epigenetic changes associated with melanoma progression to advanced and metastatic stages are still poorly understood. To shed light on the CpG methylation dynamics during melanoma development, we analyzed the methylome profiles of a four-stage cell line model of melanoma progression: non-tumorigenic melanocytes (melan-a), premalignant melanocytes (4C), non-metastatic melanoma cells (4C11−), and metastatic melanoma cells (4C11+). We identified 540 hypo- and 37 hypermethylated gene promoters that together characterized a malignancy signature, and 646 hypo- and 520 hypermethylated promoters that distinguished a metastasis signature. Differentially methylated genes from these signatures were correlated with overall survival using TCGA-SKCM methylation data. Moreover, multivariate Cox analyses with LASSO regularization identified panels of 33 and 31 CpGs, respectively, from the malignancy and metastasis signatures that predicted poor survival. We found a concordant relationship between DNA methylation and transcriptional levels for genes from the malignancy (Pyroxd2 and Ptgfrn) and metastasis (Arnt2, Igfbp4 and Ptprf) signatures, which were both also correlated with melanoma prognosis. Altogether, this study reveals novel CpGs methylation markers associated with malignancy and metastasis that collectively could improve the survival prediction of melanoma patients.

Introduction

Cutaneous melanoma is a type of skin cancer arising from the malignant transformation of melanocytes. According to the World Health Organization, the worldwide incidence of melanoma surpasses 132,000 cases per year. Despite melanoma having a low incidence compared to other types of skin cancer, its aggressive behavior and potential to develop metastasis results in a high incidence of deaths when detected in later stages [1, 2].

DNA methylation is a well-studied epigenetic mark characterized by the incorporation of a methyl group to cytosines by DNA methyltransferases (DNMTs). In cancer cells, abnormal methylation patterns can lead to genome instability, oncogene expression, and silencing of tumor suppressor genes. While cancer cells are characterized by global DNA hypomethylation, the hypermethylation of gene promoters has an essential role in gene silencing and impaired cell proliferation [6]. In melanoma, epigenetic abnormalities occur along with genetic alterations. Changes in methylation patterns have already been described, including hypermethylation of tumor suppressor genes such as PTEN and RAR-b2, and hypomethylation of repetitive elements and oncogenes, such as MAGE and GAGE families of cancer-testis genes [7,8,9]. Besides the existing DNA methylation data in melanoma, the dynamics of those alterations along tumoral progression is still widely unknown.

Studies using melanoma tumor samples have the downside of including stromal and infiltrating immune cells along with melanoma cells, which can mask the discovery of specific cancer cell biomarkers [12, 13]. Even studies conducted in recent years using single-cell technology have not yet achieved a good resolution and scalability to detect genome alterations that could be achieved only by using bulk methods [14]. To get insights into the dynamic alterations of DNA methylation occurring along with melanoma progression, we have studied epigenetic changes using a four-stage cellular model developed by our group [15, 16]. This model comprises different cell lines representing the main steps of melanoma progression: non-tumorigenic melanocytes (melan-a), premalignant melanocytes (4C), non-metastatic (4C11−), and metastatic (4C11+) melanoma cells. The in vitro characterization of the model has found that the 4C cell line resembles premalignant melanocytes (i.e., non-transformed, mesenchymal-like and undifferentiated cells), whereas 4C11− are non-metastatic, mesenchymal-like and undifferentiated melanoma cells, and 4C11+ as metastatic, highly proliferative and differentiated melanoma cells [17–22, summarized information may be found in the Table S8 in [16]]. A progressive global DNA hypomethylation was observed during melan-a malignant transformation, and both protein level and mRNA expression of DNMTs were previously shown to be altered among cell lines [21]. These results show that epigenetic marks are involved in the malignant transformation of melanocytes in this model [16]. Corroborating our findings, Preston-Alp and colleagues revealed DNA methylation as a key molecular mechanism of melanomagenesis induced by UV radiation (UVR) [23]. The authors showed that UVR directly causes stable changes in the DNA methylome and transcriptome, which affect signaling pathways with role in melanocyte biology. More importantly, these alterations correlate to methylation changes observed in melanoma.

Considering the potential of this melanoma progression model to better understand the malignant transformation of melanocytes up to intermediate and metastatic stages, we analyzed here the genome-wide methylation patterns of CpG-rich areas of the four cell lines using the method of Enhanced Reduced Representation Bisulfite Sequencing (ERRBS) [24]. Then, we evaluated CpG sites in promoter regions to discover epigenetically regulated genes associated with malignancy and metastasis. To identify these alterations, we have taken two approaches: first, we compared the CpG methylation profile of melan-a cells with those of 4C, 4C11− and 4C11+ to identify genes altered during the malignant transformation of melanocytes (malignancy signature). Second, we identified CpG markers and genes that could be specifically related to the aggressive features of metastatic cell line 4C11+ (metastasis signature). Using these signatures, we have conducted a multivariate survival analysis using the TCGA-SKCM cohort and found 140 genes presenting at least 2 CpGs as potential methylation biomarkers for melanoma patients’ prognosis. Furthermore, we discovered CpG methylation panels which provide risk scores with prognostic value for melanoma malignancy and metastasization. In summary, we have identified potential genes regulated by CpG methylation and, more importantly, CpG methylation signatures with potential prognostic value by providing risk scores associated with melanoma malignancy and metastasis.

Results

DNA methylation profiling of cell lines corresponding to different stages of melanoma progression

The DNA methylation profiling of the cell lines in the model—melan-a, non-tumorigenic melanocytes [25]; 4C, premalignant melanocytes; 4C11−, non-metastatic melanoma cells; and 4C11+, metastatic melanoma cells [15, 16, 20]—was performed using ERRBS [24], a technique developed to obtain the methylation status of CpGs at base-pair level, with a focus in CpG islands and shores. An average of 1,406,315 cytosines (with a standard deviation of 31,656) in CpG context at minimum 10X coverage was uncovered for each sample. The mapping to the bisulfite converted mouse genome (mm10 assembly) using Bismark [26] had an efficiency of 71% on average. As expected, a bimodal distribution of the methylation status was observed for all cell lines, with peaks in 0 and 100% methylation (Additional file 7: Figure S1).

A methylation average of 47.4% was obtained among CpGs for melan-a melanocytes, whereas 35.8% for 4C, 34.2% for 4C11−, and 36.8% for 4C11+ cell lines. A similar pattern was obtained when analyzing the global content of 5-methylcytosines by HPLC, as shown in Fig. 1A. Methylated cytosines are averaged in 4% of the genome for melan-a lineage and 3.8% for the other cell lines (ANOVA test F = 1.273, p-value = 0.38). The principal component analysis of these data (Fig. 1B) effectively clustered the biological replicates from the same cell lines. Also, it revealed that the DNA methylation status discriminates each cell line from our model, with pre-malignant 4C and malignant non-metastatic 4C11− cell lines displaying similar overall methylation values. This could also be observed by the patterns shown in the heatmap (Fig. 1C), which also include an unsupervised hierarchical clustering that clear stratifies the samples based on their respective cell lines, and also clusters together 4C and 4C11− cells, the two cell lines with similar undifferentiated/mesenchymal phenotype [21].

Fig. 1
figure 1

Distinct methylation patterns characterize cell lines representing melanoma progression. A. Global DNA methylation content analyzed by HPLC shows the percentage of 5-methylcytosine in the genome for two biological replicates from the cell lines melan-a, 4C, 4C11−, 4C11+. B. Principal component analysis of DNA methylation data from three biological samples of the melan-a, 4C, 4C11− and 4C11+ cell lines shows close values of Principal Components 1 and 2 for samples from the same cell line. C. A heatmap represented with light yellow for regions with low methylation and dark blue for regions with high methylation shows the top most variant regions of 1,000 bp for each sample. Unsupervised hierarchical clustering in the upper part of the heatmap groups 4C and 4C11− in a common branch closer to 4C11 + than melan-a melanocytes. melan-a: non-tumorigenic melanocyte lineage; 4C: premalignant melanocyte lineage; 4C11−: non-metastatic melanoma cell line; 4C11 + : metastatic melanoma cell line

Full size image

The distribution of CpG methylation across the genome in different stages of melanoma progression

The differential methylation analysis was conducted using the package methylKit [27] in the R environment [28], with multiple hypotheses testing correction by the SLIM method [29]. CpGs with a q-value less than 0.01 and percentage of differential methylation higher than 25% were selected. An overview of the comparative results between the four cell lines from the model can be seen in Fig. 2 and Additional file 8: Figure S2. Additional file 8: Figure S2 contains the heatmaps for pairwise comparisons between the cell lines. A scheme of melanoma model progression can be seen in Fig. 2A, depicting the stepwise transformation of the cells starting from the non-tumorigenic melanocytes melan-a. The comparison between 4C, the premalignant melanocyte lineage and melan-a, the non-tumorigenic parental melanocyte lineage, showed more hypomethylated (51,563) than hypermethylated CpGs (9,190) (Fig. 2B), both equally distributed between CpG islands and other regions, and less among CpG shores (Fig. 2C, D, E, respectively). In parallel, the comparison between the intermediate stages of the model, i.e., between non-metastatic melanoma 4C11− and premalignant 4C cell lines, which display epithelial-to-mesenchymal transition features, had the lowest number of differentially methylated cytosines (DMCs) (18,549) (Fig. 2B). DMCs in other regions of the genome not including CpG islands and shores were more frequent seen in this comparison (Fig. 2C, D, E). Finally, the pairwise comparison between the metastatic 4C11+ and the non-metastatic 4C11− melanoma cell lines revealed the highest number of DMCs (131,202) (Fig. 2B). More than half of the hypermethylated cytosines (41,912) were located in CpG islands (Fig. 2C, D, E).

Fig. 2
figure 2

The pattern of differentially methylated CpGs distributed across genomic features and transcription features. A Cell lines used in this study that represent different stages of melanoma progression. B Number of hypermethylated (red) and hypomethylated (blue) individual cytosines in CpG context per pairwise comparison among all the four cell lines. DMCs according to their location in islands, shores, and other regions are represented in C, D, and E, respectively. Annotation of all hypomethylated (F) and hypermethylated (G) CpGs by region related to transcription: promoter (salmon), exon (blue), intron (green) and intergenic (purple) for each pairwise comparison

Full size image

The comparison between melan-a and 4C11+ cells displayed 22.4% of differentially hypermethylated CpGs (21,092) (Fig. 2B), most of them located in CpG islands (14,592) (Fig. 2C, D, E). The comparison between 4 and 4C11+ cells showed similar results to those obtained in the comparison between 4C11− and 4C11+ (140,261 DMCs). An even distribution between CpGs that gained or lost methylation was observed in this comparison, with the hypermethylated CpGs concentrated in CpG islands (Fig. 2C, D, EE). Last, the pairwise comparison between the 4C11− cell line and melan-a cells revealed a similar pattern to the one found between melan-a and 4C cells. Most differentially methylated CpGs found were hypomethylated in 4C11− cells (Fig. 2B) and only ~ 11% were located in shores (Fig. 2C, D, E).

We also categorized DMCs according to their transcript annotation (promoter, exon, intron or intergenic). The sum of the percentages shown in Figs. 2F (hypomethylated CpGs) and 2G (hypermethylated CpGs) does not total 100% because there is an overlap between transcriptional regions annotated in the genome (for example, between first exons and promoters).

For hypomethylated CpGs, we observed a similar distribution of differential CpGs in the promoter, exon, intron, and intergenic regions across all comparisons, except for the comparison between 4C11− and 4C (Fig. 2F). All comparisons showed a higher frequency of differential CpGs in the promoter regions (41–43%), except for the comparison between 4C11− and 4C, in which they were more located at intergenic regions. For the hypermethylated CpGs, there was a more heterogeneous distribution of CpGs across the transcript annotations, with the most considerable portion of CpGs annotated to promoter regions in all comparisons against the metastasis-prone 4C11 + cell line. In particular, the pairwise comparison between melan-a and 4C11 + cells exhibited the most extensive annotation of differentially methylated CpGs is in promoter (59%) and exon regions (41%). These exon regions probably comprise the first exons of the related transcripts, which also include promoter regions.

Genes characterized by differential CpG promoter methylation during melanoma progression

We next decided to identify potential genes associated with melanoma progression that could be epigenetically regulated by CpG methylation in promoter regions. The CpGs located in the region from -500 to 500 nucleotides surrounding the TSS (transcription start site) and genes with more than three differentially methylated CpGs were considered in the selection of the relevant genes. The genes with differentially methylated CpGs in their promoter regions are referred hereafter as differentially methylated genes. The complete list of identified genes and CpG sites can be found in Additional file 1: Table S1.

The number of hypermethylated and hypomethylated genes in each pairwise comparison is indicated in Fig. 3A; these numbers have followed the same pattern observed in Fig. 2A across the same pairwise comparisons. In addition, the distribution of the number of CpGs in each annotated gene at each pairwise comparison can be seen in Fig. 3B. All pairwise comparisons containing the metastatic 4C11 + cell had few genes with more than 100 CpGs hyper- or hypomethylated.

Fig. 3
figure 3

Genes and CpG sites annotated as differentially methylated in each pairwise comparison. A Number of hypermethylated and hypomethylated genes per pairwise comparison. Differentially methylated genes were considered as those with differentially methylated CpG sites between 500 bp up and downstream of the TSS, and a minimum of three differentially methylated CpGs. B Frequencies of CpG count per gene for each pairwise comparison, with peaks around five cytosines having the greatest number of differentially methylated genes

Full size image

Gene signatures associated with DNA methylation changes occurring in early and late stages of melanoma progression

We sought to identify gene signatures associated with promoter methylation alterations occurring early in malignant transformation (4C pre-malignant cells) and lasting until late stages of melanoma progression (4C11− and 4C11 + melanoma cells). To accomplish that, we conducted an intersection analysis of the differentially methylated genes identified in the pairwise comparisons against melan-a to identify a malignant transformation signature. Five hundred and forty hypomethylated genes and 37 hypermethylated common genes were identified in this analysis (Fig. 4A, B). Functional enrichment analyses were conducted using the software clusterProfiler [30] and the database Gene Ontology (Biological Process), and significantly enriched functions were selected based on FDR corrected p-values. Among the hypomethylated genes, the most significantly enriched functions were related to negative regulation of cell development, nervous system development, and neurogenesis (Fig. 4C), which indicate the role of epigenetic regulation in cell dedifferentiation that occurs along melanoma progression. Moreover, the most significantly enriched functions among the hypermethylated genes were involved in ossification, cartilage development and cyclic nucleotide metabolic process (Fig. 4D).

Fig. 4
figure 4

DNA methylation signature associated with melanocyte malignant transformation (malignancy signature). Venn diagram of hypomethylated (A) and hypermethylated (B) genes in common between pairwise comparisons of 4C, 4C11− and 4C11 + versus melan-a cell line. Enrichment plots showing the pathways enriched in Gene Ontology: Biological Process for hypomethylated (C) and hypermethylated (D) genes

Full size image

In parallel, we decided to investigate hyper- and hypomethylated genes in the metastatic 4C11 + cell line compared to the other three non-metastatic cell lines (melan-a, 4C and 4C11−), in the attempt to find genes potentially related to metastasis acquisition (metastasis signature). Six hundred and forty-six hypomethylated (Fig. 5A) and 520 hypermethylated (Fig. 5B) genes were found in common between the comparisons 4C11− vs. 4C11 + , melan-a vs. 4C11 + and 4C vs. 4C11 + . Independent functional enrichment analyses were also conducted for the hypomethylated (Fig. 5C) and hypermethylated (Fig. 5D) genes. Among the hypermethylated genes, significantly enriched GO terms were associated with nervous system development: axon development, axonogenesis, axon guidance, neuron projection guidance, and neuron fate commitment. Conversely, the hypomethylated genes in the metastasis signature were related to cell–matrix adhesion, as well as to the negative regulation of nervous system development.

Fig. 5
figure 5

DNA methylation signature associated with metastasis (metastasis signature). Venn diagrams of hypomethylated (A) and hypermethylated (B) genes only in the metastatic 4C11 + melanoma cell line compared to non-metastatic melan-a, 4C and 4C11− cell lines. Enrichment plots showing the pathways enriched in Gene Ontology: Biological Process for hypomethylated (C) and hypermethylated (D) genes

Full size image

DNA methylation alterations as independent prognostic markers for melanoma survival

To further investigate the relationships between the DMCs and their potentially regulated genes identified in the murine model with melanoma patient data, we performed overall survival analyses using cutaneous melanoma data from TCGA (TCGA-SKCM). Additional file 6: Table S6 contains a list of all ERRBS CpGs annotated to mouse genes and the ortholog human gene with the respective Infinium 450 k CpG probes. Using the methSurv tool, we obtained Hazard Ratios (HR) and respective p-values for all the TCGA-SKCM's CpGs annotated to the ortholog genes to the hyper- and hypomethylated genes in the malignancy and metastasis signatures. We obtained statistically significant CpGs (adjusted p-values ≤ 0.01) associated with overall survival, resulting in 19 hyper- and 224 hypomethylated CpGs in the malignancy signature, and 330 hyper- and 253 hypomethylated CpGs in the metastasis signature (Fig. 6A, and Additional file 2: Table S2).We next analyzed the statistically significant CpGs according to their location in islands, shores, shelves and open seas. Interestingly, the number of hypermethylated CpGs correlating with patients’ survival in the metastasis signature were concentrated in islands (167 CpGs), followed by open seas (70 CpGs) and N- (53 CpGs) and S-shores. This same pattern was seen among the hypomethylated CpGs with prognostic value from the metastasis signature; which were also found enriched in islands (94 CpGs), followed by open seas (71 CpGs), and N- (45 CpGs) and S-shores (31 CpGs). In the malignancy signature, open seas and islands sites hosted the most hypermethylated CpGs (7 CpGs each). In contrast, more hypomethylated CpGs correlating with survival were found in islands (85 CpGs), open seas (65 CpGs) and N- (29 CpGs) and S-shores (33 CpGs) (Fig. 6A). Nevertheless, the differences between regions in each signature were not significant (Fisher test, malignancy: p-value = 0.944, and metastasis p-value = 0.3571).

Fig. 6
figure 6

Distribution of hyper- and hypomethylated CpGs correlated with melanoma patients’ outcomes according to their genome location and number of associated genes. A The number of hyper- and hypomethylated CpGs correlated with overall survival is shown according to their location in islands, shores, shelves and open sea regions, for malignancy and metastasis signatures. B The distribution of hyper- and hypomethylated CpGs along genomic regions is shown according to their correlation with a good or poor prognosis for malignancy and metastasis signatures. C. The number of genes from both signatures containing at least two CpGs correlating with melanoma patient outcome is shown

Full size image

We also verified if the CpGs with prognostic value would be enriched in specific genomic locations. Interestingly, in both malignancy and metastasis signatures, the hyper- and hypomethylated CpGs significantly correlated with good prognosis were located mainly in islands, whereas those correlating with poor survival were predominantly in open seas (Fisher test, hypermethylated in malignancy: p-value = 0.1837, hypomethylated in malignancy: p-value = 6.243–12, hypermethylated in metastasis: p-value = 0.0004998, and hypomethylated in metastasis: p-value = 2.333–5) (Fig. 6B).

We also explored the genes whose CpG methylation patterns were related to survival in each signature; a subset of genes containing at least 2 CpGs with significant Cox proportional-hazards models was selected in each signature (Fig. 6C). There were 140 (77 hyper- and 63 hypomethylated) genes identified in the metastasis signature and 55 (5 hyper- and 50 hypomethylated) genes identified in the malignancy signature (Additional file 3: Table S3). Then, only those differentially methylated genes present in similar signatures (malignancy or metastasis) based on gene expression levels (Additional file 5: Tables S5 and Additional file 6: Table S6 from [16]) and with an inverse relationship between expression and DNA methylation patterns were chosen for further evaluation. Finally, the analysis of the prognostic value of these genes based on expression levels was also used to guide gene selection using previous information from the Leeds Melanoma Cohort (Table S8 from [16]; European Genome-Phenome Archive accession number EGAS00001002922), which consists of molecular date from 703 treatment-naïve primary melanoma patients.

Based on the above criteria, TBC1D9 was selected from the hypermethylated subset, and four genes (P2RX7, PRDM1, PTGFRN and PYROXD2) were selected from the hypomethylated group in the malignancy signature. In the metastatic signature, six genes were selected from those hypermethylated (FAM107B, LMX1B, MKX, MYH10, PTPRF and ZBTB16), and eight (ADSSL1, ARNT2, CHN2, DIXDC1, IGFBP4, LMX1B, PCDHGB2 and TBX15) from those hypomethylated..

No genes from the malignancy signature, both from hypo- or hypermethylation subsets, presented expression levels associated with prognostic value. However, for the metastatic signature, the hypermethylated genes FAM107B, PTPRF and ZBTB16 had their low expression correlated with poor overall survival, and the hypomethylated genes ARNT2, IGFBP4 and CHN2 had their expression levels related to prognostic value. Although no gene from the malignancy signature had its expression correlated with patient survival, we decided to select PYROXD2 and PTGFRN for additional analysis (Fig. 7).

Fig. 7
figure 7

Genes from the malignancy signature presenting CpGs’ methylation status with significant prognostic value. Pyrodx2 and Ptgfrn are shown for gene expression (A) and CpG methylation (B) in melan-a, 4C, 4C11− and 4C11+ cell lines. C. Kaplan–Meier plots show the result of overall survival analysis conducted in patients from the TCGA-SKCM cohort, with samples split by CpG methylation status for CpGs of the respective genes. Gene expression is shown as raw fragment counts, with the bar representing standard error. *p ≤ 0.05, **p ≤ 0.01, ***p ≤ 0.001, ****p ≤ 0.0001 for ANOVA statistical test. CpG methylation sites for the promoter region of the genes are shown as methylated (black circle) for ≥ 70% methylation and white circles for < 70% methylation

Full size image

Among the genes in the malignancy signature, we have confirmed that Pyrodx2 and Ptgfrn have increased gene expression in all cell lines compared to melan-a (Fig. 7A), which is in line with the hypomethylation of individual CpGs in their promoter regions in all cell lines compared to melan-a (Fig. 7B). In order to evaluate the effect of DNA methylation on gene expression, the cell lines were treated with the demethylating agent, 5-aza-2’-deoxycytydine (5azaCdR). These genes, highly expressed and hypomethylated in 4C, 4C11− and 4C11− cells compared to melan-a melanocytes, were not derepressed in melan-a cells after 5azaCdR treatment (Additional file 3: Figure S3).

The survival analyses of two CpGs located in the first exons of these genes can be seen in Fig. 7C. Methylation of the PYRODX2 CpGs identified by the probe IDs cg22752023 and cg08397758 independently classify a group of patients with better survival. The same good prognosis profile is seen for the hypomethylation of PTGFRN CpGs cg12406825 and cg03752628, located, respectively, at a distance of 1500 bp from the TSS, in a CpG island, and in the gene body, in a south shore.

In the metastasis signature, Arnt2 and Igfbp4 genes, highly expressed (Fig. 8A) and hypomethylated (Fig. 8B) only in 4C11 + metastatic melanoma cells compared to melan-a, 4C and 4C11− cells, had their expression significantly increased in melan-a, 4C and 4C11− cells after 5azaCdR and 5azaCdR + TSA treatment (Fig. 8C). For both genes, 5azaCdR treatment significantly reverted the gene silencing in melan-a cells, while 5azaCdR + TSA in 4C11 + melanoma cells. In parallel, for Ptprf gene, not expressed in metastatic 4C11 + cells compared to melan-a, 4C and 4C11− cells and presenting concordant increased methylation in almost all CpGs located in the promoter region (Fig. 8A, B, respectively), the treatment with 5azaCdR was able to revert gene silencing in 4C11 + cells (Fig. 8C). Survival analyses showed that the hypermethylation of ARNT2 CpGs in the gene body region correlates with poor survival (Fig. 8D). This is in line with the finding that methylation in the body of this gene is correlated to its higher expression [5; 32–36]. For the IGFBP4 gene, the hypermethylation of two CpGs, once located in the 1500 bp distance from the TSS, in a north shore, and another located 200 bp of the TSS, inside a CpG island, predicted good patient survival outcomes (Fig. 8D). Interestingly, the methylation status of two CpGs from PTPRF exhibited different outcomes for patient survival: cg05661060, located in the gene body, in a shelf region, is hypomethylated in the patient subgroup with poor survival. In contrast, cg06796515 (located on an island) is hypermethylated in patients with poor prognosis.

Fig. 8
figure 8

Genes from the metastasis signature with CpGs methylation status correlated with melanoma patient survival. The genes Arnt2, Igfbp4, and Ptprf are shown for gene expression (A, E and I, respectively) and CpG methylation (B, F and J, respectively) in melan-a, 4C, 4C11− and 4C11+ cell lines. C, F and K. Cell lines presenting, respectively, low Arnt2, Igfpb4 and Ptprf expression were treated for 48 h only with 5azaCdR or for additional 18 h with Trichostatin A (TSA) and had their expression evaluated by RT-qPCR. D, H and L. Kaplan–Meier plots show the result of overall survival analysis conducted in patients from the TCGA-SKCM cohort, with samples split by CpG methylation status for CpGs of the respective genes. Gene expression is shown as raw fragment counts, with the bar representing standard error

Full size image

Multi-CpG methylation panels of malignancy and metastasis markers discriminate patient prognosis

After analyzing the importance of single CpGs methylation status to patient survival, we decided to investigate the combined prognostic value of a panel of CpG methylation sites and build more accurate Cox proportional-hazard multivariate models for melanoma prognosis based on the malignancy and metastasis CpG signatures. The characteristics of the patients selected for the models are described in Additional file 4: Table S4.

To determine the subset of relevant predictor variables and avoid overfitting in the Cox model, a LASSO regularization step was applied to employ only variables with non-zero coefficients in the final Cox survival model. This final Cox model was built using the survival R package with the selected CpGs as input. The HR was obtained, along with respective coefficients and p-values, corrected using the Benjamini–Hochberg method for multiple testing hypotheses. Furthermore, the significant CpGs (p-adjusted ≤ 0.05) were selected to compose the panels for malignancy and metastasis signatures, totalizing 33 CpGs in the malignancy (Table 1) and 31 CpGs in the metastasis signature (Table 2). The complete set of the CpGs and covariates used in both models can be found in the Additional file 5: Table S5.

Table 1 Significant covariables in malignancy signature Cox model
Full size table
Table 2 Significant covariables in metastasis signature Cox model
Full size table

Interestingly, there were two CpGs in common among the signatures: cg00478851 and cg01607369. The first CpG exhibited the same direction of methylation regulation, while the second showed opposite directions between the two signatures. cg00478851 is present within a distance of 1500 base pairs from the TSS of PPFIBP2, in a north shore of the closest CpG island; it has HR values of 3.6 and 2.5, respectively, in the malignancy and metastasis signatures, showing its importance for both tumor initiation and metastasization. On the other hand, cg01607369, which is located in the gene body of PPFIBP2, and an open sea region, had an HR of 0.3 in the malignancy but a HR equal of 2.9 in the metastasis signature; it may be concluded that its methylation status is associated with better prognosis during tumor initiation but with worse prognosis at the metastatic state.

The regression coefficients from the CpGs in the Cox multi-CpGs models were used as weights to their respective methylation β-values to calculate a risk score per patient using either the malignancy or metastasis signatures. The Kaplan–Meier survival curves from patients stratified on by high- and low-risk score groups based on the median value are shown in Fig. 9 (log-rank malignancy p-value = 1e-11, metastasis p-value = 2e-13), demonstrating the accuracy of both CpG panels in correctly predicting patient prognosis.

Fig. 9
figure 9

Risk scores for melanoma survival calculated from CpG methylation malignancy and metastasis-related panels accurately separate samples into low and high-risk groups. Kaplan–Meier plots show survival of high and low risk patients according to a (A) Malignancy Risk Score (HR = 3.29 [2.48–4.37]), and (B) Metastasis Risk Score (HR = 1.84 [1.56–2.16]), both obtained from the respective signatures of CpG methylation status, containing, respectively, 33 and 31 CpGs

Full size image

Discussion

In this work, we analyzed the dynamics of DNA methylation patterns in a murine cellular model of melanoma progression. In addition to genetic and chromosomal alterations occurring throughout melanoma progression, it is well accepted that epigenetic alterations play a role in tumorigenesis by either promoter hypermethylation of tumor suppressor genes or global hypomethylation. These epigenetic changes collectively induce an unstable genomic state, prone to genomic mutations and expression of transposable elements [36,37,38].

Despite different studies having described DNA methylation alterations in melanomas, few of them have evaluated these alterations in samples or cells with similar genetic backgrounds corresponding to different stages of melanoma progression. To fill this gap, we have performed a differential methylation analysis using pairwise comparisons between cell lines representing distinct stages of a linear model of melanoma progression. This analysis has provided vital information about the CpG methylation dynamics occurring in each step of the melanoma progression.

All cell lines (premalignant 4C, non-metastatic 4C11−, and metastatic 4C11 +) derived from the melan-a melanocyte lineage showed progressively more hypomethylated CpGs concentrated in gene promoters and equally distributed in either CpG islands or regions not including islands or shores (called “open sea”). This observation is in line with the findings from Ecsedi and colleagues, in which decreased promoter methylation was correlated with increased Breslow thickness in melanomas [39]. Opposingly, the hypermethylation of upstream regulatory regions of genes and the hypomethylation of repetitive genome elements have been described in melanoma cell lines corresponding to stages III and IV compared to normal melanocytes [40]. Another study classified melanomas in three groups, presenting high, intermediate, and low methylation in promoter regions. It showed that most of them were categorized in the high methylation group [41]. These contradictory results could be related to cell lines and/or tumor samples corresponding to different melanoma stages. In addition, melanoma tissue specimens, a diversity of patient treatment procedures could contribute to the heterogeneous molecular alterations observed in different studies. Moreover, adjacent tissue and lymphocyte infiltrates can influence the results of studies using tissue biopsies [11, 12].

The comparison between the metastatic melanoma cell line 4C11 + with premalignant melanocytes 4C and non-metastatic melanoma cells 4C11− showed hypermethylated CpGs concentrated mainly in islands and promoter regions. This is consistent with the findings of Wouters et al., in which higher methylation levels were also seen in CpG islands and promoters of metastatic melanomas compared to nevi and primary melanomas [42]. In our model, 1,351 hypermethylated and 848 hypomethylated genes in the comparison between metastatic 4C11 + and non-metastatic 4C11− melanoma cells, which could be related to the acquisition of metastatic capacity.

Early in melanoma progression (4C vs. melan-a), there are many hypomethylated genes, most with a low CpG count (Fig. 3B). In the next transition between premalignant (4C) and non-metastatic melanoma cell line (4C11−), the number of genes presenting alteration in promoter methylation is small. Differences of methylation per gene do not reach a high number of CpGs (Fig. 3). In the metastatic cell line (4C11 +), which represents the latest stage of progression in our model, both hypo- and hypermethylated genes containing high number of DMCs are present compared to its non-metastatic counterpart 4C11−, and also with premalignant 4C and melan-a melanocytes.

In the metastatic stage of our model, there were more genes with DMCs in their promoters. In a previous study, our group showed a higher protein and mRNA expression of the de novo DNA methyltransferase DNMT3a in 4C11 + compared to all the other cell lines. In contrast, the maintenance DNA methyltransferase DNMT1 is highly expressed in all cell lines derived from melan-a melanocytes [21]. Zhang and collaborators have shown that the ectopic expression of DNMT3a in embryonic stem cells leads to the hypermethylation of 5-methylcytosines in the DNA [43]. Therefore, it is possible that the hypermethylation state and increased number of methylated CpGs per gene observed in 4C11 + cells might be induced by the increased DNMT3a expression. Nevertheless, more studies would be necessary to confirm this hypothesis.

Hypermethylated genes in the metastasis signature were enriched with functions related to nervous system development. Both peripheral neurons and melanocytes are cells that originated from the neural crest; therefore, it is reasonable to conceive developmental genes would be regulated by epigenetic reprogramming during melanoma progression. Numerous developmental genes are upregulated in 4C and 4C11− cells in our model, such as Twist1, Snai1, Gata4, and Nanog [16]. Moreover, these cells exhibit lower expression of E-cadherin and gain of N-cadherin expression, leading to the conclusion that they undergo cell reprogramming and exhibit a mesenchymal and less differentiated phenotype.

The contribution of each discovered DMCs from the malignancy and metastasis signatures was assessed in the context of melanoma survival, using TCGA-SKCM CpG methylation data [44]. There were more CpGs whose methylation levels were related to good prognosis, both in the malignancy and metastasis signatures. They were predominantly hypomethylated in malignancy, similarly to the findings of de Unamuno Bustos et al., where the unmethylation of promoter CpG islands of multiple genes was associated with good prognosis in a cohort containing predominantly primary melanomas (140 of 170 samples) [45]. In the metastasis signature, the CpGs with good prognostic value were mainly hypermethylated, in line with the findings of Tanemura and colleagues, in which MINT31 gene CpG island methylation conferred good prognosis in stage III melanoma [46].

We have found a concordant relationship between lower promoter methylation and higher gene expression levels for Pyroxd2 and Ptgfrn from the malignancy signature, which was also associated with the poorest survival in melanoma patients (Fig. 7). This is the first time that the promoter methylation levels of Pyrodx2 and Ptgfrn have been described in melanoma, which could be further explored as novel individual biomarkers for this disease. PYROXD2 is differentially expressed in basal cell carcinoma [47] and is localized in the mitochondrial inner membrane, where it is associated with increased cell proliferation and ATP production [48], highlighting its role in melanoma survival and metabolism.

Notably, PTGFRN is overexpressed in glioblastoma multiforme and is associated with poor survival [49]. However, its methylation status has not yet been investigated. Glioblastoma and melanoma are closely related types of cancer, both aggressive and with the same developmental origin from neural crest cells.

We have also investigated the relationship between DNA methylation and gene expression levels for three genes from the metastasis signature (Arnt2, Igfbp4 and Ptprf). Arnt2 expression was increased in the metastasis-prone 4C11 + cells and diminished DNA methylation compared to the other cell lines (melan-a, 4C and 4C11−) (Fig. 8). This gene belongs to the bHLH/PAS transcription factor family, and its promoter has already been demethylated in P19 embryonic carcinoma cells [50]. Besides, its higher expression was associated with breast cancer metastasis [51] and a tumorigenic molecular signature in glioblastoma [52]. Curiously, we have found that methylation alterations in two CpGs at this gene body correlated with poor survival. Indeed, there is consistent evidence about the positive correlation between methylation at a gene body and increased gene expression [31,32,33,34,35,36].

IGFBP4, a member of the insulin-like growth factor binding protein that regulates angiogenesis [53], exhibited higher gene expression in the 4C11 + cell line and concordantly lower promoter methylation (Fig. 8). Furthermore, in TCGA-SKCM, the analysis of the methylation levels of CpGs close to the TSS was able to identify a subgroup of patients with low methylation and worst survival outcome. There is vast published literature about Igfbp4 and different types of cancer. Particularly in melanoma, one study compared primary and metastatic tumors and found increased protein levels of IGFBP4 in primary versus metastatic tumors [54]. Elevated transcriptional levels of this gene have also been found in glioblastoma biopsies [55] and lung cancer tissue [56].

Finally, the gene Ptprf was found to have a lower expression level in 4C11 + cells and concomitant hypermethylation of its promoter compared to the other cell lines (Fig. 8). Curiously, one study reported the opposite finding, showing that PTPRF was over-expressed in metastatic and primary melanoma cells in vitro and in situ compared to melanocytes [57]. In contrast, its expression was found to be downregulated in non-small cell lung cancer cell lines [58]. The survival analysis of DNA methylation within PTPRF revealed a remarkable high HR (1.84) for a CpG located in the gene body in a region classified as a CpG island: patients with lower methylation levels exhibited better survival. This particular CpG methylation shall be better evaluated in the clinic as a potential biomarker of melanoma aggressiveness.

To obtain a more accurate prediction of melanoma survival based on CpG methylation levels, we decided to build panels of CpGs using multivariate Cox regression after LASSO regularization. The final models were composed of a malignancy signature panel with 33 CpGs, and another panel from the metastasis signature with 31 CpGs (Tables 1 and 2). These panels were then applied to calculate survival risk scores and more accurately classify patients into low and high-risk groups (Fig. 9). Wouters and colleagues also identified two signatures associated with different 4-year survival outcomes for melanoma using more than 734 DNA methylation markers [42]. However, survival multivariate analyses were not evaluated for these probes. This multivariate approach has already been conducted for clear cell renal carcinoma [59], leukemia [60], esophageal squamous cell carcinoma [61], colorectal cancer [62], lung adenocarcinoma [63], among others. Future validation of the identified CpGs methylation panels in a new cohort will be deemed necessary for their use in the clinic.

Conclusions

Our findings demonstrate the usefulness of the melanoma progression cell model to unravel epigenetic changes occurring during melanoma progression. This model can further contribute to discovering molecules present in one specific disease stage, which could be clinically employed to tailor treatment for patients. Furthermore, the CpGs methylation differences identified in this study were practical to identify relevant genes that participate in the control of phenotypic switching along with melanoma progression.

The limitations of our study are related to the fact that ERRBS is a technique that evaluates CpG sites located mostly in CpG islands; therefore, even with our findings of interesting methylation differences in different regions of the genome, the exploration of differential methylation in non-island or shore regions was limited. Furthermore, we focused on the differential methylation analysis for the CpGs located in the promoter region of the genes. Consequently, further studies comparing CpG methylation differences in non-promoter regions could unveil novel results.

The DNA methylation profiling of the melanoma progression steps in our model revealed markers associated with malignancy and metastasis, which were also predictive of human melanoma survival. These markers include DNA methylation changes in Pyrodx2 and Ptgfrn, which are hypermethylated in melan-a non-tumorigenic cells, Arnt2 and Igfbp4, that are hypomethylated in malignant metastatic 4C11 + cells, and Ptprf, which is hypermethylated in 4C11 + cells. Furthermore, we established panels of CpGs using multivariate prognostic models that could accurately predict survival outcomes in human patients, and deserve to be evaluated in future studies in search of reliable biomarkers for the clinical assessment of melanoma prognosis.

Material and methods

Cell culture

The cell lines 4C, 4C11− and 4C11 + were cultured in RPMI 1640 medium (Gibco) pH 6.9 with addition of 1.0 × 105 U/L antibiotic Pen-Strep (Gibco) and supplement of 5% fetal bovine serum (FBS, Gibco). The cultures were maintained in humidified incubators containing 5% of CO2 at 37 °C. Melan-a cells were cultured under the same conditions, but with addition of 200 nM phorbol 12-myristate 13-acetate (PMA, Amresco), an activator of PKC needed for the growth of melanocytes in culture [24]. Three biological replicates of each cell line, obtained from distinct passages, were obtained in the conditions described above.

Enhanced reduced representation bisulfite sequencing (ERRBS)

To obtain pure samples, DNA was extracted from cultured cells using the Gentra Puregene Cell Kit (Qiagen), according to the fabricant recommendations, including the addition of RNase A solution (contained in the kit). Approximately 200 µL of a minimum concentration of 80 ng/µL were forwarded to the Epigenomics Core Lab of Weill Cornell Medical College, with sufficient quality of double-stranded DNA determined by fluorimetry using QuantiFluor® dsDNA System (Promega), following the Epigenomics Core recommendations. After passing the quality controls, all the samples had their libraries constructed as described by Garrett-Bakelman and colleagues [24]. The samples were multiplexed and distributed in three lanes for sequencing in Illumina HiSeq 2500.

The output files in format.bcl (base call) were demultiplexed and converted to fastq files using the software bcl2fastq. Adaptors and low-quality reads were removed, and the reads were aligned to the mm10 mouse bisulfite converted genome using Bismark [26], to generate CpG calls. These steps were conducted by the Epigenomics Core team.

Differential methylation analyses

Differential methylation analyses were conducted for each of the pairwise comparisons in the R environment [28], using the package methylKit [27]. From the alignment output file, containing minimum of 10 count coverage, counts above 99.9 percentile were filtered out, to avoid PCR bias. The data were corrected for overdispersion for a more stringent statistical analysis. P-values were corrected to q-values using the package default method: SLIM [29].

The parameters for selection of differentially methylated CpGs (DMCs) were: q-value ≤ 0.01 [29] and minimum of 25% of difference in methylation between compared lineages. Annotation of functional features (promoter, exon, intron, intergenic regions, TSSs), and of CpG classification (CpG islands, shores, and other regions), was conducted using the package Genomation [64]. Reference.bed files were downloaded from the table browser UCSC [65], containing only RefSeq curated genes for functional features, and default CpG islands classification, both based on mm10 genome assembly.

Promoter regions for annotation of genes were defined as within the range from -500 to + 500 base pairs of distance from TSS, and minimum of three DMCs were considered for accounting the gene as having its promoter differentially methylated.

From the annotated genes with differentially methylated promoters, a selection was made based on the cell lines’ characteristics. Hypomethylated and hypermethylated genes present in the intersection among the comparisons between each of the cell lines 4C, 4C11− and 4C11 + with melan-a were gathered as being related to melanocyte malignant transformation, whereas genes present in the intersection among the comparisons between all non-metastatic cell lines (melan-a, 4C and 4C11−) with the aggressive lineage 4C11 + were accounted as related to melanoma metastasis/aggressiveness in the model. Venn Diagrams were constructed using the BioVenn tool [66].

5-Methylcytosine content

Nuclear DNA (6 mg), extracted as described in the ERRBS method, was transferred to 54.5 µL of Milli-Q water containing 2.5 µL of 200 mM Tris/MgCl2 buffer (pH 7.4) and 1 µL (2.5 units) of deoxyribonuclease I (Qiagen). The samples were incubated at 37 °C for 1 h. Then, 1 µL (0.001 units) of phosphodiesterase I and 1 µL (2 units) of alkaline phosphatase were added, followed by further incubation at 37 °C for 1 h. Aliquots of 50 µL of residual volume were analyzed by a HPLC/PDA system (Shimadzu), as follows: a 250 mm × 4.6 mm i.d., 5 µm, Luna C18(2) column (Phenomenex) was eluted with a gradient of 0.1% formic acid in water (Solution A) and CH3OH:H2O (1:1, v/v; added of 0.1% formic acid; Solution B; from 0 to 25 min, 0 to 36% B; from 25 to 27 min, 36 to 0% B; from 27 to 37 min, 0% B) at a 1 mL/min flow rate and 30 °C. The PDA detector was set at 286 nm. Calibration curves were constructed at intervals of 0.5 to 8 nmol for 2′-deoxycytidine (dC) and 0.01 to 0.8 nmol for 5-mC. The percentage of global DNA methylation was calculated using the following equation:

$${5mC}_{\left(\%\right)}=\frac{{5mC}_{\left(nmol\right)}x100}{{5mC}_{\left(nmol\right)}+{dC}_{\left(nmol\right)}}$$

Pathway enrichment

The gene enrichment analysis in the database Gene Ontology: Biological Process of the overlapped genes involved in malignant transformation and melanoma aggressiveness was conducted using the R package clusterProfiler [30]. The enrichment analysis is based on hypergeometric distribution, and it calculates q-value for FDR control.

Survival analyses per CpG

The Cox proportional-hazard models for overall survival were built using TCGA-SKCM data for all the CpGs annotated to orthologs of the differentially methylated genes from the malignant transformation and melanoma metastasis signatures, using the methSurv tools developed by Modhukur et al. [44]. In methSurv, the analysis is made by individual CpG sites, and patient subgroups are separated based on methylation status into high and low methylation according to distinct parameters, such as median, mean, first quartile (q25) and fourth quartile (q75). Each model contains the covariables: age, sex, and a single CpG methylation. Kaplan–Meier plots of selected CpGs were built using the website provided by methSurv authors, which is https://biit.cs.ut.ee/methsurv/. All plots were generated splitting the samples by the best method, with log-rank p-value accounted, and adjusted for the covariates age and sex. All models with adjusted p-values ≤ 0.01 are represented in the Additional file 4: Table S4.

Survival analyses with multiple CpGs

The CpGs whose methylation patterns were significantly associated with survival in the Cox models (single CpG per model) for each signature were selected and combined in a new multivariate model with multiple CpGs. LASSO regularization was performed in the Cox proportional-hazards models for selecting variables and avoiding overfitting. The glmnet R package was used, with minimum λ obtained after 20 simulations. After regularization, a multivariate Cox regression was performed, using the survival R package and CpGs methylation, gender, age, and Breslow depth as covariates. The coefficients from the significant CpGs in the Cox models were used in to obtain a risk score for each patient: \({\alpha }_{j}=\sum {\varphi }_{i}\times {\beta }_{ij}\), where φ is the coefficient in the Cox model for the ith CpG, β is the methylation β value for the ith CpG in the jth patient, and α is the risk score for the jth patient. The samples were split based on the mean risk score value in high and low risk subgroups, and Kaplan–Meier plots were built for malignancy and metastasis signatures using survminer R package. All p-values were corrected for multiple-testing hypotheses using the Benjamini-Hochberg (False Discovery Rate) method. 345 patients from TCGA-SKCM were analyzed, after removing individuals with missing data.

RNA expression

RNA sequencing of the cell lines melan-a, 4C, 4C11− and 4C11+ was performed in triplicates in the Illumina HiSeq 1500 platform. Analysis of the data was conducted by trimming of the fragments using Trimmomatic [67], STAR alignment [68], and Rsubread [69] counting of the fragments, to assess the expression level of RNA in each of the cell lines. A complete description of the RNA-seq analysis can be found in Pessoa and colleagues [16].

Availability of data and materials

CpG methylation calls are available at the following link: https://github.com/flaviaerius/methylation-data/tree/master/input/raw-methylation-data.

References

  1. Trotter SC, Sroa N, Winkelmann RR, Olencki T, Bechtel M. A global review of melanoma follow-up guidelines. J Clin Aesthet Dermatol. 2013;6:18–26.

    PubMed  PubMed Central  Google Scholar 

  2. Matthews NH, Li W-Q, Qureshi AA, Weinstock MA, Cho E. In: Ward WH, Farma JM, editors. Cutaneous melanoma: etiology and therapy. Brisbane: Codon Publications; 2017. p. 3–17. https://doi.org/10.15586/codon.cutaneousmelanoma.2017.ch1.

    Chapter  Google Scholar 

  3. Clark J, Wallace H, Elder DE, Guerry D, Epstein MN, Greene MH, Van Horn M. A study of tumor progression: the precursor lesions of superficial spreading and nodular melanoma. Hum Pathol. 1984;15:1147–65. https://doi.org/10.1016/S0046-8177(84)80310-X.

    Article  PubMed  Google Scholar 

  4. Elder DE. Melanoma progression. Pathology. 2016;48(2):147–54. https://doi.org/10.1016/j.pathol.2015.12.002.

    Article  PubMed  Google Scholar 

  5. Jones PA. Functions of DNA methylation: Islands, start sites, gene bodies and beyond. Nat Rev Genet. 2012;13(7):484–92. https://doi.org/10.1038/nrg3230.

    Article  CAS  PubMed  Google Scholar 

  6. Pan Y, Liu G, Zhou F, Su B, Li Y. DNA methylation profiles in cancer diagnosis and therapeutics. Clin Exp Med. 2018;18:1–14. https://doi.org/10.1007/s10238-017-0467-0.

    Article  CAS  PubMed  Google Scholar 

  7. Schinke C, Mo Y, Yu Y, Amiri K, Sosman J, Greally J, Verma A. Aberrant DNA methylation in malignant melanoma. Melanoma Res. 2010;20:253–65. https://doi.org/10.1097/CMR.0b013e328338a35a.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  8. Fu S, Wu H, Zhang H, Lian CG, Lu Q. DNA methylation/hydroxymethylation in melanoma. Oncotarget. 2017;8(44):78163–73. https://doi.org/10.18632/oncotarget.18293.

    Article  PubMed  PubMed Central  Google Scholar 

  9. Charlet J, Liang G. Epigenetics of melanoma. Epigenet Dermatol. 2015;2015(3):339–54.

    Article  Google Scholar 

  10. Koroknai V, Szász I, Hernandez-Vargas H, Fernandez-Jimenez N, Cuenin C, Herceg Z, Vízkeleti L, Ádány R, Ecsedi S, Balázs M. DNA hypermethylation is associated with invasive phenotype of malignant melanoma. Exp Dermatol. 2020;29(1):39–50. https://doi.org/10.1111/exd.14047.

    Article  CAS  PubMed  Google Scholar 

  11. Rodger EJ, Almomani SN, Ludgate JL, Stockwell PA, Baguley BC, Eccles MR, Chatterjee A. Comparison of global DNA methylation patterns in human melanoma tissues and their derivative cell lines. Cancers (Basel). 2021;13(9):2123. https://doi.org/10.3390/cancers13092123.

    Article  CAS  Google Scholar 

  12. Simiczyjew A, Dratkiewicz E, Mazurkiewicz J, Ziętek M, Matkowski R, Nowak D. The influence of tumor microenvironment on immune escape of melanoma. Int J Mol Sci. 2020;21:1–28. https://doi.org/10.3390/ijms21218359.

    Article  CAS  Google Scholar 

  13. Li H, van der Leun AM, Yofe I, Lubling Y, Gelbard-Solodkin D, van Akkooi ACJ, van den Braber M, Rozeman EA, Haanen JBAJ, Blank CU, Horlings HM, David E, Baran Y, Bercovich A, Lifshitz A, Schumacher TN, Tanay A, Amit I. Dysfunctional CD8 T cells form a proliferative, dynamically regulated compartment within human melanoma. Cell. 2019;176(4):775-789.e18. https://doi.org/10.1016/j.cell.2018.11.043.

    Article  CAS  PubMed  Google Scholar 

  14. Lähnemann D, Köster J, Szczurek E, McCarthy DJ, Hicks SC, Robinson MD, Vallejos CA, Campbell KR, Beerenwinkel N, Mahfouz A, Pinello L, Skums P, Stamatakis A, Attolini CS, Aparicio S, Baaijens J, Balvert M, Barbanson B, Cappuccio A, Corleone G, Dutilh BE, Florescu M, Guryev V, Holmer R, Jahn K, Lobo TJ, Keizer EM, Khatri I, Kielbasa SM, Korbel JO, Kozlov AM, Kuo TH, Lelieveldt BPF, Mandoiu II, Marioni JC, Marschall T, Mölder F, Niknejad A, Raczkowski L, Reinders M, Ridder J, Saliba AE, Somarakis A, Stegle O, Theis FJ, Yang H, Zelikovsky A, McHardy AC, Raphael BJ, Shah SP, Schönhuth A. Eleven grand challenges in single-cell data science. Genome Biol. 2020;21(1):31. https://doi.org/10.1186/s13059-020-1926-6.

    Article  PubMed  PubMed Central  Google Scholar 

  15. Oba-Shinjo SM, Correa M, Ricca TI, Molognoni F, Pinhal MA, Neves IA, Marie SK, Sampaio LO, Nader HB, Chammas R, Jasiulionis MG. Melanocyte transformation associated with substrate adhesion impediment. Neoplasia. 2006;8(3):231–41. https://doi.org/10.1593/neo.05781.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  16. Pessoa DO, Rius FE, Papaiz DD, Ayub ALP, Morais AS, de Souza CF, da Paixão VF, Setubal JC, Newton-Bishop J, Nsengimana J, Azevedo H, Reis EM, Jasiulionis MG. Transcriptional signatures underlying dynamic phenotypic switching and novel disease biomarkers in a linear cellular model of melanoma progression. Neoplasia. 2021;23(4):439–55. https://doi.org/10.1016/j.neo.2021.03.007.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  17. Toricelli M, Melo FH, Peres GB, Silva DC, Jasiulionis MG. Timp1 interacts with beta-1 integrin and CD63 along melanoma genesis and confers anoikis resistance by activating PI3-K signaling pathway independently of Akt phosphorylation. Mol Cancer. 2013;12:22. https://doi.org/10.1186/1476-4598-12-22.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  18. Monteiro AC, Muenzner JK, Andrade F, Rius FE, Ostalecki C, Geppert CI, Agaimy A, Hartmann A, Fujita A, Schneider-Stock R, Jasiulionis MG. Gene expression and promoter methylation of angiogenic and lymphangiogenic factors as prognostic markers in melanoma. Mol Oncol. 2019;13(6):1433–49. https://doi.org/10.1002/1878-0261.12501.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  19. da Cruz AT, Hunger A, de Melo FHM, Monteiro AC, Paré GC, Lai D, Alves-Fernandes DK, Ayub ALP, Cordero EM, Filho JFDS, Schneider-Stock R, Strauss BE, Tron V, Jasiulionis MG. miR-138-5p induces aggressive traits by targeting Trp53 expression in murine melanoma cells, and correlates with poor prognosis of melanoma patients. Neoplasia. 2021;23(8):823–34. https://doi.org/10.1016/j.neo.2021.05.015.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  20. de Souza CF, Xander P, Monteiro AC, Silva AG, da Silva DC, Mai S, Bernardo V, Lopes JD, Jasiulionis MG. Mining gene expression signature for the detection of pre-malignant melanocytes and early melanomas with risk for metastasis. PLoS ONE. 2012;7(9): e44800. https://doi.org/10.1371/journal.pone.0044800.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  21. Molognoni F, Cruz AT, Meliso FM, Morais AS, Souza CF, Xander P, Bischof JM, Costa FF, Soares MB, Liang G, Jones PA, Jasiulionis MG. Epigenetic reprogramming as a key contributor to melanocyte malignant transformation. Epigenetics. 2011;6(4):450–64. https://doi.org/10.4161/epi.6.4.14917.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  22. Melo FHM, Molognoni F, Morais AS, Toricelli M, Mouro MG, Higa EM, Lopes JD, Jasiulionis MG. Endothelial nitric oxide synthase uncoupling as a key mediator of melanocyte malignant transformation associated with sustained stress conditions. Free Radic Biol Med. 2011;50(10):1263–73. https://doi.org/10.1016/j.freeradbiomed.2011.02.022.

    Article  CAS  PubMed  Google Scholar 

  23. Preston-Alp S, Jelinek J, Issa J-P, Zaidi MR. Ultraviolet radiation modulates DNA methylation in melanocytes. BioRxiv. 2021. https://doi.org/10.1101/2021.10.14.464470.

    Article  Google Scholar 

  24. Garrett-Bakelman FE, Sheridan CK, Kacmarczyk TJ, Ishii J, Betel D, Alonso A, Mason CE, Figueroa ME, Melnick AM. Enhanced reduced representation bisulfite sequencing for assessment of DNA methylation at base pair resolution. J Vis Exp. 2015;96: e52246. https://doi.org/10.3791/52246.

    Article  CAS  Google Scholar 

  25. Bennett DC, Cooper PJ, Hart IR. A line of non-tumorigenic mouse melanocytes, syngeneic with the B16 melanoma and requiring a tumour promoter for growth. Int J Cancer. 1987;39(3):414–8. https://doi.org/10.1002/ijc.2910390324.

    Article  CAS  PubMed  Google Scholar 

  26. Krueger F, Andrews SR. Bismark: a flexible aligner and methylation caller for Bisulfite-Seq applications. Bioinformatics. 2011;27(11):1571–2. https://doi.org/10.1093/bioinformatics/btr167.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  27. Akalin A, Kormaksson M, Li S, Garrett-Bakelman FE, Figueroa ME, Melnick A, Mason CE. methylKit: a comprehensive R package for the analysis of genome-wide DNA methylation profiles. Genome Biol. 2012;13(10):R87. https://doi.org/10.1186/gb-2012-13-10-r87.

    Article  PubMed  PubMed Central  Google Scholar 

  28. R Core Team. R: a language and environment for statistical computing. Vienna: R Foundation for Statistical Computing. 2017; https://www.R-project.org/

  29. Wang HQ, Tuominen LK, Tsai CJ. SLIM: A sliding linear model for estimating the proportion of true null hypotheses in datasets with dependence structures. Bioinformatics. 2011;27(2):225–31. https://doi.org/10.1093/bioinformatics/btq650.

    Article  CAS  PubMed  Google Scholar 

  30. Yu G, Wang L-G, Han Y, He Q-Y. ClusterProfiler: an R package for comparing biological themes among gene clusters. OMICS. 2012;16(5):284–7. https://doi.org/10.1089/omi.2011.0118.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  31. Jjingo D, Conley AB, Yi SV, Lunyak VV, Jordan IK. On the presence and role of human gene-body DNA methylation. Oncotarget. 2012;3(4):462–74. https://doi.org/10.18632/oncotarget.497.

    Article  PubMed  PubMed Central  Google Scholar 

  32. Fan L, Zhang A, Deng P. LMX1B mRNA expression and its gene body CpG methylation are valuable prognostic biomarkers for laryngeal squamous cell carcinoma. Biomed Pharmacother. 2019;117: 109174. https://doi.org/10.1016/j.biopha.2019.109174.

    Article  CAS  PubMed  Google Scholar 

  33. Patsouras MD, Karagianni P, Kogionou P, Vlachoyiannopoulos PG. Differential CpG methylation of the promoter of interleukin 8 and the first intron of tissue factor in Antiphospholipid syndrome. J Autoimmun. 2019;102:159–66. https://doi.org/10.1016/j.jaut.2019.05.001.

    Article  CAS  PubMed  Google Scholar 

  34. De Almeida BP, Apolónio JD, Binnie A, Castelo-Branco P. Roadmap of DNA methylation in breast cancer identifies novel prognostic biomarkers. BMC Cancer. 2019;19(1):219. https://doi.org/10.1186/s12885-019-5403-0.

    Article  PubMed  PubMed Central  Google Scholar 

  35. Hage M, Chaligné R, Viengchareun S, Villa C, Salenave S, Bouligand J, Letouzé E, Tosca L, Rouquette A, Tachdjian G, Parker F, Lombès M, Lacroix A, Gaillard S, Chanson P, Kamenický P. Hypermethylator phenotype and ectopic GIP receptor in GNAS mutation-negative somatotropinomas. J Clin Endocrinol Metab. 2019;104(5):1777–87. https://doi.org/10.1210/jc.2018-01504.

    Article  PubMed  Google Scholar 

  36. Jones PA, Baylin SB. The epigenomics of cancer. Cell. 2007;128(4):683–92. https://doi.org/10.1016/j.cell.2007.01.029.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  37. Esteller M. Cancer epigenomics: DNA methylomes and histone-modification maps. Nat Rev Genet. 2007;8(4):286–98. https://doi.org/10.1038/nrg2005.

    Article  CAS  PubMed  Google Scholar 

  38. Baylin SB, Jones PA. Epigenetic determinants of cancer. Cold Spring Harb Perspect Biol. 2016;8: a019505. https://doi.org/10.1101/cshperspect.a019505.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  39. Ecsedi S, Hernandez-Vargas H, Lima SC, Vizkeleti L, Toth R, Lazar V, Koroknai V, Kiss T, Emri G, Herceg Z, Adany R, Balazs M. DNA methylation characteristics of primary melanomas with distinct biological behaviour. PLoS ONE. 2014;9(5): e96612. https://doi.org/10.1371/journal.pone.0096612.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  40. Li J-L, Mazar J, Zhong C, Govindarajan SS, Zhang Z, Dinger ME, Meredith G, Adams C, Zhang S, Mattick JS, Ray A, Perera RJ. Genome-wide methylated CpG island profiles of melanoma cells reveal a melanoma coregulation network. Sci Rep. 2013;3:2962. https://doi.org/10.1038/srep02962.

    Article  PubMed  PubMed Central  Google Scholar 

  41. Lauss M, Ringnér M, Karlsson A, Harbst K, Busch C, Geisler J, Lønning PE, Staaf J, Jönsson G. DNA methylation subgroups in melanoma are associated with proliferative and immunological processes. BMC Med Genomics. 2015;8:73. https://doi.org/10.1186/s12920-015-0147-4.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  42. Wouters J, Vizoso M, Martinez-Cardus A, Carmona FJ, Govaere O, Laguna T, Joseph J, Dynoodt P, Aura C, Foth M, Cloots R, van den Hurk K, Balint B, Murphy IG, McDermott EW, Sheahan K, Jirström K, Nodin B, Mallya-Udupi G, van den Oord JJ, Gallagher WM, Esteller M. Comprehensive DNA methylation study identifies novel progression-related and prognostic markers for cutaneous melanoma. BMC Med. 2017;15(1):101. https://doi.org/10.1186/s12916-017-0851-3.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  43. Zhang Z-M, Lu R, Wang P, Yu Y, Chen D, Gao L, Liu S, Ji D, Rothbart SB, Wang Y, Wang GG, Song J. Structural basis for DNMT3A-mediated de novo DNA methylation. Nature. 2018;554(7692):387–91. https://doi.org/10.1038/nature25477.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  44. Modhukur V, Iljasenko T, Metsalu T, Lokk K, Laisk-Podar T, Vilo J. MethSurv: a web tool to perform multivariable survival analysis using DNA methylation data. Epigenomics. 2018;10(3):277–88. https://doi.org/10.2217/epi-2017-0118.

    Article  CAS  PubMed  Google Scholar 

  45. de Unamuno BB, Murria Estal R, Pérez Simó G, Simarro Farinos J, Pujol Marco C, Navarro Mira M, Alegrede Miquel V, Ballester Sánchez R, Sabater Marco V, Llavador Ros M, Palanca Suela S, Botella ER. Aberrant DNA methylation is associated with aggressive clinicopathological features and poor survival in cutaneous melanoma. Br J Dermatol. 2018;179(2):394–404. https://doi.org/10.1111/bjd.16254.

    Article  CAS  Google Scholar 

  46. Tanemura A, Terando AM, Sim MS, van Hoesel AQ, de Maat MF, Morton DL, Hoon DS. CpG Island methylator phenotype predicts progression of malignant melanoma. Clin Cancer Res. 2009;15(5):1801–7. https://doi.org/10.1158/1078-0432.CCR-08-1361.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  47. Dai J, Lin K, Huang Y, Lu Y, Chen WQ, Zhang XR, He BS, Pan YQ, Wang SK, Fan WX. Identification of critically carcinogenesis-related genes in basal cell carcinoma. Onco Targets Ther. 2018;11:6957–67. https://doi.org/10.2147/OTT.S170504.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  48. Wang T, Xie X, Liu H, Chen F, Du J, Wang X, Jiang XY, Yu F, Fan H. Pyridine nucleotide-disulphide oxidoreductase domain 2 (PYROXD2): role in mitochondrial function. Mitochondrion. 2019;47:114–24. https://doi.org/10.1016/j.mito.2019.05.007.

    Article  CAS  PubMed  Google Scholar 

  49. Aguila B, Morris AB, Spina R, Bar E, Schraner J, Vinkler R, Sohn JW, Welford SM. The Ig superfamily protein PTGFRN coordinates survival signaling in glioblastoma multiforme. Cancer Lett. 2019;462:33–42. https://doi.org/10.1016/j.canlet.2019.07.018.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  50. Hao N, Bhakti VLD, Peet DJ, Whitelaw ML. Reciprocal regulation of the basic helix-loop-helix/Per-Arnt-Sim partner proteins, Arnt and Arnt2, during neuronal differentiation. Nucleic Acids Res. 2013;41(11):5626–38. https://doi.org/10.1093/nar/gkt206.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  51. Li W, Liu J, Zhang B, Bie Q, Qian H, Xu W. Transcriptome analysis reveals key genes and pathways associated with metastasis in breast cancer. Onco Targets Ther. 2020;13:323–35. https://doi.org/10.2147/OTT.S226770.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  52. Bogeas A, Morvan-Dubois G, El-Habr EA, Lejeune FX, Defrance M, Narayanan A, Kuranda K, Burel-Vandenbos F, Sayd S, Delaunay V, Dubois LG, Parrinello H, Rialle S, Fabrega S, Idbaih A, Haiech J, Bièche I, Virolle T, Goodhardt M, Chneiweiss H, Junier MP. Changes in chromatin state reveal ARNT2 at a node of a tumorigenic transcription factor signature driving glioblastoma cell aggressiveness. Acta Neuropathol. 2018;135(2):267–83. https://doi.org/10.1007/s00401-017-1783-x.

    Article  CAS  PubMed  Google Scholar 

  53. Dallinga MG, Habani YI, Kayser RP, Van Noorden CJF, Klaassen I, Schlingemann RO. IGF-binding proteins 3 and 4 are regulators of sprouting angiogenesis. Mol Biol Rep. 2020;47(4):2561–72. https://doi.org/10.1007/s11033-020-05339-0.

    Article  CAS  PubMed  Google Scholar 

  54. Yu JZ, Warycha MA, Christos PJ, Darvishian F, Yee H, Kaminio H, Berman RS, Shapiro RL, Buckley MT, Liebes LF, Pavlick AC, Polsky D, Brooks PC, Osman I. Assessing the clinical utility of measuring Insulin-like Growth Factor Binding Proteins in tissues and sera of melanoma patients. J Transl Med. 2008;6:70. https://doi.org/10.1186/1479-5876-6-70.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  55. Praveen Kumar VR, Sehgal P, Thota B, Praveen Kumar VR. Insulin-like growth factor binding protein 4 promotes GBM progression and regulates key factors involved in EMT and invasion. J Neurooncol. 2014;116(3):455–64. https://doi.org/10.1007/s11060-013-1324-y.

    Article  CAS  PubMed  Google Scholar 

  56. Xiao Y, Zhu S, Yin W, Liu X, Hu Y. IGFBP-4 expression is adversely associated with lung cancer prognosis. Oncol Lett. 2017;14(6):6876–80. https://doi.org/10.3892/ol.2017.7014.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  57. Liu W, Peng Y, Tobin DJ. A new 12-gene diagnostic biomarker signature of melanoma revealed by integrated microarray analysis. PeerJ. 2013;1: e49. https://doi.org/10.7717/peerj.49.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  58. Scrima M, De Marco C, De Vita F, Fabiani F, Franco R, Pirozzi G, Rocco G, Malanga D, Viglietto G. The nonreceptor-type tyrosine phosphatase PTPN13 is a tumor suppressor gene in nonsmall cell lung cancer. Am J Pathol. 2012;180(3):1202–14. https://doi.org/10.1016/j.ajpath.2011.11.038.

    Article  CAS  PubMed  Google Scholar 

  59. Wei JH, Haddad A, Wu KJ, Zhao HW, Kapur P, Zhang ZL, Zhao LY, Chen ZH, Zhou YY, Zhou JC, Wang B, Yu YH, Cai MY, Xie D, Liao B, Li CX, Li PX, Wang ZR, Zhou FJ, Shi L, Liu QZ, Gao ZL, He DL, Chen W, Hsieh JT, Li QZ, Margulis V, Luo JH. A CpG-methylation-based assay to predict survival in clear cell renal cell carcinoma. Nat Commun. 2015;6:8699. https://doi.org/10.1038/ncomms9699.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  60. Jiang H, Ou Z, He Y, Yu M, Wu S, Li G, Zhu J, Zhang R, Wang J, Zheng L, Zhang X, Hao W, He L, Gu X, Quan Q, Zhang E, Luo H, Wei W, Li Z, Zang G, Zhang C, Poon T, Zhang D, Ziyar I, Zhang RZ, Li O, Cheng L, Shimizu T, Cui X, Zhu JK, Sun X, Zhang K. DNA methylation markers in the diagnosis and prognosis of common leukemias. Signal Transduct Target Ther. 2020;5(1):3. https://doi.org/10.1038/s41392-019-0090-5.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  61. Kuo IY, Chang JM, Jiang SS, Chen CH, Chang IS, Sheu BS, Lu PJ, Chang WL, Lai WW, Wang YC. Prognostic CpG methylation biomarkers identified by methylation array in esophageal squamous cell carcinoma patients. Int J Med Sci. 2014;11(8):779–87. https://doi.org/10.7150/ijms.7405.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  62. Chung HH, Kuo CC, Hsiao CW, Chen CY, Hu JM, Hsu CH, Chou YC, Lin YW, Shih YL. A novel prognostic DNA methylation panel for colorectal cancer. Int J Mol Sci. 2019;20(19):4672. https://doi.org/10.3390/ijms20194672.

    Article  CAS  PubMed Central  Google Scholar 

  63. Yan P, Yang X, Wang J, Wang S, Ren H. A novel CpG island methylation panel predicts survival in lung adenocarcinomas. Oncol Lett. 2019;18(2):1011–22. https://doi.org/10.3892/ol.2019.10431.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  64. Akalin A, Franke V, Vlahoviček K, Mason CE, Schübeler D. Genomation: a toolkit to summarize, annotate and visualize genomic intervals. Bioinformatics. 2015;31(7):1127–9. https://doi.org/10.1093/bioinformatics/btu775.

    Article  CAS  PubMed  Google Scholar 

  65. Raney BJ, Dreszer TR, Barber GP, Clawson H, Fujita PA, Wang T, Nguyen N, Paten B, Zweig AS, Karolchik D, Kent WJ. Track data hubs enable visualization of user-defined genome-wide annotations on the UCSC Genome Browser. Bioinformatics. 2014;30(7):1003–5. https://doi.org/10.1093/bioinformatics/btt637.

    Article  CAS  PubMed  Google Scholar 

  66. Hulsen T, de Vlieg J, Alkema W. BioVenn—a web application for the comparison and visualization of biological lists using area-proportional Venn diagrams. BMC Genomics. 2008;9:488. https://doi.org/10.1186/1471-2164-9-488.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  67. Usadellab. Trimmomatic Manual: V0.32. http://www.usadellab.org/cms/uploads/supplementary/Trimmomatic/TrimmomaticManual_V0.32.pdf.

  68. Dobin A, Davis CA, Schlesinger F, Drenkow J, Zaleski C, Jha S, Batut P, Chaisson M, Gingeras TR. STAR: ultrafast universal RNA-seq aligner. Bioinformatics. 2013;29(1):15–21. https://doi.org/10.1093/bioinformatics/bts635.

    Article  CAS  PubMed  Google Scholar 

  69. Liao Y, Smyth GK, Shi W. The R package Rsubread is easier, faster, cheaper and better for alignment and quantification of RNA sequencing reads. Nucleic Acids Res. 2019;47(8): e47. https://doi.org/10.1093/nar/gkz114.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

Download references

Acknowledgements

We acknowledge all specimen donors of TCGA-SKCM, as well as the research groups involved in the TCGA project. We also would like to thank the very contributive advice of Dr. Cem Meydan and Prof. Dr. Ari Melnick concerning the analyses of ERRBS data. Furthermore, we acknowledge Epigenomics Core for their contribution with the conduction and analyses of the ERRBS experiments.

Funding

This work was supported by the São Paulo Research Foundation (FAPESP) under Grants 2014/13663-0, 2018/20775-0, 2019/23480-4, and 2022/00322-7; Conselho Nacional de Desenvolvimento Científico e Tecnológico (CNPq) under Grant 140594/2014-0; and Coordenação de Aperfeiçoamento de Pessoal de Nível Superior (CAPES) under Grant 88881.132045/2016-01.

Author information

Authors and Affiliations

  1. Departamento de Farmacologia, Escola Paulista de Medicina, Universidade Federal de São Paulo, São Paulo, Brazil

    Flávia E. Rius, Debora D. Papaiz & Ana Luísa P. Ayub

  2. Divisão de Urologia, Departamento de Cirurgia, Universidade Federal de São Paulo, São Paulo, Brazil

    Hatylas F. Z. Azevedo

  3. Departamento de Bioquímica, Instituto de Química, Universidade de São Paulo, São Paulo, Brazil

    Diogo O. Pessoa & Eduardo M. Reis

  4. Departamento de Farmacociências, Universidade Federal de Ciências da Saúde de Porto Alegre, São Paulo, Brazil

    Tiago F. Oliveira

  5. Departamento de Análises Clínicas e Toxicológicas, Faculdade de Ciências Farmacêuticas, Universidade de São Paulo, São Paulo, Brazil

    Tiago F. Oliveira & Ana Paula M. Loureiro

  6. Bioinformatics Graduate Program, Instituto de Matemática e Estatística, Universidade de São Paulo, São Paulo, Brazil

    Fernando Andrade

  7. Department of Biology, Loyola University Chicago, Chicago, USA

    Fernando Andrade

  8. Departamento de Ciências da Computação, Instituto de Matemática e Estatística, Universidade de São Paulo, São Paulo, Brazil

    André Fujita

  9. Department of Physiology and Biophysics, Weill Cornell Medicine, New York, USA

    Christopher E. Mason

  10. Departamento de Farmacologia, Escola Paulista de Medicina, Universidade Federal de São Paulo, São Paulo, 04039-032, Brazil

    Miriam G. Jasiulionis

Authors
  1. Flávia E. Rius
    View author publications

    You can also search for this author in PubMed Google Scholar

  2. Debora D. Papaiz
    View author publications

    You can also search for this author in PubMed Google Scholar

  3. Hatylas F. Z. Azevedo
    View author publications

    You can also search for this author in PubMed Google Scholar

  4. Ana Luísa P. Ayub
    View author publications

    You can also search for this author in PubMed Google Scholar

  5. Diogo O. Pessoa
    View author publications

    You can also search for this author in PubMed Google Scholar

  6. Tiago F. Oliveira
    View author publications

    You can also search for this author in PubMed Google Scholar

  7. Ana Paula M. Loureiro
    View author publications

    You can also search for this author in PubMed Google Scholar

  8. Fernando Andrade
    View author publications

    You can also search for this author in PubMed Google Scholar

  9. André Fujita
    View author publications

    You can also search for this author in PubMed Google Scholar

  10. Eduardo M. Reis
    View author publications

    You can also search for this author in PubMed Google Scholar

  11. Christopher E. Mason
    View author publications

    You can also search for this author in PubMed Google Scholar

  12. Miriam G. Jasiulionis
    View author publications

    You can also search for this author in PubMed Google Scholar

Contributions

FER, HA, CM and MGJ conceived and designed the study. FER, HA and MGJ interpreted the data, and participated in manuscript preparation. FER performed most experiments. FER and DP drafted the manuscript. FER, ALPA, DP and MGJ prepared the figures and tables. ALPA performed epigenetic drugs treatment. TFO and APML contributed to 5mC quantification by HPLC. DOP and EMR performed the RNA sequencing. FER carried out the methylome analysis. FER, HA, FA and AF performed the survival analysis. All authors read and approved the final manuscript.

Corresponding author

Correspondence to Miriam G. Jasiulionis.

Ethics declarations

Ethics approval and consent to participate

This work was approved by the Ethics Committee of Federal University of São Paulo—Brazil (Comitê de Ética e Pesquisa—Unifesp) under the identification 7004230514. Consent to participate: Not applicable.

Consent for publication

Not applicable.

Competing interests

The authors report no competing interests.

Additional information

Publisher's Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Supplementary Information

Additional file 1: Table S1.

Complete list of differentially methylated CpGs with the respective genes for each pairwise comparison.

Additional file 2: Table S2.

Genes with differentially methylated promoters in the malignancy and metastasis signatures.

Additional file 3: Table S3.

Characteristics of the melanoma patients from TCGA-SKCM selected for the Cox proportional-hazards models.

Additional file 4: Table S4.

Cox proportional-hazards models of significant hypermethylated and hypomethylated CpGs individually from malignancy and metastasis signatures.

Additional file 5: Table S5.

Complete list of covariates included in the final Cox survival model for malignancy and metastasis signatures with their respective metrics.

Additional file 6: Table S6.

Equivalence list between mouse ERRBS CpGs location and human Infinium 450 k CpGs probes and location referenced by their closest gene.

Additional file 7: Figure S1.

Distribution of different percentages of methylation across all CpGs analyzed in a sample. Three biological replicates of each cell line have their methylation levels represented.

Additional file 8: Figure S2.

Heatmaps with differentially methylated CpGs between cell lines in the melanoma progression model. Results from pairwise comparisons between non-malignant melan-a melanocytes and non-metastatic (4C11−) and metastatic (4C11+) melanoma cell lines, between 4 and 4C11− and 4C11+, and between 4C11− and 4C11+ are shown.

Additional file 9: Figure S3.

Pyroxd2 and Ptgfrn gene silencing is not reverted in melan-a melanocytes after 5azaCdR treatment. Melan-a cells were treated with 5azaCdR for 48 h, and the expressions of Pyroxd2 and Ptgfrn were determined by RT-qPCR.

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Rius, F.E., Papaiz, D.D., Azevedo, H.F.Z. et al. Genome-wide promoter methylation profiling in a cellular model of melanoma progression reveals markers of malignancy and metastasis that predict melanoma survival. Clin Epigenet 14, 68 (2022). https://doi.org/10.1186/s13148-022-01291-x

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doi.org/10.1186/s13148-022-01291-x

Keywords

Clinical Epigenetics

ISSN: 1868-7083

Contact us