Down but Not Out: The Role of MicroRNAs in Hibernating Bats

MicroRNAs (miRNAs) regulate many physiological processes through post-transcriptional control of gene expression and are a major part of the small noncoding RNAs (snRNA). As hibernators can survive at low body temperatures (Tb) for many months without suffering tissue damage, understanding the mechanisms that enable them to do so are of medical interest. Because the brain integrates peripheral physiology and white adipose tissue (WAT) is the primary energy source during hibernation, we hypothesized that both of these organs play a crucial role in hibernation, and thus, their activity would be relatively increased during hibernation. We carried out the first genomic analysis of small RNAs, specifically miRNAs, in the brain and WAT of a hibernating bat (Myotis ricketti) by comparing deeply torpid with euthermic individual bats using high-throughput sequencing (Solexa) and qPCR validation of expression levels. A total of 196 miRNAs (including 77 novel bat-specific miRNAs) were identified, and of these, 49 miRNAs showed significant differences in expression during hibernation, including 33 in the brain and 25 in WAT (P≤0.01 &│logFC│≥1). Stem-loop qPCR confirmed the miRNA expression patterns identified by Solexa sequencing. Moreover, 31 miRNAs showed tissue- or state-specific expression, and six miRNAs with counts >100 were specifically expressed in the brain. Putative target gene prediction combined with KEGG pathway and GO annotation showed that many essential processes of both organs are significantly correlated with differentially expressed miRNAs during bat hibernation. This is especially evident with down-regulated miRNAs, indicating that many physiological pathways are altered during hibernation. Thus, our novel findings of miRNAs and Interspersed Elements in a hibernating bat suggest that brain and WAT are active with respect to the miRNA expression activity during hibernation.


Introduction
Hibernation (multiday torpor) and daily torpor in heterothermic mammals are characterized by substantial, temporal reductions of energy expenditure and body temperature (T b ). Torpor is often expressed when animals are exposed to low ambient temperature (T a ), and the decrease in energy use is crucial for survival during adverse conditions and food shortages [1]. Mammalian hibernation is an extreme example of hypometabolism with metabolic rate depression often to <5% of the euthermic rate, a minimum T b of approximately 0°C, and cell preservation during prolonged dormancy periods [2,3]. The physiological changes, such as the decrease and subsequent increase of metabolism and T b during entry and arousal from torpor, are active and adaptive processes. They are often initiated in response to changes in environmental temperature and photoperiod and require the coordination of core and peripheral tissues by the brain [4,5]. Importantly, the carbohydrate oxidation pathway shifts upon hormone activation towards lipid metabolism, and fat (mainly stored in white adipose tissue, WAT) becomes the primary energy source during the hibernation season [6]. Considering the brain's important role in integrating peripheral tissues and WAT to provide energy during hibernation, we hypothesized that to adequately perform these tasks, both must be activated during hibernation, in contrast to other tissues and organs, where metabolism is generally suppressed [7].
A detailed understanding of the molecular mechanisms of this phenomenon would improve our understanding of hibernation at the subcellular level, but may also be useful in human medicine for improving organ transplantation, preventing skeletal muscle atrophy during long-term bed rest, and suppressing carcinogenesis [7]. In the past ten years, molecules that regulate hibernation and the molecular mechanisms underlying hibernation have been identified by transcriptome analysis [8][9][10][11][12]. Thousands of genes and proteins are differentially expressed in many species during hibernation [12][13][14][15][16][17]. However, the precise and timely control of gene expression during hibernation remains unclear.
There are more than 1100 species of bats (order Chiroptera), and they are among the most geographically widespread mammals [29]. Most bat families include species that enter some form of torpor, whereas multiday torpor (i.e., hibernation) occurs in at least five families [30]. The wide taxonomic diversity of heterothermy in bats suggests that bats are a good model to explore the evolutionary history of mammalian heterothermy [31]. Importantly, unlike ground squirrels, for which most molecular work on hibernation has been conducted, bats are not strongly seasonal in their torpor expression and many use multiday torpor throughout the year [32]. Therefore, a systematic study of genes or proteins involved in bat hibernation by comparing active and hibernating individuals will help us better understand the molecular mechanisms of hibernation and to explore the evolution of mammalian heterothermy [11,31,33,34]. However, a systematic analysis of small RNAs, especially miRNAs, with regard to function and evolution of torpor in bats is currently lacking.
The latest 7× coverage genome sequence of Myotis lucifugus (http://www.ensembl.org/ Myotis_lucifugus/Info/Index) provides us with the opportunity to survey miRNA expression in hibernating bats and to conduct comparative analyses with sequences from other metazoan species. Moreover, the RNA-seq sequencing-by-synthesis technique allows to sequence millions of short cDNAs and increases sensitivity by reducing personal error in library construction [35]. Based on the Illumina/Solexa sequencing platform and combined with qPCR validation, we performed the first in-depth miRNA analysis on two key tissues (brain and WAT) of the hibernating bat Myotis ricketti, which is congeneric with M. lucifugus. Our study aimed to discover bat-specific miRNAs, illustrate their potential role during bat hibernation, and contribute to the understanding of the evolution of mammalian heterothermy.

Ethics statement
All procedures involving animals were performed by strictly following the Guidelines and Regulations for the Administration of Laboratory Animals (

Sample and RNA preparation
A total of six Myotis ricketti individuals with a body mass of 18-22 g were used in this study. Three active M. ricketti bats (T b~3 6°C) were captured using mist nets in the Fangshan area (Beijing, China; 39°48 0 N, 115°42 0 E) on October 26, 2008 (T a = 21°C). On March 6, 2009, three bats in deep hibernation (T b~9°C ) were captured in the same location (T a = 7°C). The field studies did not involve endangered or protected species, and no specific permissions were required for M. ricketti in the Fangshan area. The bats were sacrificed humanely by decapitation immediately after measuring T b and body mass. Whole brain and white intra-abdominal adipose tissue (WAT) of hibernating and euthermic bats were rapidly removed, immediately frozen in liquid nitrogen, and stored at -80°C until RNA extraction.
Total RNA was isolated from entire brains and adipose tissues (~0.2 g) of six bats using the RNAiso kit (TakaRa, Japan) according to the manufacturer's instructions, and the concentration and RNA Integrity Number (RIN) was determined by an Agilent 2100 bioanalyzer. RNA samples were pooled to construct the transcriptome sequencing library as described previously [12,21]. Briefly, total RNA from bat brains of the same stage were mixed in equal amounts (10 μg each) into two pooled samples, representing the hibernation and active state, hereafter referred to as HB (hibernating brain) and AB (active brain). Similarly, total WAT RNA from hibernating and active bats (10 μg each) were pooled and referred to as HA (hibernating WAT), or AA (active WAT). In total, four RNA pools (two tissues in two physiological states) were used to construct small RNA libraries. The RNA quality and quantity were as follows: HB: 480 ng/μl, RIN = 8. Construction and high-throughput sequencing of small RNA libraries Small RNA cDNA libraries were constructed as described previously [36]. Briefly, small RNA (18~33 nt) was purified and enriched from each total RNA sample by polyacrylamide gel electrophoresis (PAGE), and 10 μg small RNA was ligated with the proprietary adapters used for cDNA synthesis and library construction, then Solexa sequencing-by-synthesis was performed (Illumina). All small RNA data series were submitted to NCBI Gene Expression Omnibus (GEO) with the accession number GSE29053.

Filter of small RNA reads
Individual sequence reads with base quality scores were produced by Illumina/Solexa. After removing contaminant reads (adapters, low quality and redundant reads), clean unique reads were mapped onto the Myotis lucifugus genome (http://asia.ensembl.org/Myotis_lucifugus/ Info/Index) using the Bowtie program [37]. Perfectly mapped reads were scanned against the metazoan mature miRNA in the Sanger miRBase (Release 19) [38] to identify orthologs of known miRNAs using the Patscan program [39] with two mismatches allowed. Non-conserved unique reads were screened against Rfam databases (Release 10) [40] to filter the sequences of tRNA, rRNA, snoRNA, and other ncRNAs except miRNAs using Bowtie. Reads that matched to the genome more than 20 times were removed by miRDeep [41], and reads sequenced only one time were also removed. Finally, the remaining reads were considered potential miRNA reads and were used for miRNA identification.

MiRNA identification
MiRNA and its antisense strand can both be sequenced from the same precursor. Thus, to avoid repeated predictions and to reduce calculations, candidate reads whose distance in the reference genome were <200 nt were combined and examined as one genomic block. For each block, 150 nt of upstream and downstream extensions were extracted for secondary structure prediction. First, inverted repeats (IR) with stem-loop or hairpin structures were identified by Einverted of Emboss [42], with the following parameters: threshold = 30, match score = 3, mismatch score = 3, gap penalty = 6, and maximum repeat length = 240, as described previously [43]. The IR secondary structure was then predicted by RNAfold [44] with 10 nt upstream and downstream extensions and evaluated by MirCheck [43].
Predicted precursors that miRNA and antisense strand can be found in both arms were deemed miRNA candidates. Moreover, if the candidates have two or more unique reads resulted from the ±2 nt bias during cleavage and located at mature positions, they were deemed highly probable. When several length variants of the same miRNA were sequenced, only variants with the highest representation were considered. Finally, the miRNA candidates were submitted again to miRBase, and the miRNA precursors (hairpins) that passed MirCheck were manually inspected against the canonical miRNA structure to remove false predictions.

Repeat-derived siRNA detection
To screen repeat-derived small interfering RNA (siRNAs), the repeat sequences of the M. lucifugus genome were annotated by Repeatmasker (http://www.repeatmasker.org/) [45]. Reads that perfectly matched the M. lucifugus genome were aligned to repetitive elements using Bowtie, and reads that perfectly matched the repeats were considered genomic repeats-derived siRNAs.

Differential miRNA expression analysis and target prediction
Differentially expressed miRNAs that were statistically significant in relative abundance (reflected by Transcript Per Million, TPM) between hibernating and active states were identified by the edgeR function in Bioconductor [46]. Empirical Bayes estimation and exact tests based on the negative binomial distribution were used, and P0.01 & │logFC│1 considered statistically significant.
As there is no or little information about 3' untranslated region (UTR) of M. lucifugus reference genes in database of Ensembl (http://www.ensembl.org/Myotis_lucifugus/Info/Index), to predict miRNA target genes, we first identified the orthologs of M. lucifugus genes by searching against Homo sapiens mRNA (http://www.ncbi.nlm.nih.gov/RefSeq/) with tblastx (e value 1e-10 and identity 60%) [47]. Orthologs with the best hits were kept. The H. sapiens mRNAs were used to analyze the 3' UTR lengths, and 93.6% (67,594/72,204) of human mRNAs had a 3' UTR less than 3 Kb (see S1 Fig), excluding 32,559 mRNAs that did not have a 3' UTR. Thus, the 3 Kb region downstream of M. lucifugus reference gene coding sequences were extracted and aligned with the 3'UTR of H. sapiens mRNA orthologs by CLUSTALW [48] to trim the candidate sequences of 3' UTR of M. lucifugus. For M. lucifugus reference genes that did not have a human ortholog, the 3 Kb region downstream of the coding sequence was considered the 3'UTR. The miRNA target was then predicted using the PITA program based on the interaction between miRNAs and their targets with the default criterion and ΔΔG−10 kcal/mol [49]. Finally, target sequence Gene Ontology (GO) annotation was performed by Interproscan [50] and target sequences were compared to the Kyoto Encyclopedia of Genes and Genomes database (KEGG, release 50) by BLASTX at E values 1e-10 [47,51]. A Perl script was used to retrieve the KO (KEGG Orthology) information from BLAST results and to establish pathway associations between target genes and databases.

Quantitative miRNA expression by qPCR assay
Stem-loop reverse transcription (RT) qPCR was performed to quantify the expression level of 10 differentially expressed miRNAs obtained by Solexa sequencing in brains and WAT as described previously [52], and 5S rRNA was used as an endogenous control. Primers used in this study are listed in S1 Table. Each RT reaction contained 1 μg total RNA, 1× Reaction Buffer, dNTPs (5 mM each), 50 nM miRNA-specific stem-loop RT primer, 100 U RevertAid Premium Reverse Transcriptase and 0.5 U RiboLock RNase Inhibitor. RT reactions were incubated at 37°C for 30 min, 42°C for 60 min and 85°C for 5 min. PCR reactions were performed in a total of 20 μl, including 0.8 μl RT product, 10 μl 2× SYBR Green I Master Mix, 6 μM each of forward primer and reverse primer with the following program: 95°C for 10 min, followed by 40 cycles of 95°C for 15 s, 60°C for 30 s and 72°C for 15 s, then 95°C for 15 s, 60°C for 1 min and a final hold at 4°C. All reactions were performed on three biological replicates in each tissue and each physiological state, and each was run three times. Negative controls containing all reagents except template were included on each reaction plate. The 2 −ΔΔCT method was used to calculate relative expression (fold change); and data were analyzed by One-Way ANOVA using SPSS 18.0 software.

Small RNA sequencing and statistics
We isolated total RNA from hibernating and active Myotis ricketti brain and intra-abdominal adipose tissue, and the RNA concentration and RNA Integrity Number (RIN) met the requirements for Solexa small RNA library construction and sequencing.
Less than 5% of total small RNA reads from the bat libraries perfectly matched the Myotis lucifugus genome >20 times (S4 Table); thus, we removed these unique reads to avoid self-contradiction of miRNA prediction. Therefore, a total of 3.16×10 7 remained after mapping the sequences to the M. lucifugus genome (Table 1 and S2 Fig). The size distribution of perfectly matched small RNA reads is shown in Fig 1A. We identified a total of 2.47×10 7 potential miRNA reads, including 159,807 distinct reads, by searching against miRBase and filtering the ncRNAs and low-expressing reads (S5 Table). Subsequently, 1.51×10 7 (47.7% in total) miRNA reads met the fold-back structure (hairpin) and MirCheck criteria, and miRNA reads were the most abundant fraction (range from 32.7% in HA to 55.6% in AB, Table 1 and Fig 1B). Other small RNAs, such as ncRNAs, genomic repeats, transcript repeats, and unknown genomic regions, comprised 20.4%, 2.3%, 2.8% and 26.8%, respectively (Table 1 and Fig 1B). Upon further inspection of genomic repeat-derived siRNAs, we determined that Short INterspersed Element (SINE) and Long INterspersed Element (LINE) were the major contributors to the four libraries (Table 1 and Fig 2).

MiRNA identification
To better understand the length distribution of metazoan miRNA precursors, we analyzed known miRNA precursors in five out-grouped organisms (two vertebrates and three invertebrates, S3 Fig) in miRBase. The majority of pre-mature miRNAs from these five organisms were less than 150 nt with a mean length of approximately 90 nt (S3 Fig). We mapped the potential miRNA reads onto the reference genome and retrieved 509,206 block sequences with 150 nt upstream and downstream extensions; we then used these sequences for secondary structure prediction. From these, we used Einverted of Emboss [42] and identified 380,605 potential reads with inverted repeats. We identified 90 conserved hairpins in 28 families and 101 novel hairpins (S6 Table).

Target prediction and annotation
To better understand the potential role and mechanism of miRNAs in hibernation, we predicted the miRNA target genes.  Table). In contrast, further analysis showed that many pathways in both tissues were markedly affected by the differentially expressed miRNAs, especially by the miRNAs down-regulated during hibernation (S9 Table). We further examined the adipocytokine signaling pathway to determine the relationship between miRNAs and target mRNAs during hibernation, as adipocytokine signaling is co-regulated by both brain and adipose tissues. Moreover, we searched for mRNAs from our ongoing project of digital gene expression (DGE) sequencing to enrich for target genes. DGE sequencing has been conducted on the same animals used here (unpublished data). A total of 64 mRNAs involved in adipocytokine signaling pathway were identified by DGE sequencing (S10 Table, DGE unpublished data). By combining the miRNA and target mRNA expression patterns, we found that 41 of 49 differentially expressed miRNAs regulated the expression of 51 genes, and more than half of the genes in each tissue were up-regulated during hibernation.  Among these, the differential expression of three mRNAs (Insulin Receptor Substrate, IRS; Adenosine 5'-monophosphate (AMP)-activated protein kinase, AMPK; Retinoid X receptor, RXR) was significant (P0.01), including IRS over-expression. The others were down-regulated during HA/AA (Fig 4 and S10 Table).

Discussion
Our study shows that in the hibernating bat Myotis ricketti, 33 brain and 25 white adipose tissue (WAT) miRNAs were differentially expressed between euthermia and hibernation. Consequently, they are likely involved in regulating functional processes and energy metabolism during hibernation.
Our findings were possible with the ultrahigh-throughput sequencing technique (RNA-seq) Illumina/Solexa sequencing, which provides a great platform to discover and analyze small RNAs. It has been previously used in non-model species, such as the Arctic ground squirrel [Urocitellus (Spermophilus) parryii] and the sea cucumber (Apostichopus japonicus) during hypometabolism, but these species' genomic backgrounds and miRNA data are not available in miRbase [18,21]. The clear advantages of Solexa sequencing are the generation of several million small RNA sequences in each small RNA library of one run, the ability to identify novel miRNAs, the potential to overcome the drawbacks of microarrays, and the generation of expression profiles in a greater and reproducible dynamic range [53]. However, this technology can be influenced by sequencing errors and raw data processing before miRNA identification [53]. In the present study, we performed Solexa high-throughput sequencing with stem-loop qPCR validation of expression levels to analyze miRNAs involved in hibernation. Moreover, we chose a widely distributed bat species in China, congeneric with Myotis lucifugus, as the animal model to explore miRNA expression in two vital organs (brain and WAT) during hibernation. Our data show that the sequencing information from four libraries was saturated, and 88.7% unique reads were tissue/state-specifically expressed (S3 Table). These data indicate a high quality of small RNA library construction and screening and that the libraries meet the needs of small RNA annotation and analysis. After mapping the unique reads to the M. lucifugus genome database and removing the reads that matched genome 20, the small RNA annotation indicated that miRNAs (47.7% in total) are the most abundant fraction of small RNAs (see Table 1 and Fig 1). These results are similar to those of previous miRNA profiling of multi-animal species [18,21,36]. Furthermore, we found that SINE and LINE, which are retroposons widely used as phylogenetic markers [54,55], were the major source of genomic repeatderived siRNAs in all four libraries. The detailed information of SINE and LINE provided us with the ability to further explore the phylogenetic evolution of heterothermy in mammals, especially in bats.
Our data also revealed a large number of conserved and novel miRNAs. We identified a total of 196 mature miRNAs, including 119 conserved miRNAs and 77 novel miRNAs (S7 Table). Moreover, we found that 31 (22 conserved and nine novel) of 196 miRNAs were expressed in a tissue/state-specific manner, suggesting a complex mechanism of miRNA regulation during hibernation. Of those, six conserved miRNAs with expression counts >100 were specifically expressed in the brain, including miR-153, which has been shown to function in neuronal protection and tumorigenesis [56,57], miR-124, which regulates metabolic and immune processes [58,59], miR-551, which is implicated in the stress response [60], and miR-1298, which has been associated with neural and endocrinological disease [61]. However, in the absence of direct and detailed evidence between the molecular data and animal function, it is necessary to further investigate the roles of these tissue-and state-specific miRNAs during hibernation.
Of 196 miRNAs, 49 were differentially expressed (P0.01 & │logFC│1) ( Table 2 and S7  Table), including ten miRNAs that had been shown to regulate mammalian hibernation in other species (S8 Table) [18,20,[22][23][24][25]. Both microarray and qPCR have been used to validate miRNA expression in many studies, as they are more accurate than RNA-seq sequencing [18,21]. In our study, qPCR validation of 10 out of 49 differentially expressed miRNAs provided consistent results with those of Solexa sequencing (Table 3). We conclude that our data are reliable. Although we validated a few differentially expressed miRNAs by qPCR in our study, the limitations of RNA-seq, and the value of sample pooling in high-throughput sequencing, we suggest that microarray hybridization should be performed in the future to further validate the significance of differentially expressed miRNAs obtained here. Microarray analysis could potentially compensate for the drawbacks of RNA-seq sequencing [53] and estimate the biological variability.
To gain insight into the potential broader functions of the miRNAs we identified, we predicted their putative targets and classified them by KEGG pathway assay and GO annotation. Comparing the targets of 196 miRNAs, we found that miRNAs were differentially expressed in both brain and adipose tissues between euthermia and hibernation states, and this is likely a reflection of the different activities of many essential processes during these states (e.g., lipid metabolism and signal transduction, S9 Table). We found that only 49 of 196 miRNAs regulated most target gene expression (12,909 of 20,370), and thus most physiological processes during hibernation, contributing to reduced energy expenditure during hypometabolism. Moreover, both KEGG and GO assays identified the same enriched functional groups (P<0.05), and all of them were affected by differentially expressed miRNAs, particularly by down-regulated miRNAs (S9 Table). By binding to target mRNAs, miRNAs can reversibly degrade their targets at the transcriptional level and/or inhibit their translation [19]. Thus, during torpor, brain and WAT physiological processes were globally affected by down-regulated miRNAs, suggesting that these processes were activated, supporting our hypothesis. Moreover, previous studies have shown that some specific brain regions such as the hypothalamus are involved in regulating mammalian hibernation [62][63][64]. Our study was limited to the use of entire brains, but future work on specific brain regions may identify the precise regulation of miRNA expression in important brain areas. Such studies could further expand our knowledge of hibernation regulation. The brain plays an important role in controlling and synchronizing peripheral physiology to reduce energy expenditure and minimizing tissue damage and/or disease during hibernation [3,4,7]. In contrast, adipose tissue is the primary energy source during hibernation and is especially critical during periodic rewarming [6,65]. Thus, specific molecular changes are required to accommodate the specific roles of brain and adipose tissue during hibernation. Because the potential network of miRNAs and targets is highly complicated, we chose to examine the adipocytokine signaling pathway, which plays an important role in signaling transduction and energy metabolism [66,67]. In both tissues, more than half mRNAs were upregulated ( Fig 4A) in contrast to miRNA expression, indicating that this pathway was activated to some extent. According to KEGG analysis, the adipocytokine signaling pathway can be subdivided into three interrelated parts, TNFα (Tumor Necrosis Factor α)-related genes, LEP (Leptin)-related genes and ADIPO (Adiponectin)-related genes (Fig 4B). Andrews et al. [13] previously showed that energy metabolism of hibernators is shifted at the gene level from glycometabolism in the liver and/or pancreas to lipid metabolism in adipose tissue during hibernation. We found that the significant up-regulation of IRS in WAT could be indicative of glucose uptake inhibition, whereas a decrease in AMPK in WAT could promote gluconeogenesis. Our data help to illustrate how miRNAs and mRNAs control the metabolic switch by inhibiting the glycometabolic pathway to burning fat as the main energy source, supporting previous reports [65]. While the regulation of fat metabolism is crucial during hibernation, the regulation of food intake is important during euthermia. Negative AMPK expression or the inhibition of AMPK and AMPK-related gene activity in the hypothalamus result in reduced food intake [68,69]. We hypothesize that down-regulation of AMPK and RXR in WAT lead to decreased food intake and increased energy expenditure. Network construction increases our understanding of molecular mechanisms mediated by miRNAs and target genes. It also provides clues to uncover how brain and adipose tissue cooperate to regulate hibernation by decreasing food intake, inhibiting glucose metabolism and increasing fatty acid metabolism to meet the physiological needs of hibernators.

Conclusions
We identified 196 miRNAs (77 novel miRNAs) in bats by in-depth analysis of the small RNAs in the hibernating species Myotis ricketti. Differential miRNA expression suggests that some physiological pathways in the brain and adipose tissue are activated during hibernation. In contrast to other peripheral tissues that are physiologically suppressed, our findings advance the understanding of the molecular mechanisms of hibernation. However, further study is required to validate the expression pattern during hibernation vs. euthermia, estimate biological variability, and elucidate the mechanisms of specific miRNAs and their target mRNAs in hibernation. Moreover, the Interspersed Elements (LINE and SINE) identified here may help further explore the phylogenetic evolution of heterothermy in bats. With the progress of the M. lucifugus genome sequencing project and the increased accuracy of gene annotation, more detailed information on miRNAs involved in hibernation will likely be uncovered from bat transcriptome libraries. These can be used to better understand heterothermy the physiology and evolution of heterothermy.    Table. Classification of miRNA target genes. KEGG pathway analysis (A) and GO annotation (B) of target genes regulated by 196 miRNAs, 119 conserved miRNAs, 8 conserved miRNAs differentially expressed in the brain, 13 conserved miRNAs differentially expressed in adipose tissue, 77 novel miRNAs, 25 novel miRNAs differentially expressed in the brain, and 12 novel miR-NAs differentially expressed in adipose tissue. T-tests were performed to analyze the bio-activity of each physiological process in both tissues by comparing hibernation vs. active state. (XLS) S10