Experimental modification of morphology reveals the effects of the zygosphene–zygantrum joint on the range of motion of snake vertebrae Free

Trade-offs between the mobility and stability of skeletal joints are determined by many factors including bone geometry, cartilage, and soft tissues such as ligaments and muscles. The shape of joints is highly variable, both across and within species, as well as within individuals, and it depends on both evolutionary history and the function of a joint in the body. For example, canine limbs provide stability during running but are unable to pronate and supinate to the same degree as in felids, reflecting the canine hunting strategy of chasing prey rather than ambushing and grappling prey as cats do (Andersson, 2004). Depending on the joint, soft tissue and bone geometry can have different contributions to the stability and ROM of the joint. Joints with minimal bony constraints, such as the human shoulder, have a high range of motion (ROM) but also rely upon soft tissue for support, resulting in a higher injury risk (Kazár and Relovszky, 1969; Veeger and van der Helm, 2007). By contrast, joints with substantial bony constraints, such as the human hip, have moderate to minimal soft tissue support and a correspondingly lower ROM and injury rate (Anderson et al., 2001; Johnston and Smidt, 1970; Scopp and Moorman, 2001; Zakani et al., 2017).

The vertebral column is essential for protecting the spinal cord and joining the cranial and appendicular parts of the skeleton together. All extant tetrapod vertebrae articulate via the centra and pre- and post-zygapophysis joints that maintain vertebral connection and provide stability (Sumida, 1997). As in other joints, there is a trade-off between stability and flexibility, depending on the needs of the animal (Anderson et al., 2001; Johnston and Smidt, 1970; Kazár and Relovszky, 1969; Scopp and Moorman, 2001; Veeger and van der Helm, 2007; Zakani et al., 2017). Limbless species often rely on the vertebral column for propulsion, and enhanced axial flexibility probably benefits the locomotion of terrestrial limbless animals by allowing them to make contact with and conform to a wide variety of shapes and sizes of surfaces that are used to generate propulsive forces.

Snakes, the most speciose limbless tetrapods, have vertebrae with variable shape, but always possess three articulations: (1) the cotyle–condyle joint at the centrum, (2) the pre- and post-zygapophysis joint, and (3) a distinctive joint formed by the zygosphene and zygantrum (Fig. 1; Hoffstetter and Gasc, 1969; Johnson, 1955; Romer, 1956). This zygosphene–zygantrum joint is large and prominent in all snakes, diagnostic for their vertebrae, and absent or minimal in all other vertebrates (Hoffstetter and Gasc, 1969; Romer, 1956). Furthermore, snakes are the only clade of limbless squamates that have a zygosphene–zygantrum joint, despite many other independent evolutionary origins of limblessness (Wiens et al., 2006). Although the morphology of the zygosphene–zygantrum joint is well described (Auffenberg, 1963; Gasc, 1974; Hoffstetter and Gasc, 1969; Holman, 2000; Johnson, 1955), its function is poorly understood. The zygosphene has long been postulated to prevent axial torsion (Gasc, 1976; Romer, 1956), but empirical tests of the function of this distinctive joint are lacking. However, because all snakes have the zygosphene–zygantrum joint, determining the consequences of its presence or absence is not possible by comparing different snake taxa, and comparisons with other squamate vertebrae lacking a zygosphene would be confounded by other morphological differences.

Methods such as finite element analysis, computed tomography (CT) scanning and 3D printing offer the ability to experimentally alter morphology to test hypothetical alternatives, providing a new way to test the function of existing morphology. For example, caecilians have two types of skulls, one of which was hypothesized to enhance burrowing by reducing bone strain, but no extant species have an intermediate anatomy suitable for testing this hypothesis (Kleinteich et al., 2012). Thus, Kleinteich et al. (2012) digitally altered the skulls of multiple caecilians into the opposite skull type, and finite element analysis during loading revealed no significant difference between the performance of altered and unaltered skulls. This clever use of digital manipulation of morphology solved the problem of isolating the functional consequences of morphological variation not found in extant species. These methods are a way of experimentally manipulating morphology instead of relying only on natural variation and existing biological species to test functional consequences.

In this paper, we provide the first empirical test of the hypothesis that the function of the zygosphene is to prevent roll (Gasc, 1976; Romer, 1956) by experimentally manipulating the morphology of the vertebrae by digitally deleting the zygosphene of four species of snakes. Although we primarily tested the hypothesis that removing the zygosphene will increase roll (torsion), we also examined the broader consequences of the zygosphene–zygantrum joint on ROM.

To digitally reconstruct and experimentally manipulate vertebral morphology, we used specimens of four phylogenetically diverse species of snakes from three families. We dissected and cleaned vertebrae from the mid-body of the following four specimens (Fig. 2): boa constrictor [Boidae, Boa constrictor Linnaeus 1758, American Museum of Natural History AMNH-R176819, snout–vent length (SVL) 54.3 cm], corn snake [Colubridae, Pantherophis guttatus (Linnaeus 1766), AMNH-R176816, SVL 96.7 cm], brown tree snake [Colubridae, Boiga irregularis (Bechstein 1802), B.C.J. personal collection, SVL 184.0 cm] and prairie rattlesnake [Viperidae, Crotalus viridis (Rafinesque 1818), University of Michigan – no specimen number, SVL 83.4 cm].

During prior fieldwork, two authors (H.C.A. and B.C.J.) gathered data on maximal lateral and dorsoventral bending from a freshly euthanized but otherwise intact brown tree snake. We photographed the mid-body while it was maximally bent laterally, dorsally and ventrally, and then we analyzed the images in FIJI (Schindelin et al., 2012; ImageJ 1.8.0_66 64 bit, 3D Viewer, Wayne Rasband, NIH, Bethesda, MD, USA). We calculated intervertebral ROM by determining the total angular displacement over 8–14 vertebrae along an arc with a uniformly minimal radius of curvature and then divided this quantity by the number of intervening joints. For the brown tree snake, these data were obtained from the same snake and same body segments as were used for micro-CT (μCT) scanning and digital rendering (see below). Using procedures similar to those for the brown tree snake, we also determined the ROM for one intact, freshly euthanized specimen each of corn snake (SVL 107 cm) and boa constrictor (SVL 135 cm), although these were different individuals from those used to make the 3D prints of vertebrae. Hereafter we use ‘intact’ to refer to all measurements from these fresh specimens.

For each species, two sequential vertebrae from the mid-body were isolated, dissected and μCT scanned (SkyScan 1172, 80 kV, 120 μA, Bruker, Billerica, MA, USA), segmented manually with Adobe Photoshop (CC 2015 Adobe Inc., San Jose, CA, USA) and digitally rendered with FIJI (Fig. 1). Voxel size was 26.16 μm for all snakes except B. irregularis, for which it was 19.88 μm. After segmentation, the vertebrae were smoothed in MeshLab (v.2016.12 ISTI-CNR, Pisa, Italy) with a Laplacian smooth method (coefficient of three). Two copies were made of the posterior vertebra for each species, one of which was unaltered and the other which was edited digitally to remove the zygosphene (Fig. 1E,F) using Meshmixer (v.3.4.35 Autodesk, Inc., San Rafael, CA, USA), referred to as the altered vertebra henceforth. All vertebrae were then 3D printed [Lulzbot TAZ 6, Fargo Additive Manufacturing Equipment 3D, Fargo, ND, USA; layer height 0.18 mm, xy resolution 0.05 mm, acrylonitrile butadiene styrene (ABS)] at 14× their actual size (7–16 mm) to limit the effects of print resolution and the expanding plastic as the vertebrae are printed.

To determine ROM during manual manipulations of the 3D-printed vertebrae (Movie 1), we used four motion capture cameras (Flex 13, 120 images s⁻¹, NaturalPoint, Inc., Corvallis, OR, USA). We attached six or more adhesive reflective markers to the anterior vertebra and affixed four reflective spherical markers to the posterior vertebra (Movie 1). The anterior vertebra was fixed in place to prevent movement. We tracked the markers using Motive Optitrack v.2.0.2 (NaturalPoint, Inc.), and the markers on each vertebra were used to define a rigid body. The software uses a non-linear least-squares solver to reconstruct rigid body position, and computes rotations using quaternions, which are then decomposed into yaw, pitch and roll relative to the external frame of reference. Because the computation is initially in quaternions, the order of rotations that can affect Euler angles is not problematic for this methodology (Richards, 2019; Richards and Porro, 2018). Dorsal pitch is defined as dorsal motion of the vertebra due to rotation about a horizontal transverse axis (Fig. 3E), whereas ventral pitch is ventral motion of the vertebra due to rotation about the same axis (Fig. 3F). Yaw is defined as lateral motion of the vertebra due to rotation about a dorso-ventral axis, while roll is defined as rotation of the vertebrae about an antero-posterior axis. In all cases, the center of rotation was considered to be the cotyle–condyle joint. We calibrated cameras via a wanding procedure in the Motive Optitrack software and advanced to trials only if the calibration error was below 0.1 mm. Vertebrae were manually manipulated over the ROM for a minimum of 11,000 frames of data to ensure thorough coverage of the possible ROM, including intermediate values. When manipulating the vertebrae, the cotyle–condyle joint remained articulated (Fig. 3), as disarticulation may damage living snakes. In order to achieve full coverage in kinematic space, the graphed data were analyzed for gaps, and subsequent trials were performed to achieve full coverage. Certain postures at which roll was effectively unlimited would likely damage a living snake; hence, we excluded all data when roll was beyond 40 deg. Because we used natural vertebrae, we note that morphological asymmetries could be present and produce asymmetries in ROM (Fig. S1F,G), though none of these asymmetries appeared to affect yaw at zero pitch. The yaw, pitch and roll angles were not smoothed or filtered.

We used a custom-written script in MATLAB (MathWorks, Natick, MA, USA) to analyze data (Script 1). Because of the imprecision of manually centering the vertebrae, we calculated the 1st and 99th percentile values of yaw at zero pitch, and assumed that these values would be symmetrical about the neutral axis of the vertebrae in yaw. Thus, we subtracted half of the difference between these percentile values from all yaw values to center the data. Because our data were not smoothed, we used percentile values to ensure that isolated high values due to error did not unduly bias our data. We analyzed a subset of data from the fixed vertebrae to determine error for yaw, pitch and roll in all four species. The errors from the fixed vertebrae (n=11,000 data points minimum per fixed vertebra) for yaw, pitch and roll ranged from 0.6 to 1.2 deg, 1.1 to 1.9 deg and 0.4 to 1.4 deg, respectively, showing high precision and accuracy.

To determine the anatomically neutral pitch (true zero), we constructed a line between the centers of the circles of curvature of the cotyle and condyle in a transverse μCT scan (Fig. 4A). After determining this line independently for each vertebra, we locked custom, 3D-printed parts (red component, Fig. 4A) onto the condyle and projected the reference line outward (Fig. 4A). This gave us the orientation in ‘lab space’ of the true zero for pitch of the anterior (fixed) vertebra. The corresponding part was attached to the posterior (mobile) vertebra (Fig. 4C), which was then raised to maximum dorsal pitch. The resultant difference between the slopes of the true zero axes of each vertebra in world space was used to determine true maximum pitch. We used this value to re-zero the pitch of the datasets measured from Motive Optitrack. Roll was oriented in lab space by making a line with the pre- and post-zygapophyses and making this line parallel to the table line (Fig. 4B,E). Yaw was oriented by making the neural spines line up in dorsal view (Fig. 4D,E).

The hypapophyses of both the brown tree snake and prairie rattlesnake created a bony limit on ventral pitch. However, the boa constrictor and corn snake lacked hypapophyses that limited ventral pitch and thus prevented disarticulation of the cotyle–condyle. Because this is unrealistic, we used the intact dorsal and ventral pitch data from the boa constrictor and corn snake measurements to define their ventral pitch limits. Once we obtained the maximum isolated value of dorsal pitch based on the printed vertebrae, we removed points from the dataset that would have exceeded the dorsoventral pitch ROM from the intact snakes. Reported maximum and minimum yaw angles for the altered isolated vertebrae were restricted to values within the range of roll reported from normal vertebrae to prevent high roll angles influencing yaw. Values reported for ROM refer to only one side.

To quantitatively characterize and compare ROM between species (between both normal vertebrae and altered vertebrae) and between normal and altered vertebrae of the same species, we created a custom-written Matlab script to calculate overlap of the areas in 2D (yaw and pitch) kinematic space (supplementary information). First we used the ‘boundary’ function with a shrink factor of 0.8 to create a boundary that conforms to the shape of the point cloud of pitch and yaw values. These boundaries were then converted to polygons using the ‘polyshape’ function. Finally, the ‘intersect’ function was used to determine overlap. If two separate regions of overlap occurred, the data were analyzed using the above steps in two parts – one for positive yaw and one for negative yaw – and the resulting areas were then added together. Using these tools, we conducted two separate tests. First, to assess whether values of increased roll in altered vertebrae occurred at pitch and yaw combinations within the pitch and yaw ROM of normal vertebrae, we compared the overlap between the pitch and yaw ROM of normal vertebrae and the pitch and yaw values of altered vertebrae at which roll was >2.5 deg (Fig. S2). Second, in order to assess shape similarity of ROM across species, the overlap between the pitch and yaw ROM was computed for all interspecific pairs of both normal and altered vertebrae after being normalized by area for each combination of snakes (Figs S3, S4). In altered vertebrae, regions with high roll were eliminated (Figs S3, S4). To avoid being confounded by differences in the total ROM between species, we normalized the polygon areas of the ROM of the species being compared by equalizing their areas to each other. Thus, two highly overlapping shapes show high similarity while shapes that are highly different will show low overlap.

Despite the variable shapes and proportions of the vertebrae, some general patterns of yaw, pitch and roll were consistent for all four species (Table 1, Figs 5–8; Fig. S1, Movie 2). The normal isolated vertebrae from all species had a range of maximal yaw values from 13.9 to 18.5 deg, ventral pitch values from −13.5 to −8.7 deg, dorsal pitch values from 4.8 to 10.7 deg, and roll values <±2.5 deg for all species at all combinations of yaw and pitch (Table 1, Figs 5–8).

After normalizing the yaw–pitch ROM areas of the normal isolated vertebrae to better analyze ROM shape space, the values of yaw–pitch ROM overlap between pairs of species ranged from 56% to 89% (Table S1). The prairie rattlesnake consistently had the lowest overlap with other species (56–69%) as a result of asymmetries in the vertebra, whereas the other three species had quite similar overlap values with each other, between 82% and 89% (Table S1).

The altered isolated vertebrae had increased ROM for yaw, pitch and roll, with most yaw and pitch combinations showing similar values of roll to unaltered vertebrae (Table 1, Figs 5–8). Yaw ROM depended on pitch position and vice versa for normal and altered vertebrae. The normalized areas of the altered vertebrae all had very similar overlap values ranging between 79% and 89% (Table S2). The prairie rattlesnake overlap percentage increased to be similar to that of the other three snakes analyzed in the altered vertebrae compared with the normal vertebrae (Tables S1, S2).

The maximal bending in yaw of the intact vertebrae for all species agreed closely with values obtained from the normal isolated vertebrae with a difference ranging between 0.7 and 4.1 deg and an average difference of 2.4 deg (Table 1, Figs 5–8). Yaw increased for all species in the altered vertebrae (Table 1, Figs 5–8). Yaw values for normal isolated vertebrae ranged from ±13.9 to ±18.5 deg but increased to ±18.8 to ±26.9 deg for the altered vertebrae. For the altered vertebrae, the smallest increase was 27% for the prairie rattlesnake, and the largest increase was 42% for the brown tree snake. The corn snake and boa constrictor had increases of 31% and 38%, respectively.

The maximal values of bending in ventral pitch for the intact vertebrae of the brown tree snake agreed closely with values obtained from the normal isolated vertebrae with a difference of only 2.1 deg (Table 1, Fig. 5). Compared with ventral pitch of the normal isolated vertebrae, that of the altered vertebrae increased for both the brown tree snake and prairie rattlesnake (Table 1, Figs 5–8). We could not determine maximal ventral pitch angles for 3D-printed vertebrae of the boa constrictor and corn snake because they pitched ventrally indefinitely. In the isolated vertebrae of the brown tree snake and the prairie rattlesnake, ventral pitch was limited by the zygosphene–zygantrum articulation. Maximal values of ventral pitch for the normal isolated vertebrae ranged from −13.5 deg to −8.7 deg, and they increased to −14.7 deg and −11.6 deg in the altered vertebrae of the brown tree snake and prairie rattlesnake, respectively. When the zygosphene was removed, ventral pitch increased for the brown tree snake by 10% and in the prairie rattlesnake by 36%, and ventral pitch was further limited by the hypapophyses on the anterior vertebra contacting the rim of the cotyle on the posterior vertebra (Figs 3F, 5 and 6).

The maximal values of bending in dorsal pitch for the intact vertebrae agreed closely with values obtained from the normal isolated vertebrae with a difference in range between 0.6 and 2.4 deg and an average difference of 1.6 deg (Table 1, Figs 5–8). Altering the isolated vertebrae increased dorsal pitch for all species in this study except for the nearly constant values for the boa constrictor (Table 1, Figs 5–8). Among all of the species, maximal values of dorsal pitch for normal isolated vertebrae ranged from 4.8 to 10.7 deg and increased in the altered vertebrae from 7.7 to 15.9 deg. Maximal dorsal pitch increased for the brown tree snake by 67%, for the prairie rattlesnake by 58% and for the corn snake by 49%.

Roll increased for all species in the altered isolated vertebrae but by different amounts (Table 1, Figs 5–8). Roll increased from 1.7 deg to 25.3 deg in the brown tree snake and from 1.7 deg to 12.5 deg for the prairie rattlesnake. Roll for the corn snake and the boa constrictor was effectively unconstrained in the altered vertebrae, with values greater than the 40 deg cut off. Roll was limited in the brown tree snake and the prairie rattlesnake because roll caused the cotyle–condyle to twist and disarticulate, which is not possible in the corn snake and boa constrictor. Thus, we halted roll at values before disarticulation would occur.

Discernible roll only occurred in the altered vertebrae at combinations of high yaw and ventral pitch for all species. Roll ROM quickly shifted from <2.5 deg for a particular pitch–yaw combination to >30 deg over small changes in pitch–yaw position in the altered vertebra (Figs 5–8). Such increased roll only occurred in the altered vertebrae at combinations of yaw and pitch of approximately ±15 to ±20 deg and −10 deg, respectively. These extreme postures did not overlap with the normal isolated vertebrae ROM in the brown tree snake and the prairie rattlesnake. However, there was overlap of regions of increased roll in the altered vertebrae that occurred within the normal vertebrae ROM of 3% and 11% in the boa constrictor and the corn snake, respectively. The pre- and post-zygapophyses prevented roll at lower yaw angles, especially in the prairie rattlesnake, which has large pre- and post-zygapophyses compared with the other species (Fig. 2).

Snake vertebrae have highly variable shapes, but a prominent zygosphene is always present (Gasc, 1974; Johnson, 1955; Romer, 1956). Our experimental removal of the zygosphene increased roll in all snakes examined, but only at a combination of low ventral pitch and high yaw that could not occur when the zygosphene was present in two species, as shown by our overlap data (Figs S3, S4). ROM was generally similar between species with and without a zygosphene (Figs S3, S4). The overlap of the pre- and post-zygapophyses prevented roll at the other positions (Fig. 3). The base of the hypapophyses in the prairie rattlesnake and the brown tree snake limited ventral pitch because they collided with the cotyle of the adjacent vertebra. The zygosphene in the normal vertebrae and the neural arch or the pre- and post-zygapophyses in the altered vertebrae limited dorsal pitch. By contrast, the neural spines did not contact each other and therefore did not limit dorsal flexion on the 3D-printed models.

The ability of bony structures other than the zygosphene–zygantrum to constrain roll raises the question of why the novel zygosphene–zygantrum articulation evolved, rather than enlargement of existing processes. Perhaps the zygosphene acts as an osseous limit on the ROM of the vertebral column, particularly in yaw. An osseous limit is stronger than relying on soft tissue for restricting movement because of the ability to withstand higher forces without damage. Snake vertebrae can have a wide ROM during some normal behaviors (Jayne, 1988; Morinaga and Bergmann, 2019; Sharpe et al., 2015), and reliance on soft tissue limits could increase the propensity for injury as seen in rotator cuff injuries and their prevalence among humans (Yamamoto et al., 2010). This bony limit could be beneficial for multiple behaviors that rely upon both high forces and tight bending such as constriction, concertina locomotion or gripping narrow arboreal substrates. The zygosphene could also reduce the need to rely on proprioception and motor control to prevent the body from reaching postures that could damage tissues similar to the antitrochanter articulating with the femur to restrict rotation in birds (Kambic et al., 2017). Thus, a bony blocker has some benefits that expanding the zygapophyses would be unable to provide (limiting yaw), potentially outweighing the reduced ROM due to the zygosphene. Additionally, snakes can still achieve high flexibility despite the limitations of the ROM via their large total number of vertebrae compared with other vertebrates. Together, these benefits could potentially give snakes access to increased range of motion compared with legless lizards, which lack a zygosphene–zygantrum articulation (Hoffstetter and Gasc, 1969). The lack of a zygosphene in other limbless tetrapods is puzzling, and may be due to evolutionary trade-offs, constraints or contingency. Presumably, limbless lizards must limit ROM using either soft tissue structures, which are vulnerable to injury, or modifications of existing vertebral processes, which may reduce ROM (Anderson et al., 2001; Johnston and Smidt, 1970; Kazár and Relovszky, 1969; Scopp and Moorman, 2001; Veeger and van der Helm, 2007; Zakani et al., 2017). More research is needed to determine how limbless lizards limit roll and measure how their overall ROM compares with that of snakes.

Our results suggest that normal isolated snake vertebrae can roll 2 to 3 deg even with the zygosphene present. However, we suggest caution, as such small values could be due to cumulative errors from limited CT resolution, errors in segmenting and the 3D printing process. The error from the fixed vertebra showed errors in pitch, yaw and roll of between 0.4 and 1.9 deg. Any small errors due to our process will comprise a small fraction of the total ROM of pitch and yaw as a result of the large values of these variables (up to 20 deg), while small or zero values could be dominated by error. Therefore, roll in vivo could be minimal and not biologically important. A prior study (Moon, 1999) reported similar amounts of roll in snakes in vivo and in skinned body segments and suggested that roll is actively used in locomotion. However, these in vivo estimates of roll were determined from 2D dorsal-view distances between marks along the dorsal midline and the sides of the snake, although the author acknowledged that rib movements can alter this distance (Moon, 1999). Furthermore, apparent roll across multiple joints (with no actual roll between adjacent vertebrae) can arise from combined lateral and dorso-ventral flexion without departing from planar motion (Zhen et al., 2015), which could have confounded both the in vivo and skinned body segments data in Moon (1999). Without more direct methods of obtaining in vivo data on vertebral motion, it remains unclear whether substantial roll occurs in vivo in snake vertebrae.

Our ROM values of yaw are similar to some previously reported maximal values for snakes (Jayne, 1988; Sharpe et al., 2015), suggesting they are biologically relevant. Even though one might expect diverse species of snakes to have large differences in ROM, many values are actually rather similar. The maximal amount of realized yaw, however, does vary among the different modes of snake locomotion. For example, Jayne (1988) estimated yaw from marks on the mid-dorsal scales of Crotalus cerastes, Nerodia fasciata and Pantherophis obsoleta and found less yaw during lateral undulation (approximately 5 deg) than during sidewinding (approximately 7–10 deg). After commonly finding values of yaw between 15 and 16 deg in both N. fasciata and P. obsoleta performing concertina locomotion in tunnels ranging from 3% to 10% of total snake length, Jayne (1988) suggested that snakes using this mode often may approach their maximal yaw ROM. Although our study did not include the same species, our maximal yaw for a congeneric (Pantherophis) reached a similar value of 16 deg (Table 1), supporting the idea that at least some snakes do indeed use their maximal yaw ROM during concertina locomotion. Morinaga and Bergmann (2019) found Nerodia sipedon had inter-vertebral joint angles between 6 and 9 deg when using lateral undulation between different peg spacings. These values are higher than Jayne (1988) found for lateral undulation, but they are still below any of our observed values, supporting Jayne's (1988) conclusions that during lateral undulation, and many other locomotor activities, snakes do not bend maximally. Additionally, Sharpe et al. (2015) found that anesthetized Chionactis occipitalis (165–175 vertebrae) could form an average of 6.2 complete (360 deg) coils when bent by an experimenter (though whether these are maximal values remains unclear), resulting in average yaw ROM between 12.8 and 13.5 deg, which is similar to Jayne (1988) and our own results (Table 1). Thus, despite considerable phylogenetic and ecological differences among the few species studied to date, ROM appears broadly similar across snake species, suggesting that variation in overall flexibility depends primarily on variation in the total number of vertebrae.

Snake vertebral ROM values for yaw, dorsal pitch and ventral pitch are also similar to values previously reported for crocodiles (lateral, dorsal and ventral flexion: 18.6, 9.4 and 12.7 deg, respectively) (Molnar et al., 2015), skinks (lateral flexion: 13.5 deg) (Sharpe et al., 2015) and armadillos (lateral, dorsal and ventral flexion: 6, 7 and 8 deg, respectively) (Oliver et al., 2016). Collectively, the data suggest that the considerable flexibility of snakes is mainly from large numbers of vertebrae per unit rather than unusually high ROM.

Digitally altering bones to create hypothetical morphologies is a powerful tool to circumvent limitations of using naturally occurring biological morphologies (Kleinteich et al., 2012; this study). In our study system, differences in vertebral form between lizards and snakes would create a confounding factor that could not be excluded via traditional analyses. Consequently, creating hypothetical morphologies provides a way to test effects of morphology that would otherwise be impossible using only natural specimens. This approach can also be applied to paleontological specimens (Shiino et al., 2012), modeling of transitional forms throughout the fossil record, and for small bones for which testing can be difficult (this study). However, caution is warranted in interpreting these results, because these are not natural morphologies, and may not represent the ancestral form of snake vertebrae prior to the development of the zygosphene. Potential sources of discrepancies include rib interactions, partial cotyle–condyle disarticulation, and soft tissue limitations (particularly in pitch). Our results have a ventral bias, which may result from an effect of soft tissue such as synovial fluid in the cotyle–condyle joint, but resolving this requires further investigation. It is also possible that cartilage in the joints of the cotyle–condyle provides a slight increase in ROM of yaw and pitch. Our data obtained from the intact brown tree snake, boa constrictor and corn snake did not include forces generated when bending the snakes, which may differ from in vivo ranges. Nonetheless, the data from the brown tree snake, boa constrictor and corn snake matched closely with the observed intact values, supporting the validity of our results and suggesting minimal soft tissue influence in contrast to the observed influence of soft tissue in avian hip and limb joints (Baier, 2012; Kambic et al., 2017; Manafzadeh and Padian, 2018). Thus, biological validation should be used wherever possible.

ROM in joints relies on soft tissue and joint geometry for stability though the contribution of each differs between joints. Regardless of whether a system relies on soft tissue or joint geometry, the body typically has to constantly interact with that system, whether that is stiff connective tissue around the joint (Oliver et al., 2016), a deep joint cotyle or articular facets such as the pre- and post-zygapophyses. However, in snakes, the zygosphene only interacts at positions of high yaw or high dorsal pitch with no effect at most other positions. Through the interaction of multiple joints (i.e. zygosphene–zygantrum and pre- and post-zygapophyses) and a process (i.e. hypapophysis) or soft tissue, an overall functional response is achieved that limits yaw, dorsal pitch and ventral pitch even though any given component listed above may not actively interact at all joint positions. Thus, the zygosphene may allow snakes to avoid certain trade-offs between mobility and stability, while also providing an additional load-bearing structure at high yaw positions. Future work investigating other joint articulations that engage only at high ROM could provide insight into whether other species employ similar mechanisms.

We thank K. Sondereker for help with segmenting vertebrae, G. Nikolov and A. Knoll for their help with obtaining μCT scans of the vertebrae, and the American Museum of Natural History and University of Michigan for loaning us their specimens. We would also like to thank three anonymous reviewers for helping improve the manuscript.