Silencing TAK1 alters gene expression signatures in bladder cancer cells.

The aim of the present study was to identify the differentially expressed genes (DEGs) that are induced by the silencing of transforming growth factor-β-activated kinase 1 (TAK1) in bladder cancer cells and to analyze the potential biological effects. Dataset GSE52452 from mutant fibroblast growth factor receptor 3 (FGFR3) bladder cancer cells transfected with control siRNA or TAK1-specific siRNA was downloaded from Gene Expression Omnibus. The DEGs between the two groups were identified using Limma package following data pre-processing by Affy in Bioconductor. Enrichment analysis of DEGs was performed using the Database for Annotation, Visualization and Integrated Discovery, followed by functional annotation using TRANSFAC, TSGene and TAG databases. Integrated networks were constructed by Cytoscape and sub-networks were extracted employing BioNet, followed by enrichment analysis of DEGs in the sub-network. A total of 43 downregulated and 21 upregulated genes were obtained. The downregulated genes were enriched in five pathways, including NOD-like receptor signaling pathway and functions related to cellular response. The upregulated genes were associated with cellular developmental processes. Transcription factor EGR1 and 9 tumor-associated genes were screened from the DEGs. Among the DEGs, 10 hub nodes may represent important roles in the complex metabolic network, including EGFR, CYP3A5, MAP3K7, GSTA1, PTHLH, ALDH1A1, KCND2, EGR1, ARRB1 and ITPR1. Additionally, EGFR was correlated with ERBB2, GRB2 and PIK3R1, and these were enriched in ErbB signaling pathway and various cancer-associated pathways. Silencing TAK1 may decrease cellular response to chemical stimulus via downregulating CYP3A5, MAP3K7, GSTA1, ALDH1A1, ARRB1 and ITPR1; increase cancer cell development via upregulating EGFR, EGR1 and PTHLH; and regulate cancer metastasis through EGFR, ERBB2, GRB2 and PIK3R1.

Introduction

Bladder cancer may be classified into two stages: Non-muscle-invasive disease of low grade, and progressive muscle-invasive disease of high grade, which may further deteriorate into metastatic cancer (1). This type of cancer predominantly occurs as urothelial cell carcinoma and ranks as the fourth most common malignancy and the eighth most common cause of cancer-associated mortality in men (2). Overexpressed or mutant fibroblast growth factor receptor 3 (FGFR3) may contribute to the development of the transformed phenotype of urothelial carcinoma (3), and the targeted inhibition of FGFR3 may thus prevent the development of superficial bladder cancer (4). The transcription factor nuclear factor κB (NFκB) is an important mediator of the angiogenesis and metastasis of bladder cancer (5). Activating mutations in FGFR3 can promote NFκB transcriptional activity through the important mediator transforming growth factor-β-activated kinase 1 (TAK1), which exerts a positive regulatory effect on the activity of NFκB in urothelial cell carcinoma (6). However, the exact regulatory effect of TAK1 on the progression of bladder cancer remains unclear.

RNA interference may effectively aid in defining the roles of specific genes in the progression of cancer when specifically designed siRNAs are used to silence target genes (7,8). Using the microarray expression data of bladder cancer cells transfected with control siRNA or TAK1-specific siRNA which were deposited by Salazar et al (6), the present study aimed to identify the differentially expressed genes (DEGs) and to screen for tumor-associated DEGs, followed by identification of the biological processes or pathways implicated by DEGs and the hub nodes in the complex protein-protein interaction (PPI) network and sub-network, in an attempt to provide a deeper insight into the molecular mechanisms mediating the effect of silencing TAK1 on bladder cancer.

Materials and methods

Gene expression profiles

The gene expression profiles (accession number GSE52452) (6) from 6 samples of MGHU3 (Y375C) mutant FGFR3 bladder cancer cells that were transfected with control siRNA (GSM1267159-61) or TAK1-specific siRNA (GSM1267150-52), were downloaded from the public functional genomics data repository of Gene Expression Omnibus (9). These 6 samples (3 control samples vs. 3 TAK1 siRNA-treated samples) were previously researched using Affymetrix Human Gene 1.0 ST Array.

Data processing and screening of DEGs

The raw expression data were processed through RMA background correction (10), quantile normalization, log base 2 (log2) transformation and probeset summarization to obtain the gene expression matrix by employing Affy package in Bioconductor (11) and probe annotation files provided by Brain Array Lab (brainarray.mbni.med.umich.edu/). Subsequently, the lists of DEGs between control and TAK1-specific siRNA-treated cells were generated under the thresholds of |log2 fold change (FC)|≥1 (fold-change magnitude, >2) and P-value <0.05 using the empirical Bayes method offered by Limma package (12).

Functional enrichment analysis of DEGs

To provide insight into the biological functions or pathways involving the identified DEGs, the Gene Ontology (GO) Biological Process (13) and Kyoto Encyclopedia of Genes and Genomes (KEGG) (14) pathway enrichment analyses of DEGs were performed with the criterion of P<0.05 using the Database for Annotation, Visualization and Integrated Discovery (15).

Functional annotations of DEGs

To identify the specific functions of DEGs, transcription factors were screened from the DEGs based on the TRANSFAC database (16), accompanied by screening of tumor-suppressor genes (TSGs) using the TSGene (17) database and identification of oncogenic genes using the tumor-associated genes (TAG) database (18).

Construction of PPI network and screening of sub-network

In consideration of the gene-gene interactions in the complex biological systems, the identified DEGs were inputted into the Search Tool for the Retrieval of Interacting Genes (STRING) database (19) to identify the interacting pairs with a combined score >0.9 (data downloaded on May 9, 2014). The PPI network was constructed using Cytoscape software (20), followed by excavation of the sub-network using the BioNet tool (bionet.bioapps.biozentrum.uni-wuerzburg.de/). KEGG pathway enrichment analysis of DEGs in the sub-network was performed to identify the associated biological pathways.

Results

Screened DEGs

With the criteria of |log2 FC|≥1 and P<0.05, a total of 43 downregulated genes and 21 upregulated genes were identified in the mutant FGFR3 bladder cancer cells transfected with TAK1-specific siRNA in comparison with the cells transfected with control siRNA.

Functional enrichment results of DEGs

By performing a KEGG pathway enrichment analysis of the separate upregulated and downregulated genes, it was demonstrated that the downregulated genes were significantly enriched in five pathways, including the NOD-like receptor signaling pathway [P=0.014731; mitogen-activated protein kinase kinase kinase 7 (MAP3K7) and tumor necrosis factor α-induced protein 3 (TNFAIP3)] and the upregulated genes were enriched in the cytokine-cytokine receptor interaction pathway [P=0.018547; epidermal growth factor receptor (EGFR) and interleukin 18 (IL18)] (Table I). The GO Biological Process enrichment analysis revealed that the downregulated genes were associated with cellular responses, including regulation of catalytic activity [P=5.45×10−4; e.g. aldehyde dehydrogenase 1 family member A1 (ALDH1A1), β-arrestin (ARRB1), inositol trisphosphate receptor (ITPR1) and MAP3K7)] and cellular response to chemical stimulus [P=1.38×10−3; e.g. ALDH1A1, ARRB1, ITPR1, MAP3K7, cytochrome P450 family 3 subfamily A polypeptide 5 (CYP3A5) and glutathione S-transferase α1 (GSTA1)]; the upregulated genes were associated with cellular developmental processes [P=5.80×10−4; e.g. EGFR, early growth response 1 (EGR1) and parathyroid hormone-like hormone (PTHLH)] (Table II).

Table I.

Significantly enriched KEGG pathways of upregulated and downregulated genes.

KEGG pathway	P-value	Gene list
Downregulated genes
Aldosterone-regulated sodium reabsorption	0.0079126	SCNN1G, SGK1
NOD-like receptor signaling pathway	0.014731	MAP3K7, TNFAIP3
Retinol metabolism	0.0177602	ALDH1A1, CYP3A5
Metabolism of xenobiotics by cytochrome P450	0.0215999	CYP3A5, GSTA1
Drug metabolism-cytochrome P450	0.0227555	CYP3A5, GSTA1
Upregulated genes
Cytokine-cytokine receptor interaction	0.018547	EGFR, IL18

KEGG, Kyoto Encyclopedia of Genes and Genomes.

Table II.

Top 10 significantly enriched GO terms of upregulated and downregulated genes.

GO term	P-value	Gene list
Downregulated genes
GO:0006805 xenobiotic metabolic process	5.32E-04	ACSL1, ALDH1A1, CYP3A5, GSTA1
GO:0050790 regulation of catalytic activity	5.45E-04	ACER2, ACSL1, ALDH1A1, ARRB1, DUSP10, FGD3, ITPR1, MAP3K7, SGK1, SYTL2, TNFAIP3
GO:0071466 cellular response to xenobiotic stimulus	5.45E-04	ACSL1, ALDH1A1, CYP3A5, GSTA1
GO:0009410 response to xenobiotic stimulus	5.85E-04	ACSL1, ALDH1A1, CYP3A5, GSTA1
GO:0045824 negative regulation of innate immune response	5.89E-04	DUSP10, TNFAIP3
GO:0032715 negative regulation of interleukin-6 production	9.53E-04	ARRB1, TNFAIP3
GO:0065009 regulation of molecular function	1.07E-03	ACER2, ACSL1, ALDH1A1, ARRB1, BAMBI, DUSP10, FGD3, ITPR1, MAP3K7, SGK1, SYTL2, TNFAIP3
GO:0070887 cellular response to chemical stimulus	1.38E-03	ACER2, ACSL1, ALDH1A1, ARRB1, BAMBI, CYP3A5, FGD3, GDF15, GSTA1, ITPR1, MAP3K7, TNFAIP3
GO:0002819 regulation of adaptive immune response	2.02E-03	DUSP10, MAP3K7, TNFAIP3
GO:0018149 peptide cross-linking	2.08E-03	SPRR1A, SPRR3
Upregulated genes
GO:0036296 response to increased oxygen levels	1.01E-06	EGR1, IL18, PDPN
GO:0055093 response to hyperoxia	1.01E-06	EGR1, IL18, PDPN
GO:0030324 lung development	2.45E-05	HEG1, IL18, PDPN, PTHLH
GO:0030323 respiratory tube development	2.64E-05	HEG1, IL18, PDPN, PTHLH
GO:0060541 respiratory system development	4.12E-05	HEG1, IL18, PDPN, PTHLH
GO:0001945 lymph vessel development	2.81E-04	HEG1, PDPN
GO:0060571 morphogenesis of an epithelial fold	4.59E-04	EGFR, PTHLH
GO:0048869 cellular developmental process	5.80E-04	ANTXR1, EGFR, EGR1, HEG1, IL18, MDK, MEA1, PDPN, PTHLH, STRADB
GO:0009725 response to hormone stimulus	6.92E-04	EGFR, EGR1, IL18, MDK, STRADB
GO:0071320 cellular response to cAMP	7.60E-04	EGR1, IL18

GO, Gene Ontology.

Results of functional annotation of DEGs

Using TRANSFAC, TSGene and TAG databases, the functional annotation analysis screened 1 transcription factor (upregulated EGR1) and 9 TAGs, including 1 oncogene (upregulated EGFR), 5 TSGs (downregulated TNFAIP3, ras association domain family member 5 and carcinoembryonic antigen related cell adhesion molecule 7; upregulated GLI pathogenesis related 1 and EGR1) and 3 other TAGs (downregulated breast carcinoma amplified sequence 1 and activin membrane-bound inhibitor; upregulated PTHLH).

Analysis of the PPI network

Using the STRING database and Cytoscape tool, a PPI analysis of the DEGs was performed and an integrated PPI network was obtained. Based on this network, the top 10 hub nodes representing important roles in the complex metabolic networks were EGFR (degree, 178), CYP3A5 (degree, 83), MAP3K7 (degree, 58), GSTA1 (degree, 56), PTHLH (degree, 47), ALDH1A1 (degree, 45), potassium voltage-gated channel subfamily D member 2 (degree, 43), EGR1 (degree, 42), ARRB1 (degree, 36) and ITPR1 (degree, 33) (Fig. 1). Additionally, this network revealed a correlation of MAP3K7 with TNFAIP3.

Figure 1.

Protein-protein interaction network. Red represents upregulation, green represents downregulation, and yellow represents no difference.

Analysis of the screened sub-network

Using the BioNet tool, a sub-network with EGFR at the core was screened from the PPI network, in which EGFR was associated with v-erb-b2 avian erythroblastic leukemia viral oncogene homolog 2 (ERBB2), growth factor receptor-bound protein 2 (GRB2) and phosphoinositide-3-kinase, regulatory subunit 1 (PIK3R1) (Fig. 2). Furthermore, the DEGs in the sub-network were found to be enriched in various pathways associated with cancer, including the ErbB signaling pathway (P=3.07×10−7; e.g. EGFR, ERBB2, GRB2 and PIK3R1), pathways in cancer (P=5.90×10−7; e.g. EGFR, ERBB2, GRB2 and PIK3R1), prostate cancer (P=1.04×10−4; e.g. EGFR, ERBB2, GRB2 and PIK3R1), endometrial cancer (P=1.64×10−4; EGFR, ERBB2, GRB2 and PIK3R1), non-small cell lung cancer (P=1.90×10−4; EGFR, ERBB2, GRB2 and PIK3R1) and glioma (P=3.90×10−4; e.g. EGFR, GRB2 and PIK3R1) (Table III).

Figure 2.

Extracted sub-networks. Color depth represents fold-change magnitude of differential expression. Red represents upregulation and green represents downregulation. Squares represent genes with low importance in the sub-network. Circles represent genes with high importance in the sub-network.

Table III.

Top 10 significantly enriched KEGG pathways of differentially expressed genes in the sub-network.

KEGG pathway	P-value	Gene list
ErbB signaling pathway	3.07E-07	EGFR, ERBB2, GRB2, NRG1, HBEGF, ERBB3, PIK3R1, EGFR, ITGA6, PTGS2, JUP, IGF1R,
Pathways in cancer	5.90E-07	MET, FOS, IL8, ERBB2, GRB2, PIK3R1
Focal adhesion	7.34E-07	EGFR, ITGA6, CAV1, IGF1R, MET, VAV2, ERBB2, GRB2, PIK3R1
Hepatitis C	5.77E-06	EGFR, STAT2, IRF1, IL8, GRB2, OAS1, PIK3R1
Endocytosis	8.07E-05	EGFR, ARRB1, AP2M1, CAV1, IGF1R, MET, ERBB3
Prostate cancer	1.04E-04	EGFR, IGF1R, ERBB2, GRB2, PIK3R1,
Malaria	1.52E-04	MET, IL8, SDC3, SDC4
Endometrial cancer	1.64E-04	EGFR, ERBB2, GRB2, PIK3R1
Non-small cell lung cancer	1.90E-04	EGFR, ERBB2, GRB2, PIK3R1
Glioma	3.90E-04	EGFR, IGF1R, GRB2, PIK3R1

KEGG, kyoto encyclopedia of genes and genomes.

Discussion

By re-analyzing the microarray data from specific mutant FGFR3 bladder cancer cells using bioinformatic methods, the present study identified 64 genes with significantly altered expression between TAK1-specific siRNA-treated cells and the control siRNA-treated cells, including 9 TAGs. The downregulated genes were related to the NOD-like receptor signaling pathway and cellular response functions, whereas the upregulated genes were associated with cellular developmental processes. By constructing a PPI network, the present study identified 10 hub nodes that may exert the predominant effects on the network. Furthermore, the DEGs in the sub-network with EGFR at the core were associated with various types of cancer and the ErbB signaling pathway.

TAK1 siRNA vs. control siRNA samples yielded 43 downregulated genes, of which 6 were identified to be hub nodes in the PPI network, including CYP3A5, MAP3K7, GSTA1, ALDH1A1, ARRB1 and ITPR1. CYP3A5 and GSTA1 are associated with the detoxification of chemical stimuli, electrophilic compounds or other damage stimuli (21,22). They were also enriched in the term metabolism of xenobiotics by cytochrome P450, which plays important roles in detoxification (23). MAP3K7, which correlated with TNFAIP3, was enriched in the NOD-like receptor signaling pathway which is involved in sensing intracellular microbial motifs or other damage stimuli (24). ALDH1A1 and CYP3A5 were enriched in retinol metabolism; the loss of retinol acyltransferase is inversely correlated with the invasion bladder cancer (25). ARRB1 is considered to cause a specific dampening cellular response to stimuli or sensory signals (26) and ITPR1 mediates calcium release which then amplifies apoptosis in response to specific stimuli (27,28). Accordingly, the present study further demonstrated that the 6 hub nodes in the PPI network were enriched in the GO function cellular response to chemical stimulus. It is therefore possible that silencing of TAK1 using specific siRNA transfection may alter the cellular response-associated pathways or functions in response to stimuli via regulating the expression of the 6 DEGs.

TAK1 siRNA treatment induced 21 upregulated genes in bladder cancer cells, among which EGFR, PTHLH and EGR1 were identified as hub nodes in the PPI network. These 3 genes were all enriched in cellular developmental processes. As exemplified, activated EGFR, which contributes to phenotypic characteristics in various tumor types (29), is an effective therapeutic target for the treatment of bladder cancer (30,31). Increased expression of PTHLH resulting from downregulated p38MAPK signaling is associated with metastatic lesions to the liver and lung from colon cancer cells (32). EGR1 also exerts pro-tumorigenic effects by contributing to tumor infiltration, node formation and metastasis (33). Thus, the upregulation of these genes in TAK1 siRNA-treated bladder cancer cells may suggest that silencing TAK1 has the effect of promoting cancer cell developmental processes. Furthermore, the present study also predicted EGFR, EGR1 and PTHLH as TAGs, which suggests their potential use as therapeutic targets for the diagnosis and treatment of bladder cancer.

Using BioNet software, a sub-network with EGFR at the core was screened from the PPI network. According to the PPI network, EGFR was correlated with ERBB2, GRB2 and PIK3R1. These 4 DEGs were enriched in various cancer types, including prostate cancer, endometrial cancer, non-small cell lung cancer and glioma. As previously reported, the EGFR family of four receptors including EGFR and ERBB2, is implicated in the development and progression of various human cancer types (34). The activation of ERBB2 may result in resistance to cetuximab-based therapy targeting EGFR, while the inhibition of this gene can restore cetuximab sensitivity in patients with cetuximab-resistant cancers (35). GRB2 amplification is observed in esophageal squamous cell carcinoma and is significantly involved in lymph node metastases (36). Phosphatidylinositol 3-kinase (PI3K), which may promote cancer cell survival, is an important therapeutic target in cancer (37). PIK3R1, the inhibitory subunit of PI3K, often undergoes mutations in endometrial cancer (38,39). The enrichment results also revealed the these 4 DEGs were enriched in ErbB signaling pathway, which is involved in regulating cell survival and adhesion (40). Silencing TAK1 in bladder cancer cells led to the upregulation of EGFR and downregulation of ERBB2, GRB2 and PIK3R1, suggesting an inhibitory effect of silencing TAK1 on the metastasis of bladder cancer cells via modulation of the ErbB signaling pathway. However, the explicit molecular mechanisms require further research.

In summary, the present study demonstrated that silencing TAK1 induced 43 downregulated and 21 upregulated genes. Silencing TAK1 may lead to decreased cellular response to chemical stimuli via downregulating CYP3A5, MAP3K7, GSTA1, ALDH1A1, ARRB1 and ITPR1, as well as increased cancer cell developmental processes via upregulating EGFR, EGR1 and PTHLH. In addition, silencing TAK1 may exert regulatory effects on bladder cancer metastasis and other various cancer types via regulating the expression of EGFR, ERBB2, GRB2 and PIK3R1.