Heterothermy in pouched mammals – a review

Database

Data from published studies on 41 marsupial species (Australidelphia and Ameridelphia) were available for the present analysis (Table 1). Species ranged from small marsupials such as the narrow-nosed planigale Planigale tenuirostris weighing 7.0 g to larger species such as the Western quoll Dasyurus geoffroii and the Eastern quoll Dasyurus viverrinus weighing 1000 g (Table 1). We aimed to include only heterothermic marsupial species for which body mass and at least one heterothermic variable, that is T_{b min} or TMR, has been reported. Marsupial species that are likely to be heterothermic, for example by phylogenetic relationships to known heterothermic species, have not been included in our analysis. Thus, our dataset comprised 41 heterothermic marsupial species with data for body mass, 38 for normothermic body temperature (T_b), 38 for T_{b min}, 30 for TMR, 30 for maximum TBD (TBD_max) and 23 for mean TBD (TBD_mean). The T_b and BMR values were taken either from the original publications or from Clarke, Rothery & Isaac (2010) and White, Blackburn & Seymour (2009), respectively. Torpor was defined as a reduction in metabolic rate below 75% of resting metabolic rate (RMR) and TBD as the time metabolic rate remained below 75% of RMR (Geiser et al., 2006). Or, if only T_b measurements were available, as the time T_b remained below 30°C (Körtner & Geiser, 2000). TMR was defined as the lowest metabolic rate during torpor when TMR was stable and minimal. Data for hibernators and daily heterotherms were analysed separately because most physiological variables related to heterothermy, such as TBD_mean, TBD_max, T_{b min} and TMR, differ substantially between these two groups.

Phylogeny

The phylogeny was derived from a published mammalian supertree, which includes 4510 species with updated branch lengths derived from dated estimates of divergence times (Bininda-Emonds et al., 2008). The supertree for mammals was pruned to include only the species of concern for this study, that is heterothermic marsupials (n = 41), using the ‘analysis in phylogenetics and evolution’ (APE) package (Paradis, Claude & Strimmer, 2004) and the ‘analysis of evolutionary diversification’ (GEIGER) package (Harmon et al., 2008) in R (Ihaka & Gentleman, 1996). The tree is available in Supporting Information Appendix S1.

Phylogenetic signal and analysis

The method of phylogenetic generalized least squares (PGLS; Grafen, 1989; Martins & Hansen, 1997; Garland & Ives, 2000; Rohlf, 2001) was implemented for the trait data using BayesTraits (Pagel, Meade & Barker, 2004). PGLS is an extension of Felsenstein's phylogenetic independent contrasts (PIC; Felsenstein, 1985) that allows for flexibility in the underlying evolutionary assumptions (Martins & Hansen, 1997). While PIC analysis assumes a ‘Brownian motion’ model of evolution where the variation in the trait data is fully explained by the phylogeny, PGLS analysis allows more flexibility through the use of a parameter lambda (λ). The parameter λ is determined by maximum likelihood (ML) and can range between 0 [no phylogenetic signal, similar to ordinary least squares, i.e. ordinary least square (OLS) analysis] and 1 (pattern of trait data variation is fully explained by the phylogeny, similar to PIC analysis), and thus indicates how strong the phylogenetic signal is for a certain trait or the relationship between the two traits. Intermediate values of λ indicate that the trait evolution is phylogenetically correlated, but does not follow fully a Brownian motion model (White et al., 2009). A more in-depth description and further mathematical details on PGLS analysis can be found in Pagel (1999), Garland & Ives (2000) and Freckleton, Harvey & Pagel (2002).

To compare phylogenetic signals for individual traits for measured trait data in daily heterotherms, we compared λ = 0 and λ = 1 models to the PGLS model (i.e. λ = ML) using BayesTraits (Pagel et al., 2004). A likelihood ratio (LR) test was used, similar to that described in Capellini, Venditti & Barton (2010), with LR = 2 × [Lh (better-fitting model) − Lh (worse-fitting model) ] with the best-fitting model having the highest log-likelihood score (Lh). The significance was tested with a chi-square distribution with 1 degree of freedom (Pagel, 1999; Freckleton et al., 2002).

Estimation of missing values using phylogeny

For species without measured values for T_b (n = 4), T_{b min} (n = 3), TMR (n = 11), TBD_max (n = 11), TBD_mean (n = 18) and/or BMR (n = 5), we used the program PhyloPars (Bruggeman, Heringa & Brandt, 2009) to estimate these values using the marsupial phylogeny. The program operates on a known phylogeny and an incomplete feature matrix that includes the available observations. A detailed description of the web-based program and its application was published recently (Riek & Bruggeman, 2013). Missing values are estimated using phylogenetic and phenotypic covariances. In a first step, ‘the optimal phylogenetic and phenotypic covariance are combined with the tree topology to calculate the covariances between the observations and the missing values. These are subsequently used to express the estimate of each missing value as the product of all original observations and an estimate specific set of associated weights’ (Bruggeman et al., 2009) forthwith termed ‘evolutionary model’. The estimated phylogenetic covariances, which characterize the rate and direction of evolution, are then combined with the dataset to infer missing trait values. The resulting estimates are the weighted sum of known values for both traits, determined by correlations between species (i.e. phylogenetic distance) and correlations between traits (i.e. allometric relationships). Additionally, to determine the accuracy of the estimated trait values, the model results are validated by cross-validation in which each observed parameter (i.e. measured values) is excluded from the input data and re-estimated using the ML estimation-derived phylogenetic and phenotypic covariances to determine the prediction error. In other words, known trait values are omitted one by one from the dataset and then re-estimated from the remaining values using phylogeny. This procedure results in a set of residuals (one for each known trait value), which determine the quality of the estimated missing values. Furthermore, the mean of the absolute values of these residuals is used as a measure of accuracy (Riek & Bruggeman, 2013). For comparison purposes, the program also estimates the error for a simple mean model, which assumes the best estimate for a missing value is given by the mean of all observations, and for a nearest-neighbour model, which assumes the best estimate for a missing value is given by the phylogenetically closest observation (Bruggeman et al., 2009).

Phylogenetic signals for individual traits

For marsupial daily heterotherms, individually analysed traits that were measured exhibited low (TBD_mean, TBD_max), moderate (T_b) and strong (body mass, T_{b min}, TMR, BMR) phylogenetic signals, with a ML of λ values ranging from 0.34 to 1.00 (Table 2). All phylogenetic signals for individual traits in daily heterotherms differed significantly (P < 0.05) from λ = 0, confirming the need for phylogenetically informed analysis. An exception was TBD_mean, where λ = 0.34 was not significantly different from λ = 0 with P being slightly above the significance level (P = 0.051). While the phylogenetic signal for body mass and T_{b min} were not significantly different from λ = 1, the phylogenetic signals for the traits BMR, TMR, T_b, TBD_max and TBD_mean were significantly different from λ = 1, suggesting that a fraction of the variability is not due to evolutionary constraints in these traits. Contrarily, all phylogenetic signals (λ) for all analysed traits in marsupial hibernators were 0, that is no phylogenetic signal (results not shown).

Allometric relationships

We analysed allometric relationships between body mass and TMR (Fig. 2) as well as between TBD_mean and T_{b min} and TMR (Fig. 3) and between T_{b min} and the ratio of TMR/BMR (Fig. 4) separately for observed values (Table 3) as well as for observed and estimated values (Table 4) using the marsupial phylogeny in daily heterotherms and hibernators. For daily heterotherms, all allometric relationships were significant (P < 0.05) except for the relationship of measured as well as measured and estimated values combined between TBD_mean and T_b (Tables 3 and 4). The phylogenetic signals for relationships between traits in daily heterotherms ranged between 0.42 and 0.70 for measured values and between 0.49 and 1.00 for measured and estimated values, and differed significantly (P < 0.05) from 0 in all cases, emphasizing the need for phylogenetic correction of the trait data and the inappropriateness of the use of OLS for these data. In hibernators, however, the phylogenetic signal of only one allometric relationship (TBD_mean vs. T_{b min}) was significantly different from 0 (P < 0.05) for measured values as well as for measured and estimated values; for all other analysed relationships, the phylogenetic signal was 0, that is equivalent to OLS (Tables 3 and 4).

The phylogenetically controlled allometric relationship between body mass and TMR (Fig. 2) for measured values was described by the equations TMR (mL O₂ h⁻¹) = 0.13 body mass^0.57 (g) for hibernators and TMR (mL O₂ h⁻¹) = 0.53 body mass^0.77 (g) for daily heterotherms, with a significant (P < 0.01) lower exponent for hibernators than for daily heterotherms and body mass explaining 55 and 61% of the variation in TMR, respectively. Including estimated values into the regression analysis for daily heterotherms, yielded the equation TMR (mL O₂ h⁻¹) = 0.42 body mass^0.89 (g) with body mass explaining 84% of the variation in TMR. Furthermore, our results show that normothermic T_b increases with increasing TBD_mean in hibernators, but this relationship is not significant in daily heterotherms for both measured as well as measured and estimated values (Tables 3 and 4). Contrarily, T_{b min} decreases with increasing TBD_mean in both daily heterotherms and hibernators, whether measured or measured and estimated data were analysed. However, the regression for measured data in hibernators was not significant (Fig. 3a). The exponents for the corresponding PGLS models were lower for hibernators (measured: −0.45, measured and estimated: −0.89) than for daily heterotherms (measured: −0.27, measured and estimated: −0.23). Similarly, TMR decreases with increasing TBD_mean in both hibernators and daily heterotherms for measured as well as for measured and estimated data (Fig. 3b). In hibernators, however, this decrease is more pronounced as indicated by a steeper slope (measured: −0.52, measured and estimated: −0.62) compared with daily heterotherms (measured: −0.40, measured and estimated: −0.36). The T_{b min} had a significant (P < 0.001) influence on the ratio TMR/BMR (Fig. 4a) in daily heterotherms but not in hibernators, whether measured or measured and estimated values were included in the analysis.

Estimation of missing values using phylogeny

In our analysis, we estimated values for T_b (n = 4), T_{b min} (n = 3), TMR (n = 11), TBD_max (n = 11), TBD_mean (n = 18) and BMR (n = 5) for several marsupial species on the basis of their body mass and phylogeny, for which trait data were not available (Table 1), with species ranging in size from 7 g (Planigale spp.) to 1.0 kg (Dasyurus spp.). An exception was the marsupial mole Notoryctes caurinus for which TBD_mean and TBD_max could not be estimated because of the low number of species in this order (Notoryctemorphia, n = 2) and thus an undetermined model. However, for all other estimated values the evolutionary model, the model that estimates missing feature values, had the lowest cross-validation errors for all traits, compared with the mean and nearest-neighbour model, respectively (Table 5). The accuracy of inferred trait values was further determined by re-estimation of measured trait values. On average re-estimated values of measured traits were within 1.3 ± 0.13% (mean ± standard error, T_b), 3.7 ± 1.03% (T_{b min}), 1.1 ± 0.10% (TMR), 4.2 ± 1.57% (TBD_max), 8.2 ± 2.92% (TBD_mean) and 5.7 ± 1.71% (BMR) of observed values, respectively.

Phylogenetic signals

Phylogenetic signals for individual measured traits and their relationships ranged from low to high for daily heterotherms, consistent with the results of Cooper & Geiser (2008). However phylogenetic signals in hibernators were absent, with one exception, that is the relationship between TBD_mean and T_{b min}. The reason for the absence of a phylogenetic signal for most relationships between traits and individual traits in hibernators is most likely due to the low number of species analysed in this group (n = 6) and the high uniformity of the data. This is in contrast to daily heterotherms, which are physiologically diverse, have a large body mass range and belong to four different orders. Despite the lack of a phylogenetic signal in the hibernators, we decided to analyse hibernators and daily heterotherms separately, because most physiological variables related to heterothermy, such as TBD_mean, TBD_max, T_{b min} and TMR, differ substantially between these two groups.

In daily heterotherms, the PGLS models always fitted the data better than OLS models and should therefore be preferred. However, there are situations when it may be inappropriate to remove the phylogenetic signal (i.e. use PGLS models), for example if explanatory traits are also closely correlated with phylogeny. Nevertheless, our results show that the phylogenetic signals for all individual measured traits and measured trait relationships examined in daily heterotherms were significantly different from λ = 0 (with the exception of TBD_mean), and thus phylogenetic correction for these traits was warranted. Furthermore, for two trait relationships (TBD_mean vs. TMR and T_{b min} vs. TMR/BMR) and five individual traits (BMR, TMR, T_b, TBD_mean and TBD_max), the strength of the phylogenetic signal estimated by the ML λ was also significantly different from 1. This indicates that the evolution of these relationships and traits was less than expected under the Brownian motion model. Our results therefore emphasize that assuming a set value for λ, for example 0 as in OLS analysis or 1 as in PIC analysis (Felsenstein, 1985), will likely under- or overestimate the influence of shared evolution, which can lead to false conclusions as pointed out by Capellini et al. (2010).

Estimation of missing values using phylogeny

In our analysis, we estimated missing values for several traits related to heterothermy using a phylogenetic approach. Although it needs to be noted that the observed diversity in any physiological parameter will reflect sources of variation other than body mass and phylogeny, as, for example, genotypic adaptation or phenotypic plasticity. Thus, for simplicity, the procedure used here to estimate trait values assumes that body mass and phylogeny are the only sources of variation. However, in a recent study, it has been shown that estimating missing values using phylogeny is the most effective and accurate way to date in estimating missing values, and thus not only increases the number of values in the database, but also increases the confidence and reduces bias in the regression coefficients (Riek & Bruggeman, 2013), which is much more precise than estimating missing values by allometry alone. A problem often encountered in phylogenetic analyses is either the lack of trait data for a species with a known phylogenetic position or the lack of information on the phylogenetic position of this species because of the lack of, for example, molecular data. Usually, if information on either is missing (trait data or phylogenetic position), this species would be excluded from the analysis, resulting in a reduced confidence and increased bias of the regression coefficient. However, we had a fully resolved phylogeny (Bininda-Emonds et al., 2008) and missing trait values for some species. Removing missing rows from the dataset and thus not including the estimates in the analysis, for example for the relationship between body mass and TMR in daily heterotherms, had a direct effect on the allometric exponent. which changed from 0.89 (error 0.069) to 0.77 (error 0.133), reducing the coefficient of determination from 0.84 to 0.61. Importantly, however, the relationships remained significant. Thus, including estimates of missing values, derived on the basis of their phylogeny and body mass, greatly decreases the error of the regression coefficient. Garland & Ives (2000) already demonstrated that trait values can be estimated with higher accuracy and reduced uncertainty if the phylogenetic position of a species is taken into account. Furthermore, the approach of predicting missing trait values based on body mass and phylogeny will be very useful for future tests of hypotheses regarding adaptive variation.

Heterothermy in marsupials

Our comparative analysis shows that many marsupials are heterothermic, expressing either daily torpor or hibernation. Torpor is known to occur in 15% of all marsupials, well above the percentage for placental mammals or birds (Geiser & Ruf, 1995; McKechnie & Lovegrove, 2002; Lovegrove, 2012). The high proportion of marsupial heterotherms is likely linked to some extent to a disproportionally large research effort on the group, which in turn seems to be related to their perceived ‘primitive’ position that may allow a better understanding of the evolution of endothermy and heterothermy, rather than a genuinely higher proportion of heterotherms in marsupials than in placentals, of which about two-thirds are bats and rodents containing many heterothermic species.

The diversity of torpor use in marsupials is not only obvious at the species level, but also at the order and family level. Whereas in marsupials, most likely at least six of seven orders contain heterothermic species (it is highly probable that at least some of the rat-opossums, Paucituberculata, are heterothermic), only about half of the placental orders contain known heterotherms. This is reflected in the large number of marsupial families that contain heterothermic species. These are more than half of the marsupial families (n = 10), if again the rat-opossums (Caenolestidae) are included.

With regard to food of heterothermic marsupials, it is not surprising that torpor has been observed only in families that eat meat, insects, nectar or other unpredictable and fluctuating food items. In contrast, families that rely on predictable leaves and grass, or are omnivorous and are able to switch diets, seem to contain mainly or exclusively homeotherms. The other obvious difference between homeothermic and heterothermic marsupials is size, as in most placental mammals, and is explained to a large extent by relative surface area and mass-specific energy requirements and storage. Heterothermic marsupials weigh up to 1 kg, but most weigh less than 50 g. This differs from placentals at the upper end of the body mass scale, because several large carnivores weighing up to about 100 kg are heterothermic (Tøien et al., 2011).

The occurrence of torpor in certain taxonomic groups and the strong phylogenetic signal detected in some of the variables examined here raise the question regarding the evolution of heterothermy in marsupials. If heterothermy is an ancestral trait in the entire infraclass (Geiser, 2008), those taxa that are now homeothermic must have lost the ability to use torpor because it was no longer required because of size or diet. If torpor is polyphyletic in the group, it must have evolved secondarily from homeothermic ancestors in those taxa that were energetically challenged by their size, diet or habitat, but not in those that did not need to use it. Although the plesiomorphic explanation appears more plausible and is also supported by developmental evidence, recent data show that torpor expression in marsupials is not unvarying within a species, but is affected by the environmental conditions of their habitat (Geiser & Ferguson, 2001), which is likely due to selection, and also is known to occur in placentals (Dunbar & Brigham, 2010; Zervanos et al., 2010; Stawski & Geiser, 2011). Moreover, short-term exposure to different temperatures during development can modify torpor expression within one generation (Riek & Geiser, 2012). Individuals that experience low temperatures during development and growth show more frequent and longer torpor than those raised under warm conditions, demonstrating that torpor expression is subject of developmental phenotypic plasticity. These findings show that variables of torpor are not entirely species-specific constants, but rather are affected by both short-term and long-term factors that may obscure the evolutionary history of torpor.