Mammalian hibernators experience physiological extremes, e.g. ischemia, muscle disuse and hypothermia, which are lethal to non-hibernators, implying the existence of underlying mechanisms that allow hibernators to withstand these physiological extremes. Increased cell proliferation is suggested to be such a strategy, but its molecular basis remains unknown. In this study, we characterized the expression pattern of ZBED1 (zinc finger, BED-type containing 1), a transcription factor that plays a crucial role in regulating cell proliferation, in five tissues of the greater horseshoe bat (Rhinolophus ferrumequinum) during pre-hibernation, deep hibernation and post-hibernation. Moreover, we investigated the ZBED1 genetic divergence from individuals with variable hibernation phenotypes that cover all three known mtDNA lineages of the species. Expression analyses showed that ZBED1 is overexpressed only in brain and skeletal muscle, not in the other three tissues, suggesting an increased cell proliferation in these two tissues during deep hibernation. Evolutionary analyses showed that ZBED1 sequences were clustered into two well-supported clades with each one dominated by hibernating and non-hibernating individuals, respectively. Positive selection analyses further showed some positively selected sites and a divergent selection pressure among hibernating and non-hibernating groups of R. ferrumequinum. Our results suggest that ZBED1 as a potential candidate gene that regulates cell proliferation for hibernators to face physiological extremes during hibernation.

Hibernation is a strategy taken by hibernators to survive the pressures of a harsh environment such as low ambient temperature and/or scarce food availability (Wang et al., 1996). In order to survive, hibernators must suffer potential adverse effects induced by long periods of low body temperature, and repeated ischemia and perfusion, some of which are highly stressful and even lethal for non-hibernators (Carey et al., 2003). However, hibernators obviously have abilities to protect themselves from potential harm, and several underlying mechanisms are proposed, such as induced expression of stress (heat shock) proteins to combat cellular stress (Carey et al., 1999, 2000; Lei et al., 2014) and increased protein–ubiquitin conjugate concentration to help degrade damaged proteins during hibernation (Carey et al., 2000).

Among the protective mechanisms that underlie tolerance of physiological extremes is increased cell proliferation (Cerri et al., 2009). Both the number of apoptotic cells and the number of cells that undergo proliferation in the main brain areas of hibernating frogs increase compared with cell numbers in active frogs, suggesting a compensatory effect of increased cell proliferation for increased apoptotic cells during hibernation (Cerri et al., 2009). During hibernation, it is stressful for cells to be exposed to low temperatures for long periods because of its effects on processes such as protein stability, membrane function, ATP synthesis, activity of key regulatory enzymes and cytoskeletal integrity (Kandror and Goldberg, 1997; Sonna et al., 2002). Additionally, ischemia-induced oxidative stress during hibernation is considered to play a crucial role in cell death (Dave et al., 2012). Because of the high risk of cell death, increased cell proliferation might be an important protective mechanism for hibernating mammals to cope with potential cell damage during hibernation.

ZBED1 (zinc finger, BED-type containing 1), also called DREF (DNA replication-related element-binding factor), is a transcription regulatory factor playing a key role in the control of cell proliferation (Matsukage et al., 2008). ZBED1 regulates cell proliferation by controlling many cell proliferation-related genes, such as those encoding PCNA, Orc 2, DNA polymerase α 180 kDa and 73 kDa subunits, RFC140, raf, E2F, cyclin A and Skp A (Hirose et al., 1991, 1993; Ohno et al., 1996; Takahashi et al., 1996; Kandror and Goldberg, 1997; Ryu et al., 1997; Sawado et al., 1998; Phuong Thao et al., 2006; Tsuchiya et al., 2007; Matsukage et al., 2008). To date, no study about the role of ZBED1 in hibernation has been reported. Given the important role of ZBED1 in cell proliferation, it can be hypothesized to be a potential factor for helping physiological adaption in hibernators. If it is the case, an increased expression of ZBED1 in tissues vulnerable to stresses induced by physiological extremes during hibernation would be expected and adaptive evolution of the ZBED1 gene in hibernating taxa may occur compared with non-hibernating taxa.

The greater horseshoe bat Rhinolophus ferrumequinum Schreber 1774, has a wide distribution in China and is found from southern China to northern China, and three divergent mtDNA evolutionary groups (Northeast, Central-East, Southwest) are identified (Sun et al., 2013). The Northeast group, which is distributed in Northeast China, hibernates strictly from December to April during winter, whereas the Southwest group, which is distributed in Southwest China, where the weather is warm and food is available even during winter, is active year-round (Aulagnier et al., 2008). Many individuals in the Central-East group (e.g. Beijing and Jinan), also hibernate during winter. The different hibernation phenotypes therefore provide us with an ideal model system to study the potential role of ZBED1 in hibernation.

In this study, we characterized relative mRNA expression levels of the ZBED1 gene in five tissues including heart, kidney, brain, liver and skeletal muscle, from Northeast individuals of R. ferrumequinum during pre-hibernation, deep hibernation and post-hibernation, respectively. We also analyzed adaptive evolution of ZBED1 based on individuals from all three mtDNA groups. We aimed to (i) determine the expression pattern of ZBED1 in different tissues across the hibernation cycle and (ii) determine the potential genetic divergence and possible adaptive evolution of ZBED1 between hibernating and non-hibernating groups of R. ferrumequinum.

Gene expression analyses of ZBED1

Sample collection

A total of 16 female greater horseshoe bats, belonging to the Northeast group, were used in this study and these individuals were collected at different time points in the Jilin Province of Northeast China (Table 1). Three individuals were sampled in September in the Ground Cave (125°50′25″E, 41°4′8″N), referred to as pre-hibernation individuals. Three and five individuals were sampled in December and April, from the Ground Cave and Pigeon Cave (125°56′19″E, 41°28′29″N), respectively – referred to as deep hibernation individuals. Five bats were sampled in May from Pigeon Cave, referred to as post-hibernation individuals. The pre-hibernation individuals and post-hibernation individuals were active and had obvious responses to stimuli (sound, touch and light) and high body temperatures, whereas the deep hibernation individuals were in a torpid state, had no response to stimuli and body temperatures close to ambient temperature. All bats were killed by decapitation as soon as they were caught in the caves to ensure RNA profiles most representative of their active or hibernation state, as well as to minimize pain and suffering. Surgical procedures were promptly performed to protect RNA from degradation. Five tissue types (heart, kidney, brain, liver and skeletal muscle) of each animal were rapidly excised, flash-frozen in liquid nitrogen and then stored at −80°C until processed for RNA isolation. All experimental procedures carried out in this study were approved by the National Animal Research Authority of Northeast Normal University, China (approval number: NENU-20080416) and the Forestry Bureau of Jilin Province of China (approval number: [2006]178).

Table 1.

Sample information used for gene expression analyses

Sample information used for gene expression analyses
Sample information used for gene expression analyses

Total RNA isolation

Total RNA was isolated from bat tissues using TRIzol Reagent (Life Technologies, Carlsbad, CA, USA). Approximately 50 mg tissue was ground in liquid nitrogen using a pestle and mortar, and then solubilized in 1 ml ice-cold TRIzol reagent and incubated for 5 min. The cell debris was pelleted by centrifuging the homogenates at 12,000 g at 4°C for 10 min and 0.2 ml 100% chloroform was added to the supernatant, then tubes were shaken vigorously by hand for 15 s and incubated at room temperature for 5 min. Then the aqueous phase was transferred to a fresh tube after centrifuging at 12,000 g at 4°C for 15 min and an equal volume of 100% chloroform was then added. After shaking tubes vigorously by hand for 15 s, they were incubated on ice for 2–3 min and centrifuged at 12,000 g at 4°C for 10 min to extract again. Next, 0.5 ml of isopropyl alcohol was added to the aqueous phase transferred into a fresh tube and the RNA was precipitated at −80°C overnight. The next day, the supernatant was removed from the tube after centrifugation at 12,000 g at 4°C for 10 min. Then the RNA pellet was washed with 1 ml of 75% ethanol by vortexing and centrifuging at 7500 g, 4°C for 5 min. Finally, the RNA pellet was air-dried for 5–10 min, and then dissolved in RNase-free water. Non-denaturing agarose gel electrophoresis and a NanoDrop spectrophotometer (Thermo Fisher Scientific, Waltham, MA, USA) were used to assess the quality and quantity, respectively, of the isolated RNA. A260/280 values were all above 2.0 and electrophoresis of the RNA samples demonstrated that 28S and 18S rRNA were intact.

cDNA synthesis

A total of 80 RNA samples were converted to cDNA by using TransScript One-Step gDNA Removal and cDNA Synthesis SuperMix (TransGen Biotech, Beijing, China) following the manufacturer's instructions. The reaction mixture contained approximatly 1 μg RNA template, 0.5 μg Anchored Oligo (dT)18 primer and 6 μl RNase-free water and was incubated at 65°C for 5 min, then on ice for 2 min. In the next step, a mixture containing 2×TS Reaction Mix, TransScript RT/RI Enzyme Mix and gDNA Remover was added into the mixture from the first step. The final 20 μl reaction mixture was incubated at 42°C for 30 min, and finally 85°C for 5 min.

Real-time quantitative PCR

Quantitative RT-PCR was performed using StepOne Real-Time PCR System (Applied Biosystems) with an automatic threshold calculated by StepOne software v.2.1. The primer pairs of ZBED1 and reference gene β-actin specified for R. ferrumequinum are listed in Table S1, including their sequences, annealing temperatures and product lengths. The efficiencies of these two genes are all approximately 100%, which ensures that the ΔΔCt equation is used correctly. The total 10 μl PCR mixture contained 5 μl Thunderbird Sybr qPCR Mix (Toyobo), 0.2 μl 50× ROX reference dye, 1 μl cDNA template and 0.25 μmol l−1 of each primer. The PCR was performed under the following conditions: pre-denaturation at 95°C for 1 min, then 40 cycles: 95°C, 15 s; 60°C for 1 min, with data collection after each cycle. Next is melting/dissociation curve analysis.

The relative quantity was calculated by using the 2−ΔΔCt method (Livak and Schmittgen, 2001). The fold changes of ZBED1 gene in five tissue types between September (pre-hibernation) and December (deep hibernation) were calculated relative to September; between April (deep hibernation) and May (post-hibernation) were calculated relative to April. The relative expression folds were expressed as mean±s.e. and statistical significance between active and torpid samples was analyzed using an independent-samples t-test. In further analysis, the fold changes of ZBED1 in kidney, brain, liver and skeletal muscle were calculated relative to heart. The statistical significance among the five tissues was analyzed by one-way ANOVA followed by post hoc Tukey's test (SPSS v.19.0). Differences between means were statistically significant when P<0.05.

Adaptive evolution analyses of ZBED1

Samples collection, cloning and sequencing

Forty-six individuals of R. ferrumequinum were sampled from 13 different localities covering all three known mtDNA groups (Fig. 1, Table S2). One individual of Rhinolophus macrotis was also sampled to be used as an outgroup for phylogenetic analyses. Biopsy punches (3 mm diameter) were used to acquire bat wing membrane (plagiopatagium) samples as outlined by Worthington Wilmer and Barratt (1996). The 3 mm holes in each wing heal in ∼4 weeks (Worthington Wilmer and Barratt, 1996). Genomic DNA was extracted from the wing biopsies using the Animal Genomic DNA Isolation Kit (Sangon, China). Primer pairs for amplification of ZBED1 coding domain sequence (CDS) were designed according to cDNA sequence (KU312321) of ZBED1 and whole genome shotgun sequence of R. ferrumequinum (Accession numbers: AWHA01106959.1 and AWHA01106960.1) using Primer Premier 5.0. Sense. Primer sequence was 5′-GCTGTTTGTCTCACTGCTT-3′ and anti-sense primer sequence was 5′-ACCCCATCTGTCTGTCACT-3′. To reduce mutations introduced by amplification, Q5 super fidelity 2× Master Mix (Biolabs) was used according to the manual instructions. Purification of PCR product was completed by using a TaKaRa MiniBEST Agarose Gel DNA Extraction Kit. Cloning of ZBED1 CDS was performed using a pEASY Blunt Zero Cloning Kit (TransGen Biotech). Five positive clones of each individual were randomly selected to be sequenced by Shanghai Sangon Biotechnology. All CDS sequences of ZBED1 were deposited in GenBank with accession numbers KR868578-KR868649.

Fig. 1.

Distribution of R. ferrumequinum sampling sites in China. The corresponding three mtDNA groups (Northeast, Central-East and Southwest) are outlined with dashed lines. Numbers 1–13 represent 13 isolated sampling sites.

Fig. 1.

Distribution of R. ferrumequinum sampling sites in China. The corresponding three mtDNA groups (Northeast, Central-East and Southwest) are outlined with dashed lines. Numbers 1–13 represent 13 isolated sampling sites.

Data analyses

All 46 individuals sampled were sequenced successfully. Some individuals were determined to be heterozygotes and both allele sequences from the heterozygote individuals were used for subsequent analyses. Sequences were aligned using ClustalW implemented in MEGA 6 (Tamura et al., 2013). Polymorphic site detection and haplotype analysis were implemented in Dnasp 5.0 (Librado and Rozas, 2009). The sequence alignment was analyzed with RDP 3 (Martin et al., 2010), implementing several algorithms to detect recombinant sequences.

The phylogenetic relationship of the ZBED1 gene was constructed using MrBayes 3.2.0 (Ronquist and Huelsenbeck, 2003). The DNA substitution model that best fit the data was estimated using the Akaike information criterion using Modeltest 3.7 (Posada and Crandall, 2005). Bayesian inference was run with four simultaneous chains, each of 1×106 generations, sampled every 100 generations and the first 25% of trees were discarded as ‘burn-in’. The ZBED1 CDS of R. macrotis was also cloned and sequenced, which was used as the outgroup.

The adaptive evolution of the ZBED1 gene was determined by positive selection analyses. For this, the CODEML program in PAML package was used to derive maximum-likelihood estimates of the rate of nonsynonymous substitutions (dN) and synonymous substitutions (dS) (Yang, 2007). A dN/dS ratio (ω) of >1 signifies positive selection, a ratio of ∼1 signifies neutrality and a ratio of <1 signifies purifying selection. Branch models, branch-site models, site-specific models (M1a, M2a, M3, M7, M8), and clade model (model C) were implemented in the CODEML program. Each model was fit to the data multiple times with different initial omega values to help ensure local optima were avoided. Likelihood ratio tests (LRTs) were used to compare the model fit of complex models against simpler, nested models. LRTs were carried out by comparing twice the difference in ln likelihood scores of nested models against a χ2 distribution with the degrees of freedom equal to the number of extra parameters estimated by the more complex model. The species tree was constructed based on the published mtDNA phylogenetic tree of R. ferrumequinum (Sun et al., 2013), and the inner relationships within each of the three mtDNA groups were respectively constructed using MrBayes 3.2.0 and nucleotide sequences of the ZBED1 CDS. The best-fit models selected for 1st, 2nd and 3rd codon sites within the Northeast group were TIM, HKY and HKY, respectively. Similarly, the best-fit models selected for 1st, 2nd and 3rd codon sites within the Central-East and Southwest group were K81uf, TIM and HKY+I+G, and TIM, F81 and HKY+G, respectively. Based on the species tree, branch models were used to identify branches under positive selection. For this, a one-ratio model and a free-ratio model were initially employed, and then two-ratio models that assumed the positive selection of the ancestral branches of the Northeast group and the Central-East group of R. ferrumequinum were used, respectively. In addition, we also applied two-ratio models to the branches with ω>1 identified under the free-ratios model, which assumes an independent ω ratio for each branch, to identify potential positively selected branches in the phylogeny of R. ferrumequinum.

Branch-site models were applied to the ancestral branches of the Northeast and the Central-East groups to identify the certain sites under positive selection. Estimates of site-wise ω values in these foreground branches were then compared with estimates across the remaining background branches in the phylogeny. Under model A, the ω values were assigned to four predefined site classes. Site class 0 was estimated from data but constrained (0<ω0<1), site class 1 ω1=1, site class 2a ω2a could exceed 1 on the foreground but is constrained to be under purifying selection on the background, and site class 2b ω2b could exceed 1 on the foreground but not on the background. The null model fixes ω2=1. The comparison between these two models is so-called test 2, with one degree of freedom (d.f.). In addition, test 1 was also used for the analyses, which compares model A with M1a (nearly neutral), with d.f.=2.

A clade model was applied to test potential divergent selection between hibernating and non-hibernating groups of R. ferrumequinum. We compared ω ratios averaged across branches within focal clades (foreground) with ω ratios estimated for the rest of the tree (background) (Bielawski and Yang, 2004). In clade models, we applied a ‘multi-clade’ approach, which allows for more than two tree partitions to examine complex patterns of divergence in selection across the phylogeny by comparing such models against simpler, nested models with fewer tree partitions (Weadick and Chang, 2012a). Models M1a and M2a_rel were used as null models for model C in LRTs. M1a is the standard null model for model C. However, model C versus M1a LRT was reported prone to false-positive test results when faced with moderate among-site variation in selective constraint (Weadick and Chang, 2012,b). M2a_rel null model proposed by Weadick et al. (2012b) was therefore also used in this study.

Positively selected sites were analyzed using both PAML and HyPhy. With PAML, we compared site models M2a (selection) versus M1a (neutral), M7 (beta) versus M8 (beta and ω). We also compared models M3 (discrete) versus M0 (one ratio) to test change of ω ratio among sites. The Bayes empirical Bayes (BEB) was used to identify sites under positive selection. In HyPhy we used open-source software packages available under the datamonkey web-server (http://www.datamonkey.org/; Pond and Frost, 2005a,b; Pond and Muse, 2005; Murrell et al., 2012), to implement five different methods, including single-likelihood ancestor counting (SLAC), fixed-effects likelihood (FEL), random-effects method (REL), mixed-effects model of evolution (MEME) and fast unbiased Bayesian approximation (FUBAR).

Differential expression of ZBED1

The Northeast group of R. ferrumequinum is known to have a strict hibernation behavior during winter, and here we analyzed its ZBED1 expression level in five different tissues across pre-hibernation, deep hibernation and post-hibernation. Three tissues (heart, kidney and liver) consistently presented a decreased expression during deep hibernation (December and April) compared with both pre-hibernation (September) and post-hibernation (May). By contrast, the two other types of tissue, brain and skeletal muscle, showed an increased expression during deep hibernation and decreased expression in both pre-hibernation and post-hibernation (Fig. 2A,B). However, only the expression differences of ZBED1 in brain and skeletal muscle between April and May were statistically significant, with P-values of 0.011 and 0.017, respectively. The smaller sample size in September and December than in April and May may limit the statistical difference of ZBED1 expression in brain and skeletal muscle between September and December.

Fig. 2.

Relative quantification of ZBED1 gene expression in five tissues of R. ferrumequinum. mRNA expression was measured in heart, kidney, brain, liver and skeletal muscle from pre-hibernation (September) to deep hibernation (December) (A), and from deep hibernation (April) to post-hibernation (May) (B). The statistical significance was analyzed using an independent-samples t-test. Asterisks indicate significant differences between the two states (P<0.05). Error bars indicate s.e.

Fig. 2.

Relative quantification of ZBED1 gene expression in five tissues of R. ferrumequinum. mRNA expression was measured in heart, kidney, brain, liver and skeletal muscle from pre-hibernation (September) to deep hibernation (December) (A), and from deep hibernation (April) to post-hibernation (May) (B). The statistical significance was analyzed using an independent-samples t-test. Asterisks indicate significant differences between the two states (P<0.05). Error bars indicate s.e.

We further analyzed the relative expression levels of ZBED1 for each of the five tissue types and for the four time points (September, December, April and May) (Fig. 3). The results show that ZBED1 was expressed most strongly in skeletal muscle compared with the other four tissues (Fig. 3). Liver showed the least expression (Fig. 3). Expression of ZBED1 in kidney, brain and liver showed no significant difference at the three time points: May, December and April, with P-values >0.05 (Fig. 3B–D), except for expression at the active state in September, which showed significant differences among these three tissues (Fig. 3A). The differences in expression between heart and kidney were not significant at the four time points (Fig. 3).

Fig. 3.

Differential expression of the ZBED1 gene in different R. ferrumequinum tissues collected at different times. Samples were collected in September (A), May (B), December (C) and April (D). Individual mRNA expression values are normalized to the mean value in heart, and calculated as 2−ΔΔCt. Normalized mRNA expression values for multiple individuals were averaged by tissue. Different letters indicate significant differences among tissues (Tukey's test, P<0.05). Error bars in the figure are s.e.

Fig. 3.

Differential expression of the ZBED1 gene in different R. ferrumequinum tissues collected at different times. Samples were collected in September (A), May (B), December (C) and April (D). Individual mRNA expression values are normalized to the mean value in heart, and calculated as 2−ΔΔCt. Normalized mRNA expression values for multiple individuals were averaged by tissue. Different letters indicate significant differences among tissues (Tukey's test, P<0.05). Error bars in the figure are s.e.

Genetic analysis

Among 46 individuals analyzed, 26 individuals were determined to be heterozygous, and a total of 39 haplotypes were detected (Table S3). Among the 39 haplotypes, 64 polymorphic sites were found, of which 7 sites (1092, 1206, 1281, 1335, 1344, 1555 and 1626) were homozygous and identical between the Northeast group and the Central-East group but different from the Southwest group (Table 2). No evidence of recombinant sequence was found in our dataset. To determine the genetic divergence of the ZBED1 gene within the species, the Bayes tree based on the complete CDS of the ZBED1 gene was constructed (Fig. 4A). The best-fit models for 1st, 2nd and 3rd codon sites were TIM, TIM and K81uf+I+G, respectively. Accordingly, 39 haplotypes of the ZBED1 gene were clustered into two well-supported clades, with one clade containing only individuals from the Southwest group, and the other one containing individuals from both the Northeast group and the Central-East group. Unlike the mtDNA result, the haplotypes of the ZBED1 gene from the Northeast group and the Central-East group did not form two clades, but were clustered into one clade.

Table 2.

Single-nucleotide polymorphisms of ZBED1 CDS among three groups: Northeast, Central-East and Southwest

Single-nucleotide polymorphisms of ZBED1 CDS among three groups: Northeast, Central-East and Southwest
Single-nucleotide polymorphisms of ZBED1 CDS among three groups: Northeast, Central-East and Southwest
Fig. 4.

R. ferrumequinum phylogenetic trees. (A) The Bayes phylogenetic tree of the ZBED1 gene of R. ferrumequinum. The gene tree was based on complete ZBED1 CDS. Posterior probabilities over 0.50 are shown on tree branches. Squares and circles indicate haplotypes that belong to the Northeast group and the Central-East group, respectively. The haplotypes shared between the Northeast and the Central-East groups are indicated by both squares and circles. (B) The phylogenetic tree of R. ferrumequinum. The topologic relationship among the Northeast group, the Central-East group and the Southwest group is based on previous mtDNA phylogenetic results (Sun et al., 2013). The inner relationship within each group was constructed by MrBayes 3.2.0 using ZBED1 CDS in this study. Hap20 and Hap21, marked with ‘N’ and ‘C’, are haplotypes shared between the Northeast and the Central-East groups, respectively. Unlike haplotypes from the Northeast and the Central-East groups, which form one clade labeled with Central-East and Northeast in the gene tree (A), the haplotypes from the Northeast and the Central-East groups formed two separate clades in the species tree (B).

Fig. 4.

R. ferrumequinum phylogenetic trees. (A) The Bayes phylogenetic tree of the ZBED1 gene of R. ferrumequinum. The gene tree was based on complete ZBED1 CDS. Posterior probabilities over 0.50 are shown on tree branches. Squares and circles indicate haplotypes that belong to the Northeast group and the Central-East group, respectively. The haplotypes shared between the Northeast and the Central-East groups are indicated by both squares and circles. (B) The phylogenetic tree of R. ferrumequinum. The topologic relationship among the Northeast group, the Central-East group and the Southwest group is based on previous mtDNA phylogenetic results (Sun et al., 2013). The inner relationship within each group was constructed by MrBayes 3.2.0 using ZBED1 CDS in this study. Hap20 and Hap21, marked with ‘N’ and ‘C’, are haplotypes shared between the Northeast and the Central-East groups, respectively. Unlike haplotypes from the Northeast and the Central-East groups, which form one clade labeled with Central-East and Northeast in the gene tree (A), the haplotypes from the Northeast and the Central-East groups formed two separate clades in the species tree (B).

Positive selection analysis

To determine the potential adaptive evolution of the ZBED1 gene, different models implemented in PAML were employed. In branch models, the one-ratio model (M0) showed a single ω ratio (ω0=0.01926) for the entire tree (Table 3). The free-ratio model showed high variable ω ratios among branches, of which seven branches (a–g in Fig. 4B) were considerably higher than 1 (Table 3), although the LRT statistic of the free-ratio model versus the one-ratio model was not significant with 2Δl=60.90 (P=0.924 and d.f.=78). The seven branches were further designated as the foreground branches and a two-ratio branch model was then analyzed (Table 3). Among seven branches, only five branches had a significantly better fit than the one-ratio model (Table 3). Interestingly, all these five branches belonged to the Northeast group and the Central-East group. To test whether the ancestral branches of the Northeast group (branch NE) and the Central-East group (branch CE) underwent positive selection, branches NE and CE were then respectively defined as foreground branches with others being background branches, and ω ratios of these two branches were calculated to be 2.01765 and 1.48366, respectively. The LRT statistic showed that the two-ratio models were not significantly better than the one-ratio model. We further used branch-site model A to identify possible positively selected sites of NE and CE branches. Whereas the LRT statistics of both test 1 and test 2 were not significant, the ω2a and ω2b ratios for the ancestral NE and CE branches were higher than 1 (Table 4).

Table 3.

Positive selection analyses of the ZBED1 gene based on branch models

Positive selection analyses of the ZBED1 gene based on branch models
Positive selection analyses of the ZBED1 gene based on branch models
Table 4.

Positive selection analyses of the ZBED1 gene based on branch-site models

Positive selection analyses of the ZBED1 gene based on branch-site models
Positive selection analyses of the ZBED1 gene based on branch-site models

Subsequently, different site-specific models (M1a, M2a, M3, M7 and M8) were further analyzed (Table 5). M2a (positive selection) identified one amino acid (324A) under positive selection at the 72.4% probability, whereas the LRT statistic of the M1a–M2a comparison was not significant (P=0.664). M8 showed that five sites were under positive selection with ω=1.00952, and one (324A) of the five sites showed a high probability of 93.2%. The difference between M7 and M8 was statistically significant (2Δl=13.70, d.f.=2, P=0.001). The analyses based on datamonkey web-server also supported that amino acid 324A was under positive selection, because at least two methods, REL and FUBAR, significantly identified 324A to be under positive selection with ω>1 (Table S4). The discrete model (M3) with K=3 site classes suggested that 2.067% of sites were under positive selection with ω2=1.00949 and identified 10 amino acids under positive selection at 100% probability. The LRT statistic of M3(K=3)-M0 comparison was significant (2Δl=45.62, d.f.=4, P<0.001).

Table 5.

Positive selection analyses of the ZBED1 gene based on site models

Positive selection analyses of the ZBED1 gene based on site models
Positive selection analyses of the ZBED1 gene based on site models

We used clade models (branch-site model C) to test for potential divergent selection pressure among our three focal groups (Table 6). We first applied a multi-clade model C assuming three partitions: Northeast group, Central-East group and Southwest group. Under this model, which we call ‘model CCE&NE’, a small proportion of the dataset (only 1.819%) evolved under divergent selective constraint across the three partitions, with a ω ratio less than 1 along the SW branches [ω2(SW)=0.00000], above 1 along the Northeast group [ω2(NE)=2.07667] and slightly higher than 1 along the Central-East group [ω2(CE)=1.07229]. Then we used model C with either the Northeast group (‘model CNE’) or the Central-East group (‘model CCE’) as the foreground partition and all others comprised the background partition. Using model CNE, the ω ratio of the Northeast group was estimated to be 2.83291, whereas under model CCE, the ω ratio of the Central-East group was less than 1 [ω2(CE)=0.59606]. LRT statistics suggested model CCE&NE fits the data better compared with model CNE and model CCE. The comparison between model CCE&NE and M1a was not statistically significant, but the comparison between model CCE&NE and M2a_rel was statistically significant. We also compared model CCE and model CNE versus M1a and M2a_rel, respectively; neither of these comparisons were statistically significant, suggesting neither model CCE nor model CNE fit the data better than the null models. These results suggested that model CCE&NE was more suitable for the data, and the Northeast and Central-East groups were under positive selection.

Table 6.

Positive selection analyses of the ZBED1 gene based on clade models

Positive selection analyses of the ZBED1 gene based on clade models
Positive selection analyses of the ZBED1 gene based on clade models

Mammalian hibernators experience physiological extremes, e.g. ischemia, muscle disuse and hypothermia, that are lethal to non-hibernators but have no deleterious effects on hibernators. This could suggest the existence of underlying mechanisms that allow hibernators to withstand these physiological extremes. To date, several possible protective mechanisms have been proposed (Carey et al., 2000; Cerri et al., 2009; Lei et al., 2014) and among them, increased cell proliferation appears to represent a major protective mechanism to compensate for increased apoptotic cells during hibernation (Cerri et al., 2009).

In this study, using R. ferrumequinum as a model system, we evaluated the expression and adaptive evolution of an important transcription factor (ZBED1) that plays a key role in the control of cell proliferation (Matsukage et al., 2008). Gene expression analyses showed that out of the five tissues analyzed, the ZBED1 gene is overexpressed only in the brain and skeletal muscle tissue during deep hibernation when compared with pre-hibernation and post-hibernation. This tissue-specific overexpression may suggest that these tissues have a high level of cell proliferation. The increased expression of ZBED1 in the brain observed in this study is consistent with the study of Cerri et al. (2009), in which increased cell proliferation in the brain of a species of hibernating frog (Rana esculenta) is demonstrated. In addition, an increased cellular proliferation in tests in some hibernators, e.g. ground squirrels and woodchucks, is suggested to account for their full sexual activity during hibernation (Wimsatt, 1969). Increased cell proliferation in certain tissues found in these studies and our own study might suggest tissue heterogeneity, because many other specific tissues (e.g. lymphocytes, liver and hippocampus) are suggested to undergo a reduced cellular proliferation during hibernation (Shivatcheva and Hadjioloff, 1987; Kruman, 1992; Kiss et al., 2011; Popov et al., 2011). Different tissues might suffer different pressures during hibernation and hence might need variable protective mechanisms. Given the stress that the brain suffers during hibernation, such as ischemia-induced oxidative stress and thermal stress (Carey et al., 2003; Dave et al., 2012), the increased ZBED1 expression in the brain during deep hibernation observed in this study could be associated with cell proliferation and hence may be helpful to compensate for cell death and the onset of neurological damages. Further analyses showed that the expression of the gene in the skeletal muscle is consistently higher than that of the other four tissues used in this study. The overexpression of the ZBED1 gene in skeletal muscle tissue might imply a high level of cell proliferation. The possible cell proliferation might partly account for the lack of muscle atrophy in hibernators after extended periods of disuse during hibernation. Inactivity in humans, such as confined bed rest, weightlessness or limb immobilization, can lead to atrophy of skeletal muscle, loss of muscle tone and impaired strength (Harlow et al., 2001). However, in hibernating animals, such as ground squirrels and bats, although their skeletal muscle mass is reduced by 14–65% depending on muscle type and species, neither atrophy nor dysfunction is found in their skeletal muscle (Steffen et al., 1991; Wickler et al., 1991). The possible increase of cell proliferation in skeletal muscle could therefore help to maintain the function of skeletal muscle during hibernation.

Phylogenetic analyses of the ZBED1 gene showed a strong genetic divergence between the Northeast group and the Southwest group of R. ferrumequinum, which represent hibernating and non-hibernating groups, respectively. This would imply a functional divergence between them. Positive selection analyses further showed some evidence of the adaptive evolution of the ZBED1 gene. The branch model detected five positive selected branches, which were found to be restricted to the Northeast and the Central-East groups that are dominated by hibernating individuals. Under the site model, several positively selected sites were detected, of which, one site, 324A, was found to be with a high probability under positive selection. The amino acid 324A is putatively located in the ribonuclease H-like catalytic domain in the protein structure of ZBED1 (Mitchell et al., 2015). Our clade model supports the existence of divergent selection pressures of the ZBED1 gene among the three mtDNA groups. Estimates of ω were >1 for ∼2% of sites when the Northeast and Central-East groups (ωNE=2.07667, ωCE=1.07229) were defined as the foreground clades. This may indicate a strong positive selective pressure on the Northeast group and a slight positive selective pressure on the Central-East group. Nevertheless, in the Southwest group, we detected no evidence for positive selection. These findings remained unchanged when the ZBED1 phylogenetic tree (Fig. 4A) was used for our positive selection analyses (data not shown). Identification of the positively selected branches and positive selective pressure on the Northeast and Central-East groups support our prediction that adaptive evolution of ZBED1 might occur in the hibernating groups of R. ferrumequinum.

In this study, using R. ferrumequinum as a model system, we investigated the expression and adaptive evolution of the ZBED1 gene, an important cell proliferation-related transcription factor, and found a potential role of the ZBED1 gene in hibernation. We found an overexpression of the ZBED1 gene only in the brain and skeletal muscle, and detected some evidence of adaptive evolution of this gene in hibernating groups of R. ferrumequinum. Our results suggest that ZBED1 may be a potential candidate gene that regulates cell proliferation for hibernators to face the physiological extremes of hibernation.

We are thankful to Hongjun Lin and Bo Luo for their great contributions to samples collection and Yuyang Hao for her help in the experiments. We are grateful to Katie Andrea Solari for her help in language revision.

Author contributions

J.F., Y.W. and K.S. conceived and designed the experiments. Y.X., Y.W., H.W., T.J., A.L., X.H., X.Y. and L.S. collected the samples and performed the experiments. J.F. contributed reagents, materials and analysis tools. Y.X. wrote the paper. All authors read and approved the final manuscript.

Funding

This study was supported by the National Natural Science Foundation of China [grant numbers 91131003, 31270414, 31370399].

Aulagnier
,
S.
,
Hutson
,
A. M.
,
Spitzenberger
,
F.
,
Juste
,
J.
,
Karataş
,
A.
,
Palmeirim
,
J.
and
Paunovic
,
M.
(
2008
).
Rhinolophus ferrumequinum. The IUCN Red List of Threatened Species 2008: e.T19517A8947355
.
Bielawski
,
J. P.
and
Yang
,
Z.
(
2004
).
A maximum likelihood method for detecting functional divergence at individual codon sites, with application to gene family evolution
.
J. Mol. Evol.
59
,
121
-
132
.
Carey
,
H. V.
,
Sills
,
N. S.
and
Gorham
,
D. A.
(
1999
).
Stress proteins in mammalian hibernation
.
Am. Zool.
39
,
825
-
835
.
Carey
,
H. V.
,
Frank
,
C. L.
and
Seifert
,
J. P.
(
2000
).
Hibernation induces oxidative stress and activation of NK-kappaB in ground squirrel intestine
.
J. Comp. Physiol. B Biochem. Syst. Environ. Physiol.
170
,
551
-
559
.
Carey
,
H. V.
,
Andrews
,
M. T.
and
Martin
,
S. L.
(
2003
).
Mammalian hibernation: cellular and molecular responses to depressed metabolism and low temperature
.
Physiol. Rev.
83
,
1153
-
1181
.
Cerri
,
S.
,
Bottiroli
,
G.
,
Bottone
,
M. G.
,
Barni
,
S.
and
Bernocchi
,
G.
(
2009
).
Cell proliferation and death in the brain of active and hibernating frogs
.
J. Anat.
215
,
124
-
131
.
Dave
,
K. R.
,
Christian
,
S. L.
,
Perez-Pinzon
,
M. A.
and
Drew
,
K. L.
(
2012
).
Neuroprotection: lessons from hibernators
.
Comp. Biochem. Phys. B Biochem. Mol. Biol.
162
,
1
-
9
.
Harlow
,
H. J.
,
Lohuis
,
T.
,
Beck
,
T. D. I.
and
Iaizzo
,
P. A.
(
2001
).
Muscle strength in overwintering bears
.
Nature
409
,
997
.
Hirose
,
F.
,
Yamaguchi
,
M.
,
Nishida
,
Y.
,
Masutani
,
M.
,
Miyazawa
,
H.
,
Hanaoka
,
F.
and
Matsukage
,
A.
(
1991
).
Structure and expression during development of Drosophila melanogaster gene for DNA polymerase alpha
.
Nucleic Acids. Res.
19
,
4991
-
4998
.
Hirose
,
F.
,
Yamaguchi
,
M.
,
Handa
,
H.
,
Inomata
,
Y.
and
Matsukage
,
A.
(
1993
).
Novel 8-base pair sequence (Drosophila DNA replication-related element) and specific binding factor involved in the expression of Drosophila genes for DNA polymerase alpha and proliferating cell nuclear antigen
.
J. Biol. Chem.
268
,
2092
-
2099
.
Kandror
,
O.
and
Goldberg
,
A. L.
(
1997
).
Trigger factor is induced upon cold shock and enhances viability of Escherichia coli at low temperatures
.
Proc. Natl. Acad. Sci. USA
94
,
4978
-
4981
.
Kiss
,
A. J.
,
Muir
,
T. J.
,
Lee
,
R. E.
Jr
. and
Costanzo
,
J. P.
(
2011
).
Seasonal variation in the hepatoproteome of the dehydration- and freeze-tolerant wood frog, Rana sylvatica
.
Int. J. Mol. Sci.
12
,
8406
-
8414
.
Kruman
,
I. I.
(
1992
).
Comparative analysis of cell replacement in hibernators
.
Comp. Biochem. Physiol. A Physiol.
101
,
11
-
18
.
Lei
,
M.
,
Dong
,
D.
,
Mu
,
S.
,
Pan
,
Y.-H.
and
Zhang
,
S.
(
2014
).
Comparison of brain transcriptome of the greater horseshoe bats (Rhinolophus ferrumequinum) in active and torpid episodes
.
PLoS ONE
9
,
e107746
.
Librado
,
P.
and
Rozas
,
J.
(
2009
).
DnaSP v5: a software for comprehensive analysis of DNA polymorphism data
.
Bioinformatics
25
,
1451
-
1452
.
Livak
,
K. J.
and
Schmittgen
,
T. D.
(
2001
).
Analysis of relative gene expression data using real-time quantitative PCR and the 2(-Delta Delta C(T)) method
.
Methods
25
,
402
-
408
.
Martin
,
D. P.
,
Lemey
,
P.
,
Lott
,
M.
,
Moulton
,
V.
,
Posada
,
D.
and
Lefeuvre
,
P.
(
2010
).
RDP3: a flexible and fast computer program for analyzing recombination
.
Bioinformatics
26
,
2462
-
2463
.
Matsukage
,
A.
,
Hirose
,
F.
,
Yoo
,
M.
and
Yamaguchi
,
M.
(
2008
).
The DRE/DREF transcriptional regulatory system: a master key for cell proliferation
.
Biochim. Biophys. Acta
1779
,
81
-
89
.
Mitchell
,
A.
,
Chang
,
H.-Y.
,
Daugherty
,
L.
,
Fraser
,
M.
,
Hunter
,
S.
,
Lopez
,
R.
,
McAnulla
,
C.
,
McMenamin
,
C.
,
Nuka
,
G.
,
Pesseat
,
S.
, et al. 
(
2015
).
The InterPro protein families database: the classification resource after 15 years
.
Nucleic Acids. Res.
43
,
D213
-
D221
.
Murrell
,
B.
,
Wertheim
,
J. O.
,
Moola
,
S.
,
Weighill
,
T.
,
Scheffler
,
K.
and
Pond
,
S. L. K.
(
2012
).
Detecting individual sites subject to episodic diversifying selection
.
PLo.S Genet.
8
,
e1002764
.
Ohno
,
K.
,
Hirose
,
F.
,
Sakaguchi
,
K.
,
Nishida
,
Y.
and
Matsukage
,
A.
(
1996
).
Transcriptional regulation of the Drosophila CycA gene by the DNA replication-related element (DRE) and DRE binding factor (DREF)
.
Nucleic Acids Res.
24
,
3942
-
3946
.
Phuong Thao
,
D. T.
,
Ida
,
H.
,
Yoshida
,
H.
and
Yamaguchi
,
M.
(
2006
).
Identification of the Drosophila skpA gene as a novel target of the transcription factor DREF
.
Exp. Cell Res.
312
,
3641
-
3650
.
Pond
,
S. L. K.
and
Frost
,
S. D.
(
2005a
).
Not so different after all: a comparison of methods for detecting amino acid sites under selection
.
Mol. Biol. Evol.
22
,
1208
-
1222
.
Pond
,
S. L. K.
and
Frost
,
S. D. W.
(
2005b
).
Datamonkey: rapid detection of selective pressure on individual sites of codon alignments
.
Bioinformatics
21
,
2531
-
2533
.
Pond
,
S. L. K.
and
Muse
,
S. V.
(
2005
).
HyPhy: hypothesis testing using phylogenies
. In
Statistical Methods in Molecular Evolution
(
ed. R. Nielsen
), pp.
125
-
181
.
New York: Springer
.
Popov
,
V. I.
,
Kraev
,
I. V.
,
Ignat'ev
,
D. A.
and
Stewart
,
M. G.
(
2011
).
Suspension of mitotic activity in dentate gyrus of the hibernating ground squirrel
.
Neural Plast.
2011
,
867525
.
Posada
,
D.
and
Crandall
,
K.
(
2005
).
Modeltest 3.7. Program Distributed by the Author
.
Spain
:
Universidad de Vigo
.
Ronquist
,
F.
and
Huelsenbeck
,
J. P.
(
2003
).
MrBayes 3: bayesian phylogenetic inference under mixed models
.
Bioinformatics
19
,
1572
-
1574
.
Ryu
,
J.-R.
,
Choi
,
T.-Y.
,
Kwon
,
E.-J.
,
Lee
,
W.-H.
,
Nishida
,
Y.
,
Hayashi
,
Y.
,
Matsukage
,
A.
,
Yamaguchi
,
M.
and
Yoo
,
M.-A.
(
1997
).
Transcriptional regulation of the Drosophila-raf proto-oncogene by the DNA replication-related element (DRE)/DRE-binding factor (DREF) system
.
Nucleic Acids Res.
25
,
794
-
799
.
Sawado
,
T.
,
Hirose
,
F.
,
Takahashi
,
Y.
,
Sasaki
,
T.
,
Shinomiya
,
T.
,
Sakaguchi
,
K.
,
Matsukage
,
A.
and
Yamaguchi
,
M.
(
1998
).
The DNA replication-related element (DRE)/DRE-binding factor system is a transcriptional regulator of the Drosophila E2Fgene
.
J. Biol. Chem.
273
,
26042
-
26051
.
Shivatcheva
,
T. M.
and
Hadjioloff
,
A. I.
(
1987
).
Seasonal involution of gut-associated lymphoid tissue of the European ground squirrel
.
Dev. Comp. Immunol.
11
,
791
-
799
.
Sonna
,
L. A.
,
Fujita
,
J.
,
Gaffin
,
S. L.
and
Lilly
,
C. M.
(
2002
).
Invited review: effects of heat and cold stress on mammalian gene expression
.
J. Appl. Physiol.
92
,
1725
-
1742
.
Steffen
,
J. M.
,
Koebel
,
D. A.
,
Musacchia
,
X. J.
and
Milsom
,
W. K.
(
1991
).
Morphometric and metabolic indices of disuse in muscles of hibernating ground squirrels
.
Comp. Biochem. Phys. B Comp. Biochem.
99
,
815
-
819
.
Sun
,
K.
,
Luo
,
L.
,
Kimball
,
R. T.
,
Wei
,
X.
,
Jin
,
L.
,
Jiang
,
T.
,
Li
,
G.
and
Feng
,
J.
(
2013
).
Geographic variation in the acoustic traits of greater horseshoe bats: testing the importance of drift and ecological selection in evolutionary processes
.
PLoS ONE
8
,
e70368
.
Takahashi
,
Y.
,
Yamaguchi
,
M.
,
Hirose
,
F.
,
Cotterill
,
S.
,
Kobayashi
,
J.
,
Miyajima
,
S.
and
Matsukage
,
A.
(
1996
).
DNA replication-related elements cooperate to enhance promoter activity of the Drosophila DNA polymerase alpha 73-kDa subunit gene
.
J. Biol. Chem.
271
,
14541
-
14547
.
Tamura
,
K.
,
Stecher
,
G.
,
Peterson
,
D.
,
Filipski
,
A.
and
Kumar
,
S.
(
2013
).
MEGA6: molecular evolutionary genetics analysis version 6.0
.
Mol. Biol. Evol.
30
,
2725
-
2729
.
Tsuchiya
,
A.
,
Inoue
,
Y. H.
,
Ida
,
H.
,
Kawase
,
Y.
,
Okudaira
,
K.
,
Ohno
,
K.
,
Yoshida
,
H.
and
Yamaguchi
,
M.
(
2007
).
Transcriptional regulation of the Drosophila rfc1 gene by the DRE-DREF pathway
.
FEBS J.
274
,
1818
-
1832
.
Wang
,
L. C.
,
Lee
,
T.
,
Fregley
,
M.
and
Blatteis
,
C.
(
1996
).
Torpor and hibernation in mammals: metabolic, physiological, and biochemical adaptations
. In
Handbook of Physiology: Environmental Physiology
(ed.
J. R.
Pappenheimer
,
M. J.
Fregly
and
C. M.
Blatties
), pp. 507-532.
New York
:
Oxford University Press
.
Weadick
,
C. J.
and
Chang
,
B. S. W.
(
2012a
).
Complex patterns of divergence among green-sensitive (RH2a) African cichlid opsins revealed by Clade model analyses
.
BMC Evol. Biol.
12
,
206
.
Weadick
,
C. J.
and
Chang
,
B. S. W.
(
2012b
).
An improved likelihood ratio test for detecting site-specific functional divergence among clades of protein-coding genes
.
Mol. Biol. Evol.
29
,
1297
-
1300
.
Wickler
,
S. J.
,
Hoyt
,
D. F.
and
Van Breukelen
,
F.
(
1991
).
Disuse atrophy in the hibernating golden-mantled ground squirrel, Spermophilus lateralis
.
Am. J. Physiol.
261
,
R1214
-
R1217
.
Wimsatt
,
W. A.
(
1969
).
Some interrelations of reproduction and hibernation in mammals
. In
Dormancy and Survival
(ed.
A. H. W.
Woolhouse
), pp.
511
-
549
.
London
:
Cambridge University Press
.
Worthington Wilmer
,
J.
and
Barratt
,
E.
(
1996
).
A non-lethal method of tissue sampling for genetic studies of chiropterans
.
BAT Res. News
37
,
1
-
3
.
Yang
,
Z.
(
2007
).
PAML 4: phylogenetic analysis by maximum likelihood
.
Mol. Biol. Evol.
24
,
1586
-
1591
.

Competing interests

The authors declare no competing or financial interests.

Supplementary information