m6A-Related lncRNAs Are Potential Biomarkers for the Prognosis of Metastatic Skin Cutaneous Melanoma

Background: The incidence of skin cutaneous melanoma (SKCM) has risen more rapidly than any other solid tumor in the past few decades. The median survival for metastatic melanoma is only six to nine months and the 5°years survival rate of patients with conventional therapy is less than 5%. Our aim was to reveal the potential molecular mechanism in m6A modification of lncRNA and provide candidate prognostic biomarkers for metastatic SKCM. Methods: lncRNAs expression level was obtained by re-annotation in TCGA and CCLE datasets. m6A-related lncRNAs were selected though correlation analysis. Univariate cox regression analysis was used to screen out independent prognostic factors. LASSO Cox regression was performed to construct an m6A-related lncRNA model (m6A-LncM). Univariate survival analysis and ROC curve were used to assess the prognostic efficacy of this model and candidate lncRNAs. Enrichment analysis was used to explore the candidate genes’ functions. Results: We obtained 1,086 common m6A-related lncRNAs after Pearson correlation analysis in both two datasets. 130 out of the 1,086 lncRNAs are independent prognostic factors. 24 crucial lncRNAs were filtered after LASSO Cox regression analysis. All the m6A-LncM and the 24 lncRNAs were related to overall survival. Stratified survival analysis of m6A-LncM showed that the model retains its prognostic efficacy in recurrence, radiation therapy and other subgroups. Enrichment analysis also found that these lncRNAs were immune associated. Conclusion: Here, we obtained 24 crucial lncRNAs that may be potential biomarkers to predict survival of metastatic SKCM and may provide a new insight to improve the prognosis of it.


INTRODUCTION
Based on the anatomical location, melanomas could be subdivided into limbic melanoma, skin cutaneous melanoma (SKCM), and mucosal melanoma (Curtin et al., 2005;Curtin et al., 2006), of which skin cutaneous melanoma is the major subtype of melanoma in Caucasians and the proportion of SKCM is roughly 20% of the Asian population (Chang et al., 1998;Chi et al., 2011). The incidence of SKCM has risen more rapidly than any other solid tumor in the past few decades (Eggermont et al., 2014). In recent years, researchers have made significant progress in understanding the biology, genetics, and treatment of SKCM. However, the prognosis remains poor due to the high rate of invasion and metastasis (Bsirini and Smoller, 2018). SKCM can metastasize extensively to the skin, subcutaneous, lymphatic system, lungs, and other non-pulmonary organs, so the patients with metastasis have poor survival rates and metastasis is also a major obstacle to improving prognosis (Balch et al., 2009;Domingues et al., 2018;Kremenovic et al., 2020). Therapies such as neoadjuvant treatment and targeted therapy have been shown to improve the prognosis of metastatic melanoma to some extent (Bai et al., 2019). The median survival for metastatic melanoma is only 6 to 9°months and the 5°years survival rate of patients with conventional therapy is less than 5% (Agarwala, 2009). Therefore, the search for effective biomarkers for prediction of prognosis and new therapeutic targets is urgent.
M6A is a mechanism of post-transcriptional modification of RNA prevalent in eukaryotic cells . m6A modification is involved in the degradation, translation, splicing, and other processes of mRNA Chen et al., 2019;Liu and Gregory, 2019). m6A is frequently found in the 3′ UTR stop codon region and exon region, respectively (Dominissini et al., 2012;Meyer et al., 2012). The process of m6A modification is regulated by a variety of relevant factors (Jia et al., 2013). m6A regulators can be divided into three types based on previous research: 1) Writers that can recognize RNA and modify m6A, which includes KIAA1429, METTL3, RBM15, METTL14, WTAP, METTL16, and ZC3H13 (Balacco and Soller, 2019); 2) Erasers, which are primarily responsible for the removal of m6A modifications, including FTO and ALKNH5 Pan et al., 2018) and; 3) Readers, including HNRNPC, HNRNPA2B1, YTHDF1, YTHDC1, YTHDF2, YTHDC2, and YTHDF3 (Wang et al., 2018). The Readers could identify RNA methylation modifications and participate in RNA translation, degradation, and other processes . m6A affects many important life processes and is essential for cell division and proliferation, focal death and apoptosis . Numerous studies have confirmed that aberrant m6A modifications play a key role in the genesis and development of a variety of tumors, including SKCM (He et al., 2021;Yu et al., 2021). For example, M6A methyltransferase METTL3 is upregulated in melanoma and modulates melanoma cell invasiveness through MMP2 (Dahal et al., 2019). FTO can promote bladder cancer by regulating the miR-384, MALAT, MAL2 axis through m6A modifications (Tao et al., 2021). The IGF2BP2, LINC00460, and DHX9 complex could promote colorectal cancer proliferation and metastasis by mediating the stability of HMGA1 .
LncRNAs are transcripts with more than 200 nucleotides and with non-coding potential (Ramilowski et al., 2020). The function of lncRNAs remains largely unknown. LncRNA may be involved in the regulation of mRNA expression through translational regulation, histone modifications, or posttranscriptional (Huang et al., 2018). LncRNAs can influence various aspects of tumor cells including survival, proliferation, and migration though participating in gene regulation (Ramilowski et al., 2020). Aberrant expression of lncRNA is associated with tumor malignancy and has been shown to play a key role in the development of numerous cancers including SKCM. For instance, lncRNA TTN-AS1 promotes SKCM development and metastasis by maintaining TTN expression . LncRNA HCP5 inhibits the development of SKCM through regulation of miR-12 expression (Wei et al., 2019).
In this article, we identified m6A modifications related lncRNAs and explored their prognostic ability by bioinformatic analysis and finally obtained potential biomarkers which can predict SKCM prognosis. We also established an important m6A related lncRNA-mRNA regulatory network to provide a new insight to investigate the mechanism of SKCM ( Figure 1).

Ethical Compliance
Public data was used for this study and there are no ethical issues.

Data Sources
Expression data of 484 SKCM patients and corresponding clinical characteristics were obtained from the GDC data portal using gdc-client (https://gdc-portal.nci.nih.gov/). As described in the Introduction, SKCM has a high rate of metastasis and a poor survival rate of patients with metastasis. Thus, our goal is developing the prognostic biomarkers that could predict the poor survival of metastatic SKCM. Then, we filtered 221 patients who with metastatic status for further analysis ( Table 1). We also extracted the expression profile of 38 SKCM cell lines through CCLE database. All the gene expression profiles were quantified by FPKM and normalized though log2-based transformation. In addition, the expression level of 21 m6A-related genes (FTO, ALKBH5, RBM15, RBM15B, METTL3, METTL14, METTL16, WTAP, VIRMA, ZC3H13, HNRNPC, HNRNPA2B1, IGF2BP1, IGF2BP2, IGF2BP3, RBMX, YTHDC1, YTHDC2, YTHDF1, YTHDF2, and YTHDF3) (Zaccara et al., 2019) were constructed from the two datasets, respectively.

Correlation Analysis Between lncRNAs and m6A-Related Genes
Pearson analysis was used to evaluate the correlation between these genes and lncRNAs based on the expression level of lncRNAs and the 21 m6A-related genes. LncRNAs with p value < 0.05 and an absolute Pearson correlation coeffcient (PCC) ≥ 0.3 were selected as m6A-related lncRNAs. We then also used Spearman correlation to double-check the correlation (p value < 0.05, absolute correlation coefficient ≥0.3). LncRNAs that were significantly associated with at least one of the 21 m6A genes in both data sets were selected for subsequent analysis.

Predict the Interactions Between lncRNAs and m6A Regulators and Predict m6A Modification Sites of lncRNAs
RNAInter (RNA Interactome Database) was used to predict the interactions between these lncRNAs and m6A regulators based on the RNA sequence (Lin et al., 2020). SRAMP database was used to predict m6A modification sites of the lncRNAs (Zhou et al., 2016).

Univariate Cox Regression Analysis and Lasso Analysis
The univariate cox regression analysis was used to screen out the prognostic lncRNAs. Subsequently, the R package glmnet (Friedman et al., 2010) was used to construct an m6A-related FIGURE 1 | Overview of the comprehensive analysis. We first obtained 12,711 lncRNAs expression level by re-annotation the TCGA and CCLE SKCM data. Then, 1,086 lncRNAs which correlated with at least one of the 21 m6A genes in both datasets after correlation analysis were filtered. Then, 130 out of the 1,086 lncRNAs are independent prognostic factors after Univariate Cox regression and had been selected for further analysis. We used LASSO Cox regression analysis to further filter a model that contained 24 crucial lncRNAs and have an excellent prognostic efficacy. Stratified survival analysis further testified the model's prognostic efficacy. Finally, enrichment analysis, WGCNA, correlation analysis were used to explored these lncRNAs' potential function in SKCM.
Frontiers in Molecular Biosciences | www.frontiersin.org May 2021 | Volume 8 | Article 687760 lncRNA prognostic model (m6A-LncM) of SKCM patients by LASSO Cox regression. The riskscore of the LASSO regression model could be calculated as follow: where E i is the expression value of the i gene in the model, and βi is the coefficient calculated by LASSO.

Differential Expression Analysis
Based on the riskscore of m6A-LncM, we classified the 221 SKCM patients into high or low risk scores groups. Then, differential expression analysis between the high and low groups was performed using the limma R package. The differentially expressed genes (DEGs) were identified with an adjusted p value < 0.05 and an absolute log2 fold change ≥0.585 (1.5 fold change).

Survival Analysis
Kaplan-Meier analysis was used to evaluate the prognostic efficiency of m6A-LncM and candidate lncRNAs. Survival curves reflect the relationship between the survival model or lncRNA expression level and SKCM patients' survival status though the survival R package. The pROC package was used to calculate the area under the ROC curve (AUC) of the m6A-LncM to assess the prognostic efficiency of it. We also compare the overall survival information between different subgroups of SKCM patients based on the riskscore of the model. The subgroups separated by the following features: age (≤57 or >57°years), neoadjuvant treatment (Yes or No), gender (male or female), melanoma clark level (I, II, III or IV, V), Recurrence (Yes or No), pathologic stage (I, II or III, IV), neoplasm cancer status (with tumor of tumor free), radiation therapy (Yes or No), and tumor issue site (regional or distant). Caret R package was used to randomly separate the 221 samples in two self-dependent test datasets (1:1) which were used to testify the survival efficacy of these lncRNAs.

WGCNA Network Construction and Module Identification
First, the co-expression network of m6A-LncM-related DEGs was constructed by an automatic network construction function in WGCNA R package. Second, co-expression modules were detected by the hierarchical clustering function. Then, modules were associated with clinical characteristics by calculating gene significance (GS) and module membership (MM). Finally, the candidate module with key genes were selected for further analysis.

Construction of the Co-Expression Network
The co-expression network between m6A related genes, m6A-LncM and riskscore-based DEGs was constructed based on the PCC that calculated by Pearson analysis. Those dysregulated lncRNA-mRNA pairs with an absolute PCC ≥ 0.3 and p value < 0.05 were selected to construct the co-expression network.

Function Enrichment Analysis
All the DEGs in the study were extracted for further functional enrichment by using the clusterProfile R package and Metascape software. The lncRNAs' potential functions were obtained from ImmLnc database . Functions with a false discovery rate <0.05 were selected.

Screen m6A-Related lncRNAs in Metastatic SKCM Patients and Cell Lines
A total of 12,711 lncRNAs' expression levels were re-annotated in 221 SKCM patients and 38 SKCM cell lines. Then, we got the expression levels of 21 m6A-associated genes in the two datasets separately and performed Pearson analysis between the 21 genes and the 12,711 lncRNAs. The lncRNAs which exceed the threshold value (p value < 0.05 and |PCC| > 0.3) were defined as m6A-related lncRNAs. Finally, we obtained 1,479 lncRNAs that were significantly associated with m6A in 221 SKCM patients and 8,263 eligible ones in 38 SKCM cell lines. Finally, 1,086 lncRNAs that were significantly associated with at least one of the 21 m6A genes in both data sets were selected for subsequent analysis (Figure 1 and Supplementary Table S1).

Identification of Potential Prognostic lncRNAs and Construct the m6A-LncM
Combined with clinical information, we used univariate Cox regression analysis to filter independent prognostic lncRNAs from the 1,086 lncRNAs related to m6A (p value < 0.05). We obtained 130 out of the 1,086 lncRNAs which were significantly associated with overall survival (OS) of metastatic SKCM patients (Supplementary Table S2). Next, LASSO analysis was used on the 130 prognostic lncRNAs to generate an m6A-associated lncRNA model (m6A-LncM) which contains 24 lncRNAs (Figures 2A,B). The risk score of each sample in the dataset was calculated based on the coefficient of the 24 lncRNAs ( Figure 2C). The 24 lncRNAs' PCC were showed in Figure 2D. We then also used Spearman correlation to double-check the correlation. The result showed that all the 24 lncRNAs were significantly associated with at least one of the 21 m6A genes (p value < 0.05, absolute correlation coefficient ≥0.3) (Supplementary Figure S1A and Table 2).
Then, based on the median risk score, we divided the patients into high-and low-risk score subgroups. The survival curve hinted that patients whose risk scores were high had a worse survival rate ( Figure 2E). The ROC curves also hint that the m6A-LncM had a good efficacy to predict overall survival of metastatic SKCM patients ( Figure 2F). We also performed multivariate survival analysis (cox proportional hazard models by combining the expression of the 24 lncRNAs and OS) to testify the survival efficacy of these lncRNAs. The result showed that they remained significant in multivariate survival models (Supplementary Figure S1B).

The 24 lncRNAs Have m6A Modification Sites and Have Interactions with m6A Regulators
We used RNAInter (RNA Interactome Database) to predict the interactions between these lncRNAs and m6A regulators based on the RNA sequence. We found that all the lncRNAs have interactions with m6A regulators (Supplementary Table S3). On the other hand, we also used SRAMP database to predict m6A modification sites based on the RNA sequences of these lncRNAs. We found that most of the lncRNAs have potential m6A modification sites except RP11-247L20.4 (Supplementary Table S4). The results of the two tools showed that the 24 lncRNAs are regulated by m6A modification and m6A regulators. It also illustrated that the correlation analysis that we performed is credible.  Table 2). The two self-dependent test datasets also verified their survival value. Both these 24 lncRNAs have good prognostic efficacy ( Supplementary  Figures S2, S3).

Stratified Survival Analysis of the m6A-LncM
To further assess the prognostic efficacy of m6A-LncM, we analyzed how the model's risk scores differed across various clinical traits and found that this model could distinguish melanoma clark level, recurrence, neoplasm cancer status, radiation therapy, and tumor issue site ( Figure 4A). Therefore,  we performed the stratified survival analysis to explore whether the model could predict OS in these subgroups. Interestingly, the higher risk metastatic SKCM patients had worse survival rate in all the melanoma clark levels ( Figures 4B,C). Similarly, patients with high riskscores had a poor prognosis, regardless of recurrence ( Figures 4D,E). When the neoplasm cancer status of patients is with tumor, the high riskscore was also associated with a lower survival rate ( Figure 4F). We also confirmed that m6A-LncM retained its prognostic efficacy of patients with radiation therapy (Figures 4G,H). Patients with higher risk also had a worse survival rate in both the regional lymph node and distant metastasis subgroups ( Figure 4I). These data proved that the model could be a novel predictor and performed an excellent prognostic efficacy in metastatic SKCM and may help to improve the prognosis of this cancer.

Enrichment Analysis of Key lncRNAs
To further confirm the potential functions of the m6A-LncM and crucial lncRNAs in SKCM, we classified the 221 patients into high-and low-risk score groups based on the model and obtained differentially expressed genes which are regulated by it ( Figure 5A). The enrichment analysis showed that the upregulated DEGs were enriched in the Melanogenesis pathway, this further confirms that our model plays an important role in SKCM. Other cancer-related signaling pathways were also enriched, include Calcium signaling pathway, Wnt signaling pathway, Hippo signaling pathway, MAPK signaling pathway, and so on ( Figure 5B). The downregulated DEGs were enriched in numerous immune-associated pathways such as Th1 and Th2 cell differentiation, intestinal immune network for IgA production and NF−kappa B signaling pathway ( Figure 5C). Based on the result, we conjecture these lncRNAs may play a crucial role in the immune process of SKCM patients. Then, we used ImmLnc database to further testify whether the 24 lncRNAs are related to immune response. As showed in Figure 5D we confirmed these crucial lncRNAs are also enriched in multiple immune-associated pathways. All of these suggest that the m6A-LncM plays an important role in the prognosis of SKCM and the lncRNAs which we obtained also involved in immune response.

Identification of m6A-LncM-Regulated Gene Modules Associated with Clinical Traits
After differential analyses, we selected the m6A-LncM associated DEGs to construct a gene co-expression network by WGCNA. The soft threshold power was set at 6 for further analysis (Supplementary Figure S4). Next, we constructed the gene network and identified modules using the one-step network construction function of the WGCNA R package. To cluster splitting, the minimum module size was set at 30 and three modules (blue, brown and turquoize) were generated ( Figure 6A). We then analyzed the relationship between these modules. It hinted that there are two clusters over the three modules ( Figure 6B).
Subsequently, we correlated modules with patients' clinical traits to find the crucial genes that associated with clinical characteristics ( Figure 6C). The results showed that the turquoize module was significantly associated with melanoma clark level and neoplasm cancer status ( Figures 6D,F). The brown module was associated with recurrence and tumor issue site ( Figures 6E,G). The blue module was not associated with clinical phenotypes. Finally, we executed an enrichment analysis of the DEGs in the two functional modules and it showed that these key genes were also enriched in immune-associated pathways, such as activation of immune response, alpha-beta T cell activation and lymphocyte migration (Supplementary Figure S5). It may help us to further explore the potential regulatory mechanism of these critical lncRNAs in metastatic SKCM.

Construct a m6A-Regulated lncRNA-mRNA Co-Expressed Network
To better understand the m6A mediated regulation of lncRNA, we constructed a co-expressed network based on the PCC between the m6A genes, m6A-lncM, and m6A-LncM regulated DEGs. Ultimately, 21 m6Agenes, 24 lncRNAs, and 1,251 mRNAs were selected in the network (Figure 7). The colors of the nodes in the network represent their function in metastatic SKCM.

DISCUSSION
In recent years, integrate multiple omics analysis has been widely used in cancer research (Chen et al., 2020b). There are also many important studies involved in the prognosis of SKCM. Anjali et al. developed a webserver to predict the survival of SKCM patients  (Dhall et al., 2020). Li et al. also predicted the metastatic progression of melanoma based on mRNA expression status (Li et al., 2015). The machine learning model has also been performed to predict primary or metastatic SKCM patients (Bhalla et al., 2019). LncRNAs were also found to play a key role in SKCM. Several studies found that some lncRNAs could distinguish the subtypes of SKCM and predict their survival (Ma et al., 2017;Yang et al., 2018). Overexpressed lncRNA HCP5 could also decrease SKCM cell malignancy in vitro which may be upregulated by RARRES3 (Wei et al., 2019). m6A modification has been found that plays an important role in the occurrence and development of tumor. For example, sublethal heat treatment increases EGFR m6A modification near the 5′UTR region and promotes its binding to YTHDF1, thereby enhancing the translation of EGFR mRNA and promoting hepatocellular carcinoma progression (Su et al., 2021). METTL3 dependent m6A modification could regulate the differentiation of T follicle helper cells (Yao et al., 2021). WNT7B-m6A-TCF7L2 positive feedback loop could promote gastric cancer progression and metastasis . The m6A regulators have also been found by many studies to influence the progression and prognosis of cancers though regulating crucial lncRNAs. For instance, LINC00958 regulated by m6a modification can promote breast cancer tumorigenesis through miR-378a-3p and YY1 axis (Rong et al., 2021). A novel hypoxic long-stranded noncoding RNA KB-1980E6.3 maintains stemness in breast cancer stem cells by interacting with IGF2BP1 to promote the stability of c-Myc mRNA (Zhu et al., 2021). MALAT1 promotes thyroid cancer progression by binding to miR-204 and upregulating IGF2BP2, thereby affecting miR-204/IGF2BP2/m6A-MYC signaling (Ye et al., 2021). In summary, m6A modification of lncRNA is not negligible in tumor progression, however, there are no studies on m6A modification of lncRNAs in SKCM.
Here, we identified a prognostic model including 24 m6Aassociated lncRNAs from 221 metastatic SKCM patients. RP11-480I12.5 was reported to promote breast cancer growth and tumorigenesis by inhibiting mir-29c-3p mediated degradation of AKT3 and CDK6 (Lou et al., 2020). LncRNA JPX can promote cervical cancer progression by regulating the miR-25-3p/SOX4 axis (Chen et al., 2020a). These lncRNAs were subsequently found to be closely related to immunity in our study. Interestingly, most of the 24 lncRNAs are novel to cancer research. Compared to the previous lncRNA studies in SKCM, the 24 lncRNAs identified in this study are novel and have the potential for application in clinic. Therefore, we hope that our findings will help to identify potential prognostic lncRNAs regulated by m6A and thus provide ideas to improve the poor prognosis of metastatic SKCM.

Limitations of the Study
Although we screened several novel m6A-related biomarkers, their specific relationship with m6A modifications needs to be verified experimentally. Meanwhile, due to the lack of lncRNA sequencing data and corresponding clinical information for SKCM, our model was not validated using an independent validation set. Subsequent studies may focus on the development of lncRNA sequencing in SKCM to provide strong evidence for screening more potential biomarkers.

CONCLUSION
This is the first study to explore the role of m6A-modified lncRNAs in the prognosis of metastatic SKCM. In this research, we identified the stable correlated m6A-lncRNA pairs in SKCM patients and SKCM cell lines. Then, we constructed a m6A-LncM by lasso regression and found potential prognostic lncRNAs by survival analysis. Twentyfour lncRNAs with independent prognostic efficacy were obtained in the analysis. Next, we comprehensively studied the relationship between these lncRNAs and clinical traits and explored their potential functions. Finally, we delineated a network mediated by these lncRNA, m6A genes and the m6A-LncM regulated differential genes. This study helps to identify potential prognostic targets for metastatic SKCM to improve its poor prognosis.

DATA AVAILABILITY STATEMENT
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article/ Supplementary Material.