Ecomorphology of toothed whales (Cetacea, Odontoceti) as revealed by 3D skull geometry

Extant odontocetes (toothed whales) exhibit differences in body size and brain mass, biosonar mode, feeding strategies, and diving and habitat adaptations. Strong selective pressures associated with these factors have likely contributed to the morphological diversification of their skull. Here, we used 3D landmark geometric morphometric data from the skulls of 60 out of ~ 72 extant odontocete species and a well-supported phylogenetic tree to test whether size and shape variation are associated with ecological adaptations at an interspecific scale. Odontocete skull morphology exhibited a significant phylogenetic signal, with skull size showing stronger signal than shape. After accounting for phylogeny, significant associations were detected between skull size and biosonar mode, body length, brain and body mass, maximum and minimum prey size, and maximum peak frequency. Brain mass was also strongly correlated with skull shape together with surface temperature and average and minimum prey size. When asymmetric and symmetric components of shape were analysed separately, a significant correlation was detected between sea surface temperature and both symmetric and asymmetric components of skull shape, and between diving ecology and the asymmetric component. Skull shape variation of odontocetes was strongly influenced by evolutionary allometry but most of the associations with ecological variables were not supported after phylogenetic correction. This suggests that ecomorphological feeding adaptations vary more between, rather than within, odontocete families, and functional anatomical patterns across odontocete clades are canalised by size constraints.


Introduction
Cetaceans are a monophyletic group of aquatic mammals and comprise two clades: Odontoceti (toothed whales) and Mysticeti (baleen whales). Odontocetes are pursuit predators that use echolocation to locate prey, while mysticetes feed on large aggregations of smaller organisms using their sievelike baleen. Odontocetes diverged from mysticetes about 34 Mya (Marx et al. 2016) and are much more diverse; they comprise 10 extant families, including at least 72 species and 32 genera, with the families Delphinidae (oceanic dolphins) and Ziphiidae (beaked whales) showing the largest number of species (Hooker 2009). Compared to mysticetes, this group displays wide variation in body size and skull morphology related to feeding ecology (Werth 1992(Werth , 2006aCozzi et al. 2016). Skull synapomorphies of Odontoceti include the presence of a concave facial area to accommodate the melon (a fatty organ functioning as a lens for the propagation of echolocation sounds), the presence of premaxillary foramina and premaxillary sac fossa, posterior expansion of maxilla over the supraorbital region covering the frontal bones, and directional facial asymmetry, all features that are related to their highly specialized sonar system (Cranford et al. 1996;Uhen 2004;McKenna et al. 2012;Geisler et al. 2014;Marx and Fordyce 2015;Marx et al. 2016;Martínez-Cáceres et al. 2017;Churchill et al. 2018;Coombs et al. 2020Coombs et al. , 2022. Due to such a highly specialised system, the odontocete skull has received particular attention in studies of functional morphology (Rommel et al. 2009;Mourlam and Orliac 2017;Gillet et al. 2019;Park et al. 2019;Roston and Roth 2019;Coombs et al. 2020Coombs et al. , 2022. One general evolutionary trend in odontocetes since the Oligocene is the telescoping of the skull (Churchill et al. 2018). Winge (1918) introduced this term and Miller (1923) subsequently formalised it to describe the unique posterior elongation of the rostral elements relative to the backward shift of the bony nares. It is now recommended that the term be used more specifically to refer to the maxillo-occipital proximity provided by the extensive bone overlap of cranial elements in the skull (Roston and Roth 2019). Such a dramatic anatomical innovation resulted in an overlap of adjacent facial bones in a manner reminiscent of an antique folding telescope (Winge 1918(Winge , 1921Miller 1923;Rommel et al. 2009;Roston and Roth 2019). This provides improved functionality of the toothed whale skull in an aquatic environment and is accompanied by the evolution of echolocation and associated changes in the ear (Geisler et al., 2014;Mourlam and Orliac 2017), increased hydrodynamicity (Fordyce and de Muizon 2001;Marx et al. 2016), changes in thermoregulation (Manger 2006;Werth, 2007;Marino et al. 2008), and eye pigmentation (Fasick and Robinson 2016), as well as the evolution of new oxygen storage units such as different types of myoglobin and respiratory pigments that increase the oxygen carrying capacity (Noren and Williams 2000;McClellan et al. 2005).
The skull of odontocetes is also shaped by sensory, cognitive and feeding functions. The parts of the skull that are more relevant to feeding are the rostrum, the temporal fossa (origin of the temporalis muscle, the main adductor of the jaws in odontocetes), the zygomatic process of the squamosal bone (point of articulation with the mandible), the teeth, the hyoid apparatus, and the mandible (Perrin 1975). The latter has received more attention in the literature due to its importance in both feeding strategies and the transmission of sounds to the ear (Werth 2006a;Rommel et al. 2009). In addition, the presence of the melon greatly influences the external shape of the head in odontocetes. Indirectly, skull morphology can provide ecological and trophic information, as the skull assists toothed whales in exploiting surrounding resources (Marshall 2009).
Previous studies of odontocetes demonstrated that interspecific differences in skull size are associated with reproductive parameters, peak auditory frequencies, and prey size rather than dietary categories (MacLeod et al. 2007;Loy et al. 2011;Barroso et al. 2012;Guidarelli et al. 2014Guidarelli et al. , 2018McCurry et al. 2017b), and a link between skull size and deep diving abilities has been described (MacLeod et al. 2006;McCurry and Pyenson 2019). Nevertheless, a comprehensive ecomorphological analysis of odontocete skull evolution within a phylogenetic framework is still lacking.
It is expected that skull shape and size correlate with ecological factors. As ecological factors such as predation have been suggested to be one of the major selective pressures that drove the evolution of biosonar mode in toothed whales (Galatius et al. 2018;Jensen et al. 2018), a strong relationship can be predicted between skull size and biosonar mode. Since mammalian skull morphologies appear to show a strong phylogenetic signal, it is also possible that this expected association might eventually be obscured after clade evolutionary history is accounted for (Marcus et al. 2000). Comparative methods like the phylogenetic generalised least squares (PGLS) was introduced to overcome this issue (Martins and Hansen 1997), especially in conjunction with geometric morphometric (GM) data (Monteiro 2013). GM is an ideal tool to study skull morphological evolution across ecomorphologically diverse vertebrate groups since it allows separation of size from shape variation, simultaneously providing effective visualization (Rohlf and Marcus 1993;Cardini and Polly 2013). Here, we investigated 60 out of 72 species of toothed whales covering 32 living genera (~ 86% of the current diversity) to test for association between the evolution and diversification of odontocete skull shape and size and various biotic and abiotic factors including diet (MacLeod et al. 2007;Slater et al. 2010;McCurry et al. 2017a;Galatius et al. 2020;Coombs et al. 2022), biosonar mode and maximum/minimum peak auditory frequencies (kHz) (Barroso et al. 2012;Galatius et al. 2018;Jensen et al. 2018), diving ecology (Noren andWilliams 2000;Werth 2006a;Würsig 2009), prey size (minimum, maximum, andaverage) (MacLeod et al. 2006), surface temperature (ST), encephalization quotient (EQ) and brain mass (Montgomery et al. 2013). Our study was designed within a comparative framework in order to determine the relative contributions of ecological adaptations to skull size and shape and therefore to reveal any potentially adaptive ecological variation.

Specimens examined
Skull data were collected by DV from 111 individual toothed whale specimens representing 60 species (1-5 individuals per species) of nine families: Delphinidae, Iniidae, Kogiidae, Lipotidae, Monodontidae, Phocoenidae, Platanistidae, Pontoporidae, and Zaiphiidae. Only adult specimens were used for this study; adult status of delphinids was confirmed by the caudal extension of the maxilla to the nuchal crest and the visibility of the frontal bones in dorsal view (Cozzi et al. 2016), whereas non-delphinid adult specimens were selected comparing their skull size with total body length measurements in the literature and in the museum's database. The measured specimens were from the collections of the following institutions: Natural History Museum of London ((NHM), Muséum National d'Histoire Naturelle (MNHN, Florence, Italy), Natural History Museum of the University of Pisa (MSNUP), La Specola (MZUF), National Museums Liverpool (NML, Liverpool, UK) and the Natural History Museum of Denmark (NHMD). All measured specimens are listed in Online Resource 1. In addition, 3D models were downloaded from the website Phenome 10 k (Goswami 2015) to cover more individuals per species as well as the species Orcaella heinsohni.

Sampling: Photogrammetry protocol
A 3D model for each skull specimen was reconstructed from a set of photographs taken in three different orientations (ventral, dorsal, and lateral). Photographs were taken using a digital camera Canon EOS 1100D 12.2-megapixel resolution digital single-lens reflex camera with 18-55 mm lens, on a tripod. Specimens were fixed vertically on a rotating table and ~ 100-150 photographs were taken at intervals of approximatively 10 degrees. For larger specimens (i.e., Hyperoodon spp., Ziphius cavirostris, and Indopacetus pacificus) the operator (DV) moved the tripod with a mounted camera around the object placed on a pad on the floor, and the "walk-around method" was used (Mallison and Wings 2014). Millimetre scale measurements (rostrum length and bizygomatic width) were taken for further scaling reference of the virtual models. Images were imported into Agisoft PhotoScan Professional (Agisoft 2018) and photos from each chunk/orientation were aligned in order to generate three dense point clouds (ventral, dorsal and vertical orientation) that were subsequently merged together (Falkingham 2012;Katz and Friess 2014;Mallison and Wings 2014;Evin et al. 2016;Linder 2016;Muñoz-Muñoz et al. 2016;Agisoft 2018;Mallison 2018). 3D models with texture were exported as. PLY files and scaled by dividing the scaling factor identified in Meshlab by the scale measurements (in mm) (Cignoni et al. 2008a(Cignoni et al. , 2008b.

Landmark data
A total of 28 landmarks ( Fig. 1; Table 1; raw data and scripts are available at: https:// github. com/ Vicar iD/ Macro-Data. Vicari-et-al. git) were taken on skull virtual models using IDAV Landmark software (Wiley et al. 2005). The landmarks were selected based on previous studies (Churchill et al. 2018;Galatius et al. 2020;Vicari et al. 2022a, b) to describe broad anatomical skull regions relevant to the research questions, including the rostrum, the temporal region, and the facial concavity. The landmark scheme used in this study was assessed on the mean shape dataset using the Landmark Sampling Error Curve function (LaSEC) within the LaMBDA R package (Watanabe, 2018). This function subsamples the dataset by randomly selecting three landmarks, and generating a median fit obtained shape (= Procrustes) distances between the original dataset and the subsampled one for each selected number of iterations. Results showed that 11 landmarks provided a median fit value of 0.90 and 15 provided a median fit value of 0.95, confirming that our 28-landmark scheme was dense enough to achieve a fit R value of 1 that characterised both toothed whale shape and size. Size was quantified from the original 3D landmark coordinates as centroid size (CS), which is the square root of the sum of the squared distances between each landmark and the centroid of the entire configuration, while shape data were extracted using the generalised Procrustes analysis (GPA; Bookstein 1991).
Generalized Procrustes analysis is an iterative procedure where variation in size is first removed by scaling each landmark configuration so that it has a CS of 1.0; rotation and translation are taken into account by centring and rotating the original landmark coordinates in order to obtain an optimal solution that minimizes the quadratic distances between homologous points (Procrustes method). The thin plate spline (TPS) (Gunz et al. 2009) was used to reconstruct the positions of missing landmarks using individuals belonging to the same species and/or genus (when not possible within the same species) as reference specimens. Missing  Table 1 1 3 landmarks were identified in 34 specimens (Table 1), and in most cases, they were concentrated in the pterygoids. To detect their impact on subsequent analyses, sensitivity analyses were performed on two separate datasets, one including 28 landmarks and another with a smaller number of landmarks (n = 26).

Skull asymmetry and measurement error
Landmark recording was taken twice in four different sessions on 111 individuals in order to quantify skull directional asymmetry (Klingenberg et al. 2002) and, indirectly, measurement error (Fruciano 2016), repeatability, and digitizing ability of the operator (DV). A Procrustes ANOVA was performed using MorphoJ on the shape component and the percentage of variance explained by the measurement error was calculated by examining Mean Square (MS) values following the formula described in Sherratt (2015) and Fruciano (2016). The plotOutliers function in the package geomorph vs 4.0.3 was used (Adams et al. 2016) to reveal outliers, and misplaced or inaccurate landmarks (Fruciano 2016). Principal Component Analysis (PCA) in R Team (2015) using geomorph package (Adams et al. 2016), the Unweighted Pair Group Method with Arithmetic Mean (UPGMA) (Cardini 2014) in PAST 2.17 (Hammer et al. 2001), and Mantel tests in MorphoJ (Klingenberg et al. 2011) were performed to identify individuals adversely affected by measurement error (Online Resources 2 and 3).

Geometric morphometric (GM) analyses
The registered shape coordinates and log transformed CS were averaged for each species for subsequent macroevolutionary analyses (Meloro and Tamagnini 2021). Geometric Morphometrics (GM) permits partitioning of the asymmetric and symmetric components of shape variation (Klingenberg et al. 2002). As many species of toothed whales show a high degree of asymmetry in their crania (MacLeod 2002;MacLeod et al. 2007;Fahlke et al. 2011;Galatius and Goodall 2016;del Castillo et al. 2017;Huggenberger et al. 2017;Coombs et al. 2020) and as the asymmetric component is relevant to answer the intended research questions, these variables were partitioned using the function bilateral.symmetry in geomorph package (Adams et al. 2016). Analyses were also performed on the whole skull Table 1 Description of landmarks taken on odontocete skulls for our GM analysis

Area
No.

Landmark description
Facial region: area into which the maxilla and premaxilla expand during telescoping; it is bounded posteriorly by the nuchal crest and laterally by the orbitotemporal crest

1-2
Anterior tips of the right and left premaxillae 3-4 Point between the maxillary flange and the antorbital notch on the right and left lateral portions of the maxilla 5 Septum nasi osseum 6 Apex of the nuchal crest or lambdoid crest (posteriormost structure) in the midline of the skull/dorsomedial margin of the supraoccipital at the intersection of this margin and the external occipital crest Planum parietale: or temporal fossa, the origin of the temporalis muscle. This area is bounded by the parietal, squamosal, frontal, and alisphenoid bones

7-8
Junction of nuchal crest, temporal, parietal, occipital and frontal sutures on the dorsal border of the temporal fossa 9-10 Posteriormost point on the temporal crest/on the curve of the shape data, without partitioning into symmetric and asymmetric components.
A PCA was performed on the whole mean shape (without splitting the shape into symmetric and asymmetric component) of the Procrustes coordinates to identify patterns of variation between species at a macroevolutionary scale (Reznick and Ricklefs 2009).
We examined four categorical ecological variables (surface temperature, diet, dive depth, biosonar mode) and nine continuous variables (adult average body mass, adult average body length, maximum hearing frequency minimum hearing frequency, encephalization quotient, absolute brain mass and mean, maximum, and minimum prey size) (Online Resource 4). Surface temperature (ST) is defined as the surface temperature at the area of the maximum abundance for each species and included the following categories: warm, temperate, cold, and cold-temperate as defined in Wursig and Perrin (2009). Diet was defined following Slater et al. (2010) and then grouped into three categories: fish, squid, and fish/mammals eaters. Diving depth was quantified by compiling diving range depths from Galatius et al. (2020) and Würsig and Perrin (2009), and comparing these with Würsig (2009). Species were additionally grouped into three diving categories: deep (> 600 m), semipelagic (between 100 m and 500 m), and shallow (< 100 m). Biosonar mode was defined following Surlykke et al. (2014) and Jensen et al. (2018) as broad-band (BB), narrow-band highfrequency (NBHF), and frequency modulated (FM). As the sperm whale (Physeter macrocephalus Linnaeus,1758) is not present in the dataset, the multi-pulsed (MP) category was not considered. Body mass (kg) and length (m) were calculated as mature adult average measures (Online Resource 4). Frequencies at maximum and minimum energy (dB) of the echolocation clicks (kHz) were considered following Jensen et al (2018). The encephalization quotient (EQ) was included to quantify the variation in brain mass not explained by the allometric component (Montgomery et al. 2013). Absolute brain mass and EQ were compiled from (Montgomery et al. 2013). Mean, minimum, and maximum prey size were compiled from MacLeod et al. (2007).
The Procrustes ANOVA (function "procD.lm") in geomorph (Adams et al. 2016) was employed to test for allometry (with log CS as the predictor variable and Procrustes shape coordinates as the response) and to test the association between skull size, whole shape, symmetric and asymmetric shape against the variables listed above. Echolocation peak frequencies, prey size, EQ and brain mass were not available for all the species, so analyses were run on four separate datasets with 60, 56, 31, and 26 species respectively (taxa included in each dataset can be found in Table 2 and Online Resource 5).

Comparative methods
Species-level phenotypic data are not statistically independent due to their phylogenetic history and so the phylogenetic component of the interspecific variation must be estimated. To test if level of species similarity differs with respect to phenotypic traits, we quantified phylogenetic signal in skull data (both size and shape) using the K statistic (see section below). Secondly, the phylogenetic generalised least squares (PGLS) (Rohlf 2001) approach was implemented using the function "procD.pgls" (Adams and Collyer 2015). This allowed incorporation of the phylogenetic covariance matrix as an error term in the Procrustes ANOVA models.
Our analyses used the time-calibrated phylogenetic tree from McGowen et al. (2009) (Fig. 2) as it covered most species in our dataset (selected using the function "drop.tip", R package ape 5.0, Paradis and Schliep 2019), except for Orcaella heinsohni, Sousa plumbea and S. teuszii. These three missing species were added to the topology manually using Mesquite (Maddison and Maddison 2021) by breaking each branch leading to Sousa and Orcaella in half and attaching the missing species; all three Sousa species were treated as a polytomy.
The phylogenetic signal for skull size and shape data was quantified using the "K" statistic (for log CS) and its multivariate equivalent Kmult, applying the function "physignal" in the geomorph package (Blomberg et al. 2003;Adams 2014). Higher K or K mult values represent stronger phylogenetic signal in a trait or character. A value of K or K mult = 1 indicates the trait evolved under Brownian Motion (BM) (Blomberg et al. 2003), while K or K mult < (or > 1) means that closely related species resemble each other either less than or more than expected by BM, respectively (Blomberg et al. 2003).

Skull asymmetry and measurement error
Matrix correlation between replicates supported a strong positive and significant correlation for both skull size and shape in the whole dataset of 28 landmarks (r = 0.99; p < 0.0001). The repeatability index calculated on the mean square of the replicas was 0.92 (Table 3), while it was 0.99 with the function "rep_index" (Marcy et al. 2018). The Procrustes ANOVA (Table 3) showed that both the percentage of variance explained by fluctuating asymmetry (FA) was greater than that explained by directional asymmetry (DA). This is shown in the effects of between-individuals variation on shape as well as measurement "side" representing the directional asymmetry (DA), and interaction between individual and side, representing fluctuating asymmetry (FA, which are both significant (Table 3).

Overall skull shape
The first two principal component (PC) vectors together accounted for 67.5% of the total variance (Fig. 3), while  smaller species are in blue, while larger species are in red. Skull differences across major toothed whale genera are shown PC3 accounted for only 9.7%; the remaining PCs each accounted for ≤ 4% of the total variation, and so will not be further discussed. PC1 (55.0% of var.) described relative changes in the proportion of the rostrum from narrow and elongated (as in 'river dolphins', PC1 negative scores) to short and wide (PC1 positive scores, i.e., Kogiidae). The braincase relative width and length similarly loaded positively on PC1. For PC1 negative scores, the foramen magnum, characterised by 4 landmarks (LM 15,16,17,18), assumed a more circular shape and ventral position. Also, landmarks on the pterygoid hamuli, which delimited the posterior margin of the hard palate and the border of the internal bony nares, shifted more anteriorly as PC1 scores increased. The PC2 (12.6% of var.) axis described changes in the overall area of the temporal fossa and the concavity of the profile of the facial region. PC2 negative values indicated a reduction in the relative size of the temporal fossa, where the temporalis muscle (which elevate the mandible) originate. It also shows the dorsal shift of the unpaired landmarks on the nuchal crest. PC2 positive values characterise a shortening of the pterygoids and an anterior shift of the nasal area as described by landmark 5, together with the anterior shift of the landmarks describing the ventral-most point of the paroccipital process. When the categorical variables were mapped onto the morphospace, the PC2 separated species based on diet, biosonar mode, and diving ecology (Fig. 4), showing a potential pattern of association with skull shape. The lower middle part of the PC scatterplot is occupied by species who are squid-eaters and deep divers, with beaked whales (deep divers that use FM biosonar) specifically occupying the lower left part of the plot. The regression of shape coordinates versus log CS revealed a significant (p = 0.002) allometric component, with size explaining 9.0% of shape variance (Fig. 5). The lower part of Fig. 5 is occupied by small species belonging to the Phocoenidae, while the right side is occupied by larger species belonging to the Ziphiidae. The Delphinidae occupies a broad range of sizes from small dolphins, such as the common dolphin (D. delphis), to large species, such as the killer whale (Orcinus orca).

Ordinary least squares (OLS) and phylogenetic generalized least squares (PGLS)
The Procrustes ANOVA showed a significant impact of nearly all ecological variables on both skull shape and size on the complete dataset (n = 60) (Table 4). However, prey size had no impact on skull shape and size, while surface temperature (ST) had no effect on size. Among all variables, skull shape was most closely associated with biosonar mode and surface temperature, which explained ~ 20% and 15% of the whole shape variance, respectively. Skull size was most closely associated with body length (86% of variation) and body mass (77% of variation).
A significant phylogenetic signal across odontocetes species could be detected, being stronger in skull size (K mult = 0.653, p < 0.001) than skull shape (K mult = 0.565, p < 0.001). After phylogenetic correction, skull size was found to be significantly influenced by biosonar mode (10% var.), brain mass (69% var.), and maximum peak frequencies (42% var.), PGLS also confirmed only body length, brain and body mass, average and minimum prey size, and surface temperature to be significantly associated with skull Table 3 Procrustes ANOVA on skull size log centroid size (CS) and shape component. Directional (DA) and fluctuating asymmetry (FA) explained a considerably small % of total variation that was compara-ble but lower than those explained by replicas. The Rsq of variation explains the contribution of each factor to overall variation. R is the intraclass correlation index * p-values are in bold when significant (p < 0.05) shape with brain mass explaining 13% of variation and all other significant variables explaining around 5% each. Moreover, significance of each of the variables listed above are due more to one component (asymmetric or symmetric) of skull shape than whole (not-partitioned) skull shape (Online Resource 6). Analyses were also conducted using residuals of the regression of whole shape against log CS (Online Resource 7) demonstrating that diet and metric variables were not significantly associated with skull shape, and minimum and mean prey size were instead associated after phylogenetic correction.

Discussion
Foraging underwater has an enormous cost due to the challenging 3D environment and toothed whales have evolved ways of minimizing this cost, for example, evolving echolocation abilities. It is well established that odontocete skulls show extensive asymmetry (Fahlke and Hampe 2015;Cozzi et al. 2016;Coombs et al. 2020), and our results confirm previous expectations on the presence of directional asymmetry also at the interspecific scale (Fahlke and Hampe 2015). Directional asymmetry of the skull in toothed whales is driven by adaptations for high frequency sound production (Cranford et al. 1996;MacLeod et al. 2006;MacLeod et al. 2007;Geisler et al. 2014;Hirose et al. 2015;Park et al. 2019;McCurry et al. 2017b;Churchill et al. 2018;Coombs et al. 2020). The asymmetric shape and presence of specialized fats in the melon allow for the focussing of vocalizations into a highly directional sonar beam for prey echolocation (Surlykke et al. 2014). This system of sound production for echolocation has diversified into distinct forms resulting in varying degrees of skull directional asymmetry within toothed whales (Fahlke et al. 2011;Huggenberger et al. 2017;Coombs et al. 2020). Trade-offs between size, frequencies emitted, and beam directionality are known (Jensen et al. 2018), and our results confirm the correlation between skull size, biosonar mode and maximum peak frequencies. In fact, the narrow-band high frequency biosonar mode appears to have evolved four times in small distinct ecological groups of toothed whales, which also display a difference in size i.e., Kogidae, Phocoenidae, Pontoporiidae, and Lissodelphininae (Surlykke et al. 2014;Galatius and Goodall 2016;Galatius et al. 2018;Jensen et al. 2018).
Our results demonstrated that skull size is strongly associated with brain mass, and with maximum peak frequency, biosonar mode and prey size variables. This was expected as a pattern of body size constraints on sound production for prey detection is well-established (Jensen et al. 2018).
Brain mass occupies the larger portion of the toothed whale skull variation, and it is worth noting that encephalization quotient (EQ) is not correlated with both skull shape and size. Moreover, the impact of body metrics on skull size is evidently a consequence of allometric changes that equally influence skull shape where brain mass explains generally the higher proportion of shape variance.
Sound is produced by high electrical potential, and having a large brain facilitates a greater information processing (Marino et al. 2004;Dudzinski et al. 2009), increases cognitive abilities, social ecology, communication (Montgomery et al. 2013), and the ability to control their hearing to maximize the returning echo (Nachtigall and Supin 2008;Pacini et al. 2011). Echolocation clicks are generated by right and left phonic (or monkey) lips which are both able to produce clicks, then focus them though the melon that acts as an amplifier to detect the prey/ target in the 3D environment (Cranford et al. 2011). Once the target has been detected, the ability to hear the returning echo is equal to their click production mechanism: they can detect all of the sounds they have emitted (Nachtigall and Supin 2008). That also means that the animal hearing differs depending on prey absence/presence, prey size and distance from it (Nachtigall and Supin 2008). In fact, when the target is present the toothed whales increase their hearing sensitivity and can select the larger target (Nachtigall 1980). While hile sound production is active, the hearing ability is a passive mechanism which has been demonstrated that can be regulated by toothed whales (Natchgall andSupin 2013, 2015). This implies that creating a sound must be a requirement to hear the environment and detect the maximum and minimum prey size at minimum energy cost. Unsurprisingly, the evolution of biosonar modes to detect their prey has occurred in parallel with cochlear shape adaptation to deep environments (Park et al. 2019).
Furthermore, having a large body size increases the dive duration through an increase in the amount of oxygen stored in the muscles, and a decrease in the mass specific metabolic rate (Kleiber 1975). Thereby their ability to perform long dives at depth is improved, leading to the evolution of different biosonar types to enhance directional sonar beams for prey echolocation at specific depths (Surlykke et al. 2014). For instance, size and slow click rate in Ziphiidae plays an important role in foraging performance as having a large size increases their prey detection range (Jensen et al. 2018). Toothed whale species produce their diving ecology (deep in red, semipelagic in green, and shallow in blue); c. three biosonar mode groups (broad band in red, frequency modulated in green, and narrow-band-high-frequency); and d. four water surface temperature groups. Letters indicate species abbreviations listed in Fig. 3 and Online Resource 1 different echolocation frequencies and the larger the animal is, the lower the frequency produced (Surlykke et al. 2014). Interestingly this pattern is well established also in bats (Giacomini et al. 2022) suggesting that there might be a biological "rule" within echolocating mammals. PC2 shape vector represented a posteriorly compressed braincase with a more concave profile of the facial region, mainly in Ziphiidae and Kogiidae (Fig. 3). This is a character associated with sound production, directionality of sonar clicks (Galatius and Gol'din 2011;Galatius and Goodall 2016), and pelagic habitats (Cozzi et al. 2016). Along this axis, it is also possible to detect the elevation of the nuchal crest, which in ziphiids is linked to the presence of a very large melon (Bianucci et al. 2016). Changes in this area are associated with the development of premaxillary crests, the general elevation of the vertex, and the increased surface area for attachment of facial muscles, which is associated with movements of the melon to focus the echolocation sound beam (Heyning 1986;Cranford et al. 2008). This elevation may also increase the surface insertion of the muscles on the occipital plate. One of these is m. semispinalis, which originates from the dorsolateral surface of the skull and progresses in a caudal direction to the middle of the thoracic region (Cozzi et al. 2016). This muscle increases the swimming stability, and is usually associated with pelagic and deep-water ecology (Cozzi et al. 2016). However, our results show that at a macroevolutionary scale, changes in skull shape are more related to prey size and peak frequency emitted than to locomotion/diving abilities, which is consistent with a previous study (Bianucci et al. 2016;McCurry et al. 2017b).
Size-related variables also -influence shape variation at the interspecific scale, with brain mass always explaining  Table 4 Procrustes ANOVA and PGLS analyses performed on toothed whale skulls on whole shape in four datasets (60, 56, 31, 26 species) to test covariation between skull size (CS = centroid size), shape, and ecological (EQ = encephalic quotient, kHz = minimum and maximum peak frequency, ST = surface temperature) and metric variables most of the variance in the symmetric shape component, followed by body mass, length and mean and minimum prey size. On other hand, the asymmetric component of shape was impacted by minimum peak frequencies, surface temperature and diving ecology when corrected for phylogeny. Our results show that skull shape correlates with surface temperature and minimum and average (mean) prey size, while general diet categories have no impact on both skull size and shape possibly due to relatively low variation that occurs within the odontocetes clades. Clearly shape is optimised for dealing with the commonest prey items while maximum prey size might be the result of peculiar local adaptations still detectable in the morphology but not on a macroevolutionary scale.
A change in temperature is likely to change abundance, distribution and size of prey, which may induce changes in cranial shape between different temperature regimes (McCurry and Pyenson 2019). These oscillations influenced the evolution of rostral morphology in toothed whales in the Miocene and the Pliocene (McCurry and Pyenson 2019), which was linked to the emergence of different ecological feeding niches (McCurry et al. 2017b). Moreover, many studies hypothesized the evolution of dolphin ecotypes that involve subtle differences in rostral shape due to ecological pressures (Natoli et al. 2004(Natoli et al. , 2006Escorza-Trevino et al. 2005;Adams and Rosel 2006;Morin et al. 2010;Charlton-Robb et al. 2011;Amaral et al. 2012;Mendez et al. 2013;Moura et al. 2013). It indeed makes more sense that mean and minimum prey size exert selective pressure over long evolutionary time scale on skull shape irrespective of phylogenetic history. Toothed whales need to detect both minimum and average prey size with their echosonar system as: i) prey captured of a certain size will meet better their forage daily requirements, but ii) they also need to know the minimum prey size as they do not chew their prey. This can be explained by looking at their skull anatomy. In toothed whales the nasal bones are shifted backwards and this results in a skeleton and anatomical rearrangement of air and food pathways. The larynx shifts to the left and it is encircled by the palatopharyngeal sphincter so that water does not enter in the airways (Reidenberg and Laitman 2018). In fact, toothed whales are left-side breathers and right-side swallowers. This shift allows whole prey to be swallowed, and our findings testify the association with minimum prey size swallowed and skull shape and size.
A possible link between prey size and prey capture mode has previously been hypothesized (MacLeod et al. 2006;McCurry et al. 2017a). Specifically, rostrum elongation divides longirostrine (e.g., Pontoporia) from brevirostrine (e.g., Kogia spp.) species. This is a major feature that relates to capture mode (Werth 2006a). Nonetheless, a nuance is present in feeding adaptations, longirostrine species (i.e., common dolphin) have a long and slender rostrum that generally allows for the rapid capture of prey (ram feeding), and brevirostrine species have a broad, short rostrum and dental reduction (McCurry and Pyenson 2019) usually associated with suction feeding (Werth 2006a, b). Indeed, suction feeding consists of a mechanism that allows a rapid decrease the intra-oral pressure in order to facilitate the entry of the prey (Werth 2006a, b). In addition, the longirostrine species consume a broader prey size range while suctions feeders tend to consume smaller prey. The PCA plot (Fig. 3) shows how selective pressures such as feeding strategy may cause divergence and convergence of skull shape between families. Taxa originating from early diverging nodes in odontocetes evolution, such as the longirostrine Platanista and the brevirostrine Kogia, were quite divergent in morphology from all other odontocetes. This divergence had a significant influence on the identification of large sources of variation in the analysed skull dataset. In contrast, the extant paraphyletic group of 'river dolphins' (Platanista, Inia, Pontoporia, and Lipotes) provide a good example of convergence among disparate toothed whale lineages (Marshall 2009;Page and Cooper 2017). Nonetheless, Geisler et al. (2011), examined living and fossil odontocetes in a phylogenetic context and concluded that many features of the skulls of 'river dolphins' are retained from ancestors and did not necessarily evolve convergently. It is clear from our PCA of skull shape data (Fig. 3) that related toothed whales tend to resemble one another. For example: i) species known to feed on marine mammals (Orcinus orca and Pseudorca crassidens) showed a similar robust skull shape (Vicari et al. 2022a), which is also advantageous for catching and killing large prey (Galatius et al. 2020); ii) many species (e.g. belonging to Globicephalinae and Lissodelphininae) have positions near their closest relatives in the morphospace (Galatius et al. 2020) and skull features are conserved in Lissodelphininae compared to Delphininae, which occupies a larger range of PC1 scores (Galatius and Goodall 2016;Galatius et al. 2020).

Conclusions
This study found significant associations between skull size and brain size, biosonar mode, prey size, maximum peak frequency, and body metric variables across 60 species of odontocetes. Biosonar mode and maximum peak frequencies clearly had an impact on skull size evolution only and not on shape, differently from minimum peak frequencies which impact the asymmetric component of skull shape.
Size is also a strong predictor of skull shape as evidenced by the allometric signal detected in this sample. In this study, most of the ecological relationships with skull morphology were not significant after removal of phylogenetic effects, and only a correlation between skull size and biosonar mode was detected together with the expected association between skull size, body length, brain mass, prey maximum and minimum, and maximum peak frequency. Brain mass, surface temperature, body mass and length, were also important drivers of toothed whales skull 1 3 shape evolution. Restricting the sample to the species for which prey size data were available provided support for an association between both skull size and shape and prey size. This also suggests that hunting specialisation plays a key role in the evolution of skull morphology of odontocetes. This applies to average prey mass that impacts skull shape, while size correlates positively with both maximum and minimum prey size. Larger skulls might be argued to allow production of a stronger bite force necessary to catch and hold potential large prey, as seen in the killer whale (Orcinus orca), and false killer whales (Pseudorca crassidens), and to play an important role in foraging performance. Our finding that, in toothed whales, the shape of the head differs markedly from the shape of the skull demonstrates the importance of considering the orofacial morphology and the shape of the head when examining associations between ecological variables and morphological variation related to feeding.

Acknowledgements
The authors are grateful to museum staff, Christine Lefèvre, Amandine Blin, and Virginie Bouetel at MNHN, Roberto Portela Miguez at NHM, Morten Tange Olsen and Daniel Klingberg Johansson at NHMD, Paolo Agnelli at "La Specola", and Chiara Sorbini at MSNUP for providing access to the collections and infrastructures. We wish to thank Isabelle De Groote and Alessio Veneziano who provided initial guidance on photogrammetry techniques, Carmelo Fruciano and Gabriele Sansalone for their constructive feedback. We would like to thank Morgan Churchill, two anonymous reviewers, the editor Daryl Croft and the associate editor D. Rex Mitchell who provided comments that greatly helped improve this manuscript, and the BioAcoustic Summer School (SeaBASS) at Syracuse University (NY) that provided DV with additional background information on marine mammals acoustics.
Author contributions DV collected the data and performed the statistical analyses. DV together with CM, RPB, MRM, RCS, OL, GB, wrote and revised the manuscript and its analytical interpretations. RCS helped DV during data collection at NHM. CM supervised the project.
Funding This project was funded by Liverpool John Moores University (LJMU) PhD scholarship and LJMU research grant (both awarded to D. Vicari) and SYNTHESYS Project which is financed by the European Community Research Infrastructure Action (FR-TAF-6867, and DK-TAF-6759, both awarded to D. Vicari).
Data availability All data generated or analysed during this study are available at: https:// github. com/ Vicar iD/ Macro-Data. Vicari-et-al. git.

Declarations
Competing interests The authors have no competing interests to declare that are relevant to the content of this article.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.