The dynamic adaptive landscape of cetacean body size

long-term


SUMMARY
Adaptive landscapes are central to evolutionary theory, forming a conceptual bridge between micro-and macroevolution. [1][2][3][4] Evolution by natural selection across an adaptive landscape should drive lineages toward fitness peaks, shaping the distribution of phenotypic variation within and among clades over evolutionary timescales. 5 The location and breadth of these peaks in phenotypic space can also evolve, 4 but whether phylogenetic comparative methods can detect such patterns has largely remained unexplored. 6 Here, we characterize the global and local adaptive landscape for total body length in cetaceans (whales, dolphins, and relatives), a trait that spans an order of magnitude, across their $ 53-million-year evolutionary history. Using phylogenetic comparative methods, we analyze shifts in long-term mean body length 7 and directional changes in average trait values 8 for 345 living and fossil cetacean taxa. Remarkably, we find that the global macroevolutionary adaptive landscape of cetacean body length is relatively flat, with very few peak shifts occurring after cetaceans entered the oceans. Local peaks are more numerous and manifest as trends along branches linked to specific adaptations. These results contrast with previous studies using only extant taxa, 9 highlighting the vital role of fossil data for understanding macroevolution. [10][11][12] Our results indicate that adaptive peaks are dynamic and are associated with subzones of local adaptations, creating moving targets for species adaptation. In addition, we identify limits in our ability to detect some evolutionary patterns and processes and suggest that multiple approaches are required to characterize complex hierarchical patterns of adaptation in deep time.

RESULTS
Total length data imputation By using phylogenetically informed imputations and a range of proxy measurements (see STAR Methods), along with total length information from the literature and museum specimens, we generated the largest total length dataset to date for Cetacea, comprising 345 taxa (89 living and 256 fossils; Figure 1). Where total lengths were known, imputed values showed a strong correlation with them, regardless of the amount of information used for estimation (see Figure S1), meaning that total lengths could be estimated for even incomplete fossil taxa ( Figure S1), allowing their inclusion in downstream analyses. These additional data are critical to robust macroevolutionary inference, though we note that it is impossible to calculate total length values without a degree of uncertainty (see STAR Methods), and this must be considered in any interpretation of the results.
Exploring tempo and mode in body size evolution As in previous work, 9,14 we find that disparity through time in extant cetaceans is lower than would be expected under a constant-rates process (center of gravity = 0.43, p = 0.005), consistent with early burst adaptive radiation models ( Figure 2). However, this signal is entirely erased by the inclusion of fossil taxa, regardless of whether we include all fossil cetaceans or limit our sample to Neoceti only (center of gravity = 0.29, p = 0.294).
The macroevolutionary adaptive landscape is often conceptualized as an n-dimensional topographic map, the ruggedness of which depicts the relationship between morphological form and fitness 4 ( Figure S2). If the topographic relief of the landscape is sufficient, selection should cause traits to evolve rapidly toward peaks and remain there until selection relaxes or the topography of the landscape is altered, for example, due to environmental disturbance or changes in resource availability. The topography of the adaptive landscape can be approximated as an Ornstein-Uhlenbeck process, where the positions of peaks on the landscape are represented as a vector of distinct long-term means, q = (q 1 , q 2 , ., q n ), and the ruggedness of the landscape is represented by the magnitude of a parameter, a, which determines the strength of the pull toward the peaks as an exponential function of time. 15 As a/0, the ruggedness of the adaptive landscape declines and evolution reverts to a random walk. Such a result would not necessarily mean that a trait does not experience strong selection over short (e.g., generational) timescales; rather, it would imply that there are no long-term active trends in trait evolution at the clade level.
As one of the most conspicuous and fundamental features of an organism, body size can improve our understanding of changes in adaptive landscapes over macroevolutionary timescales. 16 Size influences multiple aspects of ecology, physiology, and, ultimately, fitness, [17][18][19] and has therefore become an obvious proxy for life history in analyses of evolutionary tempo and mode. 20,21 We used reversible jump Markov chain Monte Carlo (MCMC) methods 7 to map the adaptive landscape for cetacean body sizes. An advantage of this approach is that it does not require the user to specify the location of peaks a priori; rather, they are estimated from the data. Estimated peaks differed by an order of magnitude, ranging from approximately 2 up to 12 m, although the 95% credibility interval for the estimated peaks indicated that this variation could be even higher, ranging from 60 cm to 50 m (see Figure S3).
It has long been recognized that secondarily aquatic tetrapods tend to be larger than their terrestrial sister lineages. 22 Diverse explanations have been suggested for this pattern. Life in water frees the limbs from mechanical constraints associated with supporting body mass, 23 which should passively increase the upper limit on potential sizes without imposing a minimum size constraint. In addition, heat is more readily lost in water than in air and, because surface area is reduced relative to volume in larger animals, this should place constraints on the lower size limits of fully aquatic taxa. 24 Our results indicate that the adaptive landscape for cetacean body size was profoundly shaped by the thermoregulatory lower limit on size in marine environments, 22,25,26 with a dramatic shift to a new adaptive peak associated with large (posterior mean = 12 m total length) size during the major transitional stage of adapting to fully aquatic life in ambulocetids, from a small (posterior mean = 2 m) ancestral size. (C) Schematic representation of the analytical approach used in the paper. We used data from the literature, new measurements of museum specimens, and the phylogeny of Lloyd and Slater 13 to estimate total length values for species lacking these data. The images in the center of the panel correspond to the dorsal (top) and ventral (bottom) views of the skull of an example odontocete. The imputed dataset and phylogeny were then analyzed using bayou 7 and the Fabric model 8 to detect shifts in the evolutionary regimes of total length, generating insights about the adaptive landscape of this trait. See also Figure S1 for a quality assessment of the imputation approach.

OPEN ACCESS
We found limited support for a trophically defined adaptive landscape for body size in cetaceans, contra previous work on extant taxa alone. 9 The cetacean body size adaptive landscape remains remarkably stable for $20 million years after their initial transition to the marine realm. However, a dramatic decrease in long-term mean body size (from $12.5 to $2.9 m) in the Late Oligocene, along the branch leading to Delphinida (oceanic dolphins, porpoises, and non-platanistid river dolphins), with a second shift to even smaller sizes (long-term mean $ 2.1 m) in the Early Miocene, along the branch leading to total group Inioidea + Lipotidae, is consistent with a dietary model. Maneuverability is negatively correlated with body size in marine mammals, 27 and evolution about these smaller long-term means would provide a performance advantage for delphinidans when foraging for small, agile prey in complex environments. It is tempting to speculate that the decrease in long-term mean body size in the Lipotidae + Inioidea clade is an adaptation to living in freshwater river systems, given their extant distributions. However, this modern distribution is convergently derived, and many fossil taxa in the clade are found in marine sediments. 28 Functional studies of Inia indicate that its flexible body enhances maneuverability at the expense of speed. 27 Thus, members of this lineage may have become smaller to optimally maneuver in shallow coastal environments, which then served as an exaptation for life in rivers. Delphinoids, in contrast, exhibit a suite of vertebral adaptations that result in a stiffer body, 29 which sacrifices some of the maneuverability attained by iniids and pontoporiids for increased speed and swimming efficiency. 27 This likely facilitated the diversification of pelagic piscivores during the Late Miocene, which succeeded the archaic stem lineages and platanistoids 30,31 and culminated in the dramatic, rapid radiation of oceanic dolphins that today comprise almost half of cetacean diversity. 9, 13,32 Aside from these major shifts, we detected few other peaks on the adaptive landscape that affect large clades or that are associated with substantial size change. Shifts within odontocetes represent increases in body size in relation to the delphinidan mean, and were detected for Kentriodontidae and closely related stem delphinidans, and for Orcaella (Figures 3 and S4). Within the baleen whales (Mysticeti), we only observed shifts within the crown group ( Figure 3). These shifts indicate localized decreases in body size within balaenids and in ''basal thalassotherians'' (sensu Bisconti et al. 33 ), with other lineages of balaenopterids containing the remaining changes (Figures 3 and S4). We found no evidence for shifts toward gigantic sizes in crown Balaenopteridae or Balaenidae, the largest extant cetaceans, likely due to the inference that a large optimal body size is acquired early in cetacean evolutionary history. However, we did find a shift toward smaller size in the minke whale. Virtually all the shifts described above were retained across our sensitivity analyses (see STAR Methods and Burin et al. 34  Although the fit of Ornstein-Uhlenbeck models to comparative data is often interpreted in terms of estimates of q, the location of the adaptive peak(s), the parameter a that describes the topography of the landscape and thus the historical strength of selection, is of equal importance. 35 Estimates from wild populations indicate that selection can drive incredibly rapid responses in heritable phenotypes within a single generation, 36 though this may represent short-term fluctuations in response to changes in local conditions 37 that do not accumulate to yield large macroevolutionary divergences. 6 Still, for cetacean body size, we estimate an extremely low rate of adaptation (posterior mean a = 0.0077, 95% highest posterior density = 0.0001-0.0141), indicating that it would take, on average, approximately 90 million years (or 1.6 times the age of the entire clade) to move halfway toward a new long-term mean body size after a peak shift. Such low estimated rates of adaptation indicate that selection on the global adaptive landscape was not the only factor shaping the distribution of body sizes for living and extinct cetaceans.
A flat adaptive landscape for cetacean body size implies that variance should have increased through time, yielding both larger and smaller taxa in spite of an inferred mean size at the upper end of the body size distribution. However, if our models are insufficient, our low parameter estimates may reflect more complex underlying dynamics rather than weak selection. For example, the surprising lack of peak shifts identified in our bayou analyses, particularly in large-bodied extant taxa, could result from that model's inability to accommodate multiple a values, causing us to miss nested shifts associated with periods of accelerated evolution toward a new local peak, or from small sample sizes (such as where shifts occur along terminal branches) that preclude parameter identifiability. To account for these possibilities, we fitted a simpler phenomenological model, the ''Fabric model,'' 8 to our data in which phenotypic shifts along individual branches of phylogeny are treated as biased deviations from an otherwise unbiased random walk in which rates are free to vary. Any trends toward larger or smaller  7 Solid black circles represent increases in average total length (q) for a given regime; hollow circles represent decreases in average total length (q) for a given regime. The bars to the right of the phylogeny show the log 10 -transformed total length values (m) for each species. The colors in the branches and their corresponding bars identify the regime to which each lineage belongs. Silhouettes indicate cetacean families that have at least one living species. In the geological timescale, Q, is the Quaternary. See also Figures S2-S4. ll OPEN ACCESS sizes can be thought of as evolutionary change that falls outside of the range of possibilities permitted by an otherwise unbiased process, though no explicit mechanistic explanation can be invoked.
Using the phenomenological model, we identified 18 branches along which evolution was biased away from an expected change of 0, none of which were found by bayou 7 (Figure 4). 11 of these shifts were toward larger sizes. Notably, some of these involved macroraptorial predators Basilosaurus, Ankylorhiza tiedemani, Livyatan melvillei, and the modern orca Orcinus orca. A massive (175%) increase in size was also detected along the branch leading to crown group balaenopterids (excluding minkes). Additional increases were detected in crown balaenids and the deep-diving ziphiids. Seven trends toward smaller size were also recovered. For example, the ancestor of total group odontocetes underwent a substantial (35%) decrease in size relative to the neocete ancestor, while smaller decreases occurred in the common ancestor of the pygmy sperm whale family Kogiidae, several odontocete terminal lineages, and a pair of nested trends toward smaller size in the toothed mysticete lineage Aetiocetidae.

DISCUSSION
The evolution of morphological diversity demands explanation. Models of gradual evolution predict a simple increase in phenotypic variance through time, with more closely related lineages exhibiting greater similarity to one another, on average, than they do to more distantly related lineages. 38,39 Where selection acts to drive lineages toward distinct phenotypes, it should mold and shape the distribution of phenotypic variation within and across clades, yielding patterns that deviate from the expectations of simple gradualistic models. 5,40 Characterizing these ''ghosts of selection past'' is a central challenge of comparative evolutionary biology. Our analyses suggest fundamental limits in their ability to detect some evolutionary patterns and processes and that a diversity of approaches is required to characterize these complex hierarchical patterns of phenotypic adaptation.
Comparative biologists increasingly model phenotypic macroevolutionary landscapes using the mathematics of the Ornstein-Uhlenbeck process, which connects microevolutionary processes 5 to patterns of variance and covariance over phylogeny. 41 In this framework, peaks on the landscape act as phenotypic attractors while selection strength reflects the steepness of the peaks. However, these models are data-hungry, and phylogenies with massive numbers of tips are required for accurate parameter estimation. 42 This, in turn, can result in a coarse, global view of the adaptive landscape, such as the identification of a single shift of large effect associated with the transition to the marine realm, or a shift to a much smaller mean size associated with increased maneuverability in piscivorous delphinidan odontocetes. What is missed, then, are more localized changes in the topography of the adaptive landscape, which may be particularly difficult to detect when different lineages experience change in different ways. 4 These hidden dynamics are evident as multiple trends in body size evolution that appear on more local scales than the global peak shifts identified by our bayou analyses. In some cases, they occur along branches leading to individual species that can be characterized as possessing similar ecologies, such as trends toward larger size along branches leading to the macropredators (see also Coombs et al. 43 ). In other cases, trends are clade-dependent, such as a recent trend toward larger sizes in balaenopterids and, to a lesser extent, balaenids, which is associated with bulk feeding on dense prey patches, 26,[44][45][46] and ziphiids, which is associated with deep-diving behaviors. 47 Trends toward smaller sizes associated with the origin of odontocetes and kogiids may reflect similar responses to foraging on small, mobile prey relative to ancestral states. 48

Figure 4. Evolutionary trends in total length evolution of living and extinct cetaceans
Phylogeny of cetaceans showing the branches with trends toward increasing or decreasing total length detected by the Fabric model. 8 Magenta lines represent branches with evidence of positive trends in average total length; blue lines represent branches with evidence of negative trends in average total length. Silhouettes indicate the main groups for which significant trends were found. In the geological timescale, Q, Quaternary. See also Figure S2.
The incongruence between a low number of global shifts and multiple conspicuous, more taxonomically restricted local body size trends suggests a much more dynamic adaptive landscape than would be inferred from any single analysis. How could this be the case? Simpson 2,3 conceptualized adaptive zones in terms of phenotypic solutions to functional problems, suggesting that traits associated with ecological behaviors like feeding and locomotion might be associated with rugged adaptive landscapes. 49 Thus the global body size adaptive landscape may be fairly flat when averaged over time, but there could still be distinct peaks on the cranial morphology adaptive landscape, for instance, associated with bulk feeding, piscivory, etc. 43 This means that there are likely multiple phenotypic dimensions to the landscape of which we are only exploring one. However, body size is a complex and multifaceted trait, interconnected with almost every aspect of a species' biology. 17,50 This, in turn, means that body size is under multiple selection pressures, each with its own adaptive mean, resulting in a locally rugged adaptive landscape (i.e., the principle of frustration, sensu Marshall 51 ). We argue that the close proximity of these long-term means facilitates movement between peaks on the adaptive landscape in relatively small increments, which are not identified as discrete shifts in our bayou analyses but appear as trends in the analysis using the Fabric model. 8 Therefore, we describe the adaptive landscape of body size in cetaceans over macroevolutionary timescales by comparing it to a turbulent ocean surface rather than to the more stable topography of terrestrial landscapes (sensu Bell 52 ). The few shifts we identified are analogous to large waves, with the smaller peaks representing unique local conditions of adaptation, 4 more ephemeral in nature, akin to ripples on the surface of the ocean that change rapidly across both space and time. These dynamics likely further increase the difficulty of peaks being detected as discrete shifts in comparative analyses. We note that a major reason our analyses detect these complex dynamics and previous studies do not is that we include fossil taxa. Fossils flatten the adaptive landscape and reveal different peaks. Together, these results highlight the need to include fossils and use multiple comparative approaches to fully characterize adaptive landscapes in deep time, especially in traits as labile and integrated as body size.

STAR+METHODS
Detailed methods are provided in the online version of this paper and include the following:

DECLARATION OF INTERESTS
The authors declare no competing interests.

INCLUSION AND DIVERSITY
One or more of the authors of this paper self-identifies as an underrepresented ethnic minority in their field of research or within their geographical location. We support inclusive, diverse, and equitable conduct of research.

Lead contact
Further information and requests for resources should be directed to and will be fulfilled by the lead contact, Gustavo Burin (gustavo. burin@nhm.ac.uk).

Materials availability
The collated dataset generated in this study has been deposited in the NHM Data Portal. 54 Data and code availability d All original code is available at https://doi.org/10.5519/vmbrpkuq Zenodo: https://doi.org/10.5281/zenodo.7563162 and is publicly available as of the date of publication. DOIs are listed in the key resources table. All R package version details available at Github Repository: https://github.com/gburin/bodysize-evolution-cetacea. d Data are available from the NHM Data Portal 54 d Any additional information required to reanalyze the data reported in this paper is available from the lead contact upon request.

EXPERIMENTAL MODEL AND SUBJECT DETAILS
Total length data We collated total length (m) for living species primarily from stranding records in the literature (n = 1,327 individuals), and the strandings database at the Natural History Museum London (n = 1,659 individuals). If no other data were available we used data from aboriginal subsistence harvests (n = 80 individuals). We only included values within the range of adult lengths for each species, and, where possible, we recorded a minimum of 10 total length values for each species. This resulted in a total length dataset of 2,986 individuals of living species. We collated data on total lengths of fossils from the literature; some were true total length measurements, others were based on estimates (a total of n = 273). 87 species in our phylogeny (see below), two living and 85 fossil, had no available total length data. To estimate total length for these taxa we collated data for nine proxy measurements 55,56 from 913 adult specimens (451 extinct and 462 living) from the published literature, or by measuring specimens using either Mitutoyo Absolute Digimatic CD-8"AX calipers, large calipers, or tape measures depending on the size of the feature. The proxy measurements were (1) width of antorbital notches, (2) bizygomatic width, (3) exoccipital width, (4) occipital condyle breadth, (5) condylobasal length, (6) maximum narial fossa width, (7) maximum nasal width, (8) atlas articular facet width, and (9) humerus length. We obtained a minimum of two proxy measurements from each specimen, as few fossils possessed all features. Prior to the analyses we removed species with uncertain taxonomy and those not present in the phylogeny.

Phylogeny
We used the ''Safe'' analysis phylogeny from Lloyd and Slater 13 and removed all taxa without any total length proxy measurements from the phylogeny leaving 345 taxa. 60 of these remaining taxa were singletons, forming terminals with branch lengths of zero. Because our analyses required branch lengths, we added 0.001 MY to all branches to remove the zero-length branches. To ensure this did not substantially influence our results we repeated all downstream analyses (and sensitivity analyses; see supplemental information) removing the species with zero length branches from the phylogeny (n = 285 taxa; see supplemental information; Burin et al. 34 ). REAGENT

METHOD DETAILS
Total length data imputation To impute total lengths for the 87 species without available total length data, we used a phylogenetically-informed imputation approach. These methods use the eigenvectors of the phylogenetic distance matrix to estimate trait values for species without the focal variable 53,57 while accounting for intraspecific variation. 53,57 We used the phylogenetic imputation algorithm in the R package Rphylopars, 53,57 using a Brownian Motion model for trait evolution, and including specimens with at least one proxy measurement to incorporate intraspecific variation. All measurements were log 10 -transformed prior to the imputation. The final total length dataset contained 345 taxa; 89 living and 256 fossil.
To assess the quality of total length estimates we generated datasets with different percentages of missing total length data for species with available empirical total length data as follows. We randomly removed data for 10% (n = 25), 20% (n = 51) and 50% (n = 129) of all species with existing total length data, leaving 233, 207 and 129 species respectively. Note that we removed total length data for all specimens belonging to these species, so in some cases the percentage of missing total length data records is larger than the percentage of missing species. We also generated a dataset with empirical total length data from extant species only (n = 89). We then assessed the accuracy of the approach by plotting the imputed values against the empirical average values for the species whose total length data was removed from each dataset ( Figure S1). Imputed total length data were strongly correlated with the empirical values, regardless of the amount of information used for the estimation of missing values.
Finally, we calculated the Normalized Root Mean Squared Error for the main dataset (NRMSE 58 ), which represents the mean squared errors for the imputed values standardized by the range of empirical values (Equation 1; see also Table S1).
Where X represents the variable "total length", X Ã i is the imputed value of total length for species i, X i is the empirical value of total length for species i, and n is the number of taxa. Mean squared errors were calculated for estimates using all possible combinations of proxy measurements used for the imputation (ranging from using just one measurement up to all nine). NRMSE values for all four datasets were low (Table S1), and within acceptable values found in other papers. 59-61

QUANTIFICATION AND STATISTICAL ANALYSIS
Exploring tempo and mode in body size evolution We performed disparity through time analyses 62 using the dtt function in the R package geiger, 63 with 10,000 simulated datasets used to generate a null, constant-rates distribution. The centre of gravity was computed as the weighted sum of the product of standardised average subclade disparity and relative node age, with significance determined as the proportion of simulated centres of gravity that were lower than the empirical value. 49 We used bayou 7 to detect evolutionary shifts in total length evolution across cetaceans. bayou uses reversible-jump MCMC algorithms to fit multi-peak OU models of trait evolution, identifying the number, magnitude and location of shifts in q, the evolutionary regime mean of total length. We used a cumulative Poisson distribution as the prior distribution for the number of shifts, setting the prior on the expected number of shifts (l) as 2.5% of the total number of branches in the corresponding tree, i.e. l = 15. We allowed an unlimited number of shifts to occur on each branch, while setting the prior probability for each branch to contain a shift to be proportional to the branch's length. For all other parameters, we used the recommended default distributions found in (https:// github.com/uyedaj/bayou/blob/master/tutorial.md). We ran the MCMC chains for 1 million generations, sampling every 1,000 generations, yielding a posterior sample of 700 after discarding the first 30% as burn-in. We assessed the convergence of each chain using Gelman and Rubin's R statistic, and confirmed that all parameters had an effective sample size greater than 200.
To assess the relative importance of the shifts we followed Uyeda and Harmon's 7 criteria: we selected shifts that had posterior probabilities greater than 0.1 (corresponding to approximately an eightfold increase in the posterior probability compared to the prior probability), with only three of these shifts leading to singleton taxa. For the selected shifts, we calculated the posterior-to-prior ratio as a measure of how much information our data is providing on the location and magnitude of changes in evolutionary regime. All bayou analyses were performed in R. 64

Sensitivity analyses
To test the sensitivity of our results to both tree topology and the prior parameter for expected number of shifts (l), we replicated the bayou protocol described above for the original tree (Full tree) and six additional trees, with five different values of l.
The six trees were as follows: 1. Full. The most comprehensive phylogeny possible that included all species with total length data (n = 345). This is the the ''Safe'' analysis phylogeny from Lloyd and Slater 13,14 with species without total length data (empirical or imputed) removed, and is the tree used for our main analyses. 2. No Archaeoceti. The Full tree after dropping all species that are classed as archaeocetes (n = 322).