Assigning metabolic rate measurements to torpor and euthermy in heterothermic endotherms: ‘torpor’, a new package for R

Torpor is a state of controlled reduction of metabolic rate (M) and body temperature (T_b) observed in numerous mammals and birds (Ruf and Geiser, 2015). This physiological state can occur over short periods (i.e. <24 h), often referred to as ‘daily torpor’, or it can last up to many days, often referred to as ‘hibernation’, which typically involves multiday torpor bouts separated by short spontaneous arousals (McKechnie and Lovegrove, 2002; Ruf and Geiser, 2015; Schleucher and Withers, 2001). The drop in T_b may be considerable in hibernators, which often allow T_b to reach values close to 0°C, but it can also be much shallower, as is often the case in species entering daily torpor (Ruf and Geiser, 2015). It has been suggested that a continuum may exist between these types of torpor as well as between torpor and euthermic rest (Boyles et al., 2013; McKechnie and Lovegrove, 2002; van Breukelen and Martin, 2015). Torpor use has profound implications on energy expenditure and allocation (Lyman et al., 1982) and affects many biological functions (Geiser and Brigham, 2012; Nowack et al., 2017). Accordingly, a vast amount of literature describes the occurrence of torpor and its associated energy savings in a multitude of species (Geiser and Ruf, 1995; Lovegrove, 2012; Nowack et al., 2020; Ruf and Geiser, 2015).

Determining whether measurements of M can be assigned to torpor or euthermy can however be challenging. The distinction between these two states is usually straightforward in typical hibernators, due to the more than 90% reduction in M common in these species when torpid at their usually hibernation ambient temperature (T_a) (Geiser, 2004). Yet this distinction can become problematic when the difference between euthermic M (M_e) and torpid M (M_t) is small, which often occurs slightly below the thermoneutral zone (TNZ) (Geiser, 2011; Hainsworth and Wolf, 1970; Humphries et al., 2002; Speakman and Selman, 2003). Additional difficulties are encountered with species that either enter daily torpor, referred to as ‘daily heterotherms’, whose minimum T_b in torpor often lies only slightly below the euthermic T_b (Bartels et al., 1998; Bonaccorso and McNab, 1997; Genoud et al., 1990; McNab, 1980b), or who exhibit a large variability in M and T_b at rest (Coburn and Geiser, 1998; Geiser et al., 1996).

Barclay et al. (2001) reviewed the variety of criteria that have been chosen in the past to distinguish torpor from euthermy. Apart from a minority of reports that identified torpor on the basis of behavioral features (Brice et al., 2002; Geiser and Kenagy, 1988; Geiser and Masters, 1994), the vast majority of studies have used patterns of variation in T_b (or skin temperature) and/or M to separate the two states in the field or in the laboratory. Animals have been deemed to be in torpor below a threshold T_b or M (Coburn and Geiser, 1998; Geiser et al., 1996; Hosken and Withers, 1999; Kelm and von Helversen, 2007; Levesque, 2008; Mzilikazi and Lovegrove, 2002), below some threshold temperature differential between the body and air (Levesque and Lovegrove, 2014), or below a threshold percentage of the euthermic rate of metabolism (Geiser, 1988a; Hudson and Scott, 1979). This diagnostic threshold value was sometimes calculated on the basis of the parameter's variation in the euthermic state (e.g. Lovegrove and Raman, 1998; McKechnie et al., 2007), but has also been predicted by an equation based on body mass and T_a (Willis, 2007). All available techniques have their limitations and are at least partly arbitrary (Barclay et al., 2001; Boyles et al., 2011), despite the efforts made to render them more objective.

With this paper, we propose an automated method for distinguishing torpor from euthermic rest in species that enter torpor, hereafter referred to as heterotherms following the terminology of Geiser and Ruf (1995) or Lovegrove (2012). Here, assignment to an either euthermic or torpid state is based on a probabilistic approach using the variation observed among measurements of M at varying T_a. We assume that the relationship between M of resting animals and T_a follows the classical ‘Scholander-Irving model’ (Scholander et al., 1950; for a discussion see McNab, 2002), which was later extended to include torpor (Hainsworth and Wolf, 1970; Humphries et al., 2002; Speakman and Thomas, 2003; Geiser, 2011). This model (Fig. 1) predicts the vast majority of patterns observed among endotherms entering daily torpor or hibernation. We explicitly consider that additional metabolic inhibition involving mechanisms other than the abolishment or reduction of thermogenesis necessary to maintain euthermy may occur during torpor (Geiser, 1988b; 2004; Geiser and Kenagy, 1988; Guppy and Withers, 1999; Withers et al., 2016), by allowing the curve of torpor M to reach T_lc at a level equal to – or lower than the resting metabolic rate within the TNZ. Thus, over all aforementioned techniques, our novel method provides major improvements, as it assigns measurements under an explicit framework (i.e. the extended Scholander-Irving model) and does so with a probabilistic approach based on observed variation.

Our method is specifically intended to facilitate the discrimination between torpor and euthermy in laboratory experiments using respirometry. Additionally, it can also be useful to predict M at varying T_a, by modelling the different parameters of the thermoregulatory curves. Further, we have assessed the performance of our method, by applying it to previously published data on the M of various mammals, including both heterotherms and species not undergoing torpor (hereafter referred to as homeotherms). Finally, we have provided a new package running within R (R Development Core Team, 2012) to allow researchers to apply our method to their own data. This package (‘torpor’) comprises several useful functions that will improve standardization of the analyses of metabolic measurements for thermal biology.

Over all datasets, the proportion of assignments, that could be validly assigned (i.e. assignment confidence >0.80, cf. ‘Assignment confidence’ in section ‘Materials and Methods’) ranged from 0.63 to 1.00 (median=0.87, Table S1). The corroboration index, which assesses the similarity between valid assignments made by our method and by the authors of the original datasets, ranged between 0.66 and 1.00 (median=1.00, Table S1). Complete matches between the method valid assignments and the authors’ descriptions were found for 15 out 28 datasets (e.g. Nyctophilus geoffroyi; Hosken and Withers, 1999; Fig. 2A). Mismatches mostly occurred close to T_lc (e.g. Peropteryx macrotis; Genoud et al., 1990; Fig. 2C) but could also be found as far as 41.0°C below T_lc (Fig. 3). Among homeotherms, no M measurement was assigned to torpor (e.g. Sorex minutus; Sparti and Genoud, 1989; Fig. 2B).

In studies of homeotherms, where no M was originally assigned to torpor, prior and posterior distribution overlaps (PPOs) of TMR ranged from 76.80 to 95.60%. In identified heterotherms, PPOs of TMR values range from 4.60 to 64.60% (Table S1).

Modelled values of T_lc, M_TNZ and T_t and extracted values from the original studies were significantly correlated (T_lc: N=28, Pearson's coefficient: 0.65, P<0.001, M_TNZ: N=28, Pearson's coefficient: 1.00, P<0.001, T_t: N=14, Pearson's coefficient: 0.97, P<0.001).

Our method aims to facilitate the discrimination between torpor and euthermy using the variation of M measured at varying T_a. We tested it using 28 published sets of metabolic data obtained by respirometry on 26 species of small or medium-sized mammals. Selected species displayed a diversity of metabolic and thermal strategies ranging from permanent homeothermy to heterothermy including shallow, daily torpor and deep, long-term hibernation. The efficiency of our method proved satisfactory. Indeed, the corroboration index was generally high. Most conflicting assignments were a mere consequence of the difference between the estimated T_lc and that defined by the authors (Table S1; Fig. 2C and Fig. 3). In particular, all mismatches between the model's and the authors’ assignments concerning M_e and M_TNZ measurements in homeotherm species were explained by differences in estimated T_lc. The remaining mismatches could generally be explained by intra-state frequency distributions of M that did not segregate clearly. This condition was often found close to TNZ, but also occurred in species where T_t was relatively close to T_lc (e.g. McNab, 1980a, 1989). Such discrepancies illustrate the difficulty of the assignment process when torpor and euthermy have to be distinguished solely on the basis of measurements of M. While the use of a statistical method enables an objective assignment, it is nevertheless worthwhile to recall that some authors also used patterns of T_b to assign their data (e.g. McNab, 1980a, 1989). These specific cases should deserve further attention as they can highlight mechanisms decoupling M from the control of T_b (e.g. Daniels, 1984; McNab, 1988; Heldmaier et al., 2004). One advantage of our method is to allow researchers to identify statistically data points that are difficult to assign. In this study, we considered assignments with a confidence lower than 0.8 as invalid. Obviously, an increase of this threshold value would lead to fewer assignment mismatches.

Our method also models the relationship between M and T_a specific to each state, as well as several parameters describing the standard energetics of the animal(s) under study, including M_TNZ, TMR, T_lc and T_t. Predicting M at any T_a is crucial to model energy costs (e.g. Boyles et al., 2020). It should, however, be recalled that our Bayesian inference-based method models all parameters whatever the available data. It remains therefore crucial to consider the parameters’ identifiability. The PPO of T_bt and T_lc are not reliable for that purpose, because the prior distribution of those parameters largely depends on the data provided. From the PPOs of TMR obtained in studies where no M was assigned to torpor (i.e. ‘homeotherms’) and those obtained in identified heterotherms, we define a PPO higher than 75% as indicative of a parameter that is not identifiable.

The parameters estimated by our method correlated significantly with those extracted from the original studies. The number of M_e values provided was probably not always sufficient to model T_lc adequately. On one hand, a statistical estimation of T_lc (e.g. in Song et al., 1995, 1997; Bonaccorso and McNab, 1997; Genoud et al., 1997; Hosken and Withers, 1999; Willis et al., 2005; Genoud, 2014) was often lacking, leaving space for some potential misleading subjectivity by the authors of the original studies. On the other hand, an overestimation of T_lc by our method remains possible if appreciable additional metabolic inhibition occurs during torpor (see below). The estimated T_t was always close to that of the original studies. This higher precision in the estimation of T_t in comparison to that of T_lc is probably due to the lower variation of M exhibited during conforming torpor in comparison to that measured during euthermic rest. It is also not surprising that the M_TNZ obtained by our method and by the authors of the original studies were close to identical, since both were calculated as the mean of nearly the same M values (i.e. all values above T_lc).

An important aspect incorporated in our method is that we allow function for M_t to be lower than M_TNZ at T_lc. We thus explicitly allow for additional metabolic inhibition to occur in some species during torpor, other than the cessation or reduction of the thermogenesis necessary to maintain euthermy (Geiser, 2004; Guppy and Withers, 1999; Withers et al., 2016). Caution should be kept in M_r interpretation as an insufficient number of M_t values corresponding to conforming torpor leads to the unidentifiability of this parameter (e.g. Melonycteris melanops; Bonaccorso and McNab, 1997; Fig. 2D). Consequently, the mentioned additional metabolic inhibition should only be considered when sufficient M_t values corresponding to conforming torpor are provided and M_r PPO indicates adequate parameter's identifiability. Considering only identifiable parameters (i.e. M_r PPO<75%), we quantified additional metabolic inhibition as the fraction M_r/M_TNZ. Values ranged between 0.27 and 0.93 (N=18, median=0.77; Table S1). By extrapolation, we should consider the possibility that a hypometabolic state might occur within the TNZ as well (e.g. Grimpo et al., 2013; Reher et al., 2018). With our method, we did not implement an assignment to such a hypometabolic state within the TNZ, hence M_TNZ corresponds to the mean of all values from the TNZ. Our method might therefore underestimate the eumetabolic M_TNZ, and might consequently overestimate T_lc, if animals occasionally enter a hypometabolic state in the TNZ, A predicted M_t at T_lc much lower than M_TNZ should be viewed as an indication for a strong additional metabolic inhibition occurring during torpor and should invite for a careful consideration of the estimated M_TNZ and T_lc.

Throughout this paper, resting metabolic rate within the TNZ is referred to as M_TNZ rather than to the basal rate of metabolism (BMR). Estimating the BMR requires several specific criteria: animals should be post-absorptive, adult, non-reproductive, and resting during a major inactive phase of the daily cycle (McNab, 1997). If these criteria are not met, the method still remains applicable, but the resting M estimated within the TNZ (i.e. M_TNZ) will correspond to a minimal resting metabolic rate rather than to the BMR.

While the performance of our method appeared to be excellent with the selected datasets, its main limitations correspond to those of the model on which it is based. Indeed, although extremely fruitful, the Scholander-Irving model (Scholander et al., 1950) and its latter extensions (Geiser, 2011; Hainsworth and Wolf, 1970; Humphries et al., 2002; Speakman and Thomas, 2003) does not describe the thermal biology of all heterothermic species. Causes for significant divergence from the hereby model include a non-linear relationship between M_e and T_a (e.g. Cercartetus nanus: Song, Körtner, and Geiser, 1997), a difference in the slopes of the regressions of M versus T_a in regulated torpor and euthermic rest (Geiser, 2004) and a strong dependence of M_t on factors other than T_a, such as body mass (Kelm and von Helversen, 2007) and duration since the last activity period or since last meal (Grigg et al., 1992; Morris et al., 1994). Furthermore, this model also does not describe adequately circadian (e.g. the shallow ‘rest-phase hypothermy’ of many small birds; McKechnie and Lovegrove, 2002; or ultradian; Heldmaier et al., 1989) variations in M.

Our method is able to reveal the presence or the absence of distinct groups of values (euthermic rest versus torpor) but is not intended to judge whether any of the data points entered correspond to stable rates and/or to minimal values. Transitions between activity and rest or between euthermic rest and torpor as well as incomplete torpor bouts during which M never reaches a stable, minimal level (e.g. ‘test drops’; Hudson and Scott, 1979; see also Genoud, 1993; Lyman, 1982) should be excluded. It should also be recalled that long runs may be necessary for some heterothermic animals to achieve a state of steady torpor (e.g. Heldmaier, Ortmann, and Elvert, 2004). The above limitations should be kept in mind for a lucid application.

Then, individual variation is not explicitly considered in our method. Especially, body mass variation within a studied population might affect the assignment of M to torpor or euthermy in a complex manner (Genoud, 2014; Genoud et al., 2018; Glazier, 2005; Nespolo et al., 2003; Sassi and Novillo, 2015; Schleucher and Withers, 2001). However, there would be no particular difficulty to apply our method to single individuals, or repeatedly on single individuals.

Our method and the associated R-package (torpor) provide a way to standardize the analysis of respirometry data in relation to T_a. Its major strength is that it uses a probabilistic approach to assign metabolic values to torpor or euthermy, rather than assigning them on the basis of a particular threshold value. Thus, the partly arbitrary nature of the assignment process is removed. Further, it can be applied to study intraspecific as well as interspecific variation in energetics. Parameters might now be extracted from the literature for comparative analysis, largely avoiding causes of variation due to the diversity of past assignment techniques. Ultimately, this method and the associated R-package (torpor) will ease intra- and inter-specific comparative analyses of endotherm energetics.

Our present method is based on three steps and requires measurements of M of resting or torpid animals at T_a ranging from below – to within the TNZ. Animals are assumed to be in one of three ‘states’: euthermic rest below the TNZ (M_e), torpor below the TNZ (M_t) or rest within the TNZ (M_TNZ). We do not consider here the relation between M and T_a above the TNZ. Consequently, metabolic rates measured at T_a higher than the upper limit of the TNZ should be excluded. During the different steps, a mixture model based on Bayesian inference is run under varying conditions.

The method initially defines the highest possible T_a within the dataset that still underestimates T_lc (T_{lc_low}). Above T_{lc_low}, a linear regression between M and T_a should neither result in a significant negative slope nor be affected by heteroscedasticity. Beginning with the ten M values measured at the highest T_a, linear regressions are performed on sets of M values progressively including values at lower T_a. The inclusion of additional M values obtained below T_lc eventually leads to a significantly negative slope in the case of animals not entering torpor and/or to heteroscedasticity in the case of heterotherms. A significantly positive slope would reveal an absence or insufficient number of values for euthermic rest below the TNZ, hence it automatically stops the analysis as T_{lc_low} is undefinable. Heteroscedasticity is assessed with a Breusch-Pagan test with the function ‘bptest’ from package ‘lmtest’ (Zeileis and Hothorn, 2002). Significance level for the Breusch-Pagan test is set at 0.05. The significance of the negative regression is assessed with a linear regression, function ‘lm’ and significance level is set at 0.01.

Then, in order to get T_lc, the model described in the second step [see ‘Second step: M measurements pre-assignment’] is first run without data points measured above T_{lc_low}. Moreover, within this step, the parameter T_lc is also modelled. M_TNZ provided in that analysis is the mean of the M measured at T_a higher than T_{lc_low}. Once T_lc has been estimated, M_TNZ is recalculated as the mean of the M measured at T_a higher than T_lc.

The assignment of M measurements to one of the three physiological states is based on a probabilistic process. Below T_lc, a posterior categorical distribution is generated for each coupled data (T_a and M), providing the probabilities to belong either to torpor or to euthermic rest (i.e. value specific state membership probabilities). Independent of the pre-assignment, M values higher than the predicted M_e or lower than the predicted M_t are automatically assigned to M_e and M_t, respectively. This latter automatic procedure is, however, only performed if at least one M measurement has been assigned to torpor during that second step. As defined previously, M values located above T_lc are automatically assigned to M_TNZ. Automatically assigned values get a state membership probability of one.

In the third step, only M values situated between the predicted M_e and M_t are assigned. Measurements’ states that are automatically assigned during the second step are provided along with the coupled M and T_a values. During that final step, the parameters of the functions relating M to T_a below the T_lc, as well as the standard deviation for the different physiological states are modelled.

Three different Markov chains are run during 50,000 iterations starting at initial values within the range of parameter space. The initial convergence phase is excluded by dropping the first 30,000 iterations. Markov chains are thinned by a factor of 10 and the Brooks–Gelman–Rubin criterion R̂ (Brooks and Gelman, 1998) is used to assess the convergence of chains (<1.1). To ease parameters’ estimation, metabolic rate measurements are divided by their mean.

The PPO enables the evaluation of the relevance of some modelled parameters. Overlap values are obtained with the function ‘MCMCtrace’ from package ‘MCMCvis’ (Youngflesh, 2018). Truncated Gaussian prior distributions are provided with the function ‘rtruncnorm’ from package ‘truncnorm’ (Mersmann et al., 2018).

The M value is assigned to the physiological state having the highest specific state membership probability. The assignment confidence represents the product of the highest specific state membership probability and of the probability that T_a is above (for M_TNZ) or below (for M_t and M_e) T_lc. That latter probability is calculated from the variation in the estimation of T_lc (cf. ‘First step: estimation of M_TNZ and T_lc’), as the value at T_a of the T_lc cumulative posterior distribution function, or 1 minus that function, respectively. A threshold proportion is then selected to consider any assignment with a significantly lower assignment confidence as invalid. That hypothesis is tested with a binomial test, function: ‘binom.test’, with the alternative hypothesis set as ‘greater’. The significance level for the binomial test is set at 0.05. For the present study, assignments with a confidence higher than 0.8 were considered as valid.

We evaluated our method by applying it to 28 published sets of M data obtained by respirometry. Mainly, we aimed at highlighting strengths and possible issues linked to the application of our method to typical sets of measurements. Thus, we were specifically interested in divergences between assignments made by our method and by the authors. The 26 species of small or medium-sized mammals investigated illustrate a diversity of metabolic and thermal strategies ranging from permanent homeothermy to heterothermy including shallow, daily torpor and deep, long-term hibernation. The selected datasets also differ in size (Table S1). For each published set, we considered all metabolic data provided, except those made above the upper critical temperature defined by the authors (see ‘the three-steps method’ above). In 20 species (here considered ‘heterotherms’, Table S1), rates of metabolism measured at T_a’s below the described thermoneutral zone were originally assigned to either of two states, which were usually referred to as euthermy (or normothermy) and torpor. In Cynopterus brachyotis, values corresponding to particularly low T_b were identified but not explicitly referred to as torpor (McNab, 1989). Some of the values obtained from torpid Lasiurus seminolus by Genoud (1993) were originally characterized by ‘irregular fluctuations’ in M (Genoud, 1993). For simplicity, we treated these values as if the authors assigned them to torpor. The remaining six species are hereafter referred to as homeotherms (Table S1), as none of the values in the corresponding datasets were assigned to torpor in the original studies. Chosen studies reported M_TNZ (mostly referred to as BMR) and T_lc estimates, although the latter had to be extracted from a graph in two cases or calculated as the intercept between the euthermic line below thermal neutrality and M_TNZ provided by the authors in another case (Table S1). In all dataset the authors provided a sufficiently precise graph of metabolic rate values as a function of T_a and several of them also estimated T_t. Plotted data were digitized using the software Plot Digitizer (Huwaldt and Steinhorst, 2015) except for the Pipistrellus pipistrellus dataset, which was provided by MG. Sample sizes of the digitized data ranged between 27 and 204 (Table S1). We tested the correlation between the M_TNZ, T_lc and T_t values modelled with the present method and those provided by the authors with Pearson's paired-samples correlation test, function ‘cor.test’. Cases in which the estimation of T_t fell outside the range of measured T_a were excluded from this analysis as T_t was considered unidentifiable. Then, we examined whether the assignment of the metabolic values to torpor or euthermy made by the authors and provided by the method coincided, and we calculated a corroboration index as the fraction of the matched assignments. Only measures with valid assignments (i.e. assignment confidence >0.80, cf. ‘Assignment confidence’ in previous section) were considered for the correlation tests and the corroboration index calculations. Finally, in order to define a PPO range highlighting an identifiable parameter, the PPOs of TMR from studies where part of the M measurements were assigned to torpor were compared with those from studies where no M measurement was assigned to torpor. In those later studies TMR was considered unidentifiable.

We are thankful to Philippe Christe and Nicolas Salamin for their support during the elaboration of the model. We also thank four anonymous reviewers who contributed to improve the original manuscript.