Identification and Validation of Ferroptosis-Related DNA Methylation Signature for Predicting the Prognosis and Guiding the Treatment in Cutaneous Melanoma

Cutaneous melanoma (CM) is one of the most aggressive skin tumors with a poor prognosis. Ferroptosis is a newly discovered form of regulated cell death that is closely associated with cancer development and immunotherapy. The aim of this study was to establish and validate a ferroptosis-related gene (FRG) DNA methylation signature to predict the prognosis of CM patients using data from The Cancer Genome Atlas (TCGA) and the Gene Expression Omnibus (GEO) database. A reliable four-FRG DNA methylation prognostic signature was constructed via Cox regression analysis based on TCGA database. Kaplan–Meier analysis showed that patients in the high-risk group tended to have a shorter overall survival (OS) than the low-risk group in both training TCGA and validation GEO cohorts. Time-dependent receiver operating characteristic (ROC) analysis showed the areas under the curve (AUC) at 1, 3, and 5 years were 0.738, 0.730, and 0.770 in TCGA cohort and 0.773, 0.775, and 0.905 in the validation cohort, respectively. Univariate and multivariate Cox regression analyses indicated that the signature was an independent prognostic indicator of OS in patients with CM. Immunogenomic profiling showed the low-risk group of patients had a higher immunophenoscore, and most immune checkpoints were negatively associated with the risk signature. Functional enrichment analysis revealed that immune response and immune-related pathways were enriched in the low-risk group. In conclusion, we established and validated a four-FRG DNA methylation signature that independently predicts prognosis in CM patients. This signature was strongly correlated with the immune landscape, and may serve as a biomarker to guide clinicians in making more precise and personalized treatment decisions for CM patients.


Introduction
Cutaneous melanoma (CM) is one of the most aggressive malignant skin tumors that arise from mutations in melanocytes [1]. The incidence of CM has increased rapidly over the past decades [2,3]. Although medical advances have been able to alleviate symptoms and reduce the burden of tumors, they still provide limited help in improving survival, with 5-year overall survival (OS) of 15% [4,5]. An accurate assessment of CM progression is essential to improve OS [6]. Therefore, there is a need to identify novel robust biomarkers to predict patient prognosis.
Excessive ultraviolet (UV) radiation is an important environmental trigger in the pathogenesis of CM [7]. UV from sunlight can lead to excessive intracellular production of reactive oxygen species (ROS), causing oxidative stress damage to cells [8]. Mitochondria, as a major iron-rich and ROS-producing organelle, is considered to be an important site of cell ferroptosis [9], and the accumulation of ROS is one of the characteristics of ferroptosis [10]. Thus, the role of ferroptosis-related genes (FRGs) in the development and prognosis of CM has attracted our interest. Ferroptosis is a newly discovered form of programmed cell

Screening and Establishment of Prognosis-Related FRG DNA Methylation Signature
A total of 567 FRGs and corresponding 8781 DNA methylation sites of these FRGs were extracted. Through univariate Cox regression analysis, 631 DNA methylation sites that were significantly associated with patient prognosis were screened out (p < 0.01). These DNA methylation sites were subsequently used in multivariate Cox regression analyses to construct a prognostic signature. An optimum model consisting of four DNA methylation sites (cg12336709, cg23750391, cg15674193, and cg06904403) was constructed for predicting OS. Forest plots showed the association of these four DNA methylation sites with patient OS in a univariate and multivariate Cox proportional hazard model ( Figure 1). The genes corresponding to these four sites were ZEB1 (zinc finger E-box binding homeobox 1), CISD1 (CDGSH iron sulfur domain 1), LRRFIP1 (LRR binding FLII interacting protein 1), and DUOX1 (dual oxidase 1). The chromosomal locations of these four DNA methylation sites and their p values in Cox regression analysis are shown in Table 1.

Evaluation of the FRG DNA Methylation Prognosis Signature
The risk scores of CM patients were calculated according to the riskScore formula: riskScore = (−1.959 × β value of cg12336709) + (−9.761 × β value of cg2375039) + (−1.229 × β value of cg15674193) + (1.418 × β value of cg06904403). Patients were separated into high-and low-risk groups based on the median risk score (−2.137) (Figure 2A), and the number of deaths increased with the increasing risk score ( Figure 2B). Principal component analysis (PCA) showed that the contribution of the four components was sustained and there were identifiable dimensions between the high-and low-risk groups ( Figure 2C). In

Evaluation of the FRG DNA Methylation Prognosis Signature
The risk scores of CM patients were calculated according to the riskScore formula: riskScore = (−1.959 × β value of cg12336709) + (−9.761 × β value of cg2375039) + (−1.229 × β value of cg15674193) + (1.418 × β value of cg06904403). Patients were separated into highand low-risk groups based on the median risk score (−2.137) (Figure 2A), and the number  Kaplan-Meier survival analysis showed significantly lower survival probability in the high-risk group (p = 1.19 × 10 −12 ) ( Figure 3A). ROC curve analysis indicated good predictive accuracy for the four-DNA methylation signature, with area under the curve (AUC) values of 0.738 (95% CI: 0.66-0.82), 0.730 (95% CI: 0.68-0.79), and 0.770 (95% CI: 0.72-0.82) for the 1-year, 3-year, and 5-year survival, respectively ( Figure 3B). To further investigate the prognostic value of this signature, we analyzed another independent DNA methylation dataset of CM patients. Patients were divided into high-and low-risk groups according to the median risk score (−2.137) obtained from the training TCGA cohort. As expected, survival prognosis was significantly poorer in patients with high-risk scores (p = 0.0052) ( Figure 3C), with AUCs of 0.773 at 1 year, 0.775 at 3 years, and 0.905 at 5 years ( Figure 3D), confirming the good prognostic predictive role of the signature for prognosis assessment of CM patients.   methylation dataset of CM patients. Patients were divided into high-and low-risk groups according to the median risk score (−2.137) obtained from the training TCGA cohort. As expected, survival prognosis was significantly poorer in patients with high-risk scores (p = 0.0052) ( Figure 3C), with AUCs of 0.773 at 1 year, 0.775 at 3 years, and 0.905 at 5 years ( Figure 3D), confirming the good prognostic predictive role of the signature for prognosis assessment of CM patients.

Correlation between FRG DNA Methylation Signature and Clinical Characteristics
We further analyzed the value of the signature in patients stratified by different clinical factors in TCGA, and investigated whether clinical characteristics could distinguish patients' prognostic survival. The results demonstrated that male patients had higher risk scores and poorer OS ( Figure 4A), whereas there was no significant difference between patients stratified by gender in survival analysis ( Figure 4E). Patients who were older (≥60) or with thicker Breslow depth (≥2 mm) had significantly higher risk scores and worse OS ( Figure 4B,C,F,G). Stage II patients tended to have higher risk scores than stage I patients ( Figure 4D), and stratified survival analysis indicated that patients with advanced stage (stages III and IV) had poorer OS ( Figure 4J). Pathogenic variation of BRAF and NRAS genes plays a very important role in CM, about 50% of CM carry an activating BRAF mutation [26]. Here, we analyzed the alterations of NRAS and BRAF in high-and low-groups, and found that patients in the low-risk group had more BRAF mutations and fewer NRAS mutations. BRAF-mutant patients had lower risk scores than BRAF wild-type (BRAF-wt) patients ( Figure 4F), but the survival analysis showed no significant difference between BRAF-mutant and BRAF-wt groups ( Figure 4L). For the NRAS gene, there were no significant differences in risk scores or survival analyses between the NRAS-mutant and NRAS-wt groups ( Figure   To explore the independence of the signature from other clinical factors (gender, age, stage, Breslow depth, and NRAS and BRAF variations), we performed univariate and multivariate Cox regression analyses ( Figure 5A,B). The results showed that risk score, age, and stage were correlated with the OS of CM patients in univariate Cox regression, and remained significant in the multivariate Cox regression analysis, suggesting that the signature could be regarded as an independent prognostic indicator. Subsequently, we developed a nomogram for OS prediction using the risk score, age, and stage ( Figure 5C). The results indicated that the predictive nomogram for OS was useful for clinical management.

Association of FRG DNA Methylation Prognostic Signature with Immune Landscape
To further understand the potential correlation of the signature with the immune landscape of the CM samples, we compared the differences in the various immune cell components between the low-and high-risk groups. The results indicated that risk scores were negatively correlated with ESTIMATE scores, immune scores, and stromal scores ( Figure 6A). The ESTIMATE score, immune score, and interstitial score were significantly increased in the low-risk group ( Figure 6B-D). TME immune cell infiltration analysis showed that high-risk score was negatively associated with plasma cells, CD8+ T cells, activated CD4 memory T cells, Gamma Delta T cells, and M1 macrophages, while it was significantly positively correlated with follicular helper T cells, resting NK cell, M2 macrophages, and activated mast cell ( Figure 6A,E).
Considering the importance of checkpoint inhibitors in clinical treatment, we further compared the differences in the expression of 12 ICBs in the high-and low-risk groups and found substantial differences in PD-1, PD-L1, PD-L2, CTLA4, CD276, CD80, IDO1, CD4, CD8A, CD8B, and CD86 between the two groups ( Figure 7). These results indicated that the low-risk group exhibited a higher immune infiltration status and that immune checkpoint-related pathways might play an essential role in the good prognosis of the low-risk group. This signature could help to select the appropriate checkpoint inhibitors to treat CM patients. Finally, considering that immune infiltration is an important parameter for CM patients, a nomogram was constructed to predict patients' OS based on risk score, age, stage, and immune score ( Figure 8).

Functional Enrichment Analyses
To investigate the underlying mechanism in different prognostic risk groups, we identified the DEGs between the high-and low-risk groups. A total of 716 up-regulated and 76 down-regulated genes were identified in the low-risk group ( Figure 9A). Thereafter, GO, KEGG, and Reactome analyses were performed on these DEGs. GO enrichment analysis demonstrated that the up-regulated genes in the low-risk group were mainly associated with immune response, immunological synapse, T cell activation, and positive regulation of T cell proliferation ( Figure 9B and Table S1). KEGG pathway enrichment analysis revealed that immune-related pathways such as cytokine-cytokine receptor interaction, chemokine signaling pathway, NF-kappa B signaling pathway, and B cell receptor signaling pathway were significant ( Figure 9C and Table S2). Reactome analysis showed that these up-regulated genes were mainly enriched in the immune system, cytokine signaling in the immune system, the adaptive immune system, and FCGR activation ( Figure 9D and Table S3). The above results illustrated that the risk signature might be involved in immune microenvironment formation in CM patients.    CD4, CD8A, CD8B, and CD86 between the two groups ( Figure 7). These results indicat that the low-risk group exhibited a higher immune infiltration status and that immu checkpoint-related pathways might play an essential role in the good prognosis of t low-risk group. This signature could help to select the appropriate checkpoint inhibito to treat CM patients. Finally, considering that immune infiltration is an important param eter for CM patients, a nomogram was constructed to predict patients' OS based on ri score, age, stage, and immune score (Figure 8).

Functional Enrichment Analyses
To investigate the underlying mechanism in different prognostic risk groups, we identified the DEGs between the high-and low-risk groups. A total of 716 up-regulated and 76 down-regulated genes were identified in the low-risk group ( Figure 9A). Thereaf-

Discussion
CM is a significant cause of skin cancer-related death worldwide with a poor prognosis [3,27]. Despite advances in CM treatment with targeted therapies combined with immunotherapy in recent decades, cross-drug resistance still affects the prognosis of patients. Well-established risk assessment and treatment stratification tools can help to select treatment options and improve patients' outcomes [28]. In this study, we combined ferroptosis-related gene (FRG) sets and other datasets to construct and validate a ferroptosisrelated DNA methylation signature in CM. The signature may be the first risk-prediction signature based on FRGs DNA methylation with good clinical applicability in CM.
Ferroptosis refers to an iron-dependent, non-apoptotic form of cell death that plays an essential role in the tumor microenvironment (TME) [11,29]. Dysregulation of ferroptosis has been associated with the development and prognosis of several cancers [12,30]. Among the genes corresponding to DNA methylation in the prognostic signature, we identified that ZEB1 and DUOX1 are drivers of ferroptosis and that CISD1 is a suppressor of ferroptosis, which play a crucial role in cancer development [31,32]. For example, ZEB1 is a major element in the control of epithelial-to-mesenchymal transition (EMT) and is

Discussion
CM is a significant cause of skin cancer-related death worldwide with a poor prognosis [3,27]. Despite advances in CM treatment with targeted therapies combined with immunotherapy in recent decades, cross-drug resistance still affects the prognosis of patients. Well-established risk assessment and treatment stratification tools can help to select treatment options and improve patients' outcomes [28]. In this study, we combined ferroptosis-related gene (FRG) sets and other datasets to construct and validate a ferroptosisrelated DNA methylation signature in CM. The signature may be the first risk-prediction signature based on FRGs DNA methylation with good clinical applicability in CM.
Ferroptosis refers to an iron-dependent, non-apoptotic form of cell death that plays an essential role in the tumor microenvironment (TME) [11,29]. Dysregulation of ferroptosis has been associated with the development and prognosis of several cancers [12,30]. Among the genes corresponding to DNA methylation in the prognostic signature, we identified that ZEB1 and DUOX1 are drivers of ferroptosis and that CISD1 is a suppressor of ferroptosis, which play a crucial role in cancer development [31,32]. For example, ZEB1 is a major element in the control of epithelial-to-mesenchymal transition (EMT) and is closely related to ferroptosis sensitivity in cancer cells, which can influence tumorigenesis from the early steps of cancer [32][33][34]. Moreover, ZEB1 promotes immune escape in melanoma [35]. CISD1, a negative regulator of autophagy, was differentially expressed in a variety of tumors and correlated with OS [36,37]. Meanwhile, CISD1 can be served as a prognostic biomarker in patients with hepatocellular carcinoma [38]. DUOX1 is a new target for macrophage reprogramming, and is involved in the redox regulation of EGFR signaling; ATP-mediated DUOX1 activation led to the activation of the NF-kappaB pathway and production of IL-8 [39][40][41]. LRRFIP1 plays an important role in the invasion of tumor cells [42]. High expression of LRRFIP1 was found to be associated with a better response to teniposide in glioblastoma, and could be a candidate gene for tumor-targeted therapy [43]. The major biological characteristics of these four genes are summarized in Table 2. Here, we found the signature consisting of DNA methylation of these four genes was closely correlated with the prognosis of CM patients. Aberrant DNA methylation is associated with tumorigenesis of diverse malignancies, including CM [44,45]. Epigenetic changes have been shown to alter gene expression and play a crucial role in the occurrence and progression of cancer [46]. For example, PD-L1 and PD-L2 promoter methylation can regulate their mRNA expression and are associated with patient OS in CM [22]. Here, correlation analysis between DNA methylation levels and gene expression indicated that DNA methylation levels at the cg12336709 site were significantly positively correlated with ZEB1 expression (p < 0.0001), while the expression of CISD1 and DUOX1 was significantly negatively correlated with their methylation levels ( Figure S1), suggesting that these DNA methylations may regulate gene expression and, thus, affect CM development.
Tumor immunotherapy is a new therapeutic approach that has made substantial progress in both fundamental research and clinical practice in recent decades [47]. Here, we calculated immune and stromal scores of tumor samples to investigate the differences in immune microenvironment between low-and high-risk groups, and found that ESTIMATE scores, immune scores, and stromal scores were significantly higher in the low-risk group. We also evaluated the relationship between the infiltration of immune cell factors with the risk score, and found that CD8 T immune cells and activated CD4 memory T cells were negatively correlated with the risk score, while M2 macrophages and activated mast cells were positively correlated with the risk score. High fractions of activated memory CD4 T cells were known to have an anti-tumor function [48], while high fractions of activated mast cells and M2 macrophage have a pro-tumorigenic effect, and are associated with poor prognosis [49,50]. Meanwhile, most immune checkpoint genes also presented strong activity in the low-risk group, such as PD-1, PD-L1, PD-L2, CTLA4, CD4, CD8A, CD8B, and IDO1. Moreover, functional enrichment analysis demonstrated that the up-regulated genes in the low-risk group were mainly associated with immune response, immunological synapse, T cell activation, positive regulation of T cell proliferation, cytokine-cytokine receptor interaction, the chemokine signaling pathway, the NF-kappa B signaling pathway, and the B cell receptor signaling pathway, the immune system, cytokine signaling in the immune system, the adaptive immune system, and FCGR activation. The up-regulated genes in the low-risk group may be able to activate immune responses. Patients with higher immune activity have a better prognosis [51]. The inhibition of NF-kappa B was correlated with increased disease-free survival in patients with breast cancer [52]. These findings indicated that the FRG DNA methylation signature might influence the prognosis of CM patients through modulating the immune microenvironment of CM.
In addition, univariate Cox regression, Kaplan-Meier, and ROC analyses were carried out on the four individual methylation sites to further analyze their prognostic significance in CM. Kaplan-Meier analysis showed that the single DNA methylation site also could distinguish patients between high-and low-risk, and patients with higher methylation levels for cg12336709, cg23750391, and cg15674193 exhibited higher OS ( Figure S2A-C), and patients with higher cg06904403 methylation had lower OS ( Figure S2D). ROC analysis showed that the predictive performances of individual methylation sites were lower than the combination of these four DNA methylation sites ( Figure S2E,F). These results implied a combination of methylation sites can provide a better prognostic prediction for CM patients.

Dataset Source and Pre-Processing
Publicly available DNA methylation data and clinical information of CM patients were downloaded from TCGA "https://portal.gdc.cancer.gov/ (accessed on 26 September 2021)" and GEO databases "https://ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE51547 (accessed on 26 September 2021)". The somatic mutation data were obtained from TCGA database, and NRAS and BRAF mutation information was extracted. DNA methylation levels were expressed as β values ranging from 0 (no methylation) to 1 (100% methylation). FRGs were acquired from the FerrDb database "http://www.zhounan.org/ferrdb/ (accessed on 17 September 2021)" [53], and the corresponding DNA methylation site information was extracted using Perl script. Patients from TCGA database were utilized as the training cohort for the identification and construction of the prognostic signature, and other data obtained from the GEO database were used as an independent validation set validation of prognostic signature.

Prognosis-Related DNA Methylation Filtering and Risk Model Construction
The FRG DNA methylation levels in TCGA data were subjected to univariate Cox regression analysis to determine DNA methylation related to OS. Statistically significant (p < 0.01) variables were then selected for multivariate Cox regression analysis to construct models comprising all possible combinations of two to five factors, aiming at further filtering out combined biomarkers associated with OS. Hazard ratios (HRs) and 95% confidence intervals (CIs) were calculated. The prognostic risk scores were calculated for each patient by summing the DNA methylation levels and the corresponding regression coefficients. Patients were divided into high-and low-risk groups based on the median risk score.

DNA Methylation Prognostic Signature Validation and Nomogram Construction
To assess the usability of the prognostic signature, the Kaplan-Meier survival curve and Wilcoxon rank test were used to demonstrate and compare the difference in prognosis between the two groups [54]. ROC curves were applied to calculate the sensitivity and specificity of the signature and to estimate the prognostic performance [55]. Univariate and multivariate Cox regression analyses were performed to assess the independence of the signature. Nomograms for the 1-, 3-, and 5-year OS were constructed to predict survival rates by summing the points corresponding to each factor of patients.

Immunogenomic Landscape Evaluation
The infiltration levels of stromal and immune cells were calculated for each CM patient using the "ESTIMATE" R package [56]. Three kinds of immunophenoscores and infiltration of 22 immune cells for CM patients were obtained, and the differences between the high-and low-risk groups were compared. We also analyzed the correlation of FRG DNA methylation signature with immune checkpoint-related genes (ICGs), and compared the expression values of ICGs in patients between high-and low-risk groups. Correlations were calculated using Spearman correlation analysis, and the differences between the two groups were evaluated using Mann-Whitney U-test.

Functional Enrichment Analysis
Differentially expressed genes (DEGs) between low-and high-risk patients were determined using R "limma" package, and the adjusted p < 0.01, |logFC| ≥ 0.5 were considered a statistically significant difference. GO functional enrichment and KEGG pathway analysis were conducted and the biological function and signaling pathways were demonstrated by ggplot2.

Survival Analysis
Statistical analyses were performed using R language (version 4.0.3). Univariate and multivariate Cox proportional hazards analyses were performed using the R "survival" package. Time-dependent ROC analysis was conducted using R "survivalROC" package, and ROC curves were plotted for 1-, 3-, and 5-year. Nomograms were established with the R "rms" package. The differences between survival curves were detected using the generalized Wilcoxon test. In addition, the correlation between DNA methylation sites and the relative expression of the corresponding genes were examined by Spearman correlation analysis. Unless otherwise stated, p < 0.05 was considered significant.

Conclusions
In summary, we constructed and validated a ferroptosis-related four-DNA methylation signature that could be utilized as an independent and effective biomarker to predict the prognosis of CM patients. The signature was correlated with tumor immune infiltration and immune checkpoint genes, and would facilitate the selection of more effective immunotherapeutic strategies. Functional enrichment analysis revealed that the signature might be involved in immune microenvironment formation in CM patients. Despite the specific function requiring further exploration, this study provided a theoretical foundation for improving the clinical treatment and promoting the development of personalized immunotherapy and precision medicine for CM.

Supplementary Materials:
The following supporting information can be downloaded at: https://www.mdpi.com/article/10.3390/ijms232415677/s1, Figure S1: Relationship between methylation level and its gene expression level of each site was evaluated using Spearman correlation analysis. Gene expression levels are reported as log 2 (FPKM), and methylation β-values are defined by the Infinium HumanMethylation 450 BeadChip. The Reported p values are bipartite. Figure S2: Association of single DNA methylation sites with patient OS. (A-D) Kaplan-Meier analysis with Wilcoxon test was performed to estimate the differences in OS between patients with different DNA methylation levels. (E,F) The AUC value of individual methylation level in predicting 3-year and 5-year OS. Table S1: The top 50 Gene Ontology (GO) terms of differentially up-regulated genes in low-risk group. Table S2: The top 50 KEGG pathways of differentially up-regulated genes in low-risk group. Table S3: The top 50 Reactome pathways of differentially up-regulated genes in low-risk group.