Gene expression and adaptive evolution of ZBED1 in the hibernating greater horseshoe bat (Rhinolophus ferrumequinum)

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.

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

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. A_260/280 values were all above 2.0 and electrophoresis of the RNA samples demonstrated that 28S and 18S rRNA were intact.

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)₁₈ 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.

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⁻¹ 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.

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.

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×10⁶ 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 (d_N) and synonymous substitutions (d_S) (Yang, 2007). A d_N/d_S 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 χ² 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<ω₀<1), site class 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 ω₂=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).

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.

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

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.

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

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 ω₂=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).

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 C_CE&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 C_NE’) or the Central-East group (‘model C_CE’) as the foreground partition and all others comprised the background partition. Using model C_NE, the ω ratio of the Northeast group was estimated to be 2.83291, whereas under model C_CE, the ω ratio of the Central-East group was less than 1 [ω_2(CE)=0.59606]. LRT statistics suggested model C_CE&NE fits the data better compared with model C_NE and model C_CE. The comparison between model C_CE&NE and M1a was not statistically significant, but the comparison between model C_CE&NE and M2a_rel was statistically significant. We also compared model C_CE and model C_NE versus M1a and M2a_rel, respectively; neither of these comparisons were statistically significant, suggesting neither model C_CE nor model C_NE fit the data better than the null models. These results suggested that model C_CE&NE was more suitable for the data, and the Northeast and Central-East groups were under positive selection.

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.