RNA-seq expression profiling of rat MCAO model following reperfusion Orexin-A

Orexin-A is a neuropeptide with potent neuroprotective activity towards cerebral ischemia-reperfusion (I/R) injury, but few studies have attempted to elucidate the mechanism. Herein, we performed global gene expression profiling of the hippocampus following reperfusion with Orexin-A using RNA sequencing (RNA-seq). RNA-seq identified 649 differentially expressed genes (DEGs) in the Orexin-A group compared with saline controls (I/R group), of which 149 were up-regulated and 500 were down-regulated. DEGs were confirmed using qRT-PCR, their molecular functions, biological processes and molecular components were explored using Gene Ontology (GO) analysis and 206 KEGG pathways were associated with Orexin-A treatment. MAPK, chemokine and calcium signalling pathways were mainly responsible for the neuroprotective effects of Orexin-A. Hspb1, Igf2 and Ptk2b were selected for functional interaction analysis by GeneMANIA. The results suggest that Orexin-A modifies gene expression in the hippocampus, leading to neuroprotection from I/R injury. The study provides a basis for future elucidation of the molecular mechanisms underlying Orexin-A.

Similarly, Orexin-A has a prominent neuroprotective effect on cerebral I/R injury. Kitamura et al. found that the number of neurons expressing Orexin-A was significantly lower in the non-ischemic side than the ischemic side [16]. Additionally, they found Orexin-A pre-treatment significantly reduced the brain infarction volume induced by reperfusion injury, ameliorated neurologic deficit scores and reduced the infarction volume after rat cerebral I/R injury. The neuroprotective effect of Orexin-A was shown to be mediated by an increase in hypoxia-inducible factor-1 (HIF-1α) activity [17], and our lab demonstrated that Orexin-A protects SHY-5Y cells against oxidative stress by activating the PI3K/AKT signalling pathway.
Few studies have investigated the neuroprotective effects of Orexin-A in I/R injury. Therefore, in the present study we performed global gene expression profiling www.impactjournals.com/oncotarget/ Oncotarget, 2017, Vol. 8, (No. 68), pp: 113066-113081 Research Paper www.impactjournals.com/oncotarget in the hippocampus after reperfusion with Orexin-A by RNA sequencing (RNA-seq). We identified differentially expressed genes (DEGs) and signalling pathways putatively associated with Orexin-A reperfusion. The findings provide a basis for further elucidation of the molecular mechanisms underlying the neuroprotective effects of Orexin-A.

Orexin-A decreases the infarct volume after cerebral I/R injury
Infarction volumes in model rats were measured using triphenyltetrazolium chloride (TTC) staining. As shown in Figure 1, normal tissues were stained red, while infarction areas were stained white. No infarction area was identified in the sham group, while the I/R group displayed a significantly higher rate of infarction (30.620% ± 2.665%) than the sham group. In addition, the Orexin-A group displayed a rate of infarction that was 12.497% ± 1.912% lower than that displayed by the I/R group ( ** p < 0.01), indicating that pre-treatment with Orexin-A significantly decreased the infarction area, and Orexin-A therefore exhibited a clear neuroprotective effect on I/R injury.

Global analysis of RNA-seq data
Clean reads from RNA-seq were filtered, and 14,199,598, 13,499,205 and 13,116,887 total reads were acquired in the sham, I/R and Orexin-A groups, respectively. Analysis of global gene expression is shown in Figure 2. About 89% of reads were successfully mapped, 11% were unmapped, 71% were perfectly matched with reference sequences and 17-18% were mismatched ( Figure 2A). Figure 2B shows the number of common and unique genes in the sham, I/R and Orexin-A groups. A total of 12,457 genes were found to be commonly expressed in all three groups, and 221 genes were specific to the sham sample, 294 were unique to the I/R group and 206 were expressed only in the Orexin-A group.

Screening of DEGs
Genes were considered differentially expressed when FDR ≤0.001 and |log2 ratio| >1. As shown in Figure 3A, 156 genes were up-regulated and 26 genes were down-regulated in the sham group compared with the I/R group. After reperfusion with Orexin-A, 649 genes were differentially expressed compared with the I/R group, of which 149 were up-regulated and 500 were down-regulated. In addition, there were 186 upregulated genes and 353 down-regulated genes in the Orexin-A group compared with the sham group. Figure  3B shows the number of DEGs and the fold change in expression between the Orexin-A and I/R groups. The absolute fold change (log2 ratio) ranged from 1 to 17.13. Figure 3C shows the expression levels of DEGs in the I/R and Orexin-A groups. Nine genes (Glycam1, Hfe2, Hmgn5, Lect1, LOC100911576, Mybl2, Scgb1c1, Tmem72 and Tmem27) not expressed in the I/R group were significantly up-regulated after reperfusion with Orexin-A. Eight genes (Hmgn5b, LOC100362172, Upk1b, Tmprss5, Msln, Pcdha1, Pcdh11x and Hdx) expressed in the I/R group were not expressed at all after reperfusion with Orexin-A. Additionally, 632 genes were detected in both groups. A selection of DEGs between the Orexin-A and I/R groups is listed in Table 1. These DEGs are likely to be involved in the neuroprotective effects of Orexin-A, and were further investigated.

Figure 1: Orexin-A decreases the cerebral infarction volume following middle cerebral artery occlusion (MCAO) in rats. (A)
Representative images from sham and MCAO groups treated with or without Orexin-A after a 24 h reperfusion using TTC staining. Brain sections from the sham group are red, and infarction volumes are not visible. Infarction volumes in the I/R group are white and obviously increased in size compared with the sham group. Reperfusion of Orexin-A significantly decreases the infarction volume compared with the I/R group. (B) Percentage of cerebral infarction volumes to total brain volumes. Columns represent infarction volumes as a percentage of total volumes. Data are expressed as mean ± SD (n = 6). ** p < 0.01 vs. the I/R group.

qRT-PCR validation
qRT-PCR was performed to validate selected DEGs identified by RNA-seq. The relative expression levels of eight DEGs were determined using a LightCycler 480 (LC480). The sequences of primers were listed in Table 2. As shown in Figure 4, the relative expression levels determined by qRT-PCR were in good agreement with the RNA-seq data. Fives genes (Enpp2, Slc24a2, Pk2b, Cacna1e and Hspb1) were up-regulated in the I/R group compared with the sham group ( ** p < 0.01, *** p < 0.001), and all were down-regulated following pre-treatment with Orexin-A compared with the I/R group ( $ p < 0.05, $$ p < 0.01). Three genes (Fibcd1, Wfs1 and Ttr) were down-regulated at 24 h after reperfusion compared with the sham group ( ** p < 0.01), whereas Orexin-A reperfusion clearly increased their expression compared with the I/R group ( $ p < 0.05, $$ p < 0.01). The results confirmed the reliability of the RNA-seq data and further indicated the involvement of the identified DEGs in the neuroprotective effects of Orexin-A.

GO annotation of DEGs
GO annotation was used to further investigate the functions of the identified DEGs. As shown in Figure 5, DEGs were annotated into three categories: molecular function, biological process and cellular component. The molecular function category was divided into eight main subcategories, and the top four were binding (GO:0005488), catalytic activity (GO:0003824), receptor activity (GO:0004872) and transporter activity (GO:0005215; Figure 5A). The top hit in the biological process category was cellular process (GO:0009987), metabolic process (GO:0044238), biological regulation (GO:0065007) and response to stimulus (GO:0050896) were also well represented ( Figure 5B). Cell part (GO:0044464), membrane (GO:0016020), organelle (GO:0043226) and macromolecular complex (GO:0032991) were the main cellular component subcategories ( Figure 5C). A selection of DEGs classified into different subcategories is listed in Table 3.
We chose two subclasses for further exploration using AmiGO 2 ( Figure 6). DEGs related to programmed cell death could be involved in autophagic cell death, cornification, pyroptosis, mitochondria programmed cell death and death receptor activity. DEGs associated with defence responses may play an important role in the cellular defence response, the innate immune response, the inflammatory response, clearance of foreign intracellular displaying differentially expressed genes between the I/R and Orexin-A groups. The x-and y-axes show the fold change (log2 ratio) in gene expression. Yellow points represent up-regulated genes with a fold change less than 1 and p < 0.05; blue points denote down-regulated genes with a fold change greater than 1 and p < 0.05. Brown points represent genes with no significant differences. (C) Venn diagrams showing the number of commonly and uniquely differentially expressed genes between the I/R and Orexin-A groups. Nine genes were expressed only in the I/R group, and eight genes were expressed only after reperfusion with Orexin-A.   Chemokine signalling pathway (ko04062; 15 DEGs), GnRH signalling pathway (ko04912; 13 DEGs) and Neurotrophin signalling pathway (10 DEGs; Figure 7 and Table 4). These signalling pathways likely play important roles in a various manner in the neuroprotective effects of Orexin-A.

Gene interaction network analysis of DEGs
We selected three of the more interesting DEGs in the Orexin-A group compared with the I/R group for gene interaction analysis by GeneMANIA. The functional interaction network of Hspb1 ( Figure 8A) includes seven Hspb1-interacting genes that were found to be differentially expressed by RNA-seq, of which five (Hspa1b, Odf1, Cryaa, Hspb2 and Flinc) were down-regulated and two (Hspb3 and Csrp3) were up-regulated (Table 5). Meanwhile, eight down-regulated and two up-regulated genes (Igf1, Igfbp1, Igfbp4, Gpc3, Padi3, Cyp4a8, Rln3, Insl3, Igfbp2 and Ins1; Table 5) appeared to functionally interact with Igf2 ( Figure 8B). Similarly, six genes interacted with Ptk2b, of which Slc2a1, Jak2, Mcam and Epha7 were down-  Table 5).
In addition, we analysed the signalling pathways involving Hspb1, Igf2 and Ptk2b that may be involved in Orexin-A reperfusion. RNA-seq data showed that Hspb1 was mainly involved in the MAPK and VEGF signalling pathways, whereas Igf2 was most closely associated with the P53 signalling pathway. Ptk2b could potentially affect four signalling pathways: the chemokine, nicotine addiction, calcium and Jak-STAT signalling pathways. Elucidation of these DEGs and their associated signalling pathways could prove critical for understanding the neuroprotective mechanisms of Orexin-A.

DISCUSSION
RNA-seq technology has emerged as a powerful method for screening transcripts and has been applied to a broad range of fields [18,19]. Bermudez et al. used RNAseq to screen DEGs and corresponding pathways involved in glucocorticoid responsiveness in a bovine trabecular meshwork (BTM) cell strain. Sequencing results showed that 93 and 606 genes were differentially expressed in responder and non-responder BTM cells, respectively [20]. Wang et al. analysed stage-specific genes associated with gastric cancer using RNA-seq and identified 2,224 genes specifically expressed in stages I and II, and 539   genes specific to stage IV [21]. Meanwhile, Holmes et al. assessed the impact of three different intensities of transcranial direct current stimulation (tDCS) on gene expression in rat cerebral cortex by RNA-seq. Compared with the sham group, ~1,000 DEGs were identified at each treatment intensity [22]. In the current study, Orexin-A alleviated the cerebral infarction volume induced by I/R injury, consistent with the neuroprotective effects of Orexin-A on I/R injury. To more comprehensively reveal which genes and pathways are involved in the neuroprotective effects of Orexin-A, RNA-seq was carried out on sham, I/R and Orexin-A rat groups. To our knowledge, this is the first study to use RNA-seq to perform global gene expression profiling after reperfusion with Orexin-A in the hippocampus. Compared with the I/R group, ≥2-fold up-regulation of 149 genes and ≥2-fold down-regulation of 500 genes was observed following reperfusion with Orexin-A, identifying 649 genes that are potentially involved in the neuroprotective effects of Orexin-A. Go analysis based on the RNA-seq data identified DEGs involved in adhesion, metabolism, localisation, immune reaction, programmed cell death and defence responses (Figures 5 and 6; Table 3). We screened several immune response-related genes including IL11 (-2.359), Cd55/44 (-1.681/-3.4947), Ccl2 (-4.059), Ifit1 (-4.2288) and Cadm1 (-1.079) differentially expressed after reperfusion with Orexin-A. A number of genes related to cell death, including Bcl3 (-4.9307), caspase 3/4 (-0.826/-1.838), Rassf5 (-1.2599), Myc (-1.0524), Casp8ap2 (-1.307) and Bag4 (-1.0895), were down-regulated in the Orexin-A group. In addition, reperfusion induced the up-regulation of Hif1a, Stat3, Hspb2, Hspa1a, Hspa1b and Hspb1, which play an important role in stimulus or stress, but Orexin-A lowered their expression levels (log2 ratio) by - 1.1357-, -1.135-, -9.951-, -3.474-, -1.51-and -2.602fold, respectively. These differences in expression after reperfusion with Orexin-A may contribute to its neuroprotective effects. Orexin-A has been shown to activate multiple protein kinases such as PKA, PKC, MAPK/Erk and PDK1 in various cellular contexts [23]. Orexin-A activates the mTORC1 pathway via extracellular calcium influx and the lysosome pathway in recombinant cell lines, leading to the central regulation of cell growth and metabolism [24]. Orexin-A increases glucose transporter 4 expression and lipid accumulation in 3T3-L1 adipocytes via ERK1/2, JNK and p38 MAPK signalling [25,26]. Orexin-A treatment promoted the uptake of Glu via the OX1R/ PKCα/ERK1/2 pathway and protected co-cultured cells against anoxia/hypoglycaemic injury in astrocytes [27]. Orexin-A stimulates ERK 1/2 phosphorylation and facilitates cell migration via the PLC-PKCα signalling pathway in cultured rat astrocytes [28]. In the current study, KEGG pathway analysis identified 206 pathways that were closely associated with Orexin-A reperfusion, of which metabolic pathways accounted for 34 DEGs, and MAPK signalling accounted for 21 DEGs. It was not surprising to find that metabolic pathways and MAPK signalling pathways may play important roles in Orexin-A reperfusion.
Besides these two pathways, several others were also detected, including Wnt, calcium, chemokine, gonadotropin-releasing hormone (GnRH) and neurotrophin signalling pathways. Orexins play important roles in nociceptive modulation via activating intracellular calcium signalling [29]. Orexin-A and Orexin-B significantly increase the expression of neurotrophin-3 (NT-3) at the mRNA level when primary cortical neurons are cultured with 0.01, 0.1 or 1 mM Orexins [30], and Orexin-A directly stimulates GnRH neurons to secrete GnRH, implicating Orexin-A in GnRH signalling [31].
Taken together, these results implied that Orexin-A might regulate the above signalling pathways, and thereby modulate cerebral I/R injury.
Among the identified DEGs, three genes of particular interest were selected for analysis using GeneMANIA. Heat shock protein beta-1 (Hspb1), also known as Hsp27, is a member of the small HSP family and a key mediator in the cellular response to environmental and physiological stresses [32][33][34]. Upon expose to oxidative stress, HSP27 functions as an antioxidant to decrease the levels of reactive oxygen species (ROS) [35]. In benign prostatic hyperplasia (BPH), the expression of HSP27 is gradually increased and accompanied by inflammation [36]. In retinal capillary endothelial cells, HSP27 functions as an anti-apoptotic protein to preserve the endothelial barrier [37]. Others have similarly noted that HSP27 inhibits cell death via a c-jun-mediated cascade in cerebral ischemia [38]. HSP27 transgenic mice display amelioration of the infarction volume, behavioural deficits and blood brain barrier (BBB) damage in a stroke model [39]. Insulin-like growth factor 2 (Igf2) is indispensable in a number of cells because it promotes cellular proliferation and differentiation [40]. In medulloblastoma, Igf2 promotes the proliferation of precursor cells and the PZp53 cell line [41]. IGF2 also mediates the proliferation of hippocampal neural stem cells via IGF1R and AKTdependent signalling [42]. IGF2 is highly expressed in microglia and exhibits a strong protective effect against cytokine-mediated neuronal death in vitro and in vivo [43]. Protein tyrosine kinase type 2 beta (Ptk2b) was found to play a vital role in the response to diverse cellular stimuli and diseases. Ptk2b specifically induces cell death and reduces clonogenic growth in myeloma cell lines and multiple myeloma stem cells. Furthermore, Ptk2b inhibition attenuates multiple myeloma progression [44]. Ptk2b can be up-regulated by lovastatin, resulting in an anti-apoptotic effect on ox-LDL-induced cell injury in endothelial cells. The protective effect of Ptk2b is associated with the AKT signalling pathway [45]. At 24-72 h after focal cerebral ischemia, phospho-Ptk2b was rapidly induced in microglia surrounding the necrotic infarction area, and phospho-Ptk2b was confirmed to colocalise with phospho-p38 and act as an upstream mediator of the p38 signalling pathway after cerebral ischemia [46]. Consistent with these previous findings, these three genes likely play a vital role in Orexin-A-mediated neuroprotection, and their functions and molecular mechanisms will be investigated in future studies. It has found that orexin-A has multiple effects on different tissues. By comparing the mechanism of orexin-A in different tissues, it can be found that there is commonality and tissue specificity. It will provide a comprehensive understanding of Orexin-A effects on different organs and reveals whether its neuroprotective effects are tissue specific at molecular levels.
In conclusion, our global expression profiling identified genes, pathways, functions and interaction networks potentially involved in the neuroprotective effects of Orexin-A reperfusion. However, the results obtained by RNA-seq and bioinformatic analysis were general and indicate the involvement of numerous complex processes, and further experiments are needed to validate the identified DEGs and signalling pathways to clarify the molecular mechanisms involved in the process of Orexin-A reperfusion.

Animals
Male Wistar rats, aged 6-8 weeks, were obtained from Pengyue Experimental Animal Breeding Institute of Jinan, China (License No. SCXK (Lu) 2014-0007). All rats were maintained in a room at 25 ± 3°C and 50-60% humidity, with a 12 h light-dark cycle and free access to food and water. All experiments were approved by the Ethics Committee of Jining Medical University and performed in accordance with the National Health Guide for the Care and Use of Laboratory Animals.

Animal grouping
A total of 18 rats were randomly divided into three groups: (1) sham-treated group (sham), from which the right carotid artery was isolated without further processing; (2) I/R group, in which rats were subjected to a 2 h right middle cerebral artery occlusion (MCAO) followed by a 24 h reperfusion with saline (10 μl, 0.9% NaCl) after ischemia at the onset of reperfusion; (3) Orexin-A reperfusion group (Orexin-A), in which rats were injected with Orexin-A (50 ng/kg; Phoenix Pharmaceuticals, Inc., Burlingame, CA, USA) after reperfusion.

MCAO model
Rat MCAO model experiments were performed using the bolt line method. Firstly, rats were anesthetised by intraperitoneal injection of 10% chloral hydrate (0.3 mL/100 g) and swabbed with iodophor. The common carotid artery (CCA) was isolated, and a midline incision was made to expose the internal carotid artery (ICA), external carotid artery (ECA) and right CCA. A nylon suture (0.26 mm diameter) was inserted into the right CCA lumen and gently injected into the ICA at a depth of 18 mm, and the neck incision was sutured. Finally, nylon sutures were gently removed at 2 h after ischemia, and reperfusion was performed for 24 h.

TTC staining
The cerebral infarction volume was measured using TTC (Sigma, St Louis, MO, USA) staining. At 24 h after reperfusion, brains were removed, frozen at -20°C for 30 min, cut into 2 mm thick coronal slices and stained with 0.2% TTC solution in the dark for 30 min at 37°C. After fixing with 4% formalin, slices were photographed using a IX71 microscope (Olympus, Japan). Infarction zones were analysed with Image-Pro Plus v 6.0 software (Media Cybernetics, Bethesda, MD, USA). The infarction vs. normal volume percentage was calculated.

Construction of cDNA library and RNA-seq
Total RNA was extracted using TRIzol reagent (Invitrogen, Carlsbad, CA, USA) according to the manufacturer's instructions. RNA was quantified and the integrity was determined using a BioSpec-nano spectrophotometer (Shimadzu, Japan). Total RNA was treated with DNase I to remove DNA contamination, and RNA was used to generate a cDNA library. The mRNAs were multiplexed and tagged with standard Illumina tags, and fragments were amplified by PCR. Library quality and quantity were determined using an Agilent 2100 Bioanalyzer and an ABI StepOnePlus Real-Time PCR System. Library products were sequenced using an Illumina HiSeq TM 2000 platform by Shenzhen Huada Genetics Corporation (Shenzhen, China).

qRT-PCR
The same RNA samples were reservetranscribed into cDNA using a cDNA synthesis kit (Tiangen Biotech, Beijing, China) according to the manufacturer's protocol. qRT-PCR was performed in 96-well plates using a Roche LC480. Reactions (20 μL) were mixed with 10 μL of 2×SYBR SuperReal PreMix (Tiangen Biotech, Beijing, China), 2 μL of cDNA, 0.6 μL of forward primer, 0.6 μL of reverse primer and 6.8 μL of ddH 2 O. The amplification procedure was as follows: 95°C for 10 min, followed by 40 cycles of 95°C for 15 s and 60°C for 1 min. PCR primers are listed in Table 1. β-actin was amplified as an endogenous control. The relative mRNA level of each gene was calculated using the 2 -ΔΔCT method. At least three replicates were performed for each cDNA sample, and reactions were prepared in triplicate.

Statistical analysis
All data are expressed as mean ± standard error of the mean (SEM) and were analysed using GraphPad Prism one-way analysis of variance (ANOVA). A p-value less than 0.05 was considered statistically significant.