Investigation into the underlying molecular mechanisms of hypertensive nephrosclerosis using bioinformatics analyses.

Hypertensive nephrosclerosis (HNS) is a major risk factor for end-stage renal disease. However, the underlying pathogenesis of HNS remains to be fully determined. The gene expression profile of GSE20602, which consists of 14 glomeruli samples from patients with HNS and 4 normal glomeruli control samples, was obtained from the Gene Expression Omnibus database. Gene ontology (GO) and pathway enrichment analyses were performed in order to investigate the functions and pathways of differentially expressed genes (DEGs). Pathway relation and co‑expression networks were constructed in order to identify key genes and signaling pathways involved in HNS. In total, 483 DEGs were identified to be associated with HNS, including 302 upregulated genes and 181 downregulated genes. Furthermore, GO analysis revealed that DEGs were significantly enriched in the small molecule metabolic process. In addition, pathway analysis also revealed that DEGs were predominantly involved in metabolic pathways. The tricarboxylic acid (TCA) cycle was identified as the hub pathway in the pathway relation network, whereas the sorbitol dehydrogenase (SORD) and cubulin (CUBN) genes were revealed to be the hub genes in the co‑expression network. The present study revealed that the SORD, CUBN and albumin genes as well as the TCA cycle and metabolic pathways are involved in the pathogenesis of HNS. The results of the present study may contribute to the determination of the molecular mechanisms underlying HNS, and provide insight into the exploration of novel targets for the diagnosis and treatment of HNS.

Introduction

Hypertensive nephrosclerosis (HNS), also termed hypertensive nephropathy or ‘benign’ nephrosclerosis, is one of the major risk factors for end-stage renal disease (ESRD) (1). The development of HNS has previously been reported to be closely associated with hypertension, race and aging (2–4), for which the diagnosis is predominantly based on clinical manifestations (5). However, despite numerous studies investigating HNS, the underlying pathogenic mechanisms of this disease remain to be elucidated, which limits the ability to pertinently treat and improve prognosis. Therefore, clarifying the molecular mechanisms of HNS is imperative for the development of effective diagnostics and therapeutic strategies.

In a disease network, a systems biology approach can be adopted as a means to reveal complex molecular interactions, rather than single molecular components (6). The ‘omics’ analyses, which represent major cornerstones of systems biology research, are considered to be unbiased methods for the identification of biomarkers and to elucidate the pathological mechanisms of chronic kidney disease (7). These analyses are invaluable tools in nephrology research, and greatly facilitate the work of scientists (8). Bioinformatics has the potential to enable scientists to comprehensively analyze altered mRNA or microRNA expression patterns in a disease and identify key genes and pathways via construction of correlative networks (9). This analytical approach has been widely used in order to reveal the potential mechanisms of renal diseases, including lupus nephritis (10), membranous nephropathy (11) and diabetic nephropathy (12). The findings of these analyses have greatly contributed to the development of knowledge with regards to renal diseases. However, to the best of our knowledge, the bioinformatics network analysis of HNS has not previously been reported.

The present study aimed to use different bioinformatics approaches in order to determine the differentially expressed genes (DEGs) in HNS. The original GSE20602 dataset was downloaded in order to identify the DEGs between glomeruli specimens from patients with HNS and normal glomeruli specimens. This dataset was published by Neusser et al (13). Subsequently, pathway enrichment and network analysis were performed in order to identify the key genes and signaling pathways implicated in HNS. The results of the present study may improve understanding of the pathogenesis of HNS and be valuable for future studies investigating HNS.

Materials and methods

Microarray data

The gene expression profile GSE20602 was downloaded from the Gene Expression Omnibus (GEO, www.ncbi.nlm.nih.gov/geo/) database, and it was based on the platform of GPL96 [HG-U133A] Affymetrix Human Genome U133A Array. The dataset included 14 renal biopsy samples from patients with HNS and four healthy control samples. Identification of HNS samples was based upon light and electron microscopy analysis and associated predetermined clinical and histological criteria, as well as immunofluorescence examination (14–16). The control samples were derived from normal kidney tissues of patients with tumor nephrectomy (13). The patient glomeruli were microdissected from biopsy tissues for RNA isolation and further microarray experiments (13).

Data pre-processing and differential analysis

Based on the Gene-Cloud Biotechnology Information (www.gcbi.com.cn) platform, the raw CEL files were transformed into probe-level data, and subsequently converted into corresponding gene symbols. According to the algorithms described by Bolstad (17), the Robust Multi-chip Average method was used in order to compute expression levels of probes, which consisted of three steps: Background correction, data normalization and expression measure. Briefly, the background corrected intensities on every GeneChip were computed for every perfect match cell. Following this, normalization data was then acquired using the quantile normalization algorithm, among which individual values were replaced with the mean of original values if the values were the same as another perfectly matched cell. Finally, the expression level for each probe was determined via an additive linear model. The Student's t-test was used in order to calculate the P-values of genes, and Hochberg's method (18) was used to adjust the raw P-value via calculation of the false discovery rate (FDR). Only genes with |log2 fold-change (FC)|>1 and FDR<0.05 were selected as DEGs for further investigation (19). Hierarchical clustering analyses of DEGs were performed as previously described (20).

Gene ontology (GO) and pathway enrichment analyses

GO analysis may be applied for the annotation of genes of high-throughput genomic or transcriptomic data (21). It is capable of predicting the function of genes in three aspects, including biological processes, molecular function and cellular components. The Kyoto Encyclopedia of Genes and Genomes (KEGG, www.genome.jp) is a recognized pathway-associated database for the systematic analysis of gene function (22). The Database for Annotation, Visualization and Integrated Discovery (DAVID, david.ncifcrf.gov) is an online bioinformatics resource for the systematic extraction of biological function from large gene or protein lists (23). In the present study, DAVID was applied in order to conduct GO and KEGG pathway enrichment analyses of the 483 identified DEGs. Fisher's exact test was used to calculate the P-value, and the FDR was calculated to correct this. The P<0.05 and FDR<0.05 were set as the significance threshold.

Network construction and hub module identification

Hub nodes have increased complex correlativity compared with other genes within the networks, therefore the are more likely to be involved in the underlying mechanisms of disease (24). Pathway relation network analysis can simultaneously reveal the pathway that has the greatest regulatory effect on both the highest and lowest stream pathways. Gene co-expression network analysis determines the association between genes, and can aid in the search for the key gene from complex regulatory associations. The networks in the present study were constructed by Genminix Informatics Co., Ltd. (Shanghai, China) (25). The algorithms heavily reference previously published methods (26). Briefly, in networks, the nodes represent genes or pathways, and the edges indicate interactional relationships among them. The centrality of a network is represented by the central degree, which is the contribution of one gene (or pathway) to the genes (or pathways) in the vicinity and is represented by the area of nodes. The greater the degree value, the greater the area of the node. Therefore, the key genes and pathways may be distinctly identified from the networks.

Results

Clinical and histological characteristics

Of the 14 patients with HNS, 11 were males, and 3 were females. The mean age was 58 years. The mean systolic blood pressure was 143 mmHg, and the mean diastolic blood pressure was 86 mmHg. The mean creatinine, estimated glomerular filtration rate (eGFR) and proteinuria level were 2.3 mg/dl, 42 ml/min and 1.6 g/24 h, respectively. In the control group, the mean age, creatinine and eGFR level were 65 years, 1.0 mg/dl and 59 ml/min, respectively. The mean systolic and diastolic blood pressures were not available for this group. No proteinuria was detected in the control group (Table I). The clinical and histological characteristics of patients with HNS and controls were presented in a study by Neusser et al (13).

Table I.

Characteristics of patients with HNS (n=14) and controls (n=4). Data are presented as mean ± standard deviation.

Category	HNS	Control
Gender (male/female/NA)	11/3/0	0/2/2
Age (years)	58±12	65±9
BP systolic (mmHg)	143±18	NA
BP diastolic (mmHg)	86±13	NA
Creatinine (mg/dl)	2.3±1.7	1.0±0
Estimated glomerular filtration rate (ml/min)	42±24	59±1
Proteinuria (g/24 h)	1.6±1.7	0

BP, blood pressure; NA, not applicable; HNS, hypertensive nephrosclerosis.

Identification of DEGs

According to the cut-off criteria of |log2FC|>1 and FDR<0.05, a total of 483 DEGs were identified in HNS samples compared with control samples, including 302 upregulated genes and 181 downregulated genes. The results of DEG identification are presented in the form of a heat map and a volcano plot (Fig. 1), and the top 10 dysregulated genes are presented in Table II. As revealed by Fig. 1A, the 4 control samples form a sub-cluster with some of the HNS samples. The most significantly upregulated and downregulated genes are albumin (ALB; FC=9.74) and DEAD-box helicase 3, Y linked (DDX3Y; FC=−9.14) (Table III).

Figure 1.

DEGs between HNS and control samples. (A) Heat map of the DEGs, red indicates upregulated genes and green indicates downregulated genes. (B) Volcano plot of the DEGs (Yellow indicates DEGs with |log2 fold-change|>1 and false discovery rate <0.05. Blue indicates non-DEGs). HNS, hypertensive nephrosclerosis; DEG, differentially expressed genes.

Table II.

Top 10 differentially expressed genes associated with hypertensive nephrosclerosis.

Gene symbol	Gene description	Fold-change	P-value	False discovery rate
ALB	Albumin	9.74	1.42×10−3	1.24×10−3
DDX3Y	DEAD-box helicase 3, Y-linked	−9.14	2.18×10−2	2.11×10−2
EGR1	Early growth response 1	6.78	4.60×10−5	0
FOS	FBJ murine osteosarcoma viral oncogene homolog	6.06	1.90×10−4	0
HPD	4-hydroxyphenylpyruvate dioxygenase	6.04	3.95×10−3	1.83×10−3
RPS4Y1	Ribosomal protein S4, Y-linked 1	−5.92	3.23×10−2	3.17×10−2
ABP1	Amiloride binding protein 1 (amine oxidase (copper-containing))	5.62	8.38×10−3	5.04×10−3
CYP4A11	Cytochrome P450, family 4, subfamily A, polypeptide 11	5.48	2.18×10−3	1.24×10−3
ALDOB	Aldolase B, fructose-bisphosphate	5.37	7.44×10−3	5.04×10−3
BBOX1	γ-butyrobetaine hydroxylase 1	5.05	4.00×10−3	3.23×10−3

Table III.

Top 10 significant gene ontology terms of differentially expressed genes associated with hypertensive nephrosclerosis.

Gene ontology name	Differentially expressed genes	P-value	False discovery rate
Small molecule metabolic process	76	1.10×10−37	2.14×10−34
Cellular nitrogen compound metabolic process	21	6.42×10−17	6.25×10−14
Transmembrane transport	30	4.48×10−15	2.91×10−12
Xenobiotic metabolic process	16	2.80×10−13	1.36×10−10
Response to drug	20	1.72×10−12	6.71×10−10
Cellular response to calcium ion	9	7.68×10−12	2.49×10−9
Platelet degranulation	12	1.68×10−11	4.67×10−9
Platelet activation	16	1.16×10−10	2.82×10−8
Steroid metabolic process	10	1.36×10−10	2.95×10−8
Gluconeogenesis	9	3.41×10−10	6.64×10−8

GO and KEGG pathway enrichment analyses

All DEGs were uploaded to the DAVID software in order to perform GO and KEGG analyses. GO analysis revealed that 384 significant GO categories (FDR<0.05) were regulated by DEGs, and the top 10 significantly enriched GO terms are presented in Fig. 2A and Table IV. Among them, the small molecule metabolic process (P=1.10×10−37) was most significantly associated with HNS, which contained 76 DEGs. Pathway analysis demonstrated that 112 pathway categories (FDR<0.05) were affected by DEGs. Both Fig. 2B and Table IV present the top 10 significantly enriched pathways, the most prominent of which were the metabolic pathways (P=1.32×10−31), with 65 DEGs.

Figure 2.

Histogram of the top 10 significant (A) GO categories and (B) pathways of differentially expressed genes associated with hypertensive nephrosclerosis. GO, gene ontology; PPAR, peroxisome proliferator-activated receptor.

Table IV.

Top 10 significant pathways of differentially expressed genes associated with hypertensive nephrosclerosis.

Pathway name	Differentially expressed genes	P-value	False detection rate
Metabolic pathways	65	1.32×10−31	2.69×10−29
PPAR signaling pathway	13	1.29×10−13	1.32×10−11
Protein digestion and absorption	12	4.61×10−11	3.13×10−9
Cell adhesion molecules	14	1.27×10−10	6.50×10−9
Arginine and proline metabolism	10	1.64×10−10	6.68×10−9
Drug metabolism-cytochrome P450	10	2.38×10−9	8.10×10−8
Glycine, serine and threonine metabolism	7	6.81×10−8	1.98×10−6
Proximal tubule bicarbonate reclamation	6	8.44×10−8	2.15×10−6
Tryptophan metabolism	7	1.20×10−7	2.73×10−6
Carbohydrate digestion and absorption	7	2.82×10−7	5.76×10−6

Pathway association network analysis

The pathway association network analysis was conducted based on the interrelation of the aforementioned 112 enriched pathway categories from KEGG, and this analysis revealed that the 59 nodes (representing pathways) and 172 edges (representing line connections between pathways) were assigned to the network (Fig. 3). The top ten hub nodes with higher degrees included the tricarboxylic acid (TCA) cycle, glycolysis/gluconeogenesis, mitogen-activated protein kinase (MAPK) signaling pathway, pyruvate metabolism, pathways associated with cancer, propanoate metabolism, cell adhesion molecules, T-cell receptor signaling pathway, ascorbate and aldarate metabolism and the synthesis and degradation of ketone bodies (Table V).

Figure 3.

Relationship network of significant pathways associated with hypertensive nephrosclerosis. Nodes represent pathways. The size of the nodes represents the associated degree value. Lines with arrows represent interactions between pathways, whereas the pathways with arrow tails regulate pathways with the associated arrow's head. Red nodes represent upregulated pathways, blue nodes represent downregulated pathways, and yellow nodes indicate that there were genes associated with the pathway that were both up and downregulated. PPAR, peroxisome proliferator-activated receptor; MAPK, mitogen-activated protein kinase; ECM, extracellular matrix.

Table V.

Top 10 key pathways according to the degree.

Pathway name	Out degree	In degree	Total degree
Tricarboxylic acid cycle	6	13	19
Glycolysis/gluconeogenesis	4	12	16
MAPK signaling pathway	0	16	16
Pyruvate metabolism	7	8	15
Pathways in cancer	12	0	12
Propanoate metabolism	6	5	11
Cell adhesion molecules	4	7	11
T cell receptor signaling pathway	3	8	11
Ascorbate and aldarate metabolism	4	4	8
Synthesis and degradation of ketone bodies	4	4	8

Out degree, number of downstream pathways; In degree, number of upstream pathways; Total degree, sum of out degree and in degree values.

Gene co-expression network analysis

The gene co-expression network was constructed with respect to gene function associations (25) and contained 51 nodes (representing the DEGs) and 168 edges (representing line connections between nodes; Fig. 4). Based on the connectivity degree, the hub genes were revealed to be sorbitol dehydrogenase (SORD), glycine-N-acyltransferase (GLYAT), cubilin (CUBN), pipecolic acid oxidase (PIPOX), ureidopropionase β (UPB1), 4-hydroxyphenyl pyruvate dioxygenase (HPD), dihydropyrimidinase (DPYS), solute carrier family 7 member 9 (SLC7A9), kynurenine 3-monoxygenase (KMO) and serine hydroxymethyltransferase 1 (SHMT1; Table VI).

Figure 4.

Gene co-expression network of differentially expressed genes associated with hypertensive nephrosclerosis. Nodes denote genes; edges represent gene-gene interactions and the size of the nodes represents the associated degree value.

Table VI.

Top 10 key genes according to the degree value.

Gene	Description	Degree
SORD	Sorbitol dehydrogenase	22
GLYAT	Glycine-N-acyltransferase	21
CUBN	Cubilin (intrinsic factor-cobalamin receptor)	15
PIPOX	Pipecolic acid oxidase	14
UPB1	Ureidopropionase, β	14
HPD	4-hydroxyphenylpyruvate dioxygenase	14
DPYS	Dihydropyrimidinase	14
SLC7A9	Solute carrier family 7 member 9	14
KMO	Kynurenine 3-monooxygenase (kynurenine 3-hydroxylase)	13
SHMT1	Serine hydroxymethyltransferase 1 (soluble)	13

Degree, number of the associated gene interactions with the key gene in the network.

Discussion

HNS is a common disease of the kidney, which significantly impacts patient quality of life. However, the exact molecular mechanisms underpinning HNS remain to be determined. High-throughput technologies can simultaneously reveal the expression levels of thousands of molecules and thus can be used for the prediction of potential therapeutic targets for kidney diseases (27).

In the present study, gene expression data was extracted from the GSE20602 dataset in order to identify the underlying molecular mechanisms of HNS via application of numerous bioinformatics approaches. A total of 483 DEGs, including 302 upregulated and 181 downregulated genes, were selected for in the present study. ALB demonstrated the highest FC of all identified DEGs in HNS. ALB encodes albumin, the most abundant protein in human blood. Notably, urinary albumin has previously been demonstrated to be a risk factor for HNS (28). Gupta et al (29) revealed that patients with nephrotic proteinuria and a serum albumin >35 g/l suffered from HNS and had poor renal survival. Furthermore, a recent study using mice demonstrated that clinical outcomes of kidney disease were significantly improved following gene knockdown of ALB, which therefore suggested that filtered albumin is deleterious to kidney cells (30).

In the present study, GO analysis was conducted with the aim of improving current understanding of the main functions of DEGs with regards to HNS. The results of this analysis yielded 384 significant GO terms, including small molecule metabolic process, cellular nitrogen compound metabolic process, transmembrane transport and the xenobiotic metabolic process. Among them, the small molecule metabolic process was the most significant. It is one of the GO terms that belongs to the biological process domain. This gene function category includes thousands of small molecules, and GO analysis revealed that there were 76 DEGs enriched in the small molecule metabolic process term. These metabolic processes are considered to be involved in homeostasis, and individuals with abnormal metabolic statuses are at a significantly increased risk of developing chronic kidney disease (31). However, due to the differences of data preprocessing and filtering criteria, the DEGs associated with HNS revealed by the present study differ from those obtained by Neusser et al (13). In terms of GO analysis, Neusser et al (13) predominantly focused on the role of hypoxia in nephrosclerosis, thus the results of the GO analysis were hypoxia-associated biological processes, such as angiogenesis, inflammation and renal fibrosis. The present study, however, aimed to reveal the pathogenesis of HNS via global bioinformatics analysis; therefore, the GO findings of the present study are inconsistent with those determined by Neusser et al (13).

Similar to the GO analysis, pathway enrichment analysis was also performed to further investigate the DEGs in HNS. It was revealed that the DEGs were predominantly involved in metabolic pathways and the peroxisome proliferator-activated receptors (PPAR) signaling pathway. The kidney is an important metabolic organ; therefore, the metabolic pathway that included 65 DEGs was proven to be the most significant pathway in the present study. PPARs are nuclear hormone receptors and are critical for lipid metabolism (32). A previous study demonstrated that PPARs are highly expressed in the kidney (33). The PPAR signaling pathway was previously revealed to be a common pathway associated with renal dysfunction, such as hypertensive nephropathy (34). However, GO and pathway enrichment are preliminary analyses of DEGs and whether DEGs are implicated in the progression of HNS remained to be verified by network analysis.

Following pathway analyses, 59 significantly enriched pathways were selected in order to establish a pathway relation network. The TCA cycle, glycolysis/gluconeogenesis, MAPK signaling pathway and pyruvate metabolism were identified as the core pathways. The TCA cycle is a mitochondria-dependent process, and the cell bodies of podocytes contain a substantial number of mitochondria (35). Previous studies have revealed that the TCA cycle functions as a bridge in order to connect other metabolic pathways to one another (36,37). Dysfunction of the TCA cycle was previously demonstrated to be associated with kidney injury (38). In the present study, 3 DEGs phosphoenol pyruvate carboxykinase 1 (PCK1), PCK2 and oxoglutarate dehydrogenase like (OGDHL) were revealed to be involved in the TCA cycle. The PCK1 gene encodes the cytosolic isozyme of phosphoenolpyruvate carboxykinase (PEPCK), whereas PCK2 encodes the mitochondrial isozyme of PEPCK (39). PEPCK is a rate-limiting enzyme of gluconeogenesis occurring in the liver and renal cortex, and it is essential for glucose homeostasis (40). Numerous studies have revealed that PCK1 is a multi-functional gene and is implicated in physiological processes in the liver, kidney and adipose tissues (41–43). Therefore, the TCA cycle and its associated pathways and DEGs may take part in the regulation of HNS development.

Finally, a gene co-expression network analysis was conducted in order to reveal hub genes associated with HNS. A number of prominent genes were identified as a result of their degree value. The majority of these genes encoded metabolic enzymes, such as SORD, GLYAT, PIPOX, UPB1, HPD, DPYS, KMO and SHMT1. For example, SORD is the second enzyme of the polyol pathway (involved in glycolysis), which catalyzes the conversion of sorbitol to fructose and is highly expressed in the kidney (44). Due to the polymorphic variation of the SORD gene, both the accumulation and toxicity of sorbitol are associated with the development of microvascular problems (45), which may be one of underlying risk factors for HNS. Furthermore, SORD and aldolase B, another upregulated gene, are both implicated in the sorbitol pathway, which is closely associated with hyperglycemia (46). Furthermore, synthesis and degradation of ketone bodies was identified as one of the hub pathways, which is also associated with diabetes. Whether the results of the present study are also affected by hyperglycemia remains to be investigated by further studies.

In addition to metabolism-associated genes, CUBN was also identified as a hub gene (Table VI). CUBN encodes cubilin, a proximal tubular epithelial cell protein, which was also revealed to be expressed in glomerular podocytes (47). Mutations in CUBN have been associated with the susceptibility to ESRD (48). Considering the association between HNS and ESRD, it may be suggested that CUBN may be involved in the progression of HNS. Furthermore, the megalin protein, encoded by LRP2, is implicated in the facilitation of the internalization of the cubilin-albumin complex (49) and is expressed in the proximal tubule and the glomerulus (50). LRP2 was also screened as a DEG, which confirmed the potential role of CUBN in the pathogenesis of HNS. Further studies on these hub genes may contribute to the development of more effective therapeutic approaches for patients with HNS. However, due to limited experimental conditions, the results of the current study were predicted using only bioinformatics approaches. Further molecular biology studies are required in order to verify these results.

In conclusion, the present study aimed to investigate the potential underlying molecular mechanisms of HNS using bioinformatics analysis. A total of 483 DEGs were identified in HNS samples compared with control samples. Furthermore, the present study revealed that SORD, CUBN and ALB genes, as well as the TCA cycle and metabolic pathways, may be implicated in the pathogenesis of HNS. These results may prove valuable for further studies aiming to investigate novel targets for the diagnosis and treatment of HNS.