Terrestrial locomotion energy costs vary considerably between species: no evidence that this is explained by rate of leg force production or ecology

Inter-specifically, relative energy costs of terrestrial transport vary several-fold. Many pair-wise differences of locomotor costs between similarly-sized species are considerable, and are yet to be explained by morphology or gait kinematics. Foot contact time, a proxy for rate of force production, is a strong predictor of locomotor energy costs across species of different size and might predict variability between similarly sized species. We tested for a relationship between foot contact time and metabolic rate during locomotion from published data. We investigated the phylogenetic correlation between energy expenditure rate and foot contact time, conditioned on fixed effects of mass and speed. Foot contact time does not explain variance in rate of energy expenditure during locomotion, once speed and body size are accounted for. Thus, perhaps surprisingly, inter-specific differences in the mass-independent net cost of terrestrial transport (NCOT) are not explained by rates of force production. We also tested for relationships between locomotor energy costs and eco-physiological variables. NCOT did not relate to any of the tested eco-physiological variables; we thus conclude either that interspecific differences in transport cost have no influence on macroecological and macrophysiological patterns, or that NCOT is a poor indicator of animal energy expenditure beyond the treadmill.

SCiEntifiC REPORTs | (2019) 9:656 | DOI: 10.1038/s41598-018-36565-z perspective of optimal foraging 28,29 and energy allocation theories 30,31 , such vast variability seems counter to the concept that animals have evolved to economise their daily energy expenditure, within the constraints of their ecology, in order to maximise the energy they have remaining, which can be directed towards reproduction 32 . In turn, this variability calls into question the importance of energetics, or at least locomotory energy expenditure, in shaping an animal's ecology 33 . Yet there has been little research into the ecological implications of these differences in locomotion cost exhibited by similarly sized species. It has been known for a long time that penguins have a high cost of transport 34,35 . Initially, this was believed explicable by their waddling gait. However, Griffin and Kram 36 posited that penguins expend a lot of energy when walking because of their relatively short hind limbs resulting in low foot contact times, and thus high rates of force production 9 . They therefore argued that the proposed explanation for the greater relative NCOT in smaller animals might also explain the greater relative NCOT in the short-legged penguins 4 . By extension, the inference has been made that in general, species exhibiting a high NCOT for their size might have relatively short legs and, conversely, animals that display particularly good energy economies tend to have long legs (see also 11 ). Kram and Taylor's (1990) function relating metabolic rate of locomotion to foot contact time has successfully predicted gait-metabolic cost relationships in various contexts [37][38][39] . However, there has yet to be a substantial across-species analysis to test the hypothesis of a relationship between foot contact time and the cost of transport. Thus although it is often assumed to be the case, we do not know whether rate of force production describes variation in the cost of transport inter-specifically.
Limb morphology was included as an explicit determinant of locomotion energy costs in a biomechanical model (LiMb) developed by Pontzer 40 . The LiMb model recognises that energy expenditure during terrestrial locomotion, while mostly explained by the costs of generating the vertical component of the ground reaction force, also includes the costs of the horizontal component associated with cyclical braking and propulsion 41 and the costs to swing the limbs 42,43 . The model was initially applied to an intra-specific analysis of human participants and performed better at predicting running metabolic rate than did foot contact time 40 . A similar conclusion arose from intra-specific analyses of two quadruped species: dogs and goats 44 .
The LiMb model requires the input of three parameters, two of which are the energy cost of swinging the limb and the excursion angle of the limb during the stance phase. Both of these are generally only obtainable under laboratory conditions. However, the model predicts that over a wide range of body sizes the other and most easily measured input parameter, effective limb length (the functional length of the limb as a mechanical strut), will accurately estimate locomotion costs because it is the only one of the three parameters that scales with body mass. Evidence supporting this prediction came from an across-species study; log-transformed effective limb length explained 98% of the variance in log-transformed mass-specific locomotor costs (as opposed to 94% explained by log-transformed body mass) 45 , in turn supporting aforementioned previous work suggesting that the magnitude and frequency of muscle forces generated to counter gravitational acceleration is a key determinant of NCOT across species [36][37][38][39] (cf 6 ).
However, despite this impressively high coefficient of variation, similarly to the aforementioned relationships between cost of transport and body mass, the use of log-transformed data across a wide mass range masks the absolute size of the residuals; a third of the residuals have a magnitude greater than 20%. Moreover, many comparisons between relevant pairs of species within the dataset do not support the premise that longer effective limbs associate with lower transport costs 45 . This might be explained by effective limb length not relating to rate of force production, for example because it does not account for the excursion angle of the limb (which as mentioned previously is an additional variable in the LiMb model); excursion angle dominates the variation in locomotor cost at the intraspecific level 40,45 . This proposed explanation is indirectly evidenced by an across-species regression, based on the present dataset, of relative effective limb length against rate of force production represented by relative foot contact time; accounting for speed, this regression returned a comparatively low R 2 of 0.61 ( Fig. 1).
Foot contact time is mathematically more closely associated with rate of force production against the ground than is effective limb length. Assuming that rate of force production is a key driver of NCOT across species then foot contact time is perhaps the proxy likely to associate best with variations in NCOT between similarly-sized species despite the great variation in morphologies and walking gaits they exhibit. At the least, it seems reasonable to re-asses this hypothesis, drawing upon the more extensive data sets for foot contact times and metabolic rates of terrestrial animals on the move that are now available. We collated these data to investigate whether, across species, rate of force production represented by foot contact time is predictive of a species' locomotion costs. Using a phylogenetically informed approach, we focussed on birds and mammals (since data for these taxa are by far the most prevalent) to test for an interspecific relationship between metabolic rate during locomotion and foot contact time, while also recognising and accounting for the relationships between each of these variables and both body mass and locomotion speed, using phylogenetic mixed models to formally incorporate the non-independence associated with multiple measurements of a single species. We then expanded the analysis to consider additional predictors of interspecific variation in NCOT, testing for correlations between NCOT and a host of putative eco-physiological correlates, to examine the ecological relevance of NCOT beyond the laboratory. For this analysis we focussed on mammals to take advantage of the rich data sets published by Jones et al. 46 and multiple other sources 24,47-52 .

Materials and Methods
A total of 288 observations of mass-specific rate of oxygen consumption (  V O 2 , ml kg −1 min −1 ) and 126 observations of foot contact time (t c , s) during pedestrian locomotion were obtained for 21 species (combining data for large leghorn and bantam leghorn chickens, and for musculus and musculus longshanks mice; Table 1) mainly from published studies 9,15,34,36,[53][54][55][56][57][58][59][60][61][62][63][64][65][66][67][68][69][70][71] . The species in this study represent a mass range of 32 g to 467 kg. Previously unpublished kinematic data for king penguins are included as part of the present study. In most cases collation of data from the literature required data extraction from digitised versions of published figures. We used SCiEntifiC REPORTs | (2019) 9:656 | DOI:10.1038/s41598-018-36565-z data for  V O 2 (i.e. rather than NCOT) to test for the relationship between transport costs and 1/t c because the relationship between speed and 1/t c is not always linear, and it is therefore not possible to calculate species-specific speed-independent values for foot contact times. Because the y-intercept of the relationship between  V O 2 and speed is greater than zero and typically also greater than resting metabolic rates e.g. 72 , for each data set we calculated a linear regression relating  V O 2 and speed, and subtracted the y-intercept of this relationship from each value of  V O 2 in the data set, following previous studies (e.g 15 .). Our data for  V O 2 therefore represent only the incremental increases in  V O 2 associated with increases in speed. Data for NCOT were obtained from a published compilation 4 . Data for putative eco-physiological correlates of NCOT were obtained from the PanTheria database 46 , namely: home range size, geographical range, group size, terrestriality, diet breadth, trophic level and habitat breadth. Additional variables were obtained from other sources: basal metabolic rate 47 , field metabolic rate 48 , maximum aerobic metabolic rate 49 , daily movement distance 24 , maximum running speed 50 , and body fat 51,52 . Absolute aerobic scope was calculated as the difference between maximum aerobic metabolic rate and basal metabolic rate. Factorial aerobic scope was calculated as the ratio of maximum aerobic metabolic rate to basal metabolic rate. Activity metabolic rate was calculated as the difference between field metabolic rate and basal metabolic rate 73,74 .
The speed-and size-independent relationship between log 10 -transformed mass-specific rate of oxygen consumption (  V O 2 , ml kg −1 min −1 ) and log 10 -transformed 1/t c (s −1 ) was analysed using multivariate phylogenetic mixed models 75-77 . Phylogenetic mixed models were implemented in the ASReml-R v3.0 78 package of R v3.0.2 79 , with inverse relatedness matrices calculated from phylogenetic covariance matrices using the MCMCglmm package v2.21 80 . Phylogenetic mixed models were selected over the more commonly used methods of independent contrasts 81,82 and phylogenetic generalised least squares 82,83 because the former can formally incorporate phylogenetic non-independence as well as non-independence associated with multiple measurements of single species (i.e. measurements of  V O 2 and t c for a single species running at multiple speeds). Phylogenetic mixed models are an analogue of the mixed model from quantitative genetics, which partitions phenotypes of related individuals into heritable (additive genetic) and non-heritable components to estimate inter-specific variances and covariances between traits 77 .
We used multivariate phylogenetic mixed models to examine the phylogenetic correlation among the traits  V O 2 and 1/t c , accounting for either body mass or speed, or accounting for both body mass and speed (a phylogenetic correlation between two traits represents the proportion of variance that these two traits share due to phylogenetic relatedness, after accounting for any fixed effects in the model). The multivariate models included log 10 (  V O 2 ) and log 10 (1/t c ) as response variables, either or both of log 10 -transformed speed (U, m s −1 ) and log 10 -transformed body mass (M, kg) as fixed effects, and phylogenetic relatedness as a random effect (it was not possible to include species identity as a random effect in the multivariate models because there was insufficient data to fit covariance matrices for both phylogenetic relatedness and species identity). Completely parameterised (unstructured) (co) variance matrices were specified for the random effect associated with the phylogeny, as well as the residuals. The significance of phylogenetic correlations between log 10 (  V O 2 ) and log 10 (1/t c ), conditioned on the fixed effects, were inferred by determining if the phylogenetic covariance between log 10 (  V O 2 ) and log 10 (1/t c ) differed significantly from zero by using a likelihood ratio test to compare models with and without the appropriate covariances fixed at zero. Approximate standard errors for phylogenetic correlations were calculated using the R 'pin' function 84 .

Residual effective limb length
Before running the multivariate phylogenetic mixed model, we first visually verified that the relationships between log 10 (U) and both of log 10 (  V O 2 ) and log 10 (1/t c ) were approximately linear for all species (Fig. 2), and used univariate phylogenetic mixed models to test for significant interactions between the fixed effects of log 10 (U) and log 10 (M) for both log 10 (  V O 2 ) and log 10 (1/t c ). The significance of fixed effects was tested using Wald-type F-tests with conditional sums of squares and denominator degrees of freedom calculated according to Kenward and Roger 85 . Phylogenetic heritability (h 2 ), a measure of phylogenetic correlation equivalent to Pagel's 86 λ 77 , was estimated as the proportion of variance attributable to the random effect of phylogeny; the proportion of variance attributable to species identity independent of phylogeny was also calculated. The significance of phylogenetic heritability and the random effect of species identity were assessed using likelihood ratio tests to compare models with and without the random effects. Approximate standard errors for the estimate of phylogenetic heritability were calculated using the R 'pin' function 84 .
To test for correlations between NCOT and each of the putative eco-physiological correlates of NCOT, because there was only a single value per species we used the phylogenetic generalised least squares (PGLS) method 82,83 and the 'ape' v3.2 87 and 'caper' v0.5.2 88 packages of R v3.0.2 79 , with Pagel's λ 86 estimated using maximum likelihood. Body mass was included as a covariate in all analyses, and all data were log 10 -transformed for analysis. The PGLS method is equivalent to univariate phylogenetic mixed models when only one response value per species is analysed, as is the case for NCOT data.
The phylogeny used for all analyses was obtained from the open tree of life 89 using R v3.2.2 and the R package 'rotl' 90 , with branch lengths estimated using the arbitrary method of Grafen 83 in 'ape' v3.2 87 .
Parameter and variance component estimates are shown ± SE, unless otherwise noted, and α was set at 0.05 for all tests. Multiple comparison corrections were not performed 91 .  60 Hoyt and Taylor 59 15 Watson et al. 61

Results
Assessed using univariate phylogenetic mixed models, there was no interaction between log 10 (M) and log 10 (U) as predictors of log 10 (  V O 2 ) ( Table 2) or log 10 (1/t c ) ( Table 3) during treadmill locomotion, indicating that there is no need to include an interaction terms in the multivariate phylogenetic mixed model. With only the additive effects of log 10 (M) and log 10 (U) included, there was no significant phylogenetic heritability for log 10 (  V O 2 ). However, there was significant among-species variation unrelated to phylogeny ( Table 4), indicating that there are consistent mass-and speed-independent differences among-species in log 10 (  V O 2 ) that are not related to patterns of phylogenetic relatedness (i.e. related species are not more similar than unrelated species). Conversely, for log 10 (1/t c ) there was significant phylogenetic heritability but no significant among-species variation unrelated to phylogeny (Table 5); this indicates that there are significant speed-and mass-independent differences in log 10 (1/t c ) among species, but that related species are more similar than unrelated ones.
None of the eco-physiological traits were significantly related to NCOT when body mass was accounted for ( Table 6).

Discussion
Most explanations for animal scaling of mechanical efficiency during locomotion are premised on the concept that energy costs of locomotion derive primarily from the muscular force to support the body and, particularly for larger animals, to accelerate the body as it oscillates through the step cycle 6,40 . Smaller animals have higher stride frequencies due to their shorter limbs, which necessitate higher rates of force generation due to shorter foot contact times 9 . This requires the recruitment of faster, less efficient muscle fibres, which demand higher rates of cross-bridge cycling and Ca 2+ pumping, and thus use more ATP per gram of active muscle 16,92,93 . Consequently, at a higher stride frequency more energy is expended to deliver the necessary ground reaction force to support the body. Superficially, this theory appears to explain the present results, because foot contact time -a proxy for rate of force production -describes variation in the energy cost of locomotion between species of different sizes  Table 2. Parameter estimates for univariate phylogenetic mixed models describing the effect of speed (U, m s −1 ), body mass (M, kg), and their interaction on rate of oxygen consumption (  V O 2 , ml kg −1 min −1 ). Fixed = log 10 (  V O 2 ) ~ log 10 (M) + log 10 (U) + log 10 (M) * log 10 (U). The model includes random effects of phylogeny and species, which are used to determine the variance associated with the phylogeny and the variance associated with species-specific differences not associated with phylogeny, respectively, conditioned on the fixed effects. h 2 is the proportion of variance, conditioned on the fixed effects, accounted for by the random effect of phylogeny (shown ± SE), and is equivalent to Pagel's λ 77,86 ; χ 2 is the test statistic used to determine if the variance estimates for the random effects of phylogeny and species are significantly greater than zero. Estimates associated with 'Phylogeny' and 'Species' are estimates of the variances associated with each of these random effects, conditioned on the fixed effects, and the estimate associated with 'Residual' is the residual variance. Phylogenetic h 2 = 0.13 ± 0.24 (χ 2 1 = 0.69, P = 0.41). Proportion of variance accounted for by 'Species': 0.57 ± 0.21 (χ 2 1 = 14.7, P < 0.001).    locomoting at the same speed (Fig. 3B), as well as of the same size moving at different speeds (Fig. 3A). However, the among-species relationship disappears when both speed and body size are accounted for simultaneously (Fig. 3C). This finding indicates that the apparent relationships among speeds (Fig. 3A) and sizes (Fig. 3B) arise because of the independent effects of speed and size on both foot contact time ( Fig. 2A) and  V O 2 (Fig. 2B), rather than because of an among-species link between (speed-and size-independent) foot contact time and  V O 2 . Stated another way, for animals of similar size locomoting at similar speeds, there is no among-species relationship between foot contact time and  V O 2 . In turn, this suggests, surprisingly, that rate of ground force production is not a substantive influence on NCOT between species. There must therefore be overriding, alternate factors at play determining the often considerable variation in energy cost of transport between species of similar size.
Not only does NCOT not relate inter-specifically to foot contact time, it also does not relate with any of the eco-physiological variables investigated, from home range size to running speed to body fat mass. While such data may be noisy, this consistent lack of correlative evidence raises the question as to whether NCOT -a measure obtained from prescribed and highly controlled protocols -is meaningful beyond the laboratory i.e. meaningful within an ecological context.
In the remainder of the Discussion we first consider what factors other than rate of force generation might explain variation in  V O 2 between species during terrestrial locomotion, and then interpret and consider the implications for the lack of correlation between NCOT and all available eco-physiological variables.
Possible mechanisms underlying NCOT. The phylogenetically informed estimate of the scaling exponent for the relationship between log 10 (NCOT) and log 10 (body mass) is −0.28 4 . Given that NCOT in that study, as in the current one, is presented in mass-specific terms, the whole-animal scaling exponent for NCOT is 1 −0.28 = 0.72, which is similar to the phylogenetically informed estimate of the scaling exponent of BMR of 0.69 ± 0.01 47 . This similarity in turn raises the question as to whether the scaling laws governing NCOT are the same as those governing BMR. However, similarity of scaling exponents does not mean that two traits are functionally related 94 , and our analyses found no association between NCOT and BMR once the effect of body mass on BMR was accounted for. Thus the relationship between NCOT and mass is not explained by the relationship between BMR and mass.
Morphological differences between species are likely to include the amount of active muscle that is applying force during periods of foot-ground contact 1 (see also 95 ). For example, birds have longer legs than do mammals and thus longer muscle fibres and a greater muscle volume, which can explain the considerably higher energy expenditure of birds for a given rate of force production 15,96 . Effective mechanical advantage based on the ratio of the anatomical moment arm of the muscles to the load arm of the ground force vector 97 can also vary considerably between species, including when muscle fibre length has been accounted for 98 . Indeed, Pontzer et al. 98 demonstrate that step length, a strong proxy for foot contact time 9 , is only a limited predictor of NCOT across (eight) species of birds and mammals, arguing that this is because of diversity within those species of muscle fibre length and effective mechanical advantage. Similar effective limb lengths in animal species can be represented by diverse arrangements of the limb skeletal structure resulting in differing degrees of upright or crouched stance; the more upright the limb (the more erect the animal's posture) the lower the cross-sectional area of muscle required to generate sufficient force to counter gravity 37,99 . There may also be differences between species in the relative shortening velocity of active muscle. For example, due to the force-velocity relationship of skeletal muscle, in cases where the fibres are operating at lower shortening velocities a smaller cross-sectional area of muscle would need to be activated to provide the same force as muscles operating at higher shortening velocities 15,100 . Recent work by Cavagna and Legramandi 101 provides a new approach to considering the biomechanics underpinning the scaling of locomotion efficiencies with size. They demonstrated that there is variability in the hysteresis energy loss experienced by animals during the lower part of the vertical oscillation of the centre of mass. While muscle-tendon units are undergoing their stretch-shortening cycle, in heavier animals a greater role is played by more elastic tendons at the expense of less elastic muscles. Pontzer 102 also argues for a graded difference across body mass in the mechanics underlying locomotion energy costs, but through a different pathway. He presents a model which indicates that force production becomes a smaller part of the cost, and mechanical work a greater part of the cost, as animals become larger.
In conclusion, while briefer ground-foot contact times may require faster, less energetically efficient muscle fibres, this mechanism appears insufficient to explain speed-and size-independent differences in the energy cost    of locomotion between species. Indeed, an interpretation of recent work by Gutmann and Bertram (2017) 103 is that if species of a similar size have very different stride frequencies at a given locomotion speed, as is apparent from data presented by 54 , then Kram and Taylor's (1990) relationship between metabolic rate and foot contact time should not hold. Instead, the volume of muscle activated per step 44 , variations in effective mechanical advantage for instance due to posture 102 , muscle contraction efficiency in terms of mechanical power 104 , and/or bigger animals having a lower energy loss by hysteresis during each step 101 , may explain an important proportion of the variation in residual movement costs. Sufficient information across species on duty factor, limb length, musculature and moment arms along with kinematic information will be necessary to investigate the covariance of these factors with the energy cost of locomotion across species. Alternatively, the considerable diversity in the species represented in inter-specific scaling relationships of NCOT perhaps hints that, for any given species, there could be a plethora of reasons explaining its relatively high or low costs of transport and, furthermore, these could be idiosyncratic. Consequently, mechanisms influencing the energy cost of locomotion beyond mass and, for example, effective limb length may not be generalizable.
Ecological correlates with NCOT. For many terrestrial species, the amount of energy they expend to locomote substantially decreases the energy remaining that they can channel into reproduction [105][106][107] . Thus animals are expected to optimise the efficiency of their energy expenditure during movement, within certain constraints 108 . In turn, we might predict correlations between species ecology and energy costs of transport. Yet our analysis on mammals found no relationships between NCOT and a host of ecological variables once body mass, which accounts for the fact that bigger animals have a lower NCOT and can range further 102 , was accounted for. Effectively, our analyses show that for different species of the same size, home range size, daily movement distance, geographic range, group size, terrestriality, diet breadth, trophic level and habitat breadth do not relate to NCOT. Furthermore, despite the fact that animals which can run faster or more energetically cheaply perhaps have the opportunity to roam more extensively and those with greater energy output might be predicted to carry greater energy stores to buffer against short falls 109 , neither maximum running speed nor percentage body fat levels were related to NCOT. Finally, NCOT was also not related to FMR, MMR, AMR or aerobic scope. Thus NCOT appears to be disconnected from animal foraging behaviour, broader measures of energetics and varying aspects of ecology. This either calls into question the importance of locomotion energy expenditure in influencing an animal's ecology, or the relevance of NCOT for understanding the locomotion energetics of terrestrial animals beyond the realms of biomechanics or comparisons among many animals of varying size. If either or both possibilities are true, in turn it is important to understand why NCOT has little ecological relevance. Pontzer 110 argues that most extant terrestrial animals may have already evolved to be efficient foragers. Given the adaptive losses that these species would experience as a consequence of further enhancements to their locomotion efficiencies 108 , for example longer legs might reduce an animal's capacity to accelerate 97 , there is little such selective pressure. Harris and Steudel 111 report that details of prey pursuit and capture are the factors that describe hind limb length, in contrast finding no predictive power in home range size or daily movement distance. Perhaps other physical factors that influence locomotion energy efficiencies, such as leg muscle mass, for similar reasons are also under minimal selection for enhancement in terms of energy efficiencies, particularly in species for which locomotion costs are a relatively small proportion of their total energy expenditure 112 . The possibility that NCOT is simply a poor representation of an animal's energy costs of locomotion is particularly supported in the present study by the lack of relationship between NCOT and field metabolic rate or activity metabolic rate.  Table 6. Outputs from phylogenetic generalized least squares (PGLS) analysis models to predict log 10 (NCOT), mlO 2 m −1 . For all models, body mass was included as a covariate. For all models the maximum likelihood value of λ was 0. DMD = daily movement distance; BMR = basal metabolic rate; FMR = field metabolic rate; MMR = maximum metabolic rate; AMR = activity metabolic rate; FAS = factorial aerobic scope; AAS = absolute aerobic scope.
SCiEntifiC REPORTs | (2019) 9:656 | DOI:10.1038/s41598-018-36565-z NCOT, calculated as the slope of the linear fit between rate of energy expenditure and locomotion speed, tends towards the total energy costs of locomotion per unit distance at high running speeds 33 . However, the majority of movements by animals are conducted at low speeds relative to their maximum obtainable [113][114][115][116] . Under these circumstances, per unit distance the cost of transport is higher than indicated by NCOT due to fixed costs associated with locomotion 33 , most notably the 'postural cost of transport' PCOT 72 . Furthermore, animals often incur additional energy costs while walking associated with, for example, intermittent locomotion 117 , turning 118 , and negotiating various terrains e.g. Fig. 4B in 119 . NCOT therefore underestimates an animal's true energy costs to move, and the magnitude of this error presumably varies between species depending upon factors such as their typical running speeds and overall movement behaviours.
To progress our understanding of the energetic costs and constraints for terrestrial animals traversing their habitats, if and how these costs help shape their ecology, and how those costs relate to their locomotion characteristics, perhaps now it is time to step off the laboratory treadmill and measure transport costs in the wild 33,120-122 .

Data Availability
We have made the data collated for the present study available on Dryad.