Prediction of the molecular mechanisms and potential therapeutic targets for diabetic nephropathy by bioinformatics methods.

In this study, we aimed to explore the molecular mechanisms of and genetic factors influencing diabetic nephropathy (DN). Gene expression profiles associated with DN were obtained from the GEO database (Accession no. GSE20844). The differentially expressed genes (DEGs) between diabetic mice and non-diabetic mice were screened. Subsequently, the DEGs were subjected to functional and pathway analysis. The protein-protein interaction (PPI) network was constructed and the transcription factors (TFs) were screened among the DEGs. A total of 92 upregulated and 118 downregulated genes were screened. Pathway analysis revealed that the p53 signaling pathway, the transforming growth factor (TGF)-β signaling pathway and the mitogen-activated protein kinase (MAPK) signaling pathway were significantly enriched by upregulated genes. Serpine1 (also known as plasminogen activator inhibitor-1), early growth response 1 (Egr1) and Mdk were found to be significant nodes in the PPI network by three methods. A total of 12 TFs were found to be differentially expressed, of which nuclear receptor subfamily 4, group A, member 1 (Nr4a1) and peroxisome proliferator-activated receptor gamma (Pparg) were found to have multiple interactions with other DEGs. We demonstrated that the p53 signaling pathway, the TGF-β signaling pathway and the MAPK signaling pathway were dysregulated in the diabetic mice. The significant nodes (Serpine1, Egr1 and Mdk) and differentially expressed TFs (Nr4a1 and Pparg) may provide a novel avenue for the targeted therapy of DN.

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).

Figure 1

Gene expression data before and after normalization. Horizontal axis represents the sample symbol and the vertical axis represents the gene expression values. The black line in the box plot represents the median value of gene expression.

Figure 2

The heatmap of gene expression profiles in diabetic and non-diabetic samples. Green, low gene expression value; red, high gene expression value; black, no differential expression.

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.

Figure 3

Distribution information of the upregulated and downregulated genes on the respective chromosomes.

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.

Table I

The significant GO terms enriched by DEGs.

Term	Count	P-value
Upregulated genes
BP
GO:0042127 - regulation of cell proliferation	13	2.30E-05
GO:0045765 - regulation of angiogenesis	5	9.16E-05
GO:0016525 - negative regulation of angiogenesis	4	2.04E-04
GO:0008285 - negative regulation of cell proliferation	7	0.001184688
GO:0009611 - response to wounding	8	0.002365745
CC
GO:0005576 - extracellular region	20	5.60E-05
GO:0044421 - extracellular region part	11	0.001793488
GO:0046658 - anchored to plasma membrane	3	0.002476471
GO:0005615 - extracellular space	8	0.006777968
GO:0044459 - plasma membrane part	14	0.022475418
MF
GO:0005172 - vascular endothelial growth factor receptor binding	3	2.41E-04
GO:0005539 - glycosaminoglycan binding	5	0.002459879
GO:0030247 - polysaccharide binding	5	0.00373434
GO:0001871 - pattern binding	5	0.00373434
GO:0030246 - carbohydrate binding	7	0.00480543
Downregulated genes
BP
GO:0006873 - cellular ion homeostasis	7	0.002754032
GO:0055082 - cellular chemical homeostasis	7	0.003139335
GO:0048878 - chemical homeostasis	8	0.003390517
GO:0050801 - ion homeostasis	7	0.004850765
GO:0019725 - cellular homeostasis	7	0.010204447
CC
GO:0016324 - apical plasma membrane	6	3.76E-04
GO:0005576 - extracellular region	22	0.001139987
GO:0045177 - apical part of cell	6	0.001580783
GO:0005615 - extracellular space	9	0.014078693
GO:0005903 - brush border	3	0.020010741
MF
GO:0015293 - symporter activity	5	0.008223849
GO:0019807 - aspartoacylase activity	2	0.012005144
GO:0008201 - heparin binding	4	0.013672168
GO:0004046 - aminoacylase activity	2	0.023867944
GO:0031402 - sodium ion binding	4	0.024166754

GO, Gene Ontology; DEGs, differentially expressed genes; BP, biological process; CC, cellular component; MF, molecular function.

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).

Table II

Pathways significantly enriched by differentially expressed genes.

Term	Count	P-value
Upregulation
mmu04115: p53 signaling pathway	5	7.57E-04
mmu04110: Cell cycle	5	0.007243672
mmu04060: Cytokine-cytokine receptor interaction	6	0.015248477
mmu04350: TGF-β signaling pathway	4	0.015533727
mmu04010: MAPK signaling pathway	6	0.021089527
mmu05219: Bladder cancer	3	0.026718766
Downregulation
mmu00040: Pentose and glucuronate interconversions	3	0.003652816

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.

Figure 4

Protein-protein interaction (PPI) network for differentially expressed genes (DEGs). Green, downregulated genes; red, upregulated genes.

Table III

Top 15 significant nodes in the protein-protein interaction network based on degree, betweenness and subgraph centralities.

Node	Degree	Node betweenness	Betweenness	Node subgragh	Subgragh
Serpine1	16	Serpine1	2905.323	Egr1	337.665
Egr1	15	Egr1	2433.539	Serpine1	337.1419
Mdk	14	Mdk	2319.191	Mdk	285.9198
Nr4a1	12	Pparg	2105.681	Icam1	209.8053
Icam1	12	Lrp2	2064.846	Cdkn1a	187.7399
Cdkn1a	12	Cdh16	1583.724	Nr4a1	178.7159
Tlr4	11	Umod	1440.558	Tlr4	147.8937
Cdk1	10	Ltf	1421.087	Dusp1	119.8608
Bmp2	10	Ggt1	1384.699	Bmp2	117.067
Dusp1	9	Tlr4	1286.763	Cyr61	104.8729
Pparg	8	Slc12a1	1250.098	Pparg	85.67106
Ltf	7	Nr4a1	1248.253	Cybb	68.2394
Egf	7	Icam1	1223.037	Smad6	63.11027
Cybb	7	Cdk1	1119.84	Cdk1	56.22912
Tgfb3	7	Wfdc15b	1082.356	Map3k8	52.32016

Figure 5

Hierarchical clustering analysis for the gene expression profiles of 24 significant genes in the protein-protein interaction (PPI) network.

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.

Table IV

The significant pathways enriched by 24 significant nodes in the protein-protein interaction network.

Term	P-value	Genes
mmu04010: MAPK signaling pathway	0.005210429	Dusp1, Map3K8, Tgfb3, Nr4A1, Egf
mmu05200: Pathways in cancer	0.010457331	Cdkn1A, Bmp2, Pparg, Tgfb3, Egf
mmu04115: p53 signaling pathway	0.015340116	Cdk1, Cdkn1A, Serpine1
mmu04350: TGF-β signaling pathway	0.023759324	Bmp2, Smad6, Tgfb3
mmu04110: Cell cycle	0.048314969	Cdk1, Cdkn1A, Tgfb3

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.

Table V

Differentially expressed transcription factors in diabetic samples.

TF family	TF gene symbol
bHLH	Atoh8
bHLH	Hey1
HMG	Tox2
Homeobox	Meox1
Homeobox	Six4
IRF	Irf6
MH1	Smad6
Nuclear orphan receptor	Nr4a1
PPAR receptor	Pparg
Thyroid hormone receptor	Vdr
zf-C2H2	Klf6
zf-C2H2	Egr1

Figure 6

Interactions between transcription factors (TFs) and target genes. Red, upregulated genes; green, downregulated genes. Triangle represents TFs and the box plot represents target 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.