Abstract
Aims: The aim of this study was to screen several novel genes associated with renal cell carcinoma (RCC), and analyze the gene functions and signal pathways which were critical to RCCs with DNA microarray. Methods: The gene expression profile of GSE781 was downloaded from Gene Expression Omnibus database, including 9 RCC samples and 9 healthy controls. Compared with the control samples, differentially expressed genes (DEGs) of RCC was identified the by packages in R. The selected DEGs were further analyzed using bioinformatics methods. Gene ontology (GO) enrichment analysis was performed using Gene Set Analysis Toolkit and protein-protein interaction (PPI) network was constructed with prePPI. Then, pathway enrichment analysis to PPI network was performed using WebGestalt software. Results: A total of 429 DEGs were down-regulated and 418 DEGs were up-regulated in RCC samples compared to healthy controls. A total of 11 remarkable enhanced functions and 13 suppressed functions were identified. PPI nodes of high degrees, such as JAK2, IL8, BMPR2, FN1 and NCR1, were obtained. The DEGs were classified and significantly enriched in cytokine and cytokine receptor pathway. Conclusion: The hub genes we find from RCC samples are not only bio-markers, but also may provide the groundwork for a combination therapy approach for RCCs.
Introduction
Renal cell carcinoma (RCC), also known as hypernephroma, is a kidney cancer that originates in the lining of the proximal convoluted tubule, and known to be the most lethal of all the genitourinary tumors [1,2]. It is also the most common type of kidney cancer in adults, responsible for approximately 80% of cases [3]. RCC is relatively resistant to radiation therapy and chemotherapy, although some cases respond to immunotherapy. The major difficulty in RCC is the constitution that this disease is not one entity but rather a collection of different types of tumors, each possessing distinct genetic characteristics, histological features, and, to some extent, clinical phenotypes [4].
Recent research has demonstrated several key factors associated with RCC initiation and multiorgan metastasis, such as von Hippel-Lindau (VHL), BRCA1 associated protein-1 (BAP1) and polybromo 1 (PBRM1) [5,6,7,8]. Inactivation of the VHL tumor suppressor gene is an archetypical tumor-initiating event in clear cell RCC (ccRCC) that leads to the activation of hypoxia-inducible transcription factors (HIFs), which have opposing effects on clear cell carcinogenesis [9,10]. VHL-HIF response maintains an epigenetically altered which is specific to metastatic ccRCC in VHL-deficient cells [11]. Focusing on the two most prominent pro-metastatic VHL-HIF target genes, it has been shown that loss of polycomb repressive complex 2 (PRC2)-dependent histone H3 Lys27 trimethylation (H3K27me3) activates HIF-driven chemokine (C-X-C motif) receptor 4 (CXCR4) expression in support of chemotactic cell invasion, whereas loss of DNA methylation enables HIF-driven cytohesin 1 interacting protein (CYTIP) expression to protect cancer cells from death cytokine signals [12]. As a tumor suppressor gene, PBRM1 is a critical regulator of p21 and a physiologic mediator of G1 arrestment [13]. BAP1 is identified as one of several putative two-hit tumor suppressor genes in a two-hit manner though a whole-genome and exome sequencing [14].
However, the paramount regulator and distinct molecular mechanism of RCC has yet to be evaluated as the complexity of its pathogenesis. In this present paper, we aimed to screen several novel tumor suppressor genes and explore the molecular mechanism of RCC with a computational analysis. The most significant DEGs which served as potential biomarkers may provide a new sight on RCC clinical therapy.
Materials and Methods
Affymetrix microarray
The gene profile of GSE781 was downloaded from Gene Expression Omnibus (http://www.ncbi.nlm.nih.gov/geo/), which was based on the platform of Affymetrix Human Genome U133B Array condensed. This expression dataset was deposited by Lenburg ME [15]. A total of 17 samples were derived from 9 RCC patients and 8 age- and gender- matched healthy control subjects.
Data preprocessing and differential expression analysis
We first converted the probe-level data in CEL files 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. And then, we imputed missing data and performed quartile data normalization [16,17]. The multtest package in R was used to identify differentially expressed genes (DEGs) in patients with RCC compared to healthy controls [18]. To circumvent the multi-test problem which might induce too much false positive results, the Benjamini & Hochberg (BH) method [19] was used to adjust the raw P-values into false discovery rate (FDR). The FDR less than 0.05 and the absolute logFC larger than 1 were chosen as cut-off criteria.
Gene ontology enrichment analysis
To produce a dynamic, controlled vocabulary that can be applied to all eukaryotes, Gene Ontology (GO) analysis has been used frequently in functional study large scale genomic and transcriptional data [20,21]. DEGs were separated into two sets based on different expression behavior, and then the GO analysis was performed using Gene Set Analysis Toolkit suit [22,23]. A p-value less than 0.05 was considered statistically significant.
Construction and analysis of the protein - protein interaction (PPI) network
PrePPI is a database of predicted and experimentally determined protein-protein interactions (PPI) using a Bayesian framework that combines structural, functional, evolutionary and expression information for yeast and human [24]. In this analysis, the most significant up- or down-regulated DEGs were screened, and a PPI network was constructed by collecting interactions from prePPI database [25,26].
Pathway enrichment analysis
For functional annotation the DEGs, the over-represented KEGG categories in pathways was identified using KAAS [27]. The p-value less than 0.05 was chosen as cut-off criterion.
Results
Differential gene expression in RCC samples
After normalization(Figure 1), the data was performed differential expression analysis. The genes with p-value less than 0.05 and logFC absolute value over than 1, were considered as DEGs. A total of 429 DEGs were down-regulated and 418 DEGs were up-regulated in RCC samples compared to healthy controls. Among these 848 DEGs, fatty acid binding protein 7 (FABP7) was the most significantly up-regulated gene and uromodulin (UMOD) was down-regulated the most (p=3.4E-04, logFC=6.93; p=3.9E-08, logFC=-8.34, respectively). In addition, cytokine receptors, such as interleukin 15 receptor alpha (IL15RA) and colony stimulating factor 1 receptor (CSF1R), were also up-regulated (p=3.8E-04, logFC=2.21; p=4.9E-05, logFC=1.96, respectively).
Protein-protein interaction network
Interaction networks of DEGs were constructed as shown in Figure 2. According to the results, interaction relationship of 207 DEGs was obtained. A high level of connectivity was shown in the network constructed based on UMOD (Figure 2A). PPI nodes janus kinase 2 (JAK2), interleukin 8 (IL8), bone morphogenetic protein receptor, type II (BMPR2), fibronectin (FN1) and natural cytotoxicity triggering receptor 1 (NCR1) were of high degrees in the network, all of which were obtained from RCC samples. In addition, the PPI network based on FABP7 was also constructed, with imperfect connections due to the absence of several key hubs, as is shown in Figure 2B.
Immune response enhancement and hormone suppression in RCCs
We performed a functional enrichment analysis using Gene Set Analysis Toolkit and identified 11 remarkable enhanced functions and 13 suppressed functions (Figure 3, FDR<0.05). The DEGs was classified with GO analysis according to their functions. Among these enhanced functions, the most remarkable up-regulation was immune response based on cytokine and cytokine receptor enrichment (GO: 0006955, FDR= 5.8E-29). The other significant functions included defense response (GO: 0006952, FDR= 6.3E-14), response to wounding (Go: 0009611, FDR=1.5E-11), leukocyte activation (GO: 0045321, FDR= 1.2E-08), positive regulation of immune system process (GO: 0002684, FDR= 2.2E-10), inflammatory response (GO: 0006954, FDR=1.2E-05) and so on. Genes, such as toll-like receptor (TLR) family, IL super-family, chemokine (C-X-C motif) ligand 10 (CXCL10), TAP binding protein (TAPBP), C-type lectin domain family 4 member A (CLEC4A), major histocompatibility complex class I (HLA), complement component 3a receptor 1 (C3AR1) and integrin alpha L (ITGAL), were involved in GO analysis.
The most evidently downward modulation function in RCC was oxidation reduction which regulated by hormones (GO: 0055114, FDR= 3.9E-36). The other significantly suppressed function included ion transport (GO: 0006811, FDR=9.9E-07), response to organic substance (GO: 0010033, FDR=3.2E-03) and generation of precursor metabolites and energy (GO: 0006091, FDR=7.9E-13). The response to hormone stimulus was also suppressed (GO: 0009725, FDR=5.9E-03), which in turn weaken the chemical homeostasis.
Cytokine and cytokine receptor pathway enrichment in RCCs
Pathway enrichment analysis of all genes involved in PPI was performed with KAAS (Figure 4). Only cytokine and cytokine receptor pathway was enriched with a remarkable FDR less than 0.001 in PPI network of UMOD. Genes associated with cytokine and cytokine receptor pathway were chemokine, hematopoietin, PDGF family, TNF family and TGF-β family proteins genes. Five chemokines genes participated in this pathway, including IL8, CXCL9, CCL5, CCL3, CCL6 and CCL9. The primary hematopoietins were cytokine receptors, such as OSMR, CNTFR and CSF3R. TGF-β family proteins, such as BMP7, BMPR2, TGFB2, TGFB3 and their receptor TGFBR1, were also involved in this pathway.
Discussion
In the present paper, we identified 848 DEGs in RCC patient samples compared to the control samples. DEGs in RCCs were highly associated with immune response, hormone response and defense response, which may play important roles in tumor initialization and migration. Pathway enrichment demonstrates that only cytokine and cytokine receptor pathway was enhanced in RCC patients compared to healthy controls. Genes obtained from the results of pathway enrichment analysis are the subject of our investigation. Among these DEGs, UMOD and FABP7 were the most significantly down-regulated and up-regulated genes, respectively. In addition, cytokine receptors, such as IL15RA and CSF1R, were also up-regulated.
UMOD was the most abundant protein in urine. Previous studies have shown that mutation of UMOD was highly associated with nephronophthisis (NPHP), medullary cystic kidney disease (MCKD), and familial juvenile hyperuricemic nephropathy (FJHN) [28,29]. This protein severed as a receptor for binding and endocytosis for tumor necrosis factor (TNF) and cytokines [30]. In the RCC patients, UMOD was significantly down-regulated, which is accordant to that of chronic kidney disease patients in previous studies [28,31].
Recent research has shown that fatty acid binding protein 7 (FABP7), also known as B-FABP7, was 20 folds higher in RCC than that in healthy donors, and was undetectable in patients with RCC after surgical operation [32,33]. Thus, FABP7 was chosen as one of the markers for RCC [32,34]. The expression of FABP7 was controlled by OCT1, OCT6 and nuclear factor I (NFI), and supershift experiment suggested that BRN2 may bind to the promoter of FABP7 as well [35]. Further studies demonstrated that FABP7 can be regulated by protein kinase C (PKC) and mitogen-activated protein kinase (MAPK) pathways [36]. In addition, FABP7 was also high expressed in cutaneous malignant melanoma and breast cancers [37,38].
Cytokine and cytokine receptor pathway was enriched in RCC patients. However, different cytokines have distinct functions in tumors. IL 15, a proinflammatory cytokine and growth factor, was produced by macrophages and dendritic cells, and can promote nature killer T cell leukaemias [39]. Clinical treatment has applied IL15 as an ancillary drug for cancer therapy [40]. In our data, IL15RA was highly enriched in RCC patients, which suggested a function enhancement of IL15. This enrichment may help patients in promoting innate and adaptive immunoreactions, and enhancing cytotoxicities. It has been reported that granulocyte macrophage colony-stimulating factor (GM-CSF) can inhibits lymphomas and carcinomas, but macrophage colony-stimulating factor (M-CSF) may promotes breast cancer invasion, even both of them bound to the same receptors [41,42]. The enrichment of CSF1R in our study was accordant to previous study, though its function remains unclear [43].
Even though the DEGs and hub genes may provide a new sight into the thrapy approach for RCCs, however, the expression of these genes was not confirmed by real-time PCR and the function in RCCs was not evaluated. The potential biomarkers and hub genes, such as JAK2, IL8, BMPR2, FN1 and NCR1, only provided potential targets for RCC therapy, and further effort to confirm the hypothesis was still in great needed.
Conclusion
We have demonstrated 847 DEGs in RCC samples compared to the control samples. Cytokine and cytokine receptor pathway was enriched in RCC samples, which suggested a role in the development of RCC for this signal pathway. The DEGs and hub genes we find from RCC samples are not only bio-markers, but also may provide the groundwork for a combination therapy approach for RCCs.
Conflicts of Interest
The authors have no conflict of interest to state.