Autophagy-Related Long Non-coding RNA Is a Prognostic Indicator for Bladder Cancer

Bladder cancer (BC) is one of the most common malignant urinary system tumors, and its prognosis is poor. In recent years, autophagy has been closely linked to the development of BC. Therefore, we investigated the potential prognostic role of autophagy-related long non-coding RNA (lncRNA) in patients with BC. We obtained the lncRNA information and autophagy genes, respectively, from The Cancer Genome Atlas (TCGA) data set and the human autophagy database (HADb) and performed a co-expression analysis to identify autophagy gene-associated lncRNAs. Then, we divided the data into training group and testing group. In the training group, 15 autophagy-related lncRNAs were found to have a prognostic value (AC026369.3, USP30-as1, AC007991.2, AC104785.1, AC010503.4, AC037198.1, AC010331.1, AF131215.6, AC084357.2, THUMPD3-AS1, U62317.4, MAN1B1-DTt, AC024060.1, AL662844.4, and AC005229.4). The patients were divided into low-risk group and high-risk group based on the prognostic lncRNAs. The overall survival (OS) time for the high-risk group was shorter than that for the low-risk group [risk ratio (hazard ratio, HR) = 1.08, 95% CI: 1.06–1.10; p < 0.0001]. Using our model, the defined risk value can predict the prognosis of a patient. Next, the model was assessed in the TCGA testing group to further validate these results. A total of 203 patients with BC were recruited to verify the lncRNA characteristics. We divided these patients into high-risk group and low-risk group. The results of testing data set show that the survival time of high-risk patients is shorter than that of low-risk patients. In the training group, the area under the curve (AUC) was more than 0.7, indicating a high level of accuracy. The AUC for a risk model was greater than that for each clinical feature alone, indicating that the risk value of a model was the best indicator for predicting the prognosis. Further training data analysis showed that the gene set was significantly enriched in cancer-related pathways, including actin cytoskeleton regulation and gap junctions. In conclusion, our 15 autophagy-related lncRNAs have a prognostic potential for BC, and may play key roles in the biology of BC.


INTRODUCTION
Global Cancer Observatory (GCO) data (https://gco.iarc.fr/) (1) show that ∼550,000 people were diagnosed with bladder cancer (BC) in 2018, making BC the 10th common and 13th most fatal tumor worldwide. Of all patients with BC, ∼77% are male and BC is the sixth most common cancer in males worldwide. Furthermore, 53% of the patients with BC are more than 70 years old and the mortality is as high as 47%. The diagnosis of BC is mainly dependent on the findings of specimen biopsy, and the treatments for BC have generally included radical/partial cystectomy, neoadjuvant chemotherapy, adjuvant chemotherapy/radiation, and checkpoint inhibitorbased targeted therapies (2)(3)(4)(5)(6). BC is considered to be one of the tumors possessing the greatest economic burden on its treatment because of the difficulty of an accurate diagnosis and a precise treatment, inevitable frequent surveillance strategies, and the high recurrence risk (7,8). Therefore, timely diagnosis and clinical staging of BC for enabling reasonable and effective treatment measures become an urgent problem.
Autophagy is a process of elemental catabolism responsible for intracellular degradation. In the process of autophagy, excessive, FIGURE 1 | Prognostic network of long non-coding RNAs (lncRNAs) and co-expressed autophagy genes in bladder cancer (BC). In the center, the red node represents lncRNA, and the blue represents autophagy genes. The co-expression network was constructed by using CYTOSCAPE 3.7.2 software. aggregated, or damaged proteins and organelles are transferred to lysosomes for the degradation and reuse, thereby producing macromolecular molecules and energy (9,10). Ishaq and his colleagues (11) provided a good overview of the role of autophagy in tumors, describing that autophagy can exert anticancer and pro-cancer dual effects by affecting the web of inter-organelle membrane contact sites, cancer stem cells, cancer metabolism, and tumor immunity (11). Moreover, autophagy is also involved in gene repair and supports the synthesis of DNA in response to oncogenes (12,13). Few studies show that microRNA (miRNA) 516a exerts carcinogenic effects in BC by attenuating beclin-1 (BECN1)-dependent autophagy (14), the miR-139-5p/Bmi-1 axis is closely related to AMP activated protein kinase (AMPK) and mammalian target of rapamycin (mTOR) pathway-dependent autophagy in human BC cells (15), and autophagy-related gene 7 (ATG7) promotes the invasion of BC via an autophagy-mediated increase in the rho GDP dissociation inhibitor beta (ARHGDIB) messenger RNA (mRNA) stability (16). Although the role of autophagy and targeted therapy strategies based on this process have gradually become a research hotspot in BC in recent years, the involvement of specific mechanisms and checkpoints remains controversial. Since the development of high-throughput genomics, long non-coding RNA (lncRNA) has been validated as a non-protein coding RNA comprising different types of RNA polymerase II (Pol II)-transcribed molecules longer than 200 nucleotides (17,18). lncRNA participates in the occurrence and development of cancer by regulating biological processes including cell growth, cell cycle, cell metastasis, cell death, angiogenesis, metabolic disorders, drug resistance, immune escape, DNA damage, and cell stemness (19,20). Many reports have corroborated the association between lncRNA and autophagy in cancer. BLACAT1, MALAT1, XIST, SNHGs, HULC, CASC2, GAS5, H19, and HOTAIR all are lncRNAs that mediate chemoresistance in the tumor by modulating autophagy (19). The expression of lncRNA maternally expressed 3 (MEG3) leads to a reduced autophagy activation, which reduces the proliferation of BC cells (21). lncRNA urothelial cancer associated 1 (UCA1) acts as an oncogene and promotes malignant progression of BC through the UCA1-miR-582-5p-ATG7 autophagy signaling pathway (22). Therefore, autophagy-related lncRNAs may be used as prognostic factors in patients with BC and also as potential therapeutic targets. Here, we strived to establish a prognostic system for autophagy-related lncRNA in BC and to promote the targeted therapy of BC.

Screening for lncRNA Co-expressed With Autophagy-Related Genes
Autophagy-related genes were extracted from the human autophagy database (HADb, http://autophagy.lu/clustering/ index.html). The Limma R package was used to extract the autophagy-related gene expression information from the TCGA-BLCA RNA sequencing data. To identify the lncRNAs-targeting autophagy-related genes, we analyzed the Pearson correlation between the lncRNA and autophagy-related gene expression levels using the screening criteria: (| R2 | > 0.3 and p < 0.001). Finally, 1435 lncRNAs co-expressed with autophagy-related genes were identified in a TCGA-BLCA data set.

Establishment of a Risk Score Model and Building a Risk Curve
We use the survival R package to construct the lncRNA signature. The BLCA sample criteria for identifying and verifying lncRNA characteristics were: (1) complete lncRNA expression value; and (2) complete clinical characteristics (the survival time and survival status). A total of 406 BLCA samples were included for further construction of lncRNA characteristics. We randomly divided the 406 BLCA samples into the training group (n = 203) and testing group (n = 203). The training cohort was used to develop the lncRNA features. The experimental cohort was used to verify the lncRNA characteristics. Secondly, 406 lncRNAs were selected from a training cohort for univariate Cox regression analysis to identify prognosisrelated lncRNAs (p < 0.005). Univariate Cox regression analysis identified 45 prognostic lncRNAs. We used the Akaike information criterion (AIC) method to select the optimal model from 45 prognostic lncRNAs. Finally, we selected 15 lncRNAs with the lowest AIC values to build a prediction model.
Multivariate Cox regression analysis was performed on 15 prognostic lncRNAs to determine their prognostic characteristics and calculate the correlation coefficient (24,25). The risk score for the prognostic lncRNA characteristics of each patient was FIGURE 3 | Kaplan-Meier (KM) survival curve of the autophagy-related lncRNA BC risk score in TCGA training group. In the TCGA data, the 5-year survival rate of high-risk patients is lower than that of low-risk patients. calculated by using the following formula: Risk score = expression lncRNA1 × coefficient lncRNA1 The risk score of prognostic lncRNAs was calculated by using a linear combination of lncRNA expression levels weighted by a regression coefficient (β). β was calculated by a log conversion of the hazard ratio (HR) obtained from multivariate Cox regression analysis (26). Low-risk group and high-risk group were determined by using the median-risk score.

Risk Model Prediction Ability in Training Group and Testing Group
Survival and phopmap packages were used to evaluate the predictive value of the risk model in a training cohort (n = 203).  The training cohort was divided into low-risk group and highrisk group based on the median risk score. The Kaplan-Meier (KM) survival curve was drawn to compare the difference in the overall survival (OS) between the two groups and a risk curve was drawn. p < 0.05 was considered as statistically significant. Subsequently, the training group was used to sort the complete data (n = 172) for each clinical feature. These clinical features include age, gender, union for international cancer control (UICC) stage, histologic grade, and pathological TNM (pTNM) stage. Univariate and multivariate Cox regression analyses were  used to evaluate the predictive power of the risk model. Then, survivroc was used to draw the time-dependent receiver operating characteristic (ROC) curve and area under the curve (AUC) values to evaluate the predictive value of the comparative risk model and the assessed clinical features. p < 0.05 was considered as statistically significant. We then conducted a correlation analysis between the model risk score and the clinical characteristics of the training group. Each clinical feature was divided into two categories (e.g., dividing age into > 65 and ≤ 65, and dividing the grade into high and low). SPSS 25 software was used for statistical analyses.
The same method was used in the testing group to verify the stability and reliability of the risk model. The testing group (n = 203) was divided into low-risk group and high-risk group based on the median risk. The KM survival and risk curves were drawn. The testing group data with complete clinical data was selected (n FIGURE 7 | KM survival curve of the autophagy-related lncRNA risk score for BC in TCGA testing group. TCGA data shows that the 5-year survival rate of high-risk patients is lower than that of low-risk patients (log rank p < 0.01).
FIGURE 8 | Multi-factor receiver operating curve details. The AUC value of the testing group is significantly more than 0.7, which is greater than the predicted value of clinical data.
= 172), and the time-dependent ROC curve and AUC value were determined to evaluate the predictive value of the comparative risk model and these clinical features.

Gene Set Enrichment Analysis
Gene set enrichment analysis (GSEA) was used to determine the gene set functions. The aim was to use a predefined set of genes, sort the genes according to the differential expression levels in the two sample groups, and check whether the gene set was enriched at the top or bottom of this sorting table.

Statistical Analysis
Cytoscape software (version 3.7.2) was used to construct an autophagy-related lncRNA co-expression network. R language (version x64 3.6.3, survival library) was used to estimate the survival rate, generate a survival curve, and perform the Cox regression analysis. We used the R language (pheatmap library) to generate the risk value, survival time, and risk heat maps. We used the R language (survival ROC library) to draw a multiindex ROC curve. GSEA (version 4.0.3) was used to distinguish between the two sets of functional annotations. SPSS 25 was used to analyze the correlation between the risk model and clinical features. Statistical significance was set at the threshold of p < 0.05 for both tails.

Autophagy-Related Gene lncRNAs Co-expression Network Construction
We extracted the lncRNA data from the TCGA database transcriptome data set. A total of 232 autophagy-related genes were extracted from the HADb (http://autophagy.lu/clustering/ index.html). We constructed an autophagy-related lncRNA coexpression network and performed a co-expression analysis to identify autophagy-related lncRNAs. The screening criteria were (| R2 | > 0.3 and p < 0.001). A co-expression network diagram with autophagy genes was drawn by using the 15 prognostic-related lncRNAs identified by multivariate Cox analysis (Figure 1).

Selection of 15 Autophagy-Related Prognostic lncRNAs for Patients With BC
First, we used univariate Cox regression analysis based on the identified autophagy-related lncRNAs to initially screen the prognostic genes. We used a value of p = 0.005 as a cutoff value, and the lncRNAs that met this criterion were used to judge prognostic indicators. Our training set included 203 patients with BC from the TCGA data set. A total of 45 lncRNAs had a prognostic value for patients with BC (p < 0.01) ( Table 1). Then, we used the AIC method to select the best model from the 45 prognostic lncRNAs. About 15 autophagy-related lncRNAs were found to be suitable for model building. Of these lncRNAs, 5 were unfavorable prognostic factors (AC037198.  Table 2). These 15 autophagy-related lncRNAs were used as a prognostic model.

The Effect of Autophagy-Related lncRNA on BC Prognosis
The risk scoring method was used to assess the prognosis of BC in the training group (n = 203). We divided patients with BC into low-risk group and high-risk group based on the median risk score (Figure 2). The risk score can significantly predict the OS in patients with BC, with OS being longer in the low-risk group than in the high-risk group. Additionally, a KM survival curve for BC autophagy-related lncRNA risk scores was plotted (Figure 3). The KM survival curve showed that the 5-year survival rate was greater in the low-risk group than that in the high-risk group (log rank p < 0.01). We screened the training group data with complete clinical characteristics for the next analysis (n = 172) and used multivariate Cox regression analysis to explore if the risk score is an independent predictor of prognosis in patients with BC. The HR was 1.080, which indicates that the risk score model can significantly help predict the survival in patients with BC, and eliminates the effects of other factors (such as age, gender, and grade; Table 3). Next, the ROC curve of clinical data and risk score was drawn. The AUC of the ROC curve of the risk score was 0.865, which is > 0.7 and indicates that the model has good accuracy. The AUC of the model is greater than that of any other clinical characteristics, which indicate that the model is superior to other clinical traits in prognostic prediction (Figure 4).

Clinical Correlation Between Autophagy-Related lncRNAs and Patients With BC
Subsequently, we determined the clinical value of the 15 autophagy-related lncRNAs in terms of age, gender, grade, and pTNM stage ( Table 4). The risk score is closely related to grade, UICC stage, and pT stage (Figure 5), but is not strongly related to age, gender, pN stage, and pM stage. Additionally, the risk score tends to increase in higher grades, which indicates that the lncRNA signature may be related to BC progression.

Validation in the TCGA Data Set
Next, this model was validated by using the TCGA testing group. A total of 203 patients with BC were recruited to verify the lncRNA characteristics (Figure 6). We divided the patients with BC into low-risk group and high-risk group based on the median risk score. Consistent with the results obtained from the TCGA data set training group, the survival time was shorter in highrisk patients than in low-risk patients (Figure 7). Testing group data were then screened with complete clinical characteristics (n = 172), and the time-related ROC curve was drawn (Figure 8).
The AUC value was 0.833, which is >0.7, and indicates that the risk model has great predictive ability and that the predictive ability of the risk model is greater than that of other clinical features. Together, these results show that the established lncRNA signature provides reliable prognostic predictions for patients with BC.

Gene Set Enrichment Analysis
Further functional annotations were made through GSEA. The results showed that differentially expressed genes between the two groups were enriched in tumor-related pathways. A total of 11 gene sets were significantly enriched when the nominal value of p < 0.05. We selected seven pathways closely related to tumor metastasis and invasion, including actin cytoskeleton regulation, gap junctions, focal adhesion, and extracellular matrix (ECM) receptor interaction ( Table 5). These are important signaling pathways related to tumorigenesis and cancer progression, including in small cell lung cancer, which is highlighted by the relationship between the Wnt signaling pathway and glioma (Figure 9). In conclusion, the identified autophagy-related genes contribute to important cancer pathways. This may provide a new direction for the treatment of BC.

DISCUSSION
Bladder cancer is defined as malignant tumors that occur on the bladder mucosa. BC is the most common malignant tumor of the urinary system (27). There are currently some BC prognostic indicators, including tumor grade, UICC stage, and pTNM stage (28). However, additional prognostic factors are required to more accurately predict and analyze the survival rate in patients with BC. At present, a large number of lncRNAs have been analyzed in the diagnosis and prognosis of various cancers. Long non-coding RNA Erbb4-IR is downregulated in prostate cancer and can predict the prognosis (29). Overexpression of lncRNA ROR1AS1 predicts a poor prognosis and promotes cervical cancer metastasis by activating the Wnt/β-catenin/EMT signaling cascade (30). lncRNA-D16366 is a potential biomarker for the diagnosis and prognosis of hepatocellular carcinoma (31). lncRNA PVT1 can predict a prognosis in patients with prostate cancer and regulate the tumor growth (32). lncRNA FENDRR and FOXF1 are also prognostic factors for the survival in patients with lung adenocarcinoma (33). In addition, many reports confirm a link between lncRNA and autophagy in cancer.
Reduced expression of lncRNA PRRT3-AS1 can inhibit the prostate cancer cell proliferation and promote apoptosis and autophagy (34). In acute myeloid leukemia, lncRNA UCA1 promotes autophagy by targeting miR-96-5p (35). lncRNA MALAT1 activates autophagy and promotes the cell proliferation by downregulating the miRNA-204 expression in gastric cancer (36). However, to-date, there has been no systematic and comprehensive analysis of autophagy-related lncRNA in BC. This study is the first to use bioinformatics to analyze the value of autophagy-related lncRNA in the diagnosis and prognosis of BC. These results are expected to provide new ideas and research directions for the diagnosis and treatment of BC.
In this study, we collected TCGA data sets and divided them into training group and testing group. In the training group, autophagy-related lncRNAs were used to investigate the prognosis of BC in patients. We identified autophagyrelated lncRNAs by constructing a lncRNA autophagy gene coexpression network. We identified a signature of 15 autophagyrelated lncRNAs and divided the patients with BC into low-risk group and a high-risk group based on the median risk score. The 5-year survival rate of patients in the high-risk group was shorter than that of patients in the low-risk group. Using univariate and multivariate Cox regression analyses, we showed that the risk score is an independent BC prognostic prediction factor. Subsequently, a multi-factor ROC curve was made, and its AUC was 0.865, which was greater than the AUC of each clinical feature. Together, these results show that the lncRNA signature can better predict a prognosis in patients with BC than can other clinical features. Then, the model was used to predict a prognosis in the testing group, and satisfactory results were obtained, which prove the effectiveness of the model.
Previous studies have shown that THUMPD3-AS1 is associated with non-small cell lung cancer and regulates self-renewal through miR-543 and ONECUT2 (37). The Wnt/β-catenin signaling pathway promotes the progression of liver cancer through the activation of SUMO1P3 lncRNA by targeting miR-320a (38). A large number of studies also indicate that lncRNA can directly or indirectly regulate the development of cancer through the regulation of miRNA. Therefore, lncRNA-mRNA co-expression analyses are essential to fully evaluate the function of lncRNAs in cancer. Additionally, the autophagy-related functions of lncRNAs occur through the regulation of mRNA. Therefore, it is important to analyze the specific function of lncRNAs by examining the co-expression of lncRNA-mRNA (39,40). Ultimately, our results indicate that the 15 autophagy-related lncRNAs identified in this study are potential BC therapeutic targets.
By constructing an autophagy lncRNA co-expression network, we identified a risk model of 15 autophagy-related lncRNAs. Using this model, we were able to separate the group as low risk and high risk based on the median-risk score. This risk score has a prognostic value for patients with BC. In the future, this lncRNA model can be used for risk scoring to predict the prognosis of a patient. However, our research contains several limitations. Bioinformatics methods and public databases were used for this study. Although the stability of the autophagy-related lncRNA signature was successfully validated, the exact molecular mechanisms of these lncRNAs have not been experimentally investigated and their prognostic values have not been experimentally proven. Therefore, further experimental studies using a larger sample size are required to verify these results.

DATA AVAILABILITY STATEMENT
The data that support the findings of this study are available from the corresponding author upon reasonable request.

AUTHOR CONTRIBUTIONS
JW and CG: conceptualization, methodology, validation, formal analysis, investigation, data curation, writing of the original draft, and reviewing and editing. HF, ZX, and YH: validation, formal analysis, investigation, data curation, and visualization. YL: conceptualization, methodology, investigation, data curation, reviewing and editing the manuscript, visualization, supervision, and project administration. All authors: agree to be accountable for the content of the work.