Identification of hub miRNA biomarkers for bladder cancer by weighted gene coexpression network analysis

Abstract

Bladder cancer (BC) is a common urinary system tumor with high aggressiveness, and it results in relatively high mortality due to a lack of precise and suitable biomarkers. In this study, we applied the weighted gene coexpression network analysis method to miRNA expression data from BC patients, and screened for network modules associated with BC progression. Hub miRNAs were selected, followed by functional enrichment analyses of their target genes for the most closely related module. These hub miRNAs were found to be involved in several functional pathways including pathway in cancer, regulation of actin cytoskeleton, PI3K-Akt signaling pathway, mitogen-activated protein kinase (MAPK) signaling pathway, Wnt signaling pathway, proteoglycans in cancer, focal adhesion and p53 signaling pathway via regulating target genes. Finally, their prognostic significance was tested using analyses of overall survival. A few novel prognostic miRNAs were identified based on expression profiles and related survival data. In conclusion, several miRNAs that were critical in BC initiation and progression have been identified in this study. These miRNAs, which may contribute to a comprehensive understanding of the pathogenesis of BC, could serve as potential biomarkers for BC prognosis or as new therapeutic targets.

Keywords:

• bladder cancer
• miRNA expression
• functional enrichment analysis
• bioinformatics analysis

Introduction

Bladder cancer (BC) is one of the most common malignant cancers of the urinary system and is associated with an annual mortality of over 150,000.^1,2 Bladder transitional cell carcinoma accounts for more than 90% of these deaths. Approximately 70%–75% of transitional cell carcinoma cases are non-muscle invasive bladder cancer and the remaining 25%–30% are muscle invasive bladder cancer. Approximately 70% of non-muscle invasive bladder cancer patients ultimately develop relapse within 5 years after undergoing transurethral resection of BC. Meanwhile, muscle invasive bladder cancer is likely to be linked with lethal metastasis at the time of diagnosis.^3–5 To date, sensitive and specific markers for the diagnosis and treatment of BC have not been found. Therefore, it is crucial to search for biomarkers to predict the development and prognosis of BC, which can also be used as therapeutic targets.

miRNAs are constituted by a large number of naturally occurring endogenous small non-coding RNAs. miRNAs, with lengths of 22–25 nucleotides, can regulate gene expression through translational inhibition or mRNA cleavage.⁶ A variety of studies have shown that miRNAs play important regulatory roles in diverse physiologic processes, including cell proliferation, adhesion and migration.^7,8 Furthermore, miRNAs can target multiple mRNAs, which render them powerful regulators. For example, loss of miR-218, a tumor suppressor, enhances cancer cell migration and invasion in prostate cancer through direct regulation of LIM and SH3 domain protein 1 (LASP1), which is a member of the nebulin family of actin-binding proteins.⁹ miR-124/miR-9 plays a role in the neuronal development of primary neurons by regulating glycoprotein M6A (GPM6A), a transmembrane glycoprotein.¹⁰ Moreover, a study has shown that miR-1 inhibits cell proliferation, invasion, metastasis and further progression of esophageal squamous cell carcinoma by binding to its target gene, TAGLN2, which also belongs to the family of actin-binding proteins.¹¹ In addition, abnormal miRNA expression may interfere with many cellular signaling pathways; moreover, it can profoundly affect the pathogenesis and progression of cancer. Abnormal expression of miRNAs has been observed in several kinds of human tumors, including BC.^12–14 The role of miRNAs in BC is a promising area of research to better understand the mechanisms of BC initiation, development and metastasis.

The rapid growth of high-throughput microarray technology in recent years has provided new methods to identify BC progression-related genes. Previously, identification of these genes was accomplished through basic gene expression analysis methods. Molecular mechanisms can be partially understood by screening differentially expressed genes or differentially expressed miRNAs (DEMs).^15,16 Moreover, recent bioinformatics advancements allow us to apply the weighted gene coexpression network analysis (WGCNA) algorithm to interpret expression profiles, which ensure more comprehensive identification of new biomarkers.¹⁷ WGCNA can detect highly coexpressed genes that are closely connected within the network based on similarities in sample expression profiles. These subnetwork regions are called modules and are usually composed of functionally related genes.^17,18 Therefore, the WGCNA algorithm for BC microarray miRNA expression data was performed in order to construct a coexpression network based on interactions between DEMs. In the meantime, hub miRNAs potentially associated with the pathogenesis of BC were identified and validated as prognostic biomarkers. These miRNAs may further reveal the biologic mechanisms of BC initiation and progression and may also serve as potential therapeutic targets.

Materials and methods

Acquisition of gene expression data and DEMs

The raw miRNA expression data from 418 BC and 19 normal tissue samples were downloaded from the Cancer Genome Atlas (TCGA).¹⁹ DEMs between BC and normal samples were selected using the “edgeR” R package.²⁰ miRNAs with a false discovery rate of <0.05 and | log₂ fold change | of >0.585 were identified as DEMs.

WGCNA and the identification of clinically significant modules

As reported above, the WGCNA package provided an approach to constructing miRNA modules within a network based on correlations between DEM expression profiles.^21,22 First, a Pearson’s correlation matrix was computed for all pair-wise miRNAs. Subsequently, an adjacency matrix was constructed using a power adjacency function (α_mn = Power (S_mnβ) = |S_mn|^β, α_mn = adjacency between two miRNAs, S_mn = Pearson’s correlations between two miRNAs). β is a soft thresholding parameter allowing for a maximal correlation coefficient. The power of β=3 (model fitting index R²=0.83) was chosen in accordance with the scale-free topology criterion. Next, the topological overlap matrix was calculated from the adjacency matrix to assess the interconnectedness of the miRNAs.²³ Finally, the “cutreeStaticColor” function was applied to classify similar expression miRNAs into gene modules (minModuleSize =20 and mergeCutHeight =0.95). To identify the modules that were significantly associated with BC development, gene significance (GS) was calculated according to the formula: GS_i = −log(P_i), where P_i was obtained from the value of the student’s t-test (two-tailed) to identify the differential expression between two groups. GS represented the correlation between gene expression and disease progression. The module significance was defined as the mean GS in a module, while a module with a high module significance value was considered to be the key module.²¹

Detection of candidate hub DEMs

Hub miRNAs were defined based on intramodular connectivity and module membership (MM) in our study. The intramodular connectivity was calculated as follows: $K_i={\displaystyle {\sum }_{j\ne i}^n{a}_{ij}}$, where i, j belongs to Module q, K_i = intramodular connectivity of gene i and a_ij = adjacency between genes i and j.²⁴ Intramodular connectivity reflects the correlation between miRNAs in the coexpression network. The MM was computed as originally described, namely, MM_i = |cor (x(i), ME), where the module eigengene stands for the first principal component of a given module. The MM_i ranged from 0 to 1, which reflected the Pearson’s correlation coefficient between gene expression and the module eigengene in the network.²⁵ In general, miRNAs showing a high degree of connectivity and MM values within a module were defined as “hub miRNAs”, which were assumed to be functionally significant.^22,24,26 The top 25% of miRNAs with the highest connectivity and MM values in the blue module were selected based on the module sizes as hub miRNAs for further study.

Bioinformatics analysis of miRNA target genes

Potential target genes of hub DEMs were retrieved from the target prediction database TargetScan (http://www.targetscan.org/). In addition, genes with a cumulative weighted context++ score of ≤−0.4 in the software were selected.²⁷ Functional Gene Ontology (GO) enrichment analysis and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis were performed using DAVID v6.8 (http://david.abcc.ncifcrf.gov/) bioinformatics online resources in order to uncover the biologic function of target genes. A P-value cutoff of <0.01 was set, and at least five genes were selected as the enriched function of target genes.

Human tissue specimens

Human tumor tissues and adjacent normal tissue specimens were obtained from the Department of Urology at Nanjing First Hospital, Nanjing Medical University. Ethical approval was obtained from the ethics committee of Nanjing First Hospital and written informed consent was obtained from all patients. The specimens were collected during surgery and frozen in liquid nitrogen.

Western blotting

Total protein was extracted from the tumor tissues and adjacent normal tissue specimens using radioimmunoprecipitation assay lysis buffer containing a protease inhibitor cocktail (Roche, Shanghai, China). The protein concentration was quantified using a bicinchoninic acid protein assay kit (Beyotime, Shanghai, China). Equivalent amounts of proteins were separated by sodium dodecyl sulfate-polyacrylamide gel electrophoresis, transferred to polyvinylidene fluoride membranes (EMD Millipore, Billerica, MA, USA) and blocked with 5% non-fat milk in distilled water for 2 h at room temperature. Later, the membranes were incubated with diluted anti-TAGLN2 primary antibodies (Abcam, #ab121146) overnight at 4°C, followed by incubation with horseradish peroxidase-conjugated secondary antibody (YI FEI XUE BIOTECHNOLOGY, Nanjing, China) for 2 h at room temperature. β-Actin was used as an internal control.

Human protein atlas

Expression of the TAGLN2 protein in BC tissues and normal bladder tissues was determined from immunohistochemical data publicly available in the Human Protein Atlas (www.proteinatlas.org).²⁸

Survival analyses

Survival analyses in this study were performed using the TCGA database. The miRNA expression profile and survival data from 406 BC patients were used in our study after eliminating patients without available clinical information. The prognostic values of signature candidate miRNA profiles were assessed using GraphPad Prism 7. BC samples were divided into two groups according to their median proposed miRNA expression. The survival of two patient cohorts with candidate miRNAs was compared using a Kaplan–Meier survival plot. Moreover, log-rank P-values and the hazard ratio (HR) with 95% CIs were also calculated.

Results

Identification of DEMs

Using TCGA miRNA expression data, DEM expression profiles from BC and normal tissue samples were acquired by the “edgeR” package in R, with P<0.05 and log₂ (fold change)| >0.585 being selected as the criteria. A total of 480 DEMs were identified after analysis, among which 350 genes were upregulated while 130 were downregulated (Figure 1; Table S1).

Figure 1 Volcano plot of differentially expressed miRNAs.

Abbreviations: FC, fold change; FDR, false discovery rate.

Coexpression network and clinical significant modules

The miRNA coexpression network in this study was constructed using the WGCNA. In addition, DEMs with similar expression patterns were divided into modules adopting the average hierarchical clustering linkage. Modules containing at least 20 miRNAs were extracted, and 5 modules were eventually identified (Figure 2A). The blue module had the highest mean GS among these modules (Figure 2B), and the correlation coefficient (r) of this module indicated a significant correlation of intramodular connectivity with gene significance (Figure 2C). This suggested that the blue module was closely related to BC progression. The top 25% of the 54 miRNAs (Table 1) showing the highest connectivity and MM values in the blue module were chosen as candidate hub miRNAs.²⁹

Figure 2 Identification of key modules associated with the progression of BC.

Notes: (A) Results of a cluster analysis and five modules identified from DEM expression networks. (B) The blue module had the highest mean gene significance among these modules. (C) Relationship between gene significance and intramodular connectivity.
Abbreviations: BC, bladder cancer; DEMs, differentially expressed miRNAs.

Table 1 Top 25% of the 54 miRNAs showing the highest connectivity and MM values in the blue module

microRNA	Module	Connectivity	Module membership
hsa-miR-133a-1	Blue	14.52815127	0.942911887
hsa-miR-133a-2	Blue	14.13870062	0.931813002
hsa-miR-1–2	Blue	13.22743045	0.921363767
hsa-miR-1–1	Blue	13.2169924	0.920517623
hsa-miR-133b	Blue	12.77790952	0.914496915
hsa-miR-143	Blue	11.9099187	0.890564772
hsa-miR-195	Blue	11.52715541	0.868339224
hsa-miR-139	Blue	11.20058022	0.854359765
hsa-miR-6507	Blue	10.72163066	0.82522312
hsa-miR-145	Blue	10.39113951	0.801384187
hsa-miR-28	Blue	10.24439735	0.8227846
hsa-miR-3199–2	Blue	9.044028532	0.818183518
hsa-miR-548ba	Blue	8.813461621	0.812107848

Abbreviation: MM, module membership.

Hub miRNA target prediction and functional annotation

Target genes of hub miRNAs were predicted using Target-Scan V7.1. Those with cumulative weighted context++ scores of ≤−0.4 in the software were selected (Table S2).²⁷ The miRNA regulatory network was visualized using CytoScape software, and miRNAs with many interactions were graphically highlighted (Figure 3). In total, 1351 predicted target genes were selected for further analysis. Altogether, 27 enriched GO biologic process terms (Figure 4A) and 17 KEGG pathways (Figure 4B) were obtained using a threshold of P<0.01 within the online software DAVID. GO analysis results suggested that the most represented GO terms in biologic processes were related to cell proliferation and regulation of growth. On the other hand, target genes were found to be significantly enriched in cancer-related pathways, such as pathways in cancer, regulation of actin cytoskeleton, PI3K-Akt signaling pathway, mitogen-activated protein kinase (MAPK) signaling pathway, Wnt signaling pathway, proteoglycans in cancer, focal adhesion and p53 signaling pathway.

Figure 3 A network of hub miRNAs and their target genes.

Notes: The cumulative weighted context++ scores, which indicate how effectively a given miRNA may target mRNAs, are shown by the node color. The figure was generated using CytoScape.

Figure 4 Functional annotations of the target genes.

Notes: (A) The 27 enriched GO biologic process terms. (B) The 17 enriched KEGG pathways.
Abbreviations: GO, Gene Ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes.

Validation of the identified target genes

TAGLN2, the target gene with the highest repeatability, was randomly selected for further validation in order to determine the clinical relevance of target gene expression. It was a common target gene for miR-1–1/miR-1–2, miR-133a-1/miR-133a-2 and miR-133b. Western blots showed significantly higher TAGLN2 expression levels in BC tissue than in adjacent normal tissue specimens (Figure 5A), which was consistent with immunohistochemistry from the human protein map (www.proteinatlas.org, Figure 5B). The TAGLN2 mRNA level was consistently higher in BC tissues than in normal bladder tissue (Figure 5C and D), and increased TAGLN2 expression also significantly increased the risk of advanced TNM stages (Figure 5E), as shown in Compendia Biosciences (www.oncomine.org) and TCGA BC RNAseq.

Figure 5 TAGLN2 is upregulated in human BC specimens.

Notes: (A, B) TAGLN2 expression in normal bladder tissue and BC specimens. (C) TAGLN2 expression in the TCGA BC RNAseq dataset (N=408). (D) Oncomine data showing TAGLN2 expression in normal (n=48) vs tumor (n=81) bladder specimens. (E) One-way ANOVA analysis showed that TAGLN2 expression levels were significantly associated with TNM stage in the TCGA BC RNAseq dataset.
Abbreviations: ANOVA, analysis of variance; BC, bladder cancer; N, normal; TCGA, The Cancer Genome Atlas; T, tumor.

Identification of hub miRNAs with prognostic value

Hub miRNAs in the blue module associated with the overall survival of BC patients were assessed in order to validate hub miRNAs as prognostic biomarkers. miRNA expression and clinical data of BC patients obtained from TCGA were applied toward assessing prognostic significance. As shown in Table 2, the low-expression miRNA candidate group had significantly poorer overall survival, which included hsa-miR-133a-1 (P=0.0327, HR: 1.372, 95% CI: 1.026–1.834), hsa-miR-133b (P=0.0498, HR: 1.530, 95% CI: 1.145–2.045), hsa-miR-143 (P=0.0018, HR: 1.597, 95% CI: 1.195–2.134), hsa-miR-195 (P=0.0127, HR: 1.452, 95% CI: 1.086–1.94) and hsa-miR-145 (P=0.0004, HR: 1.705, 95% CI: 1.276–2.279). These results indicate that the downregulated hsa-miR-133a-1, hsa-miR-133b, hsa-miR-143, hsa-miR-195, hsa-miR-6507, hsa-miR-145 and hsa-miR-548ba miRNAs were negative factors for the survival of BC patients.

Table 2 Overall survival of the hub miRNAs in the blue module

Hub miRNA	HR (95% CI)	P-value
hsa-miR-133a-1	1.372 (1.026–1.834)	0.0327
hsa-miR-133a-2	1.284 (0.961–1.716)	0.0921
hsa-miR-1–2	1.205 (0.902–1.610)	0.2069
hsa-miR-1–1	1.222 (0.9146–1.633)	0.1759
hsa-miR-133b	1.530 (1.145–2.045)	0.0498
hsa-miR-143	1.597 (1.195–2.134)	0.0018
hsa-miR-195	1.452 (1.086–1.940)	0.0127
hsa-miR-139	1.133 (0.848–1.514)	0.3977
hsa-miR-6507	1.585 (1.186–2.118)	0.0020
hsa-miR-145	1.705 (1.276–2.279)	0.0004
hsa-miR-28	0.823 (0.6159–1.100)	0.1878
hsa-miR-3199–2	1.107 (0.828–1.479)	0.5110
hsa-miR-548ba	1.408 (1.054–1.882)	0.0212

Note: Bold type indicates statistical significance.

Abbreviations: CI, confidence interval; HR, hazard ratio.

Discussion

BC is a fatal form of cancer, and the identification of new precise and suitable targets is of vital importance. Circulating miRNAs may serve as potential targets due to their involvement in various biologic processes, including cell proliferation, differentiation, apoptosis, metabolism and tumorigenesis. The identification of miRNA biomarkers and targets was accomplished by analyzing the microarray-based miRNA expression data obtained from TCGA database. Initially, 480 DEMs were identified in each group, including 350 upregulated and 130 downregulated miRNAs. A weighted miRNA coexpression network was subsequently constructed and several miRNA modules were identified based on DEM expression profiles. The blue module of coexpressed miRNAs was found to be of more importance than other modules, suggesting that it played a crucial role in the development of BC. Meanwhile, the top 25% miRNAs with the highest connectivity and MM values were selected as hub miRNAs for further research.

Among the blue module hub miRNAs, miR-1 and miR-133a are located in the same chromosomal region, or “cluster”. Furthermore, there are two different chromosomal loci harboring miR-1 and miR-133a clusters (miR-1–1 and miR-133a-2 on 20q13.33, as well as miR-1–2 and miR-133a-1 on 18q11.2). Yamasaki et al³⁰ reported that silencing the miR-1/miR-133a cluster could regulate the carcinogenicity of prothymosin-α and purine nucleoside phosphorylase and promote the proliferation and invasion of BC cells. In another study by Kawakami et al,³¹ the expression levels of miR-1 and miR-133a were remarkably inhibited in renal cell carcinoma cell lines and specimens, while miR-1 and miR-133a restoration could induce apoptosis in renal cell carcinoma cell lines. miR-133b, which is located in chromosomal region 6p12.2, has been shown to be downregulated in various tumors. Furthermore, it is distinctly associated with invasive clinicopathologic features, tumor subtypes and lower survival rate.^32,33 Some evidence has shown that downregulating miR-195 may inhibit apoptosis and promote the proliferation and invasion of non-small-cell lung cancer cells by upregulating the oncogene CHEK1 (checkpoint kinase 1).³⁴ Moreover, it has been found that miR-143, miR-145 and miR-28 are frequently markedly downregulated in colon cancer. These miRNAs are involved in regulation of several cellular processes, including proliferation, migration and chemoresistance.^35,36

Furthermore, the enriched GO terms and KEGG pathways of hub miRNA target genes were also characterized by functional enrichment analysis. These target genes were found to be enriched in pathways in cancer, regulation of actin cytoskeleton, PI3K-Akt signaling pathway, MAPK signaling pathway, Wnt signaling pathway, proteoglycans in cancer, focal adhesion, p53 signaling pathway, cell proliferation and regulation of growth. In particular, many of these enriched functions regulated by hub miRNAs participate in the formation and progression of multiple tumors. For instance, the MAPK cascade is a critical pathway for human cancer cell survival, dissemination and resistance to drug therapy. Numerous stimuli, including DNA damage pathways and internal metabolic stress, as well as signaling from external growth factors, cell–matrix interactions and communication from other cells, can cause increased protein amplification. Moreover, these stimuli can alter the tumor microenvironment, thus over-activating the pathway and leading to tumor progression, as other studies have shown.^37,38 Studies have demonstrated that miR-1/miR-133a regulates the MAPK pathway by inhibiting the phosphorylation of extracellular regulated protein kinases, which are involved in the metastasis of colorectal carcinoma.^39,40 The regulation of actin cytoskeleton signaling pathways is identified as a cancer-related factor in multiple cancer datasets. Many key proteins are involved in this pathway, including WASP family proteins, Arp2/3 complexes, LIM-kinases, cofilin and cortactin. These proteins are regulated by miR-143 and overexpressed in multiple cancers, and can promote the migration, invasion and metastasis of cancer cells.^41–44 The p53 gene plays an important role in various functional activities, such as avoiding the accumulation of damaged DNA, maintaining genomic stability, regulating cell differentiation by inducing cell cycle arrest, promoting apoptosis and DNA repair senescence. Cell division and proliferation will be uncontrolled if the p53 gene is damaged by cancer-induced factors.^45,46 Increasing evidence has highlighted the important roles of miRNAs in mediating tumor suppression functions of p53. Sachdeva et al have reported that miR-145 is involved in the regulation of the dynamic balance between p53, a tumor suppressor, and c-Myc, a proto-oncogene, so as to maintain normal cell growth and proliferation. The downregulation of miR-145 expression will subsequently trigger cell malignancy.⁴⁷

Notably, it is found that many genes are targeted by multiple hub miRNAs. For example, GPM6A, LASP1 and TAGLN2 are targeted by at least four hub miRNAs (Table S3). GPM6A has been shown in previous studies to promote the aggressiveness of lymphoid leukemia by inducing cell migration.⁴⁸ Zhao et al demonstrated that LASP1 might promote liver, gut and mesentery metastases by using an orthotopic xenograft mouse model of pancreatic ductal adenocarcinoma, the upregulation of which was closely related to poor prognosis for pancreatic ductal adenocarcinoma patients.⁹ In this study, TAGLN2 was randomly selected for further verification. The results suggested that TAGLN2 expression was significantly higher in tumor tissue than in adjacent normal bladder tissue, and increased TAGLN2 expression also significantly increased the risk of advanced TNM stages, indicating that these genes may be closely associated with the development of BC. Thus, further studies regarding BC initiation and progression should be conducted on these genes independently from the proposed miRNAs.

Finally, the potential prognostic value of hub miRNAs was validated based on miRNA expression profiles in BC patients, as well as their survival information in the TCGA database. Seven out of the 11 hub miRNAs (including hsa-miR-133a-1, has-miR-133b, has-miR-143, has-miR-195 and hsa-miR-145) were found to be closely associated with the overall survival for BC patients.

Conclusion

To summarize, a key module significantly associated with BC progression has been confirmed in this study. Hub miRNAs within this module may play vital roles in the initiation and progression of BC, and may therefore be candidates for functional studies. Moreover, they were found to regulate a variety of cancer-related pathways at the functional level, such as pathways in cancer, regulation of actin cytoskeleton, MAPK, Wnt, proteoglycans in cancer, focal adhesion and p53-associated cellular signaling pathways, as well as cell proliferation and various cancer-related genes. In addition, some miRNAs were also found to be associated with overall survival for BC patients. Our study not only provides new insights into the initiation and progression of BC, but also identifies new potential biomarkers and therapeutic targets. Nonetheless, further studies are still needed.

Acknowledgments

This study was supported by grants from the National Natural Science Foundation of China (81370853 and 81570613) to RPJ, and Research and Innovation Program for Graduates of Jiangsu Province (CXZZ13_0583), Jiangsu Provincial Medical Youth Talent (QNRC2016073), Key Program of Science and Technology Foundation from Nanjing Medical University (2016NJMUZD042), and Natural Science Foundation of Jiangsu Province (BK20170135) to YZG.

Disclosure

The authors report no conflicts of interests in this work.