Microarray analysis reveals gene and microRNA signatures in diabetic kidney disease.

The current study aimed to identify therapeutic gene and microRNA (miRNA) biomarkers for diabetic kidney disease (DKD). The public expression profile GSE30122 was used. Following data preprocessing, the limma package was used to select differentially-expressed genes (DEGs) in DKD glomeruli samples and tubuli samples and they were compared with corresponding controls. Then overlapping DEGs in glomeruli and tubuli were identified and enriched analysis was performed. In addition, protein‑protein interaction (PPI) network analysis as well as sub‑network analysis was conducted. miRNAs of the overlapping DEGs were investigated using WebGestal. A total of 139 upregulated and 28 downregulated overlapping DEGs were selected, which were primarily associated with pathways involved in extracellular matrix (ECM)‑receptor interactions and cytokine‑cytokine receptor interactions. CD44, fibronectin 1, C‑C motif chemokine ligand 5 and C‑X‑C motif chemokine receptor 4 were four primary nodes in the PPI network. miRNA (miR)‑17‑5p, miR‑20a and miR‑106a were important and nuclear receptor subfamily 4 group A member 3 (NR4A3), protein tyrosine phosphatase, receptor type O (PTPRO) and Kruppel like factor 9 (KLF9) were all predicted as target genes of the three miRNAs in the integrated miRNA‑target network. Several genes were identified in DKD, which may be involved in pathways such as ECM‑receptor interaction and cytokine‑cytokine receptor interaction. Three miRNAs may also be used as biomarkers for therapy of DKD, including miR‑17‑5p, miR‑20a and miR‑106a, with the predicted targets of NR4A3, PTPRO and KLF9.

Introduction

Diabetic nephropathy (DN) is a complication of diabetes. It is estimated that ~33% of patients with diabetes develop end-stage kidney disease (ESKD), which may further progress to DN (1–3). Diabetic kidney disease (DKD) is the primary cause of kidney failure. DKD is frequently identified as a degenerative disease of the glomerulus and occurs due to alterations of the intrarenal metabolism and structure (4,5). Recently, the incidence rate of DKD has increased (6). In addition, type 2 diabetes (T2D)-ESKD has an unfavorable 5-year survival rate at <30% (7).

Pathogenesis of DKD is multifactorial and complicated. Inflammation and oxidative stress are identified as the two critical regulators in DKD. Nuclear factor, erythroid 2 like 2 (Nrf2) is a transcription factor that mediates the expression levels of antioxidant genes. Bardoxolone methyl, which may activate the expression of Nrf2, was determined to increase the estimated glomerular filtration rate in patients with DKD (8). Previous studies also suggested inflammation as a major pathophysiological mechanism for DKD progression (9,10).

Although clinical treatments for DKD have not been progressed, several effective agents and biological molecules have been identified. For example, metformin extracted from herbal medicines, may be a useful agent for the control of DKD, as it has hypoglycemic characteristics and renoprotective activities (11). Other herbal medicines have also been identified, such as curcumin and glycosides (12). A recent study revealed that the JAK-STAT pathway was associated with the development of DKD and gene expression of JAK2 was increased in glomerular podocytes and JAK2 inhibitor had an evident reduced effect on albuminuria, and may be used as a novel therapy of DKD (13).

MicroRNAs (miRNAs) have significant roles in the regulation of gene expression and cellular activities. A previous study determined that miR-192 was elevated in DKD via the mediation of renal fibrosis (14). Additional important miRNAs for DKD progression have also been previously identified, including miR-21, miR-29c, miR-200b, miR-216a and miR-377 (15).

However, comprehensive understanding of the pathogenesis of DKD is still unclear. A previous study established a microarray profile, GSE30122, identifying the crucial genes in DKD samples and highlighted several key molecules that were differentially expressed in glomeruli and the diabetic tubili, including cell division cycle 42 and vascular endothelial growth factor. Inflammation-associated pathways were primarily enriched for these genes (16). Although they performed expression validations of these genes, potential interactions between the differentially-expressed genes (DEGs), particularly at the protein level, were not investigated further. Therefore, Tang et al (5) used this expression profile to perform the weighted gene co-expression network analysis and identified 10 novel potential therapeutic targets for DKD, including ETS proto-oncogene 1, transcription factor, lipopolysaccharide induced TNF factor, nuclear factor, erythroid derived 2, retinoic acid receptor, γ and signal transducer and activator of transcription 5A (5). However, regulation from miRNAs was not taken into consideration. Therefore, the present study re-analyzed the dataset of GSE30122, and performed a series of bioinformatics analyses. Following the DEG selection in glomeruli and tubili, the overlapping DEGs were selected to perform enrichment analysis and protein-protein interaction (PPI) network analysis. Additionally, a sub-network was further extracted from the PPI network. It is of note that the present study identified the miRNAs that may target these DEGs and constructed an integrated miRNA-target regulatory network. Using these comprehensive analyses, the present study aimed to provide novel insight into the pathogenesis of DKD and find novel therapeutic biomarkers, particularly miRNAs.

Materials and methods

Data resources

The expression profile GSE30122 (16) was downloaded in the Gene Expression Omnibus (GEO, www.ncbi.nlm.nih.gov/geo) database, was used in the current study. The dataset consisted of 35 glomeruli and 34 tubuli tissue samples. The glomeruli samples were from 9 diabetic human kidney patients (DKD-G group) and 26 healthy controls (G-control group) and the tubuli samples were from 10 DKD patients (DKD-T group) and 24 healthy controls (T-control group). The platform used was Affymetrix Human Genome U133A 2.0 array (GPL571 HG-U133A_2; Affymetrix; Thermo Fisher Scientific, Inc., Waltham, MA, USA).

Data pretreatment and differential analysis

Raw data in the format of CEL was obtained and the Affy package (version 1.48.0) (17) in R (version 3.3.1) was used to perform background correction and quantiles normalization. After these pretreatments, differentially-expressed genes (DEGs) between DKD-G and G-control groups, and between DKD-T and T-control groups were identified, respectively, using non-paired t-test in the limma (www.bioconductor.org/packages/release/bioc/html/limma.html) package (18). The thresholds for significant DEG selection were log |fold-change|≥1 and P<0.05. In addition, the pheatmap (cran.r-project.org/web/packages/pheatmap/index.html) software was utilized to draw the clustering heatmap of these DEGs (19). Then the overlapped DEGs in the two comparisons were selected.

Enrichment analysis of the overlapped DEGs

The function and pathway enrichment analyses of the overlapped DEGs were performed using the online tool of Database for Annotation, Visualization and Integration Discovery (DAVID; david.abcc.ncifcrf.gov), in combination with the gene ontology (GO; www.geneontology.org) and the Kyoto Encyclopedia of Genes and Genomes (KEGG, http://www.genome.jp/kegg/pathway.html) databases. Statistical significance was determined based on hypergeometric test and the cut-off values for significant GO terms and pathway terms were gene count (number of enriched genes in a specific function or a pathway) was ≥2 and P<0.05.

Construction of the PPI network and sub-network

To explore potential interactions of the overlapped DEGs at the protein level, a PPI network was constructed based on the Search Tool for the Retrieval of Interacting Genes (STRING; string-db.org) database (20) and visualized using Cytoscape (version 3.2.0, http://cytoscape.org/) (21). Parameters for the network establishment were: The species was Homo sapiens and the PPI score (indicating medium confidence) was ≥0.4. A node was defined as the protein product of a DEG in the network and it was required that all nodes in the network were DEGs.

Score of a node in the network was calculated using one of the topological properties, degree centrality, which is based on the concept that nodes with the highest number of interactions were significant in the network (22). Hub nodes were those with the high degrees in the network.

Sub-networks were extracted from the PPI network to further elucidate the most significant functional modules of the DEGs, using the MCODE Cytoscape plug-in (23).

Integrated miRNA-target regulatory network construction

Potential miRNAs that may target the overlapping DEGs were identified using the WebGestal tool (version 2.0) (24). This was also visualized using Cytoscape software.

Results

DEGs identified in different DKD samples

Following preprocessing of raw data, the DEGs were identified when the DKD-G and G-control samples and DKD-T vs. T-control samples were compared based on the aforementioned criteria. A total of 680 or 724 DEGs were selected in DKD-G or DKD-T samples, respectively. A heatmap of gene expressions were presented in Fig. 1, which revealed that these DEGs may distinguish the DKD-G/DKD-T samples from the control samples.

Figure 1.

Heatmap of expression of differentially-expressed genes. (A) Glomeruli samples of patient with DKD vs. healthy control. (B) Tubuli samples of DKD patient vs. healthy control. DKD, diabetic kidney disease.

Using overlapping analysis, a total of 139 upregulated and 28 downregulated DEGs in the two sets of DEGs were identified.

Enrichment results of the overlapping DEGs

The upregulated 139 DEGs were significantly enriched in biological processes (BP), such as inflammatory/immune response, molecular functions (MF, such as peptidase inhibitor activity and 13 pathways, including extracellular matrix (ECM)-receptor interactions, including CD44, FN1 and collagen type VI α3 chain (COL6A3), primary immunodeficiency, including protein tyrosine phosphatase, receptor type, C (PTPRC), CD3d molecule and LCK (LCK proto-oncogene, Src family tyrosine kinase) and cytokine-cytokine receptor interaction, including C-X-C motif chemokine receptor 4 (CXCR4), C-C motif chemokine ligand 5 (CCL5) and chemokine (C-C motif) receptor 2 (Table I), whereas the 28 downregulated genes were primarily enriched in metabolic BP and MF, such as carbohydrate and polysaccharide binding. No significant pathways were enriched for the downregulated DEGs.

Table I.

Significant pathways for the 139 overlapping upregulated differentially-expressed genes.

ID	Name	P-value	Genes
hsa05322	Systemic lupus erythematosus	0.000136	C1QA, HLA-DQB1, C1QB, C7, C3, ACTN1, HLA-DPA1, HLA-DRA
hsa04512	ECM-receptor interaction	0.001105	CD44, TNC, COL6A3, COL1A2, THBS2, SPP1, FN1
hsa05340	Primary immunodeficiency	0.00143	PTPRC, CD3D, LCK, IL7R, BLNK
hsa04514	Cell adhesion molecules (CAMs)	0.001654	HLA-DQB1, PTPRC, SELL, CD2, VCAN, ITGB2, HLA-DPA1, HLA-DRA
hsa04510	Focal adhesion	0.006241	RAC2, TNC, COL6A3, COL1A2, PDGFRA, ACTN1, THBS2, SPP1, FN1
hsa05310	Asthma	0.006949	FCER1A, HLA-DQB1, HLA-DPA1, HLA-DRA
hsa04640	Hematopoietic cell lineage	0.007496	CD3D, CD44, CD1C, CD2, IL7R, HLA-DRA
hsa05020	Prion diseases	0.013215	C1QA, C1QB, C7, CCL5
hsa05416	Viral myocarditis	0.015647	HLA-DQB1, RAC2, ITGB2, HLA-DPA1, HLA-DRA
hsa04610	Complement and coagulation cascades	0.01646	C1QA, C1QB, C7, C3, CFB
hsa04060	Cytokine-cytokine receptor interaction	0.021393	CXCR4, IL10RA, CCR2, PDGFRA, CCL19, CXCL6, IL7R, CCL5, LTB
hsa04672	Intestinal immune network for IgA production	0.026364	HLA-DQB1, CXCR4, HLA-DPA1, HLA-DRA
hsa04062	Chemokine signaling pathway	0.036074	DOCK2, RAC2, CXCR4, CCR2, CCL19, CXCL6, CCL5

The PPI network of the overlapping DEGs

Using the aforementioned parameters, a PPI network was established, containing 115 nodes and 303 interactions. Top ten nodes with high degrees were ALB (degree = 47), PTPRC (degree = 29), LCK (degree = 18), CD44 (degree = 17), C3 (degree = 16), FN1 (degree = 16), APOE (degree = 15), CLU (degree = 13), CCL5 (degree = 11) and CXCR4 (degree = 11) (Fig. 2).

Figure 2.

Protein-protein interaction network of the overlapped differentially-expressed genes. Red circles, upregulated genes; green circles, downregulated genes.

Sub-network analysis obtained two significant module networks, A and B. Module A contained 13 nodes and 31 interactions and Module B consisted of 8 nodes and 28 interactions. It is of note that the majority of the of the top ten nodes identified in the PPI network were also present in these two modules (Fig. 3).

Figure 3.

Sub-network of the overlapped differentially-expressed genes. (A) Module A and (B) module B. Red circles, upregulated genes; green circles, downregulated genes.

Integrated miRNA-target regulatory network

The miRNAs that may regulate the DEGs were identified based on the up and downregulation expressions. As indicated in this integrated miRNA-target network (Fig. 4), only one miRNA, miR-373 was predicted to target the upregulated DEGs, such as ADAM like decysin 1, SRY-box 4, FN1, B-cell linker and SRY-box 9. Conversely, the network consisted of 11 downregulated DEGs, which were targeted by 37 miRNAs. Among them, OPCML, nuclear receptor subfamily 4 group A member 3 (NR4A3), protein tyrosine phosphatase, receptor type O (PTPRO) and Kruppel like factor 9 (KLF9) were identified as the four predominant genes with high topological scores, which suggested that they may be the targets of multiple miRNAs. It is of note that miR-17-5p, miR-20a and miR-106a were all predicted to simultaneously target the genes of NR4A3, PTPRO and KLF9.

Figure 4.

Integrated miRNA-target regulatory network. Red circles, upregulated genes; green circles, downregulated genes; triangles, miRNAs; arrows, regulatory relationship. miRNA, micro RNA.

Discussion

The present study re-analyzed the dataset of GSE30122 and identified several important genes that were differentially expressed in both of glomeruli and tubuli DKD samples, including CD44, FN1, CCL5 and CXCR4. All were predominant nodes in the PPI network and were associated with pathways like ECM-receptor interaction and cytokine-cytokine receptor interaction. It is of note that the present study identified important miRNAs, including miR-17-5p, miR-20a and miR-106a and the overlapping DEGs such as NR4A3, PTPRO and KLF9 were all predicted as target genes for the three miRNAs.

A hallmark of chronic kidney disease (CDK) is the accumulation of ECM in glomeruli and tubular interstitium, which may lead to glomerulosclerosis and tubulointerstitial fibrosis. Previous studies have demonstrated that transforming growth factor-β (TGF-β) is a key regulator in renal fibrosis (25–27) and its action is dependent on the activation of a protein complexes in ECM. Considering the close relationship of DKD and CDK, it is possible that the ECM-receptor interaction pathway is a significantly enriched pathway of the overlapped DEGs and may be important in DKD pathogenesis. CD44, also termed CD44 Antigen or ECMR-III, is a major cell-surface receptor for ECM and hyaluronate (28). A previous study reported that the protein of CD44 was downregulated in patients with kidney transplants (29). The proximal tubule (PTC) was identified to be important in renal interstitial fibrosis in DN and CD44 is the receptor of hyaluronic acid in the process of PTC hyaluronan synthesis (30). FN1 is a glycoprotein in the ECM and has been associated with DKD metabolism (31). As a long non-coding RNA (lncRNA), Pvt1 oncogene (PVT1) has been reported to affect the renal metabolism and is proposed as a regulator of DN (31). ECM degradation has an active role in promoting the glomerular fibrosis, whereas PVT1 inhibited the degradation of ECM proteins in the glomerulus. In mesangial cells, knockdown of this lncRNA may lead to an evident reduction of FN1 (31,32). A previous study confirmed the deregulation of FN1 in DN, and demonstrated the reduced expression level of this gene is regulated by the miR-1207-5p (33). These findings collectively support the results of the current study that CD44 and FN1 were important genes in DKD that enriched in the ECM-receptor interaction pathway and were highlighted in the PPI network. It is possible that the two genes may have a significant effect on DKD progression via involvement in this ECM-associated pathway.

The biology of cytokines is often complicated; they serve a variety of roles via interaction with expressed receptors, to trigger signaling pathways and subsequent cytokine-driven outcomes (34). Similar to the ECM pathway, the cytokine-cytokine receptor interaction pathway is also important in DN based on a network-centric analysis (35). In addition, this pathway is significantly enriched for upregulated genes in T2D (36). A previous study determined that mediated by Jak-STAT, this pathway may interact with TGF-β signaling, which as aforementioned, has a role in renal fibrosis (37). All these findings support the significant role of this pathway in DKD. Chemokines have been identified to promote the progression of DKD. For example, the pro-inflammatory chemokine (C-C motif) ligand 12 contributed to the glomerulosclerosis in mice with T2D and blocking its expression had a protective effect on DKD (38). The inflammatory CCL5 is expressed in various cell types, including fibroblasts and renal tubular epithelial cells. Upregulated CCL5 was previously identified in the kidney and its expression was directly associated with concentration of proteinuria in renal tubular cells (1). As one of the CXC chemokine receptors (CXCRs), CXCR4 is specific for stromal cell-derived factor-1. The CXCR4 gene was identified to be increased in injured kidneys (39), and was associated with the cytokine-cytokine receptor interaction pathway in various diseases (40,41). In the present study, CCL5 and CXCR4 were predominant in the PPI network and enriched in the cytokine-cytokine receptor interaction pathway, supporting the involvement of chemokines and the receptors in DKD progression.

In terms of the miRNAs, the miR-17~92 cluster is the focus in renal disease (42). A previous study validated that miR-17-5p and miR-20a were associated with nephron progenitor proliferation and they are upregulated in polycystic kidney disease (43). Additionally, miR-20a is upregulated by TGF-β signaling and implicated in tumorigenesis (44). As mentioned previously, TGF-β has an important role in renal fibrosis, and it may be inferred that the TGF-β-induced miR-20a may also be important in DKD pathogenesis. MiR-106 has also been previously implicated in diabetic retinopathy (45), however, the role in renal diseases has not been identified. These findings suggested the importance of the three miRNAs, miR-17-5p, miR-20a and miR-106, in DKD progression.

NR4A belongs to the nuclear transcription factor family that act as key regulators of cytokines. NR4A is suggested to have an important role in renal inflammation (46). As one of the NR4A members, NR4A3 was confirmed to be differentially expressed in the renal medulla of hypertensive patients, compared with normotensive controls (46). A previous study confirmed that in cell transformation, NR4A3 is a target of miR-17-5p (47). This provides a foundation of the targeting relationship between NR4A3 and miR-17-5p in DKD, as predicted in the integrated network of the present study. In addition, NR4A3 may be a novel target of miR-20a and miR-106, based on the predictions of the current study.

PTPRO encodes for the protein tyrosine phosphatase receptor type O, a member of the R3 subtype family. Dysregulation of this protein leads to childhood-onset nephrotic syndrome (48), and antibodies to this protein lead to an increase in the glomerular albumin permeability (49). The Wnt/β-catenin (Ctnnb1) pathway is important to DKD and genes or proteins in this pathway have been upregulated in glomeruli and podocytes. It is of note that, in heterozygous NPHS2Cre/Ctnnb1FloxE3/WT mice, expression of PTPRO was significantly reduced (50), suggesting that this gene may participate in the Ctnnb1 pathway. It has been previously established that miR-17~92 cluster members have oncogenic function in human cancers and PTPRO acts as a tumor suppressor. The evidence provided by Xu et al (51) confirmed that PTPRO was regulated by the miR-17~92 cluster in tumor (48). A previous review that focused on the role of miRNA in the pathogenesis of DN indicated that PTPRO is predicted a target of miR-25 (52). Based on our study, PTPRO may be a novel target of miR-17-5p, miR-20a and miR-106a.

KLF9 was upregulated and identified as a biomarker of DN by integrating three microarray datasets (53). However, its function in DKD remains to be elucidated. A previous study reported that KLF9 was upregulated by oxidative stress and could promote this stress induced cell death (54). It has been previously identified that as a target of miR-106 under thyroid hormone treatment (55). In the present study, KLF9 was predicted as the target gene of miR-17-5p, miR-20a and miR-106a.

Lack of the experimental validation is a limitation of the present study, particularly the validation of miRNA-target regulations these will be performed in future studies. The present study is a valuable contribution to the prediction of important miRNAs and the reveal of pathological mechanisms in DKD.

In conclusion, the present study identified several important genes in DKD, which may be involved in pathways, such as ECM-receptor and cytokine-cytokine receptor interactions. Three potential miRNAs biomarkers were identified, including miR-17-5p, miR-20a and miR-106a, with the predicted targets of NR4A3, PTPRO and KLF9. However, additional validation is required to confirm these findings.