Identification of a Novel Glycolysis-Related LncRNA Signature for Predicting Overall Survival in Patients With Bladder Cancer

Background Both lncRNAs and glycolysis are considered to be key influencing factors in the progression of bladder cancer (BCa). Studies have shown that glycolysis-related lncRNAs are an important factor affecting the overall survival and prognosis of patients with bladder cancer. In this study, a prognostic model of BCa patients was constructed based on glycolysis-related lncRNAs to provide a point of reference for clinical diagnosis and treatment decisions. Methods The transcriptome, clinical data, and glycolysis-related pathway gene sets of BCa patients were obtained from The Cancer Genome Atlas (TCGA) database and the Gene Set Enrichment Analysis (GSEA) official website. Next, differentially expressed glycolysis-related lncRNAs were screened out, glycolysis-related lncRNAs with prognostic significance were identified through LASSO regression analysis, and a risk scoring model was constructed through multivariate Cox regression analysis. Then, based on the median of the risk scores, all BCa patients were divided into either a high-risk or low-risk group. Kaplan-Meier (KM) survival analysis and the receiver operating characteristic (ROC) curve were used to evaluate the predictive power of the model. A nomogram prognostic model was then constructed based on clinical indicators and risk scores. A calibration chart, clinical decision curve, and ROC curve analysis were used to evaluate the predictive performance of the model, and the risk score of the prognostic model was verified using the TCGA data set. Finally, Gene Set Enrichment Analysis (GSEA) was performed on glycolysis-related lncRNAs. Results A total of 59 differentially expressed glycolysis-related lncRNAs were obtained from 411 bladder tumor tissues and 19 pericarcinomatous tissues, and 9 of those glycolysis-related lncRNAs (AC099850.3, AL589843.1, MAFG-DT, AC011503.2, NR2F1-AS1, AC078778.1, ZNF667-AS1, MNX1-AS1, and AC105942.1) were found to have prognostic significance. A signature was then constructed for predicting survival in BCa based on those 9 glycolysis-related lncRNAs. ROC curve analysis and a nomogram verified the accuracy of the signature. Conclusion Through this study, a novel prognostic prediction model for BCa was established based on 9 glycolysis-related lncRNAs that could effectively distinguish high-risk and low-risk BCa patients, and also provide a new point of reference for clinicians to make individualized treatment and review plans for patients with different levels of risk.


Methods:
The transcriptome, clinical data, and glycolysis-related pathway gene sets of BCa patients were obtained from The Cancer Genome Atlas (TCGA) database and the Gene Set Enrichment Analysis (GSEA) official website. Next, differentially expressed glycolysis-related lncRNAs were screened out, glycolysis-related lncRNAs with prognostic significance were identified through LASSO regression analysis, and a risk scoring model was constructed through multivariate Cox regression analysis. Then, based on the median of the risk scores, all BCa patients were divided into either a highrisk or low-risk group. Kaplan-Meier (KM) survival analysis and the receiver operating characteristic (ROC) curve were used to evaluate the predictive power of the model. A nomogram prognostic model was then constructed based on clinical indicators and risk scores. A calibration chart, clinical decision curve, and ROC curve analysis were used to evaluate the predictive performance of the model, and the risk score of the prognostic model was verified using the TCGA data set. Finally, Gene Set Enrichment Analysis (GSEA) was performed on glycolysis-related lncRNAs.

INTRODUCTION
Bladder cancer (BCa) is one of the most prevalent malignancies of the urinary system, according to statistics from 185 countries worldwide, in 2020, it ranked 9th in incidence and 13th in mortality (Sung et al., 2021). BCa is a heterogeneous disease and can be divided into two clinical types depending on whether it infiltrates the muscle layer: 25% of cases are muscle invasive BCa (MIBC), and 75% are non-muscle invasive BCa (NMIBC) (Cumberbatch et al., 2018). NMIBC has a low mortality rate, while 70% of these tumors will recur and 15% will progress in stage and grade. In contrast, the long-term survival rate of MIBC is unsatisfactory because it is prone to develop distant metastasis (Knowles and Hurst, 2015). Recent studies suggested that the poor prognosis of MIBC may also be related to the different molecular mechanisms of its occurrence and progression; some scholars even believed that MIBC and NMIBC were not the same kind of disease (Gui et al., 2011). With the increasing indepth study of molecular biological mechanisms and the rapid development of gene detection technology, molecular typing and prognostic evaluation of BCa by genetic testing is becoming a new diagnostic and therapeutic target.
Glycolysis is a metabolic pathway that provides energy to the body. Changes in the glycolytic state may signify the presence of cancer. One study by Otto Warburg found that the glycolysis of tumor cells produced more lactic acid and less ATP, which resulted in tissue hypoxia, even in an environment with a high oxygen content. This phenomenon is known as the "Warburg effect" (Warburg, 1956). Enhanced or abnormally activated glycolysis promoted tumor progression, which has been proven in breast cancer, colorectal cancer, liver cancer, etc. (Chen et al., 2019a;Weng et al., 2020;Zhang et al., 2020). In BCa, an increase of the Warburg effect has been proven to contribute to the aggressiveness of BCa tumor cells and accelerated the proliferation rate of BCa cells (Ritterson Lew et al., 2015). Therefore, changes in glycolytic status may be an emerging marker of malignancy and a potential prognostic target for patients with BCa.
LncRNAs are a class of nucleotides with a transcriptional length of more than 200 nucleotides. Although lncRNAs do not possess protein-coding ability, they can regulate gene expression through multiple layers and pathways, thus affecting the occurrence and progression of tumors, and can be used as markers for tumor diagnosis and prognosis (Evans et al., 2016).
There are more than 50,000 lncRNAs in the human genome, and they can play different roles in different tumors, acting as tumorsuppressor or tumor-promoter genes (Derrien et al., 2012). To effectively predict the prognosis of patients with BCa, the sensitivity and specificity of a single lncRNA was not sufficient; it was necessary to combine the high or low expression of multiple lncRNAs to have stronger diagnostic value.
Hyperglycolysis is a common phenomenon in the rapid proliferation of tumor cells and is also a marker of tumorigenesis. Recently, increasing amounts of data have shown that lncRNAs were a crucial regulator of glycolysis, and some studies have attempted to prevent tumor progression by inhibiting the activity of key enzymes in the tumor glycolysis pathway (Yang et al., 2014;Chen et al., 2019a;Liao et al., 2019;Wang et al., 2019a). The functional study of glycolysis-related lncRNAs has found that different lncRNAs could promote glycolysis in various ways, which in turn promoted tumorigenesis, progression, and metastasis. For example, lncRNA-SNHG7 regulated by c-MYC could increase glycolysis by inhibiting the expression of miR-34a-5p, thus promoting the progression of breast cancer , while the lncRNA LINRIS stabilized IGF2BP2 and promoted aerobic glycolysis in colorectal cancer (Wang et al., 2019b).
The functional study of glycolysis-related lncRNAs in bladder cancer has found that lncRNA-SLC16A1-AS1 could be used as a target of E2F1 and a co-activator to induce metabolic reprogramming in the progression of bladder cancer (Logotheti et al., 2020), and lncRNA UCA1 up-regulated hexokinase 2 in bladder tumor cells to promote glycolysis through the mTOR-STAT3/microRNA143 pathway (Li et al., 2014). Prognostic models constructed by differentially expressed glycolysis-related lncRNAs have been proven to have good predictive performance in gastrointestinal tumors, breast cancer, gliomas, and other tumors. In bladder cancer, there is currently no known glycolysisrelated lncRNA prognostic model.
This study aimed to clarify the relationship between glycolysisrelated lncRNAs and the prognosis of BCa. First, a total of 9 glycolysis-related lncRNAs were identified to have prognostic significance. Then, based on these 9 glycolysis-related lncRNAs, a novel prognostic prediction model for BCa was constructed and verified, which proved that the model possessed good predictive power. Finally, gene function enrichment analysis was performed on the differentially expressed glycolysis-related lncRNAs identified through this study to explore their potential functions in the progression of BCa.

Sample Sources and Processing
Transcriptome sequencing data and clinical data of BCa patients were downloaded from The Cancer Genome Atlas (TCGA) database 1 , and gene sets of glycolysis related functions and glycolysis related pathways were obtained from the Gene Set Enrichment Analysis (GSEA) official website 2 . The "edgeR" package of R software was used to normalize the entire data set, set | log2FC | > 1 and FDR < 0.05 as the threshold to construct the volcano map, and obtain the tumor tissue and pericarcinomatous differentially expressed glycolysis-related genes and lncRNA. The correlation between the lncRNAs and glycolysis-related genes was calculated using Pearson correlation. If the correlation coefficient of lncRNA | R2 | > 0.3 and P < 0.05, it was considered to be related to glycolysis and used for further analysis.

Construction of the Risk Score Model
Univariate Cox regression and Least Absolute Shrinkage and Selection Operator (LASSO) regression analysis were used to screen out 9 lncRNAs with prognostic significance. Genes with P < 0.05 were considered to be independent glycolysisrelated lncRNAs with prognostic significance. Multivariate Cox regression was used to calculate their respective coefficients (βI). Then, the formula of risk scoring model composed of βI and lncRNA expression level (EXPI) was obtained as follows: Risk score = 9 i=1 (β i * Exp i ). Then the risk score formula was used to calculate the risk score for each BCa patient. Subsequently, the median of the risk score was set as the cutoff value, and all eligible patients with bladder cancer were divided into highrisk and low-risk groups. Next, Kaplan-Meier survival analysis was performed to compare the differences in Overall Survival (OS) between the two groups of BCa patients. Univariate Cox regression and multivariate Cox regression analysis were used to evaluate whether risk score, age, gender, and TNM were independent prognostic factors for BCa. Clinical survival analysis was used to evaluate the predictive power of the risk scoring model for different clinical subgroups. ROC curve analysis was used to evaluate the efficiency of the prognostic model. These results were also verified by the validation data.

Establishment and Evaluation of the Prognostic Model
Then, a decision curve analysis (DCA) was performed to verify the predictive effect of the prognostic model.

Gene Set Enrichment Analysis
Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis was performed on the differentially expressed genes of BCa patients between the high-risk and low-risk groups. By analyzing the tumor signaling pathways and tumor progression processes involved in the differentially expressed genes, their possible roles and mechanisms in the occurrence and progression of BCa was further explored.

Differentially Expressed Glycolysis-Related Genes and LncRNAs in Bladder Cancer
The flow chart in Figure 1 shows the general step-by-step process of this study. First, the gene expression data, lncRNA sequencing data and clinical data of the corresponding BCa patients (411 tumor tissues, 19 pericarcinomatous tissues) were downloaded from the TCGA database, and glycolysis-related gene sets (GLYCOLYTIC_PRCOESS, HALLMARK_GLYCOLYSIS, REACTOME_GLYCOLYSIS, GLYCOLYSIS_GLUCONEOGENE SIS, and BIOCARTA_GLYCOLYSIS_PATHWAY) were downloaded from the GSEA database. The "edgeR" package of R software was used to normalize the entire data set, Set | log2FC | > 1 and false discovery rate (FDR) < 0.05 as the threshold to obtain glycolysis-related genes and lncRNAs that were differentially expressed between tumors and pericarcinomatous tissues (Figure 2).
Next, Pearson correlation analysis was used to analyze the correlation between the differentially expressed glycolysis-related genes and lncRNAs, and to set correlation coefficients | R2 | > 0.3 and P < 0.05 for lncRNAs that were considered to be related to glycolysis. A total of 59 glycolysis-related lncRNAs were identified as differentially expressed in BCa tumors and pericarcinomatous tissues (Supplementary Table 1). Subsequently, univariate Cox regression analysis was performed, through which 10 glycolysisrelated lncRNAs related to the prognosis of BCa were identified ( Figure 3A). Lastly, 9 key glycolysis-related lncRNAs associated with BCa prognosis were further screened by LASSO regression (Figures 3B,C).

Construction of a Risk Scoring Model for 9 Glycolysis-Related LncRNAs
The data of bladder cancer patients with a follow-up time of more than 30 days was selected to construct a risk scoring model, where their clinical parameters and pathological stages were showed in Supplementary Table 2. Patients were divided into training set and validation set at a ratio of 2:1. The risk score formula (Risk score = 9 i=1 (β i * Exp i )) was used to calculate the risk scores of all BCa patients. The median of the risk score was set as the cutoff value, and all eligible patients with BCa were divided into either high-risk or low-risk groups. Results of Kaplan-Meier survival analysis showed that the Overall Survival (OS) of the high-risk group of BCa patients was significantly lower than that of the low-risk group, both in the training set and the validation set (P = 1.546e-10, Figures 4A,B). The glycolysisrelated lncRNAs AC011503.2, AC078778.1, and ZNF667-AS1 were expressed at a lesser extent in tumor tissues of patients with BCa and were thus considered to be protective factors for the prognosis of bladder cancer, while glycolysis-related lncRNA AC099850.3, AL589843.1, MAFG -DT, NR2F1-AS1, MNX1-AS1, and AC105942.1 were highly expressed in tumor tissues and were therefore considered to be risk factors for the prognosis of BCa. ROC curve analysis proved that the model has a good predictive performance (3-year AUC = 0.753, 5-year AUC = 0.797) ( Figure 4C). We evaluated the survival times of patients in the high-and low-risk groups and found that mortality rates for patients with high-risk scores were higher than those with low-risk scores (Figures 4E,G). Heatmap analysis was performed to reveal the expression profiles of the 9 glycolysisrelated lncRNAs in high-risk and low-risk groups (Figure 4I). These results were then validated using validation set data (Figures 4B,D,F,H,J).

Construction and Evaluation of Prognostic Model of 9 Glycolysis-Related LncRNAs
A nomogram was then constructed to obtain the 3-year and 5year OS of BCa patients ( Figure 5A). This nomogram included the results of multivariate Cox regression, risk scores, and clinical information such as gender, age, and TNM stage. The calibration chart (Figures 5D,F) generated by using the "rms" package of the R software helped verify the performance of the nomogram. Next, the ROC curve ( Figure 5B) was used to obtain the area under the curve (3-year AUC = 0.781, 5-year AUC = 0.821) to evaluate the accuracy of the nomogram. These results further proved that the prognostic signatures of the 9 glycolysis-related lncRNAs were independent prognostic factors for bladder cancer. Subsequently, the C index was calculated (training set 0.79, validation set 0.724), and Decision Curve Analysis (Figures 5H,I) was conducted to verify the predictive effect of the prognostic model. Finally, the validation set data was used to verify these results (Figures 5C,E,G). Subgroup analysis results (Figure 6) showed that the model had good predictive efficacy in all age groups, in the T1-T2 and T3-T4 stages, with or without lymph node metastasis, and with or without distant metastasis (P < 0.05) (Figures 6A-I). There was no significant difference observed in the female group (P = 0.107) (Figure 6J). Possibly because of the small number of cases in the group.

Gene Function and Pathway Enrichment Analysis of 9 Glycolysis Related LncRNAs Signatures
Kyoto Encyclopedia of Genes and Genomes enrichment analysis was performed to explore the possible roles and mechanisms of the differentially expressed glycolysis-related lncRNAs in the occurrence and progression of tumors in the high-risk group and the low-risk group of BCa patients, Analysis results (Figures 7A-L) showed that these glycolysis-related FIGURE 3 | Screening glycolysis-related lncRNAs with prognostic value. (A) Risk ratio Forest plot identified 10 glycolysis-related lncRNAs with prognostic significance (P < 0.05). (B,C) LASSO regression further screened out 9 glycolysis-related lncRNAs that were significantly related to prognosis. lncRNAs were involved in multiple signaling pathways, such as: B-cell receptor, T-cell receptor, Hedgehog, MAPK, WNT, Calcium, and Chemokine. They were also involved in oxidative phosphorylation, ECM receptor interaction, and the development of basal cell carcinoma. These findings may help researchers determine the direction of further in-depth research in order to continue studying the mechanism of how glycolysis-related lncRNAs affects the progression of BCa.

DISCUSSION
Accurately predicting the prognosis of BCa patients is of great concern for clinicians in order to develop personalized treatment and review plans. For patients with a high risk of recurrence, we must be more cautious in the treatment strategy of bladder preservation, and for low-risk patients, we need to preserve the bladder as much as possible for the patients to help improve their postoperative quality of life. To determine whether a BCa patient was of a high-risk type, current clinical practice was mainly based on TNM staging, postoperative pathological diagnosis, and molecular classification, but there was no quantitative risk score or predictive model. Prognostic models constructed by differentially expressed glycolysis-related lncRNAs have been proven to have good predictive performance in gastrointestinal tumors, breast cancer, gliomas, and other tumors. Ho et al. (2021) constructed a glycolysis-related lncRNA prognosis model that could identify a subgroup of patients with different prognosis and different degrees of immune infiltration. In bladder cancer, there is currently no known glycolysis-related lncRNA prognostic model. This prediction model was established based on the following theoretical basis and previous research results: First, an increasing amount of evidence shows that changes in glucose metabolism FIGURE 4 | Construction of risk score model. Kaplan-Meier survival analysis showed that the Overall Survival (OS) of patients in the high-risk training set (A) and validation set (B) was significantly lower than that in the low-risk group. The ROC curve of the training set (C) showed that the AUCs of the 3-year and 5-year OS were 0.753 and 0.797, respectively, while the validation set (D) were 0.713 and 0.705, respectively. Survival rate and survival status of bladder cancer patients in the training set (E) and validation set (F). The distribution of 9 glycolysis-related lncRNAs risk scores for each patient in the training set (G) and validation set (H). Heatmap of 9 glycolysis-related lncRNAs in the low-risk group and the high-risk group in the training set (I) and validation set (J). are a sign of tumorigenesis. Some research has pointed out that the aerobic glycolysis of bladder cancer cells is more active than that of normal urothelial cells, and the increase in aerobic glycolysis is beneficial to the proliferation, invasion, and metastasis of bladder tumor cells (Berrondo et al., 2016). Other studies have found that the activation of the Warburg effect was positively correlated with the degree of tumor malignancy (Agnihotri and Zadeh, 2016). Secondly, previous studies have confirmed that lncRNAs were also involved in glycometabolic reprogramming of tumor cells in BCa. Different lncRNAs play different roles in the glycolysis process of bladder tumor cells.
Some lncRNAs accelerate tumor progression by promoting glycolysis, and some of them inhibit glycolysis to suppress cancer. For instance, Li et al. (2014) found that the lncRNA UCA1 could promote glycolysis by upregulation of HK2 through the mTOR-STAT3/miRNA and miR-143 pathways. Hu et al. (2017) found that the lncRNA CASC 8 may act as a tumor suppressor by reducing glycolysis in bladder tumor cells.
In this study, 9 glycolysis-related lncRNAs were screened in relation to the prognosis of BCa patients (AC099850.3, AL589843.1, MAFG-DT, AC011503.2, NR2F1-AS1, AC078778.1, ZNF667-AS1, MNX1-AS1, and AC105942.1). Among them, FIGURE 6 | The results of subgroup analysis. Analysis showed that the OS of bladder cancer patients in the high-risk group was significantly lower than that in the low-risk group based on age and TNM stage (P < 0.05). (A-I) In male patients, the OS in the high-risk group was also significantly shorter than that in the low-risk group (P < 0.001), while no significant difference was shown in female patients (P = 0.107) (J), possibly due to the relatively small number of cases.
glycolysis-related lncRNA AC011503.2, AC078778.1, and ZNF667-AS1 were protective factors for the prognosis of BCa, and glycolysis-related lncRNAs AC099850.3, AL589843.1, MAFG-DT, NR2F1-AS1, MNX1-AS1, and AC105942.1 were risk factors for the prognosis of bladder cancer. MNX1-AS1 regulates the expression of RAB1A in bladder tumor cells by competitively binding with miR-218-5p, thereby promoting the proliferation, migration, invasion, and epithelial-mesenchymal transition of bladder tumor cells, which contributes to the growth and metastasis of bladder cancer cells (Wang et al., 2020). Multivariate Cox analysis was used to calculate the regression coefficient and construct the prognostic model. According to the risk score, the median was used to divide the bladder cancer patients into high-risk and low-risk groups. Analysis showed that the OS in the low-risk group was longer than that in the high-risk group. According to the results of multivariate Cox regression, a histogram was established, including age, sex, TNM stage, and risk score. The prediction performance of the model was verified by both a calibration diagram and decision curve analysis, and the results demonstrated that the model had good predictive performance.
In terms of GSEA functional enrichment analysis, it was found that these glycolysis related lncRNAs were significantly enriched for the regulation of multiple processes, such as the B cell and T cell receptor signaling pathway, Hedgehog signaling pathway, MAPK signaling pathway, TGF-β signaling pathway, chemokine signaling pathway, EMC receptor interaction, and oxidative phosphorylation processes. Li et al. (2016) found that GALNT1mediated glycosylation and activation of Sonic Hedgehog signaling were involved in tumor-initiation and self-renewal of bladder cancer stem cells (Li et al., 2016). In another study, Teng Hu et al. found that RSPO3 promoted the aggressiveness of bladder cancer via the Wnt/β-catenin and Hedgehog signaling pathways . Transcriptional regulatory factor FIGURE 7 | Enrichment analysis. (A-L) KEGG pathway analysis indicates that these glycolysis related lncRNAs were involved in the B cell and T cell receptor signaling pathway, Hedgehog signaling pathway, MAPK signaling pathway, TGF-β signaling pathway, chemokine signaling pathway, EMC receptor interaction, and oxidative phosphorylation processes.
Frontiers in Genetics | www.frontiersin.org 10 August 2021 | Volume 12 | Article 720421 YAP was overexpressed in bladder cancer tissue and promoted the spread of bladder cancer by affecting the MAPK pathway (Qiu et al., 2020). All of these prior results further confirmed that the glycolysis associated lncRNAs screened in the present study play an important role in the occurrence and progression of BCa. It should be pointed out that this study has some limitations: First, this study was a retrospective study of patient data from the TCGA database, which may lead to selection bias. Secondly, the patient information obtained for this study lacks longer followup time; therefore, the predictive model still needs to be further verified in large-scale prospective clinical trials.

CONCLUSION
This study established a novel prognostic model of BCa based on 9 glycolysis-related lncRNAs, which can effectively distinguish high-risk and low-risk bladder cancer patients from a new perspective. It can also provide a new point of reference of ideas for clinicians to formulate individualized treatment and review plans for patients with different levels of risk. This study also found that some new glycolysis-related lncRNAs were related to the prognosis of bladder cancer, which points out the direction for the next step of verification experiments and exploration of mechanisms.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author.

AUTHOR CONTRIBUTIONS
YY designed the study and took overall control of the manuscript. ZZ analyzed and interpreted the data and drafted the manuscript. CL performed the data analysis and helped write the manuscript. WL collected relevant references, sorted out the data, and helped draft the manuscript. KM and CZ was involved in collecting the data and assisting with manuscript revision. All authors agreed to be accountable for the content of the work.

FUNDING
This work was supported by grants from Guangdong Provincial Clinical Research Center for Urological Diseases (2020B1111170006), from Guangdong Science and Technology Department (2020B1212060018).