DEGs identifified between DKD and Normal patients
The expression matrices of the 2 datasets GSE30528 and GSE30529 were normalized. The boxplots of the intensity of all samples demonstrated that the expression values of each sample were close to the same after normalization (Fig. 1). Unbiased clustering was performed by principle component analysis (PCA), followed by PCA plots and UMAP dimensionality reduction visualization. Each dot in the scatterplot represents a single sample separately(Fig. 2A-B). PCA analysis results showed that PC1 was 21% and PC2 was 16.9% and showed a clear difference. After screening with the threshold of an adjusted |log2(FC)| >1 and P value < 0.05, a total of 274 DEGs between DKD and Control patients were obtained. The 205 up-regulated genes and the 69 down-regulated genes were shown in the volcano plot (Fig. 3A). Blue dots represent the downregulated DEGs; red dots represent the upregulated DEGs. A heatmap of differentially expressed genes shows the expression levels of the top20 with upregulated and downregulated genes (red, high expression; blue, low expression) (Fig. 3B). Venn plots were also created (Fig. 4), Using the “VennDiagram”R package 45, we conducted Venn diagram analysis of DEGs and genes screened from the CTD Database associated with DKD. These diagrams show overlapping genes. Figure 4A shows the intersection of upregulated and positively relevant genes in DKD (see the attached table). Figure 4B shows the intersection of downregulated and negatively relevant genes in DKD.
GO and KEGG pathway enrichment analyses of up-regulated DEGs
To reveal deeper insights into the biological functions of DEGs, GO annotation and KEGG pathway enrichment analyses were performed using DAVID(Figure 5). A total of 772 GO terms and 44 KEGG pathways were enriched and set the criteria with P-value < 0.05 as signifificantly enriched. The GO terms were comprised of 3 parts: 662 biological process (BP), 54 cellular component (CC), and 56 molecular function (MF).The top 15 enriched GO terms and The top 10 enriched KEGG pathway were shown in Circle map(Figs. 6A) and Chord plot(Figs. 6B). The result of GO enrichment indicated that for BP, up-regulated DEGs were significantly enriched in positive regulation of cell activation, positive regulation of leukocyte activation, regulation of peptidase activity,leukocyte cell-cell adhesion, regulation of tumor necrosis factor production etc. Regarding CC, DEGs were significantly enriched in collagen-containing extracellular matrix, cytoplasmic vesicle lumen, vesicle lume,secretory granule lumen, platelet alpha granule etc. For MF, up-regulated genes were significantly enriched in glycosaminoglycan binding heparin binding, sulfur compound binding, peptidase regulator activity, extracellular matrix structural constituent etc.(Figs. 6A and Table 1 ) These DEGs were also enriched in KEGG pathways including Phagosome Hematopoietic cell lineage, Staphylococcus aureus infection, Rheumatoid arthritis, Leishmaniasis, Cell adhesion molecules, Focal adhesion etc.(Figs. 6B and Table 2).
Ontology
|
ID
|
Description
|
GeneRatio
|
pvalue
|
p.adjust
|
Table 1
GO function analysis for up-regulated DEGs. The top 5 terms were selected based upon p value rankings.
BP
|
GO:0050867
|
positive regulation of cell activation
|
28/247
|
4.56e-13
|
1.69e-09
|
BP
|
GO:0002696
|
positive regulation of leukocyte activation
|
27/247
|
1.23e-12
|
2.15e-09
|
BP
|
GO:0052547
|
regulation of peptidase activity
|
29/247
|
2.19e-12
|
2.15e-09
|
BP
|
GO:0050900
|
leukocyte migration
|
30/247
|
4.61e-12
|
2.15e-09
|
BP
|
GO:0045785
|
positive regulation of cell adhesion
|
27/247
|
4.87e-12
|
2.15e-09
|
CC
|
GO:0062023
|
collagen-containing extracellular matrix
|
29/251
|
5.86e-14
|
1.77e-11
|
CC
|
GO:0060205
|
cytoplasmic vesicle lumen
|
24/251
|
1.10e-11
|
1.19e-09
|
CC
|
GO:0031983
|
vesicle lumen
|
24/251
|
1.17e-11
|
1.19e-09
|
CC
|
GO:0034774
|
secretory granule lumen
|
23/251
|
2.52e-11
|
1.91e-09
|
CC
|
GO:0031091
|
platelet alpha granule
|
12/251
|
1.81e-09
|
1.09e-07
|
MF
|
GO:0005539
|
glycosaminoglycan binding
|
25/242
|
1.01e-15
|
4.74e-13
|
MF
|
GO:0008201
|
heparin binding
|
20/242
|
1.80e-13
|
4.21e-11
|
MF
|
GO:1901681
|
sulfur compound binding
|
21/242
|
3.52e-11
|
5.48e-09
|
MF
|
GO:0061134
|
peptidase regulator activity
|
19/242
|
1.80e-10
|
2.11e-08
|
MF
|
GO:0005201
|
extracellular matrix structural constituent
|
15/242
|
7.05e-09
|
6.59e-07
|
Ontology
|
ID
|
Description
|
GeneRatio
|
pvalue
|
p.adjust
|
Table 2
KEGG pathway enrichment analyses for up-regulated DEGs. The top 10 terms were selected based upon p value rankings.
KEGG
|
hsa04145
|
Phagosome
|
17/154
|
3.62e-09
|
4.92e-07
|
KEGG
|
hsa04640
|
Hematopoietic cell lineage
|
14/154
|
4.37e-09
|
4.92e-07
|
KEGG
|
hsa05150
|
Staphylococcus aureus infection
|
13/154
|
2.74e-08
|
2.05e-06
|
KEGG
|
hsa05323
|
Rheumatoid arthritis
|
12/154
|
1.66e-07
|
7.20e-06
|
KEGG
|
hsa05133
|
Pertussis
|
11/154
|
1.67e-07
|
7.20e-06
|
KEGG
|
hsa05140
|
Leishmaniasis
|
11/154
|
1.92e-07
|
7.20e-06
|
KEGG
|
hsa04514
|
Cell adhesion molecules
|
14/154
|
8.23e-07
|
2.64e-05
|
KEGG
|
hsa04510
|
Focal adhesion
|
16/154
|
1.24e-06
|
3.49e-05
|
KEGG
|
hsa05152
|
Tuberculosis
|
15/154
|
1.53e-06
|
3.82e-05
|
KEGG
|
hsa05144
|
Malaria
|
8/154
|
3.97e-06
|
8.91e-05
|
PPI Network Construction and Hub Gene Identifification
To screen the differential genes that were strongly linked to DKD, we performed protein interaction network analysis on the basis of 194 DEGs .Nondirectional scores were calculated for individual nodes in the protein interaction network using the Network Analyser in Cytoscape. PPI network was constructed involving 168 nodes and 1317 edges. The nodes indicate the DEGs and the edges indicate the interactions between two genes.The color depth of node from pink to rose red refers to the degree standard and represents the values; the color depth of edges from white to blue refers to the combined-score standard and represents the importance in the network.The connectivity degree is an important parameter and the high connectivity degree indicated the protein interacted with more surrounding proteins and play a more important role. The nodes and edges with higher connectivity degree were shown as deeper color with rose red or blue color. All the protein nodes were arranged using the layout of the Degree Sorted Circle Layout (Fig. 7). According to the degree value, the top 15 hub DEGs were determined using cytoHubba. These hub genes were PTPRC、ITGAM、ITGB2、TYROBP、FN1、CD44、SELL、CYBB、IL10RA、LCP2、CCL5、IRF8、CCR2、CCL2、CD53(Fig. 8A). Hence, we consider the 15 genes as the hub genes for further analysis. These hub DEGs were also enriched in KEGG pathways including Cell adhesion molecules, AGE-RAGE signaling pathway in diabetic complications, Cytokine-cytokine receptor interaction, Natural killer cell mediated cytotoxicity, NOD-like receptor signaling pathway, Chemokine signaling pathway, Rap1 signaling pathway, ECM-receptor interaction, T cell receptor signaling pathway, TNF signaling pathway etc. (Figs. 8B and Table 3).
Ontology
|
ID
|
Description
|
GeneRatio
|
pvalue
|
p.adjust
|
qvalue
|
Table 3
KEGG
|
hsa04514
|
Cell adhesion molecules
|
4/14
|
9.65e-05
|
0.003
|
0.002
|
KEGG
|
hsa04933
|
AGE-RAGE signaling pathway in diabetic complications
|
3/14
|
6.07e-04
|
0.006
|
0.003
|
KEGG
|
hsa04060
|
Cytokine-cytokine receptor interaction
|
4/14
|
0.001
|
0.008
|
0.004
|
KEGG
|
hsa04650
|
Natural killer cell mediated cytotoxicity
|
3/14
|
0.001
|
0.008
|
0.004
|
KEGG
|
hsa04621
|
NOD-like receptor signaling pathway
|
3/14
|
0.003
|
0.014
|
0.008
|
KEGG
|
hsa04062
|
Chemokine signaling pathway
|
3/14
|
0.004
|
0.015
|
0.008
|
KEGG
|
hsa04015
|
Rap1 signaling pathway
|
3/14
|
0.005
|
0.018
|
0.010
|
KEGG
|
hsa04512
|
ECM-receptor interaction
|
2/14
|
0.010
|
0.028
|
0.015
|
KEGG
|
hsa04660
|
T cell receptor signaling pathway
|
2/14
|
0.014
|
0.032
|
0.018
|
KEGG
|
hsa04668
|
TNF signaling pathway
|
2/14
|
0.016
|
0.036
|
0.019
|
Core genetic screening associated to diabetic nephropathy
Using the “VennDiagram”R package, we conducted Venn diagram analysis of 205 upregulated key DEGs and A search of 144 confirmed diabetic nephropathy-related genes in the CTD database. As shown in the Venn diagram (Fig. 9A), the five candidate genes reported in the literature were: FN1, TGFBI, CCL2, VDAC1, and TXNIP, in which FN1 and CCL2 are simultaneously overlapped with Top 15 hub genes (Fig. 9B).
Construction of an mRNA-miRNA regulatory network
To further assess the potential of miRNAs as a DKD marker to screen out important miRNA and mRNA, We used the miRWalk database to predict target miRNAs of 18 target genes (TGFBI, VDAC1, TXNIP, PTPRC, ITGAM, ITGB2, TYROBP, FN1, CD44, SELL, CYBB, IL10RA, LCP2, CCL5, IRF8, CCR2, CCL2, CD53). At the same time, combined with the miRTarBase database, we obtained 151 target miRNAs of 5 specifically expressed target genes (TXNIP, CD44, CYBB, IL10RA, VDAC1) and determined 264 mRNA-miRNA pairs. According to the prediction results, a co-expressed network of mRNAs and miRNAs, which comprised 105 nodes and 264 edges, was constructed by Cytoscape (Fig. 10).
GSEA Enrichment Analysis of Biological Functions and Pathways of Hub Genes
To further investigate the potential functions of TXNIP, CD44, CYBB, FN1, CCL2, TGFBI and VDAC1 in DKD, we performed Gene Set Enrichment Analysis (GSEA) (Fig. 11 ). From the GSEA results, the DKD-related of FN1 and TGFBI were all enriched in “core matrisome ”, “ECM glycoproteins ”, “integrin3” and “integrin1” pathways. In the KEGG enrichment analysis, FN1 and CD44 were particularly involved in “ECM receptor interaction” and “hematopoietic cell lineage”. Interestingly, this analysis revealed that the 2 pathways were also enriched in the same direction at the GSEA.“signaling by interleukins”, “IL18 signaling pathway ” and “interleukin4 and interleukin13 signaling pathway” were enriched in the FN1 and CCL2 high-expression groups. Genes in high expression groups of TXNIP和CCL2 were significantly enriched with the NOD-like receptor signaling pathway, and CCL2 was also enriched in AGE-RAGE signaling pathway in diabetic complications. In the GSEA enrichment analysis, after screening with the threshold of |NES|>1, False discovery rate (FDR) < 0.25 and p.adjust < 0.05, in high expression groups of TXNIP was enriched in REACTOME _PURINERGIC _SIGNALING _IN _LEISHMANIASIS _INFECTION (NES = 2.071, P = 0.001, FDR = 0.019);REACTOME _THE _NLRP3 _INFLAMMASOME (NES = 1.820, P = 0.046, FDR = 0.0367)༛REACTOME _INFLAMMASOMES (NES = 1.823, P = 0.036, FDR = 0.028), REACTOME _NUCLEOTIDE _BINDING _DOMAIN _LEUCINE _RICH _REPEAT _CONTAINING _RECEPTOR _NLR _SIGNALING _PATHWAYS (NES = 1.899, P = 0.025, FDR = 0.019), REACTOME _LEISHMANIA _INFECTION(NES = 1.715, P = 0.025, FDR = 0.019), (Fig. 12 ).
Expression and Screening Identification of Target Genes Associated With Diabetic Nephropathy
GSE30528 and GSE30529 were used to detect the expression of the screened target genes. The results showed that the expression of five of the seven diabetic nephropathy-associated pivotal genes (TXNIP, CD44, CYBB, FN1, CCL2, TGFBI and VDAC1) between healthy individuals (TXNIP, CD44 ,FN1, CCL2, TGFBI) was consistent with the predicted results, with all DKD groups higher than the Control group, and the differences were statistically significant (Fig. 13). To further assess the diagnostic value of target genes in diabetic nephropathy, ROC curves were established with area under the ROC curve values between 0.5 and 1. The closer the AUC was to 1, the better the diagnosis was. In the GSE30528 and GSE30529 datasets, the variables TXNIP, CD44 ,FN1, and TGFBI all had high accuracy in predicting the predictive ability of outcome in normal patients and diabetic nephropathy patients, and CCL2 had average accuracy (Fig. 14A-J). We also used the GSE1009 database for validation, which included three diabetic nephropathy patients and three normal patients, to create validation ROC curves. The results showed that GSE1009 confirmed that TXNIP, CD44 ,FN1, TGFBI, and CCL2 were also of diagnostic value (Fig. 14K-O), which was consistent with the predicted results in GSE30528 and GSE30529. Therefore, we used these five genes (FN1, TGFBI, CCL2, CD44, TXNIP) as pivotal genes for further analysis.
Changes of serum indexes, TXNIP, CD44, TGF-β1 and MicroRNA-130a in DN and healthy participants
As shown, the levels of blood glucose, urinary microalbumin, low-density lipoprotein (LDL-C), triglyceride (TG), total cholesterol (TC), urea nitrogen (BUN), and creatinine (Cr) were significantly elevated in DN patients compared to healthy controls (all P < 0.05) (Table 4 − 1). High-density lipoprotein (HDL-C) was significantly decresaed in DN patients compared to healthy controls ( P < 0.05) (Table 4 − 1). Moreover, serum TXNIP, CD44, and TGF-β1 in DN patients were markedly higher, while serum MicroRNA-130a were markedly lower than healthy controls (P < 0.05) (Table 4 − 2). Serum TXNIP, CD44, and TGF-β1 was negative correlated with MicroRNA-130a in DN patients (Fig. 15).
Table 4
− 1. Blood and urinary parameters of DN patients
Groups
|
n
|
blood glucose (mmol/L)
|
urinary microalbumin (24h, mg)
|
HDL-C
(mmol/L)
|
LDL-C (mmol/L)
|
TG (mmol/L)
|
TC
(mmol/L)
|
BUN
(mmol/L)
|
Cr
(mmol/L)
|
NC
|
50
|
4.76 ± 0.70
|
24.08 ± 4.39
|
1.60 ± 0.37
|
2.04 ± 0.57
|
0.99 ± 0.37
|
3.50 ± 1.23
|
3.02 ± 0.84
|
65.85 ± 14.92
|
DN
|
100
|
11.14 ± 3.30
|
127.97 ± 41.99
|
0.63 ± 0.29
|
3.75 ± 0.78
|
3.29 ± 0.99
|
7.35 ± 1.03
|
10.97 ± 2.79
|
131.06 ± 11.48
|
t
|
|
-13.49
|
-11.65
|
17.37
|
-13.73
|
-15.80
|
-20.19
|
-19.64
|
-29.52
|
P
|
|
< 0.01
|
< 0.01
|
< 0.01
|
< 0.01
|
< 0.01
|
< 0.01
|
< 0.01
|
< 0.01
|
Table 4
− 2. Serum parameters of DN patients
Groups
|
n
|
MicroRNA-130a
|
TXNIP(pg/mL)
|
CD44 (pg/mL)
|
TGF-β1(ng/mL)
|
NC
|
50
|
1.61 ± 0.94
|
81.83 ± 13.10
|
40.52 ± 13.27
|
11.55 ± 5.86
|
DN
|
100
|
0.62 ± 0.28
|
196.42 ± 55.73
|
114.78 ± 42.72
|
48.76 ± 23.99
|
t
|
|
9.70
|
-14.32
|
-11.99
|
-10.79
|
P
|
|
< 0.01
|
< 0.01
|
< 0.01
|
< 0.01
|
Expression levels of MicroRNA-130a, TXNIP, CD44, TGF-β1 and pathological changes of liver and heart in DN rats.
The MicroRNA-130a levels in the liver and heart were significantly decreased in diabetic nephropathy group compared to the normal control group (P < 0.05, Figure.16). The expression levels of TXNIP, CD44 and TGF-β1in the liver were significantly elevated in diabetic nephropathy group compared to the normal control group (P < 0.05, Fig. 17A-B and 18A). The levels of that in the heart were significantly higer in the diabetic nephropathy group compared to the normal control group (P < 0.05, Fig. 17C-D and 18B). Immunohistochemistry analysis showed that the protein levels in hepatocytes of TXNIP, CD44, TGF-β1 in the diabetic nephropathy group were markedly higer than those in the normal control group(Fig. 19). HE staining revealed unclear hepatic tissue structure, swollen liver cells, loose cytoplasm, fatty degeneration, and focal necrosis with inflammatory cell infiltration in tissue sections from DN rats. Hypertrophic, distorted, and disorganized myocardial cells, increased extracellular matrix, fibroblasts, and inflammatory cells were also evident(Fig. 20).