Independent evolution of baleen whale gigantism linked to Plio-Pleistocene ocean dynamics

Vertebrates have evolved to gigantic sizes repeatedly over the past 250 Myr, reaching their extreme in today's baleen whales (Mysticeti). Hypotheses for the evolution of exceptionally large size in mysticetes range from niche partitioning to predator avoidance, but there has been no quantitative examination of body size evolutionary dynamics in this clade and it remains unclear when, why or how gigantism evolved. By fitting phylogenetic macroevolutionary models to a dataset consisting of living and extinct species, we show that mysticetes underwent a clade-wide shift in their mode of body size evolution during the Plio-Pleistocene. This transition, from Brownian motion-like dynamics to a trended random walk towards larger size, is temporally linked to the onset of seasonally intensified upwelling along coastal ecosystems. High prey densities resulting from wind-driven upwelling, rather than abundant resources alone, are the primary determinant of efficient foraging in extant mysticetes and Late Pliocene changes in ocean dynamics may have provided an ecological pathway to gigantism in multiple independent lineages.


Introduction
Vertebrates have evolved to gigantic sizes repeatedly over the past 250 Myr [1,2], reaching their extreme in today's blue whale, Balaenoptera musculus, which is the largest animal to have ever lived. All living baleen whales (Mysticeti), including the blue whale, are obligate suspension feeders and they possess a complex suite of adaptations that enhance the energetic efficiency of foraging on small-bodied, low trophic-level prey [3][4][5]. That such large animals should feed on such small prey is not without precedent; the fossil record demonstrates that large-bodied suspension feeders have arisen on several occasions in diverse clades [1,6,7]. However, the ecological mechanisms and evolutionary processes that promote and maintain gigantism remain poorly understood in general [1], and we lack a comprehensive understanding of how and when mysticetes, in particular, attained such large sizes.
Macroevolutionary analyses based on phylogenies of extant cetaceans have suggested a number of possible explanations for the evolution of mysticete gigantism, including diet-related niche partitioning among the earliest representatives of crown cetacean clades [8], rapid rates of body size evolution along the mysticete stem [9] or as one expected outcome of trade-offs between the short-term fitness benefits of size increases versus long-term costs (i.e. increased extinction risk) associated with large size [10]. While most of these hypotheses predict that gigantism should have evolved relatively early in the clade's history, as in terrestrial mammals [11], the mysticete fossil record suggests a much later origin of exceptionally large size [2,[12][13][14][15] but see [16]. From their origin in the Late Eocene or Early Oligocene through the Middle Miocene, the largest mysticetes remained less than 10 m long, though several lineages appear to have independently explored the upper reaches of this size spectrum [16][17][18]. Lambert et al. [13] argued that true gigantism, defined as body lengths larger than 10 m, evolved in the early Late Miocene in response to the evolution of large macro-predatory physeteroid odontocetes and lamniform sharks. Others have suggested a later, Plio-Pleistocene origin for gigantism, either as a direct response to increased near-shore primary productivity from the Late Miocene onwards [12] or else owing to the effects of glacial cycles on habitat availability and resource distributions [14]. Despite a rich fossil record that is ideally suited for macroevolutionary inference [19], there have been no quantitative tests of when and how mysticetes achieved gigantic sizes. As such, it remains to be determined whether the evolution of blue whale-sized animals requires special explanation or whether it is simply one plausible outcome of a stochastic, constant-rates process [20].
Here, we leverage the excellent fossil record of mysticetes and a robust phylogenetic framework to provide, to our knowledge, the first formal test of when and how gigantism evolved. Using novel phylogenetic models of body size evolution and extensive simulation, we show that the evolution of exceptionally large size (more than 10 m) is a recent phenomenon that results from a fundamental clade-wide shift in the mode of body size evolution.

Material and methods (a) Body size data
Studies of mammalian body size evolution typically focus on body mass, given the well-defined anatomical and physiological responses of this trait to abiotic factors [11]. Reliable estimates of body mass are rare for cetaceans but available data indicate that total length (TL) scales with mass 1/3 [21]. We therefore used log 10 (TL) as the metric for our analyses. We took a conservative approach to defining lengths for extant species by collecting length data from museum specimens, stranding records, aerial surveys and aboriginal subsistence harvests that, unlike commercial whaling records, are not subject to minimum length restrictions [22] and so do not lead to a potentially false distinction between extant and extinct body size distributions. Data collection was restricted to adult individuals, and resulted in an average n of 13.3 individuals per extant species (range ¼ 1-57; electronic supplementary material, table S1).
For extinct mysticetes, we estimated log 10 (TL) in millimetres as log 10 (TL) ¼ 0:92 Ã [log 10 (bizyg) À 1:72] þ 2:68, ð2:1Þ where bizyg is the bizygomatic breadth of the skull, measured in millimetres [23], which we obtained either from direct measurement of fossil specimens (by N.D.P.) or from the literature. Because this approach requires well-preserved crania, we were only able to obtain a single size estimate for each fossil species.
We computed species means for all species where n . 1. To account for measurement error in comparative model fitting (see below), we computed a pooled variance over all species represented by multiple specimens. The standard error of the mean for the ith species, including fossil taxa, was then computed as the pooled standard deviation divided by the square root of n i .

(b) Phylogenetic inference
As a framework for macroevolutionary inference, we jointly estimated topology and branch lengths, in millions of years, of mysticete phylogeny from morphological and molecular data using BEAST v. 2.2.1 [24], accessed through the CIPRES Science Gateway [25]. Morphological character data for 13 extant and 63 extinct mysticetes were based on [14]. We also downloaded 11 nuclear loci and coding regions of mitochondrial genomes, where available, for all 15 extant mysticete species from Genbank (electronic supplementary material, table S1). We used the fossilized birth -death (FBD) process [26] as a prior on the distribution of branching times and branch lengths while allowing for potential ancestor-descendant relationships [27,28]. Complete details of phylogenetic analyses are provided in the electronic supplementary material.

(c) Tempo and mode of body size evolution
We first computed two measures of phylogenetic signal (Pagel's l [29] and Blomberg's K [30]) for log 10 (TL), using PHYTOOLS v. 0.5-38 [31] for R v. 3.3.1 [32]. We also evaluated trends in subclade disparity through time (DTT) [33] using a modified version of code from GEIGER v. 2.0.6 [34]. We compared mysticete DTT to a null expectation derived from 10 000 Brownian motion (BM) simulations, and computed the morphological disparity index (MDI) as the area between the median of these simulations and the observed curve.
To explicitly test hypotheses for the evolution of mysticete gigantism, we fitted a series of macroevolutionary models to our comparative dataset using maximum likelihood (ML). We fitted time-homogenous single rate BM, single peak Ornstein-Uhlenbeck (OU), biased random walk (also referred to as a trend model), and time-dependent rate (accelerating/decelerating: AC/DC) models of morphological evolution using the fitContinuous function in geiger. Variance of the means was added to the diagonal elements (variances) of the model-specific variance-covariance matrix during model fitting to account for measurement error. We used the OUwie.slice function in the OUwie package [35] to fit a variable rate model in which we allowed evolutionary rates to shift from one rate regime to another (either higher or lower) at some point in the past, with the shift point treated as a free parameter. We also tested for an effect of mean global ocean temperature, approximated as the d 18 0 curve of [36], on rates of body size evolution using the fit_t_env function [37] in the RPANDA package [38]. We generated a smooth cubic spline (d.f. ¼ 15, figure 1) from the d 18 0 data using the smooth.spline function in the stats package in R and allowed rates of size evolution to vary as a function of this curve. Use of higher degrees of freedom, allowing for a more detailed curve, did not change ML parameter estimates or likelihoods.
We finally considered an alternative model where the mode of size evolution shifts from unbiased BM to a trended random walk at some time, t shift , in the past. Under this model, the elements V ij of the phylogenetic variance-covariance matrix are identical to those of unbiased BM but elements in the vector of expected mean trait values E(x) are given by where x i is the expected value of the ith terminal taxon, t i is its occurrence time in millions of years before present, u is the root state and b is the trend parameter, which may be positive (size increases with time) or negative (size decreases). We treated u, b, t shift and evolutionary rate (s 2 ) as free parameters in our model. Simulation tests indicate appropriate false positive rates when fitting this model, and parameter identifiability is generally good (electronic supplementary material). Relative support for each model was assessed through computation of small sample corrected Akaike weights (w A ) and support for the best-fitting model over a timehomogeneous BM model was assessed via a parametric bootstrap using 1000 simulated datasets. We assessed the robustness of our results to topological and branch length variation by fitting all models to both the maximum clade credibility (MCC) tree and 1000 trees drawn at random from the post-burn in posterior sample.

(d) Effects of preservation bias
An artificial break between the body size distributions of extant and extinct mysticetes could arise if large-bodied pelagic taxa exhibit decreased fossilization or sampling probabilities. We rspb.royalsocietypublishing.org Proc. R. Soc. B 284: 20170546 used a simulation approach to determine the effects of such size-biased sampling on subsequent macroevolutionary model inference when the true model of evolution is an unbiased, constant rates process. Our simulations treated sampling probability as a logistic function of size, with sampling ranging from random (P[sampling] ¼ 0.5 for all taxa) to completely biased against large taxa (P[sampling] ¼ 0 for any fossil taxon larger than the clade mean). Full details of our simulation procedure and methods are provided in the electronic supplementary material.

(a) Body length
Extant baleen whales have a right-shifted body size distribution compared to their fossil relatives (figure 1). Notably, this is not simply owing to the largest extant mysticetes being larger than the largest fossil species; the smallest extant mysticete, Caperea marginata, is larger than the smallest fossil taxa and is more comparable to the average size of mysticetes for the Oligocene through to the Pliocene.

(b) Phylogenetic inference
The MCC tree topology is in general agreement with previous studies of both extant and fossil mysticete phylogeny [14], including the placement of Eschrictius and Megaptera within a paraphyletic Balaenoptera (figure 2). Divergence time estimates are somewhat younger than previous estimates derived from tip-dating [14], which probably results from use of the FBD tree prior [39]. Mapping body size onto the MCC tree topology (figure 2) confirms that increases in body size occur independently in multiple lineages of extant mysticetes.

(c) Tempo and mode of body size evolution
Phylogenetic signal in mysticete body size is high (Pagel's l ¼ 0.94, likelihood ratio test vs l ¼ 0 p , 0.001; Blomberg's K ¼ 0.47, p , 0.01). DTT analysis is mostly consistent with a constant rates process (MDI ¼20.063, p ¼ 0.65) but shows a pronounced pulse of increased within-clade variation at approximately 5 Ma that is inconsistent with a timehomogeneous evolutionary process ( p ¼ 0.066; figure 3).
None of the models allowing for time-dependent or temperature-dependent rates, temporal trends in the mean (e.g. Cope's rule), or a stationary optimal size provided a better fit to mysticete body size data than a simple Brownian diffusion model with a single evolutionary rate (Akaike Weight, w A , for BM ¼ 0.33; w A other models less than 0.26). However, we found strong support (w A ¼ 0.97) for the mode-shift model in which body size evolution switches from a slow unbiased random walk to an upward-biased random walk, with an inferred shift time of 0.19 Ma. Closer inspection of results revealed a saddle point in the likelihood surface for this model that is attributable to a lack of fossil data in the interval 3-0 Ma. Pleistocene glacial cycles and associated exposure and erosion of shelf area limit the preservation potential of Pleistocene cetacean fossils and little complete material is available for this interval [40]. Rather than attempting to add fragmentary Pleistocene fossil taxa to the tree, we instead reviewed the literature for fossil occurrences of extant species that could be used to truncate terminal branches and, in turn, create a Pleistocene pseudofossil record. Two fossil occurrences fulfilled our strict requirement of being identifiable to the species level and preserving sufficient morphology to confirm that they fall within the range of sizes exhibited by extant populations; a posterior cranium of a humpback whale (Megaptera novaeangliae) from the lower Kioroshi Fm Japan is dated to 0.125-0.15 Ma [41] and a grey whale (Eschrichtius robustus) skull and skeleton from the San Pedro Sand, California [42] is dated to 0.2-0.5 Ma [43]. Truncating the terminal branches for these two species to their youngest possible ages (0.125 and 0.2 Ma, respectively) resulted in a positive-definite Hessian matrix and increased support (w A ¼ 0.99; table 1) for a slightly older mode shift (0.31 Ma) with a 2-unit support range [44] of 0.13-4.5 Ma (figure 4).

(d) Effects of preservation bias
Using simulated data, we found no effect of size-biased sampling on false detection rates for the mode-shift model. This result holds even when assigning a sampling probability of zero to all fossil species larger than the clade's mean size (electronic supplementary material).  Figure 1. Mean body lengths for extant mysticetes and estimated length for fossil species (baleen-bearing mysticetes, circles; toothed mysticetes, triangles) are plotted according to their age as inferred from our phylogeny. Shaded areas correspond to 80 (white), 90 (grey) and 95% (black) quantiles of 1000 Brownian motion simulations on mysticete phylogeny and illustrate that the modern fauna is both lacking in small species (less than 5 m) and over-represented in large ones (more than 10 m), relative to the fossil record. To the right is a smooth-spline fitted to the Eocene -Present oxygen isotope curve [36], and used as a proxy in modelling temperature-dependent body size evolution.  productivity [12], and glacial disruption of near-shore habitat [14]. Additional plausible explanations can be formulated that would also be consistent with the qualitative pattern but that invoke vastly different timings or processes, such as a simple increase in variance towards the present as expected under a constant rates process [46,47] or biases in the fossil record that prevent the sampling of large pelagic taxa. By taking a phylogenetic approach to modelling body size evolution in the fossil record, we have shown that none of these explanations can explain the observed discrepancy between fossil and extant mysticete body sizes and, instead, identify a shift in evolutionary mode during the past 4.5 Myr that resulted in the largest animals to have ever lived.

(a) Time homogeneous processes and sampling bias cannot explain the late origin of gigantism
It is well known that an apparent trend towards size increase (or decrease) through time could equally be the outcome of an actual bias in the direction of size evolution or owing to an unbiased process in which variance increases through time [20,46,47]. Although mysticete body size shows high phylogenetic signal, consistent with a constant rates process [48], we found that a time homogeneous evolutionary mode is unlikely for mysticetes. Size disparity for the first 30 Myr of mysticete evolutionary history can be reasonably well approximated by simulating under BM, but the size distribution of extant species is incompatible with these expectations ( figure 1). This discordance is emphasized by dramatic deviations in subclade disparity through time plots at around 5 Ma ( figure 3). We similarly find no support for the hypothesis that taphonomic biases have filtered the mysticete fossil record to an extent that generates erroneous inference of a recent mode shift. Taphonomic size biases in the terrestrial realm typically remove small species, resulting in an over-representation of large-bodied taxa [49]. Habitat preferences play an equally influential role in the marine realm [50] and could bias against the preservation of large-bodied pelagic taxa, in turn generating the disjunct body size distributions for extant and extinct mysticetes that we recover. Live-dead comparisons of extant cetacean communities using stranding records suggest that the size spectrum of fossil communities is unlikely to have been strongly filtered by taphonomic biases [15], but our simulations show that even if such filtering did occur, it cannot explain the support we recover for a shift in evolutionary mode associated with Plio-Pleistocene gigantism (electronic supplementary material). Key to understanding this initially confusing result is the nature of the discordance between fossil and extant body size distributions; while the largest extant mysticetes are larger than expected under a constant rates process, the smallest extant species are also too large (figure 1), meaning that failure to sample large fossils is only half the problem. Although a complex suite of preservation biases may have colluded to generate the patterns we observe, we are unable to conceive of such a scenario at present. Taken together, these results argue strongly against stochastic explanations for the onset of gigantism and suggest a more mechanistic explanation is required.

(b) Timing and mechanisms
By confidently restricting the evolutionary mode shift to the Plio-Pleistocene, we are able to rule out many earlier proposed drivers of gigantism, such as the onset of an Antarctic Circumpolar Current, the initial evolution of bulk feeding or pressure from macro-predators. Because the maximum size of all described Palaeogene mysticetes is at the lower end of the Quaternary distribution [16] and well within the confidence envelope of a constant rates model, we also rule out competition among stem mysticete lineages [51] or other neocete clades [8] as a significant factor in the recent evolution of exceptional size. In fact, appending two recently described [17,18], large, stem mysticetes from the late Oligocene into our phylogeny and liberally assigning them total lengths of 10 m is insufficient to overturn strong support for a recent origin of gigantism (electronic supplementary material).
Although the Plio-Pleistocene witnessed a steep decline in average global temperatures [36], we find no support for temperature-dependent increases in macroevolutionary rate, contrary to predictions from macroecological studies [52,53] and macroevolutionary modelling of extant clades [37]. It is tempting to speculate that large size may have evolved as a directional physiological response to declining ocean temperatures rather than through increased rates that serve to increase size disparity, but the high phylogenetic signal present in our data is inconsistent with strong selection on a narrow adaptive peak [48,54]. What, then, explains the evolution of the largest animals in Earth's history?
Near-shore productivity began to increase in Late Miocene with the onset of coastal upwelling systems [55], and it seems no coincidence that mysticete palaeodiversity reached its peak during this interval [14,19,56], driven dominantly by the radiation of small-bodied and probably piscivorous [57] Cetotheriidae. Our estimated range of times for the origin of gigantism is more recent and we cannot rule out the hypothesis that disrupting effects of Pleistocene glacial cycles played a prominent role in driving the evolution of enormous size [14]. However, the world's oceans also experienced increasingly intensified, wind-driven upwelling dynamics from the Late Pliocene onwards [58][59][60]. These dynamics have been invoked as contributing to the explosive diversification of delphinine dolphins during the past 3.5 Myr [61] and are equally consistent with the timing of the shift in evolutionary mode we recover for mysticetes.
We hypothesize, based on an understanding of mysticete feeding mechanics and energenetics [5], that the origin of gigantism lies in this Late Pliocene shift in oceanography and concommitant changes in the intensity and the distribution of primary productivity. Because extant mysticetes are active bulk feeders on small-bodied prey suspended in the water column, it is high prey densities, rather than prey abundance per se, that increase the energetic efficiency of each feeding event [3,62]. Indeed, extended foraging bouts on highly dense and ephemeral resources are a defining characteristic of life history and feeding ecology in the largest living mysticetes, regardless of ram feeding versus lunge feeding behaviours or phylogenetic affinity [63,64]. Prior to the onset of modern upwelling regimes, the tropics appear to have been in a more or less permanent El Niñ o state [60], a condition that today is associated with reduced primary productivity and low densities of large mysticetes [65]. Although alkenone fluxes in sediments from the eastern equatorial Pacific suggest that productivity fluctuated with temperature during Pleistocene glacial cycles [59], the abrupt transition from homogeneous to heterogeneous, upwelling driven productivity patterns provides a direct mechanism to explain the selective advantage of large size.
The Plio-Pleistocene loss of small mysticetes (figure 1) can also be explained by these mechanisms. Our estimated timing for the origin of gigantism is coincident with a decline in global mysticete diversity from 3 Ma [14] that is driven dominantly by the near extinction of the diminutive Cetotheriidae [66], and of small, ram-feeding right whales (Balaenidae). Although intensified, wind-driven upwelling seasonally increases productivity and prey densities, the effects are localized to continental shelf breaks. Ichthyolith records from the South Pacific Gyre indicate simultaneous crashes in teleost and chondrichthyan productivity during the transition from a permanent El Niñ o state to modern oceanic conditions [67], further corroborating a redistribution of oceanic productivity. With the evolution of larger body sizes, mysticetes would have benefited from lower mass-specific metabolic rates to buffer against periods of low resource availability, lower costs of transport for efficient long distance migration between feeding locations, and high feeding performance when these ephemeral resources became available. At the same time, small-bodied lineages would probably have been out competed by their larger congeners that were more efficient at moving between and exploiting these patchily distributed ephemeral pulses of productivity in desert-like oligotrophic ocean ecosystems.
The lack of robustly identified Late Pliocene and Pleistocene mysticete fossils of extinct or extant species represents the final impediment to confirming our hypothesis. Deméré et al. [40] noted that the rich Pliocene mysticete record consists mostly of extinct but taxonomically uncertain species, but considered named fossil species from Pleistocene deposits to be nomina dubia that most likely represent extant taxa. Conversely, confirmed fossils of extant taxa are limited to a handful of species [41,42] and the evolutionary history of size in the largest living mysticetes remains enigmatic. Taxonomic and phylogenetic resolution of Pleistocene mysticete fossils is therefore a critical last step in confirming whether the origin of gigantism was a response to intensified upwelling and the associated restructuring of marine primary productivity during the Late Pliocene, and whether this shift in body size evolutionary dynamics occurred over a prolonged period or as a more or less instantaneous increase in size over multiple lineages.

Conclusion
By analysing the evolution of mysticete body size using phylogenetic approaches and palaeontological data, we confirm that gigantism in this clade is a surprisingly recent phenomenon. Although filter feeding using baleen had probably evolved by the Mid-Oligocene (25 Ma) [17,68], the ecological preyscapes that energetically favour gigantism only arose in the   Plio-Pleistocene, with the onset of seasonally intensified upwelling regimes (ca 3 Ma). As a result, we live in a time of giants; unlike any other time in geological history, modern oceans are rich with extremely large bodied suspension feeders that rely on dense but patchily distributed prey resources. Projected climate-driven changes to ocean ecosystem structure, diversity, and productivity presage a decrease in critical habitat for large-bodied baleen whales and other suspension feeding vertebrates [69], highlighting the ecological vulnerability of these giants operating on an energetic knife-edge.
Data accessibility. Additional methods and results supporting this article have been uploaded as part of the online electronic supplementary material. All datasets and scripts used in this study are available from the Dryad Digital Repository: http://dx.doi.org/10.5061/ dryad.b68g0 [70].