Comparison of Brain Transcriptome of the Greater Horseshoe Bats (Rhinolophus ferrumequinum) in Active and Torpid Episodes

Hibernation is an energy-saving strategy which is widely adopted by heterothermic mammals to survive in the harsh environment. The greater horseshoe bat (Rhinolophus ferrumequinum) can hibernate for a long period in the hibernation season. However, the global gene expression changes between hibernation and non-hibernation season in the greater horseshoe bat remain largely unknown. We herein reported a comprehensive survey of differential gene expression in the brain between winter hibernating and summer active greater horseshoe bats using next-generation sequencing technology. A total of 90,314,174 reads were generated and we identified 1,573 differentially expressed genes between active and torpid states. Interestingly, we found that differentially expressed genes are over-represented in some GO categories (such as metabolic suppression, cellular stress responses and oxidative stress), which suggests neuroprotective strategies might play an important role in hibernation control mechanisms. Our results determined to what extent the brain tissue of the greater horseshoe bats differ in gene expression between summer active and winter hibernating states and provided comprehensive insights into the adaptive mechanisms of bat hibernation.


Introduction
Hibernation is a survival strategy for some mammals (e.g., bats, bears, ground squirrels) against the food deprivation in winter [1,2]. The body temperature, metabolic rate, heart beats, and oxygen consumption of some small hibernators are remarkably reduced during hibernation [3,4], however no neural damage is observed in the brain of hibernators after arousal from torpor [2,5]. During the hibernation season, brain controls the initiation and end of the torpor phase, and is protected from damages induced by hypoxia-ischemia [6]. It is thought that brain regions play a pivotal role in the regulation of torpor-arousal cycle [7][8][9]. Up to date, the molecular basis for the regulation of torpor-arousal cycle has been extensively investigated in hibernators, such as the arctic ground squirrels and bears [10][11][12]. Evidences suggested that the global suppression of metabolism, protein synthesis, energy consumption, and neural excitotoxicity, as well as an increased antioxidant defense in the brain are strategies used by hibernators for neuroprotection [13][14][15]. The underlying molecular mechanisms by which hibernators prevent brain damages during hibernation may help us to develop the therapeutic strategies of human neurodegenerative disorders [16].
It has long been believed that selective gene expression plays a central role during the hibernation [17]. For example, the upregulation of a_2-macroglobulin in the liver tissue increases the blood clotting time in squirrels of torpid state, whereas the downregulation of prostaglandin D2 synthase in the brain tissue of hibernating squirrels is related to arousal state changes [17,18]. Thus, gene expression variations are critical for the adaptation to change environments. Over the past decades, it has been documented that genes involved in lipolysis, non-shivering thermogenesis, and neural plasticity tend to be differentially expressed between the summer active and winter hibernating states in some animals [10,19,20]. Most previous researches on mammalian hibernation have performed in squirrels [10,21], however, little is known about the molecular mechanisms in heterothermic bats.
Bats, belonging to the order Chiroptera, are one of the largest monophyletic clades in mammalian species [22][23][24]. They are the only mammals with true flight ability [25]. Most bat families contain heterothermic species, and the hibernation behavior of bats is confirmed in at least seven families [26]. Some works have been performed to screen hibernation-related genes in the brain of the greater horseshoe bat (R. ferrumequinum) using SSH library and dot blot [27]. The greater horseshoe bats are known to have prolific populations and use large body fat preserved as fuel to survive an annual hibernation period up to 6-to 8-months. Since bats fly immediately upon arousal and the brain plasticity has been suggested to play important roles in hibernation [18,20,25], it is interesting to explore the molecular mechanisms underlying brain function and protection during bat torpor. Thanks to the emergence of next-generation sequencing technology, which opens the door for us to obtain cDNA fragments from transcriptome with reasonably complete coverage in a reduced time scale and at a lower cost [28]. In this study, we performed next-generation sequencing to compare the brain transcriptomes of active and torpid greater horseshoe bats and to identify differentially expressed genes that contribute to the hibernation phenotype.

Brain tissue sample preparation and deep sequencing
All animals used in this work were carried out in strict accordance with the approval of Regulations for the Administration of Laboratory Animals ( No specific permissions were required for the greater horseshoe bats capture in the Fish cave. The surface temperatures of the active and torpid bats were approximately 25uC and 5uC, respectively. The greater horseshoe bats were euthanized by respiratory hyperanesthesia followed thoracotomy, and all efforts were made to minimize potential pain and suffering. The whole brain tissues (without cerebellum) were stored in RNAlater solution (Takara Biotechnology Co. Ltd., Dalian, China) to preserve the RNA state for use immediately after killing.
Total RNA was extracted using Trizol (Life Technologies Corp.) according to the manufacturer's protocols and treated with RNeasy Mini kit (Qiagen). RNA concentration was evaluated using a NanoDrop spectrophotometer and RNA quality was assessed by Agilent Bioanalyzer. Total RNA was submitted to

Data preprocess and de novo transcriptome assembly
To reduce sequencing errors, DynamicTrim Perl script implemented in SolexQA package [29] was performed to control the quality of raw sequencing data with the default parameters setting. Next, we de novo assembled the transcriptome sequences using Trinity package [30]. The assembled contigs with length less than 200 bp were discarded due to their low annotation rate. The program was run on 64-bit Linux system (Red Hat 6.0) with 256 GB internal memory. After transcriptome assembly, CD-Hit program [31] was then used to reduce sequence redundancy of transcriptome with default parameter setting.

Function annotation and RNA-Seq analysis
As for the functional annotation, assembled contigs were searched by BLASTX against the NCBI non-redundant nucleotide database with an E-value threshold of 1E-5. The results of the best blast hits were extracted, and coding sequences were subsequently determined. Functional annotation was implemented using the online software Blast2GO [32]. Gene ontology (GO) terms and annotations were downloaded on May 22, 2012. These identified genes were classified into different levels of GO categories.
We then estimated read counts of every gene using RSEM [33] nested in the Trinity package [30]. To compare the gene expression difference between the active and torpid states, we normalized the read counts by calculating FPKM (fragments per kilobase of exon per million fragments mapped) values. In this work, EdgeR method was used to detect differentially expressed genes [34]. After applying Benjamini-Hochberg correction for multiple test, the false discovery rate (FDR) ,0.05 and fold change greater than 2 were selected as the criteria for significant differences. In order to explore the biological implications, all differentially expressed genes were submitted to DAVID online toolkit [35].

De novo assembly and functional annotation
To obtain an overview of brain transcriptome of summer active and winter torpid greater horseshoe bat, cDNA samples were prepared from mixtures of RNA from the whole brain (without cerebellum) and then were sequenced using Illumina Genome Analyzer II. Each sample group was sequenced in two lanes, and a total of 51,755,965 and 38,558,209 paired-end reads with 75-bp length for active and torpid samples were generated ( Table 1). After cleaning of raw sequences, the trimmed reads were assembled using Trinity package [30]. As shown in Figure 1, contigs with their lengths between 200 and 350 bp were overrepresented, making up ,38% of the total number of contigs. We even identified 7,445 contigs longer than 3,050 bp. Next, we used CD-Hit program [31] to generate non-redundant contigs. A total of 11,279 contigs with 95% sequence identity were removed, and only the longest one was retained if more than one contigs correspond to a gene. Finally, a total of 11,214 genes were obtained representing about 10% of all distinct contigs. Of these annotated genes, 9,910 (93%) were assigned to at least one GO categories and small molecule metabolic process (GO:0044281) and protein binding (GO:0005515) were most represented ( Figure 2, Table S1).

Differential expression analysis
FPKM (fragments per kilobase of exon per million fragments mapped) values were used as proxies of gene expression level. To investigate the expression patterns of the active and torpid samples of the greater horseshoe bats, we evaluated the number of reads assembled from each library for each transcript and compared gene expression difference between them. In this work, edgeR method [34] was used to measure gene expression differences between summer active and winter torpid bats. A total of 1573 genes were identified to be significantly differentially expressed between these summer active and winter torpid samples with the criteria of F.D.R ,0.05 and fold change .2 (Table S2), within which 1038 genes are up-regulated in the torpid state and 535 genes are down-regulated in the torpid state. Among these differentially expressed genes, the gene with the most significant (P-value = 3.89E-202) expression difference between the active and torpid samples is heat shock 70 kDa protein 8 (HSPA8), with a ,4.46 fold higher expression in the torpor sample. HSPA8 belongs to the heat-shock cognate subgroup, which functions as an ATPase working with auxilin to remove clathrin coated vesicles [36]. Conversely, the gene with the most significant (Pvalue = 4.77E-154) expression difference and higher expression in the active sample (6.01 fold) is CKB, which encodes a cytoplasmic enzyme involved in energy homeostasis [37]. Table 2 lists the top 10 differentially expressed genes that were significantly altered between active and torpid samples.

Gene ontology enrichment analysis
In order to gain insights into the biological implications of these differentially expressed genes, we performed the GO enrichment analysis. A total of 261 GO categories belonging to molecular function and biology process sub-categories were identified. The significant overrepresented GO terms were listed in Table S3, and the enriched GO terms of significantly up-regulated or downregulated genes in the torpid brain were all listed in Table S4. In these GO categories, those responsible for energy metabolism, cellular stress responses, oxidative stress, and cytoskeleton organization were observed to be overrepresented during bat torpor. Some important functional categories are summarized below.
Cellular stress responses. GO categories of 'unfolded protein binding' and 'protein folding' are examples representing the cytoplasmic stress response. Most of the genes involving these GO categories were over-overexpressed during torpor (Figure 4). We found many heat shock proteins (HSPs) were differentially expressed. HSPs, a group of proteins induced by heat shock, play diverse roles in folding, regulation and degradation of other proteins; repair of DNA and chromatin and cell cycle control in response to outside stresses [38,39]. The differentially expressed HSP families include Hsp70 (HSPA8, HSPA5, HSPA9,  HSPA4L), Hsp90 (HSP90AA1, HSP90B1) and Hsp40 Oxidative stress. Oxidative stress is an imbalance between free radicals and cellular antioxidant defense, which can lead to cell death by damaging proteins, nucleic acids and lipids through free radicals or peroxides [40]. Genes involved in GO category of 'response to oxidative stress' are mainly down-regulated during the hibernation ( Figure 5). For example, the glutathione peroxidase 3 (GPX3) gene, which functions in the detoxification of hydrogen peroxide, has a ,8 fold lower expression in the torpid brain compared with the active brain.
Stabilization of microtubule. Next, GO term of 'regulation of cytoskeleton organization' was enriched among differentially expressed genes between active and torpid samples ( Figure 6). For example, two neural microtubules-associated protein genes (MAP1B and MAP2) which cooperatively play important roles in dendritic outgrowth and microtubule organization [41] show elevated expression during torpor. Previously, Chen et al. also found that MAP1B up-regulated in the torpid brain [27]. In addition, three cytoskeletal genes, dystonin (DST), microtubuleactin cross-linking factor 1 (MACF1) and adenomatous polyposis coli protein (APC), are also overexpressed during the hibernation. These genes might be important for the characterization of neural morphology in the torpid brain.

Discussion
Animal hibernation is an adaptive and complex process of metabolic suppression for energy conservation in winter. To understand the molecular mechanisms of mammalian hibernation, some genomic studies have been conducted using squirrels and bears as models. However, the transcriptomic changes in response to hibernation in bats remain largely unknown. In this work, we compared the brain transcriptome of greater horseshoe bats between active and torpid states. Hibernating bats can preserve fat during pre-hibernation period as a primary energy source for surviving in the hibernation season. The concentration of blood glucose is low during bat torpor [42]. The brain prefers to utilize ketones rather than glucose as the primary energy source during hibernation. [43]. Glycolysis, a key metabolic pathway of glucose catabolism, includes 10 sequential enzyme reactions to convert glucose into pyruvate and prepares for energy production. The low expression of nine genes (i.e., PFKP, PFKM, PFKL, ALDOA, ALDOC, ENO1, ENO2, GPI, and PKM) involved in glycolysis of torpid bats supports the theory of glucose conservation during mammalian torpor. In addition to the conservation of glucose, down-regulation of glycolytic pathway plays a central role in metabolic suppression during torpor. The co-downregulation of these genes in the torpid bats further  suggests a global suppression of metabolism at bat torpor, however, malate dehydrogenase 1 (MDH1), which produces oxaloacetate for gluconeogenesis, was elevated during torpor. This may be a continual pathway that prepares glucose for utilization of bats upon arousal. Most of the NADH dehydrogenase (ubiquinone) subunits were down-regulated in torpid bats suggested that electron transport is also suppressed during hibernation. Nonetheless, the protein subunits such as NDUFS1 and NDUFA5 were elevated in expression. This seems to agree with the up-regulation of some subunits of mitochondrial membrane proteins (e.g., NDUFS1, NDUFS3, NDUFV2) in brainstem of ground squirrels at earlier torpor, suggesting the common strategies adopted by small hibernators (i.e., bats and squirrels) in control of NADH dehydrogenase activity [44]. In the mitochondria, electron transport promotes ATP synthesis by generating electrochemical gradient across membranes in the aerobic environment. The results described above strongly suggest that a reduction of electron transport occurs in the brain of torpid bats to adapt low oxygen supply and to reduce the reactive oxygen species (ROS) production. Therefore, a global metabolic suppression is achieved.
Protein synthesis and proteolysis were suppressed during the torpor [3,45]. Our findings of the elevated expression of multifamily HSP (heat shock proteins) genes suggested the preservation of proteome by protein folding during bat torpor. Stresses, such as hypoxia or ischemia, induce protein mis-folding during animal hibernation, which also occur in the pathological conditions related to neurodegenerative diseases and brain injury. HSPA8, a highly abundant member of HSP70s, was an important factor for protein refolding, which reduces the neuronal loss in the neurodegenerative disease [46]. Moreover, overexpression of two mitochondrial HSPs (HSPA9, HSPD1) increases the brain resistance to oxidative stress from focal ischemia [47,48]. In addition, the overexpression of Hdj-2 (DNAJA1, HSP70 partner) could reduce the ischemia-like injury [49]. Elevated expression of HSPs in various tissue of hibernating animals shows that HSPs play specific role in protein homeostasis to maintain tissue specific functions [50][51][52]. In our study, the elevated expressions of multiply families of HSP gene may provide a potential mechanism for synergistically protecting the nervous system during bat torpor.
During the hibernation season, ischemic/reperfusion switch in the brain leads to oxidative stress [53,54]. Several proteins related to antioxidant defense are shown to be elevated in the torpid brain of bats [55]. PRNP is a prion protein that promotes the activity of superoxide dismutase 1 (SOD1). It has been documented that knock-out of this gene results in increased levels of oxidative stress in the brain [56,57]. Therefore, elevated expression level of PRNP may increase the antioxidant defense in the torpid brain. Because PRNP has diverse vital roles in nervous system and tumorigenesis [58,59], we hypothesized that up-regulation of PRNP is involved in the maintenance of neural morphology and anti-apoptosis during the torpor. Superoxide dismutase, an important antioxidant, eliminates free radicals in the body. Up-regulation of SOD2 could protect the neural cell from free radical damage in the torpor. Yan et al. also found the elevated levels of SOD2 in the brown adipose tissue of hibernating arctic ground squirrels [10]. Besides the typical antioxidant, Oliver et al. suggested that oxidation resistance protein 1 gene (OXR1) is essential for the protection of neurons from oxidative stresses [60]. Increased expressions of OXR1 enhance the sensitivity of neurons under stress conditions. In this study, we also found an elevated expression of the OXR1 gene in the torpid brain of bat.
During the hibernation season, hibernators undergo global temperature-mediated neural plasticity where dendrites, cell bodies and spines significantly retract in the torpor [61,62]. The differentially expressed cytoskeletal gene influences the stabilization of cytoskeletal network leading to the neural plasticity [20]. Although neural plasticity protects the brain, significant changes of neural structure can result in memory loss [61,63]. However, it has been shown that bats keep their long term memory of food-finding after hibernation [64]. In this study, the expressions of two microtubule stabilizers (i.e., MAP1 and MAP2) are found to be elevated in the torpid bat brain. Previous work suggested that a lower MAP2 gene expression reduces the microtubule density in dendrites and length of dendrites [65]. Arendt et al. also showed that a lower MAP2 density enhances the dendritic retraction due to reduction of the stable microtubules during the torpor [66]. These findings suggest that bats may reduce retraction of dendrites by stabilizing the microtubule in the dendrites through elevating MAP2 gene expression. Although the neuronal capacity is significantly reduced, hibernating mammals are capable of sensing periodical arousals and environmental stimulations [55]. Therefore, the highly expressed SNAP25 and SYT1 (involved in vesicle docking and synaptic transmission) identified in torpid bats suggest that the neural plasticity occurs in the brain at bat torpor.
In conclusion, we have detected hibernation-related genes in the brain of greater horseshoe bats by deep transcriptome sequencing. We have also identified differential expression of genes involved in glycolysis, respiratory electron transport chain, unfolded protein binding, protein folding, oxidative stress, and cytoskeleton organization. Results of this study will further our understanding on the molecular mechanism of mammalian hibernation.

Supporting Information
Table S1 Top 10 represented GO categories of biological process and molecular function.