Differentiation of skull morphology and cranial kinesis in common toads

We examined the cranial morphology and cranial kinesis of the common toads Bufo bufo and B. spinosus with micro-computed tomography and geometric morphometrics and compared the results with published data for related species in a phylogenetic context. The species significantly diverge in skull shape. The skull of B. spinosus is shorter and higher, with a ventral arm of the squamosal bone and the jaw articulation point positioned perpendicular to the braincase, in comparison with a more lateral position in B. bufo. In either species, females have a shorter snout and a higher and wider skull at the jaw articulation point that is positioned more posteriorly, in comparison with conspecific males. High variation in the amount of bone ossification was recorded in both species, ranging from scarcely ossified and loosely connected bones to highly ossified and firmly connected bones. We also found that skull shape and inferred kinetic properties of the skull are highly variable across the Bufonini tribe. However, sample sizes are mostly small and intraspecific variation is high, which might compromise the analyses. Overall, the results suggest that developmental plasticity produces high variation in ossification and cranial kinesis, affecting individuals’ feeding performances. At the population level, this variation supports an efficient exploitation of the habitat and may promote morphological adaptation in a changing environment.


Introduction
The vertebrate skull is an anatomically and developmentally complex skeletal structure with multiple functions. The skull form (size and shape) may be constrained by development and by descent and by a variety of mechanical demands imposed by its functional role in perception, food gathering and vocalisation as well as locomotion and defence . For several groups it has been shown that skull shape reflects selection pressures from functional requirements, such as food intake, in mammals (Galatius et al., 2020;Van Cakenberghe et al., 2002;Wroe & Milne, 2007), birds (Felice et al., 2019), snakes (Andjelković et al., 2016;Klaczko et al., 2016) and frogs (Paluh et al., 2020). Another property of the vertebrae skull that can be directly related with feeding is cranial kinesis. Cranial kinesis denotes the mobility between cranial elements, not including the articulation of the lower jaw (e.g. Kardong, 1977;Herrel et al., 2007). As here applied to anuran amphibians it refers to just the mobility of the upper jaw and the palate relative to the braincase (Iordansky, 1989).
The description of the anuran skull morphology is necessarily general, due to the remarkable diversity and variation in the development of skull skeletal elements (Hanken & Hall, 1984;Herrel et al., 2019;Smirnov, 1990Smirnov, , 1994Trueb, 1973Trueb, , 1985Trueb & Alberch, 1985;Weisbecker & Mitgutsch, 2010). Anurans have a broad, dorsoventrally flattened and fenestrated skull, with a posterior position of the jaw articulation joints. The bones of all three skull complements (neurocranium, dermatocranium and viscerocranium) are largely reduced. The dermal complement consists of the dorsal roofing bones (paired nasals and frontoparietals) and the ventral palatal bones (paired vomers, palates and pterygoids and medial parasphenoid). The neurocranium consist of the sphenethmoid and paired prootics and exoccipitals (Fig. 1). Laterally the ventral arm of squamosal bone reaches the quadrate and its upper arm (the otic ramus) reaches the otic region of the braincase. The quadrate (an element of the viscerocranium) articulates to the lower jaw and together with the pterygoid and squamosal it forms the suspensorium, i.e. the structure that connects the jaw to the braincase. Additional functional support is provided by the palatine that stabilises and connects the upper jaw to the braincase (Roček, 2003;Trueb, 1973;Trueb et al., 1993). The existence of distinct developmental stages (aquatic larvae and metamorphosed juveniles and adults) requires distinct functional adaptations. During metamorphosis an extensive tissue remodelling takes place, including the de novo formation of skull bones, and the full (adult) complement of cranial bones does not form until after metamorphosis is complete (Hanken & Gross, 2005;Kerney et al., 2012). However, the skull morphologies of adults may be affected by diet specialisation of larvae (Herrel et al., 2019, and references therein). The anuran skulls may further diversify in shape and ossification level during juvenile and adult growth (Ponssa & Vera Candioti, 2012). Moreover, the level of ossification of the skull varies from cartilaginous to largely ossified, adding to the disparity in anuran cranial morphology (Paluh et al., 2020).
In terms of cranial kinesis, the anuran skull can be classified as rhynchokinetic and/or pleurokinetic. Rhynchokinesis describes the mobility of the snout, particularly on account of the mobility of the premaxilla. In anurans this mobility allows the closing and opening of the nostrils (Iordansky, 1989). Pleurokinesis describes the mobility of the suspensorium and characterise all lissamphibians (De Villiers, 1938;Iordansky, 2000;Wake & Hanken, 1982). This mobility involves small-scale movements of the upper jaw and palate relative to the braincase and allows for greater flexibility. Conversely, ossification tightens the inter-bone connections and reduces cranial kinesis (Iordansky, 2000;Wake & Hanken, 1982). The skull shape across the family Bufonidae is largely conserved, notwithstanding some major evolutionary changes including ear loss and the reduction of the columella bone (Womack et al., 2018a, b). The genus Bufo is characterised by a well-ossified skull, which may in some species co-ossify with the superficial skin (Duellman & Trueb, 1994;Paluh et al., 2020). Cranial kinesis across the family is less well studied, but the reduction of pleurokinesis due an increase in ossification during postmetamorphic growth has been described only for Rhinella marina (also known as Bufo marinus).
Here, we explore the skull morphology and betweenbone connections of two toad species, the common toad, Bufo bufo (Linnaeus, 1758) and the spined toad, B. spinosus Daudin, 1803. Genus Bufo is positioned within the family Bufonidae (true toads), subfamily Bufoninae and within the Bufonini tribe. The Bufonini tribe holds four highly supported branches including the subtribe Bufonina, which comprise genus Bufo (Dubois, Ohler, & Pyron, 2021).
Bufo spinosus is present in North Africa, the Iberian Peninsula and the southwest of France, whereas B. bufo has a large range extending from France deep into Scandinavia, Russia and Turkey. These species are morphologically and ecologically similar, yet genetically well differentiated, with independent evolutionary histories since the Late Miocene (Recuero et al., 2012;Arntzen et al., 2013a, b). The external morphology of these toads is variable which hampers species identification (Cvetković et al., 2009;Gingras et al., 2013;Lüscher et al., 2001), although they can be told apart by the degree of posterior divergence of the parotoid glands and by the size and shape of the metatarsal tubercle, at least in some regions (Arntzen et al., 2013b(Arntzen et al., , 2018(Arntzen et al., , 2020. Recent studies described some divergences in skull shape between the sexes in B. bufo (Üzüm et al., 2021) and between B. bufo Fig. 1 The skull of Bufo bufo in dorsal, ventral and lateral view. Cranial bones are labelled for the neurocranium (in green), the dermatocranium (in grey) and the suspensorium (in yellow). Note the full set of 12 bilaterally symmetric landmarks that was used for analyses of divergences between B. bufo and B. spinosus species and sexes. A reduced set of six bilaterally symmetric landmarks (marked by solid dots) was used for the phylogeny-based analyses of 50 Bufonini species 1 3 and B. spinosus (Sanna, 2019). There is a large variation in size, especially in B. bufo (Cvetković et al., 2009) with larger body sizes in the south than in the north and with females being larger than males. Despite this marked variation in size, similar diets were recorded, not only for sexes and populations , but also between B. bufo and B. spinosus (Vallvé & Sánchez-Iglesias, 2018), indicating that they are opportunistic feeders.
To explore differences that might exist between B. bufo and B. spinosus, we gathered and analysed data with micro-CT scanning and geometric morphometrics. Micro-CT scanning allows the detailed reconstruction of between-bone interconnections that form the morphological basis to the kinetic properties of the skull (Natchev et al., 2016;Waltenberger et al., 2021). We here provide an in-depth analysis of interspecific differences and sexual dimorphism of B. bufo and B. spinosus in skull size and shape, ossification level and cranial kinesis. To estimate the amount and directions of evolutionary changes in skull between B. bufo and B. spinosus relative to other toad species, we integrate our results on common toads with those for other species in the Bufonini tribe using published three-dimensional models (Womack et al., 2018a, b) under reference to species phylogenetic relationships (Jetz & Pyron, 2018).

Materials and methods
Three-dimensional models of toad skulls were obtained for 27 adult B. bufo individuals (13 females and 14 males) and 18 B. spinosus (nine females and nine males), preserved in ethanol at the Naturalis Biodiversity Center, Leiden, The Netherlands (Supplementary Table S1). Specimens were scanned with the Carl Zeiss Xradia Versa 520 (Carl Zeiss X-Ray Microscopy, Inc., Pleasanton, CA, USA), with 80-kV source voltage and 86-88-mA intensity. The surface 3D model reconstructions of skulls were made with Avizo 9.5 software (FEI, Thermo Fisher Scientific), with some missing data in the case of incomplete skulls (see Supplementary  Table S1). To explore the divergence in skull morphology of common toads relative to other toad species, we included 87 skull models of 48 species from Bufonini tribe available from https:// doi. org/ 10. 5061/ dryad. 6tn2n (Womack et al., 2018a, b). In this sample nine species were represented by a single individual and 39 species were represented by two individuals, all with sexes unknown.
Skull size and shape analyses were based on three-dimensional landmark configurations obtained directly from individual surface models using IDAV Landmark version 3.0 (http:// graph ics. idav. ucdav is. edu/ resea rch/ EvoMo rph). To represent anatomical positions that cover the skull entirely, we applied 12 bilaterally symmetric landmarks amounting to 24 landmarks in total. For comparisons with other toad species, we used a reduced set of six bilaterally symmetrical landmarks, so 12 landmarks in total (Fig. 1). Brief anatomical descriptions of the landmarks are provided in Supplementary Table S2. Skull size was expressed as centroid size (CS) that was calculated as the square root of the sum of squared distances from the centroid. A general Procrustes analysis was performed and the symmetric component was calculated as the averages of the original and mirrored landmark configurations in order to remove the redundancy amongst the bilaterally homologous landmarks (Dryden & Mardia, 1998;Klingenberg et al., 2002;Rohlf & Slice, 1990).
A factorial analysis of variance (ANOVA) was used to test the contribution of independent factors (species and sex) and their interaction (species × sex) to differences in skull size, along with Tukey's HSD post hoc test. A principal component (PC) analysis was performed to explore the differences in skull shape between species and sexes. Because of relatively small sample sizes, the PC scores of nine principal components which describe approximately 85% of total shape variation were used as dependent variables in a multivariate factorial analysis of variance (MANCOVA) instead of the original variables. The MANCOVA was performed with species and sex as independent grouping variables and with skull size as covariate. Statistical tests were done using R version 4.0.2 (R Core Team, 2021). Principal component analyses and skull shape change visualisation were done using the geomorph package version 4.0 . We examined five between bones connections (Fig. 2) that together are taken as a measure of cranial kinesis. We scored character states as either unconnected (value 0) or connected (value 1) and summed these values to obtain an individual's cranial kinesis score (CK i , range 0-5). The Fisher's exact test for count data was performed to test for statistical significance in character state frequency distribution using R version 4.0.2 (R Core Team, 2021).
To explore the amount of ossification and cranial kinesis in a phylogenetic context, we scored the previously described between-bone connections ( Fig. 2) on the models from Womack et al. (2018b) dataset. The between-bone connections were scored binary as described above, with the value 0.5 assigned when both character states were recorded within a species. The sum of scores represents a measure of ossification and cranial kinesis at the species level (CK s ) . The phylogenetic relationships of altogether 50 representatives of the tribe Bufonini (Dubois et al., 2021) with morphological data available were represented by a tree that was constructed by pruning the amphibian Tree of Life (Jetz & Pyron, 2018). For this we used the extraction tool in https:// vertl ife. org/ phylo subse ts/. Five-thousand trees from the pseudo-posterior distribution of trees were downloaded and a maximum clade credibility consensus tree was constructed with TreeAnnotator version 1.10.1 (distributed as part of the BEAST software package; Suchard et al., 2018) and used at all downstream analyses. The species average of symmetric component of shape variation (data for one species 1 3 missing) and CK s were mapped onto the phylogenetic tree. The strength of phylogenetic signal for skull shape was calculated using Blomberg's K statistics (Blomberg et al., 2003) adapted for highly multivariate data (Adams, 2014). The strength of phylogenetic signal for CK s was tested using the phytools package version 1.0-3 (Revell, 2012). For major clades, including the Atelopodina with nine species, Phryniscities with 20 species, Bufonities with 12 species, and Stephopaedities with six species, the amount of between-species variation or morphological disparity (MD) in skull shape and CK s was calculated and tested by pairwise permutation tests. These analyses were done in R using the geomorph version 4.0  and RRPP packages version 1.3.0  for skull shape and with the bootstrap version of Levene's test in the lawstat package version 3.4 for CK s (Gastwirth et al., 2020).

Skull size and shape
A factorial ANOVA showed that the species and sexes significantly differed in skull size with a significant species × sex interaction (F species = 49.37, P < 0.0001; F sex = 56.38, P < 0.0001 and F species × sex = 17.34, P = 0.0002). The sexual dimorphism in skull size was confirmed by post hoc tests, with females possessing larger skulls than males and with a CS ratio of 1.17 in B. bufo (P = 0.0134) and 1.47 in B. spinosus (P < 0.0001). Bufo spinosus females were larger than B. bufo females with a CS ratio = 1.41 (P < 0.0001). No differences were found between males of both species (CS ratio = 1.13, P = 0.1316).
The MANCOVA test showed a significant impact of size, species, and sex on skull shape variation. A non-significant species × sex interaction and significant CS × sex and CS × species interactions were recorded, indicating that species and sexes diverge in the pattern of size-related shape changes (Table 1).
The position of individuals along the first and second PC axes indicated a difference in skull shape amongst species and sexes. Species were best separated along the first axis (that described 26.26% of the total variation) and sexes were best separated along the second axis (20.10%) (Fig. 3). This figure also illustrates the mean skull shape for species and sexes, highlighting that the skull of B. spinosus is higher, with a ventral arm of squamosal bone and the jaw articulation point perpendicular to the braincase compared to more lateral position in B. bufo. The dorsal arm of the squamosal bone (otic ramus) is positioned laterally, more distant from the mid-dorsal line in B. spinosus than in B. bufo. Females of both species have a shorter snout and slightly wider and higher skulls at the jaw articulation point that has a more posterior position compared to males (Fig. 3). Fig. 2 Variation in the degree of ossification of the parasphenoid and prootic in common toads. The arrows indicate the adjacent cranial bones for which character states ('connected' or 'unconnected') were scored as follows: prootic to exoccipital (P-E), prootic to squamosal (P-Sq), sphenethmoid to parasphenoid (Sp-Ps), sphenethmoid to palatine (Sp-Pa) and exoccipital to parasphenoid (E-Ps)

Bone connections and ossification of skull elements
The bone connection characters showed a marked variation ( Fig. 4; Table 2), with no significant differences between species and sexes (P > 0.05 in all comparisons). The cranial kinesis (CK i ) was also highly variable (Table 3) with significant differences for species and sexes (Fisher's exact test, P = 0.035). This test also indicated a significant difference between B. spinosus females and males (P = 0.032) on account of the absence of the between-bone connections in two females (CK i = 0; Table 3).

Morphological differentiation of common toads relative to other species
In the wider set of taxa the four Bufo species (apart from B. bufo and B. spinosus these are B. japonicus and B. bankorensis) had a broad skull and a wide snout (see shapes described by PC1 in Fig. 5) and represented a small amount of the variation documented for the subclade Bufonities or the entire Bufonini tribe (Fig. 5). The lowest level of disparity in skull shape was found in the Atelopodina (MD = 0.004) and the The ossification scores across the Bufonini tribe ranged from CK s = 1.5 (a kinetic skull) to CK s = 5 (a rigid skull). The lowest scores were recorded for Wolterstorffina parvipalmata in the Bufonities and Amietophrynus mauritanicus in the Stephopaedities whereas high scores were found in many species, especially amongst the Atelopodina and the Phryniscities (Fig. 6). A high level of variation in the data was suggested by the observation that amongst 39 species with a sample size of two, ten are polymorphic. The disparity in CK s varied from MD = 1.21 in the Bufonities to MD = 0.48 in the Atelopodina with no significant difference between groups (P > 0.05 in all comparisons). Blomberg's K for the phylogenetic signal in skull shape shows a slight deviation from a Brownian motion model (K = 0.9, P < 0.001) and a larger deviation for ossification scores (K = 0.5, P < 0.01), indicating a tendency for less signal than expected under a Brownian motion model.

The interspecific divergences and sexual dimorphism in skull shape and cranial kinesis in common toads
We observed significant differences in skull size and shape between species and sexes as well as a high level of variation in ossification and cranial kinesis of common toads (B. bufo and B. spinosus). In some individuals the suspensorium and braincase were poorly ossified and loosely connected, whereas in others, they were largely ossified and firmly connected forming rigid skulls (Fig. 4). Studied populations showed the entire range of cranial kinesis ( Fig. 4 and Table 3) with some species-and sex-specific differences. For example, a movable connection between the exoccipital to parasphenoid was only found in B. spinosus females. The intraspecific variation might not have been uncovered if sample sizes had been small, as is not uncommon in present-day micro-CT scanning studies (Natchev et al., 2016), especially when the number of species is large (e.g. Bardua et al., 2021;Fig. 4 The degree of bone connectedness in Bufo bufo skulls illustrated in ventral view. The individuals with different cranial kinesis scores (CK i ) are shown as follows: left CK i = 1, middle CK i = 2 and right CK i = 5. Arrows point to connected bones. The abbreviations are the same as in Fig. 2. Note that ossification increases with increasing CK i values. The three individuals shown are males from England, Leicestershire, Gaddesby Table 2 The between-bone relationships observed in males and females of Bufo bufo and B. spinosus Numbers before and after the slash refer to the connected and unconnected character state, respectively Adjoining bones B. bufo B. spinosus

Females Males Females Males
Sphenthmoid-palatine 12/0 11/3 6/3 6/3 Sphenthmoid-parasphenoid 9/3 5/9 4/5 4/5 Prootic-squamosal 3/9 4/10 1/8 5/4 Prootic-exocciptals 4/8 3/11 1/8 3/6 Exocipital-parasphenoid 12/0 14/0 4/5 9/0  0  0  0  2  0  1  0  3  2  3  2  3  6  3  1  3  5  0  1  0  4  1  3  0  3  5  3  2  1  2 1 3 Paluh et al., 2020). In spite of this variation, skull shapes of species (B. bufo and B. spinosus) and males and females within species were significantly different. Species and sexes differ most in the relative position of the jaw articulation point and the depth of the skull (Fig. 3). Even small differences in skull shape and cranial kinesis can have direct functional consequences, in particular on feeding performance (Hanken & Hall, 1983). An increase in cranial kinesis allows the processing of larger pray, whereas firmly connected, robust constructions provide better muscle support (Herrel et al., 2007). Feeding involves coordinated movements of the cranium, the lower jaw, hyoid and tongue. Toads do not chew their prey and have limited tongue protraction which is directly associated with the lower jaw movements (Deban et al., 2001, Nishikawa & Guns, 1992, 1996. Eye retraction is an important accessory mechanism that assists the primary tongue-based swallowing (Levine et al., 2004). Common toads feed mostly on small-to medium-sized prey including beetles, ants, leaches and worms (Cornish et al., 1995;, but intake may vary depending on pray availability. Occasionally larger and hard prey is devoured, such as molluscs, lizards and mice (Sinsch et al., 2009;Vallvé & Sánchez-Iglesias, 2018). In B. bufo a dietary niche partitioning in pray size was found, with females consuming medium-sized pray in larger quantities compered to predominantly smaller sized pray in males . For food generalists such as common toads, the observed high variation in cranial morphology may help to efficiently exploit their habitat. A wide head and a loose construction of the skull would increase their capability to process large prey, but might also affect feeding performance negatively, in particular bite force, due to a reduced muscle support (Herrel et al., 2007). The increase in cranial kinesis probably does not affect the eyeretraction swallowing mechanism, because the m. retractor bulbi that moves the eyeballs into the buccal cavity is attached to the generally well-ossified, posterior part of parasphenoideum (Witzmann & Werneburg, 2017). 1 3

Phylogenetic aspects of the toad skull morphology
The divergence in skull shape of Bufo bufo and B. spinosus is small in comparison with the shape changes inferred for the Bufonini tribe and the subclade Bufonities (including the genus Bufo) in particular (Figs. 5 and 6). The species also have high intraspecific variation in cranial kinesis. This similarity may be due to the species close phylogenetic relatedness as well as to The ancestral states over the tree were interpolated by a Brownian model of evolution under the likelihood criterion. Tree depth is proportional to time. The two species that we studied are indicated by an arrow their limited ecological differentiation (Arntzen et al., 2013a(Arntzen et al., , b, 2020. To this it may be noted that in ten species of the subclade Bufonities, when just two individuals were studied, variation in cranial kinesis was nevertheless observed. In amphibians, developmental trajectories may vary within species (Roček, 2003) and developmental plasticity plays an important role in the evolution of the cranial skeleton (Trueb & Alberch, 1985;Duellman & Trueb, 1994;Smirnov, 1997;Weisbecke & Mitgutsch, 2010;Bardua et al., 2021). For anurans it has been suggested that diversity in ossification sequences and loss of ossification are unrelated to changes in skull shape (Trueb, 1993). Conversely, recent studies found that ossification sequences (Bardua et al., 2021) and change in the ossification level (Paluh et al., 2020) may both affect evolutionary changes and diversity of the anuran skull. Either way, given the observed intraspecific variation, the role of developmental plasticity has to be taken into account which requires large sample sizes.
In line with other authors (Levis & Pfennig, 2019;West-Eberhard, 2003, 2005, we suggest that the high developmental plasticity which produces a wide range of cranial morphologies may be adaptive and that it would increase evolvability. Testing these hypotheses would involve a biomechanical study of individuals and species with different levels of ossification and connectedness of the skull bones. The pronounced divergences in shape and the high intraspecific and interspecific variation in the kinetic properties of the skull make toads a promising model for the study of functional morphology.

Competing interests
The authors declare no competing interests.
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/.