Bladder cancer (BC) is the second most frequent genitourinary malignancy and the sixth most common malignancy in men. BC represents a spectrum of disease ranging from superficial, well-differentiated disease, which does not significantly affect the survival of BC patients, to highly fatal tumors for which long-term survival may be dismal (So, 2016; Ghervan et al., 2017). With the aging of population, the incidence of BC is rising year by year, and BC in elder patients will become even more frequent and evolved into a public health challenge in future. For patients with superficial BC, telescopic removal of the cancer (transurethral resection of bladder tumor, TURBT) followed by instillation of chemotherapy or vaccine-based therapy into the bladder with prolonged telescopic checking of the bladder are usually recommended, and the 5-year overall survival for these patients reaches 90%, while about 40–80% of these patients will develop disease recurrence or progression (Malmström et al., 2017). For patients with invasive BC, radical cystectomy plus pelvic lymph node dissection (PLND) followed by neo-adjuvant chemotherapy is recommended as a standard of care, and once it becomes metastatic cancer, the 5-year overall survival for patients with invasive BC is a dismal 6% (Salama et al., 2016; Sargos et al., 2016).

Biomarkers are biological substances whose detection indicates a particular disease state. So far, a variety of biomarkers have been introduced in daily clinical practice, including risk assessment, screening, differential diagnosis, determination of prognosis, prediction of response to treatment, and monitoring of progression of disease. Thus, identification of biomarkers that are associated with clinical outcomes of patients with BC might be of clinical significance (Giunchi et al., 2017).

Currently, several screening algorithms based on gene expression data, including Gene Set Enrichment Analysis (GSEA), Signalling Pathway Impact Analysis (SPIA), Topology Gene-Set Analysis, and DEGraph, and in silico Pathway Activation Network Decomposition Analysis (iPANDA), etc., have been introduced for both academic and commercial purpose, and most of these algorithms are intended to identify differentially expressed genes between groups of samples (Ozerov et al., 2016). Weighted gene co-expression network analysis (WGCNA), a systems biology algorithm, can be applied to describing the correlation patterns among genes across microarray samples, finding and summarizing modules of high related genes, and relating modules to certain clinical phenotype (Zhang and Horvath, 2005; Yip and Horvath, 2007). During the network construction, highly co-expressed genes are connected in the network and grouped into modules, thus, different modules included different functionally related genes and the most central and connected genes are treated as hub genes (Zhang and Horvath, 2005; Yip and Horvath, 2007). Correlation networks facilitate network based gene screening methods that can be used to identify candidate biomarkers or therapeutic targets. Therefore, WGCNA is usually used to facilitate the screening or identification of candidate biomarkers or therapeutic targets. Tumor node metastasis (TNM) 2009 (7th edition) was recommended for the BC, and previous studies demonstrated that patients with an infiltrative pattern had better survival than those with other pattern types (PDQ Adult Treatment Editorial Board, 2002; Yaxley, 2016). In the present study, the WGCNA algorithm was applied to identify candidate biomarkers for BC based on the TNM staging of BC patients.

As mentioned above, BC is one of the most common cancers and significant progress has been made during the past decade. The prognosis of patients with BC, especially invasive BC, remains poor. The recurrence rate for superficial BC is about two-thirds, and, despite advances in management, some patients still develop stage progression. Meanwhile, approximately 30% of muscle-invasive BCs are associated with occult distant metastasis at the time of diagnosis, which led to a dismal 5-year survival of patient with muscle invasive BC (Kamat et al., 2016). Biomarker is defined as a substance found in tissue, blood, or other body fluids that might be a sign of cancer or noncancerous conditions (Mohammed et al., 2016). Previous studies suggested that biomarkers, including BC biomarkers, exhibited several features: prognostic, predictive, and pharmacodynamic. The biological functions of prognostic biomarker, predictive biomarker, and pharmacodynamic biomarker were predicting the natural course of cancers, evaluating the probable benefit of a particular treatment, and assessing the treatment effects of a drug on a tumor and determine the proper dosage of a new anticancer drug, respectively (Kojima et al., 2016; Mohammed et al., 2016; Nandagopal and Sonpavde, 2016). Nowadays, many biomarkers and their corresponding targeted agents have been determined for the diagnosis and treatment of BC. In the present study, we identified 50 hub genes that were significantly correlated with the TNM staging of BC patients using WGCNA.

As we know, TNM staging, devised by Pierre Denoix, describes the size of the original (primary) tumor and whether it has invaded nearby tissue (T), describes nearby (regional) lymph nodes that are involved (N), and describes distant metastasis (M) (Denoix, 1946). The functional enrichment analysis of hub genes that were correlated with the TNM stages indicated that these genes were enriched in cell proliferation associated GO terms (mitotic sister chromatid segregation, mitotic cell cycle and cell cycle, etc.) and KEGG pathways (cell cycle, Oocyte meiosis, etc.). Thus, we speculated that these hub genes affected the TNM staging of BC patients through promoting the proliferation of BC cells.

The PPI network analysis of the hub genes suggested that a total number of 17 hub genes interacted with ≥5 of the hub genes were at the hub of the PPI network. Previous studies demonstrated that functional partnerships and interactions that occurred between proteins were at the core of cellular processing and their systematic characterization helped to provide context in molecular systems biology. Thus, the 17 hub genes, as mentioned above, might play an important role in the biological process of BC cells.

Our survival analyses of the hub genes indicated that lower expressions of 11 hub genes were correlated with relatively shorter overall survivals of patients with BC. Meanwhile, our GSEA results indicated that BC samples with higher expression of these genes were enriched in genes sets that were correlated with biological behaviors of tumor cells (such as P53 pathway, Apical junction, mitotic spindle, G2M checkpoint, Myogenesis, etc.).

As for the 11 survival associated hub genes, we conducted a literature review of these genes. MMP11 is a member of the matrix metalloproteinase (MMP) family involved in the breakdown of extracellular matrix in normal physiological processes, such as embryonic development, reproduction, and tissue remodeling, as well as in disease processes, and several studies demonstrated that overexpression of MMP11 were correlated with many cancers (including colon cancer, laryngeal cancer, breast cancer, prostate cancer, etc.) (Deng et al., 2005; Kou et al., 2013; Roscilli et al., 2014; Li et al., 2015; Pang et al., 2016); COL5A2 encodes an alpha chain for one of the low abundance fibrillar collagens, and previous studies suggested that COL5A2 was associated with malignancy in colorectal cancer, breast cancer, osteosarcoma, etc. (Fischer et al., 2001; Vargas et al., 2012; Wu et al., 2014); CDC25B is a member of the CDC25 family of phosphatases, and previous studies suggested that elevated expression of CDC25B promoted the proliferation of gastric cancer, esophageal carcinoma, renal cell carcinoma, etc. (Yu et al., 2012; Wang et al., 2015; Leal et al., 2016); TOP2A encodes a DNA topoisomerase, and previous studies suggested that TOP2A overexpression was a poor prognostic factor in patients with breast cancer, colorectal cancer, prostate cancer, and nasopharyngeal carcinoma, etc., (de Resende et al., 2013; Lan et al., 2014; Tarpgaard et al., 2016; Zheng et al., 2016); CENPF is a component of the nuclear matrix during the G2 phase of interphase, and previous studies demonstrated that CENPF was associated with the tumor progression of patients with prostate cancer, hepatocellular carcinomas, and breast cancer, etc. (Brendle et al., 2009; Kim et al., 2012); Cell division cycle associated 3 (CDCA3) is a part of the Skp1-cullin-F-box (SCF) ubiquitin ligase and refers to a trigger of mitotic entry and mediates destruction of the mitosis inhibitory kinase, and previous studies demonstrated overexpression of CDCA3 promoted progression of several cancers (oral cancer, prostate cancer, and lung cancer, etc.) (Uchida et al., 2012; Chen et al., 2013; O'Byrne et al., 2016); TK1(thymidine kinase 1) encodes a cytosolic enzyme that catalyzes the addition of a gamma-phosphate group to thymidine, and previous studies proved that high TK1 expression is a predictor of poor survival in patients with pT1 of lung adenocarcinoma, BC, gastrointestinal cancer, etc. (Jagarlamudi et al., 2015; Du et al., 2016); Targeting protein for Xenopus kinesin-like protein 2 (TPX2) is microtubule-associated protein and impacts spindle assembly in human cells (Huang et al., 2014), and previous studies proved that TPX2 expression was associated with poor survival in gastric cancer, breast cancer, colon cancer and esophageal cancer, etc. (Wei et al., 2013; Liu et al., 2014; Yang et al., 2015; Tomii et al., 2017); The cell division cycle associated 8 (CDCA8) gene encodes a component of the chromosomal passenger complex and plays an important role in mitosis, and overexpression of CDCA8 was reported in some human cancers, demonstrating that CDCA8 was required for the growth and progression of several cancers such as breast cancer, gastric cancer, lung carcinogenesis (Hayama et al., 2007; Yan et al., 2012; Jiao et al., 2015); Several studies demonstrated the AEBP1 upregulation conferred acquired resistance to BRAF (V600E) inhibition in melanoma and was associated with bad survivals of patients with glioma (Ladha et al., 2012; Hu et al., 2013); The protein encoded by FOXM1 is phosphorylated in M phase and regulates the expression of several cell cycle genes, such as cyclin B1 and cyclin D1, and previous studies suggested that FOXM1 promoted tumor progression of multiple cancers including gastric cancer, ovarian cancer, cervical cancer, colorectal cancer and breast cancer, etc. Trichostatin A potentiates TRAIL-induced antitumor effects via inhibition of ERK/FOXM1 pathway in gastric cancer (Barger et al., 2015; Yau et al., 2015; Zheng et al., 2015; Li et al., 2016; Wang et al., 2016; Song et al., 2017). In summary, all the conclusions were consistent with the results of survival analysis and GSEA analysis that high expressions of MMP11, COL5A2, CDC25B, TOP2A, CENPF, CDCA3, TK1, TPX2, CDCA8, AEBP1, and FOXM1 were correlated with worse overall survival of BC patients and BC samples with relatively higher expression of these genes were enriched in gene sets that were associated with cell proliferation.

In conclusion, the identified 50 hub genes that were closely correlated with the TNM staging of BC patients, and 11 hub genes (MMP11, COL5A2, CDC25B, TOP2A, CENPF, CDCA3, TK1, TPX2, CDCA8, AEBP1, and FOXM1) of which were significantly correlated with the overall survival of BC, which could be candidate biomarkers for BC. Meanwhile, further in vivo and in vitro studies were needed to make clear the exact molecular mechanisms that affected the growth of BC cells.

XW designed the study. SL and XL collected, analyzed and interpreted the data. XL, TL, XM, XY, CF, DH, and HW participated in revising the manuscript. XZ and YC participated in the study design and helped to draft the manuscript. All authors read and approved the final manuscript.