Humerus shape evolved in cetaceans under relaxed selection and random drift

Abstract

Mammalian humerus shape is mostly predicted by loading. This constraint is lifted for cetaceans and, therefore, can reflect other, specific evolutionary trends. To test this hypothesis, the three-dimensional shape of the humerus for 36 taxa of extinct and extant cetaceans was analyzed. The shape variance concurs with the evolution of diverging baleen and toothed whales and is independent of allometry. It shows traits associated with aquatic locomotion: humeral head twist, diaphysis shortening and straightening, and epiphyses enlargement. Also, changes in the anatomy of the humeral head and greater tubercle, modularity and integration patterns are associated with phylogeny. We suggest relaxed negative selection and random drift as the likely main evolutionary patterns for the evolution of the humerus shape. The evolution of genes regulating the humerus shape corresponds to the anatomical transformations and shows relaxed selection, site-level positive selection and nonsense mutations in cetaceans. Overall, the observed pattern well illustrates a “fly in a tube model” and shows its importance as a framework of emerging evolutionary innovations.

Introduction

Humerus is an essential element of the forelimb in tetrapods^1,2,3, a key element enabling their locomotion and involved in other multiple functions^4,5. The shape of the humerus in non-aquatic organisms is strongly predicted and constrained by functional factors associated with terrestrial or arboreal locomotion, such as body support, walking, digging, or climbing^{6,7,8,9,10,11,12}. Its functional constraints are also well seen in flying organisms^13,14,15. On the contrary, fully aquatic vertebrates, which are not constrained by terrestrial or aerial locomotion, show the greatest disparity in humeral shape among other ecological groups^16,17. Expectedly, this constraint can be equally lifted for cetaceans (whales and dolphins) which transited to fully aquatic life about 46–42 million years ago (Ma)^18,19 and, therefore, it could reflect other specific evolutionary trends. Flippers of living cetaceans, the forelimbs stiffly encased in soft tissues²⁰, act as hydrofoils enabling stability and maneuvering during the swimming²¹ and are used in social behavior^22,23. The humerus plays some role in the limited mobility but novel functionality of the flipper, since the shoulder joint is the only joint retaining mobility in a rigid cetacean limb and the humerus bears insertions for at least six muscles²⁴. In fact, the humeri of modern whales and dolphins are diverse in shape, as well as in size. The size limits vary between 5 and 70 cm in length, and the shape may differ in proportions of the humeral head (which may be longer than a half on the whole bone), orientation (proximodistal rotation) of the humeral head in relation to the shaft (diaphysis) and extension of crests and tubercles²⁵. However, the shape component of humerus variation is hardly quantified because of the shortening, thickening, and smoothening of this bone with diminishing crests. Evolution of cetacean humerus needs quantitative description and explanation in terms of general variation, evolutionary rates of different parts and structures, phylogenetic signaling (if present), and, specifically, integration and modularity. The definition we use here is that the integration is the ability of different parts to vary jointly, whereas the modularity occurs when integration concentrates within specific parts of the structure so that within-module integration exceeds between-modules integration²⁶.

Lifting loading constraints could lead to changing the tempo of morphological evolution or the disintegration of shape organization, as it was also evidenced in some vertebrates²⁷. On the other hand, selective pressure in a dense tridimensional water environment¹ would lead to a highly specific transformation of the limb and might result in its reintegration and the rise of new forms and states of modularity²⁸. Whether evolutionary changes in cetacean humerus were due to random genetic drift (neutral evolution^29,30,31) or selection that could be in turn negative (i.e. stabilizing or purifying) or positive (leading to adaptive evolution^{32,33,34,35,36}) remained an open question. Therefore, we aimed to understand if cetacean humeral shape evolution is more driven by phylogenetic or functional signals. We hypothesized that evolution of the humerus could follow the evolutionary pattern of Brownian motion, and in this case the shape, integration and modular structure of the humerus would change non-directionally. Therefore, correlations between and within modules of the humerus might randomly decrease or increase, but in general they were not expected to change much from the ancestral state. Meanwhile, functional constraints would result in specific transformations and should be well distinguishable after correcting for phylogenetic relationships between species. For example, different shapes of the humerus might have evolved due to adaptations to slow swimming (e.g., in most baleen whales) or fast swimming (e.g., in many toothed whales). However, the exceptions in locomotion are known from both main lineages of Cetacea: slow swimmers among toothed whales (i.e. Inia, Delphinapterus, Kogia) and fast swimmers among baleen whales (e.g., Balaenoptera).

For testing these hypotheses we analyzed patterns of evolution of the shape of the humerus and regulatory genes which can be likely associated with the shape changes. We found interplay of factors under these processes and identified relaxed negative selection and random drift as the likely main evolutionary patterns for the evolution of the humerus shape which revealed themselves in phenotypic variation, modularity an integration patterns, and evolutionary rate. However, at the same time, evolution toward increasing aquatic adaptations was suggested. Also, the evolution of regulatory genes was congruent to the anatomical transformations. Interestingly, this combination of different evolutionary processes can be well illustrated by a “fly in a tube model”.

Results

We analyzed the three-dimensional shape of the humerus (Fig. 1; Supplementary Fig. S1; Supplementary Fig. S2; Supplementary Data 1) for 36 taxa of extinct and extant cetaceans which were pooled, for the aims of this study, into several major groups based on combination of their chronostratigraphy, lifestyle and phylogeny: Protocetidae, stem cetaceans which lived between 50 and 37 million years ago (Ma) and transited from semiaquatic to fully aquatic lifestyle³⁷; early members of Pelagiceti, fully aquatic cetaceans belonging to several lineages including Basilosauridae and members of both crown groups, Mysticeti and Odontoceti, which lived before the end of the Paleogene (46-23 Ma)^38,39,40,41; and two groups morphologically similar to living whales and dolphins – Neogene and extant Mysticeti, baleen whales; and Neogene and extant Odontoceti, toothed whales including dolphins (Supplementary Table S1; Fig. 2).

Fig. 1: Landmarks for geometric morphometry of a right cetacean humerus.

Landmark partition into modules: blue, proximal epiphysis; red, diaphysis; yellow, distal epiphysis. A lateral; (B), anterior; (C), medial; (D), posterior; (E), proximal; (F), distal view. dc deltopectoral crest, gt greater tubercle, hh humeral head, lt lesser tubercle. The list of the landmarks with explanation is provided in Supplementary Table S3.

Fig. 2: Phylogenetic tree and groupings of taxa for the morphological study.

Cetacean phylogeny is based on Lloyd & Slater⁹⁰ metatrees and specific taxonomic studies³⁸. Species are grouped according to their chronostratigraphy, lifestyle and phylogeny: Protocetidae (2 taxa); early members of Pelagiceti including Basilosauridae and members of Mysticeti and Odontoceti (8 taxa); Neogene and extant Mysticeti (7 taxa); and Neogene and extant Odontoceti (19 taxa).

Also, we searched for signs of selection in evolution of regulatory genes which can accompany the morphological evolution of the cetacean humerus and are promising candidates for explaining evolutionary innovations⁴². Numerous genes are involved in the formation of limbs: different expression, timing, interactions with other genes lead to the variety of limbs proportions and shapes seen in vertebrates^43,44. Alteration of the signaling pathways of the limb morphogenesis may lead to malformations. For example, early in the limb development, altered HOXA9, HOXA10 and HOXD9, HOXD10 genes produce a shortened humerus with reduced deltoid tuberosity in mice⁴⁵. At later stages, endochondral ossification of the humerus is regulated by PANX3; knocking it out in mice produces a short diaphysis and a larger deltoid tuberosity⁴⁶. Altered GLI3 and ALX4 lead to alteration in all the forelimb, including the humerus⁴⁷, and the same may be true for PTCH1 and PTCH2⁴⁸. Proximodistal patterning of the limb is also regulated with PBX family genes⁴⁹. Therefore, we chose eleven genes (HOXA9, HOXD9, HOXD10, PANX3, GLI3, ALX4, PTCH1, PTCH2, PBX1, PBX2, PBX3) likely regulating the humerus shape, and extracted them from the genome assemblies of cetaceans and outgroups (Supplementary Table S2) to check how the evolution of humerus in cetaceans might be reflected in the evolution of these genes.

Phenotypic variation and phylogenetic signal

The set of 34 landmarks (Supplementary Table S3; Supplementary Data 1) provided a good approximation for the whole shape of the humerus. The humeri examined by us represented the whole size range of cetaceans: the smallest centroid size was found in the harbor porpoise (Phocoena phocoena), and the greatest, almost 15-fold greater, was in the blue whale (Balaenoptera musculus). The shape of the cetacean humerus was poorly explained by the size (1.8% of the variation of phylogenetically corrected Procrustes aligned landmarks, F = 0.611, p = 0.814).

We assessed whether the shape of the humerus could reflect functional constraints. For this aim, we compared morphospaces built with principal components deriving from three different principal component analyses (PCA), which concentrated on specific aspects of traits’ variation. The standard PCA described total variation, the phyPCA excluded phylogenetic variation (e.g., showed ecological niches, which were independent from phylogeny), and the PaCA focused on phylogenetic variation (Fig. 3). Blomberg’s K test for phylogenetic independence of the trait values confirmed the successful exclusion of phylogenetic relatedness in the phyPCA and high phylogenetic signal (K > 1, which meant that trait values were well predicted by the phylogeny) in the standard and phylogenetically aligned analyses (K and significance values: Supplementary Table S4). The main drivers of the shape differentiation were assessed from comparisons of the morphospaces: if they were phylogenetic, the morphospace of the standard PCA would be similarly structured with that of the PaCA, and if non-phylogenetic, with the phyPCA.

Fig. 3: Multivariate analysis of the shape of the humerus.

Sources of variation can be assumed comparing morphospaces deriving from different Principal Component Analyses. Major shape changes associated with the principal components analyses (the humeri shown from the posterior and medial views) are illustrated on each axis; their 3D models are shown at the Supplementary Fig. S2. A Standard PCA with projected phylogenetic relationships (phylomorphospace). B Corrected for phylogeny phyPCA: phylogenetic signal is excluded. С Phylogenetically aligned component analysis (PaCA): variation that is the most related to the phylogenetic signal. Taxa names are coded as following: Aege – Aegicetus gehennae, Bamy – Balaena mysticetus, Baac – Balaenoptera acutorostrata, Baed – Balaenoptera edeni, Bamu – Balaenoptera musculus, Ceri – Cetotherium riabinini, Cype – Cynthiacetus peruvianus, Dele – Delphinapterus leucas, Doat – Dorudon atrox, Esro – Eschrichtius robustus, Feke – Ferecetotherium kelloggi, Glme – Globicephala melas, Grgr – Grampus griseus, Hypl – Hyperoodon planifrons, Keba – Kelloggia barbarus, Kefu – Kentriodon fuchsii, Kein – Kentriodontidae indet, Kosi – Kogia sima, Laac – Lagenorhynchus acutus, Laal – Lagenorhynchus albirostris, Main – Maiacetus inuus, MH – MHN-UABCS-Ec/Sj29/42b, Mica – Microzeuglodon caucasicum, Momo – Monodon monoceros, Olaz – Oligodelphis azerbajdzanicus, Orbr – Orcaella brevirostris, Oror – Orcinus orca, Otvi – Otradnocetus virodovi, Pasu – Pachyacanthus suessii, Phph – Phocoena phocoena, Phsi – Phocoena sinus, Phma – Physeter macrocephalus, Stco – Stenella coeruleoalba, Suda – Sulakocetus dagestanicus, Tutr – Tursiops truncatus, Zica – Ziphius cavirostris.

The morphospace of two major standard principal components was well structured and concurred to phylogeny (Fig. 3): heterogenous Paleogene cetaceans (both Protocetidae and early Pelagiceti), Neogene and extant baleen whales (Mysticeti), and Neogene and extant toothed whales (Odontoceti) occupied three different areas of the plot. The similar position of taxa was seen in PaCA. Therefore, the phylogenetic scenario of shape differentiation was supported. Also, the phylogeny might influence evolution of the ecological niche: closely related whales had similar shapes of the humerus and also occupied similar niches (and, therefore, had similar swimming styles). Thus, exclusion of the phylogenetic signal also affected the phylogenetically dependent response to ecological signal, and no structure was seen in the phyPCA morphospace.

The first principal components of the standard PCA and PaCA showed the gradual advance of general aquatic adaptation and transition to elaborated aquatic locomotion. The second principal components described the divergence between baleen and toothed whales. Also, each analysis showed some specific aspects of shape variation (Figs. 3, S2). Major shape changes described with the first principal component were similar in all the analyses (Supplementary Table S5). Its variation was associated with changes in the humeral shaft, relative size and position of the greater tubercle, humeral head and epicondyles (Supplementary Table S6). The positive end of the standard PC1 axis (Stenella coeruleoalba, Orcinus orca, Lagenorhynchus acutus and L. albirostris) was associated with a short shaft, broad epiphyses and a prominent lesser tubercle. The negative end of the standard PC1 axis (Maiacetus inuus) was associated with a slender and long shaft and well separated greater and lesser tubercles. Variation of the PC2 and PaCA2 was mostly associated with the shape of the proximal epiphysis and, in particular, the shape and proportion of humeral tubercles, and the shape of the distal epiphysis (Supplementary Table S6). Meanwhile, the main cluster of the phyPCA morphospace was composed by diverse, heterogeneous shapes, and among the outlying taxa there were so different forms, as e.g. Tursiops truncatus, Orcaella brevirostris, Monodon monoceros and Balaenoptera musculus (Fig. 3B). The phyPC2 mostly showed variation in the size of the greater tubercle (Supplementary Table S6). Also, the PaCA analysis showed that the proximal epiphysis landmarks contributed the most to the phylogenetic variation of the studied Cetacea (PaCA1 = 46%, PaCA2 = 64%, Supplementary Table S6).

Twist of the humeral head is particularly noteworthy due to the major reshaping of the bone during evolution of cetaceans, being among their remarkable apomorphic features²⁵. Our principal component analyses reflected its variation. Proximodistal rotation of the humeral head was calculated by comparing configurations for each species to the most ancient species in our sample, Maiacetus inuus (Fig. 4). The twist was the most prominent and variable in odontocetes, from 10 to 38 degrees. However, it was also found in most mysticetes (from 2 to 13 degrees), in a basilosaurid Cynthiacetus peruvianus (16 degrees) and even in a protocetid Aegicetus gehennae (27 degrees). Both standard PCA (PC1 r_s = 0.469, p = 0.004, PC2 r_s = 0.676, p < 0.001) and PaCA (PaCA1 r_s = –0.488, p = 0.003) were associated with the degree of the rotation.

Fig. 4: Twist of the humeral head (hh) in cetacean phylogeny.

The configuration of the landmarks of Maiacetus inuus is taken as the reference. The colorful bar length equals 20 million years and represents the rotation degree.

Scores of the first three principal components from the standard PCA were used for the dimensionality reduction. PC1 and PC2 were used as complex traits reflecting the main trends in the variation of the humerus. We used them to infer potential signals of adaptive evolution in the shape of the humerus. In addition to a priori grouping (presented at Figs. 2 and 3), we let the data speak for itself. The Bayesian framework was applied to formulate the hypothesis about probable adaptive shifts for the first two standard principal components (Supplementary Fig. S3). This is a univariate method, so each PC was tested separately. The main adaptive shifts in PC1 were found at divergence between Maiacetus and other whales (posterior probability, pp = 56%) and at the base of the delphinoid clade (Kentriodontidae, Delphinidae, Monodontidae, Phocoenidae, pp = 38%). The analysis of the standard PC2 revealed main adaptive shifts at the base of Mysticeti (pp = 74%) and Physeteroidea (pp = 30%) clades, which converge in the shape of humerus (Fig. 3A: see positions of Mysticeti: Eschrichtius robustus, Otradnocetus virodovi and Physeteroidea: Kogia sima, Physeter macrocephalus).

The hypothesized adaptive pattern of evolution in the humerus shape was estimated by testing different evolutionary models for the several standard principal components. We checked our a priori groups (Fig. 2) and a posteriori (i.e., after the Bayesian inference) probable adaptive shifts to two optima: Delphinoidea and stem Mysticeti + Physeteroidea. We used a multivariate approach (mvMORPH) to capture most of the variation in the humerus and test for the best-fitting evolutionary process. The most likely evolutionary model (Supplementary Table S7) was an Ornstein-Uhlenbeck process with a specific adaptive peak in Delphinoidea, and another peak in Mysticeti + Physeteroidea, vs. all the other taxa. Humerus shape adaptation to these peaks evolved relatively fast, multivariate half-life t_0.5 = 6.76% of the tree height (3.44 Ma). Strength of selection and rate of stochastic evolution were considered as constant values in the models. A more complex model with varying evolutionary parameters could not be tested on our data due to the small sample size.

Modularity

We expected that different adaptive optima may affect the modularity structure of the humerus. Landmarks were divided into modules according to the anatomical parts of the humerus: distal epiphysis, diaphysis and proximal epiphysis. The phylogenetic signal of the 34 GPA aligned landmarks was significant (multivariate Blomberg’s K = 2.016, p < 0.001), and therefore the phylogenetically corrected landmarks were used in the modularity analysis (phylogenetic relatedness was successfully excluded: multivariate Blomberg’s K = 0.175, p = 0.991).

The inner correlation structure of the humerus was less than 0.5, which we believe to indicate a low contrast system of modules (Fig. 5, Supplementary Table S8). The most likely model of module partition for the overall cetacean sample for phylogenetically corrected landmarks was one with three modules with several within- and between- module correlations (Fig. 5, Supplementary Table S7). It had the highest support among eleven models (Supplementary Table S7). Also, we analyzed the correlation structure for each of three major a priori subgroups (Neogene and extant Odontoceti; Neogene and extant Mysticeti; Paleogene Pelagiceti) and two a posteriori subgroups (Delphinoidea; All Mysticeti (including the stem and crown members) + Physeteroidea) in our sample.

Fig. 5: Correlations (ρ) and disparities of the modules for each group.

Graphs with three nodes represent maximum likelihood estimates of the correlations: values in shapes correspond to within-module correlation; values near lines are between-modules correlation; size of shapes and lines represent magnitude of correlations. Histogram bars represent distributions of disparities of the 10,000 random partitions of the landmarks; estimates of the observed disparities are presented as shapes along the x-axis of the histograms. Vertical dotted lines at histograms of the random distributions is the 0.95 quantile. High disparity is considered when the observed sample estimate lays to the right of the 0.95 quantile of the random distribution. A the total sample. B Neogene and extant Mysticeti. C Neogene and extant Odontoceti. D Paleogene Pelagiceti. E All Mysticeti and Physeteroidea. F Delphinoidea. d diaphysis, de distal epiphysis, pe proximal epiphysis.

Paleogene Pelagiceti could be best described with a two-module structure, with the proximal epiphysis as a separate module (sum of Akaike weight for two models where proximal epiphysis is considered as separate module = 0.887, Supplementary Table S8). The humeri of the Neogene and extant Odontoceti and Delphinoidea had well defined three modules (sum of weights of four three-modules models >0.9, Supplementary Table S8). For the Neogene and extant Mysticeti all the models had low probabilities (Akaike weight < 0.25). After summing weights of the models with the same landmark partitioning (these models differ in how they account for within- and between-module correlation), Neogene and extant Mysticeti had the highest support for the two-module partition with the proximal epiphysis as a separate module (sum of weights = 0.376), while the no module hypotheses was supported with the weight = 0.130. Therefore, no firm conclusion can be drawn here. However, low probabilities of the models and low contrast between modules and their connections is an indirect evidence for the higher integration of the whole humerus in Neogene and extant Mysticeti (Supplementary Table S8). The modularity for the subgroup “All Mysticeti + Physeteroidea” (which corresponds to the adaptive peak suggested by Bayesian inference) was also low, and it was best described as a two-module structure with the separated proximal epiphysis (sum of weights of two models with proximal epiphysis as a separate module = 0.808, Supplementary Table S8).

Between and within module correlations, which were evaluated by maximum likelihood approach, were used as measures of integration. The proximal epiphysis had the lowest within-module correlation, and it was also the most separated from the other parts of the humerus (Supplementary Table S8, Fig. 5).

We used the morphological disparity as the measure of the variation of shapes and the phylogenetic correction for finding modules which diversified the most independently from phylogeny. For the overall pooled sample, the diaphysis had a significantly higher morphological disparity (compared to the distribution of disparities of random combinations of landmarks), and both epiphyses showed a lower disparity (Fig. 5, Supplementary Table S9). Within each subgroup, non-random diversification was found in different modules. Neogene and extant Mysticeti had more variable epiphyses. However, a wider group of “All Mysticeti + Physeteroidea” (which corresponds to the adaptive peak suggested by Bayesian inference) had larger disparity in the proximal epiphysis and the diaphysis (Fig. 5B, E). Disparity of all three modules was higher than random in Neogene and extant Odontoceti. However, variation in the narrower group of Delphinoidea was mostly seen in the proximal epiphysis and the diaphysis (Fig. 5C, F). In Paleogene Pelagiceti, all three modules had a high disparity corrected for phylogeny (Fig. 5D).

Low within-correlations and low contrast between modules may increase variation. The proximal epiphysis is a “weak” (low within-module correlations) module that is less dependent on the other modules (low between-module correlation). Relaxed selective constraints might increase evolutionary rates and morphological disparity in some lineages.

Rate of evolution

Biological evolution described by the Ornstein–Uhlenbeck equation is a sum of two processes: selection toward an optimum and stochastic drift regardless of the location of the optimum. The former process reflects the strength of selection, and the latter shows the neutral evolutionary rate. If the evolution is neutral, the Ornstein–Uhlenbeck process simplifies toward the Brownian motion³⁶.

Earlier we have assessed different models of evolution by estimating the most likely parameters of the Ornstein–Uhlenbeck equation, morphological integration and disparity of different parts of the humerus. Here, we hypothesize about the situation where each landmark evolved independently and the selection was absent. Under these conditions, the variation could be random and proportional to the evolutionary rate of Brownian motion, which would represent neutral evolution. Otherwise, the variation could be limited by stronger negative (stabilizing) selection or by the correlative restraining influence of other traits. Also, the variation could increase due to stronger positive selection or the correlative influence from other traits. Thus, the rate of evolution plotted against the variance of morphological traits, which are represented here by positions of landmarks, shows the evolutionary scenarios for shape changes^50,51.

The highest variance and evolutionary rates were observed in several landmarks of the proximal epiphysis (Fig. 6, Supplementary Fig. S4; L1, L2, and L3). The net evolutionary rate for the landmarks of the diaphysis was significantly slower (σ² = 1.52e–05, p < 0.038) than for the other modules. The rates for the proximal epiphysis (σ² = 2.74e-05) and the distal epiphysis (σ² = 1.87e–05) were higher and statistically indifferent from each other (p > 0.05).

Fig. 6: Regression of Procrustes variance on within-landmark evolutionary rate.

The black dashed line is the expected relationship of the evolutionary rate and the mean of the simulated variances under the Brownian motion mode of evolution. The black solid line is the observed relationship for all landmarks. According to the theoretical framework by Felice et al.⁵⁰, constraints or convergences may be suggested in the evolution of the shape when observed regression lines have lower slope than the expected line. Higher slopes may indicate strong disruptive selection or facilitated morphological disparity. Colored lines and signs represent modules. d diaphysis, de distal epiphysis, pe proximal epiphysis.

All the linear models (for the whole humerus and each module) of the observed within-landmark evolutionary rates on Procrustes variances were significant (Supplementary Table S10; Supplementary Table S11). Evolution of the shape of the whole humerus significantly differed from the expected scenario of the neutral evolution: lower than expected variances were observed, showing the negative selection (p = 0.010, Fig. 6). Meanwhile, at the level of separate modules the marginal significance was found for proximal epiphysis (p = 0.095) and diaphysis (p = 0.054), also suggesting the weak negative selection, whereas in the distal epiphysis the expected and observed regression lines overlapped with a high support (p = 0.769, Supplementary Table S10) supporting the neutral evolution of that structure. The expected neutral evolution rates were estimated as a mean of 10,000 simulations. This might give too optimistic statistical inference. So, we applied graphical analysis of the regressions of all simulations.

Graphical analysis of regressions with the global rank envelope test of the simulated data showed the high variance of the estimates of the regression slopes (Supplementary Fig. S4). The observed lines of regression were within the envelope with a high support (p > 0.244, Supplementary Table S9). Only several landmarks fell outside of the boundary interval (Supplementary Fig. S3, Supplementary Table S12). Constrained evolution may be suggested for the relative length of the deltopectoral crest (L9, L11, L15) and the lateral facet of the greater tubercle (L3, L33). Contrary, facilitated evolution may be suggested for the entepicondyle of the humerus (L29, L30).

Evolution of regulatory genes: test for positive selection

In general, the main trend for selection in HOXA9, HOXD9, HOXD10, PANX3, ALX4, GLI3, PTCH1, PTCH2, PBX1, PBX2 and PBX3 was the strong purifying selection (Fig. 7). However, specific details of this trend varied across the target phylogenetic lineages for each gene. Using branch models on separate branches, we identified several lineages for which nonsynonymous to synonymous rate ratios (ω values) significantly differed from the background (Fig. 7; Supplementary Table S13). Among them, results for PTCH1 in Hyperoodon ampullatus (0.566) and Balaena mysticetus (0.886) were much higher than the background value of 0.092. However, they did not reach cutoff levels of ω = 1, characteristic for neutral selection, meaning that selection acting on these branches was still purifying (i.e., negative), but relaxed. Instead, in the Cetacea ancestral branch, selection on the PBX3 gene was intensified, as ω = 0.0001 was lower than in the background (0.377).

Fig. 7: The results of positive selection analyses, calculated in PAML.

A Tree topology used in the study. Test lineages are highlighted in red and significant results (p < 0.05) of branch models are shown. The arrows near the gene’s names indicate increased or decreased ω values. B Comparison of ω values by gene between odontocetes and mysticetes. Significantly differing values (p < 0.05) are indicated by an asterisk. Icons are from PhyloPic (www.phylopic.org): Hippopotamus amphibius by Steven Traver; Balaena mysticetus, Hyperoodon ampullatus by Xavier Giroux-Bougard; Balaenoptera physalus, Monodon monoceros by Guillaume Dera; Physeter macrocephalus by Margot Michaud; Neophocaena asiaeorientalis by Mozillian (all in Public Domain); Delphinus capensis, Sousa chinensis by Chris huh, and from Wikipedia (commons.wikimedia.org), Orcinus orca by Chris huh (all under a Creative Commons license).

The comparison of ω values between Mysticeti and Odontoceti clades showed significantly increased results in PBX3 for baleen whales (0.877 vs 0.27, Fig. 7B). Additionally, Odontoceti had higher values than Mysticeti for the GLI3 gene (Fig. 7B). Using the branch-site model in PAML, two amino acid sites under positive selection were identified for GLI3 gene in Kogia breviceps: 260 (T → G) and 1184 (S → N). All significant p-values are shown in the Supplementary Table S13.

Discussion

Our results on the assessment of the tempo of evolution are congruent for humerus morphology and genes regulating its morphogenesis. The shape of the humerus reflects shift in the forelimb muscular efficiency during the course of the adaptation to aquatic life. Also, stronger negative (purifying) selection of the limb-related genes was seen in the cetacean ancestral lineage. Therefore, relaxed purifying selection on limb-related genes was likely to increase variation of the shapes of humerus after the transition to aquatic lifestyle which is clearly seen in the Paleogene Pelagiceti. From the deep time perspective, the swimming performance evolved toward at least two adaptive optima represented by baleen and sperm whales (group “All Mysticeti and Physeteroidea”) and dolphins (group “Delphinoidea”). However, in general, landmark variances and evolutionary rates increased proportionally to each other (suggesting a neutral pattern of evolution) or variances of some landmarks were slightly lower than evolutionary rates (suggesting weak negative selection). We will discuss this in detail.

Major trends in the evolution of cetacean humerus: Phylogenetic signal vs. environmental signal

The latest common ancestor of baleen and toothed whales was fully aquatic, with flipper-shaped forelimbs, but still with some plesiomorphic traits in the humerus: well pronounced humeral tubercles, a long deltopectoral crest, and a trochlea forming the elbow synovial joint^38,52,53. Differentiation of forelimb anatomy in the lineages of Mysticeti and Odontoceti evolved well after the transition of their common ancestor to water^39,41,54. Subsequently, their modularity networks between different parts of the humerus changed divergently: integration of the whole bone increased in Mysticeti (or at least did not change from the ancestral modular structure), whereas the modularity of the specific parts of the bone increased in Odontoceti, and particularly in Delphinoidea. The humeral shape in these two lineages largely evolved in the proximal part (Fig. 6, Supplementary Fig. S4). Low within-module correlations of the proximal epiphysis in Odontoceti might yield lower amounts of correlative changes resulting in the low (relatively to the other modules) morphological disparity, and thus the average shape of the proximal epiphysis remained somewhat similar to the ancestral shape with well pronounced tubercles. High disparity of the proximal epiphysis in Mysticeti might be due to the reduction of the greater tubercle. Possibly, relaxation of constraints (low within-module correlations, high disparity) in proximal epiphysis and diaphysis boosted their variation in dolphins, leading to increase of the range of the available phenotypes and evolution of highly derived states of twist of the humeral head.

Thus, evolutionary changes in different clades varied in different directions. From the first glance, this dynamics resembled the Brownian motion pattern of evolution. Meanwhile, it also might result from the numerous shifts of the adaptive optima. Evolutionary change of different parts of the humerus occurred at different scales and directions that in general followed weak stabilizing selection (Rate of evolution, Fig. 6). The best model to describe evolutionary dynamics of the humerus shape was the model with two optima shifts to Delphinoidea and to Mysticeti. The physeteroid (basal toothed whale) shape was close to Mysticeti, and this might be evidence for retaining plesiomorphic shape or convergent evolution. Nevertheless, the main driver of the humerus shape changes was phylogeny (as seen from the results of all the principal component analyses).

However, the environmental signal seems to be equally strong. The analyses of the shape variation (both PCA and PaCA) show that major transitions in humeral shape are shared by mysticetes and odontocetes, and these aquatic adaptations increased through the deep time. These traits include shortening and straightening of the diaphysis, widening of both epiphyses, development of the lesser tubercle, reduction of the deltopectoral crest and, at least to some extent, twist of the humeral head along with its enlargement – the features mostly recorded in earlier anatomical descriptions^20,25,55,56. Our quantitative analysis also uncovered a common trend in increasing integration of the diaphysis. Also, the higher disparity of distal epiphysis in odontocetes could be explained, among other factors, based on the higher variability of their locomotion patterns which is influenced by highly diverse behavioral and feeding strategies. They, for example, include fast swimmers like Stenella and Delphinus, and slow but maneuverable Kogia and Platanista^20,21,57. Also, the highest disparity values among delphinoids could be partially due to their better sampling for our study. From an anatomical viewpoint, the differences in the distal epiphysis can be attributed to differences in muscle apparatus⁵⁶. All living mysticetes have well-developed flexor and extensor muscles and can actively change their flipper curvature. On the contrary, living toothed whales show diverse morphology: e.g., sperm whales and beaked whales have flexor and extensor muscles, whereas dolphins and porpoises lack them⁵⁶.

“A fly in a tube” model

Evolution of the cetacean humerus showed a strong phylogenetic signal in shape variation. However, humerus shape also evolved towards increasing aquatic adaptations, with similar changes in several lineages of crown baleen and toothed whales (two major shifts of the adaptive optima were suggested). Variance-on-rate regression analysis was performed on averaged simulations that could be too optimistic, so we applied graphical analysis of all simulations with the envelope test. As we found, strong negative selection for the whole humerus and its parts could not be confirmed. Less variation was accumulated than expected under Brownian motion (Fig. 6), whereas a global rank envelope test turned out to be non-significant suggesting neutral evolution.

High phylogenetic signal may result from the Brownian motion evolutionary pattern or from other various types of adaptive processes. Among them, there are weak constant stabilizing selection, variable stabilizing selection with high rate of stochastically fluctuating adaptive peak shifts and time-dependent evolution where initially high evolutionary rates and rates of adaptive peak shifts decrease over time⁵⁸. In all these processes, stochasticity is important in shaping the evolutionary dynamics.

“Fly in a tube” model⁵⁰ well explains the combination of different evolutionary processes: populations (“a fly”) evolve along their lines of least resistance (i.e., easily occupy any points of the morphospace bounded with the “tube” walls); if bounds were relatively wide and shifts to adaptive optima were quick, the overall variation might not be affected and remained proportional to the neutral evolutionary rate. Thus, diverging evolutionary events could have happened within the Brownian motion pattern of evolution.

Evolutionary dynamics of the humerus shape seems to be complex. We found two arguments in support of the neutral evolution: variation in the landmarks was too high for stabilizing selection, and the humerus morphospace was phylogenetically structured. Thus, we suggested relaxed negative selection or random drift with weak stabilizing selection as a likely main evolutionary pattern.

Deviations from the neutral evolution may be due to an interplay of phenotypic diversification and changes in within-structure integration⁵⁹. Facilitated variation or disruptive selection and constrained variation or convergences may be suggested here⁵⁰. Originally described for the evolution of the avian skull, this model shows a process simultaneously driven by Brownian motion evolutionary pattern and developmental and ecological constraints⁵⁰.

Our research, along with other recent studies^60,61, shows the importance of the “fly in a tube model” as a framework for emerging evolutionary innovations when new phenotypes may evolve as highly integrated structures despite their divergence under a Brownian motion model. Moreover, as shown by our analysis, the modular system of the humerus also evolved divergently: from the weakening of the modules in baleen whales to a deeper separation of modules in toothed whales. Meanwhile, the variation of multiple traits of the humerus shape along the PC1 shows parallelism in forming new “whale-like” phenotypes across two main lineages of modern cetaceans – baleen and toothed whales⁵³.

This puts us back to a discussion about the origins of evolutionary innovation and its relatedness to anatomical (dis)integration^{60,62,63,64,65}. It seems that multiple evolutionary innovations may evolve as parallel events under the drift coupled with the relaxed selection in some circumstances (e.g., secondarily aquatic adaptation) and also involve some increasing morphological integration. This is an interesting effect which deserves further research.

Interestingly, the evolutionary patterns observed in the humerus of cetaceans are somewhat different from other ecological groups of mammals^11,12,66. Whereas in many groups the environmental factors were found as main drivers of morphological variation of the humerus with multiple convergence events due to similar lifestyle¹² and wide divergence of phylogenetic lineages in the morphospace¹¹, the evolution of the cetacean humerus seem to be largely driven by phylogenetic signals, and the shared trends in the evolution of humerus morphology seen through deep time (e.g., seen under the standard PCA) may be possibly best described as parallel rather than convergent or divergent, since they are highly mosaic in their manifestations. Meanwhile, this only shows the specific evolutionary mechanisms for secondarily aquatic adaptations. Moreover, our findings are congruent with those reported from the evolution of bat humeri, in which the significant phylogenetic signal was found, especially in the diaphysis, whereas the epiphyses anatomy was mostly driven by environmental signals¹³. Also, the competing trends and effects in different modules or traits of the humerus is the feature uniting cetaceans with bats¹³ and, possibly, arboreal mammals¹², and this may be the evidence for the mosaic evolution as the widespread phenomenon for the evolution of mammalian humeri¹². Broadening the range of taxa to be examined, especially on the basalmost extinct forms and other deep lineages of aquatic mammals, seems beneficial for understanding the bigger picture^67,68. Also, the “fly in a tube” model may be beneficial for the description of evolution of structural elements in other groups in which the strong phylogenetic and environmental signals were simultaneously found, such as the bat humerus¹³.

Genetic aspects of evolution of the cetacean humerus

Relaxed selection in limb related regulatory genes can be associated with the lifting of constraints and contribute to the evolution of morphological traits⁶⁹. This may be equally true for the genes involved in morphogenesis of the humerus across cetacean lineages. PTCH1, which has undergone accelerated evolution in both the Balaena mysticetus and Hyperoodon ampullatus lineages, is a crucial component of the hedgehog signaling pathway. It encodes the Sonic hedgehog (SHH) receptor, the transmembrane protein Patched 1⁷⁰, and forms a negative feedback system in Hedgehog signaling⁷¹. Sun et al.⁷² also discussed its role in cetacean flipper formation. In experiments, heterozygous PTCH1-deficient mice were characterized by increased bone mass and density, with enhanced osteoblast differentiation⁷³. Moreover, another experiment in which both copies of PTCH1 were removed indicated bone overgrowth but decreased density in the bones, especially the humerus⁷⁴. GLI3, another gene known for negatively regulating SHH⁷³, showed accelerated evolution, characterized by the significantly higher ω value in odontocetes, and two codons under positive selection in Kogia breviceps. There was also a significant increase in the ω ratio for PBX3 in the Mysticeti clade. PBX3, typically in combination with other Pbx genes, regulates humerus morphogenesis, with deficient mutants displaying fusions and expansions in the proximal humerus and humeral head duplications⁴⁹. Additionally, Pbx genes are well known for their role in pelvic girdle formation⁷⁵. Interestingly, it was previously discussed in the context of evolution of cetacean hindlimbs⁷⁶ that the deformations in Pbx and GLI3 mutants are similar to the hindlimb skeletal traits of living and extinct whales, including the anatomy of the femur which is the stylopodial element, as the humerus. It is not unlikely that changes in the regulatory genes related to hindlimb transformation (which finally led to its almost full reduction) could equally affect the forelimb and cause its transformations as a byproduct of hindlimb adaptations.

In addition, a stop codon was found in the PANX3 gene sequence of Delphinapterus leucas. PANX3 is highly expressed in skeletal tissues and plays a critical role in bone development, promoting osteoblast proliferation and differentiation: its alterations severely affect the anatomy of the mammalian humerus⁴⁶. However, the impact of this mutation is not clear, as there are no obvious bone structure abnormalities in beluga whales that can be attributed to the change of function of PANX3. Also, we found no signs of relaxation or positive selection in the HOXA9, HOXD9, and HOXD10 genes. The Hox genes, in the context of flipper formation, were investigated previously by other authors on various samples^69,77,78, and no signs of positive selection were discovered for HOXA9, HOXD9, and HOXD10 in those studies.

Therefore, even a slightly relaxed simultaneous selection in several regulatory genes along with site-level positive selection could solely drive changes of a developmental framework for humerus evolution in cetaceans under the Brownian motion model^42,63,79. However, the details of this process are to be studied in depth, as well as the mechanism enabling another aspect of the humerus evolution, namely the structural integration⁸⁰.

Conclusion

Changes in humerus shape during the evolution of cetaceans were so significant that selection toward at least two adaptive optima, following the Ornstein–Uhlenbeck model, was suggested as the pattern best explaining this process. However, the within-landmark variance remained proportional to the rate of neutral evolution or slightly decreased, indicating weak negative selection. The patterns of the humerus modularity showed divergence through deep time: in crown baleen whales, modularity weakened and morphological integration increased, whereas toothed whales retained a modular structure. The main patterns of evolution in genes regulating morphogenesis of the humerus in cetaceans were divergence under a Brownian motion model that corresponded to the relaxed negative selection, site-level positive selection and point mutations. However, due to a past transition of cetaceans to an aquatic way of life their humerus shape also evolved toward increasing aquatic adaptations, with parallel changes in several lineages of crown baleen and toothed whales, seen from the Neogene to modernity. Random drift (probably combined with weak selection) is suggested as a common pattern of humerus evolution in cetaceans. This pattern well illustrates the “fly in a tube model” and shows its importance as an environment for emerging evolutionary innovations.

Material and methods

We sampled 36 cetacean taxa (Supplementary Table S1) which included Protocetidae, Basilosauridae, extinct and extant baleen whales (stem taxa, Balaenidae, Cetotheriidae, Balaenopteridae, Eschrichtiidae), and extinct and extant toothed whales (stem taxa, Physeteridae, Kogiidae, Platanistidae, Kentriodontidae, Phocoenidae, Monodontidae, Delphinidae). The taxa were grouped as amphibious Protocetidae (Maiacetus, Aegicetus), early aquatic taxa (Paleogene Pelagiceti: Basilosauridae and basal Mysticeti and Odontoceti), Neogene and extant Mysticeti and Odontoceti (Fig. 2).

Cetacean humeri from museum collections (Supplementary Table S1) were digitized using 3D surface scanners Artec Eva (resolution 0.5 mm) and Artec Space Spider (resolution 0.1 mm) and modeled using the Artec software⁸¹. Also, 3D models of humeri of three taxa (Maiacetus innuus, Aegicetus gehennae and Dorudon atrox) were downloaded from the University of Michigan Online Repository of Fossils⁸², and the Cynthiacetus peruvianus skeleton model was downloaded from the MorphoMuseum⁸³.

Statistics and reproducibility

All the data were analyzed with publicly available software, procedures and codes, including R packages fully described and referred to below.

The shape of the humerus was described by 34 landmarks (Supplementary Table S3; Fig. 1; Supplementary Data S1), digitized with the software MorphoDig⁸⁴. A Generalized Procrustes Alignment (GPA) was performed with the gpagen function of R v. 4.4.1⁸⁵ package geomorph v. 4.0.8⁸⁶. We divided landmarks into subsets to represent main anatomical parts of the humerus. Landmark partition for modularity analysis was suggested as following: proximal epiphysis (landmarks 1, 2, 3, 4, 5, 6, 7, 8, 31, 32, 33), distal epiphysis (landmarks 22, 23, 24, 25, 26, 27, 28, 29, 34) and diaphysis (landmarks 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 30). For comparison, another approach of shape digitizing was also used, based on 500 sliding pseudolandmarks, evenly distributed on 3D models surface using the R package geomorph functions buildtemplate and digitsurface (Supplementary Fig. S1). They were slid with minimum bending energy with function slider3d, argument stepsize = 0.15, from R package Morpho v. 2.12⁸⁷. The comparison between 34 and 500 landmarks was conducted with the Procrustes test, function protest of R package vegan v. 2.6–4⁸⁸. Both sets, of 34 landmarks and of 500 surface sliding landmarks, showed similar results (null hypothesis of random similarities was rejected, p < 0.004, Supplementary Fig. S1), and therefore the set of 34 landmarks was chosen as the main framework for analyses.

Proximo-distal rotation angle of the humeral head (landmarks 4, 5, 6, 7, 8, 31, 32) was estimated in all species, and the configuration of landmarks of Maiacetus inuus was taken as the reference. Raw coordinates were used. Rotation matrix (presented as cosine) in three dimensions was calculated with the function rotonto of the R package Morpho. We chose the proximodistal axis, applied the inverse cosine to these values, and converted radians to degrees. The function contMap (R package phytools v. 1.5-1⁸⁹) was used to map the rotation angles on the phylogeny.

Phylogenetically informed morphological analyses were computed in R package geomorph⁸⁶. Cetacean phylogeny based on Lloyd & Slater⁹⁰ (Cetacea Safe Median tree) was used for the analysis (Fig. 2). Ten taxa analyzed here (Aegicetus gehennae, MHN-UABCS-Ec/Sj29/42b, Ferecetotherium kelloggi, Kelloggia barbarus, Kentriodon fuchsii, Kentriodontidae indet., Microzeuglodon caucasicum, Oligodelphis azerbajdzanicus, Pachyacanthus suessii, Sulakocetus dagestanicus) were absent in the tree, and we assigned them to their presumed phylogenetic positions based on the metatree by Lloyd & Slater⁹⁰ (the only published tree, which was complete enough in the taxonomic coverage for our research) and specific taxonomic studies³⁸. This full tree was used with caution, since many taxa positions were based solely on taxonomic opinions from the Paleobiology Database (https://paleobiodb.org).

Phylogenetic signal and phenotypic variation

The phylogenetic signal was calculated for the whole set of GPA-aligned landmarks and for each module (see Modularity analysis below), and multivariate Blomberg’s K values were estimated in the geomorph package. Phylogenetically corrected GLS-centered residuals were calculated as Z matrix⁹¹ and they were further used in the analysis. The need for size correction was estimated with the function procD.lm (geomorph); a natural logarithm of the centroid size of the 34-landmarks configurations was used for the calculations. The shape of the cetacean humerus was poorly explained by the size; therefore, we did not additionally perform size correction.

Collyer & Adams⁹¹ suggested that by comparing the different ordinations (statistical technique to project multivariate dataset into low, i.e., 2D dimensional space) one can assume sources of shape variation. Thus, phenotypic variation was summarized by standard principal component analysis (PCA), phylogenetic principal component analysis (phyPCA), and phylogenetically aligned component analysis (PaCA)⁹¹. All these ordination methods described most of the phenotypic variation in their first two components. Shape deformations were visualized using the Ziphius cavirostris humerus, the specimen which was the closest to the mean shape in the GPA-aligned 34-landmark configuration. Predicted coordinates of the extreme (minimal and maximal values of the axis) shapes of the principal components analyses were used to warp 3D mesh of the mean humerus with thin-plate spline method (functions warpRefMesh and plotRefToTarget in the geomorph package).

Standard PCA based on ordinary least squares (OLS) showed overall variation irrespective of the sources of variation. Phylogenetic PCA (phyPCA) was run on the covariance matrix corrected for phylogeny with the generalized least squares method (GLS); its output was presented as transformed GLS mean-centered residuals, which should be independent of phylogeny. Thus, major phyPCA components maximized variation that was independent of phylogeny. Phylogenetically aligned component analysis (PaCA) based on OLS was similar to two-block partial least squares analysis, but one of its matrices was the phylogenetic covariance matrix; therefore, PaCA maximized variation best describing the phylogenetic signal. Comparison of the morphospaces under three analyses may help to find phylogenetic and functional signals in the data, as similar directions of eigenvectors in phyPCA and PaCA might point to a correlation between ecological and phylogenetic factors⁹¹.

Comparisons of the ordination results were made by measuring angles between eigenvectors of the components (function angleTest in Morpho). If two components showed opposite directions (the angle between them larger than 90 degrees), we rotated one of the components, and thus the significance of the test showed the probability of differences irrespective of the direction. Two components were correlated when the angle was small, while the right angle was the evidence for the absence of correlation.

Contribution of individual landmarks to the principal components was also estimated. A contribution of each variable was measured as $\frac{{\left({X}_{i}\cdot {{sdev}}_{i}\right)}^{2}}{\sum {\left({X}_{i}\cdot {{sdev}}_{i}\right)}^{2}}$, where $X$ was an eigenvector of the $i$ component; ${sdev}$ was an eigenvalue. The contribution ranged from 0 to 1 and showed the importance of the variable in the component. Contributions of respective x, y, and z coordinates were summed for each landmark.

The Spearman rank-order correlation coefficient (r_s) was used to measure contribution of the humeral head twist to the morphospace of the principal components.

Standard PC scores were also used as the complex continuous traits. Probable adaptive shifts in the PC1 and PC2 were evaluated with the reversible-jump Bayesian method of fitting multi-optima Ornstein–Uhlenbeck model (R package bayou v. 2.3.0)⁹². We used priors: dists=list(dalpha = “dhalfcauchy”, dsig2 = “dhalfcauchy”, dsb = “dsb”, dk = “cdpois”, dtheta = “dnorm”), param=list(dalpha=list(scale=0.1), dsig2=list(scale=0.1), dk=list(lambda=10, kmax=50), dsb=list(bmax=Inf, prob= tree$edge.length), dtheta=list(mean=mean(PCvalues), sd=1.5*sd(PCvalues)). The MCMC chains for each principal component were generated with the standard error of the corresponding PC. Totally, 500,000 generations were run, then we cut off 30% of the recorded runs. The effective size for each parameter was checked for being larger than 200.

Since our sample of taxa was relatively small, the evolutionary dynamics of the humerus shape was estimated with the first three standard principal components. They summarized almost 75% of shape variation. We standardized them to zero mean and unit standard deviation. Different evolutionary models were fitted (in R package mvMORPH v. 1.2.1): Brownian motion, Brownian motion with different rates for different regimes, one-peak and multi-peak Ornstein–Uhlenbeck processes. The following parameters were set for the function mvOU: root=T, decomp = “cholesky”, decompSigma = “diagonal”. We analysed different optima formulations: two-optima models as (1) Neogene and extant Odontoceti (as at Fig. 2) vs other taxa, (2) Neogene and extant Mysticeti (as at Fig. 2) vs other taxa, (3) Delphinoidea (Kentriodontidae, Delphinidae, Monodontidae and Phocoenidae) vs other taxa, (4) all Mysticeti (the stem MHN − UABCS−Ec/Sj29/42b and all other baleen whales) vs other taxa; three-peaks models as (5) combination of models 1 and 2; (6) combination of models 3 and 4 and Physeteroidea combined with Mysticeti. Groups in models 3 and 4 were formulated after the Bayesian reversible-jump multi-optima Ornstein–Uhlenbeck model (bayou analysis⁹²). The fit of the models was compared by the weights of the Akaike coefficients corrected for the number of parameters (AICc). The multivariate phylogenetic half-life⁹³ (t_0.5, the time for the trait to evolve halfway from the ancestral state to the new optima) was calculated. Multivariate t_0.5 is calculated from the eigenvalues of the selection matrix. We expressed it in terms of the tree height percent. For pure Brownian motion process t_0.5 = Inf. The faster the adaptation process is, the closer t_0.5 is to zero.

Modularity analysis

The maximum likelihood method (R package EMMLiv2 v. 0.0.3) was used for testing several partitions of landmarks and choosing the best fit model of module assignment⁹⁴. The hypothesis of no module was supported when optimal estimates of between-landmarks correlations were equal to each other, and the structure could be analyzed as a single entity. The complexity level was assessed by comparison of Akaike coefficients (AIC_c) corrected for the number of parameters. EMMLi analysis was conducted on phylogenetically corrected GLS-centered residuals (Z matrix estimated in previous section). The input correlation matrix was calculated as a matrix of congruence coefficients of the 3D landmarks dataset (paleomorph function dotcorr⁹⁵). It was not fully positive (it had a single negative eigenvalue), so it was further transformed to the nearest positive definite matrix using the function nearPD of the R package Matrix v. 1.6-0⁹⁶. Several alternative models were tested: no modules, all combinations of two module structures (each part of the humerus vs other parts), and three modules (proximal epiphysis, distal epiphysis, and diaphysis). Different parameterizations were applied to each modularity model: same or different within- and between-module correlations. Thus, two-module arrangements had two parameterizations, and the three-module arrangements were described with four variants. In total, eleven models were tested with EMMLi. Optimal, maximum likelihood estimation of the true correlation between landmarks, was estimated for each model.

The strength of partitions was evaluated with the values of within- and between-module correlations (ρ): the closer within-module ρ to one, the more the contrast. Low-contrast correlation structure was considered when within-module ρ was less than 0.5. Modularity was analyzed for total sample and for each subgroup with a sufficiently large sample: Paleogene Pelagiceti, Neogene and extant Mysticeti, Neogene and extant Odontoceti, “Delphinoidea” (Kentriodontidae, Delphinidae, Monodontidae and Phocoenidae), and “All Mysticeti (including their stem and crown members) + Physeteroidea”. Procrustes transformation and GLS-phylogenetic correction were applied for each group separately before the analysis.

Modularity and phenotypic variation could be interconnected⁵⁹. A separate Procrustes transformation was conducted to align partitions to the unit centroid for each module. Total morphological disparity for each module was calculated from only-intercept phyANOVA (function procD.pgls of R package geomorph). After that, we estimated the morphological variation as the sum of the Procrustes variances of the phylogenetically corrected landmarks; function morphol.disparity in geomorph with phyANOVA-object as input and argument transform = TRUE to remove the variation, which was related to phylogeny. To distinguish low and high disparities we compared observed values to the morphological disparities of 10000 random partitions⁵⁹. A high morphological disparity was considered if the empirical value was above the 95%-quantile of random distribution.

Evolutionary effect of integration

Module’s net evolutionary rates under the Brownian mode of evolution were compared with function compare.multi.evol.rates of R package geomorph. Under the Brownian motion model, the starting value of traits did not change over time, but the variance of the traits increased proportionally with time. And the higher the rate of evolution (σ ²), the higher was the variance of the trait. We used Brownian motion as a model for a random walk.

Constrained or facilitation scenarios of the evolutionary shape change were suggested from the analysis of morphological disparity and evolutionary rates correlation. The Procrustes’ variance and the evolutionary rate under neutral evolution were estimated for each landmark⁵⁰. The variance was calculated as the sum of x, y, z Procrustes coordinates variances. The evolutionary rate per landmark was calculated as the sum of σ² values^50,97 in hot.dots v. 0.0.0.9000 R library (https://github.com/rnfelice/hot.dots). Simulations were used to estimate difference between expected and observed rate to variance correlation. Neutral expectation variance of each landmark was calculated with sim.char function of R package geiger v. 2.0.11⁹⁸. Mean values of variances of 10000 simulations were used for the linear regression analysis. The ANCOVA type of the linear models with unique (type II) sum of squares was used to estimate differences between slopes of the observed and simulated rates on variance regression for the pooled landmarks and for each module. We considered different slopes when the interaction term was significant. The significance of regression estimates was calculated with the permutation test, function lmp of the R package lmPerm v. 2.1.0⁹⁹. Traits were not centered. The maximum number of iterations was 50000, and the criterion for the termination of iterations was “Prob”, i.e. p-values were sampled from all possible permutations, and it was stopped when the standard error of the sampled p-value fell below 0.001 fraction of its estimate.

In addition to ANCOVA, we performed a graphical analysis of the observed and simulated regression lines. We calculated a global rank envelope test with extreme rank length (ERL) measure in the R package GET v. 1.0^100,101. The approach was used in Murrell¹⁰² to analyze disparity through time. The non-parametric rank envelope test was used to compare a set of simulated curves to the observed curve, and here we took fitted regressions lines for the comparison. All of the curves were ordered and ranked. The null hypothesis (i.e., the equality of the simulated and observed curves) was rejected if the observed curve fell outside of the bounding curves (at the 0.05 significance level).

Gene analysis

Analysis for positive selection was conducted for 11 genes HOXA9, HOXD9, HOXD10, PANX3, GLI3, ALX4, PTCH1, PTCH2, PBX1, PBX2, PBX3, which regulate morphogenesis of the humerus^{46,47,79,103,104,105}. A total of 45 cetacean species and 2 outgroup mammals (Hippopotamus amphibius and Choeropsis liberiensis) were used in the analysis. Annotated coding sequences were retrieved from NCBI¹⁰⁶. For 35 non-annotated species, we performed a BLAST search¹⁰⁷ on assemblies downloaded from NCBI and DNAZOO databases. Then, the coding sequences were aligned in PRANK software¹⁰⁸. We also translated them into protein sequences to visually assess for insertions/deletions and correct possible misaligned regions.

During the data preparation and alignment, we identified a stop codon (codon 333, CGA → TGA) in the PANX3 sequence within the Delphinapterus leucas genome. To verify this finding, we examined two other published assemblies (GCA_002288925.3 and GCA_029941455.3) of the beluga whale and confirmed the presence of the same stop codon. Consequently, this species was excluded from subsequent analyses of the PANX3 gene.

We used the codeml algorithm implemented in PAML 4 software¹⁰⁹ to calculate the ratio of nonsynonymous to synonymous substitution rates also known as omega ratio (ω). The ω ratio, defined as ω = dN/dS, signifies alterations in selective pressures. Specifically, ω > 1 indicates positive selection, ω < 1 purifying selection, and ω = 1 neutral evolution. A well-supported tree topology from McGowen et al.¹¹⁰ was used as an input in PAML analysis. To assess if rapid evolution and positive selection were confined to specific lineages, we used the branch two-ratio model, which allows for different ω ratios between foreground branches and background branches. We then defined specific lineages of interest to measure if selection pressure acting on them is different from the background. First we labeled the Cetacea ancestral branch as the foreground to find if there are general trends in the selection of target genes, outlining the transition of whales to a fully aquatic lifestyle. Other test branches were chosen based on the results of our analysis of phenotypic variation (see the Phenotypic variation and phylogenetic signal section): Balaena mysticetus, Balaenoptera musculus, Monodon monoceros, Phocoena sinus, Tursiops truncatus, and Stenella coeruleoalba as the outliers in phyPCA analysis, Orcinus orca as a species that exemplified the positive end of the PC1 axis, Physeter macrocephalus, Kogia breviceps and Hyperoodon ampullatus as characteristic species of corresponding families Physeteridae, Kogiidae and Ziphiidae; Mysticeti and Odontoceti ancestral branches as summarizing different phylogenetic trends in the shape of humerus. In total 13 lineages were labeled as the foreground and tested one by one to evaluate the ω values. Next, we performed the same branch analyses but using the whole Odontoceti and Mysticeti clades as foregrounds to compare trends in evolution between these clades. Additionally, the branch-site model was used to find amino acid sites that might be under positive selection in chosen lineages. The Bayes Empirical Bayes method¹¹¹ was then applied to identify codon sites that could be under positive selection on the foreground branches.

Statistical significance was evaluated using likelihood ratio tests (LRTs); for branch-site tests, we assumed that the null distribution was a 50:50 mixture of a χ² distribution with a point mass at zero. Obtained p values were corrected with the Benjamini–Hochberg procedure to control for the false discovery rates (FDR).

Reporting summary

Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article.