Edited by Nicholas D. Pyenson, Smithsonian Institution, Washington, DC, and accepted by Editorial Board Member David Jablonski February 23, 2018 (received for review August 8, 2017)
The reasons why aquatic mammals exhibit larger average sizes than their terrestrial relatives have long been debated. Most previous hypotheses have focused on releases from terrestrial constraints on large sizes. Through the analysis of mammal size distributions, we find the aquatic realm imposes stronger constraints on body size than does the terrestrial realm, driving and confining aquatic mammals to larger sizes. Calculations of energy intake and demand as a function of body size indicate heat loss imposes a strong lower bound on size, whereas the scaling of feeding rate versus metabolic rate imposes a constraint on maximum size. Rather than freeing animals from body size constraints, living in water appears to impose stronger selective pressures than does living on land.
Four extant lineages of mammals have invaded and diversified in the water: Sirenia, Cetacea, Pinnipedia, and Lutrinae. Most of these aquatic clades are larger bodied, on average, than their closest land-dwelling relatives, but the extent to which potential ecological, biomechanical, and physiological controls contributed to this pattern remains untested quantitatively. Here, we use previously published data on the body masses of 3,859 living and 2,999 fossil mammal species to examine the evolutionary trajectories of body size in aquatic mammals through both comparative phylogenetic analysis and examination of the fossil record. Both methods indicate that the evolution of an aquatic lifestyle is driving three of the four extant aquatic mammal clades toward a size attractor at ∼500 kg. The existence of this body size attractor and the relatively rapid selection toward, and limited deviation from, this attractor rule out most hypothesized drivers of size increase. These three independent body size increases and a shared aquatic optimum size are consistent with control by differences in the scaling of energetic intake and cost functions with body size between the terrestrial and aquatic realms. Under this energetic model, thermoregulatory costs constrain minimum size, whereas limitations on feeding efficiency constrain maximum size. The optimum size occurs at an intermediate value where thermoregulatory costs are low but feeding efficiency remains high. Rather than being released from size pressures, water-dwelling mammals are driven and confined to larger body sizes by the strict energetic demands of the aquatic medium.
Most mammal species live on land, but the largest mammals inhabit the ocean (1). The same is true of reptiles and arthropods (2⇓–4). Aquatic and terrestrial habitats clearly impose differing selective pressures on body size. However, the quantitative study of body size evolution in mammals and other diverse animal clades typically focuses on either terrestrial or aquatic clades independently (e.g., refs. 5⇓⇓⇓–9). Consequently, the extent to which the rate, magnitude, and outcome of size change associated with habitat transitions are shared among clades remains unknown. This gap in knowledge leaves open the question of whether the apparently common phenomenon of size increase associated with the adoption of an aquatic lifestyle reflects idiosyncratic responses of individual clades versus a common response to shared constraints.
Identifying the causes of size increase in aquatic mammalian clades requires information not only on the differences in current mean size between terrestrial and aquatic sister groups but also on the variance about these differing means and the value of the mean size, given sufficient evolutionary time for the size distribution to reach equilibrium. Here, we quantify these parameters for four independent terrestrial-aquatic transitions in mammals within a phylogenetic context to quantify how the shift to aquatic habitat has impacted the rate, magnitude, and directionality of body size evolution.
Several hypotheses have been put forward to explain why mammals living in water are larger, on average, than their closest terrestrial relatives, including theories regarding diet, neutral buoyancy, habitat area, protein availability, and thermoregulation (10⇓⇓⇓⇓⇓⇓⇓⇓⇓⇓–21). Although each of these theories predicts that mammalian clades entering the water will increase in average size, these theories differ in their predictions for how such a size change is achieved. More specifically, they differ in their predictions both about the rate of evolution toward the new, larger average size as well as the variance of the aquatic size distribution relative to its terrestrial sister group (22).
The Ornstein–Uhlenbeck (OU) process (23) represents a model of adaptive evolution appropriate to the question of size evolution across a transition in habitat, where the change of body mass, X, over time, t, is the result of the interactions between the body mass optimum, θ, the strength of selection toward the optimum, α, the intensity of random drift away from the optimum, σ, and white noise, B:
Idealized theoretical scenarios in which body mass optimum increases. Nine different scenarios are illustrated in which the optimal body mass of a clade increases when it transitions from the terrestrial regime to the aquatic regime. These scenarios represent combinations of discrete threefold increases and/or decreases in the phylogenetic half-life (PHL) and stationary variance (SV) upon a clade entering the aquatic regime. In all scenarios, the optimum, PHL, and SV of the terrestrial regime are 3.5, 7.0, and 0.5, respectively, and the optimum of the aquatic regime is 5.5. A, D, and G show an increase in PHF; B, E, and F show no change in PHL; and C, F, and I show a decrease in PHL. A–C show a decrease in SV; D–F show no change in SV; and G–I show an increase in SV. The protein availability, neutral buoyancy, and habitat area hypotheses predict any of G, H, or I. The thermoregulation hypothesis predicts B, C, E, or F.
Each of the previously proposed hypotheses for increase in mean size maps onto only a subset of these nine possibilities, providing a means of testing among potential driving mechanisms and even for identifying contributions from multiple processes. The diet hypothesis predicts that herbivores are often larger than closely related carnivores due to a greater metabolic demand (8, 25). This potential effect of diet is particularly important in clades where the aquatic members have changed feeding mode. If, for example, a group of terrestrial carnivores evolved into aquatic herbivores, this hypothesis would predict an increase in body size based on the change in feeding mode but no particular change in stationary variance or phylogenetic half-life (Fig. 1E). The neutral buoyancy hypothesis predicts a release from terrestrial constraints that should result in an increase in the stationary variance (11⇓⇓⇓⇓–16) (Fig. 1 G–I). The protein availability hypothesis similarly predicts a release from energetic constraint on maximum size, enabling larger maximum sizes in the oceans, and ultimately increasing the stationary variance (17, 18, 24) (Fig. 1 G–I). The habitat area hypothesis predicts that larger habitats in the aquatic realm will allow mammals to expand to larger sizes, which also results in an increase in stationary variance (8, 20, 21, 25⇓⇓–28) (Fig. 1 G–I). Finally, the thermoregulation hypothesis predicts that terrestrial mammals that enter the water should encounter selection toward larger body sizes, not merely a release from smaller sizes, which should result in an equal or decreased phylogenetic half-life and/or an equal or decreased stationary variance (10) (Fig. 1 B, C, E, and F).
Here, we use previously compiled body size data for extant and fossil mammal species to perform phylogenetic comparative analyses and fossil time series analyses to estimate these body size dynamics across the four evolutionary transitions between land and water within modern mammals. We then use these estimates to test among the various previously mooted hypotheses for the large sizes of aquatic mammals.
The best-supported phylogenetic models for size evolution in Afrotheria, Artiodactyla, and Caniformia indicate separate size optima for terrestrial versus aquatic habitats (Dataset S1). For Musteloidea, the best-fitting models treat body size evolution as simple Brownian motion, with no trend toward an optimal size (Dataset S1).
Pairwise tests between all of the estimated body mass optima indicate that they are statistically distinct from one another (Mann–Whitney test, P < 0.001). However, despite the substantial differences between the terrestrial optima of Afrotheria, Artiodactyla, and Caniformia, the estimated aquatic body mass optima for these three groups with evidence for separate aquatic body mass optima are much closer to one another (Afrotheria: 5.48, Artiodactyla: 5.60, and Caniformia: 5.38) than they are to any of their respective terrestrial optima (Afrotheria: 3.18, Artiodactyla: 4.82, and Caniformia: 3.95), implying convergence on a body mass attractor near 500 kg (Fig. 2A). Furthermore, the optimum for Mysticeti is higher than that of any of the other analyzed groups. The estimated aquatic optimum for Musteloidea does not mirror the dramatic increase of the other three clades (Fig. 2A). This sustained small size of aquatic mustelids could indicate the presence of a second attractor at a smaller mass of roughly 10 kg or competitive exclusion from the 500-kg attractor. One potential interpretation is that this finding reflects the semiaquatic nature of most otters, which causes them to be under a combination of aquatic and terrestrial selective forces. Testing this possibility is not feasible, however, because there exists only one fully aquatic otter species.
Plots of estimated body mass optima (A), phylogenetic half-lives (B), and stationary variances (C) by clade and habitat type, illustrating the statistically larger optimal body masses, equal or decreased phylogenetic half-lives, and equal or decreased stationary variance for aquatic mammals relative to their terrestrial sister groups. Optima (θ) are reported as weighted means with 2σ confidence intervals of model-averaged values (●) and are accompanied by median modern body masses (▲). Phylogenetic half-lives [ln(2)/α] and stationary variances (σ2/2α) are reported as box plots of model-averaged values. Color coding refers to habitat regime (orange, terrestrial; dark blue, toothed aquatic; light blue, baleen aquatic). Symbols at the bottom of the plots represent significance of results of pairwise Mann–Whitney tests (0 < *** < 0.001 < ** < 0.01 < * < 0.05 < . < 0.1).
OU models fit to the fossil time series of aquatic members of each of the four clades of interest produce optimal body mass estimates that overlap with the estimates of the phylogenetic modeling (Fig. S3 and Dataset S1). However, the confidence intervals for the paleoTS (Analyze Paleontological Time-Series) estimates span multiple orders of magnitude for Sirenia and Pinnipedia. These large confidence intervals are most likely due to the very linear trends of these two groups (Fig. S3), causing poor inference of the future plateaus implied by an OU model. The use of binned data for the paleoTS analyses may also be contributing to this uncertainty due to the loss of time resolution. Furthermore, while the use of fossil specimens and measurements may give a more complete picture of body size changes through time, the paleoTS software is designed specifically for use with ancestor-descendant sequences and accounts for the auto-correlation of adjacent points in the time series but does not account for the varied dependence of measurements to one another in a single time bin or across time bins (29, 30). This factor could also help to explain the decrease in statistical power when using paleoTS.
Phylogenetic half-life and stationary variance estimates for the aquatic regimes were consistently less than or equal to those of the terrestrial regimes across all four clades (Fig. 2 B and C, corresponding to Fig. 1 B, C, E, and F). Together, these parameter estimates indicate equivalent or decreased body size disparity and equivalent or increased strength of selection upon the transition into the aquatic habitat. The convergence of these three aquatic mammal clades in such a manner suggests that living in the aquatic realm, more than merely enabling large size, imposes a strong selective pressure toward larger body sizes while simultaneously imposing tighter constraints on size variation than those present in the terrestrial realm, potentially for reasons that are shared across clades.
Most previously proposed explanations for mammal size increase in aquatic habitats are either incompatible with the results of this analysis or at least incomplete in explaining the pattern. Hypotheses regarding neutral buoyancy, increased protein availability, and increased habitat area all predict a release from terrestrial constraints, which should result in an increase in the stationary variance (8, 11⇓⇓⇓⇓⇓⇓–18, 20, 21, 24⇓–28) (Fig. 1 G–I). However, the results above show no support for increased stationary variance upon entering the aquatic realm. The diet hypothesis predicts that herbivores are often larger than closely related carnivores (8, 25). For the optimal size of a group to increase, this would require some or all of the aquatic members to convert from carnivory to herbivory. However, most aquatic mammals are carnivores, and the only herbivorous aquatic mammals evolved from herbivorous terrestrial ancestors. Therefore, diet shift does not serve as a unifying explanation for the convergence of these three groups. Finally, the thermoregulatory differences between land and water (10) should cause selection toward larger body sizes, not merely a release from smaller sizes, which would manifest as equal or decreased phylogenetic half-lives and/or stationary variance (Fig. 1 B, C, E, and F), both of which are indicated by our results. Therefore, selective pressures from thermoregulatory requirements can explain convergent increases in minimum size across these groups. However, the thermoregulatory hypothesis only partially explains the observed patterns. Thermoregulatory costs scale roughly to the 1/4 power of body size and should most strongly impact the smallest end of the aquatic size distribution. They do not fully explain the optimum or maximum sizes of aquatic mammals.
An expansion of the thermoregulatory hypothesis toward a more general energetic explanation of mammal body size distributions helps to account for changes not only in minimum size but also in the calculated optimum and stationary variance as well as the observed maximum size and skewness of the aquatic size distribution. This expanded energetic model is similar to the models of energy use as a function of body size that have been applied to explain size distributions of other groups (31⇓⇓⇓–35). In this model, surplus energy, E, for an organism is the energy used for growth, reproduction, and behaviors that demand energy beyond basal metabolism. We take this to be the difference between the intake of chemical energy by feeding, F, and the expenditure of this energy by basal metabolism, M, and heat loss to the environment, H:
Plot of a plausible energetics model (B) and its components (A) that predicts the estimated optimal body masses and modern distributions of aquatic mammals. Components of energetics model with respect to body mass (A) and comparison of energetics model calculation with stacked modern aquatic mammal body mass distributions and estimated OUwie optima for Sirenia, Odontoceti, and Pinnipedia, respectively (B), are shown. The surplus energy curve predicts the OUwie optima and modern distributions.
This energetic model for fitness as a function of size in marine mammals successfully explains not only the increase in mean size but also the observed increase in constraint on size relative to terrestrial sister groups. Due to the increased energetic cost of living in water below body temperature (the H component), the minimum possible size of an aquatic endotherm is more than three orders of magnitude larger than the minimum possible size of an otherwise similar endotherm on land (10). Therefore, terrestrial mammals that enter the water encounter selection toward larger body sizes, not merely a release from smaller sizes. In addition, basal metabolism exceeds the rate of feeding at larger sizes due to differences in allometric scaling exponents, imposing a selective pressure from energy metabolism against size increase at larger body masses. The similarities in the size distributions among toothed aquatic mammals suggest that among-clade differences in the mass-feeding relationship are comparatively minor.
Changes in the slope and/or intercept of the relationship between body mass and feeding rate associated with the evolution of baleen can account for the shift toward even larger sizes in baleen whales. This shift in the energetic constraints on size evolution may be the result of an increase in feeding efficiency in baleen whales compared with their toothed relatives, due to the evolution of new foraging strategies, particularly filter-feeding and lunge feeding (41⇓–43). For example, modifying the feeding function presented in Eq. 5 to an exponent of 0.78 (versus 0.71) results in a surplus energy distribution that maintains a peak around the optima of the three toothed aquatic groups but also includes the largest sizes of the baleen whales (Fig. S4). This hypothetical relationship does not explain the absence of extant baleen whales at smaller sizes, suggesting that this feeding mode also leads to poor competitiveness with toothed relatives (or other taxa) at smaller sizes or modifies the scaling of feeding with respect to body mass in more complex ways than addressed by the model above. The fact that extreme gigantism in baleen whales is a fairly recent phenomenon, within the last 10 My, relative to the evolution of baleen further suggests that the amount and spatial distribution of food resources further determine the viability of this feeding mechanism (44, 45).
Through the combination of heat loss and feeding constraints, an energetic model for body size evolution in aquatic mammals predicts equal or decreased phylogenetic half-lives and stationary variances on evolutionary time scales (Fig. 1 B, C, E, and F), both of which are indicated by the distributions of body sizes across living and fossil aquatic mammals. Because all aquatic mammals share these energetic constraints, they can explain the rapid, convergent evolution toward larger optimal body sizes across toothed Sirenia, Cetacea, and Pinnipedia, despite differences in the body size attractors displayed by their terrestrial sister groups.
Whereas the aquatic medium has been argued to release animals from constraints on body size through mechanisms such as habitat area and neutral buoyancy (10⇓⇓⇓⇓⇓⇓⇓⇓⇓⇓–21), comparative analysis of fossil and living mammals demonstrates that the aquatic realm actually imposes tighter constraints on size evolution in aquatic mammals than in their land-dwelling relatives. Three of the four mammal clades living in aquatic environments have evolved toward larger optimal body masses than their terrestrial counterparts. Moreover, these three clades have evolved independently toward a shared body mass attractor of ∼500 kg, coupled with increased rates of evolution and decreased variance. Many previous hypotheses for mammalian body size increase in aquatic habitats fail to account for these additional body size dynamics. By contrast, an energetic cost model incorporating size-dependent functions for food intake and costs from basic energy metabolism and thermoregulation accounts for the observed combination of increase in body size, decrease in variance, and increase in evolutionary rate. Parameterizing the energetic model using data from phocids also predicts the observed lower and upper bounds of the distribution as well as the observed right skewness of the distribution. Therefore, the convergence of aquatic mammal clades on a body mass optimum near 500 kg appears to reflect the shared constraints of allometric scaling of thermoregulation, basal metabolic rate, and feeding potential with body size for toothed mammals. Rather than being released from size pressures, they are driven and confined to larger body sizes by the strict demands of the aquatic medium.
We compiled body masses for 3,859 extant and 2,999 fossil mammal species from existing datasets (5⇓–7, 46⇓⇓⇓⇓⇓⇓–53; www.helsinki.fi/science/now/) that, in aggregate, cover ∼70% of extant mammalian species and 25% of accepted extinct species (Figs. S1 and S2). We replaced synonymous taxonomic names with the accepted canonical names using the Global Biodiversity Information Facility (GBIF) Backbone Taxonomy (https://www.gbif.org/), accessed using the R package rgbif (54, 55). When multiple measurements existed, we calculated the mean and SE. When only a single measurement existed, we randomly sampled a SE from the calculated SEs of the other members of the species’ family or, if none existed, of the species’ order.
We combined binary habitat assignments (coded as marine/nonmarine) from the GBIF and categorical habitat assignments (coded as terrestrial/marine/brackish/freshwater) from the World Register of Marine Species to create a habitat assignment for each species coded as aquatic/nonaquatic (https://www.gbif.org/; www.marinespecies.org). We resolved ambiguous habitats for fossil species using “life habit” assignments from the Paleobiology Database and the primary literature. Ultimately, we excluded 273 species from further analysis because their habitat affinities remained ambiguous after these searches.
We downloaded the stratigraphic ranges of all mammal species with fossil records using the Paleobiology Database Data Service v.1.1 on November 8, 2016 (https://paleobiodb.org/data1.1/taxa/list.txt?name=Mammalia&rel=all_children&show=app&rank=species&status=senior&limit=15000). We excluded all species that lacked a stage-resolved stratigraphic range. This culling procedure yielded a set of 3,099 species that could be used for time series analyses.
We used the time-calibrated supertree of extant mammalian species for all analyses requiring a phylogenetic framework (56, 57). We used a sample of 100 trees from the posterior distribution of a Bayesian birth-death polytomy resolver for all analyses (58). Although this approach is not a perfect substitute for a fully resolved tree (59), it allows us to randomly sample over the different resolved forms of the supertree for the analyses discussed below. The supertree includes 3,832 of the species in our database.
We tested the influence of habitat on body size evolution using generalized OU process modeling (23, 60, 61) (Eqs. 1–3), as has been used in previous studies (25, 62). These OU models contain three parameters describing the evolution of body size over time: the primary optimum of body size (θ), the rate of stochastic evolution away from the optimum (σ2), and the strength of selection toward the optimum (α). The hypothetical trait value that distinguishes stabilizing and directional selection from Brownian motion is θ. This value is assumed to be constant through time for a given selective regime and represents the average value reached by an infinite number of populations descended from a single ancestor, also referred to herein as the attractor. The parameter σ2 is a constant measure of the stochastic evolution within a selective regime. Under a Brownian motion model (where α = 0), σ2 is also known as the Brownian motion rate parameter. The parameter α is a measure of how quickly a lineage within a given selective regime reaches θ. If α is large, the selective pressure of living in a given selective regime has a much stronger effect on the trait value than the phylogenetic history, whereas if α is small, the phylogenetic constraints are relatively strong. Together, the parameters σ2 and α determine the trait disparity within a selective regime.
We performed all tree manipulations and statistical analyses using the R software environment (54). Using the R package OUwie (63), we fit seven different OU models in which one, two, or all three of these parameters were estimated for the entirety of a clade and separately for members of the two habitat types, terrestrial and aquatic. In each analysis, these parameters are optimized to best fit the given model’s predictions to the known body sizes of the extant species while also taking into account the known measurement error. The BM1 model fits a single Brownian motion model to the entire clade, estimating only the σ2 parameter. The BMS model fits separate Brownian motion models to terrestrial and aquatic members within the clade, estimating σ2 separately for each selective regime. The OU1 model fits a single OU model to the entire clade, estimating the single values for the θ, σ2, and α parameters over the entire phylogeny. The OUM model, in addition to estimating σ2 and α for the entire clade, estimates θ separately for each selective regime. The OUMV model builds on the OUM model, allowing σ2 to vary across selective regimes but holding α constant across the entire clade. Meanwhile, the OUMA model allows α to vary across selective regimes but holds σ2 constant across the entire clade. The most complex model, OUMVA, allows all three parameters to vary between terrestrial and aquatic selective regimes. It is assumed that all parameter space is theoretically biologically valid.
We fit these models across the Afrotheria, Artiodactyla, Caniformia, and Musteloidea clades using a set of 100 resolved Mammalia phylogenies (discussed above) to test the independent evolutions of aquatic life habit in mammals. In this modeling process, the time-scaled phylogeny is used to estimate the amount of time that exists since the divergence of a set of species (t in Eqs. 1–3). Lutrinae (otters) were excluded from the Caniformia analyses because they are the aquatic clade of interest in the Musteloidea analyses. The ancestral habitat conditions required by the OUwie R package and function were reconstructed using the ace function from the R package ape (64). The state with the highest marginal likelihood was assigned at each node. The ape package was also used for tree manipulation (e.g., tip removal, clade isolation). Through preliminary analyses, we found that dropping the estimated root parameter (θ0) (i.e., setting root.station = TRUE) within the OUwie function helped to stabilize the parameter estimates. With positive alpha values, the weight of the information in the analysis becomes skewed toward the present and the root estimate carries very little weight, causing it to be a difficult parameter to estimate and causing stabilization problems downstream. Dropping this parameter therefore assumes that the root was within the distribution of the ancestral regime (in this case, the terrestrial regime in all analyses: θ0 = θt). Furthermore, due to the key innovations of baleen and ram feeding in baleen whales and the definite impact this evolutionary innovation has on their body size (41), we have separated Cetacea into two regimes, corresponding to Mysticeti and Odontoceti. Including this split in the analyses resulted in a very small decrease in the optimum for toothed whales and improved likelihood support compared with a two-regime model.
We determined the relative support for each model for each tree using weighted Akaike information criterion (AIC) values, corrected for sample sizes, hereafter referred to as AICc weights (65, 66). In all AICc weight calculations, the sample size is taken to be the total number of tips within the clade of interest (61). To account for the various degrees of support for the different models and for the inclusion of different parameters in each of the models, we calculated model-averaged parameter estimates using the formulae presented by Burnham and Anderson (65). Because the θ parameter estimates for the BM1 and BMS models represent the estimated body size at the root, not necessarily the optimal body size, these models were not included in the model-averaging process.
To compare these phylogenetic (extant) results with fossil data, we used the fossil mammal ranges in our database to compute time series consisting of species body size means and variances for 1-My time bins for Sirenia, Pinnipedia, Odontoceti, and Lutrinae. We then fit OU models to these time series using the R package paleoTS (67⇓–69). We used the joint parameterization methods in the log-likelihood function to account for auto-correlation among samples (29).
We thank J. Beaulieu and G. Hunt for insightful conversations about their software and for accommodating feature requests. We thank S. K. Lyons for sharing fossil mammal body size estimates. We are grateful to two anonymous reviewers for providing very helpful comments on earlier versions of this manuscript.
Author contributions: W.G., C.R.M., and J.L.P. designed research; W.G. performed research; W.G. analyzed data; and W.G., C.R.M., and J.L.P. wrote the paper.
The authors declare no conflict of interest.
This article is a PNAS Direct Submission. N.D.P. is a guest editor invited by the Editorial Board.
Data deposition: The code and data for all analyses are available at BitBucket (https://bitbucket.org/wgearty/body-size-analyses/).
See Commentary on page 3995.
This article contains supporting information online at www.pnas.org/lookup/suppl/doi:10.1073/pnas.1712629115/-/DCSupplemental.
Published under the PNAS license.
See related content: