Hibernation in bats (Mammalia: Chiroptera) did not evolve through positive selection of leptin

Abstract Temperature regulation is an indispensable physiological activity critical for animal survival. However, relatively little is known about the origin of thermoregulatory regimes in a phylogenetic context, or the genetic mechanisms driving the evolution of these regimes. Using bats as a study system, we examined the evolution of three thermoregulatory regimes (hibernation, daily heterothermy, and homeothermy) in relation to the evolution of leptin, a protein implicated in regulation of torpor bouts in mammals, including bats. A threshold model was used to test for a correlation between lineages with positively selected lep, the gene encoding leptin, and the thermoregulatory regimes of those lineages. Although evidence for episodic positive selection of lep was found, positive selection was not correlated with lineages of heterothermic bats, a finding that contradicts results from previous studies. Evidence from our ancestral state reconstructions suggests that the most recent common ancestor of bats used daily heterothermy and that the presence of hibernation is highly unlikely at this node. Hibernation likely evolved independently at least four times in bats—once in the common ancestor of Vespertilionidae and Molossidae, once in the clade containing Rhinolophidae and Rhinopomatidae, and again independently in the lineages leading to Taphozous melanopogon and Mystacina tuberculata. Our reconstructions revealed that thermoregulatory regimes never transitioned directly from hibernation to homeothermy, or the reverse, in the evolutionary history of bats. This, in addition to recent evidence that heterothermy is best described along a continuum, suggests that thermoregulatory regimes in mammals are best represented as an ordered continuous trait (homeothermy ← → daily torpor ← → hibernation) rather than as the three discrete regimes that evolve in an unordered fashion. These results have important implications for methodological approaches in future physiological and evolutionary research.


| INTRODUC TI ON
In mammals, the evolution of endothermy-the use of metabolically produced heat to maintain T B (body temperature)-has enabled survival of many different taxa across a variety of ecological niches and ecosystems (Fristoe et al., 2015;Hillenius & Ruben, 2004;Scholander, Hock, Walters, & Johnson, 1950). Homeothermy is defined by maintenance of a relatively constant T B over time (Ivanov, 2006), and endothermy and homeothermy are frequently linked in mammals. Endothermic heterothermy, characterized by a combination of metabolically produced heat and the use of torpor-a physiological state characterized by the reduction in metabolism, and subsequently T b , below basal levels-has enabled the survival of mammalian taxa globally, especially in extremely cold climates, highly variable climates, or regions with periodic constraints on food and water availability (Geiser & Stawski, 2011;Tattersall et al., 2012). The duration of discrete torpor bouts is often used to define the two traditionally recognized forms of heterothermy in mammals: hibernation and daily heterothermy, also known as daily torpor (Geiser, 2004;Ruf & Geiser, 2015). Hibernators are capable of multiday torpor bouts and drastic Tb reduction, while daily heterotherms exclusively use shorter (<24 hr), shallower torpor bouts (Geiser, 2004;Ruf & Geiser, 2015). However, this binary categorization of heterothermy has been debated, with some authors arguing that hard boundaries between hibernation and daily torpor do not exist (Boyles et al., 2013;Canale, Levesque, & Lovegrove, 2012;Geiser & Ruf, 1995;Ruf & Geiser, 2015;van Breukelen & Martin, 2015). As new methods are developed, support for a continuum of heterothermic abilities has been growing but a consensus has not been reached (Boyles, Bennett, Mohammed, & Alagaili, 2017;Levesque, Nowack, & Stawski, 2016).
Although the common ancestor of mammals was traditionally presumed to be homeothermic (Crompton, Taylor, Jagger, & a, 1978), some recent authors have argued against this interpretation on physiological and/or behavioral grounds (Grigg, Beard, & Augee, 2004;Lovegrove, 2012). Grigg et al (2004) suggested that endothermic heterothermy may be ancestral to homeothermy in mammals because endothermic heterothermy would provide a reasonable intermediate along the transition from ectothermy to full endothermic homeothermy, a shift that would have required a many-fold increase in metabolism (Geiser & Stawski, 2011).
Transitions to new thermoregulatory regimes over evolutionary time require physiological changes. Leptin, a hormone encoded by the lep gene, influences thermoregulation and satiety by signaling to the brain the degree of adiposity available for energy intake (Dodd et al., 2014;Enriori, Sinnayah, & Simonds, 2011;Kaiyala, Ogimoto, & Nelson, 2015;Zhang et al., 2011). Due to the role it plays in regulating satiety and metabolic activity, leptin concentrations may also be impacted by diet (Clarke & Connor, 2014;Harvey et al., 2000;Weigle et al., 2005). Leptin function is correlated with post-hibernation weight gain in arctic ground squirrels (Boyer et al., 1997), prehibernation weight gain in male woodchucks (Concannon, Levac, & Rawson, 2001), and feeding activity of at least one homeothermic bat (Scotophilus heathii; Srivastava & Krishna, 2008). Exon 3 of lep may be particularly important because it contains 29 amino acid variants with functional significance (Yuan et al., 2011). Some hibernating bats show structural alterations in exon 3 of lep that may cause leptin to become more physiologically active (He et al., 2010).
Outside of the relationship leptin has with food intake and temperature regulation, leptin has a mechanistically independent function in initiating torpor bouts in some species (Stehling, Doring, & Ertl, 1996;Swoap, 2008). As part of a complex matrix of gene-gene interactions, sympathetic activation of white adipose tissue during fasting triggers a dramatic decrease in circulating leptin which is thought to initiate a torpor bout (Swoap, 2008). Entrance into torpor for Siberian hamsters is triggered by leptin (Freeman, Lewis, & Kauffman, 2004), and regulation of T B is impacted by leptin in homeothermic rats (Stehling et al., 1996).
Bats (Chiroptera) represent a useful model for comparative evaluation of the evolutionary genomics of thermoregulatory regimes because they are a large monophyletic group with nearly 1,400 extant species, and show considerable variation in thermoregulatory regimes across the phylogeny including homeothermy, daily torpor, and hibernation (Lyman, 1970;Simmons, 2005;Stawski, Willis, & Geiser, 2014). Bats are also among the most diverse mammalian clades with respect to diet, with various lineages specialized as insectivores, carnivores, frugivores, nectarivores, omnivores, and even sanguinivores (Fenton & Simmons., 2015;Hill & Smith., 1984).
To better understand the unique evolutionary history of thermoregulation, we examined the relationship between thermoregulatory regime of extant taxa and selection on lep exon 3. Yuan et al. (2011) found that leptin has undergone positive selection in heterothermic bat lineages, and therefore associated it with the evolution of torpor.
However, by using a lep exon 3 gene tree instead of a species tree, Yuan et al. (2011) effectively treated two nonindependent variables (branch length and lep selection) as independent. Coalescence theory suggests that the use of gene trees to calculate ω, a measure of positive selection, may bias results because gene trees do not necessarily mirror species trees, and may have alternative branching patterns and lengths compared to the species tree (Diekmann & Pereira-Leal, 2016). To reexamine the findings of Yuan et al. (2011), we hypothesized that lep exon 3 has undergone episodic positive selection in heterothermic lineages. Episodic positive selection is defined as positive selection which occurs in a subset of lineages. We TA B L E 1 Overview of the thermoregulatory regimes and diet of bat species included in this study  Bartels, Law, andGeiser (1998), Geiser (2006) Nectarivorous Nowak (1994), Hollar and Springer (1997) AMCC: CEF801

Myotis velifer
Hibernation Tinkle and Patterson (1965), Riedesel and Williams (1976) TA B L E 1 (Continued)  TA B L E 1 (Continued) tested this hypothesis in the context of a species tree for relevant bat taxa. We also reconstructed the ancestral thermoregulatory regimes of bats with revised species-level data compared to previous family-level analyses (Yuan et al., 2011). By evaluating the number of transitions between regimes across the phylogeny through stochastic character mapping, we also evaluated the validity of categorizing thermoregulatory regimes as three discrete character states. Finally, to address other known effects of leptin, we tested for correlation between lineages with positively selected lep exon 3 and the diets of those lineages.

| Taxon sampling
To test for positive selection of lep exon 3, 31 species of bats with known thermoregulatory regimes were sampled for the lep exon 3 gene, including representatives from nine families (  between the dates of September 2015 and August 2018, using keywords "bats," "thermoregulation," "thermoregulatory regimes," "temperature regulation," "hibernation," "torpor," "daily torpor," "daily heterothermy," "heterothermy," "homeothermy," "metabolism," and "Chiroptera." These data are summarized in Table 1 Sequences were assembled and edited within Geneious (Kearse et al., 2012). Ends of forward and reverse sequences were trimmed with an error set to 0.01. Trimmed regions were ignored and forward and reverse sequences were assembled using de novo assembly, as outlined by the Geneious manual (Biomatters, 2017). The forward and reverse sequences of lep for Nyctimene major did not assemble through de novo assembly. Therefore, these were mapped to the reference sequence from which the primers were built. Reads were then manually edited to maximize the coverage and identity be-

| Branch-sites under positive selection
To identify lineages that have experienced episodic positive selection of lep exon 3, we ran a mixed-effects model of evolution (MEME) to detect a subset of branch-sites under episodic positive selection (Murrell et al., 2012). MEME uses a fixed-effect model to explain the distribution and variation of ω across sites, and a random-effects model to explain variation in the distribution of ω across branches (Kosakovsky Pond et al., 2011;Murrell et al., 2012;Nielsen & Yang, 1998). MEME was used here because other tests, which usually average ω over branches and sites, often miss positive selection when it occurs in a subset of branch-sites (Yang & Nielsen, 2002).
When evaluating if a branch-site is under positive selection, MEME considers the specific model of molecular evolution, which here was the TRN93 model (Tamura & Nei, 1993), differences in codon frequencies, and ω.
The test for branch-sites under episodic positive selection, MEME, was performed within the DataMonkey web server (Delport, Poon, Frost, & Kosakovsky Pond, 2010) using the aligned lep exon 3 sequences, the automatic substitution model selection tool, and a user-specified tree from Amador et al. (2016). This 807 taxa phylogeny was calibrated using 44 key fossils, inferred using nine nuclear and mitochondrial genes, and shows support for the majority of currently recognized bat clades (Amador et al., 2016). The Amador et al. (2016) tree was chosen for this study because it represents the most genus-and species-level diversity, 90% and 64%, respectively, compared to other phylogenies (Amador et al., 2016). The tree was pruned using the ape package in R (Paradis, Claude, & Strimmer, 2004) to only include the 31 bats and three outgroup taxa in our dataset. The model of molecular evolution that best fit these data was determined to be the TRN93 model (Tamura & Nei, 1993) using the automatic function available in DataMonkey (Delport et al., 2010).
The significance threshold was set to 0.1 and a log-ratio test  Two analyses to test for gene-wide positive selection were performed. In the first analysis, we tested for selection along branches where evidence for positive selection was previously detected ( Figure 1a; Yuan et al., 2011). In the second analysis, we selected foreground branches based on lineages with evidence for positive selection from the MEME results ( Figure 1b). Here, branches were considered foreground if they met the following criteria based on the MEME output: A branch had at least one codon site with an EBF > 3 for having ω > 1, and a posterior probability >0.25 for having ω > 1. These criteria conform with the guidelines for interpreting EBF values from Kass and Raftery (1995). In addition to the LRT calculation, Akaike information criteria (AIC) scores were calculated.
AIC statistically quantifies the quality of each model by considering the optimum log likelihood (l) and the number of parameters (p) (AIC = −2l + 2p), enabling model comparison.

| Correlation between positive selection and phenotypic traits
To determine whether evolution of lep exon 3 is correlated with evolution of thermoregulatory regimes, we ran a threshold model (Felsenstein, 2012). Here, a discrete trait is presumed to change state when an underlying variable, the liability, crosses a certain threshold (Felsenstein, 2005(Felsenstein, , 2012. This liability is assumed to have a multivariate normal distribution and to evolve under Brownian motion (Felsenstein, 2005(Felsenstein, , 2012. Contrasting the commonly used continuous-time Markov (Mkn) model of discrete character evolution (Lewis, 2001;Pagel, 1994), character states under the threshold model are inherently ordered and the evolution of discrete character states is not memoryless (Felsenstein, 2005(Felsenstein, , 2012; the character state at one node is influenced by the character state at previous nodes. We categorized thermoregulatory regimes into three states: hibernation, daily torpor, and homeothermy. The presence (1) (Revell, 2012). We ran the chain for 3 million generations, thinned to 1,000 samples per chain, to account for autocorrelation, and discarded the first 500,000 as burn-in. To quantify the relationship between lep evolution and TR, a correlation coefficient was calculated between the liability and the thermoregulatory regime. Finally, the highest posterior densities (HPDs) were estimated from the correlation coefficients to determine whether the correlation between traits statistically differed from 0, indicating a statistically significant relationship between the two variables. To test the alternative hypothesis that lep exon 3 evolution is correlated with diet, we used the same methods as above but used diet as the discrete trait evolving under the liability.

| Ancestral state reconstruction of thermoregulatory regimes
To model the ancestral thermoregulatory regimes of bats, we first used the Mkn model of discrete character evolution (Lewis, 2001;Pagel, 1994) with the states hibernation, daily torpor, and homeothermy. It was important to test which transition rate matrix best described the data. Transition matrices describe the rate of transitioning from state i to state j. Here, we tested the data to Given these parameters, we ran an ancestral state reconstruction using the Ace function from the ape package in R ( Figure 2) with maximum likelihood estimation to obtain probabilities of states at interior nodes (Paradis et al., 2004). Stochastic character mapping (Bollback, 2006;Huelsenbeck, Nielsen, & Bollback, 2003) was also used to estimate states at interior nodes and to determine how well the chosen parameters matched the real data. Stochastic character maps were built using the make.simmap function in the phytools package in R (Revell, 2012). Fifty thousand simulations were run using the parameters described previously. A Q-Q plot was generated to compare the Mkn model to the stochastic character map to evaluate the goodness of fit. The stochastic character map was then used to quantify the number of transitions between states across simulations.  F I G U R E 3 Ancestral state reconstruction from stochastic character mapping. Topology of tree generated from Amador et al. (2016). Pie charts represent the posterior probability of each thermoregulatory regime at a given node the other sequences in the dataset. Once aligned, sequences were trimmed to the coding region of lep exon 3 at 357 bp across 34 taxa.

| Measuring positive selection
From the MEME results, evidence for positive selection was detected in 15 branches (Table 2) over a total of 10 codons (Table 3).
When these lineages were selected as the foreground branches in the BUSTED analyses, evidence for episodic positive selection was detected (p = 0.009; Figure 1b). Positive selection was not detected when the foreground branches were selected according to those previously found by Yuan et al. (2011) to be under positive selection (p = 0.220).

| Correlation between positive selection and phenotypic traits
No

| Ancestral state reconstruction
The scaled likelihoods from the Mkn model (Supporting Information   Table S1) and the posterior probabilities from stochastic character mapping were highly correlated (ρ = 0.985, p = 2.2 × 10 −16 ;

| D ISCUSS I ON
We found that the common ancestor of bats most likely used daily heterothermy and is very unlikely to have used hibernation. We also found that leptin evolution is not associated with the evolution of thermoregulatory regimes in bats. Twente and Twente (1964) hypothesized that the most recent common ancestor of bats was homeothermic and that heterothermy evolved secondarily as an adaptation to survive cold climates (Bieber & Ruf, 2009;Kortner & Geiser, 2000). However, daily heterothermy remains adaptive even at warmer temperatures because it increases energy savings and long-term survival (Geiser & Stawski, 2011;Stawski & Geiser, 2012). Flexible use of torpor may have enhanced the ability of some bats to survive in warm and/or tropical climates, such as the environments likely encountered throughout the Eocene when bats likely originated (Amador et al., 2016;Czenze & Dunbar, 2017;Czenze, Brigham, Hickey, & Parsons, 2017b;Meredith et al., 2011;O'Leary et al., 2013;Simmons, 2005;Simmons & Geisler, 1998;Simmons, Seymour, Habersetzer, & Gunnell, 2008;Teeling, 2005).
Supporting our results, recent work estimating the ancestral thermoregulatory regimes of bats (Yuan et al., 2011) and evidence for the commonality of heterothermy in bats (Geiser & Stawski, 2011) suggests that heterothermy was the ancestral state for Chiroptera and that homeothermy was secondarily derived (Geiser & Stawski, 2011;Yuan et al., 2011). This scenario mirrors the hypothesis that the common ancestor of all mammals was heterothermic (Grigg et al., 2004;Lovegrove, 2012). However, until now, the question of which heterothermic regime was used by the ancestor of bats-hibernation or daily heterothermy-was unresolved. Here, we show evidence that the ancestor of bats was likely a daily heterotherm.
Consistent with the lack of evolutionary advantages that a hibernator would be expected to accrue during the early Eocene, a relatively warm time period characterized by widespread tropical and subtropical conditions, our results suggest that the common ancestor of bats did not hibernate (Humphries, Thomas, & Speakman, 2002). Although our results suggested a marginal likelihood that the ancestor of bats was a homeotherm, this seems unlikely based on previous studies (e.g., Geiser & Stawski, 2011;Yuan et al., 2011).
Taken together with the high likelihood of daily heterothermy at this node, we argue that the most recent common ancestor of bats was a daily heterotherm. Our reconstruction therefore also indi- Our analyses also suggest that the MRCA of Rhinolophidae used hibernation; however, it is unclear if this is derived or ancestral.
Conservatively, we suggest that hibernation arose at least once in this group. Similarly, we suggest that hibernation arose at least once in the clade comprised of Vespertilionids and Molossids. We found no evidence for reversals in the hibernation phenotype-no lineages that lost and subsequently regained the ability to hibernate. bone development (Crespi & Denver, 2016), induction of mitosis (Gat-Yablonski & Phillip, 2008), and immune and stress responses (Ahima & Osei, 2004;Procaccini, Lourenco, Matarese, & La, 2009).
The potential for positive selection of leptin for alternative traits (Carey, Andrews, & Martin, 2003;Jastroch et al., 2016;Yang et al., 2008) makes it difficult to find correlations between positive selection of lep and a singular function. Therefore, in retrospect it is perhaps not surprising that we found no evidence for a tight correlation between leptin and chiropteran thermoregulatory regimes, and also found no correlation between leptin selection and diet despite the known influence of high-fat diets on leptin functioning (Frederich et al., 1995;Koch et al., 2014).
Complex genomic mechanisms and associated physiological alterations to organ systems and physiological functions across species with different thermoregulatory regimes suggest that the evolution of thermoregulatory regimes may intrinsically have no correlation with positive selection of a single gene (Andrews, 2004;Grabek et al., 2011;Hindle, Grabek, & Epperson, 2014;Morin & Storey, 2009;Villanueva-Cañas, Faherty, Yoder, & Albà, 2014). A non-synonymous substitution at codon site 91 in exon 3, found here to be under positive selection, was previously inferred to cause a functional difference in hibernating bats compared to homeothermic bats (He et al., 2010). However, we found that this substitution in the sequence of the hibernating bat is shared with Homo sapiens, a homeothermic species. Therefore, a direct relationship to thermoregulatory regime and sequence variation cannot be made for this site. Recent evidence suggests that thermoregulatory regimes are mostly influenced by the regulation of gene expression, rather than the sequence specificity of protein-coding genes (Geiser & Stawski, 2011;Grabek, Martin, & Hindle, 2015;Morin & Storey, 2009;Schwartz, Hampton, & Andrews, 2013;Yan, Barnes, Kohl, & Marr, 2008).
In our ancestral state reconstruction, zero direct transitions occurred between hibernation and homeothermy. In this model, we assumed a priori that character states are not ordered. Therefore, the lack of transitions between hibernation and homeothermy is not an artifact of the model. This suggests that thermoregulatory regimes may be better represented as an ordered trait (with daily torpor as a necessary intermediate between homeothermy and hibernation) rather than as an unordered trait. Recent evidence suggests that heterothermy exists along a continuum (Boyles et al., 2013;Dunbar & Brigham, 2010;Lovegrove, 2012;Wilz & Heldmaier, 2000). If true, this suggests that there is an inherent order to the evolution of thermoregulatory regimes, which our results indicate. Our results suggest that future research on mammalian thermoregulation should F I G U R E 4 Pteronotus parnellii photograph taken by Brock and Sherri Fenton in a cave in the western end of Cuba treat thermoregulatory regimes as an ordered and possibly continuous trait.
The genomics underlying thermoregulation in mammals remains largely unclear. Future research should aim to sequence whole genomes of mammals that vary in thermoregulatory regime.
Comparing these data in a phylogenetic framework would enable a more complete understanding of the genomic components involved in the evolution of thermoregulatory regimes. Our results revealed that leptin does not appear to be directly involved in the evolution of thermoregulatory regimes but many candidate genes including LEPR (Rezai-Zadeh et al., 2014), MEF2 (Tessier & Storey, 2010), and G0S2 (Jessen et al., 2016) have yet to be examined in a similar framework.
Such studies will reveal the importance these candidate genes have in the evolution of thermoregulation across diverse taxa.

CO N FLI C T O F I NTE R E S T
None declared.

AUTH O R CO NTR I B UTI O N S
M.E.L, F.T.B., and N.B.S. contributed to the concept and design of project. M.E.L. performed laboratory work, data collection, and analyses, and drafted the manuscript. All authors contributed to manuscript preparation and approval of final draft.

DATA ACCE SS I B I LIT Y
Sequence data are available on GENBANK (https://www.ncbi.nlm. nih.gov/genbank/) with accession numbers listed in Table 1