RNA Sequencing for Gene Expression Profiles in a Rat Model of Middle Cerebral Artery Occlusion

Gene expression regulatory mechanisms in models of middle cerebral artery occlusion (MCAO) have been assessed in some studies, but questions remain. The discovery of differentially expressed genes (DEGs) between MCAO and control rats analyzed by next-generation RNA sequencing is of particular interest. These DEGs may help guide the clinical diagnosis of stroke and facilitate selection of the optimal treatment strategy. Twenty rats were equally divided into control and MCAO groups. Three rats from each group were randomly selected for RNA sequencing analysis. Sequence reads were obtained from an Illumina HiSeq2500 platform and mapped onto the rat reference genome RN6 using Hisat2. We identified 1,007 significantly DEGs with p<0.05, including 785 upregulated (fold change [FC]>2) and 222 downregulated (FC<0.5) DEGs, in brain tissue from MCAO rats compared with that from control rats, and numerous immune response genes were identified. Gene ontology (GO) analysis revealed that the majority of the most enriched and meaningful biological process terms were mainly involved in immune response, inflammatory response, cell activation, leukocyte migration, cell adhesion, response to external stimulus, cell migration, and response to wounding. Also enriched were immune-related pathways and neural-related pathways. Similar to GO molecular function terms, the enriched terms were mainly related to cytokine receptor activity. Six DEGs were verified by quantitative real-time polymerase chain reaction (qRT-PCR). Protein-protein interaction analysis of these hits revealed that MCAO affects complement and coagulation cascades and chemokine signaling pathway. Our study identified novel biomarkers that could help further elucidate MCAO mechanisms and processes.


Introduction
Ischemic stroke is a leading cause of global mortality and morbidity. Immediately following the event, neurovascular reperfusion to the infarcted area can injure the tissue, resulting in more serious disability [1,2]. Middle cerebral artery occlusion (MCAO) is a widely used model for studying neuroprotective therapies. For example, it was used to test the effectiveness of thrombolytic therapy [3,4]. Thrombolysis is one of the most effective and economical treatments for ischemic stroke [5]. It is well known that brain damage can be exacerbated by reperfusion after thrombolysis [6].
Presumably, this secondary injury impacts multiple gene functions. Most existing studies have analyzed the functions of specific genes and signaling pathways. [7][8][9]. We propose that a broad analysis of gene expression changes could be used to develop better neuroprotectants for cerebral ischemia.
In recent decades, mRNA expression profiling performed by microarray or high-throughput RNA sequencing (RNAseq) has been used to uncover molecular mechanisms and explore diagnostic and predictive biomarkers, particularly in complicated diseases such as diabetes, cancer, and stroke [10][11][12]. RNA-seq was applied for global gene expression profiling of the hippocampus following reperfusion with orexin-A [13].  [14].
The main aim of this study was to identify protein-coding genes regulation networks 4 days after MCAO. RNA-seq was used to investigate the mRNA profiles in the brain tissue of rats subjected to MCAO. By a series of bioinformatics analyses, many immune-related pathways (B cell receptor signaling pathway, primary immunodeficiency, Fc epsilon RI signaling pathway, natural killer cell-mediated cytotoxicity, chemokine signaling pathway, cytokine-cytokine receptor interaction, complement and coagulation cascades, etc.) and neuralrelated pathways (neuroactive ligand-receptor interaction, neurotrophin signaling pathway) were identified, and six key genes were verified by quantitative real-time polymerase chain reaction (qRT-PCR).

. . Animal Experiments and Sample Collection. Male
Sprague Dawley rats (180-230 g, SPF grade) were obtained from the Experimental Animal Center of Anhui Medical University by Animal Experiments of Anhui university of Chinese medicine. All rats were anesthetized with sodium pentobarbital during all surgical procedures to minimize suffering, and the MCAO and sham surgeries were performed as previously described [15]. Four days after MCAO, the left hemispheres were collected and immediately frozen in liquid nitrogen.
. . RNA Extraction and Sequencing. Total RNA was extracted using mirVana6 miRNA Isolation kit (Cat# AM1561, Ambion, Foster City, CA) following the manufacturer's instructions. RNA sample quantity and quality were determined using a NanoDrop 2000 (Thermo Fisher Scientific, Waltham, MA) and an Agilent 2100 bioanalyzer (Agilent Technologies, Santa Clara, CA). TruSeq5 Stranded Total RNA Sample Preparation kits (Illumina, San Diego, CA) were used to prepare libraries following the manufacturer's instructions. Purified libraries were quantified by a Qubit5 2.0 Fluorometer (Life Technologies, Carlsbad, CA) and Agilent 2100 bioanalyzer. Clusters were generated by cBot with the library and sequenced on an Illumina HiSeq 2500 platform (San Diego, CA). All sequencing was performed at Origin-Biotech Inc. (Ao-Ji Biotech, Shanghai, China).
. . Functional Enrichment Analysis. GO and KEGG pathways were enriched by R package (v 3.5.1) to better understand the functions of the DEGs [23]. In our study, clusterProfiler was applied to analysis of GO terms and KEGG pathways, and the top 30 GOs and pathways are presented [24].

. . PPI Network Construction and Module Analysis.
STRING is a database that provides comprehensive information about interactions between proteins, including prediction and experimental interaction data [25]. In our study, the STRING tool was used to map PPIs among the DEGs considering interactions of combined scores ≥0.4. Then, Cytoscape was used to visualize the network [26]. The PPI network was used to filter modules based on the Molecular Complex Detection plug-in (MCOD) in Cytoscape [27] with standard set following degree cut-off=2, k-core=2, node score cut-off=0.2, and max depth=100. An MCODE score ≥4 and node ≥10 were considered for functional enrichment analysis of the modules. GO and KEGG enrichment for DEGs in the four modules were performed in clusterProfiler.
. . Validation of Differentially Expressed mRNAs from the qRT-PCR Sequencing Profile. qRT-PCR is the gold standard of mRNA detection and was performed to verify the RNA-seq results. Glyceraldehyde 3-phosphate dehydrogenase (GAPDH) served as the internal control. Relative mRNA expression was determined using the 2-ΔΔCT method. Six genes were analyzed: Top2a, Cdk1, Ccna2, Ccnb1, Th, and Esr1. Brain tissue was sequenced from the MCAO and control groups (n=3/group), and each experiment was performed in triplicate.

Results
. . DEG Screening. EdgeR software was used to screen DEGs with p<0.05 and |log 2 FC|>1. We identified 1,007 DEGs: 785 upregulated and 222 downregulated. Differential mRNA expression between the MCAO and control groups was represented in volcano and scatter plots (Figures 1(a) and 1(b)). Hierarchical clustering of DEGs was visualized (Figure 1(c)). The top 20 up-and downregulated DEGs are listed in Tables 1 and 2. . . GO Functional Enrichment Analysis of DEGs. GO enrichment analysis was performed with 1,007 DEGs in clusterProfiler. The DEGs between the MCAO and control groups were enriched to 54 subclasses of GOs, and the top 30 subclasses are shown in Figure 2(a). The top 10 enriched GO biological processes were immune system process, immune response, cell activation, leukocyte migration, response to external stimulus, cell migration, response to wounding, inflammatory response, cell adhesion, and      . . Pathway Enrichment Analysis of DEGs. Pathway enrichment analysis of DEGs could provide further insight into the function of genes and their interactions. The DEGs between the MCAO and control groups were enriched to 31 subclasses of pathways in 5 broad categories (cellular processes, genetic information processing, organic systems, metabolism, and environmental information processing) after these were analyzed using the KEGG database (Figure 3(a)). We performed KEGG pathway enrichment analysis for DEGs and found 237 pathway terms including 77 that were significant (p< 0.05). The top 10 pathways with the greatest enrichment were Staphylococcus aureus infection (rno05150), osteoclast differentiation (rno04380), complement and coagulation cascades (rno04610), malaria (rno05144), leishmaniasis (rno05140), chemokine signaling pathway (rno04062), Chagas disease (American trypanosomiasis) (rno05142), hematopoietic cell lineage (rno04640), amoebiasis (rno05146), and phagosome (rno04145). The top 30 enrichment pathways are presented in Figure 3(b).
. . PPI Network. Significantly altered DEGs were used to construct a PPI network based on the STRING database. The network comprises 744 nodes and 5,266 interactions ( Figure 4). More than 100 genes of connectivity degrees were >100, and the details of the top 10 genes with connectivity degrees were listed ( A total of 25 modules were obtained using default criteria in MCODE. Modules were listed in descending order by MCODE score. Four modules with MCODE score ≥4 and nodes ≥10 were named modules 1, 2, 3, and 4. These were selected for module network visualization, GO enrichment analysis, and KEGG pathway enrichment analysis (Figures 5-8). Most DEGs in these modules were  . . qRT-PCR Verification of the DEGs. Top2a, Cdk1, Ccna2, and Ccnb1 were overexpressed in the MCAO group (Figure 9). Th and Esr1 had lower expression in the MCAO for both RNA-seq and qRT-PCR. We successfully verified that the RNA-seq and qRT-PCR results were consistent, indicating that the RNA-seq results were reliable. The primers are listed in Table 4.

Discussion
We analyzed the protein-coding mRNA expression profiles in cerebral hemispheres from experimental and control rats by high-throughput RNA-seq. This identified 1,007 DEGs, which were subjected to GO, KEGG pathway, PPI, and PPI module analysis to provide a better understanding of pathological mechanisms that occur after MCAO.  (Cdk1, FC=25.09), Cdk1 phosphorylates key substrates leading the cell to G2-phase arrest, M-phase, and cytokinesis promotion. Some studies reported that Cdk1 activation is involved in multiple neuronal death by activating phosphorylation of Bad27 (a proapoptotic protein) or inhibitory phosphorylation of Bcl-xL, Bcl-2, and Mcl-1, which are antiapoptotic [28,29]. Others described similar results in transient ischemia [28,30]. Besides the Bcl-2 family, Cdk1 also phosphorylates the transcription factor FOXO1, which is also implicated in cell death [31]. We hypothesize that during MCAO, Cdk1 induces ischemic neuronal death by Bcl2 family or Foxo1 promoting signaling pathways.
. . Inflammatory Response. GO enrichment analysis revealed that more than 78 DEGs are involved in the inflammatory response (p=0.0078), including 75 that were upregulated. Several inflammatory response pathways were activated like complement and coagulation cascades, chemokine signaling pathway, natural killer cell-mediated cytotoxicity, and leukocyte transendothelial migration. Dozens of cytokines were dysregulated. The macrophage recruitment factor CCL2 [32][33][34] was upregulated, and CD68 and CD163 were dysregulated. Therefore, we hypothesize  that MCAO could lead to the recruitment and induction of macrophages, natural killer cells, and leukocytes. The specific cell subsets will require further verification.
group. We hypothesize that MCAO could induce Htr2a and Htr2b overexpression leading to excitability in the early stage of ischemia. Other studies have shown suppression of GABAA and glycine receptors in rats with MCAO, and receptor agonists could improve it [36,37]. Our results show that some subunits of GABBA (e.g., Gabrr3, Gabrq, Gabra6, Gabre, and Gabrr1) were dysregulated following MCAO.

Conclusions
By high-throughput RNA-seq, we analyzed the proteincoding mRNA expression profile in control and MCAO groups in brain tissue. And some GOs, pathways, and genes that were identified could play key roles with MCAO. These findings may help us to understand the underlying mechanism of protein-coding mRNAs with MCAO.

Data Availability
The data used to support the findings of this study are available from the corresponding author upon request.

Conflicts of Interest
The authors declare that there are no conflicts of interest regarding the publication of this paper.