Is vertebral shape variability in caecilians (Amphibia: Gymnophiona) constrained by forces experienced during burrowing?

Caecilians (Gymnophiona) are a small (about 215 currently recognized species) monophyletic group of elongate, totally limbless and annulated amphibians. Because most caecilians are strongly fossorial, inconspicuous and rarely encountered components of tropical ecosystems, many aspects of their biology remain poorly known (O'Reilly, 2000; Summers and O'Reilly, 1997; Wilkinson, 2012). Although their cranial osteology has been relatively well documented (e.g. Bardua et al., 2019; Lowie et al., 2021; Sherratt et al., 2014; Wake, 1993; Wilkinson and Nussbaum, 1997), few studies have focused upon the postcranial morphology of adult extant caecilians (but see Lowie et al., 2022b; Peter, 1894; Renous and Gasc, 1989; Renous et al., 1993; Taylor, 1977; Wake, 2003, 1980, Wiedersheim, 1879).

In direct association with their burrowing habits, it is thought that caecilian skull evolution resulted in compact and robust skulls with many bones being fused or connected through tight sutures (Nussbaum and Pfrender, 1998; Taylor, 1968; Wake, 1993; Wake and Hanken, 1982). Unexpectedly, however, previous studies on the impact of burrowing force on skull shape suggested that there is no direct relationship between the external forces experienced during burrowing and skull shape (Ducey et al., 1993; Herrel and Measey, 2010; Kleinteich et al., 2012; Lowie et al., 2021). Rather, cranial shape variation appears more constrained by the jaw adductor muscles in relation to feeding (Lowie et al., 2022a). However, the external forces encountered during burrowing are transmitted by the head to the vertebral column and could therefore impact vertebral form and function. Moreover, a recent study by Lowie et al. (2022b) demonstrated that variation in vertebral form occurs along the vertebral column, suggesting that not all vertebrae may be impacted similarly by the constraints imposed by a head-first burrowing lifestyle.

As far as is known, most caecilians are capable of moving through narrow tunnels using a combination of hydrostatic and internal concertina locomotion (Gaymer, 1971; O'Reilly et al., 1997; Herrel and Measey, 2012; Summers and O'Reilly, 1997), thanks to a partial independence between the vertebral column and the skin with the associated external muscular sheath of the body (Naylor and Nussbaum, 1980; Nussbaum and Naylor, 1982). During internal concertina locomotion, the vertebral column is flexed inside the body to provide purchase against the substrate, while the head and anterior-most vertebrae remain largely extended. The waves created by the vertebral column then progressively straighten and the animal pushes its head into the soil (Gans, 1973; Gaymer, 1971; Herrel and Measey, 2012; Summers and O'Reilly, 1997). Whereas we anticipate that the shape of the anterior vertebrae is likely strongly impacted by the reaction forces incurred while pushing the head into the substrate, the mid-body and more posterior vertebrae are likely more impacted by the axial musculature that generates the forces used to push the head forward into the soil.

In this study, we explored whether vertebral shape is associated with the peak push forces measured in vivo for 13 species of caecilian amphibians (Lowie et al., 2021). We predicted that, globally, the higher push forces produced by active burrowers will be associated with more robust vertebrae. We also expected to find a stronger relationship between burrowing force and vertebral shape in the anterior part of the body, as forces are likely dissipated further down the vertebral column. Alternatively, relationships between the external forces and vertebral shape may be more indirect, with the muscles generating the force constraining vertebral shape. In that case, relationships between force and shape would be expected to be strongest for the mid-body vertebrae.

On the basis of vertebral anatomy alone, only a cervical region can be identified unambiguously in all caecilians; post-cervical body regions cannot be precisely delineated solely on the basis of vertebral characters (Wake, 1980, 2003). Caecilians generally have an atlas followed by 67 (Hypogeophis pti; Maddock et al., 2017) to 306 (Oscaecilia cf. bassleri; M.W., personal observation) trunk vertebrae with no sacrum and either a short tail or no tail at all (Dunn, 1942; Maddock et al., 2017; Nussbaum and Wilkinson, 1989). To be able to compare vertebral shape across species that differ in vertebral number, we followed Wake (1980) in selecting six vertebrae: the atlas (V1), the second vertebra (V2), the third vertebra (V3), and the vertebrae at 20% (hereafter referred to as V20%), 60% (hereafter referred to as V60%) and 90% (hereafter referred to as V90%) of the total number of vertebrae. This selection was made under the assumption that all caecilians included in our study have a sufficiently similar vertebral organization. For the atlas (V1), shape was quantified for 53 individuals from 13 species belonging to eight of the 10 currently recognized families. For the five other vertebrae of interest, the dataset consisted of 40 individuals from 13 species based on the availability of whole-body micro-computed tomography (µCT) scans of high resolution (Table 1). 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 sex-specific variation (Sherratt et al., 2014). Specimens were obtained primarily from our personal collections and completed with specimens from museum collections (Table S1).

About half of the µCT scans used were generated 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., 2008, for scanner settings), the Royal Museum of Central Africa (75 kV, 1440 projections) and the Natural History Museum, London (mostly 100 kV, 3142 projections; see Table S1). The isotropic voxel sizes of all scans are listed in Table S1. All the µCT scans were processed using both automatic thresholding and manual segmentation to reconstruct the vertebrae 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 might have interfered with the placement of landmarks.

All anatomical landmarks were placed by the same person (A.L.) using Stratovan Checkpoint (Stratovan corporation, v.2020.10.13.0859). Nineteen homologous landmarks were placed on each atlas, whereas 22 homologous landmarks were placed on the other vertebrae included in our analyses (see Lowie et al., 2022b, for more information). Generalized Procrustes analyses (GPA) were performed on vertebrae using the ‘gpagen’ function in the geomorph R package v 3.3.1 (https://CRAN.R-project.org/package=geomorph). Finally, prior to the analyses, asymmetry was removed from the datasets by extracting the symmetric component of shape variation using the ‘bilat.symmetry’ function of the geomorph package.

In vivo burrowing forces were measured in the lab or the field for 120 specimens belonging to 13 species (see Lowie et al., 2021; Table 2). 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. After measurement, animals captured in the field were released at the exact locations where they were found. Push forces were measured using a custom piezoelectric force platform (Kistler Squirrel force plate, ±0.1 N; Kistler, Zurich, Switzerland) as described previously (Fig. 1; Herrel et al., 2021; Le Guilloux et al., 2020; Lowie et al., 2021; Vanhooydonck et al., 2011). The force plate was mounted on a purpose-built metal base and connected to a Kistler charge amplifier (type 9865). A Perspex block with 1 cm deep holes of different diameter was mounted on the force plate, level with the front edge. One hole with a diameter corresponding to the body diameter of the animal was filled with substrate from the container of the animal being tested. A tunnel with a diameter equal to the maximum body width of the animal was positioned on the metal base in front of (but not touching) the soil-filled hole in the Perspex block. Then, a caecilian was introduced into the tunnel and allowed to move towards the Perspex block. Next, the animal was stimulated to push into the soil-filled hole by gently tapping the end of its body. Forces were recorded in three dimensions using Bioware software (Kistler) during a 60 s recording session at 500 Hz. Each animal was tested at least 3 times, with an interval of at least 30 min between trials. For each trial, we extracted the forces in the X-, Y- and Z-direction of the best push and calculated the highest peak resultant force. We used the highest peak resultant force across all trials for an animal as an estimate of its maximal push force.

Because vertebral data for species are expected to be phylogenetically structured and thus not statistically independent, 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 (Fig. S1) was computed using the ‘maxCladeCred’ function in the phangorn package in R (https://CRAN.R-project.org/package=phangorn).

All the statistical analyses were performed in R version 4.0.3 (http://www.R-project.org/). The significance threshold (Type I error rate) was set at α=0.05. Forces were transformed logarithmically (log₁₀) to fulfil assumptions of normality and homoscedasticity.

To estimate the degree of similarity due to shared ancestry, a multivariate K-statistic (Adams, 2014) was calculated on the mean Procrustes coordinates of each vertebrae and the external measurements using the ‘physignal’ function in the geomorph package. A K-statistic was calculated on the forces using the ‘phylosig’ function in the Phytools package (https://CRAN.R-project.org/package=phytools). 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.

Phylogenetic generalized least squares (PGLS) regressions were performed using the ‘procD.pgls’ function from the geomorph package in order to assess (1) the relationship between the total number of vertebrae and their shape, and (2) the relationship between maximum resultant push force and the shape of each vertebra.

None of the measurements described in this paper (force measurements) 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 push 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 vertebral shape was generally significant but the signal was moderate (V1: K_mult=0.72, P=0.001; V2: K_mult=0.78, P=0.01; V3: K_mult=0.74, P=0.01; V20%: K_mult=0.64, P=0.01; V60%: K_mult=0.69, P=0.01). However, no significant phylogenetic signal was detected for the shape of V90% (K_mult=0.54, P=0.07) or the maximal push force (K=0.66, P=0.53).

Our PGLS regressions showed no significant relationship between vertebral shape and vertebral number (Table 3). However, the PGLS regressions did show a significant relationship between shape and burrowing force for V3, V20% and V60%, but not for the two most anterior vertebrae (V1 and V2) or V90% (Table 4). Overall, for the three vertebrae showing a significant relationship with push force, high forces were associated with vertebral shapes that were antero-posteriorly shorter and taller, combined with a widening and elevation of the posterior aspect of the neural arch. Additionally, high forces were associated with round anterior cotyles whereas low forces were associated with elliptical cotyles (Fig. 2).

Among the species examined, Boulengerula fischeri had the most elongated vertebrae and the lowest forces (Fig. 2). Additionally, the aquatic Typhlonectes compressicauda also showed relatively elongated V60% vertebrae with elliptical cotyles in association with low forces (Fig. 2C). For V3, Ichthyophis kohtaoensis, Siphonops annulatus and Dermophis mexicanus were the species that showed high burrowing forces associated with short and robust vertebrae (Fig. 2A). For V20% and V60%, high forces were mainly associated with short vertebrae and a large neural arch as found respectively in D. mexicanus and S. annulatus (Fig. 2B,C).

One of the assumptions of our study is that burrowing forces generated and experienced by caecilians are relatively high, as suggested by previous studies (O'Reilly et al., 1997). Although until recently, comparative performance data from other head-first burrowing vertebrates were scarce, over the past few years, data on push forces obtained using similar protocols have been published for snakes (Herrel et al., 2021) and skinks (Le Guilloux et al., 2020). When comparing these taxa, caecilians appear to produce, on average, lower forces for a given body diameter than either burrowing snakes or skinks (Fig. 3; Table S2). Caution should be made when interpreting these results as our dataset only comprises one caecilian with a true tail – Rhinatrema bivittatum (Nussbaum, 1977; Nussbaum and Wilkinson, 1989; Wake, 2003) – whereas all the skinks and snakes in our study have tails. During the experiments, some specimens were observed to bend the terminal end of their body to likely gain more grip, potentially impacting forward pushing force. As such, the tail could provide a force advantage in skinks and snakes. However, the only tailed terrestrial caecilian included in our analyses had the lowest force among the caecilians tested (Lowie et al., 2021), suggesting that tails may not be that important in generating forward force. Yet, despite these relatively low forces, our study shows that some of the variation in vertebral shape is correlated with maximal in vivo push forces. As no relationship was observed between the shape of the first two vertebrae and the maximal in vivo forces recorded, this suggests that external forces (reaction forces from the soil, encountered while extending the vertebral column) may not be high enough to significantly impact the shape of the anterior-most vertebrae. Overall, the first two vertebrae in caecilians are rather distinct compared with the other vertebrae (Lowie et al., 2022b; Peter, 1894; Taylor, 1977; Wake, 2003). Although rotation between the two first vertebrae seems limited in caecilians (Wake, 1980), they do appear to be modified for craniocervical movements (e.g. Hoffstetter and Gasc, 1969; Čeranský and Stanley, 2019). For instance, during initial soil penetration, the head angle is important in positioning the head for optimal soil penetration (Kleinteich et al., 2012). The head and the anterior part of the body are also implied in food capture and processing. Apart from rotational feeding (Measey and Herrel, 2006), caecilians also use head movements to tear apart their prey on the surrounding walls of their burrows (Herrel and Measey, 2012). Consequently, other functional demands may be more important in driving the evolution of the shape of the first two vertebrae in caecilians. Interestingly, no link seems to exist between the total number of vertebrae and the shape of any given vertebra, suggesting that species with a different number of vertebrae can safely be compared with one another.

In contrast to the first two vertebrae, the mid-body vertebrae appear to be shaped by constraints incurred during locomotion (Lowie et al., 2022b; Wake, 1980). Our results show that the shape of V3, V20% and V60% is indeed related to the peak push force measured in vivo. However, external forces explain only a relatively small proportion (17–25%) of the variation in vertebral shape, with the shape of the more posterior mid-body vertebrae being more tightly associated with push force. Given (1) the relatively low push forces in caecilians, and (2) the absence of a relationship between push force and the shape of the two first vertebrae, we suggest that the observed relationship between vertebral shape and burrowing force reflects a functional link between axial muscle use during burrowing and vertebral shape, rather than the need for the vertebrae to withstand substrate reaction forces. To test this hypothesis, however, quantitative data on the axial musculature such as the volume or physiological cross-sectional area of the axial muscles obtained through dissection or contrast-enhanced µCT scans are needed.

As suggested previously, variation in stoutness of the body and vertebrae (with stouter vertebrae in species producing higher forces, e.g. D. mexicanus or S. annulatus) might be associated with differences in locomotor type and substrates used (Renous and Gasc, 1989; Renous et al., 1993). Our results show that species generally considered active burrowers such as D. mexicanus generate higher push forces and have short and bulky vertebrae. Moreover, as observed in S. annulatus, the articulatory facets of the zygapophyses are more obliquely oriented, likely increasing the imbrication between adjacent vertebrae, and thus the resistance of the vertebral column to torsion. In contrast, more surface-active species with low burrowing forces appear to have relatively longer vertebrae with more horizontally oriented processes.

Interestingly, the posterior-most vertebra included in our study (V90%) shows no relationship with peak push force. Although little difference in overall shape is observed between V60% and V90% (Lowie et al., 2022b), posterior vertebrae may have additional functional roles. Indeed, the correct positioning of the posterior region of the body is required during copulation (e.g. Kupfer et al., 2006). Moreover, during initial soil penetration, the posterior end of the body is also used to press against the substrate (M.W., personal observation), seemingly to gain purchase as animals push their heads into the soil. Such behaviors likely impose significantly different constraints on mobility and shape. Additionally, in aquatic species such as T. compressicauda, the posterior part of the body may act as a rudder and the increased rigidity provided by longer vertebrae could increase the efficiency of force transfer from the animal to the fluid.

While the vertebrae in the anterior- and posterior-most parts of the vertebral column are likely shaped by multiple functions, and this shows no association with maximal forces measured in vivo, the rest of the vertebrae analyzed do show relationships with maximal push forces. Yet, these relationships are likely due to co-evolution of the axial musculature that attaches to these vertebrae, and the vertebral shape. Further anatomical studies on the axial musculature are needed to fully understand the selective pressures acting on the post-cranial skeleton. Our results further suggest that isolated vertebrae recovered in the fossil record (e.g. Estes and Wake, 1972; Evans and Sigogneau-Russel, 2001; Rage, 1986) may be informative on the burrowing capacity of these species. This information may potentially allow for inference of the lifestyle of fossil caecilians and shed further light on the adaptive evolution of these secretive animals.

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 Museum of Zoology (University of Michigan), the Amphibian & Reptile Diversity Research Centre (University of Texas Arlington), the Royal Museum of Central Africa (Brussels) and 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 M. Segall for her helpful discussions on the statistical analyses, and also M. H. Wake for the gift of Dermophis specimens.