Dysregulated genes targeted by microRNAs and metabolic pathways in bladder cancer revealed by bioinformatics methods.

The present study aimed to identify bladder cancer-associated microRNAs (miRNAs) and target genes, and further analyze the potential molecular mechanisms involved in bladder cancer. The mRNA and miRNA expression profiling dataset GSE40355 was downloaded from the Gene Expression Omnibus database. The Limma package in R was used to identify differential expression levels. The Human microRNA Disease Database was used to identify bladder cancer-associated miRNAs and Target prediction programs were used to screen for miRNA target genes. Enrichment analysis was performed to identify biological functions. The Database for Annotation, Visualization and Integration Discovery was used to perform OMIM_DISEASE analysis, and then protein-protein interaction (PPI) analysis was performed to identify hubs with biological essentiality. ClusterONE plugins in cytoscape were used to screen modules and the InterPro database was used to perform protein domain enrichment analysis. A group of 573 disease dysregulated genes were identified in the present study. Enrichment analysis indicated that the muscle organ development and vascular smooth muscle contraction pathways were significantly enriched in terms of disease dysregulated genes. miRNAs targets (frizzled class receptor 8, EYA transcriptional coactivator and phosphatase 4, sacsin molecular chaperone, calcium voltage-gated channel auxiliary subunit β2, peptidase inhibitor 15 and catenin α2) were mostly associated with bladder cancer. PPI analysis revealed that calmodulin 1 (CALM1), Jun proto-oncogene, AP-1 transcription factor subunit (JUN) and insulin like growth factor 1 (IGF1) were the important hub nodes. Additionally, protein domain enrichment analysis indicated that the serine/threonine protein kinase active site was enriched in module 1 extracted from the PPI network. Overall, the results suggested that the IGF signaling pathway and RAS/MEK/extracellular signal-regulated kinase transduction signaling may exert vital molecular mechanisms in bladder cancer, and that CALM1, JUN and IGF1 may be used as novel potential therapeutic targets.

Introduction

Bladder cancer is reported to be one of the four most commonly occurring genitourinary tumors in China (1), and is one of the most expensive types of cancer to treat due to the essential life-long surveillance involving upper tract imaging, urinary cytology and cystoscopy (2). The disease usually presents in two distinct forms: Non-muscle-invasive tumors in clinical stages Ta and T1 and muscle-invasive cancer occurring in clinical stages T2-4 (3). As bladder cancer may result from numerous processes with accumulation of genetic and epigenetic changes, high-throughput technologies have been used to generate numerous genetic and genomic datasets to uncover disease causal genes and their actions involved in the initiation and development of bladder cancer (4,5). A number of studies have documented a strong link between microRNA (miRNA/miR) function and cancer pathogenesis, including bladder cancer (6–8).

miRNAs represent a class of naturally non-coding RNAs (21–23 nucleotides in length), which act as post-transcriptional silencing modulators of target genes either by negatively regulating specific target mRNAs or by inhibiting target protein synthesis (9,10). A previous study reported a potential diagnostic role of miRNAs, as the miRNAs-profiling revealed a more accurate diagnostic effect compared with the mRNAs-classifiers (7). miRNAs have received increasing attention in cancer genomic research and the novel oncogenic or tumor suppression functions of miRNAs have been revealed in a clinical study (6). Certain studies have provided a detailed insight into regulatory networks involving bladder cancer by integrated analysis of miRNA and mRNA data (5,7). Altered miRNA expression levels have been reported in bladder cancer (3,11); however, the target genes of these miRNAs have not yet been fully elucidated, particularly in bladder cancer, and further analysis is required to elucidate the subsequent processes.

Biological interpretation of large gene lists derived from microarray analysis may require researchers to select the most useful genes for further investigation. The Database for Annotation, Visualization and Integration Discovery (DAVID) has emerged as a publicly available high-throughput functional annotation tool (12). Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis (13) and Gene Ontology (GO) term analysis (14) are able to identify the main functional and metabolic signaling pathways of differentially expressed genes (DEGs). The present study used KEGG pathways (http://www.genome.jp/kegg/pathway.html) and GO term (http://geneontology.org/) analysis for enrichment analysis.

In addition, genes with similar functions have been reported to interact with each other closely, as presented in the protein-protein interaction (PPI) network, which provides a global picture used to understand molecular mechanisms and biological processes of diseases, in particular to analyze cancer (15,16). The PPI network has been used to systematically analyze and compare the disease genes that would otherwise not be identified by single gene analysis (17), in order to gain insights into distinct topological features of cancer genes using cancer-associated sub-networks. Hub nodes are known to be the most important nodes in the PPI network, as the corresponding proteins are important proteins in metabolic networks. In the present study, disease dysregulated genes were mapped into the PPI network, where proteins are represented as nodes and interactions as edges, in order to understand key biological processes involved in bladder cancer in a global sense.

miRNAs exert their function through translational repression or degradation of mRNA targets, and of particular note is the fact that no study has yet systematically analyzed the significant mechanisms underlying miRNA targets contributing to bladder cancer incidence. In the present study, numerous computational methods were used to analyze microarray data of miRNAs and mRNAs between bladder cancer and normal samples, in an attempt to describe the variation of miRNA targets and the potential mechanisms involved in the development of bladder cancer.

Materials and methods

miRNA and mRNA expression profiling data

Gene Expression Omnibus (GEO; http://www.ncbi.nlm.nih.gov/geo/), which is established at the National Center for Biotechnology Information, provides high-throughput gene expression data from the research community worldwide. The microarray expression profiling data of miRNAs and mRNAs in bladder urothelial carcinoma (‘case’) and normal bladder tissue samples (‘normal’) were downloaded from the GEO database (dataset GSE40355). Detailed messages were accessible as follows: Microarray data of miRNAs included 8 normal (GSM881439-GSM881446) and 8 case chips (GSM881455-GSM881462); and microarray data of mRNAs included 8 normal (GSM991923-GSM991930) and 8 case chips (GSM991939-GSM991946). The expression levels of miRNAs and mRNAs were previously studied in these samples on two platforms: GPL8227 Agilent-019118 Human miRNA Microarray 2.0 G4470B (miRNA ID version) and GPL13497 Agilent-026652 Whole Human Genome Microarray 4×44K v2 (Probe Name version) (both Agilent Technologies, Inc., Santa Clara, CA, USA) (18).

Identification of differentially expressed miRNAs and genes

The raw microarray expression data of miRNAs and mRNAs downloaded as Series Matrix files from the GEO database were mapped to the corresponding genes according to the SOFT formatted family files from GEO database. For each sample, the expression values of all probes for the same miRNA or mRNA were reduced to a single value by calculating the average expression value. Bayesian method was used for multiple comparisons under the control of false discovery rate (FDR) (19). Limma package in R software (version 3.24.15, https://www.bioconductor.org/packages/3.1/bioc/html/limma.html; Bioconductor; Fred Hutchinson Cancer Reearch Center, Seattle, WA, USA) (20), a linear regression model, was used to identify differentially expressed miRNAs and genes with the cut-off criteria of |log fold change (FC)|>1, P<0.05 and FDR<0.01 (21).

Identification of disease dysregulated genes

To identify the disease miRNAs associated with bladder cancer, the selected differentially expressed miRNAs were mapped to the online software, Human microRNA Disease Database (HMDD; version 2.0; 2014.06.14 update; http://www.cuilab.cn/hmdd). As the functions of miRNAs were dependent on its target genes (9), the target genes of disease miRNAs were screened using target prediction programs. A single miRNA can regulate numerous mRNAs and a single mRNA may be targeted by different miRNAs. In order to ensure a high level of accuracy, five target prediction programs, including miRanda http://www.microrna.org), MirTarget2 (http://mirdb.org/miRDB), PicTar (http://www.pictar.org/), PITA (http://genie.weizmann.ac.il/pubs/mir07/mir07_data.html) and TargetScan (http://www.targetscan.org/) were used to screen potential targets of miRNAs (22). To decrease the number of false-positive results, putative target genes predicted by at least 3 programs of the five were accepted as the disease miRNAs targets.

Finally, to identify which differentially expressed miRNAs targets were dysregulated in bladder cancer, the present study selected the intersection of disease miRNA targets and DEGs via the mRNAs chip data obtained above. The common genes were considered as the disease dysregulated genes for further analysis.

Enrichment analysis and cancer-associated gene screening

To further understand gene functions and identifying pathways closely associated with disease dysregulated genes, the functional annotation of associated genes was performed using the online software, DAVID. The present study used KEGG pathway and GO term analysis for the following enrichment analysis of disease dysregulated genes to identify the function and pathway alterations involved in bladder cancer, with P<0.05 used as the cut-off criterion. Furthermore, DAVID was used to perform OMIM_DISEASE analysis, in order to distinguish cancer-associated genes from the disease dysregulated genes.

PPI network and miRNA-gene regulatory network

Highly complex cellular functions are the results of tightly regulated interactions of groups of proteins encoded by genes. Interactions between proteins possess a crucial role in modern biology, thus a rapidly growing number of technologies have been developed for the global charting of the PPI network. The Search Tool for the Retrieval of Interacting Genes (STRING), an online database resource (http://string-db.org/), provides uniquely comprehensive information for assembling, evaluating and disseminating PPI in a user-friendly manner (23).

In the present study, STRING was used to screen the PPI of disease dysregulated genes, which was then visualized using Cytoscape software (http://cytoscape.org/; version 3.2.1; National Institute of General Medical Sciences, Seattle, WA, USA). The total connectivity degree of each node in the network was calculated using Cytoscape software and only the nodes with ≥30 degrees were selected as the important nodes (Hubs) of PPI. Additionally, the integrated miRNA-disease dysregulated genes regulatory network was also constructed using Cytoscape software, followed by network module screening using the threshold of P<0.05 by employing the clusterONE plugins in Cytoscape (24). Subsequently, protein domain enrichment analysis was performed on the screened modules using the InterPro database (25).

Results

Differentially expressed clustered miRNAs and genes

Based on the Limma package in R language, under cut-off criteria of |logFC|>1, P<0.05 and FDR<0.01, 112 miRNAs and 2,290 genes were identified to be differentially expressed in bladder cancer samples compared with the control groups. There were 56 upregulated and 56 downregulated miRNAs; as well as 1,444 upregulated and 846 downregulated DEGs.

Target screening of the differentially expressed miRNAs

Firstly, differentially expressed miRNAs associated with bladder cancer were screened from the 112 miRNAs using the HMDD database and 71 bladder cancer-associated miRNAs were obtained from this procedure, including hsa-miR-143, hsa-miR-205, hsa-miR-21, hsa-miR-93, hsa-miR-145 and hsa-miR-143. Subsequently, targets of the 71 disease miRNAs were identified using five target prediction programs, including miRanda, MirTarget2, PicTar, PITA and Target. For each target selection, the present study selected the target genes of miRNAs coexisting in at least three of the databases.

A total of 573 disease dysregulated genes associated with bladder cancer were identified by taking the intersection of disease miRNA targets and 2,290 differentially expressed genes corresponding to mRNA datasets.

Enrichment analysis of disease dysregulated genes

To analyze the biological functions of disease dysregulated genes in bladder cancer, GO and KEGG pathway enrichment analyses were performed and the results are presented in Tables I and II, respectively. Following combination of the results from the GO and KEGG pathway analyses, it was revealed that the most marked function and pathway involved in bladder cancer were muscle organ development and vascular smooth muscle contraction, respectively. Furthermore, DAVID was used to perform OMIM_DISEASE analysis, and a set of six target genes [frizzled class receptor 8 (FZD8), sacsin molecular chaperone (SACS), calcium voltage-gated channel auxiliary subunit β2 (CACNB2), peptidase inhibitor 15 (PI15) and catenin α2 (CTNNA2)] were identified to closely associate with the functional term of disease of cancer, suggesting a potential association with bladder cancer.

Table I.

Top five enriched GO terms.

Category	Term	Count	P-value
GOTERM_BP_FAT	GO:0007517~muscle organ development	34	3.83×10−14
GOTERM_BP_FAT	GO:0042692~muscle cell differentiation	22	2.54×10−10
GOTERM_BP_FAT	GO:0014706~striated muscle tissue development	21	1.23×10−9
GOTERM_BP_FAT	GO:0060537~muscle tissue development	21	3.00×10−9
GOTERM_BP_FAT	GO:0007167~enzyme-linked receptor protein signaling pathway	35	5.26×10−9

Count represents the number of genes which were enriched in the corresponding functional category. GO, Gene Ontology; GOTERM_BP_FAT, Gene Ontology term-biological process.

Table II.

Top five enriched KEGG pathways.

Category	Term	Count	P-value
KEGG_PATHWAY	hsa04270: Vascular smooth muscle contraction	15	2.45×10−5
KEGG_PATHWAY	hsa04360: Axon guidance	15	1.19×10−4
KEGG_PATHWAY	hsa04510: Focal adhesion	19	1.61×10−4
KEGG_PATHWAY	hsa04310: Wnt signaling pathway	15	6.22×10−4
KEGG_PATHWAY	hsa04020: Calcium signaling pathway	15	2.74×10−3

Count represents the number of genes which were enriched in the corresponding functional category. KEGG, Kyoto Encyclopedia of Genes and Genomes.

PPI network and integrated miRNA-target network construction

In order to gain further insights into the changes of biological pathways in bladder cancer, the present study used Cytoscape to identify the PPI of disease dysregulated genes (Fig. 1). The upregulated calmodulin 1 (CALM1), jun proto-oncogene (JUN) and insulin-like factor 1 (IGF1) were key hub nodes in the PPI network. Cytoscape was also used to perform regulatory network construction of disease dysregulated genes and the corresponding miRNAs, as shown in Fig. 2.

Figure 1.

Protein-protein interaction network of disease dysregulated genes. The red nodes represent upregulated genes and the green nodes represent downregulated genes. Rhombuses represent hub nodes.

Figure 2.

Integrated target miRNA interaction network. Red represents upregulation and green represents downregulation. Circles represent disease dysregulated genes and rhombuses represent corresponding miRNAs. miRNAs, microRNAs.

Module screening and subsequent protein domain enrichment analysis

ClusterONE plugins were used to perform module clustering analysis in the PPI network, with a criterion of P<0.05. A total of three clustering module were identified: Module 1 (density, 0.561; quality, 0.732; P=1.582×10−7), module 2 (density, 0.714; quality, 0.625; P=0.002) and module 3 (density, 0.526; quality, 0.512; P=0.004) (Fig. 3). The protein domain enrichment analysis revealed that there were no significant enrichment results in module 2 and module 3, whereas in module 1 the contained disease dysregulated genes, including SHC SH2-domain binding protein 1 (SHCBP1), were mainly enriched at IPR008271: Serine/threonine protein kinase active site (Table III).

Figure 3.

Corresponding clustering network modules extracted from the protein-protein interaction network.

Table III.

Protein Domain enrichment analysis on module 1 using the InterPro database.

Category	Term	Count	P-value
INTERPRO	IPR008271: Serine/threonine protein kinase, active site	5	7.42×10−4
INTERPRO	IPR017442: Serine/threonine protein kinase-related	5	7.82×10−4
INTERPRO	IPR017441: Protein kinase, ATP binding site	5	1.88×10−3
INTERPRO	IPR000961: AGC-kinase, C-terminal	3	1.96×10−3
INTERPRO	IPR000719: Protein kinase, core	5	2.22×10−3
INTERPRO	IPR002290: Serine/threonine protein kinase	4	3.48×10−3
INTERPRO	IPR015633: E2F Family	2	9.57×10−3
INTERPRO	IPR003316: Transcription factor E2F/dimerization partner	2	1.31×10−2
INTERPRO	IPR014400: Cyclin, A/B/D/E	2	1.43×10−2
INTERPRO	IPR004367: Cyclin, C-terminal	2	1.67×10−2
INTERPRO	IPR006671: Cyclin, N-terminal	2	3.89×10−2
INTERPRO	IPR013763: Cyclin-related	2	4.35×10−2
bINTERPRO	IPR006670: Cyclin	2	4.70×10−2
INTERPRO	IPR017892: Protein kinase, C-terminal	2	4.70×10−2
INTERPRO	IPR001752: Kinesin, motor region	2	4.93×10−2
INTERPRO	IPR019821: Kinesin, motor region, conserved site	2	4.93×10−2

Count represents the number of genes which were enriched in the corresponding functional category. ATP, adenosine triphosphate.

Discussion

Recently, emerging bioinformatic methods have accelerated the progress made to identify the mechanisms involved in bladder cancer development at the molecular level. miRNAs were reported to regulate the activity of ~30% of all protein-coding genes in the human genome (26), whereas the deregulation of miRNAs and their targeted mRNAs may serve an important function in the development of cancer. It was previously revealed that aberrantly expressed miRNAs served significant roles in the initiation, development and metastasis of bladder cancer (3,27). Previous studies have identified the large set of differentially expressed miRNAs between bladder cancer and normal samples (11,27,28). The present study further analyzed the expression profile data of miRNAs and mRNAs, and identified a total of 573 differentially expressed targets of the differentially expressed miRNAs in bladder cancer. These genes provided a global perspective for understanding the mechanism underlying bladder tumorigenesis. In addition, the results of functional enrichment analysis indicated that the genes were most significantly associated with the GO term of muscle organ development and vascular smooth muscle contraction pathway. Genes, including FZD8, EYA transcriptional coactivator and phosphatase 4 (EYA4), SACS, CACNB2, PI15 and CTNNA2 in the PPI network, were demonstrated to be potentially associated with bladder cancer by OMIM_DISEASE analysis. Furthermore, CALM1, JUN and IGF1 were important hub nodes in the PPI network with high degrees of connectivity.

Previous evidence suggested that there were two distinct and overlapping stages of bladder cancer: Non-muscle-invasive and progressive muscle-invasive tumors (3). GO term and KEGG pathway enrichment analysis of disease dysregulated genes revealed that muscle organ development and vascular smooth muscle contraction, respectively, were involved in bladder cancer, which suggested a vital role of normal muscle function in preventing bladder cancer progression. On the molecular level, a previous study reported the participation of tumor-initiating cells, i.e., cancer stem cells, in urothelial cancer (29). FZD8, EYA4, SACS, CACNB2, PI15 and CTNNA2 were demonstrated to be the cancer-associated target genes of differential miRNAs. Tumor suppression miRNAs targeting FZDs were reported to participate in the suppression of cell growth and metastasis via the Wnt signaling pathway, which served important roles in tumorigenesis (30). EYA4, PI15 and CTNNA2 functioned as tumor suppressors in the tumorigenesis of Barrett's esophageal cancer (31), glioblastoma (32) and laryngeal carcinomas (33), respectively. CACNB2 was downregulated in breast cancer cells, thereby causing poor signaling by negatively affecting the calcium metabolism (34). SACS targeted by miR-223 was reported to be upregulated in urothelial carcinoma in comparison to that in the normal mucosa (11). Notably, in the present study, the expression levels of 6 miRNA targets were revealed to be markedly increased in bladder cancer compared with those in the normal samples, which implied the potential involvement in bladder cancer development.

In order to further understand the oncogenic mechanisms underlying the screened miRNA targets involved in bladder cancer, the present study performed PPI analysis and identified three hub nodes: CALM1, JUN and IGF1. PPI networks are typically scale-free, in which hubs are ‘highly connected’ and tend to correspond to essential biological genes (35). IGF signaling, possessing an important effect on cell proliferation and differentiation, was implicated in numerous types of cancer and was thus considered to be a therapeutic target in cancer treatment (36). Dysregulated miRNAs targets were demonstrated to be involved in IGF1R signaling, as well as in its downstream antipoptotic signaling pathways, including phophoinisitide-3 kinase (PI3K)/protein kinase B (Akt)/mechanistic target of rapamycin (mTOR), in bladder urothelial carcinoma (37). Notably, the activation of the PI3K/Akt/mTOR signaling pathway was reported to be a candidate driver of the muscle-invasive phenotype of bladder cancer (38). In the present study, the upregulation of IGF1 suggested possible active IGF1R signaling and subsequent cascade dysregulation in bladder cancer. CALM, a ubiquitous Ca2+ receptor protein interacting with hundreds of different target proteins, revealed upregulated expression levels in numerous tumor cells in comparison with that in cells from normal tissues. CALM antagonists may inhibit tumor cells growth of different origins by affecting the CALM-dependent mechanisms (39). JUN, which is important for cell proliferation, survival and apoptosis, was reported to be a crucial contributing factor for tumorigenesis due to its overexpression in numerous types of human cancer (40). The results of the present study revealed that the increased expression levels of IGF1, CALM1 and JUN in bladder cancer may be potentially effective therapeutic targets for bladder cancer treatment.

In addition, the serine/threonine protein kinase active site was revealed to be enriched in genes contained in module 1, which was extracted from the PPI network. The human protein kinase family, consisting of 518 genes, was classified as protein-serine/threonine kinases, protein-tyrosine kinases and tyrosine-kinase-like proteins based on the nature of the phosphorylated-hydroxy group (41,42). The protein-serine/threonine kinases that participated in the RAS/RAF/MEK/extracellular signal-regulated kinase (ERK) signaling pathway were reported to be associated with retroviral oncogenes and were thereby attractive cancer drug targets (43–45). SHCBP1, a dysregulated gene in module 1, was previously demonstrated to regulate the expression of activated ERK1/2 and thus be necessary for the RAS/MEK/ERK signaling pathway (46). In the present study, the dysregulation of the serine/threonine protein kinase presented in module 1 may suggest abnormal RAS/MEK/ERK transduction signaling involved in bladder cancer progression.

In summary, the findings of the present study suggested the potential therapeutic ability of the IGF signaling pathway and RAS/MEK/ERK transduction signaling in bladder cancer, as well as bladder cancer-associated oncogenes, including FZD8, EYA4, SACS, CACNB2, PI15 and CTNNA2. Notably, CALM1, JUN and IGF1, which were identified as hubs in the PPI network, may be important therapeutic targets for bladder cancer intervention. However, further research, both epidemiological and mechanistic, is required to further validate the results and clarify the underlying molecular mechanisms.