A Novel Autophagy-Related lncRNA Gene Signature to Improve the Prognosis of Patients with Melanoma

Objective Autophagy and long noncoding RNAs (lncRNAs) have been the focus of research on the pathogenesis of melanoma. However, the autophagy network of lncRNAs in melanoma has not been reported. The purpose of this study was to investigate the lncRNA prognostic markers related to melanoma autophagy and predict the prognosis of patients with melanoma. Methods We downloaded RNA sequencing data and clinical information of melanoma from the Cancer Genome Atlas. The coexpression of autophagy-related genes (ARGs) and lncRNAs was analyzed. The risk model of autophagy-related lncRNAs was established by univariate and multivariate Cox regression analyses, and the best prognostic index was evaluated combined with clinical data. Finally, gene set enrichment analysis was performed on patients in the high- and low-risk groups. Results According to the results of the univariate Cox analysis, only the overexpression of LINC00520 was associated with poor overall survival, unlike HLA-DQB1-AS1, USP30-AS1, AL645929, AL365361, LINC00324, and AC055822. The results of the multivariate Cox analysis showed that the overall survival of patients in the high-risk group was shorter than that recorded in the low-risk group (p < 0.001). Moreover, in the receiver operating characteristic curve of the risk model we constructed, the area under the curve (AUC) was 0.734, while the AUC of T and N was 0.707 and 0.658, respectively. The Gene Ontology was mainly enriched with the positive regulation of autophagy and the activation of the immune system. The results of the Kyoto Encyclopedia of Genes and Genomes enrichment were mostly related to autophagy, immunity, and melanin metabolism. Conclusion The positive regulation of autophagy may slow the transition from low-risk patients to high-risk patients in melanoma. Furthermore, compared with clinical information, the autophagy-related lncRNA risk model may better predict the prognosis of patients with melanoma and provide new treatment ideas.


Background
Melanoma is a highly malignant tumor characterized by strong invasiveness and metastasis [1]. Although it accounts for a low proportion of skin cancer cases, it is associated with a high mortality rate compared with other skin cancers [2]. Its occurrence is mainly caused by long-term exposure to ultraviolet radiation, sunburn, physical stimulation, and other factors [3]. These risk factors can promote the malignant transformation of melanocytes. At present, the treatment of melanoma mainly involves surgery, radiotherapy, chemotherapy, and other means to kill tumor cells and promote tumor cell apoptosis [4][5][6].
However, the mechanism of melanoma is regulated by a complex molecular regulatory network; hence, the effectiveness of currently available treatments is limited. In addition, according to the literature, autophagy-related genes (ARGs) and TNM staging can predict the prognosis of patients with melanoma [7,8]. However, due to the heterogeneity of tumors, TNM staging cannot accurately predict the prognosis of patients. Therefore, early diagnosis and treatment of patients with melanoma are particularly important.
Autophagy is a catabolic process capturing proteins or organelles that need to be degraded and maintaining cell stability through lysosomes [9]. At present, studies have shown that autophagy plays an important role in the formation of tumors [10]. In the early stages of cancer, autophagy can inhibit tumor growth [11]. However, with the development of cancer, cancer cells show a high degree of autophagy dependence, and autophagy plays a role in promoting the growth of such cells [12]. In addition, autophagy can be used as a prognostic marker for a variety of cancers, while inhibition or activation of autophagy can destroy tumor cells [13][14][15]. Therefore, autophagy also plays an important role in the diagnosis and treatment of cancer.
Long non-coding RNAs (lncRNAs) are non-coding protein RNAs, with a length >200 bp. In recent years, lncRNAs have been identi ed as important factors in the regulation of cancer, autoimmune disease, cardiovascular disease, etc [16]. They can affect gene expression by regulating biological processes, including the binding of transcription factors, regulation of chromatin structure, induction of miRNA, etc., and subsequently regulate the biological function of cancer cells [17]. Current studies have shown that lncRNAs account for a large proportion of melanoma genomes compared with genes that encode proteins [18]; however, the mechanism of melanoma regulation remains unclear. Moreover, lncRNAs can also be used as diagnostic and prognostic markers to predict the prognosis of patients with melanoma [19]. Therefore, the role of lncRNA in melanoma is self-evident.
At present, numerous studies have con rmed the role of autophagy in melanoma [20,21]. The increase of autophagy promotes the progression of melanoma [21], but it can also inhibit its growth [22]. Moreover, these studies are based on ARGs to diagnose or predict the prognosis of patients with this disease. In addition, recent studies have found that lncRNAs can directly or indirectly regulate autophagy [23]. They can directly activate the protein that initiates autophagy or indirectly regulate the expression of ARGs through the competing endogenous RNA mechanism [24]. Therefore, the network between lncRNAs and ARGs may be the key to elucidating the autophagy mechanism of melanoma. However, there are few studies on autophagy-related lncRNAs in this setting, and fewer autophagy-related lncRNAs have been used to evaluate the prognosis of patients with melanoma. Therefore, it is important to clarify the autophagy mechanism of melanoma related to lncRNAs and identify new molecular therapeutic targets.
In this study, we downloaded melanoma expression data from The Cancer Genome Atlas (TCGA) and obtained seven lncRNAs related to autophagy. Further multivariate COX regression analysis, construction of autophagy-related co-expression network, and gene set enrichment analysis (GSEA) were performed.
Most importantly, we constructed a risk model that can more accurately diagnose and predict the prognosis of patients than clinical information, suggesting that the autophagy-related lncRNA signature is a reliable predictor for patients with melanoma.

Data preparation
We downloaded the FPKM data and related clinical information of 471 melanoma samples from TCGA (https://www.cancergenome.nih.gov/). Furthermore, the mRNA matrix and lncRNA matrix of melanoma were obtained. The TCGA database is a public database containing 33 tumor samples and matching normal samples. At present, the research on ARGs is relatively perfect, and the Human Autophagydedicated Database (HDAb; (http://www.autophagy.lu/) is the rst public database dedicated to human autophagy [25]. Hence, we downloaded ARGs from the HDAb.
Co-expression analysis of autophagy genes After obtaining the autophagy genes from the HDAb website, we compared the mRNA matrix of melanoma with this gene set, and further extracted the autophagy-related mRNA matrix of melanoma. Finally, the co-expression network of the mRNA matrix and lncRNA matrix related to melanoma autophagy was analyzed using the limma package of the R software (version 4.0.2). |Correlation coe cient| ≥5 and p-value <0.001 were the best screening criteria.

Survival analysis of lncRNAs associated with autophagy
Prior to the multivariate COX analysis of lncRNAs, we evaluated the signi cance of lncRNAs in survival through single-factor COX analysis, and conducted a univariate COX analysis of the subsequent lncRNAs. We combined the survival time and survival status of the clinical information with lncRNAs, and subsequently performed Kaplan-Meier survival analysis of the autophagy-related lncRNAs through the survival package of the R software. P-values <0.05 denoted statistically signi cant differences.

Cox regression model analysis
We further screened lncRNAs that can be used for the prediction of prognosis of patients. We used the survival package of the R software to perform a stepwise COX regression analysis of lncRNAs with survival signi cance and obtain a riskScore. The patients were divided into high-and low-risk groups according to their riskScore. At present, the combination of multiple genes has demonstrated potential for predicting the prognosis of patients. Thus, the lncRNAs were screened through survival analysis and receiver operating characteristic (ROC) curve analysis using the survivalROC package. Finally, based on the median value of the riskScore, we drew the riskScore, survival diagram, and heatmap using the pheatmap package. To better predict the prognosis of patients, we used a survival package combined with clinical information of patients and the risk model constructed in this experiment for univariate and multivariate COX analyses, and obtained the related forest plots. Furthermore, the ROC curve related to clinical information was drawn to evaluate the prognostic value of the model. P-value < 0.05.

Construction of a co-expression network
The Cytoscape software (version 3.7.2) was used to visualize the lncRNA-related autophagy coexpression networks [26]. The Sankey diagram can directly show the autophagy-related co-expression network of the high-and low-risk groups; hence, this diagram was constructed using the ggalluvial package of the R software. Hazard ratios (HR) >1 and <1 indicated high and low risk, respectively. GSEA GSEA enables researchers to better understand the potential pathogenic mechanism in the high-and lowrisk groups. Therefore, we analyzed the Gene Ontology (GO) and the Kyoto Encyclopedia of Genes and Genomes (KEGG) of the high-and low-risk groups through GESA (version 4.0.3) [27]. Normalized enrichment score ≥1, nominal p-value ≤0.05, and false discovery rate q-value ≤0.25 denoted statistical signi cance.

Co-expression analysis of autophagy genes
We downloaded a total of 232 ARGs from the HADb database. In addition, by screening the mRNA matrix of TCGA, we obtained the mRNA expression matrix related to melanoma autophagy. Furthermore, we set |Correlation coe cient| ≥5 and p-value <0.001 as the cut-off standards, and used the limma package of the R software to analyze the co-expression of the mRNA matrix and lncRNA matrix of melanoma. Finally, we screened 89 mRNAs and 100 lncRNAs in patients with melanoma (Supplementary Figure S1).

Survival analysis of lncRNAs associated with autophagy
Prior to the multivariate COX regression analysis, we performed a univariate Cox survival analysis on 100 lncRNAs. Finally, the Kaplan-Meier survival analysis showed that 39 lncRNAs were closely related to the prognosis of patients with melanoma. The high expression of only two lncRNAs was signi cantly correlated with shorter overall survival (OS) in patients with melanoma, whereas the remaining 37 lncRNAs did not show such a relationship.

Multivariate Cox regression model analysis
The combination of multiple genes can better predict the prognosis of patients. Thus, we conducted a stepwise COX regression analysis of 37 lncRNA with survival signi cance, and nally obtained a gene combination ( Figure 1, Table 1) composed of seven lncRNAs (HLA-DQB1 antisense RNA 1 [HLA-DQB1-AS1], USP30 antisense RNA 1 [USP30-AS1], AL645929, AL365361, long intergenic non-protein coding RNA 520 [LINC00520], LINC00324, and AC055822). In the multivariate COX regression analysis, patients were divided into high-and low-risk groups according to the median value of their riskScore. According to the riskScore, we used the survivalROC package to draw the survival curve ( Figure 2A) and ROC curve ( Figure 2B) of this gene combination. Survival results showed that the high expression of this gene combination was signi cantly associated with poor OS in patients at high risk. The AUC of 0.742 suggests that this gene combination may be used to predict the prognosis of patients with melanoma. The riskScore results showed that the high-risk group had a higher risk coe cient than the low-risk group, rising from left to right as illustrated in Figure 2C. Compared with the low-risk group, the survival diagram results showed that the high-risk group was characterized by more deaths and shorter survival time ( Figure 2D). Finally, the expression of seven genes in the high-and low-risk groups was visualized using a heat map ( Figure 2E).
In addition, to verify the higher accuracy of the risk model for predicting the prognosis of patients versus clinical information, we rst performed a univariate COX analysis based on clinical information, including sex, stage, T, M, and N. The forest plot results of the univariate analysis showed that stage, T, M, N, and riskScore could be used to predict the prognosis of patients ( Figure 3A). However, we used a multivariate COX regression analysis and ROC curve to draw the clinical information and determine the index that may be used as the best independent prognostic factor. The forest plot results of the multivariate COX analysis showed that T, N, and riskScore could be used as independent prognostic factors ( Figure 3B).
However, the ROC results showed that the AUC of riskScore, T, and N was 0.734, 0.707, and 0.658, respectively ( Figure 3C). Therefore, the risk model we constructed is more accurate than T and N in predicting patient survival, and is better than other traits ( Table 2).

Construction of the co-expression network
We constructed an autophagy-related co-expression network of seven lncRNAs and 16 mRNAs using the Cytoscape software ( Figure 4A, Table 3). HR >1 indicated a risk factor, whereas HR <1 indicated a protective factor. We found that, in the co-expression network, only LINC00520 was a risk factor, whereas the other six lincRNAs were protective factors. Of note, the Sankey diagram yielded consistent results ( Figure 4B).

GSEA in melanoma
We used GSEA to analyze the GO and KEGG, and study the potential pathogenic mechanism of patients at high and low risks of melanoma. The results of the GO enrichment analysis showed that the BP of the high-risk group was mainly related to the metabolism of melanin. The GO of patients at low risk was related to the positive regulation of autophagy and the activation of the immune response ( Figure 5A). The results of KEGG were related to the regulation of autophagy, antigen processing and presentation, and T-cell receptor signaling pathway ( Figure 5B). Therefore, we hypothesized that the potential pathogenesis of melanoma in patients at high and low risks is related to the regulation of autophagy and the immune system. It was further demonstrated that the immune response plays an important role in the autophagy of melanoma.

Discussion
Melanoma is a disease associated with high mortality and di cult to diagnose in the early stage [28]. Therefore, there is an urgent need for the identi cation of molecular markers for early diagnosis and treatment. At present, studies have found that autophagy is closely related to the occurrence of melanoma [29]. However, most studies tend to evaluate the prognostic and therapeutic role of ARGs in melanoma. Thus far, there are no investigations on the diagnostic and prognostic value of autophagy-related lncRNA in melanoma. Hence, it is urgent to clarify the molecular mechanism of autophagy-related lncRNAs in melanoma and identify new molecular therapeutic targets.
In this study, we constructed an autophagy-related co-expression network composed of seven lncRNA and 16 mRNA genes. Following GSEA, we found that the GO in the low-risk group was more enriched in the activation of the immune system and the positive regulation of autophagy, while the enrichment results of KEGG were similar. Therefore, we hypothesized that autophagy and immune system regulation may affect the transition from low risk to high risk in melanoma. Univariate Cox regression showed that, among the seven examined lncRNAs, only LINC00520 was a risk factor, and its high expression was signi cantly correlated with poor OS in patients at high risk of melanoma. The remaining six lncRNAs were protective factors. Further multivariate COX analysis showed that the gene set composed of the seven lncRNAs had high accuracy in predicting the prognosis of patients in the high-and low-risk groups.
According to the literature, the expression of mRNA and protein of autophagy genes decreased with the development of melanoma compared with normal melanocyte [30]. Beclin 1 (BECN1) is a gene that activates autophagy and programmatically inhibits tumor proliferation. Light chain 3 (LC3) is closely related to the formation of the autophagy membrane. Miracco et al. found that the expression of these two autophagy genes decreased with the progression of melanoma compared with a benign nevus [31].
The expression of these two genes decreased with the increase in tumor grade and lymph node metastasis. A study conducted by Ding et al. on invasive liver cancer yielded similar ndings [32]. Bioinformatics analysis performed by Chen et al. revealed that the biological processes of the clear cell renal cell carcinoma (ccRCC) low-risk group was related to the regulation of autophagy [33]. Li et al. also found that the positive regulation of autophagy and enrichment of immune activation pathways were more associated with patients in the low-risk group than those in the high-risk group [14]. Moreover, in our study, activation of the immune system also occurred in patients in the low-risk group. In addition, autophagy can promote the development of T and B cells in cancer and be used to present antigens to major histocompatibility complex class II molecules, which can be recognized by the immune system to play an anti-tumor role. A good prognosis in the low-risk group may be related to autophagy and activation of the immune system. Therefore, the reduction of autophagy may increase the risk of melanoma.
The WD repeat domain, phosphoinositide interacting 1 (WIPI1), is a member of the WIPI family [34]. It is involved in the formation of autophagosomes and the fusion of autophagy lysosomes. Studies have reported that WIPI1 was highly expressed in both prostate cancer and melanoma [35]. Arcangelo et al. pointed out that WIPI1 was highly expressed in melanoma cells compared with normal melanocytes, and tumor proliferation can be inhibited by inhibiting its expression [36]. In our study, WIPI1 and LINC00520 were co-expressed. As an oncogene, LINC00520 was highly expressed in numerous types of cancers. In breast cancer, LINC00520 promotes the proliferation of tumor cells by regulating STAT3 [37]. LINC00520 was also highly expressed in melanoma. Overexpression of LINC00520 can promote the proliferation, invasion, and migration of melanoma. At present, there is no report on LINC00520 and autophagy in any type of tumor [38]. Therefore, further experimental veri cation is warranted. In our autophagy co-expression network, LINC00324 was co-expressed with C-C motif chemokine receptor 2 (CCR2), CASP8 and FADD like apoptosis regulator (CFLAR), and protein kinase C theta (PRKCQ). CCR2 is a speci c chemokine for monocyte chemotaxis, playing a role in resisting many types of tumors in the human immune system [39]. According to the literature, CCR2 can enhance the monitoring of the immune system, and play a role in tumor resistance by triggering the TH1 response and recruiting CD8+ cells [40]. CFLAR is an anti-apoptotic protein, and its function is mostly related to the autophagy of T cells. In the body, a decrease in this autophagy gene can affect the development of peripheral T cells and decrease the proliferation ability of peripheral T cells [41]. Dunkle et al. found that the rate of apoptosis of peripheral T cells with CFLAR de ciency increased [42]. Following bioinformatics analysis of ARGs in ccRCC, Chen et al. found that PRKCQ was associated with poor OS in patients with ccRCC, and may be used to evaluate their prognosis [33]. In summary, these genes are often associated with immune response and autophagy in tumors. According to our study, co-expression of CCR2, CFLAR, PRKCQ, and LINC00324 was associated with autophagy in melanoma. In our risk model, LINC00324 was a protective factor in patients at low risk. This observation may be related to activation of the immune response and positive regulation of autophagy in patients at low risk. In addition, in our study, HLA-DQB1-AS1, USP30-AS1, AL645929, AL365361, and AC055822 were also protective factors in patients at low risk. Moreover, in the coexpression network, they were closely related to CCR2, CFLAR, and PRKCQ. Surprisingly, except for the high expression of LINC00520 in melanoma, thus far, we have not observed the expression and mechanism of the other six lncRNAs in melanoma. Moreover, we did not nd any data reporting the role of these seven lncRNAs in the autophagy mechanism of melanoma. Therefore, our study provides a new idea for revealing the autophagy mechanism of this disease.
Autophagy was closely related to the prognosis of patients with melanoma. Hence, it is crucial to identify accurate prognostic markers for patients with melanoma. By constructing the co-expression network of ARGs and lncRNAs in melanoma, we performed univariate and multivariate COX regression analyses of lncRNAs, and eventually included seven lncRNAs in the risk model. To verify the speci city and accuracy of our risk model, we also performed univariate and multivariate COX analyses using the clinical information of patients with melanoma. The nal results showed that, compared with other clinical information, our risk model can more accurately predict the prognosis of patients with melanoma and provide new ideas for the mechanism and treatment of this disease.

Conclusion
we constructed an autophagy-related co-expression network using autophagy-related mRNAs and lncRNAs. This will provide a basis for our subsequent research on the autophagy mechanism of melanoma. In addition, we constructed a risk model, which may provide new insight into the diagnosis and prognosis of patients with melanoma in the future versus clinical information. Further investigations are warranted to verify the present ndings.