The relationship between head shape, head musculature and bite force in caecilians (Amphibia: Gymnophiona)

Caecilians (Gymnophiona) are a small (>210 currently recognized species) monophyletic group of elongate, totally limbless and annulated amphibians (Pough et al., 1998; Taylor, 1968). Because most caecilians are fossorial, inconspicuous and rarely encountered components of tropical ecosystems, many aspects of their biology are poorly known (O'Reilly, 2000; Summers and O'Reilly, 1997; Wilkinson, 2012). Of the 10 currently recognized families (Wilkinson et al., 2011; Kamei et al., 2012), only Typhlonectidae includes fully aquatic species, species of other families being primarily terrestrial as adults. With a few exceptions, caecilians are head-first burrowers (Taylor, 1968), which imposes specific constraints on the cranial system of limbless vertebrates (O'Reilly, 2000; Wake, 1993). Their typically compact, robust crania, with several composite (fused) bones and tight sutures, have been interpreted as adaptations for head-first burrowing (e.g. Wake, 1993; Wake and Hanken, 1982; Wilkinson, 2012). Although some features of the skull, such as the bony covering of the eye, the position of the mouth and the fenestration of the skull, are probably adaptations to head-first digging, no relationship between skull shape and burrowing force has been demonstrated to date (Ducey et al., 1993; Herrel and Measey, 2010; Kleinteich et al., 2012; Lowie et al., 2021).

Given that the cranium also plays roles in multiple other functions, such as feeding, lung ventilation and housing the brain and major sensory organs (Wake, 1993), it is thus expected to be shaped by many, possibly competing, functional demands. Whereas frogs (Anura) and salamanders (Caudata) show a variety of feeding mechanisms, all terrestrial caecilians studied use jaw prehension, while aquatic caecilians also use suction feeding (O'Reilly, 2000). Although observations of in vivo feeding behavior are scarce, caecilians are thought to be generalist and opportunistic predators, feeding mostly on earthworms and subterranean arthropods (Taylor, 1968; Wake, 1980; Verdade et al., 2000; Delêtre and Measey, 2004; Gaborieau and Measey, 2004; Measey et al., 2004; Kupfer et al., 2005; Herrel and Measey, 2012; Kouete and Blackburn, 2020; but see Moll and Smith, 1967; Presswell et al., 2002; Govindappa and Parwar, 2016, for occurrences of lizards, snakes and other caecilians in the diet of caecilians). The high modularity of the temporal region (Bardua et al., 2019b) and streptostyly observed in caecilian skulls (Kleinteich et al., 2008b; Summers and Wake, 2005; Wake and Hanken, 1982) does, however, suggest a potential relationship between variation in skull shape and feeding mechanisms (Summers and Wake, 2005).

Because of the totally (stegokrotaphic) or partially (zygokrotaphic) closed temporal region and the severe constraints that head-first burrowing puts on the maximal head diameter (Bemis et al., 1983; Gans, 1974; O'Reilly, 2000), the laterally placed external jaw adductor muscles are constrained in size by bones and strongly reduced compared with those of other amphibians (Bemis et al., 1983; Nussbaum, 1977, 1983; O'Reilly, 2000). However, caecilians possess a unique jaw-closing system involving the large and well-developed m. interhyoideus posterior (Bemis et al., 1983; Nussbaum, 1977, 1983). This pennate-fibered muscle is positioned in such a way that its physiological cross-section can be increased without a corresponding increase in head diameter (Nussbaum, 1977, 1983; Bemis et al., 1983; O'Reilly, 2000; Herrel et al., 2019). These morphological features suggest that caecilians can generate considerable bite forces (Herrel et al., 2019; Kleinteich et al., 2008a,b; Measey and Herrel, 2006; Summers and Wake, 2005). Preliminary bite force data presented by Herrel et al. (2019) suggest that the jaw-closing system may be different across caecilians, yet little is known about variation in the jaw muscles and its impact on bite force.

Using 3D geometric morphometrics and quantitative data on cranial muscles [mass, fiber length and physiological cross-sectional area (PCSA)], we investigated the relationships between skull shape, jaw musculature and bite force of 28 species of caecilian amphibians. As shape variation in the cranium and the mandible has already been described by Lowie et al. (2021), only shape covariation was assessed in this study. Whereas the cranium plays a vital role in many activities, the lower jaw is mainly involved in biting and respiration (Carrier and Wake, 1995; Wake, 1993). Consequently, we predicted a stronger covariation between the PCSA of the jaw muscles and mandibular shape than between the PCSA of the jaw muscles and cranial shape. However, we expected the shape of the cranium to covary more with the muscle volume because it reflects architectural constraints imposed by the partial or complete closure of the temporal area. We also investigated the relationships between in vivo bite forces and both skull shape and jaw muscle architecture in caecilians. We predicted that species with a more pronounced and more curved retroarticular process, and thus a greater input lever for the jaw-closing m. interhyoideus posterior, will produce more force (Summers and Wake, 2005). We further predicted that variation in the m. interhyoideus posterior will be the principal driver of variation in bite force as it is the largest jaw-closing muscle. Finally, we predicted that aquatic species using suction feeding will show lower bite forces, yet may possess longer muscle fibers as suction feeding probably puts strong demands on jaw-closing velocity.

We quantified the shape of the cranium and mandible of 81 individuals from 28 species belonging to nine out of the 10 currently recognized families (Table 1), thus capturing a broad diversity in cranial osteology, phylogeny and ecology. Our sample was restricted to adults and included both males and females. Although some sexual dimorphism is present in caecilians (e.g. Kupfer, 2009; Maerker et al., 2016), interspecific variation largely exceeds the sex-specific variation (Sherratt et al., 2014). Specimens were obtained primarily from our personal collections and completed with specimens from museum collections (Tables S1 and S2).

In vivo bite forces were measured in the field for 139 specimens belonging to 11 species (Table 2). Wild-caught specimens were maintained for a maximum of 24 h in large containers filled with substrate collected in the field at sites where the animals were found. Some species (Dermophis mexicanus, Ichthyophis sp., Schistometopum thomense, Geotrypetes seraphini, Herpele squalostoma and Boulengerula taitanus) were maintained under laboratory conditions for longer time periods (months or years), in large tubs with soil and fed with earthworms and crickets. After measurements, field-caught animals were released at the exact locations where they were found. Bite forces were measured using an isometric Kistler force transducer (type 9203, ±0.01 N; Kistler, Zurich, Switzerland) mounted on a purpose-built holder and connected to a Kistler charge amplifier (type 5995A, Kistler; see Herrel et al., 1999, 2001, for a detailed description of the set-up; Fig. 1). The free end of the holder was placed between the jaws of the animal, which immediately resulted in a fierce and prolonged biting. The place of application of the force and gape angle (roughly 20 deg), were standardized. Each animal was tested at least 5 times, with an interval of at least 30 min between trials. The maximal value recorded was considered as the maximal in vivo bite force of the animal, which was then used for subsequent analyses.

For this study, CT scans of different species were used (Table S1). About half of these scans were performed at the Centre for X-Ray Tomography at Ghent University, Belgium (UGCT, www.ugct.ugent.be) using the HECTOR µCT scanner (Masschaele et al., 2013). The scanner settings were sample dependent. The tube voltage varied between 100 and 120 kV and the number of X-ray projections taken over 360 deg was typically about 2000 per scan. Additional µCT scans were obtained from the online repository Morphosource (morphosource.org), the Zoological Museum Hamburg (see Kleinteich et al., 2008a, for scanner settings), the Royal Museum of Central Africa (75 kV, 1440 projections), and the personal collection of M.W. (100 kV, 3142 projections; see Table S1). The isotropic voxel size of all scans is listed in Table S1. All the µCT scans were processed using both automatic thresholding and manual segmentation to reconstruct the cranium and mandible in 3D using Amira 2019.3 (Visage Imaging, San Diego, CA, USA). Using Geomagic Wrap (3D systems), surfaces were prepared by removing highly creased edges and spikes that may interfere with the placement of landmarks. Next, they were decimated to a maximum of approximately 700,000 faces to reduce computational demands without compromising details. For the mandible, only the left hemimandible was used. The ‘mirror’ function in Geomagic was used for the specimens where the left hemimandible was damaged.

We also used 11 specimens from different species for soft tissue visualization. Specimens were stained using either a 2.5% phosphomolybdic acid (PMA) solution or a 6% Lugol's iodine (I₂KI) solution (Descamps et al., 2014; Gignac et al., 2016; Table S2). The staining time varied from 14 to 21 days depending on the size of the specimen. All specimens were then scanned with the HECTOR µCT scanner (100 kV, 2400 projections; Table S2). After a fully manual segmentation in Amira 2019.3, muscles volume was computed using the ‘Material Statistics’ module. The contrast threshold was then manually lowered for each muscle in order to highlight muscles fibers. The fiber length of all the muscles and pennation angle of the MIHP and the MDM were measured using the standard measure tool in Amira. Average fiber length and pennation angle were calculated based on at least 10 fibers per muscle. The PCSA was then calculated by dividing muscle volume by muscle length and multiplied by the cosine of the pennation angle when applicable.

In addition to anatomical landmarks, we also used 3D sliding semi-landmarks on curves (Bardua et al., 2019a; Bookstein, 1991; Buser et al., 2018; Fabre et al., 2018, 2015; Gunz et al., 2005). Both kinds of landmark were placed manually onto the crania and mandibles, by the same person (A.L.) using Stratovan Checkpoint (Stratovan corporation, v. 2020.10.13.0859). Nineteen homologous landmarks and two curves were placed on the left hemimandible and 88 homologous landmarks and four curves on the cranium (see Lowie et al., 2021, for a detailed and illustrated description of the landmarking process). Curves were resampled (see Botton-Divet et al., 2016, for a detailed description of the method) and all the sliding semi-landmarks were slid while minimizing the bending energy using the ‘slider3d’ function from the Morpho R package v. 2.8. Finally, a generalized Procrustes analysis (GPA) was performed using the ‘gpagen’ function from the Geomorph R package v. 3.3.1.

External measurements were collected on all specimens used in the biting trials to characterize their external morphology. Head length from the tip of the snout to the back of the parietal bone, head width at the widest point of the head just anterior to the jaw articulation, head height at the tallest point of the head, lower jaw length measured from the tip of the jaw to the back of the retroarticular process, and body width at mid-body were measured using a digital caliper (Mitutoyo, precision ±0.01 mm). Emaciated specimens were not included. The snout–vent length was measured by stretching the animals along a ruler (±1 mm). Individuals were also weighed using an electronic balance or a spring scale (Ohaus, precision ±0.01 g) (Table 2).

Because species are not independent data points, their phylogeny was taken into account in our comparative analyses (Felsenstein, 1985). The phylogenetic tree of Jetz and Pyron (2018) was pruned to only include the species used in our study. Using 10,000 trees from VertLife.org, the maximum credibility tree was computed using the ‘maxCladeCred’ function in the Phangorn package in R (http://www.R-project.org/).

All the statistical analyses were performed in R v. 4.0.3 (http://www.R-project.org/). The significance threshold was set at α=0.05. External measurements, muscle data and bite forces were transformed logarithmically (log₁₀) to fulfill assumptions of normality and homoscedasticity.

To assess the impact of size on shape, we performed a Procrustes regression on the GPA coordinates using the ‘procD.lm’ function from the Geomorph package. The log₁₀ centroid size was used as a proxy for size. Residuals for both cranium and mandible were then computed and further referred to as allometry-free shapes, in order to examine shape variation not attributable to allometry.

To visualize the evolutionary patterns of shape variation in the cranium and mandible, we performed a principal component analysis (PCA) on the mean of the allometry-corrected shapes for each species using the ‘gm.prcomp’ function from the Geomorph package. Then, we projected the phylogeny onto the morphospace. The results of shape allometry and the principal axes of shape variation of the cranium and the mandible are described in Lowie et al. (2021).

To estimate the degree of similarity due to shared ancestry, a multivariate K-statistic (Adams, 2014) was calculated on the mean Procrustes coordinates of the cranium, the mandible, the muscles (length, volume and PCSA) and the external measurements using the ‘physignal’ function in the Geomorph package. A univariate K-statistic was calculated on the forces using the ‘phylosig’ function in the Phytools package. The phylogenetic signal was calculated under the assumption of Brownian motion (Blomberg et al., 2003). The higher the K-value, the stronger the phylogenetic signal. Values of K>1.0 describe data with a greater phylogenetic signal than expected from Brownian motion alone.

A phylogenetic two-block partial least-squares (2B-PLS) analysis, with the allometry-free cranial shapes in a first block and the allometry-free mandibular shapes in a second block, was performed in order to assess the covariation between the cranium and mandible in caecilians.

Phylogenetic generalized least-squares (PGLS) regressions were performed in order to assess the relationships between external measurements and bite force using the ‘procD.pgls’ function from the Geomorph package. As head width was highly correlated with bite force (see Results; Table 3) and is biomechanically relevant for bite force (e.g. Le Guilloux et al., 2021; Vanhooydonck et al., 2011), residuals of the skull shape (cranium and mandible), bite force and muscles (length, volume and PCSA) were computed using PGLS regressions with head width as co-factor and further referred to as residuals.

Next, the relationships between the residuals of cranial and mandibular shape and the residuals of bite force were assessed using PGLS regressions.

To quantify the covariation between bite force and the jaw-closing muscles (length, volume and PCSA of the MAML, MAMI, MAMA, IHP and MPt), phylogenetic 2B-PLS analyses were performed using the ‘phylo.integration’ function from the Geomorph package, with the residuals of bite force as one block and the residuals of the muscle parameters (length, volume or PCSA) as the second block.

Finally, the covariation between skull shape (cranium and mandible) and the head muscles (length, volume and PCSA of the MAML, MAMI, MAMA, IHP, MPt and MDM) was quantified using 2B-PLS, with the residuals of the shape (cranium or mandible) in one block and the residuals of the muscles (length, volume or PCSA) in the second block.

None of the measurements described in this paper (force measurements or external morphometrics) are considered procedures requiring ethics approval under European law. Furthermore, no permits are needed to maintain caecilians in captivity in Europe. All wild-caught animals were maintained for one night and one day, checked for signs of stress and released at their exact site of capture (marked by GPS) the next night. Captive animals were maintained individually in large tubs of soil in a climate-controlled room (24°C) and fed with earthworms and crickets twice weekly. Animals were checked for signs of stress and injury after bite force measurements and monitored for signs of weight loss during the following weeks. None of the animals were harmed, or showed any signs of stress or weight loss after measurements, and all measurements were approved by local animal care and use committees.

The multivariate K-statistic calculated for the cranium (K_mult=0.83, P=0.001) and mandible (K_mult=0.88, P=0.001) was significant, but the signal was moderate (K<1). However, no significant phylogenetic signal was detected for the muscles (length: K_mult=0.63, P=0.55; volume: K_mult=0.65, P=0.54; PCSA: K_mult=0.69, P=0.44), the external measurements (K_mult=0.62, P=0.57), or bite force (K=0.65, P=0.53).

The phylogenetic 2B-PLS analysis performed on the allometry-free shapes showed a strong and significant covariation between the cranium and the mandible (P=0.001, r-PLS=0.91; Fig. 3). The positive extremes of the scatterplot, mostly driven by the aquatic Atretochoana eiselti, represent species with long and posteriorly projected quadrates that covary with slender mandibles with a proportionately small and dorso-medially curved retroarticular process and relatively posteriorly positioned articulation. The negative part of the plot represents species for which small quadrates covary with mandibles having a longer retroarticular process associated with the shift of the articulation toward the front of the mandible (Fig. 3).

Absolute maximal bite force across the 11 species included in the study ranged from 0.12±0.02 N for Boulengerula fischeri, the smallest species in our dataset, to 11.36±2.79 N for Dermophis mexicanus, one of the biggest species in our dataset (see Table 2). The PGLS regressions showed a significant relationship between maximum bite force and the external head measurements (Table 3). For a given head width, Boulengerula taitanus and Ichthyophis kohtaoensis had the highest bite forces while the aquatic Typhlonectes compressicauda had the lowest bite force, closely followed by B. fischeri (Fig. 4).

Results of the phylogenetic 2B-PLS analysis showed no covariation between residual bite force and residual muscle fiber length (P=0.07, r-PLS=0.76). However, a significant covariation was observed for the residuals of muscle volume (P=0.05, r-PLS=0.74). Species present on the negative part of the plot have bigger MAMI, MIHP, MAML and MAMA and slightly smaller MPt covarying with higher bite forces (Fig. 5A). Muscle PCSA showed the strongest covariation with bite force (P=0.02, r-PLS=0.82). For both volume and PCSA, MAMI and MIHP were the principal drivers of the covariations. Species present on the negative part of the plot have larger PCSA values (except for MPt which negatively covaries with bite force) and higher bite forces (Fig. 5B).

The results of the PGLS analyses showed no correlation between residual skull shape and residual bite force (cranium: P=0.9, R²=0.06; mandible: P=0.8, R²=0.07).

The results of the phylogenetic 2B-PLS analysis showed a significant covariation between the residual cranial shape and the residual muscle fiber length (P=0.017, r-PLS=0.94). The main muscles contributing to the covariation with cranial shape are the MAMI, MPt and MDM. Toward the positive extreme, the shape covariation observed corresponds mainly to a shortening of the os basale and the maxillary bones associated with shorter fibers in the MAMI and MDM but longer MPt fibers (Fig. 6). A slight widening of the braincase, and thus closing of the temporal region, was also observed for species occupying the positive part of the plot (Fig. 6). Additionally, species on the positive extreme have broader interdental plates of the maxillary. Interestingly, no significant covariation was observed between cranial shape and muscle volume (P=0.07, r-PLS=0.90) or muscle PCSA (P=0.27, r-PLS=0.86).

For the mandible, a significant covariation was found between residual mandibular shape and residual muscle volume (P=0.04, r-PLS= 0.90). The main muscles contributing to the covariation with mandibular shape were the MAMI, MAMA and MPt. Toward the positive extreme, the shape covariation observed corresponds to slender and straighter mandibles with short retroarticular processes associated with a large MAMI, MAMA and MPt, but a smaller MIHP (Fig. 7). The negative part shows more medially curved and bulky mandibles with longer retroarticular processes associated with larger MIHPs. No significant co-variation was observed between residual mandible shape and the length (P=0.1, r-PLS=0.88) or PCSA (P=0.19, r-PLS=0.85) of the muscles.

Our results show a strong covariation between cranial and mandibular allometry-free shapes. One extreme is driven by the unusual aquatic caecilian Atretochoana eiselti (see Wilkinson and Nussbaum, 1997). The quadrate bones in A. eiselti are caudally and ventrally elongated, unlike in other caecilians (Wilkinson and Nussbaum, 1997). This quadrate elongation covaries with a posterior displacement of the jaw articulation, leading to a short retroarticular process and a very long pseudodentary. As mentioned by Wilkinson and Nussbaum (1997), this shift of the articulation presumably allows for an increased gape and the capture of relatively large prey. A shorter retroarticular process probably also corresponds to a shift in the use of the adductors as the primary mechanism of jaw closing. This is not unsurprising as aquatic species are not feeding in confined tunnels and thus are freed from the need to reduce the width of their head. This is supported by the relatively low bite forces observed in the aquatic suction feeder Typhlonectes compressicauda. Species with a longer retroarticular process may switch toward the use of the MIHP as the primary mode of jaw closing as the moment arm provided by the retroarticular process becomes greater with an increase in the length of the retroarticular process.

Several studies have highlighted the variability of the quadrate–squamosal complex and have suggested a potential relationship between skull shape and feeding mechanisms (Bardua et al., 2019b; Kleinteich et al., 2008b; Lowie et al., 2021; Summers and Wake, 2005; Wake and Hanken, 1982). Although the main axes of shape variation observed in caecilian skulls are described by changes occurring in the tooth-bearing bones, the temporal region and the quadrate–squamosal complex (Bardua et al., 2019b; Lowie et al., 2021; Sherratt et al., 2014), no relationship between in vivo bite force and cranium shape was found. This suggests that the shape of the skull in caecilians is constrained by factors other than demands for bite force generation.

As observed in previous studies (e.g. Wake, 2003) and confirmed by Lowie et al. (2021), the main axis of shape variation in the mandible is the length of the retroarticular process and the angle it forms with the pseudodentary. Modeling studies focusing on the curvature and length of the retroarticular process, along with the potential movement of the quadrate–squamosal complex (streptostyly), showed that these features have consequences for the biomechanics of jaw closure and can amplify bite force in caecilians (Summers and Wake, 2005). Moreover, as previously shown by Herrel et al. (2019) and hypothesized by O'Reilly (2000) and Summers and Wake (2005), bite force is also highly dependent on gape distance. Although some species appear to be able to maintain bite force across different gape distances, a decrease in bite force was observed with increasing gape in others (Herrel et al., 2019). Species with a dorsoventrally curved retroarticular process have a better mechanical advantage for the MIHP at low gape (O'Reilly, 2000; Summers and Wake, 2005). In contrast, species with a retroarticular process in line with the pseudodentary are expected to perform poorly at low gape. Although the only aquatic species included in our bite force analysis, Typhlonectes compressicauda, did indeed show the lowest bite force for a given head width, no general relationship between in vivo bite force and mandible shape was found. Future studies on bite force at different gape distances are essential to fully understand the jaw-closing mechanisms in caecilians.

To date, no study has assessed the relationship between cranial musculature and skull shape in caecilians. However, this information is crucial to fully understand the evolution of the musculoskeletal system and the significance of morphological traits in driving variation in cranial function. As previously observed, the external adductors are architecturally constrained by size in the cranium. A complete coverage of the temporal region (stegokrotaphy) may limit the expansion of the muscles between the squamosal and the maxillopalatine (Bemis et al., 1983; Nussbaum, 1977, 1983; O'Reilly, 2000). As mentioned previously (Kleinteich et al., 2008b; Nussbaum, 1983; O'Reilly, 2000), the proportions of the different jaw-closing muscles are variable across species, suggesting that different feeding mechanisms may exist. We thus hypothesized that cranial shape would covary with the volume of the different head muscles in caecilians. However, only muscle length covaried with skull shape, suggesting that the organization of the jaw muscles is constrained mostly in length rather than in volume. Muscle length is important and will impact jaw-closing velocity and as such it would be of interest to better understand whether species differ in jaw-closing velocity (Herrel and Measey, 2012). The only study of feeding in these animals in their burrows (Herrel and Measey, 2012) suggested rapid jaw closure, possibly putting a premium on muscle fiber length, especially of the MAMI, which is ideally positioned to generate rapid jaw closure (Summers and Wake, 2005).

Because the large MIHP inserts on the retroarticular process (Bemis et al., 1983; Nussbaum, 1977, 1983), a covariation between mandibular shape and muscle PCSA was expected. However, only the volume of the muscles covaried with mandibular shape, suggesting that rather than the shape of the mandible being constrained by the forces exerted, it is the shape of the mandible that constrains the volume of the muscles that can insert upon it. Interestingly, species with a longer and more robust retroarticular process (e.g. Boulengerula fischeri) have a larger MIHP than species with a small retroarticular process (e.g. Rhinatrema bivittatum), the latter relying more on the jaw adductors for generating bite force. As the MIHP inserts on the retroarticular process, it makes sense that a longer retroarticular process provides a larger insertion area for the MIHP. Moreover, species with a longer retroarticular process appear to invest less in the traditional jaw adductors, which is reflected in a relatively shorter pseudoangular. The covariation observed does indeed suggest a continuum between species that invest in a large MIHP and that have a long retroarticular process versus those that invest in the ‘traditional’ jaw adductors and that have a rather short retroarticular process and a longer pseudoangular and deeper mandibular fossa.

The head musculature in caecilians is very different from that of other amphibians. Apart from their unique and transformed MIHP, caecilians also possess muscles that are not present in other amphibians, i.e. a true MPt and, in some, the m. levator quadrati (Kleinteich and Haas, 2007). Because caecilians also use rotational feeding to reduce prey (Measey and Herrel, 2006), these muscles could also play a role in stabilizing the jaw joint to avoid potential dislocation of the jaw during long-axis body rotation inside tunnels. Clearly, a deeper understanding of the variation in the jaw muscles in caecilians is needed. Moreover, quantitative data on the architecture of the jaw muscles would permit the creation of accurate multibody dynamic models allowing the functional role of the different muscle groups to be tested.

The streptostyly and the unique jaw-closing system in caecilians suggest that these animals can generate relatively high bite forces at low gape distances (Kleinteich et al., 2008b; Summers and Wake, 2005). Large bite forces are usually associated with the consumption of large or hard prey (e.g. Herrel et al., 2002, 2004; Edwards et al., 2013). However, only a few studies reported large prey in the diet of caecilians such as lizards and snakes, among others (Govindappa and Parwar, 2016; Moll and Smith, 1967; Presswell et al., 2002). Moreover, the use of rotational feeding would further suggest no requirement for very large bite forces. Indeed, the interlocking rows of sharp teeth probably maintain the jaws in a closed and locked position when feeding on soft large prey.

However, as suggested by previous studies (Herrel et al., 2019; Measey and Herrel, 2006), our results show that caecilians can generate considerable bite forces at relatively low gape angles for a given head width. Our study shows that higher bite forces are correlated with the overall size of the animal and the size of its head, a pattern also observed in other burrowing tetrapods (e.g. Le Guilloux et al., 2021; Vanhooydonck et al., 2011). Our results further show that bite force covaries with the volume and cross-sectional area of the MAMI and the MIHP. Whereas the important role of the MIHP was expected, that of the MAMI was less so (but see Kleinteich et al., 2008b). However, the MIHP is involved in multiple functions, whereas the adductors are exclusively involved in generating bite force. Indeed, the MIHP is also used in other functions such as buccal oscillation, ventilation and burrowing (Bemis et al., 1983; Nussbaum, 1983; Herrel et al., 2019). The strong covariation of bite force with the adductors suggests that the space available for this muscle may be one of the principal drivers of variation in bite force. Consequently, even small increases in the volume of the muscle may have a disproportionate impact on bite force. As shown by Herrel and Measey (2012), lunges and jaw closing are relatively fast. The longer and parallel-fibered muscles, such as the adductors, may shorten more rapidly and thus increase the velocity of the jaws. In contrast, bipennate muscles composed of shorter fibers such as the MIHP are known to produce more force to the detriment of velocity (Nussbaum, 1983; Summers and Wake, 2005). The combination of the two systems could allow caecilians to bite either forcefully or rapidly, providing them with a versatile jaw system (Summers and Wake, 2005).

Although bite force is correlated with external head measurements, our results suggest that the shape of the skull in caecilians is constrained more by other factors than by demands for bite force generation. However, the strong covariation between cranial and mandibular shapes and their covariation with respective muscle length and muscle volume suggests a continuum between species with long retroarticular processes that invest in large and forceful MIHP versus species with short retroarticular processes that invest in the traditional, but faster, long and parallel-fibered jaw adductors. Moreover, bite force mainly covaries with the volume and the PCSA of the MAMI, suggesting that the space available for this muscle might be driving bite force variation in caecilians. These results highlight the complexity of the jaw system and demonstrate that further studies including other parameters such as jaw velocity or gape distance are needed to fully understand feeding mechanics in caecilians.

We thank I. Josipovic and the people at Centre for X-Ray Tomography at Ghent University for their help with CT scanning. We thank the Natural History Museum (London), Museum of Zoology (University of Michigan), the Amphibian & Reptile Diversity Research Centre (University of Texas Arlington), the Royal Museum of Central Africa (Brussels), the Zoological Museum (Hamburg), A. Kupfer and the Staatliches Museum für Naturkunde Stuttgart and all the curators in these institutions for the loan of some key specimens. We also thank MorphoSource for making scans available. A.L. thanks A.-C. Fabre and M. Zelditch for their helpful discussions on the statistical analyses. A.L. also thanks M. H. Wake for the gift of Dermophis specimens.