Investigation of hub genes involved in diabetic nephropathy using biological informatics methods.

BACKGROUND: The aim of this study was to find genes with significantly aberrant expression in diabetic nephropathy (DN) and determine their underlying mechanisms.
METHODS: GSE30528 and GSE1009 were obtained by querying the Gene Expression Omnibus (GEO) database. The difference in target gene expression between normal renal tissues and kidney tissues in patients with DN was screened by using the GEO2R tool. Using the Database for Annotation, Visualization, and Integrated Discovery (DAVID) database, differentially expressed genes (DEGs) were analysed by Gene Ontology (GO) annotation and the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment. Then, the protein-protein interactions (PPIs) of DEGs were analyzed by Cytoscape with the Search Tool for the Retrieval of Interacting Genes/Proteins (STRING) database, and the hub genes in this PPI network were recognized by centrality analysis.
RESULTS: There were 110 genes with significant expression differences between normal and DN tissues. The differences in gene expression involved many functions and expression pathways, such as the formation of the extracellular matrix and the construction of the extracellular domain. The correlation analysis and subgroup analysis of 14 hub genes and the clinical characteristics of DN showed that CTGF, ALB, PDPN, FLT1, IGF1, WT1, GJA1, IGFBP2, FGF9, BMP2, FGF1, BMP7, VEGFA, and TGFBR3 may be involved in the progression of DN.
CONCLUSIONS: We confirmed the differentially expressed hub genes and other genes which may be the novel biomarker and target candidates in DN. 2020 Annals of Translational Medicine. All rights reserved.

Introduction

Diabetic nephropathy (DN) is a serious impairment of renal function caused by poor blood glucose control in diabetic patients. The clinical incidence of DN is high, treatment is difficult, and prognosis is poor (1). The typical clinical manifestations of DN include a decreased glomerular filtration rate and persistent proteinuria. Patients with severe DN will progress to uremia within 2–3 years of onset, and then many serious complications and death will occur (2). Pathophysiological studies have found that the progression of DN is closely related to genetic factors, metabolic changes, and hemodynamic changes (3).

The microarray technology and bioinformatic analysis allow for the screening of genetic alterations on the genome level. However, independent microarray analysis always lead to false‐positive rates. In this study, we re-analyzed two microarray datasets, GSE30528 and GSE1009, with a different statistical significance detection method which resulted in a dissimilar number of differentially expressed genes (DEGs). Key biomarkers with high specificity and sensitivity were identified by detecting genes with significant differences in expression between DN and normal kidney tissues. The pathogenesis of DN was studied by Gene Ontology (GO), Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment. and protein­-protein interaction (PPI) network analysis. The Nephroseq v5 platform was used to analyse correlations and to perform subgroup analysis among the hub genes and clinical manifestations of DN to further explore the pathogenesis, pathophysiological mechanisms, and key genes involved in DN.

We present the following article in accordance with the MDAR reporting checklist (available at http://dx.doi.org/10.21037/atm-20-5647).

Methods

Microarray data and preprocessing

From the Gene Expression Omnibus (GEO) database, established by Baelde et al., we successfully obtained GSE1009 free microarray data (4). The data in the dataset were derived from three pairs of normally functioning glomerulus and diabetic glomerulus samples that were obtained from autopsy. At the same time, to facilitate data comparisons, we obtained GSE30528 microarray data from this database. These data were first recorded by Woroniecka et al. (5). A total of 22 samples were included, 9 of which were from DN patients, and 13 of which were from normal human kidney tissues. The study was conducted in accordance with the Declaration of Helsinki (as revised in 2013).

Data processing of DEGs

DEGs in kidney tissues of DN and normal subjects were detected by GEO2R online tool (6) with |logFC| >1.5 and P<0.05. GEO2R is an interactive web tool that allows users to compare two or more datasets in a GEO series to identify DEGs across experimental conditions. Subsequently, the original data recorded in the TXT format were validated online by Venn to find significant DEGs in the two datasets; when the test results indicated log FC <0, the expression of the gene was deemed to be downregulated; if log FC >0, gene expression was deemed to be upregulated.

GO and pathway enrichment analysis

The purpose of GO is to unify biological factors while integrating specific definitions, clear structures, and controlled vocabulary into a database (7). Meanwhile, KEGG is a commonly used database for gene research. The KEGG database can conveniently classify all kinds of genes with system path management (8). The Database for Annotation, Visualization, and Integrated Discovery (DAVID) is also a commonly used tool in bioinformatics research. The advantages of DAVID include the detailed analysis and classification of gene and protein functions (9). In our study, we used DAVID to envisage the biological process (BP), molecular function (MF), pathway and cell component (CC), and other DEGs to collect more information in our research.

PPI network analysis and hub gene identification

The PPI network we used was based on Cytoscape software version 3.7.1, which can accurately predict protein-related experimental data and PPI results in the target database (10,11). The results express the interaction between different proteins as a combination fraction. After DEG information was evaluated in the PPI network, if the comprehensive score was more than 0.4, these DEGs were known as an important protein pair. This network system can use the Coexpedia (12) online platform to analyze the hub genes and related gene expression network system. After testing different gene pairs in Clusterone 1.0, PPI network modules with a minimum size of 3 and a minimum density of 0.5 were obtained. A module was considered to have a significant clustering effect when P<0.01.

Statistical analysis

Pearson’s correlation analysis between hub genes, glomerular filtration rate (GFR), and proteinuria in patients with DN was performed using Nephroseq v5 online platform (http://v5.nephroseq.org). Comparisons between two groups were performed using an unpaired Student’s t-test. All tests were two-tailed, with a P value <0.05 considered statistically significant.

Results

Identification of DEGs in DN

The free NCBI-GEO database was used to search the gene expression profiles of GSE30528 and GSE1009 in kidney tissues of DN patients. Among these datasets, the GSE30528 microarray data package contained 22 samples of human glomerular tissue. The study samples were taken from healthy human transplant donors, diagnostic renal biopsies, and post-nephrectomy specimens. Further analysis of the specific composition showed that it included 9 DN kidney tissue samples and 13 control group kidney tissue samples. The GSE1009 microarray data included 6 human glomerular samples, including DN samples and 3 samples in the control group.

The data obtained from the database were analysed by GEO2R. The cut-off criteria were [log FC] >1.5 and P<0.05. We found 274 and 1133 differentially expressed loci in GSE30528 and GSE1009 gene data. Then, Venn diagram software was used to identify the common DEGs in these two data sets. A total of 110 DEGs were detected, 96 of which were downregulated, and 14 of which were upregulated (log FC <0) ( and ).

Table 1

110 differentially expressed genes (DEGs) were identified from two profile datasets

DGEs	Genes name
Upregulated	ALOX5, FRY, CCL19, GRIK2, C1orf21, TAC1, CORO1A, IGLC1, THBS2, IL7R, UMOD, TDO2, TRBC1, IGK///IGKC
Downregulated	PLCE1, MYO1B, IGF1, KLK6, SEPT11, COL4A5, CXADR, GJA1, CDC14A, PALLD, GRK5, GHR, EVI2A, FOXC1, MGAT5, PARVA, CDC42EP3, WT1-AS, FGF1, CHI3L1, CORO2B, NEBL, PLA2R1, COL4A3, NPHS1, WT1, TCF21, TRAM2, KCND3, ST3GAL6, TNNT2, BMP7, DUSP1, KBTBD11, GALC, GPRC5A, CD200, IGFBP2, MEGF9, BMP2, CLIC5, CAND2, HPGD, VEGFA, MME, TYRO3, AIF1SRGAP2C///SRGAP2B///SRGAP2, MIR6872///SEMA3B, CTGF, CMAHP, DPP6, PTGER4, CR1, PDPN, PRKAR2B, FLT1, PTPRD, NFASC, PTGDS, CTNNBIP1, KLK7, USP46, THSD7A, MAGI2, EXPH5, LOC101930416///LOC101929792///LOC100996724///PDE4DIP, FGF9, SPOCK2, PAMR1, DDN, APOD, F2R, HTRA1, CACNB2, NTNG1, PDE4B, CTDSPL, PTPRO, CDKN1C, GAS1, ZNF185, MAP1B, LOX, SULF1, MYL9, TGFBR3, LPL, F5, MCM6, SEMA5A, ALB, F3, ITGB5, DPYSL3, SPOCK1

DN, diabetic nephropathy.

Figure 1

Identification of 110 commonly changed DEGs from the two cohort profile datasets (GSE30528 and GSE1009). Different color areas represent different datasets. The cross areas indicate commonly changed DEGs. DEGs, differentially expressed genes.

DEG GO analysis in DN

DAVID software and GO analysis were used to further study the 110 DEGs that were identified. These DEGs can be divided into three basic types: molecular structure and function, cell component (CC), and BP. MF analysis showed that the main functions of these genes were the regulation of heparin binding and protein tyrosine phosphatase activity. Genes in CC mainly involved nuclear membrane, collagen IV trimer, lamellipodium, postsynaptic density, actin filament, focal adhesion, extracellular matrix of proteins, extracellular region, extracellular exosome, and extracellular space. In BP, the functions of these genes were mainly involved in regulation of cell differentiation, phosphatidylinositol 3-kinase signaling, prostaglandin metabolic process, osteoblast differentiation, angiogenesis, ERK1 and ERK2 cascade, MAPK activity, tachykinin in receptor signaling pathway, prostatic bud formation, cysteine-type endopeptidase activity involved in the apoptotic process, angiogenesis, epithelial cell proliferation, cell migration, response to prostaglandin E, FGF receptor signaling pathway, angiogenesis, and glomerular basement membrane development ().

Table 2

GO enrichment analysis of DEGs in DN

Category	Term	Count	P value
GOTERM_BP_DIRECT	GO:0030154~cell differentiation	3	3.88E-02
GOTERM_BP_DIRECT	GO:0014068~positive regulation of phosphatidylinositol 3-kinase signaling	3	3.88E-02
GOTERM_BP_DIRECT	GO:0006693~prostaglandin metabolic process	2	3.42E-02
GOTERM_BP_DIRECT	GO:0045669~positive regulation of osteoblast differentiation	3	3.13E-02
GOTERM_BP_DIRECT	GO:0016525~negative regulation of angiogenesis	3	3.13E-02
GOTERM_BP_DIRECT	GO:0070374~positive regulation of ERK1 and ERK2 cascade	4	2.98E-02
GOTERM_BP_DIRECT	GO:0000187~activation of MAPK activity	3	2.85E-02
GOTERM_BP_DIRECT	GO:0007217~tachykinin receptor signaling pathway	2	2.75E-02
GOTERM_BP_DIRECT	GO:0060686~negative regulation of prostatic bud formation	2	2.75E-02
GOTERM_BP_DIRECT	GO:0006919~activation of cysteine-type endopeptidase activity involved in apoptotic process	3	2.71E-02
GOTERM_BP_DIRECT	GO:0001525~angiogenesis	4	2.52E-02
GOTERM_BP_DIRECT	GO:0050679~positive regulation of epithelial cell proliferation	3	2.44E-02
GOTERM_BP_DIRECT	GO:0030335~positive regulation of cell migration	4	2.25E-02
GOTERM_BP_DIRECT	GO:0034695~response to prostaglandin E	2	2.07E-02
GOTERM_BP_DIRECT	GO:0008543~fibroblast growth factor receptor signaling pathway	3	1.82E-02
GOTERM_BP_DIRECT	GO:0030324~lung development	4	4.44E-03
GOTERM_BP_DIRECT	GO:0045766~positive regulation of angiogenesis	5	1.39E-03
GOTERM_BP_DIRECT	GO:0032836~glomerular basement membrane development	3	9.76E-04
GOTERM_MF_DIRECT	GO:0004725~protein tyrosine phosphatase activity	4	1.07E-02
GOTERM_MF_DIRECT	GO:0008201~heparin binding	5	1.11E-03
GOTERM_CC_DIRECT	GO:0031965~nuclear membrane	4	3.73E-02
GOTERM_CC_DIRECT	GO:0005587~collagen type IV trimer	2	3.08E-02
GOTERM_CC_DIRECT	GO:0030027~lamellipodium	4	2.32E-02
GOTERM_CC_DIRECT	GO:0014069~postsynaptic density	3	2.11E-02
GOTERM_CC_DIRECT	GO:0005884~actin filament	3	1.59E-02
GOTERM_CC_DIRECT	GO:0005925~focal adhesion	7	7.97E-03
GOTERM_CC_DIRECT	GO:0005578~proteinaceous extracellular matrix	5	6.46E-03
GOTERM_CC_DIRECT	GO:0005576~extracellular region	8	5.65E-03
GOTERM_CC_DIRECT	GO:0070062~extracellular exosome	28	4.23E-05
GOTERM_CC_DIRECT	GO:0005615~extracellular space	18	5.89E-07

GO, Gene Ontology; DEGs, differentially expressed genes; DN, diabetic nephropathy.

Signalling pathway enrichment analysis

KEGG pathway analysis was used to verify DEG functions and signalling pathway enrichment. In all, 9 important signalling pathways were identified, including the regulation of actin cytoskeleton, ECM-receptor interaction, dilated cardiomyopathy, hypertrophic cardiomyopathy (HCM), Ras signaling pathway, pathways in cancer, Rap1 signaling pathway, PI3K-Akt signaling pathway, and focal adhesion ().

Table 3

KEGG pathway enrichment analysis of DEGs in DN

Term	Count	P value
has04810: regulation of actin cytoskeleton	5	4.84E-02
has04512: ECM-receptor interaction	4	2.13E-02
has05414: dilated cardiomyopathy	4	1.66E-02
has05410: HCM	4	1.56E-02
has04014: Ras signaling pathway	6	1.53E-02
has05200: pathways in cancer	9	3.93E-03
has04015: Rap1 signaling pathway	8	3.91E-04
has04510: focal adhesion	9	6.26E-05
has04151: PI3K-Akt signaling pathway	12	1.38E-05

KEGG, Kyoto Encyclopedia of Genes and Genomes; DEGs, differentially expressed genes; DN, diabetic nephropathy; HCM, hypertrophic cardiomyopathy.

PPI network analysis

The construction of the PPI network was based on Search Tool for the Retrieval of Interacting Genes/Proteins (STRING) online analysis (). A total of 110 DEG data points were imported into the PPI network for analysis, covering 188 edges and 79 key nodes (). Subsequently, 14 central nodes were identified by cell-type steroid optimization analysis: ALB, BMP2, FLT1, CTGF, PDPN, IGF1, GIA1, FGF9, BMP7, FGF1, VEGFA, TGFBR3, WT1, and IGFBP2 ().

Figure 2

Common DEGs PPI network constructed by STRING online database and Module analysis. (A) There were a total of 110 DEGs in the DEG PPI network complex. The nodes represent proteins, and the edges represent the interaction of proteins. Green circles represent downregulated DEGs, and red circles represent upregulated DEGs. (B) Module analysis for PPI network of gene signatures via Cytoscape software. Green represents a downregulated gene. DEGs, differentially expressed genes; PPI, protein-protein interaction; STRING, Search Tool for the Retrieval of Interacting Genes/Proteins.

Association between the hub genes and clinical features of DN

Of the common genes identified in this study, 14 were recognized as hub genes (). In , the network of hub genes and related co-expressed genes by Coexpedia is displayed.

Table 4

Gene symbol and full name of 14 hub genes in diabetic nephropathy

Gene symbol	Full name
CTGF	Connective tissue growth factor
ALB	Albumin
PDPN	Podoplanin
FLT1	Fms-related tyrosine kinase 1
IGF1	Insulin-like growth factor 1
WT1	Wilms tumor 1
GJA1	Gap junction protein alpha 1
IGFBP2	Insulin-like growth factor binding protein 2
FGF9	Fibroblast growth factor 9
FGF1	Fibroblast growth factor 1
BMP2	Bone morphogenetic protein 2
BMP7	Bone morphogenetic protein 7
VEGFA	Vascular endothelial growth factor A
TGFBR3	Transforming growth factor beta receptor 3

Figure 3

Hub genes and their co-expressed genes were analyzed using Coexpedia. Nodes with red background represent the co-expressed genes.

Subsequently, we focused on the changes in the hub genes in kidney tissue samples from DN patients. With the help of Nephroseq v5, the correlation between the unknown hub genes and GFR was determined. The expression of BRMP 2, CTGF, FGF9, and iGFBP2 in kidney tissue samples of DN patients was negatively correlated with GFR. Thus, the expression changes of these four genes directly affect the occurrence and development of DN. The expression of the ALB, BMP7, TGFBR3, and VEGFA genes was positively correlated with GFR. Thus, the expression of these four genes contributes to the maintenance and repair of renal function ().

Figure 4

Correlation between mRNA expression of upregulated hub genes in GFR in DN patients. GFR, glomerular filtration rate; DN, diabetic nephropathy.

The Nephroseq v5 results indicated that the expression of the GJA1 gene was positively correlated with proteinuria in DN patients (), indicating that the high expression of GJA1 promoted the progression and deterioration of DN. However, the level of proteinuria in DN patients decreased significantly with the increase of the expression of genes such as VEGFA, CTGF, FLT1, IGF1, and TGFBR 3. There was a negative correlation between proteinuria and the expression of these genes, which indicated that the expression of these genes contributed to the maintenance and repair of renal function.

Figure 5

Association between the expression of unexplored hub genes and proteinuria in DN patients. DN, diabetic nephropathy.

Discussion

In developed countries, the incidence of type 2 diabetes mellitus has risen dramatically, with renal dysfunction and failure having become the main cause of diabetes-related death (13). The occurrence of diabetes mellitus and the occurrence of the later stages of renal disease are influenced by both genetic and environmental factors (14). In our research, we successfully analyzed the expression pattern of target genes from the GEO database and further identified 110 genes with significant differences between DN and normal populations. Bioinformatics tools such as GO and KEGG helped us to examine the interaction between DEGs. The results showed that genes with significant expression differences were mainly involved in cell maturation and differentiation, cell migration, extracellular area, extracellular matrix components, exons, and other processes (15).

KEGG enrichment analysis of DEGs showed that many of these genes were related to regulation of actin cytoskeleton, ECM-receptor interaction, Ras signaling pathway, focal adhesion, and PI3K-Akt signaling pathway. From a pathological point of view, glomerular hypertrophy and proliferation of the renal basement membrane are often the beginning of glomerular injury in patients with DN. Extracellular bodies also play an important role in the progression of DN. The appearance of extracellular vesicles in renal parenchymal cells can be used as a marker of DN (16). The deposition of a large number of proteins in the basement membrane will eventually lead to glomerulosclerosis. The collagen binding and the interaction between ECM receptors are closely related to the degree and progression of glomerulosclerosis.

Fourteen genes, ALB, BMP2, FLT1, CTGF, PDPN, IGF1, GIA1, FGF9, BMP7, FGF1, VEGFA, TGFBR3, WT1, and IGFBP2, showed significant differences in expression and were identified as key genes. Among these genes, the albumin (ALB) gene encodes serum ALB. Serum ALB, known as the “intravascular porter”, can combine with various ions, hormone molecules, and drug molecules. Serum ALB also has many physiological functions, such as stabilizing the colloidal osmotic pressure, and anti-inflammatory and anti-oxidation functions (17). The results of this study showed that the increase in ALB expression led to a significant increase in GFR expression. A study conducted by Zhao et al. showed that serum ALB levels were negatively correlated with cholesterol histopathological changes and urinary protein content (18). Baelde et al. used bioinformatics tools to measure the expression of the VEGF-A gene in glomeruli of type 2 diabetes mellitus patients and found that the expression of this gene was significantly decreased (4). Some studies have shown that this gene is also inhibited in the renal tubules of glomerular interstitial tissue, and the decrease in renal vessels, cortical atrophy, and urinary protein is further aggravated in patients with decreased expression of the VEGFA gene (19,20). Podoplanin (PDPN) can be detected mainly on the glomerular podocyte membrane. The physiological functions of PDPN include the stability of normal cell morphology. The expression of PDPN in the kidney tissue of streptozotocin-induced diabetic mice was found to be severely inhibited during the progression of DN (21). Soluble Flt1 fms-related tyrosine kinase 1 (sFlt1) is an alternatively spliced soluble form of VEGF receptor 1 (VEGFR-1)/Flt-1, and works as a decoy ligand to kinase insert domain receptor (KDR). With the decrease in FLT1 expression, glomerular podocytes change the cytoskeletal structure of glomerular parenchymal cells, destroying the filtration barrier, and then producing a large number of urinary proteins, leading to renal function damage and failure (22). There results showed that there was a positive correlation between the expression of GFR and the expression of vascular endothelial growth factor receptor (VEGFRA) in DN patients, and there was a negative correlation between the expression of GFR and proteinuria. The expression of FIT1 was negatively correlated with proteinuria. The results showed that the expression of VEGFA and FIT1 played an important role in the maintenance and repair of renal function and had a renal protective effect (23). In another study, an Affymetrix microarray was used to analyse the expression of IGF-1 in the kidneys of DN patients. The results showed that the increased expression of IGF-1 in the kidneys of DN patients accompanied the decrease in proteinuria, suggesting that IGF-1 has a protective effect on renal function.

IGFBP-2 expression was detected in glomerular cells, the loop of Henle, and the medullary collecting ducts. This gene could have a role in the growth and maturation of podocytes (24). A study on the IGFBP-2 gene by Narayana et al. showed that the high expression of the IGFBP-2 gene was accompanied by the inhibition of GRF expression and an increase of proteinuria (25). CTGF, a 38-kDa secretory protein that is rich in cysteine, was initially found in the umbilical vein endothelial cells of patients with diabetic complications (26). Further studies showed that protein levels of CTGF were increased in glomeruli, tubules, and interstitial cells of diabetic kidneys and were negatively correlated with the expression of GFR. Previous studies have suggested that WT-1 is very important for the development of the kidneys and is a typical protective factor for renal function (27). When kidney tissue is damaged or renal function is damaged, the expression of BMP-7 is absent. Wang et al.’s research confirmed that the increased expression of BMP-7 can damage renal function in patients with DN (28). At the same time, the microarray analysis indicated that BMP-7 expression was positively correlated with GRF. The main functions of fibroblast growth factor-1 include the regulation of mitosis and insulin sensitization. We also observed a significant reduction in the renal expression of FGF1 in patients with DN by microarray analysis. Wu et al. found that when fibroblast growth factor-1 expression increased, pathological phenomena such as the overexpression of alpha-SMA and the accumulation of collagen were reversed, improving renal function (29). Therefore, fibroblast growth factor-1 also potentially has a protective effect on renal function.

Thus far, few if any studies have mentioned the interaction of GJA-1, BMP-2, TGFBR-3, and FGF-9 hub genes. In fact, the change in BMP-2 expression in the kidney tissue of DN patients is significantly and negatively correlated with the expression of GFR. Some scholars have reported that a BMP-2 binding site (30) could be detected in renal tubular epithelial tissue. Presently, there is abundant evidence that BMP-2 expression is closely related to renal lesions and that BMP-2 expression is negatively correlated with GFR in DN patients. Fibroblast growth factor-9 and BMP-7 can be found in nasopharyngeal carcinoma tissues (31). The expression of fibroblast growth factor-9 was negatively correlated with GFR in the kidney tissue samples of DN patients. Gap junction protein-alpha 1 (GJA-1) is one of components in cell junction structure. Some molecular biology studies have found that these proteins play an important role in the transport of small molecules between cells (32). There was a positive correlation between GJA-1 and proteinuria in DN patients. We speculate that GJA-1, BMP-2, and FGF-9 are closely related to the progression of renal injury in DN patients. In some studies, TGFBR-3 had effects on both tumorigenesis and tumor suppression (33). The results showed that TGFBR-3 in the renal tissue of DN patients was positively correlated with GFR and was negatively correlated with proteinuria. Therefore, TGFBR-3 may also be an effective protective factor for renal function.

Conclusions

A total of 110 DEGs and 14 hub genes were identified as biomarkers for the clinical diagnosis of DN and as potential targets for new therapies. However, the specific molecular mechanism and biological functions of these genes still require further study.

The article’s supplementary files as