A quantitative genetics perspective on the body-mass scaling of metabolic rate

Since the 1800s, the law-like nature of the BM scaling of MR has attracted the attention of many kinds of scientists seeking a general understanding of how life works. Traditionally, metabolic scaling has been thought to follow a 2/3 or 3/4 power law that has been explained in terms of fundamental geometric and physical principles, e.g. surface area to volume relationships (e.g. Sarrus and Rameaux, 1839; Rubner, 1883; Roberts et al., 2010; and other references cited in Glazier, 2014a, 2018a,b) and the geometry and physics of internal resource–transport networks (e.g. von Hoesslin, 1888; Kleiber, 1961; Sernetz et al., 1989; Spatz, 1991; West et al., 1997; Banavar et al., 2010; and other references cited in Glazier, 2014a, 2018b). However, recently, many studies have documented extensive variation in metabolic scaling slopes (b) that is systematically related to various intrinsic and extrinsic factors, including taxonomic affiliation, physiological state, body shape and composition, activity level, developmental stage, mode of thermoregulation, ecological life style, and numerous abiotic and biotic environmental factors (reviewed by Glazier, 2005, 2010, 2014a,b, 2018a, 2020; White and Kearney, 2013, 2014). Several mechanisms have been proposed to explain metabolic scaling and its variation, many of which can be classified into four major categories: (1) geometrical constraints on the exchange of resources and metabolic wastes, including heat, across body surfaces with variable relationships to body volume; (2) physically constrained optimization of internal transport of resources to metabolizing cells via variable branching vascular networks; (3) variable metabolic demand of various whole-body, volume-related processes, such as growth and locomotor activity; and (4) variable size-related changes in the relative amounts of tissues with different MRs (reviewed by Glazier, 2014a, 2018b).

Although most of the mechanisms proposed to explain metabolic scaling (such as those just mentioned) have focused on functional and structural properties of whole organisms, some explanations have invoked mechanisms operating at the cellular and biochemical levels (reviewed by Kleiber, 1961; Glazier, 2005, 2014a, 2015, 2018b; Agutter and Tuszynski, 2011). Many of these hypothetical suborganismal mechanisms involve body-size-related changes in the surface area or composition of cellular or subcellular membranes and (or) rates of flow of protons, ions, nutrients and metabolic wastes across them or within subcellular networks. For example, different cellular modes of organismal growth, involving increasing cell size or number or both, determine the total amount of cellular surface area available for uptake of resources and release of metabolic wastes in organisms of different size, with potential consequences for whole-organism metabolic scaling (Davison, 1955; Kozłowski et al., 2003, 2020; Schramm et al., 2021).

Although a subject of ongoing debate, both suborganismal and whole-organism systemic effects likely contribute to metabolic scaling (Glazier, 2014a,b, 2015, 2018b; Kozłowski et al., 2020). Accordingly, several investigators have advocated a hierarchical view of metabolic scaling (e.g. Darveau et al., 2002; Suarez and Darveau, 2005; Suarez, 2012; Glazier, 2014a). Glazier (2014a) has further suggested that a comprehensive theory of metabolic scaling should consider both upward and downward causation, as well as both proximate (functional) and ultimate (evolutionary) causes. Although molecular and cellular models of metabolic scaling provide potentially useful proximate mechanisms based on upward causation, they appear to be insufficient to explain the full diversity of metabolic scaling exponents (b) that have been observed. To explain this diversity, mechanisms involving downward causation via whole-organism systemic effects should also be considered (Glazier, 2014a). In addition, although many kinds of functional mechanisms have been proposed to explain metabolic scaling, ecologically mediated evolutionary mechanisms (ultimate causes) have largely been neglected.

Studies on the evolution of metabolic scaling have only recently begun, and show promise for causing a paradigm shift in our understanding of the relationships among MR, body size, and various other covarying biological and ecological factors. These studies have included interspecific phylogenetic analyses (Uyeda et al., 2017; White et al., 2019), experimental and comparative analyses of conspecific genetic strains or populations exposed to different selection regimes (Czarnołęski et al., 2008; Glazier et al., 2011, 2020), tests for the presence of correlational selection on MR and BM (Artacho et al., 2015; White et al., 2019; Videlier et al., 2021a), comparative analyses of species with markedly different ecological life styles and associated vulnerability to predatory mortality (e.g. pelagic versus benthic invertebrates; Glazier, 2005, 2006), and application of quantitative genetic analyses at the population level to interspecific metabolic scaling patterns (White et al., 2019). In particular, Czarnołęski et al. (2008) have shown how artificial selection on growth rate significantly alters metabolic scaling of the land snail Helix aspersa. Similarly, Glazier et al. (2011, 2020) have provided evidence for effects of size-selective fish predation on the BM scaling of MR and other associated traits, such as gill surface area and rates of feeding and growth, in multiple populations of a freshwater amphipod crustacean. In addition, Videlier et al. (2021a) used a novel multivariate approach to test whether natural selection has acted on both the variance and covariance of MR and BM in the fruit fly Drosophila melanogaster.

Ultimately, allometric scaling patterns within and among species are the result of the correlated evolution of traits (Lande, 1979; Ketola and Kotiaho, 2012). Therefore, a comprehensive understanding of the evolution of metabolic scaling requires analyses of not only what drives evolutionary changes in the phenotypic covariation of MR and BM (e.g. specific kinds of natural selection), but also the physical and genetic factors that permit or constrain these changes. Although several potential physical constraints have been identified, as already noted, hardly anything is known about how the quantitative genetic architecture of MR and BM may affect their coevolution. Furthermore, hardly anything is known about how short-term micro-evolutionary processes at the population level (e.g. natural selection, mutation and genetic drift) result in long-term macro-evolutionary patterns at the phylogenetic level (e.g. interspecific scaling relationships).

Our focus in this study is on how a quantitative genetic perspective and (co)variance partitioning approach (e.g. Hagmayer et al., 2020) can increase our understanding of the evolution of metabolic scaling. Several quantitative genetic studies have already revealed that both MR and BM exhibit significant additive genetic variance (V_A; i.e. variance caused by the additive action of genes within and between loci, that is responsible for the resemblance between parent and offspring) and narrow-sense heritability (h²; i.e. the ratio of V_A to the total phenotypic variance), and thus have the potential to evolve (Konarzewski and Książek, 2013; White and Kearney, 2013; Pettersen et al., 2018; Baškiera and Gvoždík, 2021), which has been corroborated by various artificial selection experiments (e.g. Falconer, 1955; Konarzewski et al., 2005). We also know that MR and BM share a large proportion of their V_A (White et al., 2019), as indicated by several reports of significant genetic correlations (r_A; i.e. the proportion of V_A shared between two traits due to pleiotropy and linkage disequilibrium) between these traits (Nilsson et al., 2009; Bushuev et al., 2012; Zub et al., 2012; Boratyński et al., 2013; Careau et al., 2011; Rønning et al., 2007; Mathot et al., 2013; see also Munday et al., 2017; Nespolo et al., 2014; Tieleman et al., 2009). Although a significant r_A may indicate a limited potential for the independent evolution of MR and BM, it by itself provides little direct genetic information about the quantitative nature of their scaling relationship. Therefore, a new genetically based estimate of the scaling slope between MR and BM is needed. Hagmayer et al. (2020) used a (co)variance partitioning approach to estimate the mass-scaling exponent of MR at the among- versus within-individual levels, but did not have the data needed to further partition the mass-scaling slope at the permanent environment and (additive) genetic levels. To our knowledge, only one study has estimated a mass-scaling exponent of MR based on the additive genetic variance and covariance of MR and BM, herein called the additive genetic scaling exponent (b_A). Videlier et al. (2021b) found that, in D. melanogaster, b_A (±s.e.) was 0.66 (±0.16) in females and 0.58 (±0.32) in males. We clearly need more b_A estimates in other taxa to evaluate the generality of these findings and the concordance in metabolic scaling at the micro-evolutionary scale [i.e. b_A, a manifestation of the genetic (co)variance in BM and MR] versus macro-evolutionary scale (i.e. the widely described negative metabolic allometry at the interspecific level).

Using a quantitative genetics approach to study the evolution of metabolic scaling requires considering b_A as a property of a population (just like h² and r_A), which might reveal not only how the evolution of the metabolic scaling of populations and species may be channelled genetically, but also how the genetic covariance between MR and BM may change through space and time via various evolutionary processes like mutation, genetic drift and natural selection (Ketola and Kotiaho, 2012). Genetic correlations may not only facilitate or constrain evolution, but they themselves may also ‘evolve’ as mutation, drift and selection alter allele frequencies and patterns of genetic pleiotropy and linkage (see e.g. Armbruster and Schwaegerle, 1996). Accordingly, quantitative geneticists have devoted a considerable amount of effort to understanding the factors and mechanisms that can cause changes in the additive genetic variance–covariance matrix (G; Arnold et al., 2008; Careau et al., 2015). The G matrix for BM and MR contains V_A for these traits along the main diagonal, in addition to the additive genetic covariance (Cov_A) along the off diagonal, which represents the additive genetic variance shared by two traits. Quantifying G is the first step to estimate b_A, which is calculated as the Cov_A divided by the V_A in BM. Therefore, changes in G will cause changes in b_A, which then sets the orientation of the ‘genetic line of least resistance’ along which BM and MR are most likely to evolve as populations and species diverge through mutation, genetic drift and natural selection (Schluter, 1996). Hence, we can employ tools and approaches from the field of quantitative genetics to better understand how G and b_A evolve by looking at how they differ in various taxa, evolutionary regimes and ecological conditions, and potentially link micro- to macro-evolutionary patterns of metabolic scaling.

This study has four major aims. First, we use bivariate animal models to reanalyse several published datasets and present new estimates of b_A in multiple animal species. Second, we examine the extent to which the population-specific b_A estimates are similar to the phenotypic metabolic scaling exponents (b_P). Third, we took advantage of the multiple G matrices quantified using scale-independent, log-transformed traits to compare the evolutionary potential in BM versus whole-animal MR. Finally, we discuss how the estimation of b_A and G may increase our understanding of metabolic scaling from an evolutionary perspective.

We searched the literature for published studies that have quantified V_A and/or h² in MR and BM. We then contacted the authors and asked them to send the raw data and pedigree information necessary for estimating V_A in BM and MR and their Cov_A (i.e. G). Note that the presence of some V_A is required for Cov_A, and therefore we did not consider datasets from publications that reported no significant V_A in MR (e.g. Bacigalupe et al., 2004; Mattila and Hanski, 2014; Nespolo et al., 2014; Baškiera and Gvoždík, 2021). The lack of V_A in MR simply implies that there is no additive genetic basis to its mass-scaling relationship. In practice, estimating b_A in a dataset for which there is little V_A in either BM or MR will generate model convergence issues and highly uncertain b_A estimates (i.e. large s.e.), and as such would be potentially problematic and uninformative.

Whenever BM and MR measurements are made for a set of relatives, a bivariate animal model allows the partitioning of the phenotypic (co)variance into genetic variance and other sources of variance depending on the structure of the data and pedigree (Wilson et al., 2010). Each dataset was analysed separately, following a simple workflow analysis (see Supplementary Materials and Methods). For each dataset, we first log₁₀ transformed BM and MR to allow properly scaled comparisons of variation in different-sized animals and of different traits (Glazier, 2021). Moreover, the (co)variance estimates can be directly used to calculate the slope of the log–log relationship between BM and MR (see below). We then included log₁₀-transformed BM and MR as response variables in a bivariate animal model fitted using ASReml-R version 4 (Butler et al., 2018).

Each model included a different set of fixed effects depending on the design of the study and relevant covariates, such as age, sex, generation, habitat, diet quality, season and temperature, and nuisance variables, such as measurement block, cage, location, time of day, etc. (see Table 1 for a list of covariates and original publications for an explanation of the different variables). We used the same set of fixed effects as used in each original publication of the data, fitted to both BM and MR such that their (co)variance was modelled independently from (i.e. conditioned on) the set of covariates. For example, in Zub et al. (2012), weasels were captured in three habitats (river valleys, meadows and forest); adding habitat as a fixed effect yielded (co)variance estimates that can be interpreted as representing the (co)variance structure in a homogeneous environment [see Wilson, 2018 for the importance of considering fixed effects when modelling (co)variances].

After accounting for fixed effects, the remaining phenotypic variance (V_P) in BM and MR was partitioned into various variance components using random effects. In all cases, V_A was estimated by including a random effect of the identity of the animal linked to the pedigree. In most studies, multiple offspring were measured from a given full-sib family, in which case dam identity was included as a random effect to capture common environmental variance (V_C; i.e. environmental effects that are shared by full-sib offspring from a given dam). In one case, dam identity was not known, but repeated measures were taken on multiple individuals, and therefore individual identity was included as a random effect to capture permanent environmental effects (V_PE; i.e. environmental effects that are specific to each individual). Note that V_C and V_PE may also capture variation from non-additive genetic effects (e.g. dominance variance). The residual variance (V_e) captured variation due to measurement error and specific environmental effects that are unique to each measurement (i.e. micro-environmental effects – any intra-individual variance due to variability in environmental factors other than those included as fixed effects in the model).

Whenever possible, variances in BM and MR were modelled within an unstructured (co)variance matrix structure, thus effectively capturing the Cov_A, common environmental covariance (Cov_C), permanent environmental covariance (COV_PE), and residual covariance (Cov_e) between BM and MR. The b_A between BM and MR was calculated as their Cov_A divided by the V_A in BM. Similarly, the residual scaling slope (b_e) was calculated as their Cov_e divided by the V_e in BM. Where applicable, the common environment scaling slope (b_C) was calculated as their Cov_C divided by the V_C in BM. The b_P between BM and MR was calculated as their Cov_P (sum of Cov_A, Cov_C, Cov_PE and Cov_e) divided by the V_P in BM (sum of V_A, V_C, V_PE and V_e). The s.e. for b_A, b_e, b_C and b_P was calculated with the delta method, and 95% confidence intervals (CI) were calculated as b_A±1.96×s.e. We used the R package ‘metafor’ (Viechtbauer, 2010) to calculate weighted meta-analytical mean estimates (i.e. for b_A, b_e, b_C and b_P) and conduct a formal test for heterogeneity variance across the b_A estimates.

The variance components from the fitted models also allow comparing evolutionary potential in BM and MR, defined as the expected evolutionary response to selection per strength of selection. We therefore calculated h² by dividing V_A by V_P (the sum of the variance components V_A, V_C, V_PE and V_e). Using h² to infer evolutionary potential, however, can be problematic, and Houle (1992) suggested using mean-scaled measures of evolvabilities. Given that BM and MR were log₁₀ transformed, we can compare their V_A directly because the variance of the log-transformed data approximately equals the mean-scaled variance on the original scale (i.e. the V_A of the log of a trait can be used as an estimate of evolutionary potential for the trait on the original scale; Hansen et al., 2011).

We obtained a total of 11 datasets suitable for estimating V_A in BM and MR, their Cov_A and, subsequently, their b_A (Table 1, Fig. 1). Overall, eight species are included, with two ectotherm species (D. melanogaster and Nauphoeta cinerea) and six endotherm species (two birds and four mammals). The average sample size is N=612 individuals measured, but considerable variation exists, with studies ranging from approximately 100 to more than 1000 individuals measured (Table 1).

Bivariate animal models revealed a significant and positive Cov_A between BM and MR in most datasets (Tables S1–S10). The resulting b_A estimates exhibited considerable variation, ranging from 0.187 to 0.972 (Table 2). Most b_A estimates were between 0.6 and 0.8, but each estimate had large uncertainty, and therefore none were significantly different from the theoretical value of 2/3, and only one estimate was significantly lower than 3/4 (Table 1, Fig. 1).

The meta-analytical average (±s.e.) b_A was 0.667 (±0.050), with 95% CI that exclude both 0 and 1, but include both 2/3 and 3/4 (Table 2, Fig. 2). By contrast, the meta-analytical average (±s.e.) b_P was 0.596 (±0.051), with 95% CI that include 2/3, but exclude 0 and 3/4 (Table 2; lower 95% CI=0.496; upper 95% CI=0.696). Mean b_A and b_P were not significantly different, as indicated by the 95% CI of each overlapping the mean of the other. In addition, the population/species estimates for b_A and b_P were positively correlated (r=0.797, N=11, P=0.003; Fig. 3A).

The meta-analytical mean residual scaling exponent was b_e=0.503±0.065 (Fig. 2), with 95% CI that exclude 0, 2/3, and 3/4 (lower 95% CI=0.375; upper 95% CI=0.630). Across the 11 estimates, b_e and b_P were strongly and positively correlated (r=0.835; Fig. 3C). There were only four valid estimates for b_C because V_C in BM or MR was frequently too small to allow fitting a Cov_C (see Tables S1–S10). The meta-analytical mean common environment scaling exponent was b_C=0.880±0.077 (Fig. 2, Fig. 3B).

The meta-analytical average h² (±s.e.) for BM (h²=0.396±0.059) was slightly higher than h² for whole-animal MR (h²=0.316±0.053; Table 3, Fig. 4A), but the difference was not significant when considering the uncertainty around the means. Similarly, the meta-analytical average V_A (±s.e.) was slightly higher for BM (V_A=0.00100±0.00018) than for whole-animal MR (V_A=0.00080±0.00013; Table 3, Fig. 4B), but the difference was not significant when considering the uncertainty around the means.

We compiled multiple datasets that were originally collected to estimate V_A for whole-organism or mass-adjusted MR, and re-analysed them to quantify the G matrix for BM and MR, which is required to calculate the b_A. Interestingly, the interspecific variation in b_A (range, 0.187–0.972) reported in Table 2 is within the range of b_P estimates (b_P≈0 to >1) that have been observed in animals for both within- and across-species relationships (see Glazier, 2005, 2010, 2014a,b; White and Kearney, 2014). This diversity is consistent with growing evidence that phenotypic metabolic scaling is highly diverse, but unfortunately the large uncertainty (95% CI) around the genetic metabolic scaling exponents prevents any significant interspecific differences from being discerned. In addition, all but one of these exponents are not significantly different from the classical theoretical values of 2/3 or 3/4; two values are not significantly different from 0, and six are not significantly from 1. The only exception is for one of the datasets with a large sample size (Sadowska et al., 2005), in which b_A was significantly lower than the theoretical 3/4 exponent (lower 95% CI=0.746). Therefore, extremely large sample sizes are required to estimate b_A with enough precision to test metabolic theory that predicts specific exponents such as 2/3 or 3/4.

The large error terms for b_A also prevent any clear-cut matching of the genetic metabolic scaling exponents with phenotypic metabolic scaling exponents reported in the literature. For example, no clear differences between the genetic scaling exponents of the ectothermic insects and endothermic birds and mammals can be seen (Table 2), despite the existence of significant differences in the phenotypic scaling exponents of ectothermic and endothermic vertebrate animals (reviewed in Glazier, 2005, 2010; White and Kearney, 2014). Moreover, existing data are insufficient to permit a consistent, clear-cut matching of genetic and phenotypic metabolic scaling exponents for individual species. For example, although the genetic metabolic scaling exponent for the fruit fly D. melanogaster (0.523±0.232 95% CI) is significantly different from 0, Van Voorhies et al. (2004) found no significant phenotypic correlation between MR and BM in this species. A limited BM range (<2-fold) may have prevented a phenotypic correlation between MR and BM from being detected. By contrast, Ellenby (1945, 1953) found significant positive correlations based on samples with larger body-size ranges of ∼2- and 3-fold, respectively. Furthermore, the phenotypic scaling exponent of 0.554 based on data in Ellenby (1945) is within the 95% CI of the genetic scaling exponent, although the scaling exponent of 0.772 based on data in Ellenby (1953) is not. In addition, the genetic metabolic scaling exponent for the laboratory mouse Mus musculus (0.826±0.323) is not significantly different from the phenotypic scaling exponent (0.748±0.289) reported by Glazier (2022), based on data from 15 studies collated by Genoud et al. (2018). These examples thus suggest that, with sufficient data, a match between genetic and phenotypic metabolic scaling exponents may be found.

In comparison to b_A and b_e, the meta-analytical b_C tended to be higher (b_C=0.880±0.077), raising questions about sources of common environmental covariance in BM and MR. In our models, V_C and Cov_C captured environmental effects (other than those specified as fixed effects) that applied to groups of full sibs. Because relatives (especially full sibs) are clustered in time and space (e.g. litter or nest effects), common and maternal environmental effects are usually confounded with the pedigree structure, and therefore may bias V_A and Cov_A estimates if not controlled for (Wilson et al., 2010). Because b_C tends to be higher than b_A, a failure to properly separate common environment from additive genetic effects will introduce an upward bias in b_A. Interestingly, V_C and Cov_C can also capture maternal genetic effects and non-additive genetic variance such as dominance variance (Wilson et al., 2010). Maternal genetic effects have been widely reported for growth-related traits and body size (McAdam et al., 2014), such that V_C (and Cov_C) may contain some source of maternal genetic (co)variance. Moreover, owing to the maternal route of inheritance of mitochondria, genetic variance in the function, structure and intracellular density of mitochondria should end up in the maternal effect (unless modelled explicitly; see Evans et al., 2014). This raises the possibility that V_C and Cov_C contain some cytoplasmically inherited genetic (co)variance in BM and MR that scales more steeply than b_A. In any case, future quantitative genetic studies of metabolic scaling should pay attention to common environmental and maternal effects, and if possible calculate and report b_C to add to the few (n=4) currently available estimates.

Overall, our collective data on genetic metabolic scaling exponents (Table 2) are consistent with the frequent observation that intraspecific and interspecific metabolic scaling exponents are negatively allometric in animals (b<1). The weighted average b_A across all datasets is 0.667±0.050, which is significantly less than 1 (see Table 2). Furthermore, all b_A estimates are <1 rather than >1, which is highly improbable by chance alone (P=0.5¹⁰=0.001). Although there was large variation in the b_A estimates, when considering the differences in light of the standard errors that come with each estimate, the possibility emerges that sampling noise alone may be responsible for all the variation that we see. The meta-analytical procedure yielded an I² (total heterogeneity/total variability) of 0, which means that there are no differences between the b_A estimates that cannot already be explained by sampling noise alone. Hence, one could argue that extreme values like the low value from Saxicola is fully expected from sampling noise alone (owing to its large uncertainty).

The b_A represents the direction of greatest genetic variance in BM and MR, and therefore divergent species are expected to be constrained along this general orientation (Schluter, 1996). Just as we calculated b_A from the additive genetic (co)variance matrix, one could use a bivariate phylogenetic mixed model (Hadfield and Nakagawa, 2010) to estimate the phylogenetic (co)variance in mass and MR and quantify the mass scaling exponent at the phylogenetic level (b_D), which represents the scaling exponent along which BM and MR co-varied as species diverged. A correspondence between b_A and b_D may be explained in two non-exclusive ways. First, b_A may channel the evolution of the phenotypic covariance of MR with BM, and thus b_D, along the genetic line of least resistance (Schluter, 1996). In other words, the distribution of the genetic variation in BM and MR in an ancestral population or species should influence the orientation (slope) along which BM and MR will vary among descendant populations or species. Second, natural selection may favour specific phenotypic and corresponding genetic associations between MR and BM. Indeed, selection on particular combinations of BM and MR could both (1) modify patterns of genetic (co)variance, and (2) cause adaptive evolution in BM and MR along a specific orientation.

Note that the above genetically based evolutionary explanations do not mention possible influences of physical constraints on metabolic scaling, as posited by traditional theory (see Introduction). Nevertheless, we do not mean to imply that physical constraints are not at all involved in metabolic scaling. For example, they could set boundary limits on the range of metabolic scaling exponents that can evolve (see e.g. Glazier, 2005, 2010, 2014b). However, at present, no firm conclusions about cause and effect relationships among b_A, b_P and b_D can yet be made. More interdisciplinary, multi-level work is needed to determine how genetic patterns of metabolic scaling at the intraspecific level are linked to interspecific metabolic scaling patterns (see also White et al., 2019). Ultimately, there can be no genetic variation outside of the fundamental physical and geometrical constraints of organisms, such that genetic and physical constraints should theoretically reflect two sides of the same coin (just like genetic correlations and physiological mechanisms linking complex traits; Careau and Garland, 2012).

Currently available quantitative genetic data are not yet sufficient to test metabolic scaling theory in a rigorous way. Although the meta-analytical weighted average b_A of 0.667±0.050 appears to support the view that additive genetic variance in BM and MR is distributed along a metabolic scaling exponent of 2/3, the 95% CI (0.567–0.767) also include 3/4. Moreover, the uncertainty around each b_A estimate was so large that none of the individual estimates clearly support any specific theoretical metabolic exponent, such as 2/3 or 3/4. They are also insufficient to test theory that predicts diverse metabolic scaling exponents, such as the metabolic-level boundaries hypothesis (MLBH: Glazier, 2010, 2014b). However, it is intriguing that, as predicted by the MLBH, the birds and mammals with relatively high basal metabolic rates (BMRs) [Saxicola torquata, Taeniopygia guttata, Myodes glareolus and Mustela nivalis, with mean BMR values of 110.7%, 117.9%, 141.7% and 155.1%, respectively, greater than that predicted from the eutherian mammal log-linear scaling relationship BMR=0.455+0.726*BM reported by Genoud et al., 2018] tend to have lower genetic metabolic scaling exponents (0.187–0.813) than those species with substantially lower MRs (M. musculus and Peromyscus maniculatus, with b_A values of 0.826 and 0.972, and mean BMR values of 23.5% and 47.2%, respectively, greater than that predicted by the eutherian scaling relationship). The MLBH predicts that the scaling exponent for resting (or basal) metabolic rate (RMR or BMR) should decrease as overall metabolic level (L) increases (Glazier, 2010, 2014b). Nevertheless, the negative relationship between b_A and L observed here is based on only six species with b_A values that have large error terms, and therefore is unsurprisingly not statistically significant. The MLBH also predicts that the scaling exponent for active maximal metabolic rate (MMR) should be higher than that for RMR, which, however, is not seen for either b_P or b_A in the house mouse M. musculus (0.738±0.102 and 0.880±0.022, respectively, for MMR, and 0.826±0.165 and 0.933±0.041, respectively, for RMR; Wone et al., 2009). However, activity has been shown to increase b_P within other animal species, e.g. the laboratory rat (Refinetti, 1989), humans (Rogers et al., 1995; Batterham and Jackson, 2003) and various ectothermic animals (Glazier, 2009). More data are obviously needed.

Available data on V_A given in Table 3 do permit a useful test of two different theoretical approaches regarding the evolution of metabolic scaling that make opposite predictions. The first theory posits that metabolic scaling is chiefly the result of evolutionary optimization of body size with secondary effects on MR (Kozłowski and Weiner, 1997; Kozłowski et al., 2020). According to this theory, as applied by Ketola and Kotiaho (2012), if natural selection acts more intensely on body size than MR, then V_A (and h²) should be lower for body size than that for MR. The second approach, based on dynamic energy budget theory, proposes that evolutionary optimization is focused primarily on MR with secondary effects on body size (Lika et al., 2019). If correct, then V_A (and h²) should be lower for MR than body size. However, comparisons of h² and V_A estimates from our models (Table 3, Fig. 4; see also Tables S1–S10) do not support the predictions of either of these theories, as there were no differences in the average h² and V_A estimates for MR versus BM (note that V_A estimated using log-transformed traits allows directly comparable estimates of mean-scaled evolvabilities; Hansen et al., 2011). Therefore, although natural selection may act on MR and BM to varying degrees in different environments to cause evolutionary changes in metabolic scaling (see also Witting, 2017), existing data suggest that there is no general difference in the intensity of selection on – and evolutionary potential in – either trait. In agreement, Videlier et al. (2021a) used a multivariate approach based on the variance–covariance matrix of Darwinian fitness, standard metabolic rate (SMR) and BM to show that linear selection was acting on both SMR and BM in male, but not female, D. melanogaster. However, Boratyński and Koteja (2010) found that BMR, but not BM, was related to relative fitness (reproductive success) in the bank vole M. glareolus. More studies of this kind are needed, as well as those that specifically test for effects of selection on the covariance of MR and BM, and thereby metabolic scaling.

In some investigations, the metabolic scaling exponent (b) is considered to be a trait in and of itself. For example, Fossen et al. (2019) did so in a study of the V_A and genotype-by-environment (G×E) interaction of b using clones from a parthenogenetic species (Daphnia magna). Although it may be useful to consider b as a trait when working with genetically identical clones or strains (e.g. Glazier and Calow, 1992; Yashchenko et al., 2016; Matoo et al., 2019), or ontogenetic b values estimated at the individual level (Norin and Gamperl, 2018; Beaman et al., 2021 preprint), such an approach can be impractical for studies that estimate b among genetically different conspecific individuals. Moreover, the two quantitative genetic studies conducted so far on b (considered as a trait) found no detectable V_A in b (Fossen et al., 2019; Beaman et al., 2021 preprint). For most studies of metabolic scaling, one could reasonably argue that b should not be considered a ‘trait’ per se, but a type of ‘trait covariance’, resulting from pleiotropy, gene linkage and/or selective covariance (Armbruster and Schwaegerle, 1996). Accordingly, genes and natural selection affect b only indirectly via their direct effects on the individual traits of MR and BM, the covariances of which determine metabolic scaling. We can still study how b_A ‘evolves’, in the sense that G (and b_A) may change through space and time as a result of changes in the frequencies of the genes responsible for BM and MR (Arnold et al., 2008). Otherwise, considering b as a trait makes it problematic to study metabolic scaling from a genetically based adaptive perspective, given that b represents the covariance of two complex traits (MR and BM) subject to multiple environmental and genetic effects and variable selection regimes (Arnold et al., 2021).

In this study, we hope to have shown the value of quantitative genetic studies in increasing our understanding of the evolution of metabolic scaling. Other useful methods have been proposed to examine genetic effects on metabolic scaling, including quantitative trait locus analyses (Wu et al., 2002; Long et al., 2006; Li et al., 2007; Vasseur et al., 2012; Zhang et al., 2021), inbreeding and crossbreeding experiments (Ketola and Kotiaho, 2012; Vasseur et al., 2012; Ketola et al., 2013; Boratyński et al., 2016), and artificial selection experiments (Czarnołęski et al., 2008; Sadowska et al., 2015; Malerba and Marshall, 2019). In addition, genetic effects on metabolic scaling can be studied through the comparison of metabolic scaling relationships observed in genetically different clones, strains or populations (Glazier and Calow, 1992; Glazier et al., 2011; Yashchenko et al., 2016; Matoo et al., 2019). We encourage researchers interested in the evolution of metabolic scaling to adopt a quantitative genetics approach that examines how b_A may be influenced by evolutionary processes acting on BM and MR, in accordance with how evolutionary biologists have traditionally evaluated changes in G and genetic correlations under artificial or natural selection (e.g. Careau et al., 2015; Delahaie et al., 2017).

We thank all authors from the original studies for conducting all of the work to collect the data and for sharing the datasets, and Wolfgang Forstmeier for guidance on the meta-analytical approach.