Identification of Hub Genes and Pathways in a Rat Model of Renal Ischemia-Reperfusion Injury Using Bioinformatics Analysis of the Gene Expression Omnibus (GEO) Dataset and Integration of Gene Expression Profiles.

BACKGROUND This study aimed to identify hub genes and pathways in a rat model of renal ischemia-reperfusion injury (IRI) using bioinformatics analysis of the Gene Expression Omnibus (GEO) microarray dataset and integration of gene expression profiles. MATERIAL AND METHODS GEO software and the GEO2R calculation method were used to analyze two mRNA profiles, including GSE 39548 and GSE 108195. The co-expression of differentially expressed genes (DEGs) were identified and searched in the DAVID and STRING databases for pathway and protein-protein interaction (PPI) analysis. Cytoscape was used to draw the PPI network. DEGs were also analyzed using the Molecular Complex Detection (MCODE) algorithm. Cytoscape and cytoHubba were used to analyze the hub genes and visualize the molecular interaction networks. Rats (n=20) included the IRI model group (n=10) and a control group (n=10). Quantitative reverse transcription-polymerase chain reaction (qRT-PCR) was used to measure and compare the expression of the identified genes in rat renal tissue in the IRI model and the control group. RESULTS Ten hub genes were identified, STAT3, CD44, ITGAM, CCL2, TIMP1, MYC, THBS1, IGF1, SOCS3, and CD14. Apart from IGF1, qRT-PCR showed that expression of these genes was significantly increased in renal tissue in the rat model of IRI. The HIF-1alpha signaling pathway was involved in IRI in the rat model, which was supported by MCODE analysis. CONCLUSIONS In a rat model of renal IRI, bioinformatics analysis of the GEO dataset and integration of gene expression profiles identified involvement of HIF-1alpha signaling and the STAT3 hub gene.

Background

Ischemia-reperfusion injury (IRI) occurs when previously ischemic tissue is reperfused with oxygenated blood and is a condition that results in cell death and tissue damage [1]. The restoration of blood flow to ischemic tissue and organs results in increased vascular permeability, tissue edema, endothelial cell injury, local and systemic inflammation, followed by healing and fibrosis. Renal IRI is a cause of acute renal injury, and most commonly occurs following trauma or surgery associated with excessive blood loss, and shock. Renal ischemia can result in transient injury, but the consequences of renal IRI can be severe, resulting in renal failure, which can have a poor clinical prognosis [2].

Recently, the molecular mechanisms associated with IRI have been studied, mainly in cardiac IRI and in animal models. Molecular and bioinformatics studies identified several genes that were upregulated in IRI, including SPRR2F, SPRR1A, and MMP-10 [3], MTA_TC0600002772.mm, MTA_TC1300002394.mm, U7 small nuclear RNA (Rnu7), and RGD7543256_1 [4]. These upregulated genes expressed in tissues following IRI may be associated with molecular and signaling pathways associated with tissue damage, apoptosis, tissue repair, and inflammation [3,4]. However, there have been few studies on the expression genes associated with renal IRI. Studies using animal models of renal IRI would facilitate the analysis of gene expression and molecular pathways that may increase the understanding of the etiology and possible ways to protect the kidney from the effects of IRI [5].

Therefore, this study aimed to identify hub genes and pathways in a rat model of renal IRI using bioinformatics analysis of microarray datasets from Gene Expression Omnibus (GEO) and the integration of gene expression profiles. This study design allows the screening and identification of differentially expressed genes (DEGs) in the renal (IRI) rat model and comparison with normal controls [6]. The Gene Ontology (GO) tool was used to identify the connections with the hub genes and the protein-protein interaction (PPI) network.

Material and Methods

Study design and the identification of gene expression profiles

This study was approved by the Ethics Committee of the Second Hospital of Jilin University. The two gene expression profiles, GSE 39548 and GSE 108195, were identified following a search of the Gene Expression Omnibus (GEO) (GEO query 2.40.0), which contained microarrays of renal tissue samples from a rat model of renal ischemia-reperfusion injury (IRI) and normal rat renal tissue samples. The data in these two profiles were obtained from the GPL96 (Affymetrix Human Genome U133A array) platform and were initially analyzed with GEO2R. R software (version R 3.2.3) was used to process all data using the limma package (limma 3.26.8) and Biobase 2.30.0 [7,8]. The log transformation was applied to the original data, followed by the t-test. GEO2R was used for comparisons of original submitter-supplied data. A fold-change was regarded as significant with P <0.05 and with threshold values >2.0 or <−2.0.

The rat model of renal ischemia-reperfusion injury (IRI)

Male rats, aged 8–10 weeks, were used in the rat model of renal ischemia-reperfusion injury (IRI). The rats (n=20) were randomly divided into the renal IRI model group (n=10) and the normal control group (n=10). Rats in the renal IRI group were anesthetized with isoflurane using 5% for induction and 1.5–2.5% for maintenance of anesthesia. The right kidney was resected, and the renal arteries supplying the left kidney were clamped for 45 minutes. Six hours after reperfusion, the rats were euthanized. The rat kidneys were removed and stored at −80°C for gene expression and protein analysis. The 10 rats in the normal group were anesthetized, and the kidneys were resected for further analysis.

Bioinformatics analysis and quantitative reverse transcription-polymerase chain reaction (qRT-PCR)

The co-expression of differentially expressed genes (DEGs) were identified and searched in the DAVID and STRING databases for pathway and protein-protein interaction (PPI) analysis and Gene Ontology (GO) enrichment analysis, including molecular function (MF), cellular component (CC) and biological process (BP) [9]. Pathway enrichment analysis was performed using the Kyoto Encyclopedia of Genes and Genomes (KEGG) database and DAVID. The inclusion criteria of three categories in the GO enrichment and KEGG pathways statistical significance (p<0.05).

Co-expressed DEGs, including upregulated and down-regulated genes, were searched using the ‘multiple proteins’ query in the STRING database [10], to identify the relevant pathways and functions of the DEGs. The networks were developed using Cytoscape with the basic network setting included the network edges, where the line thickness indicated the strength of the association, and the active interaction sources, including text mining the databases for gene coexpression, neighboring genes, and gene fusion. The relevant information was selected with a confidence score of >0.4 and presented using Cytoscape software [11]. Molecular Complex Detection (MCODE) [12] and cytoHubba were used to predict the hub genes. The parameters of MCODE used in this study were as follows: the degree of cut-off, 2; cluster finding, haircut; node score cut-off, 0.2; k-core, 2; and the maximum depth was 100. The parameters of cytoHubba used in this study were as follows: the top ten nodes ranked by degree; display options, the first-stage nodes; the shortest pathway, and the expanded subnetworks were displayed.

Quantitative reverse transcription-polymerase chain reaction (qRT-PCR)

The top ten genes identified by bioinformatics tools were validated on the rat renal tissues in the model of renal IRI and compared with the normal rat renal tissues by qRT-PCR using the QuantStudio 7 Flex real-time PCR system (Thermo Fisher Scientific, Waltham, MA, USA) [13,14]. The primer sequences used are shown in Table 1. The 2−ΔΔCt method was used, and the PCR reactions were normalized to the GAPDH gene.

Table 1

The primers used in this study.

Gene	Forward primer	Reverse primer
STAT3	CACCTTGGATTGAGAGTCAAGAC	AGGAATCGGCTATATTGCTGGT
CD44	TCGATTTGAATGTAACCTGCCG	CAGTCCGGGAGATACTGTAGC
ITGAM	CCATGACCTTCCAAGAGAATGC	ACCGGCTTGTGCTGTAGTC
CCL2	TAAAAACCTGGATCGGAACCAAA	GCATTAGCTTCAGATTTACGGGT
TIMP1	CGAGACCACCTTATACCAGCG	ATGACTGGGGTGTAGGCGTA
MYC	ATGCCCCTCAACGTGAACTTC	GTCGCAGATGAAATAGGGCTG
THBS1	CCTGCCAGGGAAGCAACAA	ACAGTCTATGTAGAGTTGAGCCC
IGF1	CACATCATGTCGTCTTCACACC	GGAAGCAACACTCATCCACAATG
SOCS3	TGCGCCTCAAGACCTTCAG	GCTCCAGTAGAATCCGCTCTC
CD14	ACTTCTCAGATCCGAAGCCAG	CCGCCGTACAATTCCACAT
GAPDH	AGGTCGGTGTGAACGGATTTG	GGGGTCGTTGATGGCAACA

Results

Expression of differentially expressed genes (DEGs) in the rat model of renal ischemia-reperfusion injury (IRI)

Kidney tissue samples from six rats in the model of renal IRI and kidney tissue from six normal rats were analyzed in the Gene Expression Omnibus (GEO) database, and the results were compared for GSE 39548 and GSE 108195. There were 3,489 DEGs in the GSE 39548 mRNA profile, and 889 DEGs in the GSE 108195 mRNA profile. After integrating the analysis, 232 significantly co-expressed DEGs were identified, including 170 upregulated DEGs and 62 down-regulated DEGs when the IRI rat kidney tissues were compared with matched normal kidney tissue samples.

Bioinformatics analysis of genes and functional pathways in the rat model of renal IRI

Figure 1 and Table 2 show the functions of the identified DEGs. Upregulated DEGs of biological processes were involved in inflammatory responses, and leukocyte migration, the integrin-mediated signaling pathway, and the innate immune response. Down-regulated DEGs of biological processes were involved in neurogenesis, the activation of MAPK activity, and the N-acetylglucosamine metabolic process. The inflammatory responses were the most significant biological processes for the functions of the identified DEGs, and the integrin-mediated signaling pathway had the most gene counts.

Figure 1

Bioinformatics analysis of the interactions between the differentially expressed genes (DEGs) in the rat model of renal ischemia-reperfusion injury (IRI).

Table 2

The significant enriched outcomes in renal ischemia-reperfusion injury (IRI) in the rat model.

Expression	Category	Pathway	Description	Gene count	P-Value
Up-DEGs	BP	Inflammatory response	GO: 0006954	5	2.17E-07
	BP	Leukocyte migration involved in inflammatory response	GO: 0002523	8	3.24E-06
	BP	Integrin-mediated signaling pathway	GO: 0007229	14	1.34E-05
	BP	Positive regulation of gene expression	GO: 0010628	14	3.32E-05
	BP	Innate immune response	GO: 0045087	5	3.41E-05
	CC	Focal adhesion	GO: 0005925	16	6.40E-07
	CC	Extracellular matrix	GO: 0031012	14	7.66E-07
	CC	Cytoplasm	GO: 0005737	82	3.86E-06
	CC	Extracellular space	GO: 0005615	30	8.59E-06
	CC	Extracellular exosome	GO: 0070062	43	1.07E-05
	MF	Integrin-binding	GO: 0005178	9	1.93E-06
	MF	Laminin-binding	GO: 0043236	6	2.19E-06
	MF	Structural constituent of cytoskeleton	GO: 0005200	8	4.00E-06
	MF	Protein-binding	GO: 0005515	55	1.58E-04
	MF	Calcium ion-binding	GO: 0005509	16	7.76E-04
Down-DEGs	BP	Fatty acid metabolic process	GO: 0006631	4	8.20E-03
	BP	Neurogenesis	GO: 0022008	3	1.29E-02
	BP	Activation of MAPK activity	GO: 0000187	3	1.74E-02
	BP	Insulin-like growth factor receptor signaling pathway	GO: 0048009	2	3.14E-02
	BP	N-acetylglucosamine metabolic process	GO: 0006044	2	3.65E-02
	CC	Integral component of plasma membrane	GO: 0005887	9	1.01E-02
	CC	Integral component of membrane	GO: 0016021	27	2.41E-02
	CC	Membrane	GO: 0016020	27	3.01E-02
	CC	Neuronal cell body	GO: 0043025	5	5.55E-02

Upregulated DEGs of the cellular component were involved in cell adhesion and cell cytoplasm, and deregulated DEGs were involved in integral components of the plasma membrane. Focal adhesion of upregulated DEGs was the most significant cellular component, and involvement in the cell cytoplasm had the most gene counts.

The upregulated DEGs of molecular function were involved in protein integrin-binding and laminin binding. Integrin-binding was the most significant molecular function, and protein-binding had the most gene counts. Pathway enrichment analysis showed that the upregulated genes were associated with extracellular matrix (ECM) receptor interactions, proteoglycans in cancer, tumor necrosis factor (TNF) signaling pathways, hematopoietic cell lineages, and the function of phagosomes (Figure 2). Proteoglycans in cancer was the most significant pathway in IRI. However, functional roles for these genes and pathways in IRI were not evaluated, and further studies on the combinations and co-expression DEGs in IRI are required.

Figure 2

Bioinformatics analysis of the gene pathways that were enriched and the upregulated genes in the rat model of renal ischemia-reperfusion injury (IRI). The intensity of the color and the size of the nodes indicate the mean P-value and gene number. Pathway enrichment analysis showed that the upregulated genes were associated with extracellular matrix (ECM) receptor interactions, proteoglycans in cancer, tumor necrosis factor (TNF) signaling pathways, hematopoietic cell lineage, and the function of phagosomes.

Predictive genes in renal IRI in the rat model

STRING database analysis identified 226 nodes and 459 edges that included upregulated and down-regulated DEGs (Figure 3). In STRING analysis of all the co-expressed DEGs, the biological processes, including upregulated and down-regulated expression of DEGs that were involved in response to stimuli. The gene count of 126 of 6,616 (P: 4.67e-12) was similar to the inflammatory response functions identified in DAVID. The molecular function protein-binding gene count was 103 of 6,454 (P: 3.87e-05), and the integrin-binding gene count was 10 of 114 (P: 8.77e-05), which was also similar to the protein integrin-binding identified in DAVID.

Figure 3

Protein-protein interaction (PPI) analysis in the rat model of renal ischemia-reperfusion injury (IRI). STRING database analysis identified 226 nodes and 459 edges that included upregulated and down-regulated differentially expressed genes (DEGs). The red and green nodes represent hub genes. The degree value is shown by an increased size of the hub nodes.

The nodes and edges of the top ten hub genes analyzed by Cytoscape

Cytoscape and cytoHubba analysis identified ten top hub genes, including STAT3, CD44, ITGAM, CCL2, TIMP1, MYC, THBS1, IGF1, SOCS3, and CD14. The STAT3 gene had the highest degree score of 36. DEGs were analyzed using the Molecular Complex Detection (MCODE) algorithm, which showed that MCODE 1 contained 11 gene nodes, including ITGAM, IGF1, MYC, THBS1, TIMP1, LOX, MMP3, CCL2, STAT3, CXCl1, and CD44 with 53 edges. All 11 predictive genes in module 1 underwent Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis, which identified involvement of the HIF-1α signaling pathway (Figure 4).

Figure 4

(A–C) The most significance predictive genes for renal ischemia-reperfusion injury (IRI) in the rat model. The Molecular Complex Detection (MCODE) algorithm showed that MCODE 1 contained 11 gene nodes, including ITGAM, IGF1, MYC, THBS1, TIMP1, LOX, MMP3, CCL2, STAT3, CXCl1, and CD44 with 53 edges. All 11 predictive genes in module 1 were analyzed by Gene Ontology (GO) and the Kyoto Encyclopedia of Genes and Genomes (KEGG) database, which identified the HIF-1α signaling pathway. Differentially expressed genes (DEGs) were analyzed using the Molecular Complex Detection (MCODE) algorithm. MCODE 1: 10.600. MCODE 2: 5.250. MCODE 3: 5.000.

Quantitative reverse-transcription polymerase chain reaction (qRT-PCR) was used to measure and compare the expression of the identified genes in rat renal tissue in the IRI model and the control group (Figure 5). Nine hub genes were significantly increased in the rat model of renal IRI, including STAT3, CD44, ITGAM, CCL2, TIMP1, MYC, THBS1, SOCS3, and CD14, and IGF1 expression was significantly reduced when compared with normal controls (p<0.05) (Figure 5).

Figure 5

Quantitative reverse transcription-polymerase chain reaction (qRT-PCR) for the expression of the hub genes in the renal tissue in the rat model of ischemia-reperfusion injury (IRI) and the controls. Nine hub genes were significantly increased, including STAT3, CD44, ITGAM, CCL2, TIMP1, MYC, THBS1, SOCS3, and CD14, and IGF1 expression was significantly decreased in the renal tissue in the rat model of ischemia-reperfusion injury (IRI) when compared with normal control renal tissue (p<0.05). ** P<0.01, * P<0.05.

Discussion

Ischemia-reperfusion injury (IRI) is a common condition that affects major organs that include the heart and kidney following the restoration of previously reduced or absent circulatory perfusion. IRI following surgery is a factor that can affect patient prognosis and can be associated with patient morbidity and mortality [15]. In the early stage of IRI, lack of blood supply and reduced oxygenation results in early changes of cell necrosis, apoptosis, inflammation, and microvascular endothelial cell damage [16]. When blood flow is restored, the tissue of the organ undergoes secondary injury, which can be more severe and is also potentially irreversible [17]. Renal IRI occurs most commonly in following hemorrhage and shock, and also in renal transplantation when the kidney is harvested, and also during cardiopulmonary bypass surgery. The pathogenesis of renal IRI involves inflammation, apoptosis, and autophagy. Improving the understanding of the biological and molecular mechanisms associated with renal IRI may improve the understanding of IRI, and may also result in predictive, preventive, and therapeutic biomarkers [18].

Previous studies have identified mRNA expression profiles in renal IRI, including for the SPRR2F, SPRR1A, and MMP-10 genes. However, these genes were unlikely to represent biomarkers for renal IRI due to their lack of disease-associated specificity. Previously identified differentially expressed genes (DEGs) are unlikely to represent biomarkers that are specific for renal IRI, as the identified hub genes have shown different profiles in previous studies [19]. This variability is due to the variable characteristics of the patients studied, including gender, age, and ethnic differences. Therefore, integration of the different gene profiles may remove the irrelevant information and focus on the significant DEGs in the different profiles.

The findings from the present study identified 232 DEGs that were co-expressed in the rat model of IRI, including 170 upregulated DEGs and 62 down-regulated DEGs. Using the Kyoto Encyclopedia of Genes and Genomes (KEGG) database, several relevant pathways associated with renal IRI were identified that were consistent with the findings from previous studies [20,21]. The use of protein-protein interaction (PPI) network modular analysis identified the top ten hub genes, which included STAT3, CD44, ITGAM, CCL2, TIMP1, MYC, THBS1, IGF1, SOCS3, and CD14. The HIF-1α signaling pathway was found to have a key role in IRI in the rat model. Also, the use of the Cytoscape plug-in algorithm, Molecular Complex Detection (MCODE), showed the HIF-1α signaling pathway as the most relevant pathway. Therefore, this study included the analysis of the hub gene STAT3 and the HIF-1α signaling pathway.

For the past few decades, the signal transducers and activators of transcription 3 (STAT3) gene have been known to participate in innate immune responses [22]. The STAT family has transcriptional roles and non-transcriptional activities, including cellular respiration, autophagy, and metabolism [22]. Previous studies had reported that STAT3 has a core role in IRI. Yu et al. [23] showed that STAT3 was upregulated in hepatic IRI. Tang et al. [24] reported that STAT3 was expressed in cerebral IRI. A recently reported study that used a mouse model if IRI with a genetic deletion of STAT3 was used to study the role of STAT3 signaling in renal IRI [25]. Increased severity of renal IRI was associated with increased leukocyte adhesion to the endothelium of the renal microvasculature and increased F4/80-positive macrophages [25]. The authors concluded that renal endothelial STAT3 limited renal kidney dysfunction in IRI and that selective therapeutic activation of STAT3 signaling in renal endothelium should be studied further as a potential therapeutic approach in renal IRI [25]. Also, ursolic acid had been reported as a potential therapeutic approach to reduce IRI through modulation of STAT3 activity and reduce tissue damage due to inflammation [26]. In 2018, Han et al. [27] showed that STAT3 was upregulated in hepatic IRI, and was associated with increased autophagy, and proposed that STAT3 had a central role in hepatic IRI by activating autophagy. These previous studies highlight the importance of studies in animal models and support the findings from the present study.

In this study, the use of multiple bioinformatics tools identified key hub genes, which were then verified by quantitative reverse transcription-polymerase chain reaction (qRT-PCR) in the renal tissue of the rat model of renal IRI. The expression of STAT3 was identified in the rat model, which is a finding supported by previous studies [26,27]. With the use of the cytoHubba and MCODE analysis, as well as STAT3, the hypoxia-inducible factor-1α (HIF-1α) gene was also identified in the rat model of IRI. HIF-1α is a transcriptional regulator that participates in hypoxic adaptation. HIF-1α is the main component of three HIF-α family proteins, which have essential roles in the regulation of tumor microenvironment, the inflammatory response, tumor metabolism, and cancer progression and metastasis [28]. Recently, Cheng et al. [29] studied the role of pre-reatment with Angelica sinensis before the onset of cerebral ischemia, and showed that the HIF-1α signaling pathway regulated cerebral IRI and had angiogenic and anti-apoptotic effects. Activation of the HIF-1α signaling pathway has also been shown to have lung and cardioprotective effects against IRI injury [30,31]. Conde et al. [32] performed in vitro and in vivo studies that showed that the HIF-1α signaling pathway participated in renal IRI in renal tissue used for transplantation. These authors also used a rat model of IRl and found that the induction of HIF-1α during IRI had a protective effect following IRI [33].

HIF-1α and STAT3 have functional interactions in biological processes. HIF-1α may have a central downstream role in theSTAT3 pathway to regulate complex signaling networks. Aghazadeh et al. [34] showed that the STAT3/HIF-1α/Hes-1 axis promoted trastuzumab resistance in breast cancer cells in vitro. The relationship between HIF-1α and STAT3 in renal disease had also been previously studied. Dodd et al. [35] showed that STAT3 could function as a transcription factor in renal adenoma cells in vitro to regulate the HIF-1α signaling pathway. However, the function of STAT3/HIF-1α had not been previously studied in renal IRI. Therefore, the findings from the present study, although this was a preliminary study in an animal model, provided additional support that the STAT3/HIF-1α signaling pathway may have a role in renal IRI.

This study had several limitations. In this preliminary study, a rat model was used to study renal IRI. Therefore, the molecular findings should be interpreted with caution and validated in further in vivo models. In particular, studies are required to validate the role of the HIF-1α pathway in renal IRI and the role of the STAT3 hub gene and the HIF-1α pathway.

Conclusions

This study aimed to identify hub genes and pathways in a rat model of renal ischemia-reperfusion injury (IRI) using bioinformatics analysis of the Gene Expression Omnibus (GEO) microarray datasets and the integration of gene expression profiles. In the rat model of renal IRI, bioinformatics analysis and integration of gene expression profiles identified the involvement of HIF-1α signaling and STAT3 as a hub gene, supporting the potential role of the STAT3/HIF-1α signaling pathway in this model.