Identification of m6A-related lncRNAs-based signature for predicting the prognosis of patients with skin cutaneous melanoma

Background: Skin cutaneous melanoma (SKCM) is one of the fastest developing malignancies with strong aggressive ability and no proper curative treatments. Numerous studies illustrated the importance of N6-methyladenosine (m6A) RNA modification to tumorigenesis. The aim of this study was to identify novel prognostic signature by using m6A-related lncRNAs, thus to improve the survival for SKCM patients and guide SKCM therapy. Methods: We downloaded the Presentational Matrix data from The Cancer Genome Atlas (TCGA) and analyzed all the expressed lncRNAs among 468 SKCM samples. Pearson correlation analysis was performed to assess the correlations between lncRNAs and 29 m6A-related genes. Least absolute shrinkage and selection operator (LASSO), univariate and multivariate Cox regression analysis were performed to construct m6A-related lncRNAs prognostic signature (m6A-LPS). The accuracy and prognostic value of this signature were validated by using receiver operating characteristic (ROC) curves, Kaplan-Meier (K-M) survival analysis, univariate COX or multi-variate COX analyses. After calculating risk scores, patients were divided into low-and high-risk subgroups by the median value of risk scores. Results: A total of 2973 lncRNAs were found expressed among SKCM tissues. Prognostic analysis showed that 98 lncRNAs had a significant effect on the survival of SKCM patients. The m6A-LPS was validated using K-M and ROC analysis and the predictive accuracy of the risk score was also high according to the AUC of the ROC curve in training and testing sets. A nomogram based on tumor stage, gender and risk score that had a strong ability to forecast the 1, 2-, 3-, 5-year OS of SKCM patients confirmed by calibrations. Enrichment analysis indicated that malignancy-associated biological processes and pathways were more common in the high-risk subgroup. Conclusion: Collectively, m6A-related lncRNAs exert as potential biomarkers for prognostic stratification of SKCM patients and may assist clinicians achieving individualized treatment for SKCM.


Introduction
Skin cutaneous melanoma (SKCM) is the least common malignant accounting for 2% of total skin tumors, but is an aggressive skin cancer, causing about 72% of deaths in skin carcinoma [1,2].Nowadays, surgical resection is usually preferred to treat primary melanoma patients [3].Although the five-year relative survival rate for early stage melanoma is more than 95% [4], the reported survival for advanced stages of melanoma is rarely longer than a year.To make things worse, there are no proper effective therapies to treat advanced stages of SKCM [5], which are much more aggressive and not sensitive to radiotherapy and chemotherapy [6].The key to improved survival remains early detection.Although extensive study has explored the mechanism of recurrence and metastasis, the tumorigenesis of cutaneous melanoma remains unclear [1].Exploring the tumorigenesis mechanism may help identify prognostic biomarkers that could serve to guide cancer therapy.
Over the last decade, it has been reported that long noncoding RNAs (lncRNAs) can mediate carcinogenesis, which consist of more than 200 nucleotides and have no protein coding potential [7].One of the major functions of lncRNAs is to regulate gene expression by interacting with microRNAs (miRNAs), thus preventing the latter from inducing translational repression and mRNA degradation via interacting with mRNA 5ʹ or 3ʹ UTRs with the competing endogenous RNA (ceRNA) mechanism [8].By regulating the expression of oncogenes or cancer suppressive genes, lncRNA also participates in oncogenesis and cancer progression and has garnered significant attention in the area of cancer research [9].These evidences make lncRNAs not only as potential biomarkers for diagnosis or prognosis but also as possible therapeutic targets for tumor treatment.Nowadays, several cancer-promotive or -suppressive lncRNAs have been reported in SKCM, such as the cancer-suppressive lncRNA HCP5 [10] and the tumor-promoter LINC00173 [11].Research on the lncRNAs in SKCM is still at its dawn, in-depth exploration of lncRNAs' functions in the formation and progression of SKCM is essential for identifying new targets for the diagnosis and management of this fatal disease.
Methylation occurs at the sixth position of nitrogen atoms of adenosine at the post-transcriptional level with S-adenosylmethionine serving as the methyl donor for N6-methyladenosine (m6A) formation, which is termed m6A modification [12].The reversible regulation of m6A modification is performed by the so-called 'reader', 'writer' and 'eraser' [13].Existing studies have identified m6A-related genes YTHDF1 and HNRNPA2B1 were altered in melanoma and might influence the development of the disease through signaling pathways such as p53 [14], and Wnt and PI3K-Akt signaling [15].Much evidence suggests that the roles of m6A modification in the dysregulation of lncRNAs in cancers, including colorectal cancer [16], liver cancer [17] and glioblastoma [18].However, few studies have been performed to explore the mechanisms underlying how m6A modifications contribute to lncRNA-dependent SKCM occurrence and development.The m6A-related lncRNAs may also have prognostic value for SKCM.
In the present study, we analyzed all expressed lncRNAs among SKCM samples from TCGA dataset to screen 2973 lncRNAs.Our results mainly identified 98 lncRNAs which were found closely related to overall survival (OS) in SKCM patients and 22 prognostic lncRNAs were closely correlated with 20 known m6A-genes.Subsequently, based on the result of univariate and LASSO analysis, we constructed and validated a m6A-related lncRNA risk prognostic signature by using six selected m6A-related lncRNAs, which showed a good performance to stratify the prognosis of SKCM patients.Furthermore, SKCM patients in low-and high-risk subgroups had different prognosis and tumor-related biological processes which were more common in the high-risk subgroup.Additionally, a ceRNA network was built to search the target miRNAs and mRNAs of five selected m6A-related prognostic lncRNAs.Finally, an accurate nomogram was constructed to predict 1-, 2-, 3-, and 5-year OS of SKCM patients and we validated the accuracy to illustrate their prognostic value.

Identification of differently expressed lncRNAs
In our study, gene expression profiles of 468 SKCM patients with clinical information were downloaded from the TCGA data portal (https: //tcga-data.nci.nih.gov/tcga/).All the expressed lncRNAs in the SKCM samples were identified but removed the lncRNA with too low expression (log2(TPM+1) mean <0.5), and K-M survival (p<0.05) and univariate Cox regression (p<0.01)analyses were implemented to screen the prognostic lncRNAs.

Construction of m6A-related lncRNA prognostic signatures (m6A-LPS) and validation
The entire TCGA-SKCM set (n= 468) was randomized as a training set (n= 300) and the testing set (n= 168).The training set was utilized to construct an m6A-related lncRNA model, and the entire set and testing set were applied to validate this established model.After performing LASSO Cox regression analysis to conduct least absolute shrinkage using the R package "glmmet" [20], a multivariate Cox regression analysis was conducted and the corresponding coefficient of these m6A-related prognostic lncRNAs were finally determined in the training set.The m6A-related lncRNA prognostic signature (m6A-LPS) gave patients risk scores which was calculated by a linear combination of Cox-coefficient and m6A-related lncRNA expression (FPKM value): Risk score = ∑ n i=1 (Coef i *Expi).And SKCM patients were classified into low-risk and high-risk group according to the median of risk scores.The K-M method and Univariate COX regression analyzes were utilized to analyze the overall survival (OS) by "survival" R package, and value the prognostic significance of m6A-related lncRNAs in m6A-LPS and risk score.The Receiver operating characteristic (ROC) analysis was conducted to evaluate the prediction efficiency of the m6A-LPS with the R package "pROC".To predict the binding sites of m6A in the lncRNAs, RMBase v2.0 (http://rna.sysu.edu.cn/rmbase/index.php)was used by comparing the coding sequence with its genomic sequence of lncRNAs.To identify the conserved motifs, the MEME program (http://meme-suit e.org/) was performed to count the base content of the m6A potential binding sites on the prognostic lncRNAs.

Gene ontology and pathway enrichment analysis
Among the TCGA cohort, DEGs between the high-risk subgroup and the low-risk subgroup were identified based on the standards of |log2 (FC)| > 1 and p < 0.05 using the R package "limma" [21].The top significantly gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes Pathway (KEGG pathway) on the DEGs were performed through "ggplot2" R package [22].To further explore the underlying pathway of the DEGs in SKCM, we performed gene set enrichment analysis (GSEA) of the upregulated genes in the high-risk subgroup.

Construction of ceRNA networks
The DEGs between 471 SKCM samples and 1 normal skin samples were identified as cancer-specific gene sets (cSEG) based on the "limma" R package [21] with the standard of |logFC| > 1 and P-value < 0.05.The target miRNAs of the m6A-related lncRNAs from m6A-LPS was predicted in lncBase v2.0 dataset and the shared target mRNAs of these target miRNAs was also predicted in the miRbase [23] and TargetScan databases [24], of which these shared target mRNAs were included in the DEGs between high-risk and low-risk subgroups.A ceRNA network was plotted using the software of "Cytoscape" [25].
We performed multiscale embedded gene co-expression network analysis (MEGENA) to study functionally orchestrating clusters of the targeted mRNAs in the ceRNA network with the form of functional networks, which capture well the hub genes of m6A-related lncRNAs fuctions.Furthermore, we used the GSE15605 database to verify the correlation between these hub genes and tumor metastases of SKCM patients.
W. Lin et al.

Independence of the m6A-related lncRNA model
Multivariate and univariate Cox regression analyses were conducted to test whether the prognostic pattern was an independent variable considering other clinical characteristics (tumor stage, gender, age, and BMI) in the patients with SKCM.Subsequently, a nomogram was used to predict cancer prognosis using the "rms" package software.In the TCGA datasets, risk score and other clinic predictors (tumor stage, age, gender and BMI) were included to generate the nomogram which can investigate the chancer of 1-, 2-, 3-and 5-year OS of patients with SKCM.The calibration curves for 1-, 2-, 3-and 5-year OS were performed to evaluate the predictive accuracy of the nomogram.

Human tissue specimens and quantitative real-time PCR (qRT-PCR) and WB analyses
The human SKCM tissues used in this study were approved by Ethnics Committee of Shenzhen People's Hospital.Overall, 8 pairs of surgically resected melanoma tissues and nevi were collected from our hospital.Total RNA was isolated from 5 pairs of melanoma tissues and nevi using the EZpress RNA purification kit (B0004), and cDNA was generated using the Pri-meScript RT Reagent Kit (Takara).Quantitative RT-PCR was performed using SYBR Green PCR master mix (Life Technologies) and a RT-PCR system (Applied Biosystems).GAPDH was used as the control of lncRNAs for PCR product quantification and normalization.
For WB analysis, three pairs of melanoma tissues and nevi tissue extracts were prepared with RIPA lysis buffer (Beyotime Biotechnology).Protein samples were separated by 10% SDS-PAGE and transferred to PVDF membranes (Millipore).After blocking with 5% milk for 1 h at room temperature, the membrane was incubated with the antibodies (including STAP1, CLEC17A, ZC3H12D, RHOH, and GAPDH) in 5% BSA overnight at 4 • C. The membrane was then incubated with a secondary antibody conjugated to a fluorescent tag (Invitrogen).The band signals were visualized in a FluorChem R system (ProteinSimple).Densitometry calculations are calculated using ImageJ.

Statistical analyses
All data in the present study were analyzed using the R statistical package (R version 4.0.3).The R packages includes limma, pheatmap, survival, glmnet, pROC, rms, and ggplot2.Survival status was assessed by the Cox regression analysis.OS was generated by the Kaplan-Meier method and evaluated by the log-rank test.The correlation was determined by Pearson correlation analysis.A two-tailed p < 0.05 was considered statistically significant.

Identification of differentially expressed m6A-related lncRNAs in SKCM
The process flow of this study is shown in Fig. 1.First, we downloaded the transcriptome profiling data of the SKCM project from the TCGA database, of which 468 SKCM samples were included.A total of 2973 lncRNAs were identified among SKCM.Then, 98 lncRNAs were significantly correlated with the OS of SKCM patients through K-M and Univariate COX regression analyzes.We first analyzed and displayed the association of these 98 prognostic related-lncRNAs with 29 m6A-related genes using the Pearson correlation analysis (Figure S1A), of which 22 prognostic lncRNAs had been proved significantly correlated (|Pearson R| > 0.3 and p < 0.001) with 29 m6A-related genes (Fig. 2A).These 22 lncRNAs whose expression value was correlated with one or more of the 29 m6A-related genes was defined as the m6A-related lncRNA in SKCM.Furthermore, these 22 m6A-related lncRNAs were differentially expressed in SKCM samples at different stages (Figure S1B).

Construction of m6A-LPS for SKCM patients
To better predict the clinic outcomes for SKCM by these m6A-related lncRNAs, six of 22 m6A-related lncRNAs were determined using LASSO regression and finally included in the prognostic model, including LINC02528, SEMA6A-AS1, EBLN3P, MIR155HG, LYRM4-AS1, HLA-DQB1-AS1, formed a m6A-related lncRNAs prognostic signature (m6A-LPS) (Fig. 2B and 2C).A prognostic indicator, risk score, was calculated based on the coefficient and expression values of each m6A-related lncRNA (Fig. 2D).Furthermore, as shown in Fig. 2E, these motif-logo showed the base content of the m6A potential binding sites on the three lncRNAs (EBLN3P, MIR155HG, and HLA-DQB1-AS1) was counted using MEME program, respectively.However, the remaining three m6Arelated lncRNAs did not retrieve m6A binding sites, which may be due to the presence of sites that have not been found at present, and the indirect regulation of them by some other genes modified by m6A.

Validation of m6A-LPS risk score in SKCM patients
In order to evaluate the prognostic capability of the m6A-LPS, we calculated risk scores for every patient in the training set and testing set using the uniform formula.The SKCM patients in every set were divided into high-risk subgroup and low-risk subgroup based on the median risk score calculated above, and the OS between two subgroups was compared.The distributions of risk scores and survival time of all SKCM patients in training set (Fig. 3A) and testing set (Fig. 3B) were displayed.These findings suggested that the clinical outcome of patients in the lowrisk group was superior to that of patients in the high-risk group.K-M analysis showed that patients in the low-risk subgroup had a significantly better survival than those in the high-risk subgroup both in training set (p<0.0001, Fig. 3C) and testing set (p<0.0001, Fig. 3D).Subsequently, the Univariate Cox regression analysis was firstly conducted to identify the association of m6A-LPS with OS of SKCM patients  in the entire set.The results exhibited that six m6A-related lncRNAs and risk score were significantly associated with OS (p < 0.001), among which four m6A-related lncRNAs (LINC02528, EBLN3P, MIR155HG and HLA-DQB1-AS1) were protective factors with HR < 1 whereas two m6Arelated lncRNAs (SEMA6A-AS1 and LYRM4-AS1) as well as risk score acted as risky factor with HR > 1 in SKCM patients (Fig. 3E).Notably, only two m6A-related lncRNAs (LINC02528 and MIR155HG) were significantly different in SKCM at different stages (Fig. 3F), highlighting that the two m6A-related lncRNAs could indicate the tumor stages of SKCM.However, due to its strong heterogeneity, these differentially expressed lncrnas do not increase with the increase of tumor stage.The K-M survival curves also confirmed that higher expression of LINC02528, EBLN3P, MIR155HG and HLA-DQB1-AS1 (Fig. 4A), while low expression of SEMA6A-AS1 and LYRM4-AS1 (Fig. 4B) were associated with better OS.To further confirm the potential clinical implications of m6A-related lncRNAs in melanoma immunotherapy, qRT-PCR analysis was conducted to determine the expression levels of the six m6A-related lncRNAs in a cohort of 5 melanoma tissues.The mRNA level of LYRM4-AS1 was prominently increased in tumors, but EBLN3P decreased in tumor when compared with nevi (Fig. 4C).However, the expression levels of the other four lncRNAs were not significantly different between tumors and nevi groups.To validate the optimality, we plotted the 1-, 2-, 3-, and 5-year ROC curves in TCGA entire set, training set and testing set, which suggested that all AUC values were over 0.6 (Fig. 4D).Collectively, the six m6A-related lncRNAs were identified as the independent prognostic factor for SKCM and the m6A-LPS risk score had the distinct prognostic value.

Enrichment analysis identifies biological functions and signaling pathways affected by m6A-related lncRNA in SKCM
For investigating the potential biological process and pathway involving in the molecular heterogeneity between the low-and high-risk subgroups, we identified 1923 (DEGs) [|log2 (FC)|> 1 and p < 0.05] between the low-and high-risk subgroups in the TCGA cohort (Fig. 5A).The result of GO and KEGG analyses showed that the 1332 downregulated DEGs (high-risk vs. low-risk subgroups) had an accumulation of immune-related GO items, including immune response-activating cell surface receptor signaling pathway, immune response-activating signaling transduction, regulation of immune effector process, positive regulation of cell activation, positive regulation of leukocyte activation; in addition, the KEGG pathways hematopoietic cell lineage, cytokinecytokine receptor interaction, chemokine signaling pathway were also enriched these down-regulated DEGs (Fig. 5B).Meanwhile, these 591 upregulated DEGs (high-risk vs. low-risk subgroups) were involved in mitochondrial translation, ncRNA metabolic process, mitochondrial gene expression, cellular respiration and ncRNA processing (GO terms), as well as several KEGG pathways, including oxidative phosphorylation, Alzheimer disease, huntington disease, Pakinson disease, fanconi anemia pathway (Fig. 5C).GSEA revealed that several tumor hallmarks and immune infiltration were enriched in the high-risk subgroup, such as apoptosis, B cell receptor signaling pathway, T cell receptor signaling pathway, and natural killer cell mediated cytotoxicity (Fig. 5D) and so on.These results may give us some insights into the cellular biological effects related to the m6A-LPS in SKCM.

Construction of the ceRNA network and the prognostic nomogram, and evaluation of predictive effectiveness
To further understand how these m6A-related prognostic lncRNAs regulate mRNA expression by sponging miRNAs in SKCM, we constructed a ceRNA network based on these m6A-related lncRNAs.Firstly, a total of 772 DEGs were identified between normal and SKCM samples according the criteria of p < 0.01 (Fig. 6A).Then, 53 target miRNAs were identified based on the lncBase and totally 8176 mRNAs were predicted through two databases (miRBase and TargetScan).As shown in Fig. 6B, Venny showed the intersection of all targeted mRNAs and top 300 mRNAs from DEGs.Ultimately, five lncRNAs, 53 miRNAs and 41 mRNAs were included in the ceRNA network (Fig. 6C).These data may provide us some directions for finding the potential functions of these m6A-related prognostic lncRNAs in SKCM.Furthermore, MEGENA program was performed to analysis the hub mRNAs among the ceRNA regulatory network, and found four key mRNAs (STAP1, CLEC17A, ZC3H12D and RHOH) (Fig. 6D).Interestingly, in the GSE15605 dataset, RhOH and STAP1 were found to have a significant relationship with tumor metastasis of SKCM patients (Fig. 6E).Notably, these four key mRNAs (STAP1, CLEC17A, ZC3H12D and RHOH) were generally highly expressed in tumor tissues comparing to nevi tissues (Fig. 6F).Thus, these findings indicate that the m6A-related lncRNA model may have great prognostic significance for SKCM patients.

Evaluation of the prognostic risk model of m6A-related lncRNAs and clinical features of SKCM, as well as the prognostic nomogram
We performed univariate and multivariate Cox regression analyses to evaluate whether the risk model of six m6A-related lncRNAs had independent prognostic characteristics for SKCM.The HR of the risk score was 3.81 (p < 0.001) in univariate Cox regression analysis.In multivariate Cox regression analysis, the HR was 2.83 (p < 0.001) (Fig. 7A), indicating that the risk model of the six m6A-related lncRNAs was unrelated to clinicopathological parameters, such as tumor stage, gender, age, and BMI.A prognostic nomogram whose indicators include tumor stage, gender and risk score was constructed to predict OS of individual patients with SKCM (Fig. 7B), and the calibration plots showed that the observed vs. predicted rates of 1-, 2-, 3-and 5-year OS showed perfect concordance and the nomogram was of excellent predictive effect, especially predicting short-term survival (1-and 2-year) for SKCM patients (Fig. 7C).These data indicated that the nomogram has a robust and stable ability to predictive the OS for SKCM patients.
Given that these m6A-related lncRNAs were related to tumor immunity, we also performed correction analysis for risk score with infiltration of immune cell subtypes in SKCM using the data from the TIMER database.As shown in Fig. 8A, the heatmap indicated risk score was closely related to the top 10 infiltration of immune cell types in SKCM.The correlation coefficients of risk score between B lineage, monocytic lineage, cytotoxic lymphocytes, T cells, NK cells, myeloid dendritic cells, CD8+T cells and fibroblasts were − 0.49, − 0.65, − 0.54, − 0.6, − 0.47, − 0.48, − 0.59 and − 0.18 (p < 0.001) (Fig. 8B).The above results suggested that the infiltration degree of these immune cell subtypes was significantly negatively correlated with the risk scores based on seven m6A-related lncRNAs in SKCM.

Discussion
The incidence and mortality of SKCM is increasing around the world and ranks 19th among the most common malignant cancers [26].Even though some sophisticated therapeutic regimens have evolved from surgery, targeted therapy, chemotherapy, and immunotherapy, the prognosis of patients with SKCM has not improved significantly [27][28][29][30].Emerging evidences have indicated that m6A RNA modifications can mediate various malignancy-related process, such as onco-transcript expression, tumorigenesis and invasion in lung cancer [31], tumor proliferation in bladder cancer [32], and metastasis in colorectal cancer [33].On the contrary, tumor suppressive role of m6A methylation have also been reported by various studies.For instance, Liu et al. reported that reduced m6A methylation regulates AKT activity to promote the proliferation and tumorigenicity of endometrial cancer [34].Besides, lncRNAs were identified as novel regulators of tumorigenesis and tumor progression [35].A variety of lncRNAs are abnormally expressed and are closely related to the prognosis of patients with SKCM [36,37].However, the prognostic value of m6A-related lncRNAs on tumor survival remains unclear.Here, we sought to identify novel biomarkers which might serve as a critical target for SKCM prognosis and optimizing therapeutic methods.
Firstly, we identified 2973 lncRNAs in SKCM samples.Based on K-M & univariate Cox regression and Pearson correlation analyses, we extracted 22 m6A-related prognostic lncRNA.Next, univariate Cox regression and LASSO analyses were used to identify six m6A-related lncRNAs, including LINC02528, SEMA6A-AS1, EBLN3P, MIR155HG, LYRM4-AS1, HLA-DQB1-AS1.In addition, SEMA6A-AS1 and LYRM4-AS1 were found to be negatively associated with OS time of SKCM patients, indicating that these m6A-related lncRNAs might be oncogenes.In contrast, LINC02528, EBLN3P, MIR155HG and HLA-DQB1-AS1 were positively correlated with patient outcomes, revealing their vital role in inhibiting the progression of SKCM.Low expression of SEMA6A-AS1 was associated with a poor prognosis in patients with HBV-related hepatocellular carcinoma [38].LYRM4-AS1 could regulate the growth of IL-1β-induced chondrocytes and played important role in cervical cancer [39,40].EBLN3P, as an immune-related lncRNA, was proved to be related with immune response suppression and progression in breast cancer [41].MIR155HG was associated with immune infiltration and immune checkpoints in multiple cancers as a prognostic biomarker, includes SKCM [42].HLA-DQB1-AS1 was identified as the immune-related lncRNA and autophagy-related lncRNA, which was used to construct a risk signature for lung adenocarcinoma patients [43,44].Most of the identified six m6A-related lncRNAs in our study were previously significantly associated with the progression of several cancers.However, there are no studies reported their expression and function in SKCM, as well as systematically analyzing the specific prognostic value of m6A-related lncRNAs in SKCM.Therefore, our findings firstly provide some insights into the relationship between m6A-related lncRNAs and SKCM progression, identifying valuable m6A-associated biomarkers for personalized treatment.
The ROC curves and survival analyses confirmed the advanced biological implications of our model to predict the outcomes of SKCM patients.In addition, our risk score model also is successful in predicting SKCM prognosis and has shown to have improved predictive accuracy compared to other clinical features such as tumor stage and gender.Furthermore, Cox regressions evidenced that our selected six m6Arelated lncRNAs and the risk score were independent prognostic factors for SKCM patients.Nomograms constructed based on the clinic risk signatures revealed the credibility of our risk score model in predicting the 1-/2-/3-/5-year OS time of SKCM patients.
GO enrichment analysis indicated that downregulated DEGs between high-and low-risk subgroups were mainly involved in immune response-activating cell surface receptor signaling pathway, immune response-activating signaling transduction, regulation of immune effector process, positive regulation of cell activation, positive regulation of leukocyte activation; meanwhile, KEGG pathway analysis suggested that these DEGs highly expressed in low-risk group could have a significant impact on Hematopoietic cell lineage, Cytokine-cytokine receptor interaction, osteoclast differentiation, chemokine signaling pathway (Fig. 5B).Besides, the upregulated genes in high-risk group were enriched in mitochondrial translation, ncRNA metabolic process, mitochondrial gene expression, cellular respiration and ncRNA processing.KEGG enrichment analysis portrayed that these selected upregulated DEGs were enriched in oxidative phosphorylation, Alzheimer disease, huntington disease, Pakinson disease, fanconi anemia pathway (Fig. 5C).Biological processes involved in differentially expressed genes between high-and low-risk subgroups distinct by the m6A-related LPS constructed by six m6A-related lncRNAs might give us some insights into new targets may be developed to improve therapeutic efficacy for SKCM.
Nonetheless, there are some limitations in this study despite of these findings.External validation by other datasets would be beneficial, and the biological mechanism of m6A-related lncRNAs has not been fully elucidated.First, this study was designed as a retrospective analysis; more prospective research should be performed to verify our results.In addition, our results lack in vitro or in vivo exploration to confirm the reliability of our mechanism analysis.Therefore, we need to conduct several further experiments to prove the mechanistic connections between these selected m6A-lncRNAs and the progression of SKCM.

Conclusion
Accordingly, we constructed a m6A-related lncRNAs prognostic signature and performed a series of bioinformatic analyses to identify six m6A-related lncRNAs significantly associated with the overall survival of patients with SKCM.We also successfully constructed a prognostic risk score model with powerful predictive effects.To the best of our knowledge, this is the first m6A-associated prognostic model to predict the prognosis of SKCM.This work identified a new method for understanding the specific roles of m6A-related lncRNAs in SKCM and highlighted the potential of m6A modification to elucidate clinical prognosis in SKCM patients.We believe that our study makes a significant contribution to the literature as our prognostic model may provide new insights into the pathogenesis, prognosis, and immunotherapy in patients with SKCM.A more in-depth understanding of the pathological mechanism of SKCM might help overcoming this obstacle in SKCM management and improving patients' welfare.

Fig. 2 .
Fig. 2. Identification of m6A-related lncRNAs in SKCM.(A) Heat map of the correlation between 22 differently expressed prognostic lncRNAs and 29 m6A-related genes.(B) The tuning parameters (log l) of OS-related proteins were selected to cross-verify the error curve.According to the minimal criterion and 1-se criterion, perpendicular imaginary lines were drawn at the optimal value.(C) The LASSO coefficient profile of 22 lncRNAs and perpendicular imaginary line were drawn at the value chosen by 10-fold cross-validation.(D) The coefficients of the selected six m6A-related lncRNAs.(E) Conserved binding sites of m6A among the selected three m6A-related lncRNAs (EBLN3P, HLA-DGB1-AS1 and MIR155HG) were predicted by tool MEME.

Fig. 3 .
Fig. 3. Risk score assessment model for prognosis prediction in training and testing sets.(A) The distributions of survival outcome (upper) and risk scores (bottom) of each SKCM case in training set.(B) The distributions of survival outcome (upper) and risk scores (bottom) of each SKCM case in testing set.(C) Kaplan-Meier survival curves of the OS of patients in the high-and low-risk groups for the training set.(D) Patients in the low-risk subgroup experienced a longer survival time tested by the Kaplan-Meier analysis in testing set.(E) A forest map showed six m6A-related lncRNAs and risk score identified by Univariant Cox proportional hazard regression in the stepwise method.(F) Expressions of six m6A-related lncRNAs in SKCM samples at different stages.

Fig. 4 .
Fig. 4. Prognostic value of m6A-related lncRNAs.(A-B) Kaplan-Meier curves showing that patients with different expression levels of the selected m6A-related lncRNAs had different overall survival.Validation of six m6A-related lncRNAs in melanoma tissues and nevi.(D) The 1-, 2-, 3-, and 5-year ROC curves of the optimal model suggested that all AUC values were over 0.60 in the entire, training and testing sets.

Fig. 5 .
Fig. 5. Risk score involved in definite biological functions in SKCM.(A) Volcano plot of the differentially expressed genes (DEGs) between the low-and high-risk subgroups.(B) Bubble chart of top 8 classes of Gene Ontology (GO) enrichment terms and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways about downregulated DEGs.(C) Bubble chart of top 8 classes of GO enrichment terms and KEGG pathways about up-regulated DEGs.(D) Gene set enrichment analysis (GSEA) indicating that tumor hallmarks were enriched in the high-risk subgroup.

Fig. 6 .
Fig. 6.Assessment of the prognostic risk score of m6A-LPS.(A) Volcano plot of the differentially expressed genes (DEGs) between the 471 SKCM samples and 1 normal sample from TCGA dataset.(B) Venny showed the intersections of top 300 DEGs and 8176 targeted mRNAs.(C) A ceRNA network of the selected five m6Arelated lncRNAs (blue) and 53 target miRNAs (red) and 41 mRNAs (green).(D) MEGENA network showed the hub genes from the ceRNA network.(E) Different expressions of four hub genes in metastasis or not SKCM patient from the GSE15605 dataset.(F) The expression of four mRNA (STAP1, CLEC17A, ZC3H12D and RHOH) in melanoma tissues and nevi.

Fig. 7 .
Fig. 7. Construction and evaluation of a prognostic nomogram (A) Univariate and multivariate analyses of the clinical characteristics and risk score with the OS.(B) A prognostic nomogram based on tumor stage, gender, age, BMI and the risk score.(C) The calibration assessment by calibration curves of the nomogram predicts the probability of the 1-, 2-, 3-, and 5-year OS.

Fig. 8 .
Fig. 8. Correlation between six m6A-related lncRNAs and infiltration of various immune cell subtypes (A) Heatmaps of abundance for top 10 immune cell subtypes between low-risk group and high-risk group in the entire sets.Warm colors represented high abundance, while cold colors represented low abundance.(B) The correlation between the risk scores and 10 immune cell subtypes.