2.1 Screening and establishment of prognosis-related CRG DNA methylation signature

A total of 96 CRGs were collected and 1374 DNA methylation sites located in CRGs were extracted. Epigenetic modifications have been demonstrated to affect RNA content and in turn plays a key role in the development of tumour. Therefore, we calculated the spearman correlation coefficient between DNA methylation sites and corresponding genes, and screened 695 methylation sites corresponding to and significantly associated with CRGs (p < 0.05). These DNA methylation sites were used in a multivariate Cox regression analysis to construct all possible models of 1–5 sites, a CRG DNA methylation signature was selected as the best prognostic model for predicting survival through ROC analyses. The p value of the likelihood ratio test was 4e−09 and the p value of the cox proportional hazards hypothesis test was 0.066 for the model. The final model for predicting OS consisted of five DNA methylation sites (cg00668852, cg01130192, cg18645035, cg22134611, cg25522181). The genes corresponding to these five sites are AOC1 (Amine oxidase copper-containing 1), ACO1 (Aconitase 1), DPYD (Dihydropyrimidine Dehydrogenas), ABCB8 (ATP Binding Cassette Subfamily B Member 8), and CIAO1 (Cytosolic Iron-Sulfur Assembly Component 1). Relevant information about these sites is shown in Table 1.

Table 1 Information of five CRG-located DNA-methylation sites2.2 Evaluation of the CRG DNA methylation prognosis signature

The risk scores of CM patients were calculated according to the model: \(} = \left( }\;}\;}00668852} \right) + \left( }\;}\;}01130192} \right)\) \(+ \left( }\;}\;}18645035} \right) + \left( }\;}\;}22134611 } \right)\) \(+ \left( }\;}\;}25522181} \right).\) CM patients were separated into high- and low-risk groups based on the median risk score.

In the TCGA cohort, Kaplan–Meier survival analysis revealed a significantly lower survival rate in the high-risk group (p < 0.001, Fig. 1A). The receiver operating characteristic (ROC) curve analysis exhibited that the signature had a good predictive accuracy with an area under the curve (AUC) value of 0.742 and a Confidence interval (CI) of 0.687–0.798 (Fig. 1B).

Fig. 1

Identification of CRG DNA methylation prognosis signature in CM patients. A, C Kaplan–Meier curves for high-risk and low-risk groups in TCGA cohort and GEO cohort. B, D ROC curves of CRG DNA methylation prognosis signature in TCGA cohort and GEO cohort

To further explore the prognostic value of this signature, we discovered in the independent GEO cohort that patients with a high-risk score had a significantly worse survival prognosis (p = 0.007, Fig. 1C), with an AUC of 0.733 (Fig. 1D). This finding validates the excellent predictive performance of this prognostic signature for CM patients.

In addition, univariate Cox regression analysis, Kaplan–Meier and ROC analyses were performed on the five DNA methylation sites in TCGA (Figure S1) and GEO (Figure S2) cohorts to further investigate their prognostic significance in CM. Kaplan–Meier analysis revealed that individual DNA methylation sites could also differentiate high-risk from low-risk patients. ROC analysis indicated that the prognostic performance of DNA methylation site combinations was superior to that of individual DNA methylation sites alone. These findings suggest that DNA methylation site combinations can provide a better prognostic prediction for CM patients.

2.3 Correlation between CRG DNA methylation signature and clinical characteristics

To assess the link between the signature and clinical characteristics, we further analyzed the patients stratified by different characteristics and investigated whether characteristics could differentiate patients’ survival in the TCGA cohort. The findings show that the two patient groups’ age (Fig. 2A, B), gender (Fig. 2C, D) and tumor stage (Fig. 2E, F) did not significantly differ from one another, suggesting that the signature was not influenced by these clinical factors and was independent in guiding the prognosis of patients. However, Breslow depth exhibited significant differences, reflecting the consistency between risk scores and Breslow depth (Fig. 2G, H). The pathogenic variants of the BRAF and NRAS are important for CM. So, we assessed the alterations in NRAS and BRAF within the groups with the highest and lowest risk. The analysis revealed no significantly different results between two groups, indicating that BRAF mutation (BRAF p.V600, Fig. 2I, J) and NRAS mutation (NRAS p.Q61, Fig. 2K–M) likewise did not affect the prognostic effect of the signature.

Fig. 2

Distribution of risk ratios and scores for patients with different age (A, B), gender (C, D), tumor stage (E, F), Breslow depth (G, H), BRAF (I, J) and NRAS (K–M) mutation in the TCGA cohort. **p < 0.01, *p < 0.05, ns: no significance

2.4 Analysis of independence and nomogram construction

We performed Cox regression to examine the signature’s independence from additional clinical characteristics (age, breslow depth, gender, tumor stage, and BRAF and NRAS mutations) in the TCGA cohort (Fig. 3A, B). The findings indicated that the signature could be used as an independent prognostic indicator because risk score, age, stage, and Breslow depth were all strongly related to overall survival (OS) in univariate-Cox-regression and remained so in multivariate-Cox -regression. Next, we used age, stage, breslow depth, and risk score to create nomograms for OS (Fig. 3C). Each predictor with a given value can be mapped to the Points axis, and the sum of these points can be referred to in the Total Points axis. Then the probability of 1-, 3- and 5-year survival for patient can be obtained from corresponding axis. The ROC curve analysis exhibited that the nomograms had a good predictive accuracy, and the AUC value of 1-, 3- and 5-year were 0.86, 0.84 and 0.80, respectively (Figure S3), indicating that the OS prognostic nomograms can guide the clinical management of patient.

Fig. 3

Univariate- (A) and multivariate-Cox-regression (B) of characteristics (age, gender, stage, and Breslow depth, BRAF- and NRAS-mutant, risk score) in the TCGA cohort. C The nomogram of characteristics in the TCGA cohort

2.5 Assessment regarding the CRG DNA methylation signature grouping by clinical characteristics

Different tumor locations can impact patient prognosis (Table 2). For the primary tumor site, Kaplan–Meier survival analysis demonstrates that there is no significant difference in survival between the high- and low-risk groups in the TCGA cohort (p = 0.165, Figure S4A). ROC curve analysis indicates that the predictive model holds certain prognostic value, with an AUC value of 0.725 (Figure S4B). For the remaining clinical characteristics included Distant Metastasis, Regional Cutaneous or Subcutaneous, and Regional Lymph, Kaplan–Meier survival analyses showed that patients in the high-risk group all had significantly lower survival rates than those in the low-risk group (Figure S4C-H). The signature demonstrates high sensitivity and specificity in predicting OS associated with these clinical characteristics, with AUC values exceeding 0.72 for all.

Table 2 Kaplan–Meier and ROC results of different regrouping

Different ulceration indicators also effects on the prognosis of patients (Table 2). For patients without ulceration (ulceration NO) in the TCGA cohort, Kaplan–Meier survival analysis demonstrates that there is no significant difference in survival between the high- and low-risk groups (p = 0.106, Figure S5A). ROC curve analysis reveals that the predictive model has some predictive value with an AUC of 0.703 (Figure S5B). For CM patients with ulceration (ulceration YES), Kaplan–Meier survival analysis indicated that patients in the high-risk group had significantly lower OS than those in the low-risk group (p = 0.002, Figure S5C). ROC curve analysis suggested that the AUC values are both greater than 0.7. This suggests that for melanoma patients with ulceration, the predictive accuracy is 0.730 (Figure S5D).

Some patients with CM receive different therapy in the TCGA cohort (Table 2). For these patients who did undergo radiation therapy (Radiation Therapy YES), Kaplan–Meier survival analysis showed that there was no significant difference between the high- and the low-risk group (p = 0.096, Figure S6A). ROC curve analysis indicated that this predictive model has some predictive value, with an AUC of 0.848 (Figure S6B). For patients who did not undergo radiation therapy (Radiation Therapy NO), shorter OS for high-risk patients compared to low-risk patients (p < 0.001, Figure S6C) and the AUC of the ROC curve was 0.818 (Figure S6D). In addition, Kaplan–Meier survival analysis showed significant differences between high- and low-risk groups for both patients with immunotherapy (Immunotherapy YES, p = 0.003) and non-immunotherapy (Immunotherapy NO, p < 0.001), with AUC of 0.846 and 0.727 (Figure S7). Similarly, there were significant differences between high- and low-risk groups of chemotherapy patients (Chemotherapy YES, p = 0.003) and non-chemotherapy patients (chemotherapy NO, p < 0.001), with AUC of 0.822 and 0.723, respectively (Figure S8). These results suggest that this signature still has high accuracy in predicting OS after regrouping by the above clinical characteristics.

2.6 Association of CRG DNA methylation prognostic signature with immune landscape

To further understand the potential relationship between the signature and the immune landscape of CM patients, we analyzed the correlation of CRG DNA methylation signature with immune and compared differences in immune cell infiltration among the high- and low-risk groups in the TCGA cohort. The results indicated a negative correlation between risk scores and ESTIMATE scores (tumor purity), immune scores (the infiltration level of immune cells in tumor tissues), and stromal scores (the level of stromal cells present in tumor tissue) (Fig. 4). This result suggests that the lower the three scores, the higher the patient risk. The risk scores were negatively correlated with Plasma cells, T cells CD8, T cells CD4 memory activated, T cells gamma delta, B cells memory and macrophages M1, while positively correlated with Neutrophils, Monocytes, NK cells resting, Macrophages M0, Macrophages M2, and T cells CD4 memory resting (Fig. 4). In addition, the immune cell infiltration level of B cells memory, Macrophages M0, Macrophages M1, Macrophages M2, Monocytes, NK cells resting, Plasma cells, T cells CD4 memory activated, T cells CD4 memory resting, T cells CD8,T cells gamma delta have significant differences between high- and low-risk groups (Fig. 5A). Meanwhile, the analysis confirmed differential ESTIMATE scores, immune scores, and stromal scores between the high- and low-risk groups, with all three scores significantly increased in the low-risk group (Fig. 5B). To investigate tissue-specific and the impact of immunotherapy, we explored immune infiltration with primary tumor, metastatic, immunotherapy, and non-immunotherapy patients. The ESTIMATE scores, immune scores, the immune cell infiltration level of T cells CD4 memory activated, T cells CD4 memory resting, T cells CD8, Macrophages M0, Macrophages M2, and Neutrophils have significant differences between high- and low-risk groups in primary tumor patients (Figure S9). The ESTIMATE scores, immune scores, stromal scores, the immune cell infiltration level of B cells naive, T cells CD4 memory activated, T cells CD4 memory resting, T cells CD8, Plasma cells, T cells gamma delta, Macrophages M1, Macrophages M2, Monocytes, NK cells resting, and Plasma cells have significant differences between two risk groups in metastatic patients (Figure S10). On the non-immunotherapy patients, the ESTIMATE scores, immune scores, stromal scores, the immune cell infiltration level of B cells memory, Macrophages M0, Macrophages M1, Macrophages M2, Monocytes, NK cells resting, Plasma cells, T cells CD4 memory activated, T cells CD4 memory resting, T cells CD8, and T cells gamma delta have significant differences between different groups (Figure S11). On patients with immunotherapy, only immune scores, the immune cell infiltration level of T cells CD4 memory activated and T cells gamma delta have significant differences (Figure S12). Considering the significance of immune checkpoint-associated genes (ICGs) in clinical treatment, we further compared the expression differences of ICGs between the high- and low-risk groups. The ICGs included CD44, CD276, CD4, CD8A, IDO1, CD86, PD-1, CTLA4, PD-L2, PD-L1, and CD80. The results indicated that there was a direct correlation between the risk scores and ICGs expression (Fig. 4). ICGs expression was significantly different between the high- and low-risk group (Fig. 5C).

Fig. 4

Correlation of CRG DNA methylation signature with 22 types of immune cell infiltration, stromal scores (the level of stromal cells present in tumor tissue), immune scores (the infiltration level of immune cells in tumor tissues), ESTIMATE scores (tumor purity), and ICGs in the TCGA cohort. A Correlation p value of CRG DNA methylation characteristics with 22 kinds of immune cell infiltration, immune-related scores and ICGs. The size of the dot represents the p value. The larger the dot, the smaller the p value. The color of the dot represents the category, blue represents immune cell infiltration, green represents immune-related scores, and red represents ICGs. B Correlation coefficient of CRG DNA methylation characteristics with 22 kinds of immune cell infiltration, immune-related scores and ICGs

Fig. 5

The differences in 22 types of immune infiltration (A), stromal scores (the level of stromal cells present in tumor tissue), immune scores (the infiltration level of immune cells in tumor tissues), ESTIMATE scores (tumor purity) (B), and expression of ICGs (C) between high- and low-risk CM patients in the TCGA cohort

2.7 Functional enrichment analyses

To investigate the potential mechanisms underlying different prognosis risk groups, we identified differentially expressed genes (DEGs) between the high- and low-risk groups in the TCGA cohort (Table S1). Thereafter, Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analyses were performed on these DEGs. KEGG enrichment analysis revealed significant enrichment in the cytokine-cytokine receptor interaction, viral protein interaction with cytokine and cytokine receptor, hematopoietic cell lineage, neuroactive ligand-receptor interaction, etc. (Fig. 6A, Table S2). GO functional enrichment analysis results indicated that the DEGs were mainly enriched in T-cell receptor complex, plasma membrane signaling receptor complex, and production of molecular mediator of immune response (Fig. 6B, Table S3). The above results suggest that the risk signature potentially contributes to the formation of immune microenvironment in CM patients.

Fig. 6

The KEGG enrichment analysis (A) and GO enrichment analysis (B) of DEGs in high- and low-risk groups in the TCGA cohort

Furthermore, we performed GSVA analyses to explore the differences in enriched signaling pathways between the high- and low-risk groups (Figure S13). The results showed that significant differences in natural killer cell mediated cytotoxicity, toll like receptor signaling pathway, T cell receptor signaling pathway, JAK stat signaling pathway and oxidative phosphorylation between the two risk groups. At the same time, there are differences in GO processes such as interferon gamma mediated signaling pathway, regulation of NK T cell activation, regulation of response to tumor cell, immune response to tumor cell, regulation of response to tumor cell, immune response to tumor cell, positive regulation of natural killer cell activation.

2.8 Association of CRG DNA methylation prognostic signature with gene expression

To examine the tissue characteristics of DNA methylation sites, we analyzed the differences in methylation sites between patients with metastatic and primary tissues in the TCGA cohort. Our analysis revealed the methylation levels of cg01130192 and cg18645035 were significantly higher in the metastatic tissues. However, the regression coefficient of cg01130192 and cg18645035 in the CRG DNA methylation prognostic signature is negative, indicating that they are negatively correlated with patients’ survival risk. The higher the methylation level of sites, the lower the survival risk. These results suggest that the metastasis and primary tumor are not necessarily consistent with the survival risk of patients. In contrast, the methylation levels of cg00668852, cg22134611, and cg25522181 did not differ between the two tissue types (Figure S14).

To investigate the relationship between DNA methylation sites and gene expression, we conducted Spearman correlation analyses between DNA methylation sites and their corresponding gene expressions. Our findings showed that the methylation levels of cg00668852 and cg18645035 were significantly positively correlated with the expression of their respective genes. In addition, we observed a significant negative correlation between the methylation levels of cg22134611, cg25522181, and cg01130192 and their corresponding gene expressions (Figure S15).