Bioinformatics prediction and experimental verification of key biomarkers for diabetic kidney disease based on transcriptome sequencing in mice

Abstract

Background: Diabetic kidney disease (DKD) is the leading cause of death in people with type 2 diabetes mellitus (T2DM). The main objective of this study is to find the potential biomarkers for DKD.

Materials and methods: Two datasets (GSE86300 and GSE184836) retrieved from Gene Expression Omnibus (GEO) database were used, combined with our RNA sequencing (RNA-seq) results of DKD mice (C57 BLKS-32w db/db) and non-diabetic (db/m) mice for further analysis. After processing the expression matrix of the three sets of data using R software "Limma", differential expression analysis was performed. The significantly differentially expressed genes (DEGs) (-logFC- > 1, p-value < 0.05) were visualized by heatmaps and volcano plots respectively. Next, the co-expression genes expressed in the three groups of DEGs were obtained by constructing a Venn diagram. In addition, Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis were further analyzed the related functions and enrichment pathways of these co-expression genes. Then, qRT-PCR was used to verify the expression levels of co-expression genes in the kidney of DKD and control mice. Finally, protein-protein interaction network (PPI), GO, KEGG analysis and Pearson correlation test were performed on the experimentally validated genes, in order to clarify the possible mechanism of them in DKD.

Results: Our RNA-seq results identified a total of 125 DEGs, including 59 up-regulated and 66 down-regulated DEGs. At the same time, 183 up-regulated and 153 down-regulated DEGs were obtained in GEO database GSE86300, and 76 up-regulated and 117 down-regulated DEGs were obtained in GSE184836. Venn diagram showed that 13 co-expression DEGs among the three groups of DEGs. GO analysis showed that biological processes (BP) were mainly enriched inresponse to stilbenoid, response to fatty acid, response to nutrient, positive regulation of macrophage derived foam cell differentiation, triglyceride metabolic process. KEGG pathway analysis showed that the three major enriched pathways were cholesterol metabolism, drug metabolism-cytochrome P450, PPAR signaling pathway. After qRT-PCR validation, we obtained 11 genes that were significant differentially expressed in the kidney tissues of DKD mice compared with control mice. (The mRNA expression levels of Aacs, Cpe, Cd36, Slc22a7, Slc1a4, Lpl, Cyp7b1, Akr1c14 and Apoh were declined, whereas Abcc4 and Gsta2 were elevated).

Conclusion: Our study, based on RNA-seq results, GEO databases and qRT-PCR, identified 11 significant dysregulated DEGs, which play an important role in lipid metabolism and the PPAR signaling pathway, which provide novel targets for diagnosis and treatment of DKD.

Keywords: Bioinformatics analysis; Biomarker; Diabetic kidney disease (DKD); Differentially expressed genes (DEGs); RNA-seq.

Figure 1. The workflow of this study.

Figure 2. Biochemical parameters of blood and urine in DKD mice model of different week ages.

(A) Comparison of serum creatinine in mice; (B) comparison of blood urea nitrogen in mice; (C) comparison of total cholesterol in mice; (D) comparison of triglycerides in mice; (E) comparison of low-density lipoprotein in mice; (F) comparison of high-density lipoprotein in mice; (G) weight comparison of mice; (H) comparison of blood glucose in mice; (I) comparison of urinary albumin excretion rate (UAER) in mice. * p < 0.05, ** p < 0.01, *** p < 0.001.

Figure 3. Differential expression analysis of our RNA-seq and two GEO datasets (GSE86300 and GSE184836).

(A) Heatmap of DEGs in our RNA-seq. (B) Volcano map of our RNA-seq. A total of 59 upregulated and 66 downregulated DEGs were identified between the DKD and the Control group. (C) Heatmap of DEGs in GSE86300. (D) Volcano map of DEGs in GSE86300. A total of 183 up-regulated and 153 down-regulated DEGs were identified between the DKD and the Control group. (E) Heatmap of DEGs in GSE184836. (F) Volcano map of DEGs in GSE184836. A total of 76 up-regulated and 117 down-regulated DEGs were identified between the DKD and the control group. (G) Venn diagram of three DEGs groups. A total of 13 co-expression genes were obtained. Volcano map exhibit significantly differentially expressed genes, in volcano map, red bubbles mean up-regulated genes, blue bubbles mean down-regulated genes, and gray bubbles mean non-significant genes. The dots in the area above the horizontal dotted line have a P-value <0.05. The dots outside the two vertical dotted lines have a —log2FC— >1. Based on gene expression matrix, clustering analysis was shown in heatmap, in heatmap, red mean up-regulated genes, blue mean down-regulated genes. (—log2 FC— >1 and p-value<0.05).

Figure 4. Bar graph of GO Annotation and dot plot of KEGG pathway enrichment analysis of DEGs.

(A, B) Our RNA-seq (Top 10). (C, D) GSE86300 (Top 10). (E, F) GSE184836 (Top 10). Bar graph shows that DEGs of the three groups are enriched in several biological processes (BP), cell components (CC), molecular functions (MF). In the bar graph, we sorted the top 10 of BP, CC and MF by p-value and visualize them. In the dot plot, the color represents the p-value, and the size of the spots represents the gene number.

Figure 5. The relative mRNA expression of 13 overlapping genes in DKD and control group determined by qRT-PCR.

* p < 0.05, ** p < 0.01, *** p < 0.001.

Figure 6. Functional analysis of key genes and their interactions.

(A) Dot plot of KEGG pathway enrichment analysis of verified 11 genes, in the dot plot, the color represents the p-value, and the size of the spots represents the gene number. (B) Bar plot of GO enrichment analysis of verified 11 genes. (C) PPI network of verified 11 DEGs, which also showed that other genes were involved. In the PPI analysis, green represents the monocarboxylic acid metabolic process, red represents the PPAR signaling pathway and blue represents the lipid metabolic process, white has no meaning. (D) The protein networks of the three key genes and several related genes. (E) Pearson correlation test between 11 DKD related genes, blue represents negative correlation, while red represents positive correlation.