Abstract

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.

1. 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. [79]. 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 (RNA-seq) has been used to uncover molecular mechanisms and explore diagnostic and predictive biomarkers, particularly in complicated diseases such as diabetes, cancer, and stroke [1012]. RNA-seq was applied for global gene expression profiling of the hippocampus following reperfusion with orexin-A [13]. Gaowa et al. performed RNA-seq to investigate transcriptome expression in a rat stroke model treated with the traditional Mongolian medicine Eerdun Wurile [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 neural-related 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).

2. Material and Methods

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

2.2. RNA Extraction and Sequencing

Total RNA was extracted using mirVana™ 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). TruSeq® 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 Qubit® 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).

2.3. DEG Analysis

FastQC was conducted for Quality control (QC) of RNA-seq reads (v. 0.11.3) (http://www.bioinformatics.babraham.ac.uk/projects/fastqc). Trimming was performed by seqtk for known Illumina TruSeq adapter sequences, poor reads, and ribosome RNA reads (https://github.com/lh3/seqtk). The trimmed reads were then mapped to the Rattus norvegicus reference genome (rno6) by the Hisat2 (version: 2.0.4) [16, 17]. StringTie (version: 1.3.0) was performed for each gene count from trimmed reads [17, 18]. Gene counts were normalized by trimmed mean of M-values (TMM) [19], and fragments per kilobase of transcript per million mapped reads (FPKM) in Perl script [20]. edgeR was performed for determining differential expression genes [21] and threshold with p<0.05 and absolute values of log2 (fold change) >1 [22].

2.4. 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].

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

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

3. Results

3.1. DEG Screening

EdgeR software was used to screen DEGs with p<0.05 and |log2FC|>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.

3.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 regulation of immune system process. The top 10 enriched cellular components were extracellular space, chromosome, cell surface, extracellular region, nucleosome, vesicle, external side of plasma membrane, kinetochore, condensed chromosome kinetochore, and protein-DNA complex. The enriched GO molecular functions were protein binding (GO:0005515), protein complex binding (GO:0032403), protein heterodimerization activity (GO:0046982), oxygen transporter activity (GO:0005344), lipid binding (GO:0008289), protein dimerization activity (GO:0046983), oxygen binding (GO:0019825), glycosaminoglycan binding (GO:0005539), cytokine receptor activity (GO:0004896), and iron ion binding (GO:0005506). These data are shown in Figure 2.

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

3.4. 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 (Table 3), including Top2a (topoisomerase [DNA] II alpha, degree=121, FC=38.80), Cdk1 (cyclin-dependent kinase 1, degree=104, FC=25.09), Ccna2 (cyclin A2, degree=98, FC= 28.73), Ccnb1 (G2/mitotic-specific cyclin-B1, degree=96, FC=8.072), Th (tyrosine hydroxylase, degree=46, FC=0.03147), and Esr1 (estrogen receptor 1, degree=18, FC=0.2192).

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 58). Most DEGs in these modules were upregulated in the MCAO group, especially in modules 1 and 2.

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

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.

4.1. Cell Cycle

Pathway enrichment analysis identified 21 DEGs enriched for the cell cycle (P value=1.25E-05), namely, Mcm6, Ttk, Bub1b, Cdkn2c, Ccnb1, Bub1, Plk1, Orc1, Tgfb1, Mcm3, Cdc6, Mcm5, Pttg1, Cdc20, Chek2, Ccnd1, Mcm2, Espl1, Rbl1, Cdk1, and Ccna2. Most were upregulated in MCAO. Following Cyclin A (Ccna2, FC= 28.73) or B (Ccnb1, FC=8.07) binding and activating cyclin-dependent kinase 1 (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.

4.2. 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 [3234] 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.

4.3. Neuroactive Ligand-Receptor Interaction

Pathway enrichment analysis showed that 27 DEGs were involved in the neuroactive ligand-receptor interaction (p=0.00832); 10 and 17 were up- and downregulated, respectively. Upregulated receptors in MCAO were Htr2b, Gabrr1, C5ar1, Cnr2, P2ry6, Htr2a, Galr1, C3ar1, Glra1, and Mc3r. Downregulated receptors in MCAO were Gabrr3, Gabrq, Brs3, Gabra6, Cckar, Trhr, Ptafr, Glra2, Chrnb4, Gpr50, Gabre, Tspo, P2rx2, Ntsr1, Chrna6, Calcr, and Drd3. Lan et al. found decreased serotonin receptor expression in MCAO rats, which was improved with treadmill exercise [35]. We found that Htr2a and Htr2b were upregulated in MCAO 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.

5. Conclusions

By high-throughput RNA-seq, we analyzed the protein-coding 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.

Acknowledgments

The authors gratefully acknowledge financial support from the National Key Research and Development Plan (grant nos. 2017YFC1701600 and 2017YFC1701601), National Natural Science Foundation of China (grant no. 81473387), Anhui Provincial Natural Science Foundation of China (grant no. 1508085QH191), and Key Project of the National Science Fund of Anhui Province (KJ2017A282 and KJ2017A284). We thank Mr. Qiang Fan (Ao-Ji Bio-tech Co., Ltd., Shanghai, China) for assisting with the data analysis.