Identification of key pathways and genes influencing prognosis in bladder urothelial carcinoma

Abstract

Background

Genomic profiling can be used to identify the predictive effect of genomic subsets for determining prognosis in bladder urothelial carcinoma (BUC) after radical cystectomy. This study aimed to investigate potential gene and pathway markers associated with prognosis in BUC.

Methods

A microarray dataset of BUC was obtained from The Cancer Genome Atlas database. Differentially expressed genes (DEGs) were identified by DESeq of the R platform. Kaplan–Meier analysis was applied for prognostic markers. Key pathways and genes were identified using bioinformatics tools, such as gene set enrichment analysis, gene ontology, the Kyoto Encyclopedia of Genes and Genomes, gene multiple association network integration algorithm (GeneMANIA), Search Tool for the Retrieval of Interacting Genes/Proteins, and Molecular Complex Detection.

Results

A comparative gene set enrichment analysis of tumor and adjacent normal tissues suggested BUC tumorigenesis resulted mainly from enrichment of cell cycle and DNA damage and repair-related biological processes and pathways, including TP53 and mitotic recombination. Two hundred and fifty-six genes were identified as potential prognosis-related DEGs. Gene ontology and Kyoto Encyclopedia of Genes and Genomes analyses showed that the potential prognosis-related DEGs were enriched in angiogenesis, including the cyclic adenosine monophosphate biosynthetic process, cyclic guanosine monophosphate-protein kinase G, mitogen-activated protein kinase, Rap1, and phosphoinositide-3-kinase-AKT signaling pathway. Nine hub genes, TAGLN, ACTA2, MYH11, CALD1, MYLK, GEM, PRELP, TPM2, and OGN, were identified from the intersection of protein–protein interaction and GeneMANIA networks. Module analysis of protein–protein interaction and GeneMANIA networks mainly showed enrichment of the cyclic guanosine monophosphate-protein kinase G signaling pathway, angiogenesis, cell proliferation, and differentiation, which are associated with tumor angiogenesis and cancer prognosis.

Conclusion

Genes and pathways related to cell cycle and DNA damage and repair may play a crucial role in BUC pathogenesis, whereas those pertaining to tumor angiogenesis may be key factors in influencing BUC prognosis, especially in advanced disease stages.

Keywords:

• bioinformatics analytical tools
• bladder urothelial carcinoma
• microarray
• differentially expressed gene
• prognosis

Introduction

Bladder urothelial carcinoma (BUC) is most prevalent in developed countries, with highest incidence rates in Europe, and occurs predominantly in males. In 2012, the worldwide incidence of BUC was ~429,800, with 165,100 BUC-associated deaths.¹ It was estimated that in 2015, China alone accounted for ~80,500 new BUC cases and 32,900 BUC-associated deaths. From 2000 to 2011, BUC incidence and mortality rates have increased markedly in the Chinese male population.² At least two types of BUC exist: low-grade nonmuscle invasive bladder cancers (NMIBC) and high-grade muscle-invasive bladder cancers (MIBC). Prognosis of 5-year survival in NMIBC and MIBC are ~90% and <50%, respectively, with progression, most often, to metastasis in the latter.^3,4 Molecular characteristics of MIBC and NMIBC are highly distinctive, as revealed in recent studies.^5–7

Genetic predisposition and environmental exposure are major risk factors that possibly affect the prognosis of BUC, as with other cancers. Smoking is a proven major risk factor for BUC and is associated with unfavorable outcomes.^8,9 Earlier genome-wide association studies have identified genetic variations at 8q24.21, 3q28, 8q24.3, 4p16.3, 22q13.1, 19q12, 2q37.1, and 5p15.33 as novel susceptibility loci associated with BUC.^10–12 Furthermore, the glutathione S-transferase M1 null and N-acetyltransferase-2 slow acetylator genotypes are strongly associated with BUC.^13,14

Genomic profiling can be used to identify the predictive effect of genomic subsets for BUC prognosis after radical cystectomy.¹⁵ A study in Chinese patients with BUC showed comparatively favorable cancer prognosis – >60% age-standardized 5-year survival rate.¹⁶ Given the high prevalence of smoking and the upward trends in BUC incidence and mortality rates, research into the risk and prognostic factors of BUC is expediently required.

Using bioinformatics, this study was conducted to map the gene expression profile of BUC patients with corresponding survival profiles from The Cancer Genome Atlas (TCGA; https://tcga-data.nci.nih.gov/; accessed October 20, 2016) to investigate potential key genes and pathways that influence BUC prognosis.

Materials and methods

Microarray data

A raw microarray dataset of BUC including human BUC mRNA expression and corresponding survival profiles was obtained from TCGA. There were 408 patients with BUC and 414 tumor tissues and 19 adjacent normal tissues microarray data in the TCGA dataset; information on overall survival (OS) as well as the status of events was available for 406 of these patients. A total of 406 BUC patients from the TCGA with complete follow-up profiles were analyzed in the survival study.

Identification of DEGs

DESeq, an R package for transcriptome chip-based expression profile analysis,¹⁷ was used according to the manufacturer’s instructions to test for differential expression of RNA transcript levels between BUC and adjacent normal tissues. A volcano plot was drawn using the gplots package. DEGs were identified with a | log₂ fold change (FC) | ≥2; an adjusted P-value cutoff of <0.05 was considered indicative of statistically significant differences.

Gene set enrichment analysis (GSEA)

Differences in gene expression levels of biological pathways in cancer and adjacent normal tissues were analyzed by GSEA v2-2.2.2 (http://software.broadinstitute.org/gsea/index.jsp; accessed October 23, 2016), with reference gene sets from the Molecular Signatures Database (MSigDB) of c2 (Kyoto Encyclopedia of Genes and Genomes [KEGG] gene sets: c2.cp.kegg.v5.2.symbols.gmt) and c5 (gene ontology [GO] gene sets: c5.bp.v5.2.symbols.gmt; c5.cc.v5.2.symbols.gmt; and c5.mf.v5.2.symbols.gmt).¹⁸ The number of permutations was set at 1,000. Enrichment results satisfying a nominal significance cutoff of <0.05 with an FDR <0.25 were considered statistically significant.

Identification of potential prognosis-related gene markers

Survival analyses were conducted on 406 patients with normalized mRNA expression and OS profiles. Patients were classified into two subgroups by gene expression levels: a high gene expression group (expression levels greater than the median value) and a low gene expression group (expression levels equal to or less than the median value). Kaplan–Meier analysis with the log-rank test was applied to estimate survival in the study sample by using a survival package in the R platform. A log-rank P-value cutoff of <0.05 was considered statistically significant for identifying potential prognosis-related gene markers.

Functional annotation and pathway enrichment analysis of potential prognosis-related DEGs

Potential prognosis-related DEGs were defined as the intersection of DEGs and the potential prognosis-related gene markers. The Database for Annotation, Visualization, and Integrated Discovery (DAVID) v.6.8 (https://david.ncifcrf.gov/tools.jsp; accessed October 23, 2016) was used to annotate input genes, classify gene functions, identify gene conversions, and carry out GO term analysis.¹⁹ To identify the functional annotation of prognosis-related DEGs, we analyzed GO terms, and KEGG pathway enrichment with DAVID, while specifying an enrichment P-value <0.05 for statistical significance.

Gene and protein networks and module analysis

GeneMANIA (http://genemania.org/; accessed October 25, 2016), a real-time multiple association network integration algorithm for predicting gene function,²⁰ was used for analyzing gene–gene interactions in the study. The Search Tool for the Retrieval of Interacting Genes/Proteins (STRING) database (http://string.embl.de/; accessed October 25, 2016), which provides a critical assessment and integration of protein–protein interactions (PPI),²¹ was used to assess direct (physical) and indirect (functional) associations of potential prognosis-related DEGs. Interactions of potential prognosis-related DEGs with a combined score >0.4 were considered statistically significant. Both gene–gene and protein–protein interactions were used for module screening by Molecular Complex Detection (MCODE; scores >3 and nodes >4) in Cytoscape, a bioinformatics integration platform. Moreover, we analyzed the functional annotation and pathway for potential prognosis-related DEGs in the modules.

Statistical method

The Kaplan–Meier method with a log-rank test was used to calculate clinical outcomes between different gene expression groups. Hazard ratios (HRs) and their 95% confidence intervals (CIs) were calculated from the Cox proportional hazards regression model. The adjusted P-value and FDR in DESeq and GSEA, respectively, were adjusted for multiple testing with the Benjamini–Hochberg procedure to control false discovery rate.^22–24 A value of P<0.05 was considered statistically significant. All statistical analyses were conducted with SPSS 20.0 (IBM Corporation, Armonk, NY, USA) and R 3.3.0.

Results

Screening of DEGs

In total, 414 BUC tumor and 19 adjacent normal specimens were screened for DEGs, and the microarray datasets of these samples were compared using DESeq. Specified criteria, |log₂ FC| ≥2, and P-value <0.05 were met for 766 genes, including 667 and 99 genes that were down- and upregulated, respectively (Figure 1, Table S1).

Figure 1 Volcano plot of the differentially expressed genes.

Notes: Red: upregulation; green: downregulation; black: non-differentially expressed genes.
Abbreviation: P-adj, adjusted P-value.

GSEA

Comparative gene expression studies of BUC tumor and adjacent normal samples in the TCGA dataset was done with GSEA (Figure 2, Table S2). In the GSEA of the KEGG pathway, tumor tissue gene expression was associated with the cell cycle, base excision repair, and the tumor protein p53 (TP53) signaling pathway. In GO enrichment analysis, biological processes pertaining to mitotic recombination, sister chromatid segregation, and mitotic sister chromatid segregation were significantly enriched, whereas molecular function in deoxyribonuclease, endodeoxyribonuclease, DNA-dependent ATPase, and DNA N-glycosylase activities was found to be enriched. Furthermore, the cellular component was enriched for condensed chromosome, sex chromosome, and condensed chromosomal centromeric region.

Figure 2 GSEA results of tumor and adjacent normal tissue samples.

Note: GSEA result of bladder urothelial carcinoma tumor tissue using c2 (A, B) and c5 (C–H) reference gene sets.
Abbreviations: ES, enrichment score; FDR, false discovery rate; GSEA, gene set enrichment analysis.

Screening of potential prognosis-related gene markers

Patients (n=406) with complete follow-up data were included in further survival analyses. Using Kaplan–Meier analysis with a log-rank test P-value <0.05, we isolated 5,019 genes as potential prognosis-related gene markers of BUC (Table S3).

GO term and KEGG pathway enrichment analysis

Potential prognosis-related DEGs were defined as the intersection of DEGs and potential prognosis-related gene markers. Two hundred and fifty-six genes that met this criterion were defined as potential prognosis-related DEGs (Table S4), and subjected to further GO and KEGG pathway analyses with DAVID. Heat maps of potential prognosis-related DEGs are shown in Figure 3. The GO analysis suggested (Table 1) that potential prognosis-related DEGs are significantly enriched in processes such as cell adhesion; angiogenesis; and negative regulation of cyclic adenosine monophosphate (cAMP) biosynthesis, protein phosphorylation, and apoptotic processes; however, with regard to molecular function, DEGs were significantly enriched in heparin binding, 3′,5′-cyclic nucleotide phosphodiesterase activity, and calcium ion binding. Further, GO cell component analysis showed that potential prognosis-related DEGs were enriched in the extracellular matrix and plasma membrane.

Figure 3 Heat map of the 256 prognosis-related differentially expressed genes.

Notes: Red: upregulation; green: downregulation.

Table 1 GO term analysis of potential prognosis-DEGs associated with BUC prognosis

Category	GO ID	Term	Count	%	P-value
GOTERM_BP_DIRECT	GO:0007155	Cell adhesion	22	5.64	9.15E-07
GOTERM_BP_DIRECT	GO:0001525	Angiogenesis	12	3.07	2.16E-04
GOTERM_BP_DIRECT	GO:0030818	Negative regulation of cAMP biosynthetic process	4	1.02	0.00114
GOTERM_BP_DIRECT	GO:0001933	Negative regulation of protein phosphorylation	6	1.54	0.001286
GOTERM_BP_DIRECT	GO:0060548	Negative regulation of cell death	5	1.28	0.006908
GOTERM_BP_DIRECT	GO:0045766	Positive regulation of angiogenesis	7	1.79	0.004536
GOTERM_MF_DIRECT	GO:0008201	Heparin binding	18	4.61	3.44E-11
GOTERM_MF_DIRECT	GO:0005509	Calcium ion binding	24	6.15	5.92E-05
GOTERM_MF_DIRECT	GO:0004114	3′,5′-Cyclic nucleotide phosphodiesterase activity	5	1.28	1.66E-04
GOTERM_MF_DIRECT	GO:0005516	Calmodulin binding	10	2.56	8.05E-04
GOTERM_MF_DIRECT	GO:0030553	cGMP binding	4	1.02	0.001269
GOTERM_CC_DIRECT	GO:0031012	Extracellular matrix	23	5.89	7.06E-11
GOTERM_CC_DIRECT	GO:0005886	Plasma membrane	89	22.80	2.59E-07
GOTERM_CC_DIRECT	GO:0005925	Focal adhesion	19	4.87	4.79E-06
GOTERM_CC_DIRECT	GO:0009986	Cell surface	22	5.64	1.57E-05
GOTERM_CC_DIRECT	GO:0005576	Extracellular region	40	10.25	1.19E-04

Abbreviations: GO, gene ontology; DEGs, differentially expressed genes; BUC, bladder urothelial carcinoma; cAMP, cyclic adenosine monophosphate; cGMP, cyclic guanosine monophosphate.

On KEGG pathway analysis, potential prognosis-related DEGs were significantly enriched in cyclic guanosine monophosphate-protein kinase G (cGMP-PKG), calcium, cAMP, cancer, and mitogen-activated protein kinases (MAPKs) signaling pathways that were associated with cancer prognosis (Table 2).

Table 2 KEGG analysis of potential prognosis-DEGs associated with BUC prognosis

Pathway ID	Name	Count	%	P-value	Genes
hsa04022	cGMP-PKG signaling pathway	14	3.59	3.09E-07	KCNMA1, ADCY5, MRVI1, PRKG1, EDNRA, ADRB3, AGTR1, ATP2B4, PDE2A, RGS2, PLN, PDE5A, MYLK, ADRA1D
hsa04020	Calcium signaling pathway	13	3.33	4.74E-06	PTGER3, PTGFR, EDNRA, AGTR1, ADRB3, ATP2B4, CHRM2, PDE1C, PLN, PDE1A, CAMK2A, MYLK, ADRA1D
hsa04024	cAMP signaling pathway	10	2.56	0.001408	EDNRA, PTGER3, ATP2B4, RAC3, CHRM2, ADCY5, JUN, PLN, PDE4B, CAMK2A
hsa05200	Pathways in cancer	13	3.33	0.006683	LAMA2, EDNRA, AGTR1, FGFR1, VEGFD, FGF7, PTGER3, RAC3, JUN, ADCY5, RUNX1T1, FGF10, FGF2
hsa04010	MAPK signaling pathway	10	2.56	0.007592	FGFR1, CACNA2D1, DUSP2, FGF7, RAC3, JUN, FGF10, FLNC, FGF2, FLNA
hsa04514	Cell adhesion molecules	7	1.79	0.012562	NCAM1, CD34, CNTN1, NRXN1, JAM2, NEGR1, JAM3
hsa05205	Proteoglycans in cancer	8	2.05	0.018787	FGFR1, HSPG2, DCN, THBS1, FLNC, CAMK2A, FGF2, FLNA
hsa04015	Rap1 signaling pathway	8	2.05	0.023834	FGFR1, VEGFD, FGF7, RAC3, ADCY5, FGF10, THBS1, FGF2
hsa04151	PI3K-AKT signaling pathway	10	2.56	0.044252	LAMA2, FGFR1, VEGFD, FGF7, CHRM2, ITGA7, FGF10, THBS1, FGF2, GHR

Abbreviations: KEGG, Kyoto Encyclopedia of Genes and Genomes; DEGs, differentially expressed genes; BUC, bladder urothelial carcinoma; cGMP, cyclic guanosine monophosphate; PKG, protein kinase G; cAMP, cyclic adenosine monophosphate; MAPK, mitogen-activated protein kinases; Rap1, Ras-associated protein 1; PI3K-AKT, phosphoinositide-3-kinase-AKT.

Gene and protein interactome networks

To investigate the interaction and hub genes of potential prognosis-related DEGs, gene–gene and protein–protein interactomes were constructed using GeneMANIA and STRING, respectively (Figures 4 and 5). The top 30-degree gene nodes at the intersection of two interactome networks were called hub genes. Nine genes were defined as hub genes: transgelin (TAGLN), alpha-smooth muscle actin (ACTA2), myosin heavy chain 11 (MYH11), caldesmon 1 (CALD1), myosin light chain kinase (MYLK), GTP-binding protein overexpressed in skeletal muscle (GEM), proline and arginine rich end leucine-rich repeat protein (PRELP), tropomyosin 2 (TPM2), and osteoglycin (OGN). GeneMANIA (Figure 4) and PPI analyses (Figure 5) showed gene nodes with the highest degree were TAGLN with 373 and ACTA2 with 21 connections. Module genes of GeneMANIA (Figure 6) and PPI analyses (Figure 7) were identified by MCODE in Cytoscape. Functional annotation and pathway analyses, analysis by DAVID, showed that module genes were related mainly to the cGMP-PKG signaling pathway, angiogenesis, regulation cell differentiation, and proliferation, all of which were associated with tumor biology and affect survival and prognosis in patients with BUC.

Figure 4 GeneMANIA interaction networks.

Abbreviation: GeneMANIA, gene multiple association network integration algorithm.

Figure 5 Protein–protein interaction networks.

Notes: Black text: prognosis-related differentially expressed genes; Red text: nine hub genes of the prognosis-related differentially expressed genes in the protein–protein interaction networks.

Figure 6 Top three modules from the GeneMANIA interaction networks.

Notes: A: Module 1; B: the enriched GO term and pathways of module 1; C: module 2; D: the enriched GO term and pathways of module 2; E: module 3; F: the enriched GO term and pathways of module 3.

Abbreviations: GeneMANIA, gene multiple association network integration algorithm; GO, gene ontology; cGMP-PKG, cyclic guanosine monophosphate-protein kinase G.

Figure 7 Top 3 modules from the PPI interaction networks.

Notes: A: module 1; B: the enriched GO term and pathways of module 1; C: module 2; D: the enriched GO term and pathways of module 2; E: module 3; F: the enriched GO term and pathways of module 3.

Abbreviations: PPI, protein–protein interaction; GO, gene ontology; cGMP-PKG, cyclic guanosine monophosphate-protein kinase G.

Stratified and joint-effects analysis between hub genes and advanced tumor stage on OS in BUC

Associations between the expression of nine hub genes and BUC survival were further evaluated in stratified and joint-effects analysis by tumor stage. The majority of BUC patients in TCGA were in advanced tumor stages, and only two patients had stage I disease. Therefore, patients with stage I disease and those without a tumor-staging report were excluded from further stratified and joint-effects analysis. A total of 402 BUC patients with advanced tumor stage (stages II, III, and IV) were included into stratified and joint-effects analysis. The outcomes in different tumor stages were significantly different (Figure 8); a higher stage was associated with poor OS. Further, we observed that the level of mRNA expression for hub genes significantly differed among different tumor stages, with higher stages seeming to have a higher mRNA expression (Figure S1). Stratified analyses are shown in Table S5. We observed that high GEM (P=0.021, HR =1.732, 95% CI =1.085–2.766) and OGN (P=0.048, HR =1.621, 95% CI =1.004–2.618) expressions were significantly associated with an increased risk of death in BUC stage IV. Similar results were observed with high OGN (P=0.05, HR =1.765, 95% CI =1.0–3.114, Table S5) in BUC stage III. However, significant results were not found for other hub genes. Joint-effects analysis showed that there was a statistically significant interaction between nine hub genes and tumor stage (Table S5, all interaction P-values <0.001). Compared with the high gene expression with BUC stage IV, patients with other stages and different gene expression levels had a significantly lower risk of death.

Figure 8 Kaplan–Meier survival curves for BUC patients.

Note: OS stratified by tumor stage.
Abbreviations: BUC, bladder urothelial carcinoma; OS, overall survival.

Discussion

The etiology of BUC is attributed to gene–environment interactions. In the multifactorial etiopathogenesis of BUC, smoking has been proven to play an important role as an environmental modulator of the risk of BUC and associated mortality. However, the genetic predictors of prognosis in BUC are unknown. The present study compared gene expression profiles of BUC tumor and adjacent normal tissues, from the TCGA, to identify DEGs and, thereafter, used GSEA to investigate and interpret differences of gene expression. GSEA results of KEGG suggested that differences of gene expression between tumor and adjacent normal tissues were greater in cell cycle and tumor protein p53 (TP53) pathways. Our results are consistent with three earlier studies that elicited the predictive value of cell cycle regulatory protein expression in the prognosis and progression of BUC.^25–27 Alterations of TP53 are predictive for BUC recurrence and are markedly associated with an unfavorable prognosis after radical cystectomy.^28,29 GO enrichment analysis for biological processes demonstrated that BUC tumors were associated with mitotic recombination, sister chromatid segregation, and mitotic sister chromatid segregation, all of which are related to cell cycle regulation.^30–32 Our results, moreover, are supported by the study of Guo et al,³³ which identified that genes involved in sister chromatid segregation were frequently altered in bladder cancer via whole-genome and whole-exome sequencing. However, in BUC, molecular function was enriched in deoxyribonuclease, endodeoxyribonuclease, DNA-dependent ATPase, and DNA N-glycosylase activities, whereas the cellular component was enriched in condensed chromosome, sex chromosome, and condensed chromosome centromeric region. Li et al³⁴ showed that DNA-dependent ATPase activity was associated with DNA damage, whereas Wu et al³⁵ reported the involvement of DNA N-glycosylase activity in DNA repair. Further, deoxyribonuclease and endodeoxyribonuclease activities caused DNA fragmentation in nuclei and, thereby, DNA damage. In the present study, the enrichments in molecular function pertaining to DNA damage, taken together with the GSEA enrichment of KEGG, indicate that TP53 – corresponding to DNA repair and damage – may play a vital role in BUC tumorigenesis. Therefore, cell cycle and DNA damage may potentially play a major role in bladder carcinogenesis and may be associated with genetic risk factors.

On univariate analysis of genome-wide survival, 5,019 genes were found to be significantly associated with BUC prognosis. Then, 256 genes, overlapping in DEGs and potential prognosis-related genes, qualified as potential prognosis-related DEGs. Both KEGG and GO analyses of potential prognosis-related DEGs were undertaken to identify pathways and genes associated with BUC survival. GO analysis showed that potential prognosis-related DEGs were enriched in angiogenesis, cAMP biosynthetic processes, cell apoptosis, cGMP binding, etc. Studies have implicated tumor angiogenesis and angiogenesis-related molecular markers as predictors of BUC prognosis.^36–39 Moreover, a previous study among TCGA MIBC patients revealed that genes from the ErbB family that were related to angiogenesis were frequently altered.⁴⁰ Another large-scale transcriptomic data of MIBC patients suggested that basal-like MIBC activated the epidermal growth factor receptor (EGFR) pathway, connected with frequent EGFR gains and EGFR autocrine loop activation; in addition, they demonstrated that basal-like MIBC cell lines were sensitive to anti-EGFR therapy.⁴¹ Caspase-3, an active apoptosis biomarker, is correlated with a favorable prognosis, whereas modulation of multiple apoptosis biomarkers synergistically predicts BUC recurrence and mortality.^42,43 Earlier studies have reported that the cAMP-PKA-cAMP response element-binding (CREB) signaling pathway modulates vascular endothelial growth factor (VEGF), whereas loss of the dimethylarginine dimethylaminohydrolase 1 (DDAH1) effect on the NO–cGMP–PKG pathway results in decreased angiogenesis.^44,45

KEGG enrichment analysis showed that potential prognosis-related DEGs were enriched in cGMP-PKG, cAMP, MAPK, Ras-associated protein 1 (Rap1), PI3K-AKT signaling pathways, and so on. As mentioned earlier, the cAMP and cGMP signal pathways are associated with angiogenesis. In addition, long noncoding RNA urothelial carcinoma associated 1-regulated cell cycle via CREB protein and the PI3K-AKT-dependent pathway in bladder cancer, wherein CREB plays a tumor suppressor role in BUC.^46,47 Sheta et al⁴⁸ reported that the MAPK signaling pathway mediates the regulation of tumor VEGF expression. Rap1 is a small GTPase with a role in cell adhesion and is a positive regulator of angiogenesis. A recent study reported that activation of the MAPK signaling pathway will reduce the response to VEGF stimulation in Rap1b-deficient endothelial cells.⁴⁹ Rap1 plays a key role in human angiogenesis by promoting integrin, VEGFR2, and FGF2 activation, and is a prerequisite for angiogenesis.^50–52 Several studies have shown VEGF to be regulated by PI3K-AKT signaling, which promotes angiogenesis. Toll-like receptor 4 induces VEGF expression via activation of PI3K-AKT signaling, to enhance angiogenesis in pancreatic cancer.⁵³ Similarly, a study of colorectal cancer showed that arginine ADP-ribosyltransferase 1 promotes expression of hypoxia-inducible factor 1-α via the PI3K-AKT signaling pathway, whereas upregulation of VEGF and basic fibroblast growth factor promotes cancer angiogenesis.⁵⁴ TCGA researchers identified potential therapeutic targets in BUC, with 45% of targets in the receptor tyrosine kinase–mitogen-activated protein kinase (RTK–MAPK) pathway.⁶ Our results are consistent with those reported, which suggested the MAPK pathway may play an important role in BUC prognosis. However, inhibition of p38 MAPK decreased BUC proliferation, growth, and cell invasiveness of bladder cancer via the MAPK pathway.^55,56 Previous studies have demonstrated that activation of the PI3K-AKT–mechanistic target of rapamycin (mTOR) pathway promotes bladder carcinogenesis and, further, induces tumor growth and chemoresistance by nicotine.^57,58 Moreover, activation of the PI3K-AKT–mTOR pathway is correlated with advanced tumor progression and unfavorable survival outcomes.⁵⁹ A previous study in the TCGA MIBC patient population revealed that the PI3K-mTOR signaling pathway was commonly altered and may serve as a potential therapeutic target.⁴⁰ These studies complement the results of the GO and KEGG pathway enrichment analyses that focus mainly on angiogenesis-related genes and pathways. Therefore, as angiogenesis-related genes and pathways may play a crucial role in BUC prognosis, they can be potential therapeutic targets.

Furthermore, we constructed the gene–protein interactome networks with potential prognosis-related DEGs and identified nine hub genes, including TAGLN, ACTA2, MYH11, CALD1, MYLK, GEM, PRELP, TPM2, and OGN, with TAGLN and ACTA2 exhibiting the highest degree of connectivity in gene and protein interactome networks, respectively. TAGLN is a member of the calponin family of actin-binding proteins that participates in cell motility and migration; TAGLN overexpression may lead to poor clinical outcome as it is an independent prognostic factor in oral squamous cell carcinoma, lung adenocarcinoma, and pancreatic cancer.^60–62 Moreover, TAGLN promotes tumor progression.^61,62 ACTA2 belongs to the actin family of proteins and plays a role in cell motility, structure, and integrity. Lee et al⁶³ reported that ACTA2 regulates c-MET and FAK expression, then positively influences metastatic potential of lung adenocarcinoma, and affects prognosis. CALD1 encodes a calmodulin- and actin-binding protein that regulates smooth muscle and nonmuscle contraction. Lee et al⁶⁴ reported that CALD1 promotes cell migration and invasiveness of bladder cancer, and its overexpression is significantly associated with an unfavorable prognosis. MYLK is a member of the immunoglobulin gene superfamily, and has been linked to the proliferative ability of breast cancer via extracellular signal-regulated kinase (ERK1/2) and the P38 pathway; moreover, it is associated with modulation of tumor invasiveness and metastasis in breast cancer.^65–67 TPM2 encodes beta-tropomyosin, a member of the actin filament binding protein family that is poorly expressed in high-grade, relapsed, and metastatic prostate tumors and may be a potential prognostic biomarker in prostate cancer.⁶⁸ MYH11 encodes MYH11 protein – a smooth muscle myosin protein belonging to the myosin heavy chain family – associated with the composition of the oncogenic fusion gene CBFB/MYH11 in acute myeloid leukemia; a mutation of MYH11 is implicated in human intestinal tumorigenesis.^69,70 OGN, also known as the osteoinductive factor, plays a potential role in the development of ovarian cancer and bone metastasis, as identified by bioinformatics analysis. The GEM protein belongs to the RAD/GEM family of GTP-binding proteins and is a regulatory protein in receptor-mediated signal transduction, whereas PRELP is a leucine-rich repeat protein present in the extracellular matrix in connective tissue.⁷¹ However, studies of the GEM and PRELP genes have not elicited any implications for cancer, and their biological role in BUC prognosis is unclear. Therefore, further studies are needed to determine the role of GEM and PRELP in BUC prognosis.

Module analysis of the GeneMANIA network indicated that BUC prognosis refers to muscular contraction, the cGMP-PKG signaling pathway, extracellular matrix, extracellular matrix organization, angiogenesis, cell proliferation, and positive regulation of cell differentiation, whereas analysis of the PPI network module revealed an association of BUC prognosis with cGMP binding, cGMP-PKG signaling pathway, actin binding, positive regulation of GTPase activity, positive regulation of cell differentiation, and pathways in cancer. The cGMP-PKG signaling pathway is associated with angiogenesis, and tumor angiogenesis modulates prognosis in cancer. Smith et al⁷² revealed that the Ral GTPase pathway is essential for key phenotypes in bladder cancer progression models as well as for regulation of the expression of key molecules such as prognostic markers. In addition, downregulation of DOC-2⁄DAB2 interactive protein, another Ras GTPase-activating protein family member, could promote cell proliferation, migration, and invasion, and predicted poor disease-free survival and overall survival in patients with BUC.⁷³ The degree of tumor differentiation determines the biological behavior of the tumor and also affects its prognosis. Antunes et al⁷⁴ reported that squamous differentiation was an independent prognostic factor for cancer-specific mortality in patients with bladder cancer who underwent radical cystectomy. Similarly, several studies have demonstrated that increased tumor proliferation is significantly associated with invasion and recurrence, metastasis, and poor prognosis. Moreover, determination of proliferation-related molecular markers could improve the prediction of recurrence and cancer-specific mortality of BUC following radical cystectomy.⁷⁵

In BUC patients with advanced tumor stages, our study demonstrated that the mRNA level of hub genes had a significant association with advanced tumor stage, and showed an upward trend in accordance with tumor stage progression. Moreover, we observed strong interactions between mRNA expression of hub genes and advanced tumor stage on prognosis prediction, which suggests that there might be potential modulating roles of tumor progression in the biological functions of these hub genes. In invasive BUC, the most important prognostic factor is tumor stage, which is based on the depth of tumor invasion and metastasis.⁷⁶ Further, the survival analysis in our study suggests that advanced tumor stage is associated with poor prognosis. These evidences may help establish a potential strategy for advanced BUC survival prediction and facilitate individualized therapeutic strategies.

This study had certain limitations in that the clinical information from the TCGA database was not comprehensive; therefore, multivariate survival analysis for confounding factors affecting prognosis in patients with BUC was not used. As a result, potential prognosis-related genes were screened using the Kaplan–Meier analysis in the current study.

Conclusion

On the basis of bioinformatics analysis of potential prognosis-related DEGs, we identified biomarkers that are related to the prognosis of BUC. Results from the present study provide a cluster of potential prognosis-related genes and pathways for future investigation. Furthermore, comparative analysis of tumor and adjacent normal tissues via GSEA revealed that cell cycle, together with DNA damage and repair-related genes and pathways, may play a crucial role in BUC pathogenesis, especially for advanced stage BUC. Further, bioinformatics analysis of potential prognosis-related DEGs indicated that tumor angiogenesis-related genes and pathways play a key role in prognosis. These results would promote an understanding of advanced BUC pathogenesis, and provide a number of valuable potential genes and pathways for further investigation of prognostic markers in advanced BUC. However, further research into molecular biology of BUC is required to confirm these results.

Acknowledgments

The authors thank the contributors of the Cancer Genome Atlas for sharing the BUC transcriptome chip data on open access.

Disclosure

The authors report no conflicts of interest in this work.