Autophagy-Related Long Non-coding RNA Signature as Indicators for the Prognosis of Uveal Melanoma

This study aimed to develop an autophagy-associated long non-coding RNA (lncRNA) signature to predict the prognostic outcomes of uveal melanoma (UM). The data of UM from The Cancer Genome Atlas (TCGA) were enrolled to obtain differentially expressed genes (DEGs) between metastasizing and non-metastasizing UM patients. A total of 13 differentially expressed autophagy genes were identified and validated in Gene Expression Omnibus, and 11 autophagy-related lncRNAs were found to be associated with overall survival. Through performing least absolute shrinkage and selection operator regression analyses, a six-autophagy-related lncRNA signature was built, and its efficacy was confirmed by receiver-operating characteristic, Kaplan–Meier analysis, and univariate and multivariate Cox regression analyses. A comprehensive nomogram was established and its clinical net benefit was validated by decision curve analysis. GSEA revealed that several biological processes and signaling pathways including Toll-like receptor signaling pathway, natural killer cell-mediated cytotoxicity, and B- and T-cell receptor signaling pathway were enriched in the high-risk group. CIBERSORT results showed that the signature was related to the immune response especially HLA expression. This signature could be deployed to assist clinicians to identify high-risk UM patients and help scientists to explore the molecular mechanism of autophagy-related lncRNAs in UM pathogenesis.


INTRODUCTION
Uveal melanoma (UM) is the most common and aggressive malignant neoplasm of the eye in adults, and up to 50-62% of UM patients will suffer from metastasis (Kujala et al., 2003;Amaro et al., 2017;Kaliki and Shields, 2017). Local interventions including transpupillary thermotherapy (Mashayekhi et al., 2015), photodynamic thermotherapy (Cerman and Cekic, 2015), brachytherapy (Collaborative Ocular Melanoma Study Group, 2006), and surgical enucleation  are the main strategies to treat the primary tumor, but the overall prognosis for metastatic UM remains poor. Therefore, it is imperative to identify high-risk UM patients with poor prognosis. Clinical features including patient age, tumor apical height, involvement of ciliary body, and largest basal diameter (LBD) of the tumor have been identified as prognostic factors for metastasis (Chew et al., 2015;Correa and Augsburger, 2016). Recently, scientists have proposed certain molecular signatures as predictors of high-risk UM, such as the gene expression profile (GEP) model (Zhang et al., 2014), DNA methylation-based model , microRNA (miRNA) model (Xin et al., 2019), and alternative splicing event model (Wan et al., 2020). However, the prognostic value of long non-coding RNAs (lncRNAs) in UM has not been fully explored.
LncRNAs are a set of RNAs with a transcript length of over 200 nucleotides and have limited or no protein-coding capacity. Previous data has shown that abnormal activation of lncRNAs participates in various biological processes and human tumors including UM. A number of studies have shown that aberrant lncRNA expressions like MALAT1, FTH1P3, ZNNT1, and CANT1 played important roles in the initiation and progression of UM (Xing et al., 2017;Zheng et al., 2017;Li et al., 2020;Wu et al., 2020). Recent evidence highlights the involvement of lncRNA ZNNT1 in the tumorigenesis of UM by regulating key autophagy gene expression (Li et al., 2020). Therefore, there exist some interactions between lncRNAs and autophagy processes that could contribute to the pathogenesis of UM.
Autophagy is a highly conserved intracellular degradation system required for the maintenance of homeostasis and promotion of survival when facing cellular stress. It plays critical roles in age-related diseases including cancers (Mizushima, 2007). According to previous studies, autophagy was proved to be upregulated in UM and to be associated with poor prognosis (Giatromanolaki et al., 2011). Specific autophagy inhibitors might serve as a potential approach to suppress UM (Zhang et al., 2016). On the other hand, increased autophagy could be protective and cooperate with apoptosis to induce cell death in UM (Ambrosini et al., 2013). Hence, exploring autophagy-related gene signatures may deepen the understanding of the underlying molecular mechanism of UM and provide insights into its prognosis.
Based on the above background, we postulated that autophagy-related lncRNAs might help us identify highrisk UM patients and customize their personalized targeted treatment. This study aimed to analyze the correlation between the differentiated expressed lncRNAs and autophagy-related genes in UM and to screen for the most survival-relevant autophagy-related lncRNA signature by least absolute shrinkage and selection operator (LASSO) regression. Furthermore, univariate and multivariate Cox regression analyses were used to confirm its prognostic value in UM. This study may provide new insights into prognostic markers and deepen the understanding of autophagy-related lncRNA in UM.

Data Resource and Preprocessing
The RNA expression profiles of 80 UM patients were downloaded in FPKM format from The Cancer Genome Atlas (TCGA) database using the GDC Data Transfer Tool 1 , and the mRNA and lncRNA expression data were also acquired. The relevant clinicopathological characteristics including sex, age, and cancer stage were downloaded from the UCSC Xena website 2 . The detailed clinical characteristics are shown in Table 1. The flowchart of the analysis procedure is shown in Figure 1. The gene expression data of GSE22138 that contains 63 tissue samples of UM patients were downloaded from the Gene Expression Omnibus (GEO) as a validation dataset (Laurent et al., 2011).

Identification of Differentially Expressed Genes and Autophagy-Related LncRNAs
In order to analyze the metastasis-related gene characteristics of UM patients, the differentially expressed genes (DEGs) between the metastasizing and non-metastasizing groups were calculated using the limma R package (Ritchie et al., 2015). DEGs with absolute log (fold change) > 0.05 and false discovery rate (FDR) < 0.05 were considered to be included for subsequent analysis. The result of differential analysis was presented in a heat map. Then, the autophagy-related gene set was downloaded from the MSigDB database 3 (Subramanian et al., 2005), and the metastasis-related autophagy DEGs were extracted by the VennDiagram R package (Chen and Boutros, 2011). The expression correlations between those DEGs were calculated, and the protein-protein interaction (PPI) network was constructed by the Search Tool for the Retrieval of Interacting Genes (STRING) database 4 . The autophagy-related lncRNAs were identified based on the criteria that the correlation coefficient | R 2 | > 0.3 and P < 0.05.

Construction of Prognostic Autophagy-Related LncRNA Signature
The LASSO method was used to further analyze the prognostic autophagy-related lncRNAs for dimensionality reduction and feature selection. The coefficients obtained from the LASSO regression were weighted by each normalized gene expression value to construct the prognostic signature model, and the risk score of each UM patient was calculated by the formula as follows: Based on the median risk score, UM patients in the TCGA cohort were divided into high-risk and low-risk subgroups. Then, univariate and multivariate Cox regression analyses were performed to analyze the predictive ability of the signature model and patients' clinical-pathological parameters on overall survival (OS) and progression-free interval (PFI). A prognostic nomogram integrating risk scores and clinical-pathological characteristics was constructed to analyze the probable 1-, 3-, and 5-year OS of the UM patients. Harrell's concordance index (C-index) was applied to evaluate the predictive accuracy of the nomogram. To determine the clinical predictive validity of the nomogram, decision curve analysis (DCA) was used to quantify the net benefits at threshold probabilities in the dataset.

Comparison of Immune Cell Infiltration Levels and Immune Scores Between the High-Risk and Low-Risk Subgroups of UM Patients
To quantify the proportion of different immune cells in the UM samples, the CIBERSORT algorithm and the LM22 gene set were used to analyze 22 human immune cell phenotypes in the tumor immune microenvironment (TME) (Newman et al., 2015). CIBERSORT is an algorithm for calculating the cell-type proportions in tumor samples, which uses a set of reference gene expression values as a representation (consisting of 547 marker   Frontiers in Genetics | www.frontiersin.org genes). The Mann-Whitney U test was used to compare the infiltration level of immune cells between the two groups.
ESTIMATE (Estimation of STromal and Immune cells in MAlignant Tumor tissues using Expression data) is an algorithm for quantifying the immune infiltration level in a tumor sample based on gene expression data (Yoshihara et al., 2013). We calculated the immune and stromal scores of each UM sample using the ESTIMATE R package and compared the expression differences of human leukocyte antigen (HLA) family genes between the two groups.

Survival Analysis
The effects of the signature model on OS and PFI were compared. The Kaplan-Meier survival curve by survival package was used to perform survival analysis, and the log-rank test was used to assess the significance of the difference in survival time between the two groups.

Statistical Analysis
All data processing and analysis were performed using the R software (version 3.6.2). For the comparison of the two sets of continuous variables, the statistical significance of normally distributed variables was estimated by independent Student's t-test, and differences between non-normally distributed variables were analyzed by Mann-Whitney U test. Categorical variables in the two groups were compared and analyzed using the chi-squared test or Fisher's exact test. Pearson correlation analysis was used to calculate correlation coefficients between different genes. Univariate and multivariate Cox regression analyses were performed to identify independent prognostic factors. The pROC R package was used to visualize receiveroperating characteristic (ROC) curves, and the area under the curve (AUC) was calculated to assess the prognostic accuracy of the risk score (Robin et al., 2011). P < 0.05 was considered to be statistically significant.

Identification of Metastasis-Associated Autophagy DEGs
To obtain metastasis-associated signature genes in UM patients, the gene expression profiles of metastasizing and non-metastasizing UM patients were analyzed and 1,107 DEGs were identified [| log (fold change)| > 0.5 and FDR < 0.05 filters, Figure 2A]. Meanwhile, the Venn diagram showed that there were 13 autophagy DEGs (Figure 2B). GO and KEGG analyses of all DEGs are shown in Supplementary Figure 1. GO and KEGG analyses of autophagy DEGs revealed that these DEGs were mainly enriched in cellular responses to starvation and extracellular stimulus and in relevant autophagy processes (Figures 2C,D).
We further analyzed the interrelationship between metastasisassociated autophagy DEGs. As shown in Figures 3A,B, all 13 autophagy DEGs were differentially expressed with good concordance both in the TCGA and GEO cohorts. The correlation analysis showed significant correlations between the expression of most of these genes, with the highest correlation coefficient between RAF1 and PIK3R4 ( Figure 3C). Meanwhile, the PPI network further revealed the interactions between 13 autophagy DEGs ( Figure 3D).

Identification of Prognostic Autophagy-Related LncRNAs in UM
As shown in Figure 4A and Supplementary Figure 2, 27 autophagy-related lncRNAs were identified. The univariate Cox regression analysis showed that 14 lncRNAs had statistical significance, among which five lncRNAs were considered as protective factors and nine lncRNAs were considered as risk factors (Figure 4B). In order to confirm the impacts of these lncRNAs on the UM prognosis, all UM patients were ranked according to each lncRNA expression and subsequently divided into low and high groups (each group n = 40) according to media value. The survival prognostic analysis showed that 11 lncRNAs had statistically significant impacts on the prognosis of UM patient ( Figure 4C).

Construction of Prognostic Autophagy-Related LncRNA Signature
To quantify the impact of autophagy-related lncRNAs on the prognosis of each UM patient, we constructed a prognostic signature by LASSO regression analysis. As shown in Figure 5A, six hub autophagy-related lncRNAs were included in the prognostic signature model. The risk score for each UM patient was calculated using the coefficients of each gene obtained from the LASSO algorithm. The distribution and scatter plot of risk scores, survival curves, and gene expression plots for each patient are shown in Figure 5B.
Kaplan-Meier analysis showed that the OS and PFI of patients with high-risk scores were significantly shorter than those with low-risk scores (Figure 5C). The ROC results validated the prognostic accuracy of the signature model ( Figure 5D). The area under the ROC curves (AUCs) of 1, 3, and 5 years were 0.919, 0.845, and 0.724, which confirmed the accuracy of the signature in predicting survival times of UM patients. The principal component analysis (PCA) showed that the signature model could separate the distribution patterns between high-and low-risk groups in terms of the six autophagy-related lncRNAs and all genes ( Figure 5E).

Clinical Correlation Analysis of the Signature and Nomogram Building
The clinical correlation analysis showed that high-risk patients demonstrated a higher risk of metastasis and poorer stage compared with low-risk patients, whereas no significant correlation with patient's age was found (Figures 6A-C). The univariate and multivariate Cox regression analyses showed that the lncRNA signature was an independent prognostic indicator for UM patients (Table 2 and Figures 6D,E). Then we constructed a nomogram that integrates clinicopathological characteristics and the signature to predict OS of UM patient ( Figure 6F). C-index was used to calculate the discriminative power of the nomogram, which showed a high degree of discrimination [0.93 (0.975-0.885)]. The DCA showed a prognostic benefit in 18-95% of patients using the signature model ( Figure 6G).

Gene Set Enrichment Analysis and Immune Infiltration Analysis of the Signature Model
GSEA results showed that the proteasome, protein transport, metabolism, and apoptosis pathways were significantly enriched in the high-risk group compared with those in the low-risk group (Table 3 and Figure 7A). Meanwhile, some biologically relevant pathways, such as EMT, DNA replication, and cell cycle regulatory factors (Figure 7B), were significantly different between the high-and low-risk groups (P < 0.05).
GSEA results also showed that the change in DEGs was significantly enriched in immune-related signaling pathways, including Toll-like signaling pathway, natural killer cell signaling pathway, etc ( Figure 8A). Therefore, we further analyzed the immune infiltration between the high-and low-risk groups. As shown in Figures 8B,C, the high-risk patients had significantly higher immune scores and stromal scores compared with the low-risk patients. To further analyze the infiltration of immune cells in UM patients, the gene expression data of UM patients were downloaded from the TCGA cohort and 22 different specific immune cell infiltration in each patient were obtained using the CIBERSORT algorithm ( Figure 8D). As shown in Figure 8E, most of the immune cells were significantly different between the two groups ( Figure 8E). Furthermore, the gene expression of the HLA family was significantly different between the two groups ( Figure 8F).

DISCUSSION
The mortality rate of UM patients with advanced stage is higher than 95% within 5 years, and the survival of patients with metastasis UM is below 1 year (Kujala et al., 2003), so finding novel and effective prognostic biomarkers is critical for UM. Although some epigenetic signatures like miRNA signature and DNA methylation signature have been established to precisely identify high-risk UM, the lncRNA signature for UM has not been reported. Experimental studies have shown that lncRNAs play a critical role in UM by regulating many cellular processes including autophagy. Therefore, it may be helpful to build an autophagy-related lncRNA signature to predict the prognosis of UM and improve therapeutic approaches.   In the present study, two datasets (TGGA and GSE22138) were enrolled to construct novel prognostic autophagy-related lncRNA signatures for UM. After comparing DEGs between metastasizing and non-metastasizing UM patients, we screened out 13 autophagy genes and further identified 27 autophagyrelated lncRNAs through establishing the lncRNA-autophagy gene co-expression network. Moreover, we utilized the LASSO regression mode to build a six-hub autophagy-related lncRNA model that could differentiate high-and low-risk UM patients. The survival analyses, ROC curve, and PCA analyses confirmed the accuracy of this model. Through univariate and multivariate Cox regression analyses, we also revealed that the risk score of this model was an independent factor in predicting prognosis of UM patients. This six-lncRNA model may guide in the selection of more precise and individualized treatment strategies for UM patients.
Through bioinformatics analysis, we identified 13 metastasisrelated autophagy DEGs, namely CDKN1A, DLC1, FOXO3, GABARAPL1, GRID1, ITGA6, MTMR14, PIK3R4, PRKCD, RAF1, SPHK1, TUSC1, and WIPI1. Among them, FOXO3, which acts as a transcription factor to correct autophagy perturbations, was reported to attenuate proliferation, invasion, and migration and promote cellular apoptosis of UM (Yan et al., 2016(Yan et al., , 2018. PRKCD was shown to suppress autophagy and implicated as a link to activate MAPK in GNAQ mutant UM (Chen et al., 2017;Zhang et al., 2017). Meanwhile, we found that RAF1 and PIK3R4 showed the highest correlation coefficient. RAF1 acts as a regulatory link between the membrane-associated Ras GTPases and the MAPK/ERK pathway, and this critical regulatory link can lead to inhibition of cell proliferation and induction of apoptosis in UM cells (Samadi et al., 2012). PIK3R4 belongs to the regulatory subunit of the PI3K complex that is involved in the maturation of autophagosomes and endocytosis. A previous study revealed that somatic mutation of PIK3R4 contributed to melanoma metastasis (Shull et al., 2012), and its functions in UM progress are worthy of further study. It is worth mentioning that only 13 DEGs out of 1,107 genes in our results were found to be autophagyrelated genes, and only 11 autophagy-related lncRNAs were associated with overall survival for UM. However, the autophagyrelated lncRNA signature we built was proven to be highly effective for predicting the overall survival for UM. These results may be explained by the fact that autophagy takes part in a certain but crucial phase of the metastatic cascade, but other complicated mechanisms take up more roles in regulating UM metastasis.
Autophagy is a lysosomal catabolic process that is regulated by non-coding RNAs such as miRNAs and lncRNAs. Experimental and clinical studies revealed that autophagy-related lncRNAs had significant diagnostic and therapeutic values in various cancers including UM. For example, a recent research reported that lncRNA ZNNT1 could act as a tumor inhibitor in UM by upregulating ATG12 gene expression, indicating delicate interactions between lncRNAs and autophagy (Li et al., 2020). Our model enrolled six autophagy-related lncRNAs including PVT1 that predicted overall survival in UM. PVT1 has been identified as a prognostic biomarker and played a critical role in the development of pancreatic ductal adenocarcinoma via triggering cytoprotective autophagy (Huang et al., 2018), and it affected EMT and cell proliferation and migration via (F) Comparisons of various HLA family genes between the high-risk and low-risk groups. *P < 0.05, **P < 0.01, ***P < 0.001. modulating CDKN1A expression in breast cancer . Furthermore, PVT1 also took part in drug resistance of pancreatic cancer by interacting with ATG14 and promoting PtdIns3K activity (Zhou et al., 2020). In UM, silencing PVT1 could suppress tumorigenesis by activating the miR-17-3pdependent p53 signaling pathway (Wu et al., 2019). Further indepth studies are needed to investigate the molecular mechanism between PVT1 and autophagy in UM.
In order to explore the biological functions of the prognostic six-lncRNA model in UM, GSEA analysis was performed to compare the low-and high-risk groups. The results demonstrated that these lncRNAs might exert influence on UM progression through the immune-related pathway. LncRNAs are involved in immunological processes in cancers by mediating the functional state of immune cells and related pathways. Based on this lncRNA model, we found that the immune scores and stromal scores in high-risk UM were higher than in low-risk UM, and various members of the HLA family showed statistical difference between the two groups of UM. Several studies have shown an abnormal expression of the HLA family in UM and confirmed that a high expression of HLA class I is a risk factor for the development of metastasis (Gezgin et al., 2017;Souri et al., 2019;Wierenga et al., 2019). Recent studies showed that the lncRNA MIR155 host gene (MIR155HG) can be a potential prognostic marker for predicting the curative effect of immune checkpoint inhibitor therapy in multiple cancers including UM (Peng et al., 2019). Our results, combined with these publications, indicated that the six autophagy-related lncRNAs might play crucial roles in the pathogenesis of UM by regulating related genes to affect the immune responses, and there existed close relationships between autophagy and immune regulation in UM.
This study was not without limitations. All results of this study are based on bioinformatics analysis, and the prognostic signature enrolled some autophagy-related lncRNAs that no research confirmed their roles in UM so far. Thus, further experiments combining biochemical experiments and clinical prognostic data are required to confirm how these lncRNAs affect the prognosis of UM patients through autophagy exactly.

CONCLUSION
In conclusion, this study screened and identified 13 dysregulated autophagy-related genes in UM, and a prognostic risk signature enrolling the six dysregulated autophagy-related lncRNAs was constructed, which was proven to be a useful marker for predicting UM patients' prognosis with high efficacy. Based on the prognostic lncRNA model, dysregulated immune-related pathway, especially HLA expression, may contribute to the poor prognosis in high-risk UM. This signature model could be deployed to assist clinicians to identify high-risk UM patients and help scientists to explore the molecular mechanism of autophagyrelated lncRNA in UM pathogenesis.

DATA AVAILABILITY STATEMENT
All datasets generated for this study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author/s.