Ontogenetic allometry underlies trophic diversity in sea turtles (Chelonioidea)

Despite only comprising seven species, extant sea turtles (Cheloniidae and Dermochelyidae) display great ecological diversity, with most species inhabiting a unique dietary niche as adults. This adult diversity is remarkable given that all species share the same dietary niche as juveniles. These ontogenetic shifts in diet, as well as a dramatic increase in body size, make sea turtles an excellent group to examine how morphological diversity arises by allometric processes and life habit specialisation. Using three-dimensional geometric morphometrics, we characterise ontogenetic allometry in the skulls of all seven species and evaluate variation in the context of phylogenetic history and diet. Among the sample, the olive ridley (Lepidochelys olivacea) has a seemingly average sea turtle skull shape and generalised diet, whereas the green (Chelonia mydas) and hawksbill (Eretmochelys imbricata) show different extremes of snout shape associated with their modes of food gathering (grazing vs. grasping, respectively). Our ontogenetic findings corroborate previous suggestions that the skull of the leatherback (Dermochelys coriacea) is paedomorphic, having similar skull proportions to hatchlings of other sea turtle species and retaining a hatchling-like diet of relatively soft bodied organisms. The flatback sea turtle (Natator depressus) shows a similar but less extreme pattern. By contrast, the loggerhead sea turtle (Caretta caretta) shows a peramorphic signal associated with increased jaw muscle volumes that allow predation on hard shelled prey. The Kemp’s ridley (Lepidochelys kempii) has a peramorphic skull shape compared to its sister species the olive ridley, and a diet that includes harder prey items such as crabs. We suggest that diet may be a significant factor in driving skull shape differences among species. Although the small number of species limits statistical power, differences among skull shape, size, and diet are consistent with the hypothesis that shifts in allometric trajectory facilitated diversification in skull shape as observed in an increasing number of vertebrate groups.


Introduction
Among vertebrates, changing shape as a consequence of changing size (allometry) has been shown to be a primary mechanism for generating morphological diversity (Klingenberg 1998;Erickson et al. 2003;Tokita et al. 2017;Gray et al. 2019;Sherratt et al. 2019). Studies of allometry and ontogeny among animals often focus on the effects of changes in growth and body size on ecology (Urošević et al. 2013;Esquerré et al. 2017;Morris et al. 2019;Gray et al. 2019), as well as their evolutionary and functional consequences (Mitteroecker et al. 2005;Wilson and Sanchez-Villagra 2011;Piras et al. 2011;Bhullar et al. 2016;Esquerré et al. 2017;Morris et al. 2019). Changes in ontogenetic allometry through altered developmental timing (heterochrony) have been demonstrated to be an effective mechanism for dietary adaptations (Denoël et al. 2004;Esquerré et al. 2017;Sherratt et al. 2019), often resulting in differences in skull shape that permit access to new feeding niches (Denoël et al. 2004;Frederich et al. 2008;Morris et al. 2019).
Sea turtles represent an excellent group to study the effects of ontogenetic and evolutionary allometry on dietary habits. They are geographically widespread, monophyletic (Evers et al. 2019), exhibit a range of ecological roles and body sizes (Pritchard and Trebbau 1984;Bjorndal et al. 1997), and have a fossil record with the potential to trace macroecological patterns across deep time (Parham and Pyenson 2010). The seven extant species belong to two families, the monotypic Dermochelyidae and the more speciose Cheloniidae ( Fig. 1) (Naro-Maciel et al. 2008;Duchene et al. 2012). These families likely diverged during the Late Cretaceous (Duchene et al. 2012; Thomson et al. 2021), with the crown of Cheloniidae diverging during the late Oligocene or early Miocene (Thomson et al. 2021). There is dramatic size variation among species (average adult weights between 35 and 400 kg: Pritchard and Trebbau 1984;Dodd 1988;Zug and Parham 1996) and within species, with Dermochelys coriacea (leatherback) increasing in size by three orders of magnitude during growth (Pritchard and Trebbau 1984;Jones et al. 2011). This increase in size is one of the greatest known among extant amniotes, only matched by some of the largest crocodilians (Brandt 1991;Leach et al. 2009). In contrast to other giant testudines such as the giant tortoises of the Galapagos and Aldabra (Chelonoidis nigra, Aldabrachelys gigantea), most sea turtles undergo significant changes in diet during ontogeny (Gibson 1983;Fowler de Neira and Johnson 1985;Furrer et al. 2004). Among them, heterochrony has been previously inferred for one species, D. coriacea, on the basis of postcranial features associated with juveniles found in the adult skeleton, such as cartilaginous epiphyses and a lack of fusion between cervical neural arches and their corresponding centra (Nick 1911;Pritchard and Trebbau 1984).
Despite their low taxonomic diversity, modern sea turtles display a remarkable ecological breadth, with most species inhabiting a unique dietary niche as adults (Bjorndal et al. 1997) (Fig. 1). In contrast, all juvenile sea turtles share a similar diet of plankton and small pelagic cnidarians (Bjorndal et al. 1997;Bolten 2003). Given the functional roles of the vertebrate skull in food acquisition and processing (Pritchard and Trebbau 1984;Claude et al. 2004;Parham and Pyenson 2010;Jones et al. 2012), the relationship between turtle skull structure and diet has been previously investigated in sea turtles but generally with a focus on one or two species. Examinations of skull development of Chelonia mydas (green) and Caretta caretta (loggerhead) indicate that dietary shifts are associated with morphological differences between ontogenetic stages (Nishizawa et al. 2010;Coelho et al. 2018;Lunardon et al. 2020). However, the relationship between skull shape and size among extant turtles remains poorly known, which limits our ability to distinguish ontogenetic and phylogenetic shape differences (Jones et al. 2012) which could inhibit interpretations of turtle evolution. For example, heterochrony has been suggested to be a significant factor in the diversification of other reptile groups, contributing to recent morphological disparity Morris et al. 2019;Gray et al. 2019) as well as early divergences (Esquerré et al. 2017;Morris et al. 2019). Therefore, the understanding of differing allometric patterns can potentially give insights into the evolutionary processes underlying sea turtle diversity.
There have been several morphometric analyses on the skulls of sea turtles. Previous morphometric analyses of sea turtle skulls have either involved linear measurements  Naro-Maciel et al. (2008) and Duchene et al. (2012). Specimens scaled to relative size, based on available literature (Pritchard and Trebbau 1984;Dodd 1988). Top two rows of symbols represent the dietary categories of juveniles (bottom) and adult (top) members of each species: jellyfish = Pelagic; crab = Durophage; sponge = Sponge; sea cucumber = Soft; sea grass = Herbivore; sea cucumber + crab + jellyfish = General. See Table 1 for further details. Silhouettes redrawn from National Aquarium Baltimore(https:// www. aqua. org/ blog/ 2015/ April/ oceans-seven) 1 3 Matsui 1995, 1997;Kamezaki et al. 2003), landmarks (Nishizawa et al. 2010;Coelho et al. 2018), or both (Lunardon et al. 2020). However, most of these studies focused on variation within one species: Ca. caretta (Kamezaki and Matsui 1997;Lunardon et al. 2020) or Ch. mydas (Kamezaki and Matsui 1995;Nishizawa et al. 2010;Coelho et al. 2018). Although some of these studies use large sample sizes they do not include all stages of ontogeny: a study of Ca. caretta used 80 individuals but none were younger than 10 years old (Lunardon et al. 2020), whereas the study of Ch. mydas used 145 individuals all between 3 and 5 years old (Kamezaki and Matsui 1995). Only two studies included multiple species and neither examined growth trajectories or ecological differences. Kamezaki et al. (2003) sampled three species and focused on skull characteristics of Ca. caretta. Myers (2007) included six extant species but found Lepidochelys olivacea (olive ridley) to be notably different from other Cheloniidae, which appears at odds with several qualitative anatomical comparisons (Gaffney 1979;Pritchard and Trebbau 1984;Jones et al. 2012). Other morphometric studies of turtle skull shape have been broader in scope, examining multiple families or Testudinata as a whole (Claude et al. 2004;Foth et al. 2017). As such, the relationship between skull shape, diet, and size among all seven extant sea turtle species remains unknown.
To fill this gap in knowledge, here we use geometric morphometrics to characterise ontogenetic allometry in the skulls of modern sea turtles. We sample all seven species across a wide size range encompassing hatchlings to large adults. While characterising patterns of ontogenetic allometry in each species, we attempt to distinguish among the effects of diet and phylogenetic history on the relationship between skull shape and size. We aim to determine if heterochrony has been a significant mechanism in the evolution of skull shape within Chelonioidea. This study represents the first systematic examination of sea turtle skull shape with three-dimensional landmarks and the most comprehensive study of allometry in the group to date.

Specimens
We sampled 63 specimens from museum collections representing all seven species of extant sea turtle, choosing as broad a size range as possible (Table 1). All hatchlings were ethanol-preserved as well as one large adult specimen of Natator depressus (flatback); all other specimens were dry skulls. Immature sea turtles in their pelagic stage are naturally rarer in museums due to their lower frequency of being washed up on beaches and collected. This factor limited specimen availability, although Caretta caretta and Chelonia mydas have samples representing all stages of ontogeny and total samples exceeding 10 (15 and 11 respectively). The sample of Lepidochelys kempii (Kemp's ridley) lacks hatchlings and large adults whereas Dermochelys coriacea lacks intermediate size animals. Otherwise, all species have at least one specimen for each size category (see below).
Specimens were scanned using X-ray computed tomography (CT). Some of the larger specimens were scanned in medical CT machines, whereas others were scanned using X-ray micro-CT at various facilities (Appendix 1 Table 1). The reconstructed data sets had voxel dimensions between 9 and 500 μm. The resulting tiff stacks were processed in AVIZO 9.0 Lite software (Visualisation Science Group, SAS). The cranium was isolated in an associated label file using the threshold and brush tool. Surface models of the skull were exported as PLY files for measurement (details below). Specimens were divided into three age classes that broadly reflect ontogenetic differences: hatchling, intermediate, and adult (Table 1). Individuals were categorised as hatchlings (n = 14) if they exhibited fontanelles and lacked ossified basicranial elements, individuals were categorised as adults (n = 24) if they had a skull length within two standard deviations of the average reported adult skull size for the species (Pritchard and Trebbau 1984;Dodd 1988;Zangerl et al. 1988;Nishizawa et al. 2010;Lunardon et al. 2020). The remaining individuals were categorised as intermediates (n = 25), as they had closed fontanelles and ossified basicranial elements but were two standard deviations smaller than the skull size reported for an adult of that species.
The following age classes and species combinations were used to determine a total of four datasets: 1. Chelonioidea all age classes: effectively a complete dataset including all seven species and all age classes (n = 63). 2. Cheloniidae all age classes: no Dermochelys coriacea (n = 55). 3. Chelonioidea adults only: only adult specimens (n = 21). 4. Cheloniidae adults only: only adult specimens and no D. coriacea (n = 17).
Dataset 1 provides a holistic overview of the entire sample, dataset 3 serves to examine evolutionary allometry, whereas datasets 2 and 4 provide an understanding of the samples without potential skew from the highly specialised and taxonomically isolated D. coriacea. The four datasets also facilitate comparisons to previous studies that did not include hatchlings or excluded Dermochelys coriacea from some analyses (e.g., Myers 2007). Although D. coriacea is the sister taxon to Cheloniidae, their lineages diverged in the Cretaceous and have been evolutionarily isolated for over 65 Ma (Duchene et al. 2012).
Missing landmarks were estimated using estimate.missing in the package geomorph. Only four specimens of Lepidochelys kempii were missing landmarks, these were from the anterior portion of the snout (landmarks 0, 1, 24, 25) and the tip of the supraoccipital (landmark 5). This function estimates missing landmarks based on the mean coordinates for each landmark across the entire dataset using the regression method of estimation.

Shape analysis
Statistical analyses were performed in the R packages geomorph v. 3.3.3 ) and RRPP v.0.6.2 ). We performed a principal component analysis (PCA) of the Procrustes aligned coordinates, using the function gm.prcomp to determine the main components of shape variation. To interpret the variation described by the major axes we used a combination of thin-plate spline deformation grids, vector analysis, and warping a mesh of the mean shape to the extremes (minimum and maximum PC scores) of each axis. These were implemented using the warpRefMesh function.
To test for a relationship between size and shape we used the phylogenetic generalised least squared method (procD.pgls). This function performs an ANOVA within a phylogenetic framework assuming a Brownian model of evolution (Adams 2014a;Adams and Collyer 2015). However, this technique normally relies on species being either represented by a single set of coordinates or an aggregate mean for each phylogenetic tip (Prevosti et al. 2011;Püschel and Sellers 2016;Wang et al. 2021). Since our data are multiple individuals aligned along ontogenetic trajectories of individual species, implementations currently available need a small modification to be operable. We created a tree file based on the phylogenetic tree found in Duchene et al. (2012) where each species was represented by a soft polytomy consisting of all their individual specimens, such that the phylogenetic relatedness can be considered while retaining the data structure (e.g., Sanger et al. 2013). It should be noted that the validity of phylogenetic comparative methods on groups with a small number of taxa, as is the case in sea turtles, is still unclear (Blomberg et al. 2003;Adams 2014b). Therefore, we also conducted a Procrustes ANOVA which does not account for phylogeny (using procD.lm) for all analyses. In both cases, we evaluated whether a common ontogenetic allometry model or a species unique allometry model better explained the observed patterns for each dataset. This evaluation was achieved by comparing ANOVAs of each model with the null assumption of a common (shared) allometry. The relationship between shape, size, and diet was also assessed using the above methods.
The high number of dietary categories relative to number of species makes it is difficult to disentangle differences in shape associated with diet from those associated with phylogenetic inheritance. However, this problem is somewhat mitigated because dietary overlap does occur between different age classes of some species, e.g. the intermediate age class of Caretta caretta and Lepidochelys kempii exhibit a general diet like the adults of Lepidochelys olivacea (Table 2). Thus, analysing the differences in shape and diet within species may help to inform our understanding of the differences in shape and diet between species.
To assess the strength of evolutionary allometry, the relationship between shape and size of the adult datasets were tested using a PGLS and the original tree (procD.pgls). Species pairwise comparisons were performed using the pairwise function from the package RPPP. This comparison was done to identify differences in slope vector length and orientation, where slope length is magnitude of shape change per unit size and orientation is the direction of shape change per unit size.
To assess whether any of the morphological variation characterised by individual PC axes are significantly correlated with size, we calculated a Pearson's product moment correlation coefficient (cor.test) comparing PC scores to centroid size, assuming a normal distribution. This was done towards understanding which particular aspects of morphological variation are most associated with allometry.
We assessed the morphological disparity between ontogenetic groups using the morphol.disparity function of geomorph. This function estimates disparity of a group (hatchling, intermediate, and adult) as their Procrustes variance using residuals of a linear model fit, in our case using the species unique allometry model.
A residual randomisation procedure with 10,000 iterations was used to assess statistical significance for all tests.
To assess phylogenetic signal, we used the physignal function, which assumes a Brownian motion model of evolution, and the phylogenetic tree recovered in Naro-Maciel et al. (2008) and Duchene et al. (2012) (Fig. 1). This function estimates phylogenetic signal using a generalisation of Blomberg's K statistic for high dimensional multivariate data (K mult ; Adams 2014a). We tested phylogenetic signal for the adults and hatchlings separately, using the mean adult shape and mean hatchling shape for each species, as demonstrated in Gray et al. (2019). As there is a high degree of variability among dated phylogenies of the group (e.g. Duchene et al. 2012;Thomson et al. 2021), we did not use a time calibrated tree. Instead branch length was standardised using the compute.brlen function with the "Grafen" method in the R package ape v. 5.4 (Paradis and Schlip 2019). This comparison was visualised by projecting the phylogenetic tree into morphospace. Phylogenetic signal was assessed to determine how important phylogenetic relatedness is in determining skull shape among living chelonioids, as well as to determine if there is a change in phylogenetic signal across ontogeny which might suggest an increased adaptive signal.

Chelonioidea all age classes (total group: n = 63)
The first four PC axes account for 67.9% of the total cranial shape variation and the remaining axes each account for less than 5%. PC scores vary among species (Figs. 3, 4) reflecting differences in skull shape. Overall PC1, and to a lesser extent PC2, are broadly associated with size and ontogeny whereas the other PC axes appear to reflect species differences and individual variation.
The greatest axis of variation (PC1; 32.3%) represents differences in the relative anteroposterior length of the posterior of the skull, in particular the length of the supraoccipital  PC3 and PC4 of the cranial morphospace of Chelonioidea, with different colours representing different species. Skulls show lateral and dorsal views of the mean shape mesh warped to the coordinates of the extremes for each axis. All points scaled to centroid size crest (Landmarks 5-6), as well as the relative size of the orbit (Fig. 3). High PC1 scores represent relatively large orbits, an anteroposteriorly short posterior part of the skull, and a proportionately short supraoccipital crest whereas low PC1 scores represent relatively small orbits, an anteroposteriorly elongate posterior part of the skull, and a proportionately enlarged supraoccipital crest (Fig. 3). Hatchling skulls have high PC1 scores, juveniles have moderate PC1 scores, and adults tend to have low PC1 scores. However, for Dermochelys coriacea, hatchlings have particularly high PC1 scores, and its adults have similar PC1 scores to the hatchlings of other turtles such as Caretta caretta. The adult specimens of Natator depressus have notably high PC1 scores in comparison to the similarly sized specimens of other cheloniids.
PC2 describes 16.4% of cranial shape variation. High PC2 scores represent a circular orbit, a small ventral projection of the mandibular condyle and a tapering pointed rostrum, whereas low PC2 scores represent an increasing ventral projection of the quadrate and jaw joint, an ovoid orbit, and laterally broader anterior rostrum (Fig. 3). The skulls of larger specimens tend to plot with lower PC2 scores. Skulls with the lowest PC2 scores are exclusively adult D. coriacea.
PC3 describes 10.3% of cranial variation (Fig. 4). High PC3 scores represent skulls that have laterally narrow posterior sections of the skull and have a dorsoventrally shallow lateral profile whereas low PC3 scores represent skulls that have laterally wide posterior sections of the skull and are dorsoventrally deep (Fig. 4). High PC3 scores characterise both Chelonia mydas and Eretmochelys imbricata whereas low PC3 scores characterise the wide cranium of Ca. caretta.
PC4 describes 8.5% of cranial variation. High PC4 scores represent skulls that have an elongated rostrum whereas low PC4 scores represent skulls with a shorter, blunt rostrum ( Fig. 3B). High PC4 scores characterise the long snouted E. imbricata whereas low PC4 scores characterise the blunt snouted Ch. mydas (Fig. 4).

Cheloniidae all age classes (no Dermochelys coriacea: n = 55)
When D. coriacea is excluded from the sample, PC1 (33.7%) is still characterised by the same variation seen in the total group dataset (Fig. 5A). High PC1 scores represent relatively large orbits and small posterior part of the skulls whereas low PC1 scores are characterised by relatively small orbits and large posterior part of the skulls. The characterisation of PC2 (15.7%) is similar to that of PC3 of the total group dataset and describes the relative width of the posterior part of the skull with high PC2 scores having relatively wide skulls and low PC2 scores having relatively narrow skulls.

Chelonioidea adults (n = 21)
PC1 for the total group adult is largely similar to PC1 for the total group, with orbit size and posterior part of the skull size characterising the axis (Fig. 5B). PC2 for the adult subset is similar to PC3 for the total group, characterised by the relative width of the posterior part of the skull.
Cheloniidae adults (no Dermochelys coriacea: n = 17) PC1 (30.3%) for the cheloniid adult subset is similar to PC3 for the total group, characterised by the relative width of the posterior part of the skull. PC2 (19%) is characterised by the shape of the squamosal and jugal (Fig. 5C). High PC2 scores are characterised by a more curved squamosal with the ventroposterior contact with the quadrate being more anteriorly located, and a more anteriorly located contact between the jugal and quadratojugal. Low PC2 scores are characterised by a less curved squamosal with the ventroposterior contact with the quadrate being more posteriorly located, and a more posteriorly located contact between the jugal and quadratojugal.

Chelonioidea all age classes (total group: n = 63)
A significant relationship between shape and size is found using both phylogenetic as well as non-phylogenetic comparative methods (p < 0.001; Table 3; Fig. 6). Both PC1 and PC2 are significantly related to size (both p < 0.001), whereas the other PC axes are not. This result suggests that orbit and posterior skull size both scale allometrically. Species and diet were found to have a significant interaction with size, as well as with each other. For both methods a model using unique species-specific allometries was a significant improvement over a common one (p < 0.001). There was no significant difference in model strength between size and diet or size and species, but both are significantly better than size alone (p <0.001; Table 4). Adults occupy a significantly (p < 0.01) more disparate morphospace than the other two ontogenetic groups. There is not a significant difference in disparity between intermediate and hatchling specimens. Adults occupied a larger area of morphospace (Procrustes variance = 0.0131) than hatchlings (Procrustes variance = 0.008) or intermediate specimens (Procrustes variance = 0.007). This result suggests that the crania of hatchling sea turtles more closely resemble each other than do adult crania.
Most species pairs had significantly different slope angles (Table 5). However, Natator depressus was only different to either D. coriacea or L. olivacea. Lepidochelys kempii is also not significantly different to any other species in slope angle, but this result may be due to the small size range. Caretta caretta and L. kempii have the steepest slopes, whereas L. olivacea and N. depressus having the shallowest ones.
There are few differences in slope length among species (Table 6). N. depressus differs significantly from Ch. mydas, L. olivacea, and E. imbricata. Also D. coriacea significantly differs in slope length from Ca. caretta, Ch. mydas, and E. imbricata (Table 6). L. kempii is not significantly different to any species.

Cheloniidae (n = 55)
Both phylogenetic and non-phylogenetic comparative methods recover a relatively weak but significant relationship between shape and size (p < 0.001; Fig. 7A, B; Table 7). When using phylogenetic comparative methods there was a stronger correlation between size, species, and shape (R 2 = 0.13) than when using a non-phylogenetic method (R 2 = 0.08). There is a significant relationship (Table 8) between shape and diet using both methods (p < 0.001) and both found a significant (p < 0.003) but weak interaction between diet, shape (R 2 = 0.07). Both phylogenetic and non-phylogenetic comparative methods found that a species unique allometry was significantly stronger than a common allometry model (p < 0.001; Table 4). Both methods found that there was no significant difference in model strength between diet and shape or species and shape as possible explanators for shape, but both are significantly better than size alone (p < 0. 001). The Pearson correlation tests found that only PC1 was significantly related to size (p < 0.001).

Chelonioidea adults (n = 21)
Both phylogenetic (p = 0.015) and non-phylogenetic comparative (p < 0.001) methods found a significant relationship between size and shape but neither with particularly strong correlations (R 2 = 0.09 and 0.12 respectively; Table 9; Fig. 7C). However, only the phylogenetic comparative ANOVA found a significant interaction between size, shape, and species (p = 0.005) with a relatively strong correlation (R 2 = 0.42). When using phylogenetic comparative methods a species unique allometry model was significantly stronger than a common allometry model (p = 0.005). Using the non-phylogenetic methods we found that a species unique model was not significantly stronger than a common allometric model (p = 0.167). No individual PC axis was significantly correlated with size.

Cheloniidae adults only (n = 17)
A significant relationship between size and cranial shape was also recovered for both models when D. coriacea was removed from the adult sample ( Fig. 7D; Table 10). Both phylogenetic (p = 0.016) and non-phylogenetic (p < 0.001) comparative methods found a significant relationship between size and shape However, only the phylogenetic comparative ANOVA found a significant interaction between size, shape, and species with a relatively strong correlation (R 2 = 0.47, p = 0.005), with a species unique allometry being the favoured model (p = 0.005). Using the non-phylogenetic methods, we found that a species unique model was not significantly stronger than a common allometric model (p = 0.065). The Pearson correlation test found PC1 and PC2 were both significantly related to size (PC1 p = 0.049; PC2 p = 0.009). No other PC axis was found to have a significant relationship with size.
Phylomorphospace For the mean shapes of hatchlings (Fig. 8A), PC1 represents 47.9% of the shape variation and PC2 represents 20.4% of the shape variation (68.3% for both combined).
For the mean shapes of the adults (Fig. 8B), PC1 represents 50.5% of the shape variation and PC2 represents 21.6% of the shape variation (72.1% for both combined). The phylogenetic Fig. 6 A Multivariate regression of skull shape based on a common allometry against size plotted on log centroid size. B The predicted skull shape based on the regression of skull shape on size of all species based on species allometry. Different colours represent different species, different shapes represent different diets. All points scaled by centroid size signal for the adult specimens was not significant (K = 0.62; p = 0.077), however the signal for the hatchlings was (K = 0.72; p = 0.011). There is also a visible difference between the shape of the phylomorphospace when plotted separately for hatchlings and adults. The hatchlings plot in a manner consistent with the accepted phylogenetic hypothesis for sea turtles (Fig. 1): D. coriacea is relatively isolated in phylomorphospace with a high PC1 score, N. depressus and Ch. mydas plot close to one another with moderate PC1 scores and low PC2 scores, whereas E. imbricata, and Ca. caretta + L. olivacea plot with high PC2 scores and low PC1 scores. The adults plot with less adherence to the phylogenetic tree. Species which are more distantly related (i.e., N. depressus and L. olivacea, and Ch. mydas and E. imbricata) plot with similar scores for PC1 and PC2, whereas the closely related Ca. caretta and L. olivacea, plot relatively far apart. When only the adults are plotted, PC1 represents a similar aspects of variation to PC1  . 7 Allometric regressions for three of the datasets. A shows log size compared to the regression score for shape for Cheloniidae. B shows the predicted shape values compared to size based on a unique species allometry model for Cheloniidae. C shows the log size compared to the regression score for shape for the adults of Chelonioidea. D shows the log size compared to the regression score for shape for the adults of Cheloniidae all points scaled by centroid size   . 8 A Phylomorphospace of mean shape of hatchling specimens. B Phylomorphospace of mean shape of adult specimens. Different colours represent different species. White circles represent shape of hypothetical most recent common ancestor of the clade based on maximum likelihood estimation for the entire dataset (relative size of the orbits and posterior part of the skull) whereas while PC2 is similar to PC3 of the complete dataset (posterior part of the skull width and depth).

Discussion
Our results show that much of the diversity of the skull of modern sea turtles can be explained by variation in size, including both evolutionary and ontogenetic patterns of allometry. We find that hatchling sea turtles have skull shapes that reflect the phylogenetic relationships identified and established by molecular studies (Naro-Maciel et al. 2008;Duchene et al. 2012). From this starting point, differential growth appears to be an important factor in the determination of the adult skull shape. We found that shape and phylogenetic relationships become less obviously linked and diet becomes at least as important a predictor of adult skull shape as phylogeny. This result does not have high statistical power because the number of diet categories and species are similar, but a pattern is evident. In contrast to Myers (2007), we do not find that Lepidochelys olivacea has a highly divergent skull shape from other sea turtles. Instead, our results are consistent with previous studies on individual sea turtle species that infer ontogenetic changes in diet are associated with changes in skull shape (Nishizawa et al. 2010;Lunardon et al. 2020). Variations in ontogenetic changes in skull shape are associated with variation in diet such that multiple instances of paedomorphosis and peramorphosis may be due to divergences in adult diet (Fig. 9). We find that diversity in adult skull shape corresponds to diversity in diet. The relationship appears to be stronger than found in some other groups of amniotes (Maestri et al. 2016;Bright et al. 2016;Gray et al. 2019). The two sea turtle species with the most divergent diets, Chelonia mydas and Eretmochelys imbricata, show relatively clear adaptations associated with their dietary category of herbivore and sponge, respectively. The short rounded snout of Ch. mydas is similar to other aquatic grazing amniotes such as sirenians (Marshall et al. 2012;Aragones et al. 2012) and marine iguanas (Wikeleski and Trillmich 1994). In contrast, the lack of other spongivorous amniotes makes comparisons with E. imbricata difficult. However, its long and narrow snout appears as a likely adaptation to reaching into narrow crevices in coral reefs where sponges often grow (Hill 1998;Figgener et al. 2019), analogous to selective terrestrial browsers which tend to have narrow snouts (Solounias and Moelleken 1993;Dompierre and Churcher 1996). Neither of these morphotypes are associated with size and instead appear to be adaptations to highly specialised diets. The other carnivorous cheloniids have skull shapes characterised by features strongly associated with ontogenetic allometry.
Shape variation accompanying increased body size in sea turtles is most strongly associated with a relative expansion of structures associated with jaw musculature (Jones et al. 2012). The clearest example of this is Caretta caretta, in which allometric skull changes include enlarged temporal regions and jaw muscles that provide greater bite force allowing for diets consisting largely of hard shelled invertebrates (Claude et al. 2004;Huber et al. 2005;Jones et al. 2012;Marshall et al. 2012;Figueirido et al. 2013;Figueirido et al. 2013). This pattern of increased size and development of the jaw muscles and associated structures in comparison to the similarly sized Ch. mydas (Pritchard and Trebbau 1984) is consistent with "acceleration" peramorphosis (Klingenberg 1998), as seen in multiple other vertebrate groups (Denoel et al. 2004;Herrel and O'Reilly 2006;Chemisquy 2015;Vita et al. 2020). Within turtles, this phenotype (with large adductor chambers and deep jaws) 1 3 is also consistently associated with durophagy among unrelated species such as Malayemys subtrijuga and Sternotherus odoratus (Claude et al. 2004;Bever 2009;Parham and Pyenson 2010;Ferreira et al. 2015;Lunardon et al. 2020). The great abundance of extinct cheloniids that have independently acquired durophagous traits, such as those noted in this study as well as robust mandibles with an expanded triturating surface (Gaffney 1979;Parham and Pyenson 2010;Weems and Brown 2017) might be related to a clade-wide ability to pursue peramorphy with relative ease. The two species of Lepidochelys have similar adult sizes but different skull shapes and different diets (Pritchard and Trebbau 1984;Bjorndal et al. 1997). At similar sizes L. kempii has a "more mature" skull shape than L. olivacea, meaning smaller orbits, a larger temporal region, and a longer supraoccipital crest. As in Ca. caretta, these differences are associated with a durophagous diet and may suggest peramorphy. At the same time, it could be equally parsimonious to suggest that L. olivacea is paedomorphic in relation to L. kempii. The adult skull shape of L. olivacea is similar to that of younger individuals of Ca. caretta, which also have a general diet that is less reliant on benthic organisms and better suited to an oceanic rather than coastal lifestyle (Bowen et al. 1997). At present, the small sample of extant taxa and lack of available fossil evidence make it difficult to clearly state which heterochronic patterns explain the shape differences between the two species.
We find that Dermochelys coriacea and Natator depressus have a skull shape consistent with paedomorphosis, characterised by a small supraoccipital crest and relatively large orbits associated with retention of a juvenile-like diet. However, unlike the other species of sea turtles there does not appear to be a strong functional link between skull shape and diet category, except for the large orbits of D. coriacea which may help deep water foraging in low light (Houghton et al. 2008;Horch et al. 2008;Hall 2008;Veileux and Kirk 2014). Both Nick (1912) and Pritchard and Trebbau (1984) found that the postcranial skeleton of D. coriacea also shows evidence of extensive paedomorphosis, suggesting that the paedomorphic skull may be part of a general skeletal pattern rather than the result of selection directly related to diet and the skull itself. Species develop paedomorphic morphologies in multiple evolutionary contexts, including retention of a juvenile life habit (Kordikova 2002), miniaturisation (Rieppel and Crumly 1997;Bright et al. 2016;Esquerré et al. 2017), diet (Denoël et al. 2004;Esquerré et al. 2017;Sherratt et al. 2019), and other complex factors which may be unclear (Bright et al. 2016;Morris et al. 2019;Bardua et al. 2021). Unfortunately, the selective pressures for D. coriacea to develop this phenotype are still unkonwn. The high degree of shell reduction as a consequence of paedomorphosis may assist in achieving a large size given that extensive shell reduction is common among other large sea turtles (Hirayama 1994;Cadena and Parham 2015). There may also be a link between paedomorphy and deep diving, as other deep diving amniotes such the southern elephant seal (Mirounga leonina) share some paedomorphic traits with D. coriacea such as lack of suture closure Pritchard and Trebbau 1984;Goswami et al. 2013). The paedomorphic skull of D. coriacea is likely an example of multiple factors influencing shape including behaviour, habitat, as well as a generally paedomorphic condition rather than just selective pressure on feeding mechanics (Figgener et al. 2019). In N. depressus, paedomorphosis appears to be restricted to the skull (Zangerl et al. 1988), perhaps making paedomorphosis more likely associated with diet in this species. However, this hypothesis requires further investigation.
The oldest (D. coriacea-Cheloniidae) and the youngest phylogenetic divergence (L. kempii-L. olivacea) effectively demonstrate how changes in ontogenetic allometry might accompany speciation, or even act as a mechanism for morphological divergence. Dermochelyids likely diverged from cheloniids in the Late Cretaceous (Duchene et al. 2012;Thomson et al. 2021). Fossil representatives of stem chelonioids from the Cretaceous, such as Toxochelys, do not exhibit obvious evidence of paedomorphy (Zangerl 1980;Gentry 2017;Gentry et al. 2019). The fossil dermochelyid Eosphargis breineri, from the Eocene of northern continents, provides a glimpse into early dermochelyid skull shape (Nielsen 1959). The high degree of similarity between the skull morphology of Eo. breineri (Nielsen 1959;Hirayama 1994) to the modern D. coriacea and the lack of other similar stem chelonioids suggests dermochelyids achieved a paedomorphic skull shape very early in their evolution.
The two species of Lepidochelys may be another example of allometric differences associated with different diets because they are sister taxa and a similar absolute size but have different skull shapes and diets. They are estimated to have diverged 2.5-3.5 Ma during the formation of the isthmus of Panama (Bowen et al. 1997;Naro-Maciel et al. 2008;Duchene et al. 2012): the most recent divergence between species of extant sea turtles (Duchene et al. 2012). As these two species appear to have distinct allometries despite diverging so recently, we suggest that heterochrony was likely a mechanism in their morphological diversification. More fossils representing possible ancestors of the Caretta + Lepidochelys clade are needed to better understand the exact evolutionary processes leading to the diversity we see today. Candidate fossils are available, however they are incomplete and require more extensive study (Zug and Parham 1996).
As we suggest for sea turtles, heterochrony appears to be an important mechanism for phenotypic diversification for much of Reptilia. Within closely related groups, shifts in allometric trajectories have frequently been found to be of major importance in generating adaptive diversity, whether the phenotype in question is body size (Denoel and Joly 2000; Kon and Tetsuo 2002;Sander et al. 2004;McNamara and Long 2012;Esquerré et al. 2017) or shape change (Leiberman et al. 2007;Adams and Nistri 2010;Piras et al. 2011;Morris et al. 2019: Gray et al. 2019. In other reptiles, heterochrony linked to dietary diversity is also seen in Crocodylia where shifts in allometric trajectories (i.e. deceleration, pre-displacement, acceleration) are perhaps the most significant determinants of skull shape, correlating with dietary niches (Morris et al. 2019). Similarly, in pythons, heterochrony was found to be the main mechanism in initial divergence of phenotype, with a strong adaptive link between ontogenetic allometry and ecology (Esquerré et al. 2017). Though the exact genetic mechanisms for various heterochronic patterns are varied depending on taxon and body region (Brugmann et al. 2006;Tokita et al. 2017;Bhullar et al. 2015;Ahi 2016;Morris et al. 2019), and are yet to be established for sea turtles, our study has shown that changes to the timing and extent of development have played a major role in generating sea turtle diversity.
Though our data are highly suggestive of allometry being an important factor in phenotypic diversification in modern sea turtles, the high dietary but low taxonomic diversity as well as small sample size make our results less certain. The almost one to one ratio of dietary specialisation to species is one of the most interesting aspects of sea turtle evolution. However, this condition make establishing any ties between skull shape and diet with any statistical power difficult (Blomberg et al. 2003;Adams 2014b). The challenges of collecting large sample size for some species, in particular L. kempii, also means that our findings are preliminary and require more extensive investigation.

Conclusion
Heterochrony appears to play a significant role in shaping the sea turtle skull. We infer there to be at least two instances of paedomorphy and one instance of peramorphy among extant sea turtles. These instances are associated with differences in diet, suggesting that shifts in dietary niche may be an agent of selection: robust peramorphic skulls with large muscle chambers are associated with hard shelled prey whereas paedomorphic skulls with large orbits are associated with soft bodied prey. However, the paedomorphic skull shape of D. coriacea is part of an overall paedomorphic body plan. Therefore, it may not reflect selective pressures to feed on soft bodied prey per se. Regardless of whether dietary differences are a driver of skull shape change in sea turtles or an incidental consequence, heterochrony is likely to be an important mechanism for facilitating such changes. Our study provides another example among Craniata of the importance of ontogenetic allometry as a mechanism for generating morphological variation.