Identification of potential autophagy-associated lncRNA in prostate cancer

Background: Long non-coding RNAs (lncRNAs) have been linked to autophagy. It is urgent to identify and assess the hub autophagy-associated lncRNA in prostate cancer. Methods: Differentially expressed lncRNAs associated with autophagy were identified in prostate cancer based on The Cancer Genome Atlas Prostate Adenocarcinoma (TCGA-PRAD) data. An autophagy-mediated competing endogenous RNA network was constructed to screen for autophagy-associated lncRNA, and the preselected lncRNAs were further validated using Gene Expression Omnibus (GEO) datasets. Furthermore, a prognostic lncRNA signature was established and assessed. Additionally, Gene Set Enrichment Analysis (GSEA) revealed the underlying molecular mechanisms. Results: Using a competing endogenous RNA network, 66 differentially expressed lncRNAs associated with autophagy were identified, and the differential expression of 7 lncRNAs were verified using the TCGA-PRAD, GSE21034, and GSE94767 datasets. Additionally, a lncRNA signature associated with autophagy, including MKNK1-AS1 and INE1, was identified as an independent indicator of survival with a C-index of 0.882. The GSEA analysis indicated that several autophagy-related signaling pathways were enriched in different risk groups. Conclusions: The lncRNAs associated with autophagy were identified, and a prediction model was developed that could be used as a prognostic predictor for prostate cancer, indicating the critical role of lncRNA in the regulation of prostate cancer autophagy regulation.


INTRODUCTION
AGING of cellular product [5]. Autophagy is disrupted in several diseases, including neurodegenerative diseases, inflammatory diseases, cardiovascular diseases, and cancer [6]. Autophagy has been implicated in cancer in numerous studies, and the effective targeting of autophagy may be an underlying therapeutic strategy in advanced cancer [7]. Death-associated protein kinase 3 (DAPK3), a novel autophagy regulator, may phosphorylate ULK1 and enable it to suppress gastric cancer progression [8]. CircMUC16 induced autophagy and accelerated epithelial ovarian cancer progression by interacting directly with miR-199a and ATG13 [9]. In addition, a growing number of studies have demonstrated that impaired autophagy is linked to prostate cancer [10]. Phospholipase C epsilon (PLCɛ) modulates the AR signaling activities in PCa cell lines through the degradation mechanism that is dependent on autophagy [11]. MicroRNA-381 reduces the RELN-mediated PI3K/AKT pathway activity, which inhibits cell proliferation and promotes apoptosis and autophagy in PCa cells [12]. These findings support the correlation between autophagy and PCa. However, the potential molecular mechanisms of action associated with autophagy in PCa are still unknown and unexplored.
Numerous studies have examined the critical functions of lncRNA in the initiation and progression of cancer [13]. Additionally, some lncRNAs are involved in the occurrence and progression of PCa [14]. For example, dysregulations of KCNQ1OT1 have been discovered in PCa cells and tissues. Enhanced lncRNA KCNQ1OT1 expression sponged miR-15a to facilitate immune escape and progression by improving the expression of PD-L1 in PCa [15]. LINC00173 functioned as an oncogene in PCa via the LINC00173/miR-338-3p/Rab25 pathways [16]. These findings revealed a strong connection between lncRNA and PCa. However, the study of the lncRNA in PCa is insufficient, and it is urgent to discover lncRNAs that may play a vital role in PCa.
In this study, we aim to establish a network of ceRNA to screen for lncRNA associated with autophagy as novel biomarkers in PCa and to construct a lncRNA signature associated with autophagy to predict the prognosis of PCa patients. It is valuable to provide new information about the clinical application of the lncRNA associated with autophagy in PCa.

Verification of the identified lncRNAs in the GEO database
We further selected the GSE21034 and GSE94767 datasets from the GEO database to verify the differential expression of 66 autophagy-linked lncRNAs filtered by the autophagy-mediated ceRNA network. As illustrated by Figure 3A, 3B, only 7 (ADAMTS9-AS1, ADAMTS9-AS2, MAGI2-AS3, PCA3, PCAT1, PVT1, SNHG3) of the 66 autophagy-related lncRNAs in the TCGA, GSE21034, and GSE94767 datasets showed abnormal expression levels in PCa tissues compared to normal prostate tissue.
Identification of the autophagy-associated lncRNA signature 232 autophagy-related genes were acquired from HADb and 1357 autophagy-related lncRNAs were identified using the Pearson correlation coefficients between the expression of lncRNAs and autophagy-related genes with |R| > 0.5 and P < 0.05 as screening criteria. Additionally, univariate Cox regression and Kaplan-Meier analysis were used to obtain AC008760.1, AC134775.1, MKNK1-AS1, and INE1. LASSO regression was further used to identify three autophagyassociated lncRNAs (AC008760.1, MKNK1-AS1, and INE1) in Figure 4. Subsequently, multivariate Cox AGING regression revealed that two of the three autophagyassociated lncRNAs, MKNK1-AS1 and INE1, were candidates for the prognostic signature derived from the lowest Akaike information criterion (AIC) value (Table  1). MKNK1-AS1 and INE1, two autophagy-associated lncRNAs in the prognostic signature, were both found to be adverse prognostic factors with a hazard ratio (HR)>1. In the lncRNA signature associated with autophagy, the following risk score formula for patients with PCa was as follows: risk score = 0.1277 * expression level of MKNK1-AS1 +0.0186 * expression level of INE1.

Assessment of the autophagy-associated lncRNA signature
PCa patients were divided into low-risk and high-risk groups for autophagy-associated lncRNA signature analysis based on the median risk score. Furthermore, the survival curve indicated that patients with high-risk scores had a shorter survival time ( Figure 5A). As shown in Figure 5B, the time-dependent ROC curve analysis exhibited that the AUC value of the ROC curve for the autophagy-associated lncRNA signature corresponding to three years and five years of survival AGING was 0.861 and 0.799, respectively. The results indicated that the autophagy-related lncRNA signature might be used to predict prognostic value. Additionally, the Cindex of the prognostic signature was 0.882 (95% CI: 0.800-1.055), indicating that the risk score had a reliable and accurate diagnostic impact on PCa diagnosis. As shown in Figure 5C, the risk score distributions for PCa patients were categorized and  AGING plotted based on the risk score. Besides, the scatter plot exhibited the relationships between risk score, survival status, and survival time in PCa patients ( Figure 5D). Moreover, the heatmap revealed two autophagyassociated lncRNAs were expressed differently in the high-risk and low-risk PCa patients ( Figure 5E). These findings suggested that the autophagy-associated lncRNA signature appeared to be useful in the diagnosis and prognosis of PCa.

Independent predictive factor in PCa
Univariate Cox regression exhibited that the T stage (P =0.049) and risk score (P < 0.001) were significantly associated with survival in the PCa, with the HR of risk score being 1.173 (95% CI: 1.079-1.275) in Table  2. Furthermore, multivariate Cox regression further revealed that the risk score (P =0.010) was significantly associated with PCa survival, with a HR of 1.145 (95% CI: 1.033-1.270) ( Table 2). Additionally, the ROC curve analysis revealed that the AUC value of the autophagy-associated lncRNA signature was 0.903, which was higher than the AUC values for age (AUC = 0.587), T stage (AUC= 0.628), N stage (AUC=0.572), Gleason score (AUC=0.668), PSA value (AUC= 0.680), biochemical recurrence (AUC=0.638) and Race (AUC=0.428) in Figure 6. In addition, we further analyzed the relationship between autophagy-associated lncRNA signature and clinical features (Table 3), and the results revealed that the Gleason score was significantly related to the autophagy-associated lncRNA signature (P =0.008). Overall, these results indicated that the autophagyrelated lncRNA signature was an independent prognostic factor for PCa patients.

Gene set enrichment analysis
To identify the potential enhanced signaling pathway between the high-risk and low-risk groups, we performed the Gene Set Enrichment Analysis (GSEA). As shown in Figure 7, peroxisome (NES=--2.191, P=0.002), lysosome (NES=--2.110, P<0.001), RNA polymerase (NES=-2.062, P<0.001), proteasome (NES=-1.962, P<0.001), and protein export (NES=-1.856, P=0.002) were all significantly altered in the low-risk group, indicating that autophagy was strongly related to PCa. These results provided a new horizon for the development of personalized treatments for PCa patients with different risk scores.

DISCUSSION
The occurrence rate of prostate cancer, an epithelial malignancy that occurs in the prostate gland, has remained elevated in recent years [17]. Prostate cancer is one of the malignant tumors with the highest fiveyear survival rate. Early patients encounter challenges  in diagnosis, and advanced patients who are predisposed to relapse have little medical choices [18]. Besides, autophagy is involved in the malignant progression of a variety of tumors, including PCa [19]. As a result, autophagy-related biomarkers may be used as diagnostic and therapeutic targets for PCa patients. Autophagy-related biomarkers have been implicated in the diagnosis and treatment of PCa [20].
LncRNAs are engaged in numerous biological processes, including epigenetic modulation, imprinting, and transcription [21]. Some lncRNAs can affect posttranscriptional regulation by interfering with microRNA. Besides, some lncRNAs can control gene expression at the epigenetic, transcriptional, and posttranscriptional levels [22][23][24][25]. Moreover, several lncRNAs have been shown to be effective biomarkers for PCa risk prediction with high sensitivity and specificity [26]. For example, a four-lncRNA signature was developed and used to predict biochemical recurrence-free survival in PCa patients, which will assist in the identification of high-risk patients who  require more aggressive therapeutic interventions [27]. Importantly, lncRNAs can modulate autophagy in the malignant tumor by regulating targeted gene expression, and lncRNAs that modulate autophagy are vital in signal transduction and the underlying molecular mechanism in cancer [28]. Therefore, autophagy-related lncRNAs may be used as novel diagnostic and therapeutic markers for the PCa. In the current study, we focused on screening for lncRNAs that might play a pivotal role in the PCa by establishing a ceRNA regulatory network and constructing a signature of autophagy-related lncRNA for the diagnosis and prognosis of PCa patients.

AGING
The ceRNAs, acting as miRNA sponges, compete for miRNA binding sites and competitively regulate mRNAs using shared miRNA response elements in ceRNA networks. Moreover, the ceRNAs are all transcripts that may become miRNA targets, including non-coding long RNA (lncRNA), pseudogenic RNA, Figure 6. Independent predictive factor in prostate cancer patients. Receiver operating characteristic (ROC) curve analysis shows the prognostic accuracy of clinicopathological parameters such as age, T stage, N stage, Gleason score, PSA value, biochemical recurrence, Race, and autophagy-related lncRNA prognostic risk score.  AGING and circular RNA. The competing endogenous RNAs (ceRNA) hypothesis indicated that ceRNAs could indirectly modulate gene expression through competitive miRNA binding [29,30]. However, studies on lncRNA-mediated ceRNA regulatory networks in PCa are insufficient. In this study, the analysis of DEmRNA, DEmiRNA, and DElncRNA between PCa and paracancerous tissues used the TCGA-PRAD data. Moreover, we have established a PCa-specific ceRNA regulatory network in PCa by integrating the interaction between autophagy-related DEmRNAs and DEmiRNAs or DEmiRNAs and DElncRNAs. Additionally, GEO data confirmed that the seven lncRNAs (ADAMTS9-AS1, ADAMTS9-AS2, MAGI2-AS3, PCA3, PCAT1, PVT1, and SNHG3) were expressed differentially in the ceRNA network, and these seven lncRNAs are related to autophagy in PCa. PCA3, as one of the selected autophagy-related DElncRNAs, expressed highly in PCa cells and tissues and has been identified and applied in clinical practice as a molecular marker of prostate cancer [31]. PCAT1 regulated the expression of CENPF, ID1, and ID3 in the cell cycle and proliferation, and promoted tumorigenesis in PCa by modulating FSCN1 via miR-145-5p [32,33]. ADAMTS9-AS1 functioned as a ceRNA, a sponge for the hsa-mir-96, and a supportive regulator of PRDM16 expression in PCa [34]. MAGI2-AS3 was used in a nine-RNA signature for PCa prognosis through coexpression network analysis [35]. Then, the abnormally expressed PVT1 functioned as an oncogene in PCa, contributing to tumor growth [36]. Additionally, SNHG3 interactd with miR-577 to enhance SMURF1 expression in PCa [37]. Besides, ADAMTS9-AS2 is a new marker with no previous studies in PCa. Most of the lncRNAs examined in this study have only a small amount of research, and this study has laid a strong foundation for further exploring the function of lncRNA in autophagy PCa.

AGING
The majority of previous studies have concentrated on the function of genes involved in autophagy [38,39]. There are currently no research that use autophagyassociated lncRNA signatures to predict PCa patient survival outcomes. Hence, it is necessary to establish autophagy-related lncRNA signatures to predict the prognosis of PCa patients. In this study, we obtained an autophagy-related lncRNA signature via Lasso regression and Cox regression, which is useful in the diagnosis and prognosis of PCa. Besides, the autophagy-associated lncRNA signature is an independent predictor of PCa, suggesting that the autophagy-associated lncRNA signature may be a prognostic biomarker for PCa patients. Additionally, there is currently a lack of research on the two autophagy-related lncRNAs (INE1 and MKNK1-AS1) in PCa. More research is therefore needed to explore how these lncRNAs affect the prognosis of PCa patients through autophagy.
This research contains certain drawbacks. Firstly, our data is primarily derived from publicly available databases, and the clinical implementation of our results necessitates further investigation. Secondly, we have not further experimentally verified the function of the screened lncRNAs. Finally, the signature requires external validation.
In conclusion, we developed a ceRNA regulatory network to investigate the lncRNAs that could play a key role in PCa autophagy, and we tested the differential expression of the selected lncRNAs using multiple GEO datasets. In addition, we successfully constructed an autophagy-related lncRNA signature for the diagnosis and prognosis of PCa patients, which could help with clinical decisions.

Patient data extraction
Data on RNA-sequencing (RNA-seq) and corresponding clinical information were derived from The Cancer Genome Atlas Prostate Adenocarcinoma (TCGA-PRAD, https://portal.gdc.cancer.gov/) [40]. The criteria for excluding cases were developed as follows: (1) patient follow-up lasted less than 30 days; (2) patients received medication before resection; (3) the histological type of the patients did not belong to prostate adenocarcinoma acinar type; (4) Patient clinical information or RNA-seq data is not available.

Establishment of a ceRNA regulatory network
The interaction between autophagy-associated DEmRNAs and DEmiRNAs was predicted using the miRTarBase [41], miRDB [42], and the TargetScan database [43]. In addition, the interaction between DEmiRNAs and DElncRNAs was matched by the miRcode database [44]. A ceRNA regulatory network was constructed to visualize the interaction between the selected miRNAs and mRNAs or lncRNAs AGING associated with autophagy by Cytoscape (version 3.7.2) [45].

Selection of autophagy-associated lncRNAs
Two hundred and thirty-two genes linked to autophagy were derived from HADb (Human Autophagy Database, http://www.autophagy.lu/clustering/) [47]. Autophagy-related lncRNAs were also determined based on Pearson correlation coefficients between the expression of lncRNAs and autophagy-related genes. The absolute value of the correlation coefficient > 0.5 and the P-value < 0.05 were applied to select the lncRNA related to autophagy.

Identification of the autophagy-associated lncRNA signature
Univariate Cox regression and the Kaplan-Meier approach were also utilized to screen for autophagyassociated lncRNAs with P <0.05. Then, the obtained lncRNAs were included in the least absolute shrinkage and selection operator (LASSO) regression. Subsequently, a multivariate Cox regression was further used to construct an autophagy-related lncRNA signature to predict prognosis in PCs [48]. Finally, a survival analysis was performed to estimate the prognostic value of the screened lncRNAs using the survival package in R 4.0.2.

Assessment of the autophagy-associated lncRNA signature
Based on the multivariate Cox regression's prognostic risk score, PCa patients were divided into high-risk and low-risk groups. The Kaplan-Meier survival curve and log-rank tests were applied to compare the differences in survival between high-risk and low-risk groups. The time-dependent receiver operating characteristic (ROC) curve and C-index calculation based on risk scores were conducted to assess the prognostic signature by the timeROC package and the survival package in R 4.0.2, respectively. Besides, the univariate and multivariate Cox regression analysis based on the prognostic signature and clinicopathological characteristics were further utilized to assess the independent prognostic factors for survival. Subsequently, the ROC curve was obtained by survivalROC package in R 4.0.2 to evaluate the relationship between clinicopathological features and the prognostic signature.

Gene set enrichment analysis
The TCGA-PRAD samples were divided into high-risk and low-risk groups based on the risk scores of the lncRNA signature. Gene set enrichment analysis (GSEA, http://www.broadinstitute.org/gsea/index.jsp) was conducted to identify the signaling pathways and biological processes [49]. The nominal (NOM) p-value <0.05 and false discovery rate (FDR) <25% were deemed noteworthy.

Statistical analysis
The data were processed and analyzed using the Perl data language (Version 5.30.2, http://www.perl.org) and R software (Version 4.0.2, https://www.r-project.org/). The survival curves were determined using the Kaplan-Meier method. The log-rank test was conducted to evaluate statistical differences between the high-risk and low-risk groups. A chi-squared test was used to investigate the relationship between risk scores and clinical features. The statistical tests were bilateral, and P <0.05 was considered statistically significant for all analyses.