This article is mentioned in:

Introduction

Bladder cancer (BC) is estimated to be among the top five most common types of cancer in western countries and ranks number 13 in terms of cancer-associated mortality worldwide (1,2). Histologically, BC may be classified based on the depth of invasion: pTa, papillary; pT1, lamina propria invasion; pT2, muscle invasive; pT3, invasion to peri-vesical fat; and pT4, locally advanced (3,4). BC results from long-term exposure to contaminants or other environmental factors involving gene mutations and progressive cellular damage. Parkin et al (5) demonstrated that the incidence of BC is almost four times higher in men than in women and inducing factors include tobacco smoke, prolonged exposure to chemical substances and race (6,7).

Genetic mutations in gene expression may lead to the malignant transformation of bladder cells. High-throughput DNA microarray analyses have identified multiple DNA mutations and alterations in the genesis of BC; these include genes encoding for B cell lymphoma-2, p53, H-Ras and fibroblast growth factor receptor 3 (FGFR3) (8,9). FGFR3 is a tyrosine kinase receptor that regulates fundamental developmental pathways and triggers a range of cellular processes, including proliferation, differentiation, migration and apoptosis (10). The FGFR3 gene is located on chromosome region 4p16.3 (11) and it is made up of 9 exons and 18 introns (12). The structure of FGFR3 is composed of an extracellular domain consisting of two or three immunoglobulin-like domains, a hydrophobic transmembrane domain and an intracellular tyrosine kinase domain (13). Upon ligand binding, FGFR3 forms dimers and activates the intracellular kinase domain, resulting in autophosphorylation of this domain. The phosphorylated residues are the binding targets of the adaptor proteins and their binding results in the activation of several signal transduction pathways, including the Ras mitogen-activated protein kinase (MAPK) signaling pathway and phosphoinositide 3-kinase (PI3K) Akt mammalian target of rapamycin (mTOR) pathway (14).

There are two mechanisms to explain the abnormal activation of FGFR3: Overexpression or activating mutations. FGFR3 mutations have been identified in multiple dwarfisms (15), such as hypochondroplasia, and in multiple types of cancer, including prostate cancer (16), cervical cancer (17) and BC. FGFR3 mutations were reported in BC for the first time by Cappellen et al (18). There is evidence to suggest that codons 248, 249 and 375 are the major mutation hot spots in BC (19); in addition, activating mutations of FGFR3 have been revealed primarily in pTa (60–70%) and in pT1-4 (16–20%) (19). Overexpression of FGFR3 has been frequently identified in BC; furthermore, Jebar et al (20) demonstrated that the expression of FGFR3 was higher in low stage BC. Du et al (21) identified a gene pathway linking FGFR3 with sterol and lipid metabolism through transcriptional profiling of BC cells subjected to short hairpin (sh)RNA knockdown of FGFR3 (21). FGFR3 has been demonstrated to be a promising therapeutic target for BC (22,23). However, the molecular mechanisms of FGFR3 activation, via overexpression or activating mutation, in BC remain to be elucidated.

The present study aimed to analyze microarray data in order to investigate the changes in gene expression profiles that occur following loss of FGFR3; in addition, the current study aimed to explore the target genes and molecular mechanisms of FGFR3 The genes that were differentially expressed in FGFR3-deleted cell lines as compared with the control cell lines were considered to be potential transcriptional targets of over-expressed FGFR3 in bladder cancer. Furthermore, a protein-protein interaction (PPI) network was constructed and the disturbed biological pathways were identified following FGFR3 knockdown in order to explore the pathogenesis and occurrence of BC associated with FGFR3.

Materials and methods

Messenger RNA expression profile data of BC

The transcription profile dataset of BC was obtained from National Center of Biotechnology Information Gene Expression Omnibus database (http://www.ncbi.nlm.nih.gov/geo/). The accession number was GSE41035 and the dataset consisted of a total of 24 mRNA samples, including 18 experimental samples collected from RT112 cell lines, with FGFR3 shRNA 2–4, FGFR3 shRNA 4-1 or FGFR3 shRNA 6–16, as well as 6 control enhanced green fluorescent protein (EGFP) shRNA samples. The platform used was GPL570 Affymetrix Human Genome U133 Plus 2.0 array (Affymetrix, Inc., Santa Clara, CA, USA). The original CEL files and the annotations file were downloaded based on this platform.

Identification of differentially expressed genes (DEGs)

Probe-level data in the CEL files were first converted into expression measures. For each sample, the expression values of all probes for a given gene were reduced to a single value by taking the average expression value. Subsequently, missing data was imputed and quartile data normalization was performed by robust multichip averaging using Affy package in R software (version 3.1; http://www.bioconductor.org/packages/release/bioc/html/affy.html) (24). The Limma package version 3.24.2 (http://www.bioconductor.org/packages/release/bioc/html/limma.html) (25) in R language with multiple testing correction was then used according to the Benjamini & Hochberg method (26) in order to identify DEGs between BC samples and normal controls. P<0.05 and |log(fold change; FC)|>1 were defined as the thresholds.

Gene ontology (GO) enrichment analysis

In order to investigate DEGs at the molecular and functional level, the online biological tool, Database for Annotation, Visualization and Integrated Discovery (DAVID) version 6.7 (http://david.abcc.Ncifcrf.gov/), was used for GO term enrichment and genes were clustered according to GO. GO is a collection of controlled vocabularies, which include molecular function, cellular component and biological process, to describe the biology of a gene product in any organism. P<0.05 was selected as the cut-off criterion during the analysis.

Pathway enrichment analysis

The theoretical principle for enrichment analysis is that associated functional genes are more likely to be selected in the abnormal biological process by the high-through screening technologies (27). Based on the selected genes, researchers are able to correctly identify the biological processes involved. In order to identify the enriched pathways of DEGs, DAVID was used with P<0.05 as the threshold. The pathways used as DAVID input for cluster analysis were from Kyoto Encyclopedia of Genes and Genomes (KEGG; http://www.genome.jp/kegg/) and BIOCARTA (http://www.biocarta.com/).

PPI network construction

PPIs are crucial for all biological processes. In the present study, the PPI network was constructed based on the Protein Interaction Network Analysis platform (PINA2) database (http://cbg.garvan.unsw.edu.au/pina/). PINA2 (28) is a database containing known and predicted associations of protein interaction. The interactions include direct (physical) and indirect (functional) associations. Of note, the protein names in the Universal Protein Resource database (http://www.ebi.ac.uk/uniprot/remotingAPI/), which correspond to the DEGs, were submitted to construct the PPI network. Here Cytoscape software (version 3.2.1; http://cytoscape.org/) (29) was used to visualize the PPI network to further observe the associations between genes.

Analysis of co-expressed genes

DEGs only explain a limited number of mechanisms of FGFR3-shRNA in BC. However, analysis of differential co-expression genes may reveal two or several similar genes with similar expression patterns across a set of samples. Co-expression genes were hypothesized to have a functional association, such as physical interaction between the encoded proteins (30,31). In order to further explore the pathogenesis of BC in the present study, the differential coexpression enrichment (DCe) function in DGCL package (32) version 2.1.2 (http://cran.r-project.org/web/packages/DCGL/index.html) in R language was used and the parameters in the function were set to default values. P<0.05 and the maximum absolute correlation coefficient >1.5 were set as thresholds.

Results

Identification of DEGs

The Limma package was used to analyze the transcription profile data between the experimental and control samples. P<0.05 and |logFC|>1 were used as the significant thresholds for DEGs. Based on these criteria, a total of 196 DEGs were identified, among which 101 were downregulated and 95 were upregulated.

GO analysis

Functional classification was performed using the online biological tool DAVID, with a threshold of P<0.05. Table I demonstrates the top ten significantly enriched GO terms when these DEGs were classified according to biological process. This analysis revealed that the most enriched functions detected in FGFR3 knockdown samples compared with control samples were the biosynthesis of sterol (P=8.87×10−8)and metabolism of sterol (P=1.58×10−7) or steroid (P=3.52×10−7). In addition, oxidation reduction, extracellular region part and identical protein process were also demonstrated to be enriched.

Table I. Clustering of differentially expressed genes based on biological process.

Table I.

Clustering of differentially expressed genes based on biological process.

GO ID	GO name	Gene numbera	P-value
GO:0016126	Sterol biosynthetic process	8	8.87×10−8
GO:0016125	Sterol metabolic process	11	1.58×10−7
GO:0008202	Steroid metabolic process	14	3.52×10−7
GO:0008203	Cholesterol metabolic process	10	7.47×10−7
GO:0055114	Oxidation reduction	22	6.80×10−6
GO:0006695	Cholesterol biosynthetic process	6	8.27×10−6
GO:0044421	Extracellular region part	27	3.78×10−5
GO:0006694	Steroid biosynthetic process	8	4.11×10−5
GO:0042802	Identical protein binding	20	4.16×10−5
GO:0008299	Isoprenoid biosynthetic process	5	5.97×10−5

a Number of differentially expressed genes enriched in this GO term. GO, gene ontology.

Analysis of the biological pathways of DEGs

The gene transcription profile was significantly altered in experimental samples compared with control samples. These DEGs were selected for KEGG and BIOCATRA pathway enrichment analysis. As shown in Table II, according to the threshold of P<0.05, six biological pathways were significantly enriched. Consistent with the results of the GO enrichment analysis, these pathways were primarily associated with biosynthesis and metabolism, including steroid biosynthesis (P=4.80×10−5) and arachidonic acid metabolism (P=0.03285).

Table II. Biological pathways in bladder cancer cells.

Table II.

Biological pathways in bladder cancer cells.

Category	Pathway name	P-value	Gene name
KEGG	Steroid biosynthesis	4.80×10−5	CYP51A1, SQLE, DHCR7, FDFT1, SC4MOL
KEGG	Terpenoid backbone biosynthesis	7.77×10−4	HMGCR, FDPS, IDI1, ACAT2
BIOCARTA	Fibrinolysis pathway	0.009393	SERPINB2, PLAU, F2R
KEGG	Arachidonic acid metabolism	0.03285	PLA2G4A, PLA2G10, GPX3, CYP4F3
KEGG	Complement and coagulation cascades	0.045395	F5, CFD, PLAU, F2R
BIOCARTA	Nuclear receptors in lipid metabolism and toxicity	0.047437	CYP24A1, ABCG1, CYP4B1

[i] KEGG, Kyoto Encyclopedia of Genes and Genomes.

PPI network construction

DEGs were mapped to a PINA2 database and a PPI network was constructed. The PPI network identified 1,865 genes, with 2,482 interactions between them (Fig. 1). Genes that had more interactions with other genes may have a critical role in BC. Table III demonstrates the DEGs with the ten highest degrees of gene interactions identified in these samples. The gene SETX was reported to have the highest degree (degree, 159), which indicated that it may have an important role in FGFR3-regulated BC.

Figure 1. Protein-protein interaction network for products of DEGs. A total of 1,865 gene nodes and 2,482 interaction associations were identified. Triangular nodes, 146 DEGs; circular nodes, 1,719 non-DEGs. DEGs, differentially expressed genes. For each node, the degree represents the number of interactions with other nodes and node size is proportional to the degree of each node.

Table III. Ten highest degrees of interaction and the corresponding differentially-expressed genes.

Table III.

Ten highest degrees of interaction and the corresponding differentially-expressed genes.

Gene name	Degree
SETX	159
CFD	127
HIST2H2BE	110
NME2	88
SERPINB2	85
UBC	80
CSK	77
DHRS2	75
TRIP6	72
MCM5	60

Analysis of co-expressed genes

Co-expressed genes were identified using DCe function in R language; a gene co-expression network was then constructed and visualized using Cytoscape. The network involved 1,935 co-expression associations between 140 differentially co-expressed genes and 798 non-differentially co-expressed genes (Fig. 2). Table IV lists the differentially co-expressed genes with the ten highest degrees of gene interactions, among which CCNB1 was the top hub gene with the highest degree (degree, 39).

Figure 2. Gene co-expression network. A total of 938 gene nodes and 1,935 associations were identified. Triangular nodes, 140 differentially co-expressed genes; circular nodes, 798 non-differentially co-expressed genes. For each node, the degree represents the number of interactions with other nodes and node size is proportional to the degree of each node.

Table IV. Ten highest degrees of interaction and the corresponding differentially co-expressed genes.

Table IV.

Ten highest degrees of interaction and the corresponding differentially co-expressed genes.

Gene name	Degree
CCNB1	39
MCM3	26
TMEM97	25
MCM5	25
H2AFZ	25
PPIL1	25
UBE2T	22
MCM2	21
CTSF	19
ZIC2	19

Discussion

BC is the fourth most common type of solid cancer in men and the seventh most common in women, worldwide. Although mutations and overexpression of FGFR3 have been associated with BC, the biological mechanisms underlying its pathogenesis remain to be fully elucidated. The present study aimed to analyze DEGs based on the transcription profile data of experimental FGFR3 knockdown samples and control EGFP shRNA samples. A total of 196 genes were identified to be differentially expressed in FGFR3 knockdown bladder cancer cell lines in comparison with control cell lines. GO and pathway enrichment analysis were conducted; in addition, a PPI network of DEGs and a differential co-expression network were constructed.

In the present study, the results of the GO analysis as well as the KEGG and BIOCATRA pathway enrichment analysis revealed that the primary biological process in which these DEGs were involved was sterol biosynthesis and metabolism. These results were consistent with those of a previous study, which demonstrated that FGFR3 may affect BC through the biosynthesis of sterol and lipids (21). Activation of FGFR3 was reported to promote the accumulation of mature sterol-regulatory element binding protein (SREBP)-1 through the PI3K-mTOR complex 1 pathway (21). SREBPs belong to the basic helix-loop-helix-leucine zipper family of transcription factors, which include SREBP-1a, SREBP-1C and SREBP-2. These transcription factors have been reported to be regulators of the activation or expression of enzymes in lipid and cholesterol homeostasis (33). Such enzymes include lanosterol-14α-demethylase, squalene epoxidase monooxygenase, sterol-Δ7-reductase, farnesyl-diphosphate farnesyl transferase 1 and sterol-C4-methyl oxidase, encoded for by CYP51A1, SQLE, DHCR7, FDFT1 and SC4MOL, respectively, which were genes identified in the present study to be involved in the top biological process affected by FGFR3 knockdown in BC, steroid biosynthesis. When sterol concentrations are low, SREBPs bind with SREBP cleavage-activating protein and mature SREBPs activate the biosynthesis of sterol (34). Wu et al (35) and Degener et al (36) revealed an association between the metabolism of steroid or bile acid and bladder cancer. Therefore, it was hypothesized that FGFR3 may regulate the corresponding target genes, which in turn affects the metabolism and biosynthesis of steroid substances and subsequently affects BC.

The genes with the ten highest degrees of interaction in the PPI network constructed in the present study were SETX, CFD, HIST2H2BE, NME2, SERPINB2, UBC, CKS, DHRS2, TRIP6 and MCM5. NME2, also known as NM23-H2, encodes nucleoside diphosphate kinase (NDPK)-B, which catalyzes the transposition of γ-phosphate between nucleosides (37). In addition, NME2 is known to be a motility and metastasis suppressor (38), which may inhibit cancer, as cells must survive and proliferate to become overt metastases. NME1, also known as NM23-H2, is a paralog of NME2, with 88% amino acid identity, and was first discovered to be a metastasis suppressor by Steeg et al (39) in 1988. Yong et al (40) studied the differential expression of NM23-H1 in BC and normal bladder cases using the immunohistochemical technique streptavidin-peroxidase procedure and demonstrated that the positive expression rates of NM23-H1 were 62.3 and 100.0% in BC and normal bladder, respectively. With disease progression, the positive expression rate decreased indicating its important role in BC (40). NM23 phosphorylates kinase suppressor of Ras and prevents downstream activation of the MAPK pathway (41). Therefore, FGFR3 expression may lead to the activation of the MAPK pathway. Furthermore, numerous studies have provided evidence to suggest that the NME2 gene may be associated with cancer (42,43).

Co-expression genes were hypothesized to have functional associations, such as physical interactions, between the encoded proteins. In the present study, the gene co-expression network containing 168 differentially co-expressed genes and 1,935 associations were built. The genes with the highest degrees of interaction were CCNB1, MCM3, TMEM97, MCM5, H2AFZ, PPIL1, UBE2T, MCM2, CTSF and ZIC2. CCNB1, encoding cyclin B1, was reported to contribute to the regulation of G2-M-phase transition, which is essential for DNA synthesis and cell proliferation (44). Dysregulated expression of CCNB1 may therefore result in uncontrolled growth and malignant transformation (45). Of note, CCNB1 is one of the 11 genes to predict outcome in several types of cancer, including BC (46). Yuan et al (47) reported that the specific downregulation of CCNB1 may lead to tumor regression through preventing the progression of cells in G2 phase and triggering cell death. H2AFZ encodes H2A histone family, member Z, which is a variant of histone H2A and is responsible for the thermosensory response and regulating euchromatin-heterochromatin transition (48). H2AFZ may have a role in high-grade cancer due to its ability to regulate a large numbers of genes (49). Dong et al (50) revealed that H2AFZ was overexpressed in BC and may be applied to the diagnosis of BC.

In conclusion, according to the expression profile data of FGFR3 knockdown in BC, the present study identified 196 DEGs. GO analysis as well as KEGG and BIOCATRA pathway enrichment analysis revealed that the primary biological process in which these DEGs were involved was sterol biosynthesis and metabolism. In addition, PPI networks of DEGs and co-expressed genes were constructed and revealed the information flow of PPIs. This comprehensive expression profile data of BC provided novel insight into the pathogenesis and occurrence of BC associated with FGFR3.

References