A signature of glucose metabolism and DNA methylation improves prognostic prediction for urinary bladder cancer

Background: Glucose metabolism and DNA methylation have been previously proved to play important roles in cancers. However, few studies have explored the potential connections between abnormal glucose metabolism and DNA methylation in tumor progression. In this study, we aimed to screen out glucose metabolism-related genes that were DNA methylation associated in order to establish a promising signature for the prognosis of bladder cancer (BLCA). Methods : With transcriptome data of BLCA samples from The Cancer Genome Atlas (TCGA) and methylation data downloaded from TCGA data base 450K microarray, glucose metabolism-related genes that were simultaneously methylation modulated and survival associated were obtained and a prognostic signature was established. Gene set enrichment analysis (GSEA) and weighted gene co-expression network analysis (WGCNA) were performed for identification of gene and pathway enrichment. Two most methylation-related risky genes were selected to conduct functional validations. Results : 18 target genes were identified and the samples were survival-significantly divided into two groups by using the signature. Via multivariate Cox analysis, we considered the signature an independent prognostic factor. Several biological processes and signaling pathways were enriched in the high-risk group according to GSEA results and three modules of genes were identified by WGCNA. We demonstrated that UCHL1 and PYCR1 promoted proliferation, migration and invasion ability of bladder cancer cells. Conclusions : The signature based on the 18 target genes could serve as an independent prognostic factor for urinary bladder cancer patients.

predictive and prognostic biomarkers, which has been a new challenging strategy for researchers and urologists.
Dramatic changes in energy metabolism and nutrient uptake are commonly observed in various cancer cells with enhanced cell proliferation and growth to support. As one of the most important nutrients for cancer cells uptake to maintain unrestricted proliferation, not only can glucose be utilized to synthesize ATP but also to provide cancer cells with carbon source [5][6]. One of the most well-known glucose metabolic reprograming forms is "the Warburg effect" [5,7], referring to the phenomenon that most cancer cells generate ATP by glycolysis even under aerobic conditions while it takes place only in nonaerobic conditions for normal cells. Altered expressions of many rate-limiting enzymes participate in activation of Warburg effect [5] and a large number of oncogenes and tumorsuppressors are involved in a variety forms of energy reprogramming processes [8][9][10][11][12].However, despite the progress previous studies have achieved, we still lack comprehensive analysis and summary of glucose metabolism-related genes especially those correlated to cancer prognosis, on the basis of which potential new diagnosis and treatment targets are likely to be established.
In previous studies, abnormal DNA methylation has been proved to be related to cancers and other diseases by regulating genomic transcription. Usually detected in highly and moderate repeated DNA sequences, hypomethylation can result in chromosome instability leading to activated oncogenes [13][14][15], while DNA hypermethylation in the promotor region can lead to silencing of some tumor suppressing genes [13,16]. For some genes, methylation has been found to be responsible for early carcinogenesis [17]. Thus, it has become a trend for DNA methylation to act as early diagnostic and prognostic biomarkers for tumors [18][19][20] because of its relative stability [21] and easy access to detection in blood. For instance, several studies have proved that hepatocellular carcinoma patients could be discriminated from healthy individuals with high sensitivities by detecting methylation of CpGs [22][23][24]. An FDA approved project aimed at early diagnosis of colorectal cancer by detecting promotor methylation status of SEPT9 in blood and the sensitivity and specificity were found to be 36.6 to 95.6% and 81.5 to 99.0% [25]. Considering the critical role glucose metabolism plays in cancer progressions, we intended to establish a new prognostic indicator for bladder cancer patients by exploring the associations between DNA methylations and glucose metabolism-related genes.
In our study, TCGA BLCA cohort was used as study object. We screened out 46 glucose metabolismrelated genes that were linked to overall survival, among which 18 genes' expressions were significantly correlated to their corresponding promotor methylation strengths. A signature based on the 18 genes was established to predict the BLCA prognosis. The enrichments of co-expressive genes, biological functions and signaling pathways were constructed and prove the accuracy of our bioinformatic analyses and to provide experimental bases for further research.

Materials And Methods
Glucose metabolism-related gene set From the Molecular Signatures Database (http://www.broad.mit.edu/gsea/msigdb/), we selected all the pathways related to glucose metabolism with relevant genes included.

Patient Samples
The RNA-seq transcriptome data of Bladder Urothelial Carcinoma (BLCA) cohort was downloaded from The Cancer Genome Atlas (TCGA) data base (https://portal.gdc.cancer.gov/) along with corresponding clinical information.

Consensus Clustering Analysis
Cancer samples from TCGA BLCA cohort included were clustered into groups by consensus expression of 30 glucose metabolism-related genes with the largest values of median absolute deviation (MAD) by using "ConsensusClusterPlus" in R program 3.3.2 (www.rproject.org). Chi-square test was used to compare the distribution of age, tumor grade, pathology stage, lymph node positivity, metastasis positivity and neoplasm grade between different clusters.
Identification Of Survival-related Differentially Expressed Genes (degs) 300 and 100 samples from TCGA BLCA cohort were randomly divided into training set and validation set, respectively. All the steps of gene screening were carried out in training set. Differentially expressed genes (DEGs) were screened out from all the glucose metabolism-related genes by comparing bladder cancer patients with overall survival time less than 1 year with those who lived longer than 3 years. P < 0.05 was the considered the criterion. With the use of a univariate Cox proportional hazard regression analysis, we analyzed all the glucose metabolism-related genes to obtain genes remarkably correlated with patients' survival time with p < 0.05 as the main cutoff point.
The overlapping genes gained from the two steps above were considered survival-related DEGs.

Analysis Of Dna Methylation Status
Methylation data was downloaded from TCGA database 450K microarray. The correlations between the expressions of survival-related DEGs and the strengths of promoter methylation probes of corresponding genes were analyzed to divide the probes into significant probes and non-significant probes after z-score normalization with P < 0.05 as criterion. Obtained genes were considered target genes for further steps of our study.

Analysis Of Prognostic Signature
The risk score of each patient in the training set and the validation set was calculated by using gene expression values and regression coefficients (β) which were generated from the univariate Cox proportional hazard regression analysis, according to the formula below: risk score = β(gene 1) × gene was performed in TCGA BLCA dataset. The parameter of permutations was at 1,000 and FWER adjusted P-value < 0.02 was considered the significance threshold. Hub gene identification and coexpression analysis were performed by co-expression network analysis WGCNA [26], which was visualized by Cytoscape [27].

Cell Lines And Cell Transfection
Bladder cancer cell lines, UM-UC-3 and T24, were purchased from the Chinese Academy of Sciences Cell Bank (China). At 37 °C in an atmosphere of 5% CO2, the two cell lines were cultured in RPMI 1640

Cell Proliferation Assay
Cell proliferation ability was measured by using a BeyoClick TM EdU-488 cell proliferation kit (Beyotime Biotechnology, China). Transfected cells were co-incubated with EDU working solution (1:1000) in 6well plates at 37 °C in an atmosphere of 5% CO 2 for 2 hours, after which the culture medium was removed and 4% paraformaldehyde was added for fixation at room temperature for 20 minutes. Then cells were incubated at room temperature for 15 minutes with 0.3% Triton X-100 in PBS, followed by incubation with click reaction solution according to the manufacture's protocol. Images were obtained by A fluorescence microscope (Olympus Corporation, Japan). Matrigel were used to measure cell migration. After incubated at 37 °C in an atmosphere of 5% CO 2 for 24 hours, cells remaining in the upper membranes were removed while cells that had migrated to the lower side were fixed by 4% paraformaldehyde and stained by 1.0% crystal violet. Pictures were obtained by using microscope and cell counting was conducted by ImageJ.

Immunoblotting Assay
Total protein was extracted from cells by using RIPA lysis buffer containing 1% PMSF. Bicinchoninic acid assay (Beyotime Institute of Biotechnology) was applied to determine protein concentrations.

Statistical analysis
All statistical analyses for bioinformatic analysis were performed by the R software (R version 3.3.2).
As for functional validations, GraphPad Prism version 7.0 (La Jolla, CA, USA) was used for statistical analyses. Student's t test was performed for comparisons. Each of the experiments were conducted for at least three times and each P value < 0.05 was considered to indicate significantly statistically.

BLCA data acquisition
The flowchart of the study procedure was shown in Fig. 1. Twelve pathways linked to glucose metabolism were chosen with a total of 606 glucose metabolism-related genes included for following study. The RNA-seq transcriptome data of 405 bladder cancer samples and 19 normal bladder samples from TCGA BLCA cohort was collected and corresponding DNA methylation data was obtained from TCGA 450 k microarray.
Consensus clustering of glucose metabolism-related genes identified two main clusters of BLCA samples With the increase of k value, the 405 bladder cancer samples were gradually and steadily divided into two main large clusters by 30 glucose metabolism-related genes with the largest MAD values, the light green cluster and the light blue cluster (Fig. 2a-j), which suggested that glucose metabolism could play a potential role as an index in classifying bladder cancer patients. Then we tried to categorize the 30 genes into two clusters according to their expression trends, namely glucose cluster 1 and glucose cluster 2. The result showed 405 samples were obviously divided into two clusters, cluster A and cluster B, which were characterized with significant differences on the expressions of cluster 1 genes and cluster 2 genes (Fig. 2k). By measuring the associations between the clustering and clinicopathological features, we detected significant difference between cluster 1 and cluster 2 for neoplasm grade (P = 2.999e-06), while no obvious differences were found for all of the other parameters (Table. 1), which indicated that compared to current indexes, it's a relatively independent but effective method to categorize subtypes of bladder cancer samples by using glucose metabolismrelated genes and it's worth further exploring.

Confirmation Of Glucose Metabolism-related Genes Associated With Survival
In the training set, by comparing patients with survival time less than 1 year and patients who lived longer than 3 years, 70 DEGs were screened out with p < 0.05 as the threshold (see Additional file 1).
Then, 101 genes significantly associated to patients' survival status were obtained from 606 glucose metabolism-related genes by using univariate COX regression analysis in the second step (see Additional file 2). The 46 overlapping genes from the two steps above were survival-related DEGs. In the third step, by analyzing promoter methylation levels of the 46 survival-related DEGs, we detected 18 genes of which the expression levels were significantly correlated to corresponding promotor methylation strengths, regarding the numbers of significant probes/the numbers of non-significant probes > 2/3 as the threshold (see Additional file 3). Thus, via the three steps above, a total of 18 target genes (FBP1, HIP1R, USF1, SLC45A3, SARS, LRRC59, TRAF5, PICK1, UCHL1, NUP188, COPS5, MGST2, SSTR5, ZAP70, SLC44A4, PYCR1, TXNIP and ELF3) associated with both patients' prognosis and DNA methylation status were selected as predictors for following further study (Fig. 3a). The regression coefficients (β) of the 18 genes were obtained from the univariate Cox proportional hazard regression analysis (Fig. 3b). Genes with regression coefficient values higher than 1 were identified as risky factors while those with coefficient values lower than 1 were protective factors. After calculating the risk score of each patient according to their 18-mRNA expression signature in the training set, we divided patients into a high-risk group (n = 150) and a low-risk group (n = 150) with the median risk score as the cutoff (Fig. 4a). Each patient's survival time and ending point are shown in (Fig. 4b), which reveals remarkable worse prognosis in the high-risk group than in the low-risk group. The 18 selected genes' expression levels of each patient are shown in (Fig. 4c) and we found that risky and protective genes showed opposite expression patterns according to the 18-mRNA risk score. Higher levels of 5 genes (UCHL1, PYCR1, SARS, NUP188, LRRC59) were more likely to be expressed in the high-risk group than in the low-risk group while the other genes (FBP1, HIP1R,USF1, SLC45A3, TRAF5, PICK1, COPS5, MGST2, SSTR5, ZAP70, SLC44A4, TXNIP, ELF3) were in the opposite situation, which was similar to our previous result (Fig. 3b).
We then carried out Kaplan-Meier analysis of overall survival for patients in the training set, which suggested great performance of the risk model in dividing patients into high-risk and low-risk (p = 0.0015) (Fig. 4d). In order to confirm the predictive ability of this risk formula, we calculated the 18-mRNA risk score of each patient in the validation set and stratified the patients into the same subgroups by using the median risk score as the cutoff. Similar to the result in the training set, the overall survival rate of the high-risk group was significantly lower than that of the low-risk group in the validation set (p = 0.00024) (Fig. 4e).
ROC analyses for 1-year, 3-year and 5-year survival were respectively carried out in the training set and the validation set, respectively. The results were shown in Fig. 5a-f, which exhibited satisfying survival predictive ability of the 18-mRNA risk score signature. Moreover, the sensitivity and specificity were even higher for the validation set than for the training set, which provides us with more confidence in the prognostic signature we built.
The 18-mRNA risk score signature was an independent prognostic factor for bladder cancer patients Multivariate Cox regression analysis was performed in the validation set to determine whether the prognostic signature we had previously built could function as an independent prognostic indicator.
When we included age, neoplasm grade, pathological T grade, pathological N grade, pathological M grade and diagnostic tumor stage together with the 18-mRNA risk score into the Multivariate Cox regression model, the age (P < 0.0001, HR = 1.033), pathological N stage (P = 0.022, HR = 1.195) and the 18-mRNA risk score signature (P = 0.004, HR = 1.167) showed significance while the other clinicopathological factors showed no statistical significance (Table. 2). Thus, it was indicated that the signature we established might be an independent prognostic factor for BLCA patients. GSEA analysis between high-risk patients and low-risk patients GSEA analysis result showed pathways and functions including wnt/beta-catenin-independent signaling pathway were significantly enriched in the high-risk group (Fig. 6a-c), which might provide us with insights into cancer progression and potential molecular mechanisms.

Weighted Gene Co-expression Network Analysis (wgcna)
With 24286 genes were adopted in WGCNA network construction, including mRNAs, miRNAs and lncRNAs, 34 modules in total were obtained. Among the 34 modules, 9 genes from the 18 glucose metabolism-related genes were identified as hub genes which were linked to 3 modules of coexpressive genes (Fig. 6d).

Functional characteristics of UCHL1 and PYCR1
In order to further validate our bioinformatic analyses, we tended to select two risky genes from the 18 genes to conduct cell functional experiments by using siRNAs to silence their expressions. Via the process of analyzing DNA methylation status of survival-related DEGs, we found the promotor methylation strengths of UCHL1 and PYCR1 were extraordinarily related to their own expressions. At the same time, no studies have so far reported the functions that these two genes play on bladder cancer cells. With great interest, we performed functional validations to determine the impacts that the two genes might have on bladder cancer cell lines, UM-UC-3 and T24. qRT-PCR results confirmed the knockdowns of UCHL1 and PYCR1 by siRNAs at mRNA level, respectively (Fig. 7A-B). Then, the proliferation ability, migration and invasiveness of BLCA cells were proved to be obviously inhibited after downregulating the expression of UCHL1 or of PYCR1 (Fig. 7C-D). Via performing western blotting, we noticed that silencing UCHL1 or PYVR1 resulted in downregulation of the expressions N-cadherin and CyclinD1 at protein level, which suggested that the two genes might promote epithelialmesenchymal transition (EMT)of bladder cancer cells and accelerate proliferation by regulating cell cycle (Fig. 7E). Meanwhile, no significant changes were detected in beta-catenin protein level, which was accordant to the GSEA result (Fig. 7E). In general, the results above provided us with not only experimental basis to partially validate our bioinformatic analyses but also with further insights to molecular functions and mechanisms for further research.

Discussion
Recently, several studies have applied gene expression signature to predict characteristics or prognosis for urinary bladder cancer. For example, a 20-gene signature strongly linked to bladder cancer development was established [28] and a 5-gene signature was developed for prediction of progression of T1G3 bladder cancer [29]; researchers previously reported a four-gene signature to directly predict OS of urinary bladder cancer patients [30]. However, a majority of these genes belong to various categories and owns different properties. On the other hand, it's necessary and urgent to explore new schemes for cancer prognosis with an easier access to gene abnormality detection in individuals, which will enhance the velocity and efficiency of translation from basic research to clinical application. A large number of studies have provided deep insights into molecular mechanisms of glucose metabolism alterations in cancers [8][9][10][11][12], by which we are more convinced of the critical role glucose metabolism plays in malignant tumors. Meanwhile, with abnormal alterations widely appearing in different cancers, DNA methylation has been recently analyzed by computational models to discover potential biomarkers for prediction of early diagnosis or prognosis of cancers [18][19][20].
Here, we selected genes significantly correlated to their promotor methylation status as well as to the OS of bladder cancer after conducting full-scale analyses of glucose metabolism-related genes.
In our study, from 606 glucose metabolism-related genes in 12 pathways, we first screened out 18 target genes, of which expressions were significantly correlated both to promotor methylation strengths and bladder patients' survival. Then, a signature based on the expressions of these 18 genes were constructed, which showed great survival predictive ability both in training and validation sets and was believed to act as a potential independent indicator for prognosis of bladder patients. subtypes of bladder cancer for its participation in glycolysis flux and gluconeogenesis [31]; MGST2 could be induced in bladder cancer cells by Bacillus Calmette-Guérin (BCG) internalization and thus it might act as a novel biomarker of the response to BCG immunotherapy for bladder cancer [32]; TXNIP negatively regulates bladder carcinogenesis by attenuating SDF-1-CXCR4-induced ERK activation [33]; ELF3 can protects bladder cancer cells from epithelial-mesenchymal transition [34]. These above are highly in accordance to our results, supporting the validity of our study strategy. As for the other The role of PYCR1 as an oncogene has been demonstrated in several other types of cancers via different molecular mechanisms including regulating p38 MAPK and NF-kB signaling pathway [35][36][37][38].
Besides, PYCR1 can be regulated by microRNAs in some cancers [37][38]. What appeals to us is that despite participation in glucose metabolism, PYCR1 also plays a role in proline metabolism [39] and thus it's worth more comprehensive insights into its molecular functions. Different from PYCR1, whether UCHL1 is an oncogene or a tumor-suppressor depends on types of cancer [40][41][42][43][44]. However, similar to our findings, an article indicated that prognosis of renal cell carcinoma was associated with alterations of UCHL1 promoter methylation [44]. Moreover, according to the GSEA results, the highrisk patients significantly enriched genes involved in wnt/beta-catenin-independent signaling pathway, PCP pathway and protein folding, which were demonstrated to be linked to a variety of cancers [45][46][47][48][49][50]. Via experiments, we proved that downregulation of UCHL1 or PYCR1 indeed resulted in no obvious changes in beta-catenin protein level. Three modules of co-expression genes aimed at 9 hub genes were identified by WGCNA and along with GSEA results and it's worth further and deeper explorations for molecular mechanism demonstration and clinical application. Generally speaking, with functional experiment validations of the two genes, we have more confidence not only in UCHL1 and PYCR1 as potential biomarkers but also in the signature we constructed as a prognosis predictor for bladder cancer patients.
In our study, the 18 survival-related genes were all characterized by glucose metabolism participation and were promotor methylation associated, which hasn't been reported before. Currently, most cases of bladder cancer are diagnosed at advanced stages so that patients fail to receive curative surgery.
Chemotherapy has been the main option for bladder cancer treatment [51], but the impact it has on overall survival is still limited [2][3][4]. Our study provides urologists and researchers with a potential new method for bladder cancer prognosis evaluation by DNA methylation monitoring with its advantage of stable existence in plasma; on the other hand, a rationale for treating bladder cancer patients by targeting metabolism is worth further consideration. However, a few limitations of our study should be emphasized. First, our bioinformatic analyses were limited within TCGA database and we only conducted functional experiments aimed at two risky genes that were outstandingly methylation-associated. Second, despite retrospective validations, our results need to be validated in further prospective investigations.
In summary, we developed a glucose metabolism gene-related signature that could potentially act as an independent prognosis predictor for bladder cancer patients. Further explorations into molecular mechanisms will reveal how these genes regulate glucose metabolism in cancer progression and how the expressions of the genes are modulated by promotor methylations, which can provide prospective and inspiration for future strategies for bladder cancer diagnosis, progression monitoring and even treatments.

Conclusions
In this study, we were convinced that glucose metabolism played a critical role in urinary bladder cancer. 18 glucose metabolism-related genes were filtered, which were significantly linked to bladder cancer patients' survival status and regulated by DNA methylation. A prognostic signature based on the expressions of the 18 genes was identified as a potential independent factor for predicting the prognosis of bladder cancer patients. Two risky genes from the 18 genes, UCHL1 and PYCR1, were found to be able to promote bladder cancer cells' proliferation, migration and evasion via cell functional experiments, which could further verify our bioinformatic analyses. In general, the 18 generelated prognostic signature has a potential to expand the vision of prognosis prediction and even treatment for urinary bladder cancer patients.         as the mean ± standard deviation (SD). d Proliferation of UM-UC-3 and T24 after respectively silencing the expressions of UCHL1 and PYCR1. Data were ana-lyzed by T-test and presented as the mean ± standard deviation (SD). e Beta-catenin, cyclin D1 and Ncadherin proteins of UM-UC-3 and T24 were detected after respectively silencing the expressions of UCHL1 and PYCR1. *P < 0.05; **P < 0.01; ***P < 0.001.