Identification and validation of a ferroptosis‐related gene signature for predicting survival in skin cutaneous melanoma

Abstract Purpose Ferroptosis plays a crucial role in the initiation and progression of melanoma. This study developed a robust signature with ferroptosis‐related genes (FRGs) and assessed the ability of this signature to predict OS in patients with skin cutaneous melanoma (SKCM). Methods RNA‐sequencing data and clinical information of melanoma patients were extracted from TCGA, GEO, and GTEx. Univariate, multivariate, and LASSO regression analyses were conducted to identify the gene signature. A 10 FRG signature was an independent and strong predictor of survival. The predictive performance was assessed using ROC curve. The functions of this gene signature were assessed by GO and KEGG analysis. The statuses of low‐risk and high‐risk groups according to the gene signature were compared by GSEA. In addition, we investigated the possible relationship of FRGs with immunotherapy efficacy. Results A prognostic signature with 10 FRGs (CYBB, IFNG, FBXW7, ARNTL, PROM2, GPX2, JDP2, SLC7A5, TUBE1, and HAMP) was identified by Cox regression analysis. This signature had a higher prediction efficiency than clinicopathological features (AUC = 0.70). The enrichment analyses of DEGs indicated that ferroptosis‐related immune pathways were largely enriched. Furthermore, GSEA showed that ferroptosis was associated with immunosuppression in the high‐risk group. Finally, immune checkpoints such as PDCD‐1 (PD‐1), CTLA4, CD274 (PD‐L1), and LAG3 were also differential expression in two risk groups. Conclusions The 10 FRGs signature were a strong predictor of OS in SKCM and could be used to predict therapeutic targets for melanoma.


| INTRODUCTION
Malignant melanoma, which is increasing in incidence worldwide, is the most aggressive form of skin cancer. There were approximately 106,110 new cases of melanoma and 7180 deaths from melanoma in the United States in 2021. 1 Over the last two decades, the range of available systemic treatment options for melanoma has expanded and include immunotherapy (CTLA-4, PD-1, LAG-3, and PD-L1 inhibitors), MEK, and targeted therapy (BRAF inhibitors). [2][3][4][5][6] However, the prognosis of advanced or metastatic melanoma is poor, 7 and the demand for novel biomarkers to improve melanoma early diagnosis and prognosis has increased in the recent years. 8 Moreover, despite improvements in staging, the prognosis of advanced melanoma is heterogeneous, with high variability in the overall estimate among stages III and IV. 9,10 Thus, the development of novel signatures that can significantly improve melanoma diagnosis and increase the accuracy of predicting melanom prognosis is needed.
Ferroptosis is a form of programmed cell death caused by iron-dependent lipid peroxidation 11 and has a crucial role in inhibiting tumor cell proliferation, invasion, and metastasis. 12 Ferroptosis can trigger immune responses, especially in malignant tumors resistant to conventional therapies. [13][14][15] Ferroptosis has a dual role in cancer because ferroptotic tumor cells release some signaling molecules that either promote or inhibit the tumor proliferation [16][17][18] ; however, the role of these molecules in tumor is incompletely understood. 16 A current research has demonstrated that lymph protects melanoma from ferroptosis, increases melanoma cell viability during metastasis, and increases the formation of distant metastasis, providing new perspectives for studying metastasis. 19 Therefore, identifying ferroptosis-related biomarkers for the prognosis of skin cutaneous melanoma (SKCM) is essential.
Some studies identified prognostic signatures based on FRGs in tumors based on public databases. 20,21 In this respect, a 19-gene signature predicted glioma cell death and prognosis, 20 and a 10-gene signature predicted OS in hepatocellular carcinoma. 21 However, no studies have identified ferroptosis-related signature that can predict OS in patients with SKCM. To identify and validate a prognostic signature and improve the diagnosis and therapy of SKCM, TCGA, GTEx, and GEO databases were comprehensively searched for novel and previously identified FRGs. The functions of this signature in the tumor microenvironment were assessed through enrichment analysis.

| Data acquisition and processing
The RNA-seq expression profiles of SKCM and normal skin tissues were obtained from TCGA and GTEx databases, respectively. 22 The GSE65094 gene expression dataset and follow-up clinical information were downloaded from the GEO. Read counts were normalized through scale. The TCGA-SKCM dataset (471 samples) and GSE65094 dataset (214 samples) were selected as the training and validation sets, respectively. Then, 253 FRGs were obtained from the FerrDb database. 23

| Analysis of DEGs
To identify significant DEGs in melanoma, RNA-seq data extracted from TCGA (471 SKCM samples and one normal skin sample) and the GTEx database (812 normal skin samples) were normalized by log 2 -transformed FPKM and merged into one dataset for subsequent analysis. DEGs were determined by "limma" package in R.

| Signature construction and validation
The correlation between FRGs and OS in TCGA training set was analyzed by univariate Cox regression analysis using the "survival" package. FRGs with p < 0.01 were considered to have prognostic significance. Differentially expressed FRGs were identified by overlapping prognostic FRGs and DEGs. A gene signature with 10 FRGs was selected by LASSO regression. 24,25 LASSO estimates were based on 1000-fold cross-validation to avoid overfitting. In this model, risk scores were calculated via multiplying the expression values of 10 FRGs by their regression coefficients (A), as follows: risk score = (expression of gene 1 × A 1 of gene 1) + (expression of gene 2 × a A 2 of gene 2) + (expression of gene N × A N of gene N).
Patients were divided into a low-and a high-risk group based upon median risk score. Survival status were displayed in the Kaplan-Meier curve. The model predictive ability was assessed using heat maps, forest plots, risk score maps, OS curves, and ROC curves. Based on prognostic ferroptosisrelated DEGs, PCA, and t-SNE were generated using the "stats" and "Rtsne" package of R. External validation of the signature was performed using the GSE65904 dataset. This dataset was divided into a high-and low-risk group based upon the cutoff value of the signature. OS curves were compared between these two groups using Kaplan-Meier analysis to validate the signature.

| Gene set enrichment analysis
GO and KEGG analyses of FRGs differentially expressed were performed using "clusterProfiler" package in R. GSEA of GO and KEGG gene sets was performed using the GSEA software version 4.1.0. The number of random sample permutations in each analysis was set to 1000.

| Immune correlation analysis
The ssGSEA was performed to assess the tumorinfiltrating immune cell subsets and immune-related functions. 21 The expression of immune checkpoint gene might predict therapeutic effect of immune checkpoint inhibitors (ICIs). 26 Thus, we investigated four ICIs: PD-1 and PD-L1, CTLA-4, and LAG-3 in melanoma. 27 We also estimated the relationship between ICIs and risk score through the spearman correlation analysis, which intended to assess potential application of FRGs signature in immunotherapy.

| Validation
A study of single cell malignant melanoma transcriptomes defined two main transcriptional states of melanoma cells: the MITF and AXL gene programs. 28 We choose the A2058 (MITF) and A375 (AXL) cell lines for our study. Co.,Ltd. These three cell lines were cultured in DMEM (high-glucose) medium (Gibco) containing 10% FBS (Gibco) at 37°C with 5% CO 2 in an incubator. Total RNA of HaCaT, A2058, and A375 cells were extracted using RNA-easy kit (Vazyme). The primers used in this study were synthesized by GENECREATE (WUHAN GENECREATE BIOLOGICAL ENGINEERING, LTD) ( Table 1). Following, the reverse transcription was conducted with the HiFiScript cDNA synthesis kit (Vazyme) to generate cDNA. The qPCR was performed using an LineGene 9600 Plus instrument (Bioer Technology) and 2× SYBR Green Qpcr MasterMix (SEVEN BIOTECH).
The CT values were normalized to the expression of the endogenous housekeeping gene GAPDH, and the 2(−ΔΔCt) values were calculated for relative quantification. The reactions were performed in triplicate. The comparisons among multiple groups were conducted by one-way ANOVA. Statistical analyses were carried out using GraphPad Prism 9.0.0 software.

| Identification of DEGs in SKCM
The total 104 FRGs were identified to match the RNA-seq data from the TCGA. Using "limma" package with an absolute log 2 -fold change (FC) >1 and an adjusted p < 0.05,

| Screening of differentially expressed FRGs
The comparison of RNA-seq data of melanoma samples (n = 471) and normal skin samples (n = 813) in the TCGA and GTEx databases identified 44 DEGs with an absolute log 2 -FC > 1 and an adjusted p < 0.05 (FC ≥ 1, FDR ≤ 0.05) (Figure 1). The total 253 FRGs were obtained from FerrDb database. The univariate Cox proportional hazard regression analysis showed that 55 FRGs were significant predictors of OS (p <0.01) ( Figure 2A). Of these, 13 FRGs were differentially expressed (four upregulated and nine downregulated) in tumor tissues ( Figures 2B-D). The relationship in these genes is displayed in Figure 2E.

| Construction of a prognostic model using the TCGA dataset
A prognostic signature with 10 FRGs (CYBB, IFNG The patients were divided into a high-risk and a lowrisk group according to the median cutoff value. The Kaplan-Meier analysis indicated that prognosis was significantly worse in the high-risk group than in the low-risk group in the training dataset ( Figure 3A, p < 0.01). The risk score distribution plot and survival status curves indicated that the high-risk group from this dataset had lower survival ( Figure 3B) and higher risk scores ( Figure 3C). Moreover, as the risk score increased, the expression of protective genes (CYBB, FBXW7, IFNG, ARNTL, JDP2, TUBE1, and HAMP) decreased, whereas the expression of risk genes (PROM2, GPX2, and SLC7A5) increased in this dataset ( Figure 3D). Similar results were obtained in the validation dataset ( Figures 3E-H). The AUC value of the risk score (0.754) was higher than those of clinical indicators, such as gender, age, and TNM stage ( Figure 4A). PCA and t-SNE analyses revealed that samples were distributed in two principal components (Figures 4B,C). The FRGs genes were differentially expressed between different T stages ( Figure 4D, p < 0.001).

| Functional enrichment analysis
GO and KEGG enrichment analyses of risk genes differentially expressed between the high-risk and lowrisk groups were performed to assess the biological functions and pathways correlated with risk scores.
GO analysis indicated that the DEGs were enriched in immune-related biological processes, including immune response-activating cell surface receptor signaling pathway, and immune response-activating signal transduction; and immune-related molecular functions, such as antigen binding, and immunoglobulin receptor binding; and immune-related cell components, such as external side of plasma membrane and immunoglobulin complex ( Figure 6A). KEGG pathway analysis showed that the DEGs were enriched in ferroptosis-related immune pathways, such as cytokine-cytokine receptor interaction, cell adhesion molecules, chemokine signaling pathway, and phagosome ( Figure 6B). GSEA indicated that genes in low-risk group were significantly enriched in immune pathways, including Toll-like receptor signaling pathway, chemokine signaling pathway, natural killer cell-mediated cytotoxicity, and antigen processing and presentation ( Figure 7A-D). These results suggest a ferroptosis-related immunosuppressive status in highrisk group.

| Immunity and gene expression
The results of ssGSEA suggested that almost immune cell subpopulations and functions were significantly different among the low-and high-risk groups ( Figure 8A,B). Given the importance of ICI in melanoma, we explored the difference in the expression of immune target genes between the two groups. We found the expression of PDCD-1 (PD-1), CD247 (PD-L1), CTLA4, and LAG3 were significantly higher in low-risk group compared with those in high-risk group ( Figure 8C-F). The results of the spearman correlation analysis showed that the risk score was negatively correlated with PDCD-1 (R = −0.59, p < 0.001, Figure 8G), CD247 (R = −0.67, p < 0.001, Figure 8H), CTLA4 (R = −0.49, p <0.001, Figure 8I), and LAG3 (R = −0.61, p < 0.001, Figure 8J).

| The qPCR analyses
In order to further verify the results of 10 differential expressed genes in the signature, qPCR was used to detect these genes expression at the mRNA level in vitro. These results of qPCR confirmed that the expression of JDP2, TUBE1, PROM2, GPX2, FBXW7, and ARNTL genes in A2058 and A375 cell lines (SKCM tumor cells) were significantly lower than that in the HaCaT cell line (normal skin cells) (Figure 9A-F). Conversely, the expression of HAMP, SLC7A5, IFNG, and CYBB genes in A2058 and A375 cell lines were significantly higher than that in the HaCaT cell line ( Figure 9G-J). These findings were consistent with the differential expression of FRGs genes between normal and tumors tissues, further validating the accuracy and reliability of this signature model.

| DISCUSSION
This study identified FRGs by comparing RNA-seq data from normal samples of GTEx and tumor samples of TCGA. Thirteen differentially expressed FRGs were identified in TCGA-SKCM cohort by univariate Cox regression analysis, and 10 OS-related FRGs in SKCM were selected by LASSO. There were significant differences in risk scores and survival between the high-risk and low-risk groups in the training and validation cohorts, demonstrating the effectiveness and accuracy of the signature in SKCM. GO, KEGG, and GSEA analyses demonstrated that these FRGs were significantly enriched in immune-related function and immune-related pathways. Compared with previous prognostic signature of SKCM, this study is the first time to use FRGs for training and validated in an independent cohort. This study aimed to provide more information for future SKCM studies.
Ferroptosis is an iron-dependent regulatory cell death induced by lipid peroxidation. 29 Ferroptosis F I G U R E 7 Analysis of enriched pathways. KEGG analysis (A-D) of Gene Set Enrichment Analysis in the low-risk groups in skin cutaneous melanoma differed from autophagy, apoptosis, and necrosis regarding cellular, molecular, and biochemical mechanisms and involves the plasma membrane integrity, cytoplasmic and organelle swelling, and moderate chromatin condensation. 30,31 Ferroptosis has been demonstrated to be involved in the pathophysiological processes of many diseases, such as ischemia-reperfusion injury, neurological disorders, blood diseases, kidney injury, and tumors. 11 Ferroptosis is used in cancer therapeutics 32 and regulates signaling pathways in melanoma cells. 33 In this respect, miR-9 inhibited RSL3-and erastin-induced ferroptosis, and miR-9 silencing-induced ferroptosis in these cells. 34 miR-137-regulated necroptosis in melanoma cells, and miR-137 inactivation enhanced the sensitivity of these cells to RSL3-and erastin-induced ferroptosis. 35 In immunotherapy, activated CD8 + T cells enhanced ferroptosis-mediated lipid peroxidation in tumor cells by downregulating the expression of SLC7A11and SLC3A2 through the release of γ-IFN. 36 Melanoma cells depended on oleic acid for freedom from acsl3-mediated ferroptosis in lymph and increased cell survival and metastatic potential. 19 These data demonstrate that ferroptosis has great potential of in SKCM treatment and prognosis. Nonetheless, the physiological mechanisms of ferroptosis involved in tumorigenesis and metastasis are unclear.
Some studies identified prognostic FRG signatures by searching public databases. For instance, a signature improved diagnosis and accuracy of prognosis prediction of hepatocellular carcinoma. 37 A 19-gene signature accurately predicted the outcome of glioma patients. 20 However, an FRG-based prognostic model for SKCM has not been established. The 10-FRG signature developed in the present study was strongly correlated with the prognosis of SKCM.
The signature contained the genes CYBB, IFNG, FBXW7, ARNTL, PROM2, GPX2, JDP2, SLC7A5, TUBE1, and HAMP and was a strong predictor of OS in training and validation cohort. Differential expression of these genes was validated using qPCR. The genes of HAMP, SLC7A5, CYBB, and IFNG were highly expressed in melanoma cell lines, while JDP2, TUBE1, PROM2, GPX2, FBXW7, and ARNTL were lowly expressed in melanoma cell lines. PROM2, GPX2, PROM2, GPX2, and SLC7A5 were risk factors for OS, while other genes were protective factors. The CYBB and NADPH oxidase genes exhibited sex differential expression in multiple sclerosis. 38 CYBB was associated with poorer disease-free survival and OS in melanoma. 39 IFNG predicted survival and the response to ICIs in melanoma. 40 FBXW7 is a commonly mutated and inactivated tumor suppressor and was shown to increase resistance to anti-PD-1 and improve the clinical response to ICIs in cancer patients. [41][42][43] The clock gene ARNTL inhibited melanoma cell growth and enhanced immunotherapeutic efficacy by improving effector functions, macrophage mitochondrial metabolism, and redox homeostasis. 44,45 PROM2 promoted ferroptosis resistance in tumor cells through stimulating the production of ferritin-containing exosomes and multivesicular bodies and increasing ferritin export. 46 GPX2 suppressed ferroptosis by silencing lipoxygenases and reducing lipid hydroperoxides. 47 Although the prognostic value of JDP2, SLC7A5, TUBE1, and HAMP in SKCM is currently unknown, the importance of these four genes should not be underestimated. Our study found that these four genes affected the prognosis of SKCM, which might not be directly administered but was determined by the upstream or downstream of these genes. There is no research on the impact of these four genes on SKCM prognosis, which makes our findings exciting and certainly worthy of further research.
Although ferroptosis has recently been a hot topic of study in tumor, the association between ferroptosis and tumor immunity is incompletely understood. KEGG and GO analyses of DEGs showed that many immune-related pathways and biological processes were enriched, indicating that ferroptosis associated with tumor immunity. Significant differences in immune response-activating signal transduction and immune response-activating cell surface receptor signaling pathway were observed between the low-and high-risk groups. This implied that ferroptotic F I G U R E 9 The differential expression of ferroptosis-related genes was detected by qPCR. (A-F) Compared with the HaCaT cell line, the mRNA of JDP2, TUBE1, PROM2, GPX2, FBXW7, and ARNTL were significantly lower in the A2058 and A375 cell lines, (G-J) while the mRNA of HAMP, SLC7A5, CYBB, and IFNG were significantly higher in the A2058 and A375 cell lines. *p < 0.05, **p < 0.01, ***p < 0.001, ****p < 0.0001 cells release certain signaling molecules such as lipid mediators to attract T and B cells. GSEA in the HALLMARK collection showed that gene sets for antigen processing and presentation, chemokine signaling pathway, natural killer cell-mediated cytotoxic city, and Toll-like receptor signaling pathway were enriched, demonstrating that ferroptosis is involved in cell death, host defense, and innate and adaptive immunity. 48 Cancers are currently being treated using immunotherapy and iron nanoparticles. 49,50 Nanoparticles induce ferroptosis by regulating iron and ROS levels and are a novel strategy for cancer therapy. 50 Immunotherapy has dramatically improved the outcomes of advanced-stage melanoma patients. 51 A recent research found that ferroptosis enhanced the antitumor effect of ICIs, even in ICIs-resistant tumors. 48 Therefore, a novel FRGs signature were constructed to investigate the potential relationship between ferroptosis and ICIs, and to predict immune-checkpoint blockade immunotherapy responses. In this study, the FRGs were revealed to be associated with ICIs (PD-1, PD-L1, CTLA-4, and LAG-3), which indicated that the FRGs signature might be used to predict the response to ICB therapy. In the meantime, the expression levels of these ICIs in low-risk group were higher compared with high-risk group. This suggested that FRGs signature could be applied to guide immunotherapy based on predicted expression level of ICIs.
Our study has several limitations. First, this prognostic signature was constructed using public databases. Therefore, multicenter clinical trials are needed to validate the clinical utility of this model. Second, other genes may be good predictors of SKCM prognosis. Third, the association between risk scores and immune responses has not been experimentally demonstrated.

| CONCLUSIONS
Our study identified a signature with 10 ferroptosisrelated genes. The signature was independently correlated with OS in both the training and validation cohorts, improving the accuracy of predicting the prognosis of SKCM. Nonetheless, additional studies are warranted to determine the mechanism underlying the correlation between ferroptosis-related genes and immunity in SKCM.