Weighted gene co-expression network analysis to investigate the key genes implicated in global brain ischemia/reperfusion injury in rats

Material and methods. GSE82146 was extracted from the Gene Expression Omnibus, consisting of 15 complete global brain ischemia (CGBI) reperfusion hippocampus samples and 12 non-ischemic control (NIC) hippocampus samples. The differentially expressed genes (DEGs) between the CGBI and NIC samples were selected using LIMMA package, and were then analyzed with weighted gene co-expression network analysis (WGCNA). Using DAVID software, the DEGs in significant modules were run through enrichment analysis. The DEGs in significant modules were merged, and then a protein–protein interaction (PPI) network was built for them using Cytoscape software. After miRNAs and transcription factors (TFs) were predicted for the DEGs using the WebGestalt tool, a TF-miRNA-target regulatory network was built using Cytoscape software. Furthermore, quantitative real-time polymerase chain reaction (qRT-PCR) analysis was conducted to detect the levels of key genes.


Introduction
Ischemia/reperfusion (I/R) refers to a situation where blood is perfused into tissues experiencing ischemia or hypoxia. 1 Although I/R promotes the repair of damage and the recovery of functioning in most cases, it can also lead to an inflammatory response and oxidative injury by inducing oxidative stress. 2 Ischemia/reperfusion is usually related to microvascular injury, and the imbalance of reactive oxygen species (ROS) and nitric oxide (NO) produced by activated endothelial cells is responsible for the subsequent inflammatory response. 3,4 The development of I/R injury is influenced by ischemia time, aerobic degree, collateral circulation, and reperfusion conditions. The I/R injury has a powerful influence on the ischemic cascade of the brain, involving brain trauma and stroke. 5 Hence, the molecular mechanisms of I/R injury need to be investigated to better alleviate its adverse effects in clinical practice.
By inhibiting nuclear factor-κB (NF-κB), ginkgolide B (GB) possesses anti-apoptotic and anti-inflammatory effects and has demonstrated neuroprotective roles in mice with ischemia-induced brain injury. 6,7 The inhibition of P2X7 receptors (P2X7Rs) protects rats from cerebral I/R injury by decreasing inflammatory response and may serve as a novel therapeutic approach for transient global cerebral I/R injury. By increasing B-cell leukemia-2 (Bcl-2) expression and reducing Bcl-2-associated X protein (Bax) expression, propofol functions as a neuroprotective agent in I/R rats. 8 Oxymatrine can protect the brain of stroke rats from focal I/R injury, and the activation of the nuclear factor erythroid 2-related factor 2 (Nrf2)/hemeoxygenase-1 (HO-1) pathway may promote the neuroprotective effects of oxymatrine in the focal brain I/R rat model. 9 MiR-124 mediates the expression of Ku autoantigen 70 (Ku70) and helps to reduce the neuronal death and brain dysfunction caused by I/R. MiR-134 downregulation relieves cerebral ischemic injury through regulating the enhancing cyclic AMP (cAMP) response element-binding protein (CREB) and downstream genes, which provides a potential therapeutic target for the injury. 10 Nevertheless, the genes and miRNAs affecting brain I/R injury have not been thoroughly explored.
In 2016, Wang et al. performed differential expression analysis on I/R in hippocampus CA1 and CA3, and found that CA3 is better at handling ischemic stress. However, the pathogenesis of I/R injury was not comprehensively researched by them. 11 To further identify the key genes and miRNAs involved in I/R injury, within this study, a series of bioinformatics analyses was carried out (such as differential expression analysis, weighted gene co-expression network analysis (WGCNA), enrichment analysis, and network analysis) on the expression profile data uploaded by Wang et al. 11 In addition, quantitative real-time polymerase chain reaction (qRT-PCR) analysis was conducted to confirm the key genes.

Data source
The microarray dataset GSE82146 (species: Rattus norvegicus) from the Gene Expression Omnibus database (GEO, http://www.ncbi.nih.gov/geo) was extracted, which was determined on the platform of GPL17117 ( RaGene-2_0-st) Affymetrix Rat Gene 2.0 ST Array (transcript (gene) version). There were 15 complete global brain ischemia (CGBI) reperfusion hippocampus samples and 12 non-ischemic control (NIC) hippocampus samples in GSE82146. In Long Evans rats (male, 275-300 g), CGBI was induced with the two-vessel bilateral carotid artery occlusion and hypovolemic hypotension model, 12 as previously described. [13][14][15] Our research was approved by the ethics committee of Ningbo No. 9 Hospital, China.

Differential expression analysis
The original data in CEL format was downloaded and preprocessed (including format conversion, filling in of the missing data, background correction (MicroArray Suite method), and data standardization (quartile method)) using the R oligo package (v. 1.36.1; http://www.bioconductor.org/packages/release/bioc/html/oligo.html). 16 Next, the corresponding genes of the probes were annotated based on the annotation platform. For genes with multiple expression values (mapped to multiple probes), the average value was calculated as the unique expression value.

WGCNA to identify disease-associated modules and genes
The WGCNA is a typical algorithm in system biology for constructing a gene co-expression network which can be used to identify modules of the relevant genes. 18 To screen the CGBI-associated modules and genes, the expression values of the DEGs in each group were determined with WGCNA. 18 The detailed processes of network building and module identification included consistency analysis between the datasets, the definition of gene co-expressioncorrelated matrix (the correlation coefficient between gene m and gene n was S mn = |cor (m,n) |), the definition of the adjacent function (the adjacent function was a mn = power (S mn,β ) ), the determination of the parameters for the adjacent function (the weighting coefficient β ≥ 0.8), the measurement of the degree of dissimilarity between the nodes, the identification of gene modules (the number of genes in the module ≥30), and the determination of the correlation between network module and disease state.

Enrichment analysis
Using DAVID software (v. 6.8; https://david-d.ncifcrf. gov/), 19 Gene Ontology (GO) 20 and Kyoto Encyclopedia of Genes and Genomes (KEGG), 21 enrichment analyses were conducted for the DEGs in significant modules. The thresholds for selecting significant terms were p-values <0.05 and counts of involved genes ≥2.

Protein-protein interaction (PPI) network analysis
After the DEGs in significant modules were merged, they were input into the STRING database (v. 10.0; http:// www.string-db.org/) 22 to predict the PPIs among them. The species was rat and the parameter PPI score was set at 0.4. The PPI results were downloaded in TSV format and Cytoscape (v. 3.2.0; http://www.cytoscape.org/) 23 was used to construct the PPI network.

Transcription factor-miRNA-target regulatory network analysis
Using the WebGestalt tool (http://www.webgestalt. org/), 24 miRNAs and transcription factors (TFs) were predicted for the DEGs involved in the PPI network using Overrepresentation Enrichment Analysis (ORA). The species was rat and the reference background was the Affymetrix Rat Gene v. 2.0 ST platform. According to significance levels, the top 10 results of miRNA-target and TF-target were obtained and integrated. Then, a TF-miRNA-target regulatory network was built using Cytoscape software. 23

qRT-PCR analysis
The total RNA of 7 brain I/R tissues and 7 control tissues were isolated using a Trizol total RNA extraction kit (Invitrogen, Shanghai, China) following the manufacturer's instructions. The purity and integrity of RNA were evaluated separately by spectrophotometer (Merinton, Beijing, China) and 2% agarose gel electrophoresis. The primer sequences of key genes were designed for qRT-PCR experiments (Table 1), and were then produced by Sangon Biotech Co., Ltd. (Shanghai, China). The qRT-PCR experiments were carried out using SYBR Green master mix kit (Applied Biosystems, Foster City, USA). The 20-µL PCR amplification system consisted of 10 µL of SYBR Premix Ex Taq (×2), 8 µL of cDNA template (keeping a consistent level after being diluted with ddH 2 O), 1 µL of forward primer (10 µM), and 1 µL of reverse primer (10 µM). The reaction conditions were 40 cycles of 50°C for 3 min, 95°C for 3 min, 95°C for 10 s, and 60°C for 30 s. Afterwards, a melt curve was created. All experiments were repeated 3 times, and glyceraldehyde-3-phosphate dehydrogenase (GAPDH) was utilized as the reference gene.

Statistical analysis
The expression levels of the key genes were analyzed using the 2 −ΔΔCt method. 25 All data are presented as mean ± standard error of the mean (SEM). SPSS v. 22.0 software (SPSS Inc., Chicago, USA) was used to perform statistical analysis, with p < 0.05 serving as the threshold of statistical significance.

Differential expression analysis
There were 390 DEGs in the CGBI samples compared with the NIC samples, including 330 upregulated genes (such as heat shock protein B1 (HSPB1) and heme oxygenase (decycling) 1 (HMOX1)), and 60 downregulated genes (such as nuclear receptor subfamily 4 group A member 2 (NR4A2)). The clustering heatmap suggested that the samples could be obviously differentiated by the DEGs (Fig. 1).

WGCNA analysis and enrichment analysis
In order to meet the prerequisite of scale-free network distribution, the adjacency matrix weighting parameter "power" was explored. As a result, a "power" value of 19 was selected when the square of correlation coefficient first reached 0.8 (Fig. 2).
A co-expression network was constructed based on this "power" value. Firstly, the dissimilarity coefficients among the DEGs were calculated. Using the dissimilarity matrix, hierarchical clustering was performed to obtain a system clustering tree of the DEGs (Fig. 3A). According to the standards of a hybrid dynamic shear tree, the lowest gene number of each network module was set at 30. After network modules were identified using the dynamic shear method, the feature vector "eigengenes" was calculated for each module. Subsequently, the modules underwent clustering analysis, and the modules with close clustering relationships were merged into a new module. Finally, the DEGs were divided into 5 network modules (Fig. 3B), and the grey module was a set of the DEGs that could not be gathered into other modules.
To identify CGBI-associated modules, the feature vector of each module and CGBI was conducted with correlation analysis. The absolute values of correlation coefficients for color key row Z-score height the modules were sorted, and it was found that the values for brown and turquoise modules were higher than 0.8 (Table 2). Meanwhile, the absolute values of gene significance for the modules were also calculated in order to screen CGBI-associated modules (Fig. 4). In addition, the heatmaps for the genes in the brown and turquoise modules are shown separately in Fig. 5A and 5B.

Functional and pathway enrichment analysis
The DEGs in the brown and turquoise modules were conducted with enrichment analysis. The DEGs in the brown module were mainly implicated in the positive regulation of gene expression (GO; p-value = 9.61E-09) and protein processing in the endoplasmic reticulum (KEGG; p-value = 2.86E-07). Also, the DEGs in the turquoise module were mainly involved in protein folding (GO; p-value = 3.33E-07) and aminoacyl-tRNA biosynthesis (KEGG; p-value = 4.58E-07) ( Table 3).

PPI network analysis
The DEGs in the brown and turquoise modules were merged, and a total of 259 DEGs were obtained (including 115 upregulated genes in the brown module, 115 upregulated genes in the turquoise module, 18 downregulated genes in the brown module, and 11 downregulated genes in the turquoise module). Next, the PPI network was built; it had 164 nodes and 620 edges (Fig. 6A). According to the degree values of the network nodes, heat shock protein 90 alpha family class A member 1 (HSP90AA1; degree = 47) and heat shock protein 5 (HSPA5; degree = 25) were key nodes. Additionally, HSP90AA1 and HSPA5 interacted in the PPI network.

TF-miRNA-target regulatory network analysis
A total of 228 TF-miRNA-target regulatory relationships were obtained, involving 6 TFs, 10 miRNAs, and 99 targets (including 46 upregulated genes in the brown module, 45 upregulated genes in the turquoise module, 5 downregulated genes in the brown module, and 3 downregulated genes in the turquoise module). The TF-miRNA-target regulatory network is shown in Fig. 6B. In the regulatory network, the myelocytomatosis oncogene (MYC; TF, degree = 38), heat shock transcription factor 1 (HSF1; TF, degree = 25), and miR-22 (degree = 12) had higher degree values. Importantly, MYC could target HSPA5 in the regulatory network.

qRT-PCR analysis
The levels of HSPB1, HMOX1 and NR4A2 in the brain I/R tissues and the control tissues were detected using qRT-PCR experiments. HSPB1 (p < 0.001; Fig. 7A) and HMOX1 (p < 0.001; Fig. 7B) were significantly upregulated in the brain I/R tissues compared with the control tissues, while NR4A2 (p < 0.001; Fig. 7C) was significantly downregulated. These findings were consistent with the results of differential expression analysis.

Discussion
In this study, 390 DEGs (including upregulated HSPB1 and HMOX1, as well as downregulated NR4A2) were identified in the CGBI samples. Through WGCNA analysis, the brown and turquoise modules were screened as CGBIassociated modules. After the DEGs in the brown and turquoise modules were merged, a PPI network was built for them. In the PPI network, HSP90AA1 and HSPA5 were the key nodes. Moreover, MYC, HSF1 and miR-22 had higher degree values in the TF-miRNA-target regulatory network. Additionally, the qRT-PCR experiments confirmed upregulated HSPB1 and HMOX1 and downregulated NR4A2.
Oxidative stress can induce the phosphorylation of HSPB1 and HSPB5, which play neuroprotective roles   Red triangles and green hexagons represent miRNAs and TFs, respectively. Brown and blue represent genes in the brown and turquoise modules, respectively. The higher the degree value of a node, the larger the size of the node is. Arrows indicate the directions of regulation in hippocampal neurons. Kupffer cells, which are the main expression sites of hepatic HMOX1, have anti-inflammatory effects and can resist the oxidative injury induced by I/R. 26 NURR1 (also named NR4A2) contributes to intestinal regeneration following I/R injury by suppressing p21 expression, which may provide new approaches for the therapy of intestinal I/R injury. 27 These findings support the thesis that HSPB1, HMOX1 and NR4A2 are related to the development and progression of I/R injury. The mRNA expression of HSP90AA1 is reduced following I/R and may be promoted by miR-1 inhibition during myocardial I/R. 28,29 A high protein expression of HSPA5 can exert neuroprotective effects and stop neural ischemic injury by attenuating endoplasmic reticulum (ER) stressinduced apoptosis. 30 By negatively mediating ubiquitin carboxyl-terminal hydrolase isozyme L1 (UCHL1) and HSPA5 protein levels, miR-181b downregulation protects mice from ischemic injury and provides a therapeutic strategy for ischemic stroke. 31 HSP90AA1 and HSPA5 interacted in the PPI network, suggesting that HSP90AA1 and HSPA5 might play roles in I/R injury through interaction with each other.
The MYC expression is upregulated after acute I/R injury and may promote the low expression of the anti-apoptotic N-myc downstream-regulated gene 2 (NDRG2), which may be associated with myocardial apoptosis in I/R rats. 32,33 By weakening NF-κB activation and reducing MYC expression, copper/zinc-superoxide dismutase (SOD1) overexpression helps to decrease ischemic damage. 34 Granulocyte colony-stimulating factor (G-CSF) increases HSF1 expression by promoting phosphorylation and the interaction of the signal transducer and activator of transcription-3 (STAT3) with HSF1, which possesses cardio-protective effects in I/R mice. 35 HSF1 prevents the death of cardiomyocytes following I/R partly by activating Akt and inactivating caspase 3 and Jun N-terminal kinase. 36 These reports declared that MYC and HSF1 might also be implicated in the mechanisms of I/R injury. MYC could target HSPA5 in the regulatory network, indicating a role of MYC in I/R injury through mediation of HSPA5.
The miR-22 plays a neuroprotective role by reducing inflammation and apoptosis, indicating that miR-22 can be applied in the treatment of cerebral I/R injury. 37 It can suppress the apoptosis of cardiomyocytes by targeting CREB binding protein (CBP); therefore, miR-22 may serve as a novel target for preventing myocardial I/R injury. 38,39 The miR-22 inhibition helps to keep cardiac mitochondrial function, and thus has therapeutic potential for acute myocardial I/R injury. 40 miR-22 decreases caveolin 3 (Cav3) expression and repairs endothelial nitric oxide synthase (eNOS) activity and NO production, inhibiting cardiac injury after I/R. 41 Therefore, miR-22 might be associated with the pathogenesis of I/R injury by regulating the DEGs.

Conclusions
In conclusion, 390 DEGs were identified between CGBI and NIC samples. Also, HSPB1, HMOX1 and NR4A2 were the key genes associated with I/R injury. Moreover, HSP90AA1, HSPA5, MYC, HSF1, and miR-22 might be implicated in the pathogenesis of I/R injury. However, the experimental study was insufficient and our results must still be further confirmed in subsequent studies.