Integrated Analysis of a Competing Endogenous RNA Network Reveals a Prognostic lncRNA Signature in Bladder Cancer

Long non-coding RNAs (lncRNAs) act as competing endogenous RNAs (ceRNAs) to regulate mRNA expression through sponging microRNA in tumorigenesis and progression. However, following the discovery of new RNA interaction, the differentially expressed RNAs and ceRNA regulatory network are required to update. Our study comprehensively analyzed the differentially expressed RNA and corresponding ceRNA network and thus constructed a potentially predictive tool for prognosis. “DESeq2” was used to perform differential expression analysis. Two hundred and six differentially expressed (DE) lncRNAs, 222 DE miRNAs, and 2,463 DE mRNAs were found in this study. The lncRNA-mRNA interactions in the miRcode database and the miRNA-mRNA interactions in the starBase, miRcode, and mirTarBase databases were searched, and a competing endogenous RNA (ceRNA) network with 186 nodes and 836 interactions was subsequently constructed. Aberrant expression patterns of lncRNA NR2F1-AS1 and lncRNA AC010168.2 were evaluated in two datasets (GSE89006, GSE31684), and real-time polymerase chain reaction was also performed to validate the expression pattern. Furthermore, we found that these two lncRNAs were independent prognostic biomarkers to generate a prognostic lncRNA signature by univariate and multivariate Cox analyses. According to the lncRNA signature, patients in the high-risk group were associated with a poor prognosis and validated by an external dataset. A novel genomic-clinicopathologic nomogram to improve prognosis prediction of bladder cancer was further plotted and calibrated. Our study deepens the understanding of the regulatory ceRNA network and provides an easy-to-do genomic-clinicopathological nomogram to predict the prognosis in patients with bladder cancer.


INTRODUCTION
Bladder cancer (BC) is the most common urological cancer and the major cause of cancer-related death globally, contributing to 165,100 deaths per year and resulting in significant morbidity and mortality (1). The occurrence and progression of bladder cancer are affected by a variety of intrinsic and extrinsic factors. The most common risk factors include smoking, occupational exposure, and the effects of environmental carcinogens (2). Although a large proportion of bladder cancer is non-muscle invasive bladder cancer (NMIBC) presenting at the first diagnosis, the risk of recurrence is nearly 30%, thus requiring a second transurethral resection of bladder tumor (3). Around 10% to 30% of patients with NMIBC will eventually progress into invasive disease (4). Radical cystectomy after neoadjuvant chemotherapy is a standard treatment for muscle invasive bladder cancer. The 5year overall survival rate is low (38.6%) among cT3-4aN0 patients after neoadjuvant chemotherapy and radical cystectomy (5). Recently, great progress has been made in the field of immunotherapy by PD-1/PD-L1 inhibitors, which are approved for treating patients with metastatic urothelial carcinoma. However, only one third of the late-stage patients have an initial response to immunotherapy, which may lead to remission (6). The effective options for treating patients with metastatic bladder cancer are extremely limited and required. Therefore, it is of great importance to investigate the molecular mechanisms underlying bladder cancer progression, and thus discover new molecular biomarkers of diagnosis and prognosis for patients with bladder cancer (7).
Long non-coding RNAs (lncRNAs) are a class of RNA defined as transcripts > 200 nt in length, which are transcribed from different regions of the genome (8). Increasing evidence showed that the expression of lncRNAs can not only serve as biomarkers of gene regulation in space, time, developmental stage but also regulate gene expression at the transcriptional, posttranscriptional, and chromosomal levels. lncRNAs are associated with a wide range of biological functions, such as transcriptional regulation, cell growth, and tumorigenesis (9,10). lncRNA, functioning as competitive endogenous RNA (ceRNA), competes for shared miRNA by cross-regulation among mRNA, circRNA, and pseudogenes (11,12). This miRNA-regulated lncRNA and mRNA network is a part of the "ceRNA hypothesis" (13). A large number of studies have shown that aberrant expression of lncRNA affects mRNA expression through sponge of miRNA, which is involved in tumorigenesis and cancer progression (8,14,15). The development of bladder cancer is accompanied by alterations in the expression patterns of many non-coding RNAs (ncRNAs) (16,17). Therefore, a better understanding of the regulatory lncRNA network is urgently required for a better functional investigation of the lncRNA.
In this study, we performed an integrated analysis of the transcriptome expression profiles and miRNA-seq data from multiple datasets and identified differentially expressed lncRNA, miRNA, and mRNA to generate a regulatory ceRNA network. We found that lncRNA NR2F1-AS1 and lncRNA AC010168.2 were independent prognostic biomarkers and used for constructing a lncRNA signature to predict the prognosis of bladder cancer. Moreover, a novel genomic-clinicopathological nomogram for predicting survival probability was generated. Our study would help to better understand the regulatory ceRNA network and provide an easy-to-do genomic-clinicopathological nomogram to predict the prognosis in patients with bladder cancer.

Data Collection and Preprocessing
RNA-seq transcriptomic data and miRNA-seq files were obtained from the TCGA database by using the "GDCRNATools" R package (18). Four hundred thirty-three RNA-seq files and 434 miRNA-seq files were downloaded on June 28, 2019. At the same time, clinical information and survival data (last follow-up time and all causes of death) of 412 patients were downloaded. A total of 403 patients were included in the study by eliminating cases with missing survival information and non-primary bladder cancer. Expression values were normalized by TMM method implemented in edgeR and further transformed by the voom method provided in limma. Low abundance genes (logcpm < 1 in more than half of the samples) were filtered out.
Datasets (GSE89006 and GSE55433) were downloaded from the Gene Expression Omnibus database (https://www.ncbi.nlm. nih.gov/gds/). GSE89006 contained four pairs of adjacent matched bladder tissue and bladder cancer and GSE55433 consisted of 26 bladder tissue samples from healthy controls and 57 bladder cancer samples. GSE89006 and GSE55433 were used as validation cohorts to exhibit the differential expression of lncRNA between bladder cancer and normal bladder tissue. GSE31684, which consisted of 93 bladder cancer patients with survival information, was used as a validation cohort for the prognostic value of lncRNA signature.

Preparation of Samples and Real-Time PCR
Three adjacent bladder tissues and six bladder cancer samples were collected from the Department of Urology, Second Xiangya Hospital of Central South University. Written informed consent was obtained from all patients, and this study was approved by the Clinical Research Ethics Committee of the Second Xiangya Hospital of Central South University. All methods were carried out in accordance with relevant guidelines and regulations (Declaration of Helsinki). Patient samples harvested by the surgical operations were frozen in liquid nitrogen. Real-time PCR was performed as previously described (19). Reverse transcription reaction was performed using Servicebio RT First Strand cDNA Synthesis Kit (G3330, Servicebio, China). Real-time PCR was cycled using 2× SYBR Green qPCR Master Mix (High ROX) (G3322, Servicebio, China). Primers are followed from 5′ to 3′: AC010168.2: CAGAAATGGAGGTGATGTGGC (F), CTCTCAGAG TTCTCAAGGCGTG (R); NR2F1-AS1: GTATTGACAGAGCAG GTAGATGAAAC (F), TTCTATTGCCAAAGCTCCCC(R); GAPDH: GGAAGCTTGTCATCAATGGAAATC (F), TGATGACCCTTTTGGCTCCC (R). This experiment was repeated three times.

Identification of Differentially Expressed mRNAs, miRNAs, and lncRNAs
Based on CRCh38.90 mapping from ENSEMBLE, mRNA and lncRNA were classified into two different files. Differentially expressed mRNA and lncRNA were analyzed to compare bladder cancer and adjacent bladder tissue by using the R package "DESeq2" with |logFC| ≥2 and FDR < 0.01 as the cutoff value. A volcano plot was used to visualize differentially expressed RNA analysis outcomes. The heatmap was generated by the "ggplot2" package. Regarding those patients with a different risk score, the differentially expressed genes were identified between the low-risk and high-risk groups using R package "limma" with "|logFC| ≥1.5 and FDR < 0.05" as the cutoff value.
Construction of the lncRNA-miRNA-mRNA ceRNA Network Differentially expressed lncRNA (DElncRNA), miRNA (DEmiRNA), and mRNA (DEmRNA) were used for further analysis, and lncRNA-miRNA interactions and mRNA-miRNA interactions were analyzed to search for shared miRNAs. As the official description of "GDCRNATools" R package, there are three criteria for determining the ceRNA between lncRNA and mRNA. First, lncRNA and mRNA share a certain number of miRNAs. Second, the correlation between lncRNA and mRNA expression was positive. Third, those shared miRNAs should have similar effects to regulate the expression levels of lncRNA and mRNA (http://bioconductor.org/packages/devel/bioc/vignettes/ GDCRNATools/inst/doc/GDCRNATools.html#competingendogenous-rnas-network-analysis). The lncRNA-miRNA pairs were investigated by miRcode (20) database, and the miRNA-mRNA pairs were used to identify three databases: starBase v2.0 (21), miRcode (20), and miRTarBase v7.0 (22). The Pearson correlation coefficient was evaluated to assess the linear association between lncRNA and mRNA. The included criteria were as follows: p-value of hypergeometric test and correlation was less than 0.01, and the absolute value of correlation was more than 0.4. Modulation of similarity and sensitivity correlations was used to calculate the possible regulatory effects of consensus miRNAs on lncRNA and mRNA. The function of "gdcExportNetwork" in this package was used to analyze the lncRNA-miRNA-mRNA interactions and the "edges" and "nodes" files were sequentially generated. Only the overlapped miRNA-mRNA interactions in these three databases were retained and the "edges" and "nodes" files were updated. Subsequently, all files were imported into Cytoscape (version 3.6.1) for the visualization of ceRNA network.

Functional Enrichment Analysis
Gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways were analyzed to investigate the function of differentially expressed mRNAs in the ceRNA network based on the Database for Annotation, Visualization, and Integrated Discovery (DAVID) (version 6.8, https://david. ncifcrf.gov/) (23). Differentially expressed genes between the lowrisk and the high-risk groups were analyzed using R packages "clusterProfiler" (24) and "DOSE" (25). These mRNAs were enriched in molecular function, cellular component, biological process, and signaling pathways, which was considered significant when p<0.05.

PPI Network Analysis for Protein-Protein Interactions
Differentially expressed mRNAs involved in the ceRNA network were subjected to the searching tool for the retrieval of interacting genes (STRING) for protein-protein interaction network analysis (version 11.0, https://string-db.org/) (26).
The database scored to predict the interactions of proteins. All settings were the default values. TSV file was downloaded, and then PPI network was generated by Cytoscape according to the PPI correlations. This PPI network was further analyzed to search clusters by an MCODE plug-in in the Cytoscape. Advanced options were as follows: degree cutoff: 2, haircut, nodes score cutoff: 0.2, K-core: 2, Max. Depth: 100.
Construction of lncRNA Signature Based on the ceRNA Network Univariate Cox proportional hazards regression analysis was performed to screen differentially expressed lncRNAs, which were associated with overall survival. Differentially expressed lncRNAs were selected for further analysis when p-value was less than 0.05. Multivariate Cox proportional hazards regression analysis was subsequently performed, and lncRNAs were considered as the significant difference when p value was less than 0.1. Based on these independent prognostic lncRNAs, the risk score of each patient with bladder cancer was calculated. The calculation formula of risk score equals to coef (LNC1) *exp (LNC1) +coef (LNC2) *exp (LNC2) +…+coef (LNCn)*exp (LNCn), in which coef (LNC) is a coefficient obtained by calculation of multivariate Cox analysis and exp (LNC) is normalized expression value. The median risk score was set as a cutoff, and patients were divided into two groups: low-risk group and high-risk group. Kaplan-Meier curve was plotted to predict the overall survival between the low-risk and high-risk groups. ROC curve was analyzed to evaluate the predictive power of this lncRNA signature for overall survival.

Gene Set Enrichment Analysis
According to patients with a risk score, Gene Set Enrichment Analysis (GSEA) (version 4.0.3, https://www.gsea-msigdb.org/ gsea/index.jsp) was performed to determine DEGs that were enriched in gene lists extracted from MSigDB (27). KEGG pathways that were enriched were determined by gene sets from the curated C2,cp.kegg.v7.0.symbols.gmt collections. The number of permutations was 1,000.

Estimate
Based on the expression data of TCGA, the level of stromal cells was presented, and the infiltration levels of immune cells in bladder cancer were calculated and downloaded from Estimation of STromal and Immune cells in MAlignant Tumor tissues using Expression data (ESTIMATE) (28) (version 1.0.13, https:// bioinformatics.mdanderson.org/estimate/). The differences in the stromal score and the immune score were analyzed between the low-risk and high-risk groups.

Clinical Application of lncRNA Signature
According to clinical characteristics, stratified survival analysis was performed, and Kaplan-Meier curves of clinical characteristics were plotted. Then, univariate and multivariate Cox regression proportional hazards regression analyses were performed to find independent prognostic factors of bladder cancer. In addition, a novel nomogram to predict the possibility of 1,3,5-year survival was established after stepwise regression analysis according to clinical characteristics by the package of rms (29). The nomogram was composed of four components: points line, risk factors line, survival probability line, and total points line. The length of each risk factor line reflects the regression coefficient estimated by multiple logistic regression analyses. The calibration was performed to evaluate the predictive value of the nomogram. To our knowledge, a calibration plot is drawn to verify the degree of calibration and an ideal line in the calibration plot is a 45-degree angle. The proposed nomogram has an accurate predictive ability when the actual graph is overlapped with the ideal line. All analyses were performed using R studio software version 1.2.5019 and R software version 3.6.2.

RESULTS
Identification of Differentially Expressed lncRNA, miRNA, and mRNA Two hundred and six differentially expressed lncRNA, 222 differentially expressed miRNA, and 2,463 differentially expressed mRNA were obtained by "DESeq2" analysis with cutoff (|logFC| ≥2 and p < 0.01). The volcano plot of mRNA and lncRNA was shown in Figure 1A, and the volcano plot of miRNA was displayed in Figure 1B. Heatmaps of lncRNA, miRNA, and mRNA are shown in Figures 1C-E, respectively. The clinical characteristics of patients with bladder cancer, including age, gender, race, pathologic stage, pathologic T, pathologic N, pathologic M, and vital status, were listed in Table S1.
Construction of a lncRNA-miRNA-mRNA ceRNA Network Based on miRcode database for lncRNA-miRNA interactions and starBase, miRcode, and mirTarBase databases for mRNA-miRNA interactions, 60 DElncRNA (44 upregulated lncRNA and 16 downregulated lncRNA), 66 DEmiRNA (42 upregulated miRNAs and 24 downregulated miRNA), and 60 DEmRNA (12 upregulated mRNA and 48 downregulated mRNA) were obtained as nodes, and 836 interactions were displayed as edges. A lncRNA-miRNA-mRNA ceRNA network was subsequently constructed by "Cytoscape" (Figure 2A). In this ceRNA network, we found that TIMP2 was a potential target of hsa-miR-200b/c-3p, serving as a  natural inhibitor of the matrix metalloproteinases and a potential non-invasive biomarker in kidney cancer and lung cancer. NR2F1-AS1 was also a potential target of hsa-miR-200b/c-3p and may act as a ceRNA in regulating the expression level of TIMP2 through competitive binding to hsa-miR-200b/c-3p.

GO and KEGG Pathway Analysis Revealed
That mRNAs in the ceRNA Network Were Enriched in Cancer-Associated Signaling Pathways The results from DAVID online tool (23) revealed that mRNAs in the ceRNA network were enriched in biological process, including "negative regulation of transcription from RNA polymerase II promoter", "cartilage development", "positive regulation of phospholipase C activity", "transcription, DNAtemplated", "lung development", "extracellular matrix organization", "integrin-mediated signaling pathway", and "xenophagy". The enriched cellular components were involved in "nucleus". The enriched molecular functions were involved in "sequence-specific DNA binding", "DNA-binding transcription factor activity", and "transcription corepressor activity". The enriched KEGG signaling pathways were involved in "microRNAs in cancer", "regulation of actin cytoskeleton", "pathways in cancer", "transcriptional misregulation in cancer", and "melanoma" (Table S2).

PPI Network Indicated That mRNAs in the ceRNA Network Were Further Clustered in Cancer-Associated Signaling Pathways
The protein-protein interactions network (26) was reconstructed by Cytoscape, and protein-coding mRNAs in the ceRNA network showed strong interactions ( Figure 2B). To further investigate clusters (highly interconnected regions), plug-in of MCODE was used in Cytoscape, and eight proteins (IL6ST, ERBB3, MITF, ZEB1, CTGF, TIMP2, ITGA5, and PLAU) were obtained in one cluster ( Figure 2C). Four proteins (PLAU, ITGA5, ZEB1, and ERBB3) were involved in "microRNAs in cancer" signaling pathway. Three proteins (PLAU, ZEB1, and MITF) were involved in the "transcriptional misregulation in cancer" (Figure 2C).

Validation of Differentially Expressed lncRNAs of NR2F1-AS1 and AC010168.2 as Independent Prognostic Biomarkers in Bladder Cancer
Among these differentially expressed lncRNAs in this ceRNA network, 27 lncRNAs were associated with the prognosis of overall survival through univariate Cox regression (p<0.05) ( Table S3). Subsequently, multivariate Cox analysis was performed for these 27 lncRNAs. Finally, two lncRNAs (NR2F1-AS1 and AC010168.2) were obtained for the construction of a prognostic lncRNA signature ( Figure 3A) ( Table S4). High expression of NR2F1-AS1 was associated with a poorer prognosis (HR, 1.16; p = 1.017e-02), whereas high expression of AC010168.2 was associated with a better prognosis (HR, 0.75; p = 1.494e-04) ( Figure 3B). Training datasets (TCGA-BLCA) and validation datasets (GSE89006, GSE31684) indicated that AC010168.2 expression was elevated in the high-risk group, whereas NR2F1-AS1 expression was decreased in the high-risk group, this was further validated by real-time PCR (Figures 3C and S1). The clinical information of training and validation cohorts are shown in Table 1.

Identification of a 2-lncRNA Prognostic Signature Based on the ceRNA Network
Based on the coefficient of NR2F1-AS1 and AC010168.2 of multivariate Cox analysis, a prognostic lncRNA signature was constructed. The formula of this signature was as follows: risk score= −0.2863566*exp (AC010168.2) + 0.1498769*exp (NR2F1-AS1). The lncRNA signature was associated with poorer overall survival in the high-risk group (p = 1.394e-04) ( Figure 4A). To investigate the sensitivity and specificity of this prognostic risk model, the ROC curve was plotted, and the area under curve (AUC) was 0.664 ( Figure 4B). Based on the median risk score as cutoff, 403 bladder cancer patients were divided into two groups ( Figure 4C), and more deaths occurred in the high-risk group ( Figure 4D). The expression level of NR2F1-AS1 was gradually elevated in pace with an increased risk score, but the expression level of AC010168.2 was decreased ( Figure 4E).

The lncRNA Signature Was Possibly Related to ECM Receptor Interaction in the Tumor Microenvironment
According to the lncRNA signature, risk-associated differentially expressed genes were analyzed to investigate GO function and KEGG signaling pathways. Five hundred seventy-two differentially expressed genes were obtained and used for further analysis. The main aspect of GO enrichment was related to the extracellular matrix. KEGG pathways enrichment analysis showed that lncRNA signature was involved in PI3K-Akt signaling pathway, focal adhesion, and ECM-receptor interaction ( Figure 5A). Furthermore, GSEA indicated that ECM-receptor interaction, cell adhesion molecules cams, focal adhesion and cytokine-cytokine receptor interaction were enriched in the high-risk group ( Figure 5B). Based on the heatmap of characterized signatures, our risk model could separate luminal (low risk) and basal (high risk) subtypes. Most of the signatures of ECM and smooth muscle, EMT and Claudin markers were displayed in the high-risk group ( Figure S2).

The lncRNA Signature Was Related to Cancer-Associated Fibroblast in Bladder Cancer
Estimation of STromal and Immune cells in MAlignant Tumor tissues using Expression data (ESTIMATE) analysis revealed that stromal score and immune score were both higher in the high-risk group than in the low-risk group ( Figure 6A). The tumor microenvironment is a complex mixture of tumor cells in the extracellular matrix (ECM), including diffusible growth factors, cytokines, and various types of stromal cells. It is believed that cancer-related fibroblasts (CAFs) are actively involved in the metastasis and invasion of tumor cells and CAFs are closely related to the poor prognosis of a variety of cancers (30). All CAFs expressed certain mesenchymal markers, including a-SMA, vimentin, and FAP. We found that the expression of ACTA2 (gene coded a-SMA), FAP, and ASPN were upregulated in the high-risk group, as compared to the low-risk group ( Figure 6B).

Validation of Prognostic Value of lncRNA Signature and Correlation Between lncRNA Signature and Biomarkers of Cancer-Associated Fibroblasts
GSE31684 consists of 93 bladder cancer patients managed by radical cystectomy (31). The detailed clinical information of GSE31684 is shown in Table 1. The survival information of patients and expression values of NR2F1-AS1, AC010168.2, ACTA2, FAP, and ASPN were extracted. According to the correlation coefficient, the lncRNA signature in GSE31684 was recalculated, and Kaplan-Meier curve of overall survival was displayed, indicating that the high-risk subgroup was associated with a poor prognosis (p = 1.667e-2) ( Figure 7A). The AUC was 0.573 ( Figure 7B). Risk score of lncRNA signature was positively correlated with ACTA2 (R=0.26, p=0.012), FAP (R=0.43, p=2.2e-05) and ASPN (R=0.30, p=0.0033) ( Figure 7C).  Stage  pTa  0  26  5  pTis  0  2  pT1  2  9  10  pT2-T4  399  19  78  pT2  129  17  pT3  138  42  pT4 132 19 The Prognostic lncRNA Signature Is an Independent Prognostic Factor in Bladder Cancer Stratified analysis of overall survival indicated that the high-risk score was associated with a poor prognosis in the subgroups of gender (male), pathologic stage (III-IV), pathologic N (N1-N3), pathologic M (M0) according to lncRNA signature ( Figure 8). Subsequently, univariate Cox analysis was performed. Age, pathologic stage, pathologic T, pathologic M, pathologic N, and lncRNA signature were significantly different and associated with prognosis. Multivariate Cox analysis revealed that only the lncRNA signature was an independent prognostic factor in bladder cancer ( Figure 9A).

A Novel Nomogram Was Generated for Clinical Application to Predict Overall Survival of Patients With Bladder Cancer
Based on the clinical characteristics of bladder cancer, we found four clinical characteristics (age, pathologic T, pathologic N, and risk signature) after stepwise Cox analysis for the establishment  of a novel genomic-clinicopathological nomogram to predict a 1-, 3-, 5-year survivals ( Figure 9B). The mortality score for each patient was calculated based on the scoring system at the bottom of the nomogram. A calibration of 1-, 3-, 5-year survival curves based on self-service resampling verification was consistent with nomogram predictions and was observed survival possibility ( Figure 9C).

DISCUSSION
Accurate predictions of overall survival could be extremely beneficial for patients with bladder cancer. Recently, many studies found varieties of diagnostic and prognostic biomarkers for cancer (32)(33)(34), few of which were utilized in clinical administration. A good biomarker can help identify cancers in the early stage, risk stratification, prognosis improvement, outcome prediction, variant histology (35), and targeted therapy (36). In this context, it is important to identify biomarkers that enable accurate prediction of recurrence, progression, and clinical outcome in bladder cancer. Thus, it is of great interest that we found prognostic lncRNA biomarkers and generated a lncRNA signature to evaluate 1-, 3-, 5year overall survivals in bladder cancer. A genomic-clinicopathological nomogram was established, and this could help physicians to predict the overall survival of patients with bladder cancer. lncRNA is a type of novel biomarker for clinical procedures involving the diagnosis and prognosis of bladder cancer (37). In many studies, lncRNA is also considered a potential target for cancer treatment (38). Although it was found that some lncRNAs were associated with a good prognosis in bladder cancer, the role and clinical significance of lncRNAs in bladder cancer remained largely unknown. The hypothesis is that ceRNA provides a new approach to study lncRNAs. ceRNA was studied in many biological processes. It has been shown that dysregulated ceRNA networks are closely related to tumorigenesis (39). lncRNAs can sponge with miRNAs to reduce the interaction of miRNA with mRNA. The ceRNA hypothesis provides a better understanding of the relationship among lncRNA, miRNA, and mRNA. In our study, we obtained 60 DElncRNA, 66 DEmiRNA, and 60 DEmRNA for generating a regulatory network of ceRNA. This ceRNA network could, thus, provide a new perspective to investigate the function of lncRNA. Based on the ceRNA network, two lncRNAs (NR2F1-AS1 and AC010168.2) were found and served as prognostic biomarkers. Previous studies have been shown that the expression of AC010168.2, also called RP11-174G6.5, may differ between the normal bladder and para-carcinoma tissue (40). More studies are focusing on NR2F1-AS1, which promotes the proliferation and migration in thyroid cancer by regulating the miRNA-338/CCND1 axis (41). NR2F1-AS1 targets SOX1 via sponging miR-363 to regulate the progression of endometrial cancer (42) and induces oxaliplatin resistance through miR-363 to target ABCC1 in hepatocellular carcinoma (43). We found that NR2F1-AS1 was positively correlated with ZEB1, ITGA5 (integrin subunit gene), and TIMP2 (gene of natural inhibitors of the matrix metalloproteinases), respectively. As the principal component of lncRNA signature, NR2F1-AS1 may regulate the integrin and degradation of extracellular matrix. It is suggested that lncRNA signature may play an important role in the tumor microenvironment through extracellular matrix, which is secreted by cancer-associated fibroblasts in bladder cancer. The results from KEGG pathways and GSEA also suggested a strong correlation between lncRNA signature and extracellular matrix. Furthermore, there is a cross-talk among TGF-b signaling pathways, integrins, and the extracellular matrix (44). TGF-b-associated extracellular matrix genes link cancer-associated fibroblasts to immune evasion and immunotherapy failure (45). TGF-b blockade  enhances specific responses to the immune-checkpoint blockade (46), which further provides an evidence that patients with this lncRNA signature could benefit from a potential treatment of TGFb blockade. Higher immune score and stromal score in the high-risk group implied therapeutic potential of combination treatments of immune checkpoint inhibitor and TGF-b blockade in metastatic bladder cancer. Cancer-specific survival was associated with a variety of clinical features, including age, higher tumor grade, higher pathological stage, lymph node metastasis, lymphovascular invasion, and soft tissue surgical margin (47). According to the lncRNA signature, a stratified analysis of overall survival indicated that a high-risk score was associated with a poor prognosis in the subgroups of gender (male), pathologic stage (III-IV), pathologic N (N1-N3), pathologic M (M0). The lncRNA signature is an independent prognostic factor in bladder cancer. Furthermore, different from other nomograms (48,49), we generated a novel genomicclinicopathological nomogram with four variables to predict mortality of patients with bladder cancer. The nomogram in this study may be a novel useful tool for clinical application.
There were some limitations in this study. First, further experiments were required for verifying the role of AC010168.2 and NR2F1-AS1 in affecting the prognosis of patients with bladder cancer. Second, following with more discovered interactions among lncRNA, miRNA, and mRNA, the lncRNA-miRNA-mRNA regulatory ceRNA network is dynamically updated.
In conclusion, we constructed a regulatory lncRNA-miRNA-mRNA ceRNA network and generated a novel genomicclinicopathological nomogram to predict survival rate. Our study may deepen a novel understanding of a regulatory ceRNA network and provide an easy-to-do genomic-clinicopathological nomogram to predict the prognosis of patients with bladder cancer.

DATA AVAILABILITY STATEMENT
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article/Supplementary Material.

ETHICS STATEMENT
The studies involving human participants were reviewed and approved by Clinical Research Ethics Committee of the Second Xiangya Hospital of Central South University. The patients/ participants provided their written informed consent to participate in this study.

AUTHOR CONTRIBUTIONS
MP and YW designed the project. MP and XC acquired data and wrote the manuscript. LY and YW performed surgical operations to obtain bladder cancer samples. WX performed real-time PCR. WX, LY, YW, and MP reviewed the manuscript. MP and YW supervised the study. All authors contributed to the article and approved the submitted version.