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/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-1α 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-1α 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.

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 -DDCt method was used, and the PCR reactions were normalized to the GAPDH gene.

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.
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. 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.
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).

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

8406
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-1a 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-1a signaling pathway as the most relevant pathway. Therefore, this study included the analysis of the hub gene STAT3 and the HIF-1a 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-1a (HIF-1a) gene was also identified in the rat model of IRI. HIF-1a is a transcriptional regulator that participates in hypoxic adaptation. HIF-1a is the main component of three HIF-a 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-1a signaling pathway regulated cerebral IRI and had angiogenic and  anti-apoptotic effects. Activation of the HIF-1a 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-1a 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-1a during IRI had a protective effect following IRI [33].
HIF-1a and STAT3 have functional interactions in biological processes. HIF-1a may have a central downstream role in theSTAT3 pathway to regulate complex signaling networks. Aghazadeh et al. [34] showed that the STAT3/HIF-1a/Hes-1 axis promoted trastuzumab resistance in breast cancer cells in vitro. The relationship between HIF-1a 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-1a signaling pathway. However, the function of STAT3/HIF-1a 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-1a 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-1a pathway in renal IRI and the role of the STAT3 hub gene and the HIF-1a pathway.