Introduction
Diabetic nephropathy (DN) is a chronic kidney disease and is a serious complication of long-term diabetes mellitus (1). DN develops in 30–40% of patients with type I and II diabetes mellitus and is a risk factor for increased mortality in patients with cardiovascular disease (2). Although the signs of early DN are not evident, the clinical evidence for DN is the presence of proteinuria, glomerular hypertrophy, decreased glomerular filtration and a decline in renal function (1). DN continues to present a health concern worldwide.
Accumulating experimental in vivo and in vitro evidence has indicated that multiple pathways and cytokines play a role in the pathogenesis of DN. For example, a recent study suggested a cardinal role of inflammatory molecular and pathways in the pathogenesis of DN (3). The activation of the innate immune response associated with various inflammatory molecules, such as interleukin (IL)-1, IL-18 and tumor necrosis factor (TNF) has also been shown to contribute to the renal injury observed in patients with DN (4). Furthermore, it has been reported that the nuclear factor (NF)-κB signaling pathway induces the expression of inflammatory genes during the progression of DN, and these effects are modulated by the Ras homolog gene family, member A (RhoA)/Rho-associated protein kinase (ROCK) signaling pathway (5). A good understanding of the molecular mechanisms responsible for the disease may aid in the development of effective therapies. However, the molecular mechanisms of DN have not yet been fully clarified.
Microarray data have been widely used to connect genes and molecules to diseases (6). Reiniger et al proved the target role of receptor for advanced glycation end-products (RAGE) in the treatment of DN based on microarray data (Accession no. GSE20844) (7). In the present study, we downloaded the same microarray data from the Gene Expression Omnibus (GEO) database. Subsequently, based on the gene expression profiles, the differentially expressed genes (DEGs) were analyzed and the DEG-related functions and pathways were predicted. The aim of the present study was to elucidate the mechanisms of DN pathogenesis and to identify associated significant genes.
Data collection methods
Data acquisition and preprocessing
Whole-genome micro-array gene expression data for glomeruli from diabetic male OVE26 mice (diabetic group, n=4) and glomeruli from non-diabetic male FVB mice (control group, n=3) have been deposited in the GEO archive database (Accession no. GSE20844) (7). We downloaded the raw Affymetrix CEL files based on the platform of Affymetrix Mouse Genome 430 2.0 Array.
The raw data underwent pre-processing, including background correction, quantile normalization and probe summarization with the application of bioconductor package 'affy', as previousy described (8).
Analysis of DEGs
The DEGs in the diabetic group compared with the non-diabetic controls were analyzed using the Bioconductor package 'limma', as previousy described (9). The P-value for each gene was calculated using the Student's t-test. Genes with differences in expression denoted by values of p<0.05 and |log2FC (fold change)|≥0.58, screened as DEGs. In order to compare the differences in the profiles of DEGs between the diabetic and control samples, the gene expression data was clustered using R gplots software package (http://cran.r-project.org/web/packages/gplots/index.html.). Subsequently, the chromosomal location of the DEGs was explored based on the chip annotation information.
Gene Ontology (GO) and pathway analysis
The upregulated and downregulated genes were subjected to GO and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis. The genes were enriched in three GO categories, such as biological process (BP), molecular function (MF) and cellular component (CC). The enrichment analysis based on the hypergeometric distribution was performed using the Database for Annotation, Visualization and Integrated Discovery (DAVID) online tool (10). The cut-off value for a significant GO term and pathway was set to P<0.05 and count ≥2.
Protein-protein interaction (PPI) network
The functional protein interactions among the DEGs and the encoding proteins were predicted using the Search Tool for the Retrieval of Interacting Genes (STRING) (11). The PPI score was set as 0.4 and other parameters were set as the default value. Cytoscape was used to visualize the PPI network.
Subsequently, with the application of CytoNCA (12), the hub nodes were measured based on the degree centrality (13), betweenness centrality (14) and subgraph centrality (15). The hub genes were then subjected to cluster analysis using the R gplots software package.
Transcription factor (TF) analysis
TFs encoded by DEGs were explored combined with the mouse TF information recorded in The Animal Transcription Factor DataBase (AnimalTFDB) (http://www.bioguo.org/AnimalTFDB/index.php) (16). Based on the TRANSFAC database, the TF-DEG interactions were predicted according to the information provided in the TRANSFAC database using the cytoscape plugin termed iRegulon (17).
Results
Identification of DEGs
After preprocessing, the gene expression data were normalized (Fig. 1). According to values of p<0.05 and |log2FC(fold change)|≥0.58, a total of 210 genes were found to be differentially expressed in the diabetic group, including 92 upregulated and 118 downregulated genes. Hierarchical clustering analysis revealed that the upregulated and downregulated genes were relatively distinguished between the different groups (Fig. 2).
As shown in Fig. 3, the majority of upregulated genes were located on chromosome 3, whereas chromosomes 9 and 11 had the majority of downregulated genes.
GO and pathway analysis
The over-represented GO terms for the upregulated and downregulated genes are listed in Table I. Significantly enriched GO terms for upregulated genes included the regulation of angiogenesis, extracellular region and vascular endothelial growth factor receptor binding. For genes that were downregulated, the significantly enriched GO terms included cellular ion homeostasis, apical plasma membrane and symporter activity.
Of the 7 significantly enriched pathways, Pentose and glucuronate interconversions were found to be closely associated with the downregulated genes. Other pathways, such as the p53 signaling pathway, cytokine-cytokine receptor interaction, the transforming growth factor (TGF)-β signaling pathway and the mitogen-activated protein kinase (MAPK) signaling pathway were significantly enriched by upregulated genes (Table II).
PPI network
With a PPI score >0.4, a PPI network with 123 nodes and 219 edges was constructed, as shown in Fig. 4. The top 15 nodes based on the degree, betweenness and subgraph centralities were screened (Table III). Top 3 nodes such as early growth response 1 (Egr1), Serpine1 [also known as plasminogen activator inhibitor-1 (PAI-1)] and Mdk were shared based on the degree, betweenness and subgraph centralities. Merging the overlapping genes, we obtained 24 significant genes. Hierarchical clustering analysis revealed that the diabetic and non-diabetic samples were distinguished based on the gene expression profiles of the 24 significant genes (Fig. 5), suggesting that these genes were feature genes in diabetic samples.
In order to analyze the pathways associated with these feature genes, we performed KEGG pathway analysis. As shown in Table IV, the significant genes were closely associated with the MAPK signaling pathway, the p53 signaling pathway and the TGF-β signaling pathway.
Analysis of TFs
Combined with the TF information recorded in TFDB, we obtained 12 differentially expressed TFs from 9 TF families (Table V). The interactions between differentially expressed TFs and DEGs predicted by the TRANSFAC database are shown in Fig. 6. The TFs, nuclear receptor subfamily 4, group A, member 1 (Nr4a1) and peroxisome proliferator-activated receptor gamma (Pparg), were shown to have interactions with multiple genes.
Discussion
DN is a chronic kidney disease and is more prevalent in patients with diabetes mellitus (18). ND increases the risk factor of cardiovascular disease and mortality in diabetic patients (19). There is thus a need for the development of more effective treatments for patients with DN. In this study, we attempted to explore the potential molecular mechanisms of DN based on the bioinformatics methods and to provide a prospective novel therapeutic target. In the present study, we screened out 92 upregulated and 118 downregulated genes. All the genes with changes in expression were proven to be significant according to hierarchical clustering analysis.
KEGG pathway analysis for both DEGs and significant nodes in the PPI network revealed that the p53 signaling pathway, the TGF-β signaling pathway and the MAPK signaling pathway were the significantly enriched pathways. It has been reported that the overexpression of p53 is associated with the progression of DN. The expression of p53 and TGF-β was shown to be overexpressed in the renal cortex of diabetic mice. The crosstalk between p53 and miR-192, which is mediated by TGF-β was shown to be involved in the pathogenesis of DN (20). The MAPK signaling pathway is a regulator of the expression of pro-inflammatory molecules in DN. Targeted therapy for inhibiting the p38 MAPK signaling pathway has shown preventive effects on streptozotocin-induced DN (21). The p38 MAPK signaling pathway also plays a partial role in fibrosis associated with DN (22). These findings suggest that our findings are significant.
In our study, the PPI network showed that Egr1, Serpine1 and Mdk were the top 3 nodes based on the degree of centrality, betweenness centrality and subgraph centrality. Serpine1, also known as PAI-1, is a serine protease inhibitor and a key regulator of extracellular matrix (ECM). PAI-1 has been widely investigated in many diseases including the kidney disease (23,24). It is reported that ECM accumulation is implicated in the development and progression of DN (25). The expression of PAI-1 contributes to the fibrosis of kidney by inhibiting ECM degradation (26). PAI-1 has been proposed to be the potential target in renal fibrogenesis (26). PAI-1 was found to be overexpressed in the kidney of diabetic mice and its deficiency prevents glomerular injury of diabetic mice (27). Therefore, PAI-1 contributes to the progression of DN and PAI-1 knockdown may prove to be an effective therapeutic strategy for the treatment of DN. Additionally, midkine encoded by the Mdk gene has been proven to play a physiological role in kidney disease, including DN. Mdk plays a role in the occurrence and progression of acute kidney injury and contributes to the development of DN (28). Recent evidence has indicated that the diverse role of Mdk may open a new avenues for targeted therapies for DN (28).
Furthermore, Egr1 is a TF and plays a role in inducing the overexpression of heparanase in DN. The upregulation of heparanase is closely associated with albuminuria and renal damage in diabetic mice (29). The inhibition of Egr1 may be an effective strategy for preventing DN in diabetes. In the present study, other TFs, such as Nr4a1 and Pparg were found to play regulatory roles in the differential expression of genes. Nr4a1 is a member of the nuclear orphan receptor family of TFs. Nr4a1 has been found to play a significant role in atherosclerosis, psoriasis and other chronic inflammatory diseases (30). As we all know, hypertension is closely related with the progression of DN. Nr4a1 has been found to be differentially expressed in the kidneys of hypertensive patients (31). Nr4a1 is a susceptible gene in chronic kidney disease and a deficiency in Nr4a1 has been shown to be involved in kidney injury and renal function disorder (32). Otherwise, the other TF, Pparg, from the TF PPAR receptor family has been found to be associated with DN. Previous studies have demonstrated that Pparg variation contributes to DN development in type 2 diabetes (33). Thus, the significant genes identified in our study may be candidate therapeutic targets in DN.
Although we predicted the significant genes and pathways involved in the pathogenesis of DN, the lack of experimental validation was a limitation in our study. Experimental studies need be conducted to validate the differential gene expression profile and gene interaction pairs in the future.
In conclusion, we proved the significant role of the p53 signaling pathway, the TGF-β signaling pathway and the MAPK signaling pathway in the progression of DN. The significant genes, such as Egr1, Serpine1, Mdk, Nr4a1 and Pparg may prove to be potential therapeutic targets for the treatment of DN. Our findings need be further validated by experimental evidence.