Vertebral bending mechanics and xenarthrous morphology in the nine-banded armadillo (Dasypus novemcinctus)

Throughout evolutionary history, the vertebral column has adapted to support the disparate modes of locomotion found in vertebrates, ranging from axial swimming in fish, to high-speed running in several land animals, to powered flight in birds and bats (e.g. Boszczyk et al., 2001; Buchholtz, 2001; Cullinane and Bertram, 2000; Filler, 2007; Gál, 1993a,b; Gaudin and Biewener, 1992; Hebrank et al., 1990; Hildebrand, 1959; Jones and German, 2014; Jones and Pierce, 2016; Molnar et al., 2014; Pierce et al., 2011; Rockwell et al., 1938; Schilling, 2011; Shapiro, 1995; Slijper, 1946; Ward and Mehta, 2014). Vertebral evolution is facilitated by regionalization, the Hox-mediated division of the column into morphological and functional units (Head and Polly, 2015; Schilling, 2011; Wellik, 2007). In mammals, the thoracic and lumbar regions are foremost implicated in locomotion, prompting investigations into the functional link between thoracolumbar morphology and mechanical properties of intervertebral joints (e.g. Boszczyk et al., 2001; Gál, 1993a,b; Gaudin and Biewener, 1992; Granatosky et al., 2014; Jones, 2016; Long et al., 1997; Nyakatura and Fischer, 2010; Pierce et al., 2011; Russo, 2010; Shapiro, 1995, 2007; Slijper, 1946). These studies illustrate the morphological and mechanical diversity of thoracolumbar vertebrae within and between mammals, and demonstrate that morphological predictors of mechanics can vary according to size, taxon and mode of locomotion. Thus emphasized is the need for mechanical investigations into animals with varied morphologies and modes of locomotion, as well as the different roles that the thoracic and lumbar vertebrae play during movement.

Xenarthrans (armadillos, anteaters and sloths) are a fascinating group of animals in which to study vertebral mechanics, yet they have received limited attention. Xenarthrans are hypothesized to have evolved as ancestrally fossorial mammals (Simpson, 1931). Several traits characterizing this major mammalian lineage are used as evidence for this hypothesis, including the limb morphology of modern and fossil xenarthrans, the lack of color vision in living xenarthrans, and the xenarthrous articulation (Emerling and Springer, 2015; Frechkop, 1949; Jenkins, 1970; Nyakatura and Fischer, 2011; Olson et al., 2016; Vizcaíno and Milne, 2002). Xenarthrous articulations, or xenarthrae, are ancillary intervertebral articulations found across the posterior thoracic and lumbar vertebrae, spanning the post-diaphragmatic region of xenarthrans (e.g. Gaudin, 1999). These articulations are formed between the enlarged anapophysis of the cranial vertebra, and the enlarged metapophysis of the caudal vertebra in an articulating pair of vertebrae (Fig. 1A). Although enlarged metapophyses and anapophyses are not uncommon across modern mammals (Filler, 2007), the two only articulate this way in xenarthrans. In the fossil record, however, xenarthrous articulations are not limited to xenarthrans, and have evolved convergently in unrelated Mesozoic mammals (Luo and Wible, 2005; Martin et al., 2015). The condition of xenarthry is found in all known modern and fossil xenarthrans with two exceptions: modern sloths and extinct glyptodonts. Modern sloths have secondarily lost the articulation, and the corresponding vertebrae in glyptodonts are instead fused to each other and to the massive carapace (Gaudin, 1999).

Xenarthrous articulations are thought to stiffen the vertebral column (Frechkop, 1949; Gaudin, 1999; Gaudin and Biewener, 1992; Jenkins, 1970). By stabilizing the trunk, xenarthrae are hypothesized to help brace the body, thereby freeing and providing leverage to the forelimbs when digging (Frechkop, 1949; Jenkins, 1970). Additionally, xenarthrae have been suggested to increase the vertebral surface area available to resist compression, resulting in the safer transmission of axial forces generated by digging (Gaudin and Biewener, 1992). Gaudin and Biewener (1992) investigated the involvement of xenarthry in axial stability with a biomechanical study of the armadillo vertebral column. They found an increase in lateral and dorsal stiffness, as well as shear, in the xenarthrous region of the nine-banded armadillo (Dasypus novemcinctus) relative to the corresponding region in the opossum, their choice of ancestral mammal, thus corroborating Frechkop and Jenkins' hypotheses (Frechkop, 1949; Jenkins, 1970). They were, however, unable to directly link this increased stiffness to the xenarthrous morphology itself.

Here we further explore the mechanical implications of xenarthrous articulations by conducting a detailed examination of the functional association between xenarthrous morphology and vertebral bending mechanics in the nine-banded armadillo, D. novemcinctus Linnaeus 1758. To do this, we quantify multi-axis bending mechanics of excised thoracolumbar joints spanning the pre-xenarthrous and xenarthrous regions, and correlate this with the bony morphology of individual vertebrae. By examining both pre- and xenarthrous regions, we are able to add context to the findings of Gaudin and Biewener (1992), who focused on the mechanics of xenarthrous vertebrae alone. In addition to allowing a direct comparison with the mechanical work of Gaudin and Biewener (1992), we chose to study D. novemcinctus because it displays a generalized plesiomorphic xenarthrous morphology, is considered a good model for the evolution of xenarthry, and is readily accessible because of its wide ecological range and large population sizes (Gaudin, 1999; Loughry and McDonough, 2013). Our approach provides the opportunity to link xenarthrous vertebrae to the mechanical and morphological regionalization of the D. novemcinctus vertebral column, thereby affording more insight into the role of these articulations during movement.

We examined eight wild-caught adult D. novemcinctus specimens (Table 1) that were packed and frozen as soon as possible after collection. Armadillo carcasses were collected by W.J.L. under Georgia Department of Natural Resources Scientific Collecting Permit CN91349 and by Frank M. Knight under Arkansas Game and Fish Scientific Collecting Permit 112020121, and donated to the Museum of Comparative Zoology and Department of Organismic and Evolutionary Biology, Harvard University, for experimentation and future curation.

Specimens were thawed over two nights at temperatures just above freezing before dissection. Vertebral columns were excised and the thoracolumbar region was divided into seven ‘motion segments’ that each included two vertebrae connected by an intervertebral joint (Fig. 1B). The thoracic (T) region was divided into five motion segments (T1–2, T3–4, T5–6, T7–8 and T9–10; Fig. 1B), while the lumbar (L) region was divided into two motion segments (L1–2 and L3–4; Fig. 1B). Following division into motion segments, remaining musculature was removed from each motion segment while leaving the intervertebral ligaments and joint capsule intact. Motion segments displaying any damage to the vertebrae or intervertebral joint were not analyzed.

The first three thoracic motion segments articulate via pre-diaphragmatic zygapophyseal articulations, and are also defined here as pre-xenarthrous, as they do not participate in xenarthrous articulations. In D. novemcinctus, the diaphragmatic vertebra is T7, which is also the most cranial vertebra to participate in a xenarthrous articulation (grey vertebra, Fig. 1B). It is not uncommon, however, to find small, non-articulating anapophyses and metapophyses on T5 or T6. All motion segments caudal to T5–6 are post-diaphragmatic and xenarthrous, indicating that the constituent vertebrae are joined by xenarthrous articulations in addition to post-diaphragmatic zygapophyseal articulations. In D. novemcinctus, the component bones of the sacrum are fully fused, making it difficult to isolate the motion segment formed by L5 and the first sacral vertebra. Moreover, visualization of movement within the joint was hindered by the large ilia. We therefore did not include the lumbosacral motion segment in our quantitative analyses, and will instead describe lumbosacral mobility qualitatively.

Prior to experimentation, stainless steel screws were fastened into holes drilled into the cranial aspect of the anterior vertebral centra of all but the first motion segment (T1–2), and into the caudal aspect of the posterior vertebral centra of all motion segments following the protocol detailed by Molnar et al. (2014). Drill bits and screws were marked so as to prevent penetration of the joint capsule. Digital X-ray images were taken of the joints following experimentation to ensure that all screws were positioned appropriately. The centrum of T1 in D. novemcinctus is too dorsoventrally compressed to support screws, so locking surgical forceps were used instead to clamp the vertebra and secure its immobility during experimentation. In all motion segments, the cranial screw (or forceps) was fastened by a clamp while leaving the caudal vertebra free to move about the intervertebral joint during the application of load. Physiological saline was applied to the motion segments throughout the preparation and experimentation processes.

To record angular deflection of the caudal vertebra under load, pins were secured in the neural spines of both vertebrae forming the motion segment during ventroflexion, in the ventral aspect of the centra during dorsiflexion, and in the right diapophyses or transverse processes during lateroflexion (Fig. 2A,B). Intervertebral deflection was recorded with a camera positioned in plane with the vertebrae, mounted on a tripod to ensure height equality between the motion segment and the camera.

The caudal vertebra of each motion segment was loaded incrementally with increasing metric masses (g) in ventroflexion, dorsiflexion and lateroflexion (see Table S1 for applied masses). Prior to loading, a calibration image was taken of the prepared unloaded motion segment in neutral position with a scale bar oriented relative to gravity using an inbuilt level (Fig. 2A). The caudal vertebra was loaded until no further substantial deflection was detected, as judged by the experimenter. Immediately upon loading, a photograph was taken (Fig. 2B), the mass was quickly removed and the vertebra was prodded gently back into neutral position in preparation for the next load.

Following the procedure of Molnar et al. (Fig. 2A,B; fig. 4 in Molnar et al., 2014), eight points on each of the images were digitized to measure angular deflection of and applied moment on the caudal vertebra. An additional four points were digitized on each calibration image to orient the joint relative to gravity, and to approximate the center of rotation. Center of rotation was approximated as the midpoint of the intervertebral disc. The same individual digitized all images. Digitization error was calculated as the standard deviation of deflection angles in 10 repeated digitization trials of a single randomly selected image (s.d. error=0.040 deg).

where the angle of deflection (Δθ) is the angle between the two vectors created between the tips of each pin to the center of rotation in each loaded photograph (θ′_ab) minus the angle between the two vectors in the unloaded photograph (θ_ab) (as shown in Fig. 2A,B).

where Δθ is the angle of deflection in degrees, and M is the bending moment in newton meters (N m). Range of motion (R_m) is the asymptote at which the curve plateaus and is used as a measure of maximum range of motion. Compliance constant (C) is the rate constant that scales the curve, and is a measure of the rate at which the curve reaches its asymptote, R_m. C is thus proportional to compliance of the motion segment, or the inverse of stiffness. R_m and C were used to quantitatively describe bending mechanics across motion segments. For calculated R_m and C data, see Table S2.

Following experimentation, all motion segments were μCT-scanned on a SkyScan 1173 in the Museum of Comparative Zoology, Harvard University, reconstructed using NRecon version 1.6.6.0, and segmented and 3D volume rendered in Mimics Materialise software version 17.0. Thirteen linear measurements and seven angular measurements were recorded on each vertebra in Mimics by the same individual (Fig. 3). Measures were selected according to established determinants of stiffness and range of motion in other animals, and additional xenarthry-specific measures (metapophyseal and anapophyseal measures) were included in order to address the implications of xenarthry in compliance and range of motion.

Centrum height, neural spine height and pre-zygapophyseal angle have been suggested to be positively correlated with dorsoventral stiffness, centrum width and diapophysis/transverse process width positively correlated with lateral stiffness, and lamina width, pre-zygapophyseal width and transverse process/diapophysis angle positively correlated with stiffness in both planes. Centrum length has been suggested to be negatively correlated with dorsoventral stiffness, pre-zygapophyseal angle negatively correlated with lateral stiffness, and neural spine angle negatively correlated with stiffness in both planes. Dorsoventral transverse process/diapophysis angle and inter-zygapophyseal length have both been positively correlated with dorsoventral range of motion, and pre-zygapophyseal angle positively correlated with range of motion in lateroflexion. Centrum length, neural spine length, craniocaudal transverse process/diapophysis angle and neural spine angle have also been shown to be negatively correlated with range of motion in the sagittal plane (Buchholtz and Schur, 2004; Granatosky et al., 2014; Long et al., 1997; Molnar et al., 2014, 2015; Pierce et al., 2011).

Metapophyses and anapophyses and their corresponding linear and angular measures are, with some exception, specific to the xenarthrous region, and their influence on intervertebral mechanics is the focus of the present study. To address the almost total lack of metapophyses and anapophyses in pre-xenarthrous thoracic vertebrae, we followed the precedent set by Head and Polly (2015) and by Klingenberg (2008) in measuring the evolution of morphological novelty using homologous landmarks. Klingenberg (2008) suggested assigning landmarks to the fully formed structure of interest, and in its absence, to the position on the surface from which it would emerge. This method was coopted specifically for vertebrae by Head and Polly (2015), who were able to assign landmarks to vertebral morphologies not present in all measured vertebrae. This protocol was adapted here in pre-xenarthrous vertebrae, where measures concerning the metapophysis and anapophysis are set to zero.

where M_adj is the size-adjusted measurement, M is the original measurement, L_S is the mean nose–rump length in all specimens, L₀ is the nose–rump length of an individual and b is the slope of the regression of log(M) on log(L₀) for each measurement in all specimens (Elliott et al., 1995). Normalized, size-removed measures were used in all subsequent analyses. For raw measurement data, see Tables S3 and S4.

where a significant interaction term indicates that craniocaudal patterns vary between bending directions. Compliance constant data were log transformed so as to increase the normality of their distribution. After an initial exploration of the data, a combined 64 out of 306 compliance constant and range of motion data points were removed because of low R² values of a motion segment (R²<0.75), or errors in data collection. The ANOVAs were used to identify significant differences in compliance constant and range of motion between motion segments along the column, and between bending directions (independent variables). To locate differences, post hoc pairwise comparisons were performed with Bonferroni corrections for multiple comparisons.

A principal component analysis was performed on all normalized anatomical measures in order to identify which changes in morphology were most heavily implicated in shape variation and also to recover a series of independent and uncorrelated variables for comparison with the bending data. Following this, the first three principal components were regressed onto the compliance constant and range of motion data in each bending direction in a series of stepwise linear regressions. The first three principal components were chosen as independent variables because they each explained a minimum of 5% of the morphological variation across the thoracolumbar region and each represents a major transition across the vertebral column (see Results). Principal components were included in regression models if they contributed significantly to the R² value of said models. Further, each measure was regressed independently onto compliance constant and range of motion data to assess which individual measurements were most predictive of mechanical properties. All statistical analyses were carried out with SPSS version 22.0.

Changes in compliance constant (C) across the column followed a similar pattern in all three bending directions: the first motion segment exhibited low compliance, and was followed caudally by an increase, and then a decrease in compliance until T7–8, at which point the compliance constant steadied to a low plateau (Fig. 4A, see Table S2). This plateau of low compliance in all three directions spanned the entire measured post-diaphragmatic, xenarthrous portion of the thoracolumbar region. Although not quantified, visual inspection revealed that the lumbosacral joint formed by the fifth lumbar vertebra and the sacrum departed from this trend, and was very compliant in dorsiflexion and ventroflexion. In contrast, its mobility was totally restricted in lateroflexion by the ilia.

Spanning the measured motion segments, a significant difference in compliance constant was found between bending directions (Table 2, see Table S5). T3–4 was significantly less compliant in dorsiflexion than in ventroflexion (P=0.008), and T1–2, T5–6 and T7–8 were significantly less compliant in dorsiflexion than in both lateroflexion (P=0.033, 0.001, 0.008) and ventroflexion (P<0.001, <0.001, =0.004; Table 2). Compliance in all three directions was statistically indistinguishable in the region spanning T9–L4 (Table 2).

Compliance constants were significantly different across motion segments in each direction (Table 3, see Table S5). In dorsiflexion, T3–4 was significantly more compliant than T7–8 (P<0.001), L1–2 (P=0.001) and L3–4 (P=0.015; Table 3). In lateroflexion, T1–2 was significantly more compliant than L3–4 (P=0.036; Table 3). Caudally, T3–4 and T5–6 were significantly more compliant than the region spanning T7–L4 (Table 3). In ventroflexion, with the exception of significantly indistinguishable compliance constants in T5–6 and T9–10 (P=0.083), the region spanning T1–T6 was significantly more compliant than the region spanning T7–L4 (Table 3). Additionally, compliance in ventroflexion was significantly higher in T9–10 than in L3–4 (P=0.037; Table 3).

As with compliance, range of motion (R_m) followed a similar pattern in all three bending directions (Fig. 4B, see Table S2). R_m was consistently lowest in T1–2, and increased caudally to reach a plateau at T9–10. A local peak in R_m was found at T3–4 in all three directions, and was especially pronounced in ventroflexion (Fig. 4B).

Caudal to T1–2, range of motion was almost consistently higher in ventroflexion than in dorsiflexion and lateroflexion (Fig. 4B). R_m was higher in dorsiflexion than lateroflexion in all motion segments except for T3–4 and T5–6 (Fig. 4B). R_m in ventroflexion was significantly higher than in dorsiflexion and lateroflexion in T3–4 and L1–2 (Table 2). In T5–6, R_m in ventroflexion was significantly higher than in dorsiflexion (P=0.004), and in T9–10, R_m in ventroflexion was significantly higher than in lateroflexion (P<0.001; Table 2).

Range of motion differed significantly across the column in all three bending directions (Table 3, see Table S5). In dorsiflexion, the region spanning T1–T6 was characterized by a significantly lower R_m than that of T9–L4, and R_m in T7–8 was significantly lower than in T9–10 (Table 3). In lateroflexion, R_m in T1–2 was significantly lower than all motion segments except T5–6 (Table 3). R_m was significantly lower in T5–6 than in the motion segments spanning T9–L4 (Table 3). In ventroflexion, R_m in T1–2 was significantly lower than in all caudad motion segments, and R_m in T3–4 was significantly higher than in the two motion segments caudad (Table 3). R_m in T5–8 was significantly lower than in T9–L2 (Table 3). A significant interaction was found between motion segment and direction (see Table S5).

A principal component analysis conducted on the morphological measures found that the first three components explained 82.7% of the variance in the thoracolumbar vertebral morphology of D. novemcinctus (Table 4, Fig. 5). The first three components were the only components with corresponding eigenvalues greater than 1, and that explained more than 5% of the variance (Table 4, see Table S6 for scores).

The first component (PC1), which explained 53.0% of the variance, effectively distinguished between pre-diaphragmatic pre-xenarthrous, and post-diaphragmatic xenarthrous vertebrae (Table 4, Fig. 5). PC1 was most heavily weighted (≥±0.75) in six anapophyseal, metapophyseal and zygapophyseal measurements, with all measures showing substantial changes in value between T6 and T7, or between T7 and T8 (Table 4, Fig. 6A–C). As the cranialmost anapophysis and metapophysis are typically found in T7, with exceptional presence in T5 and T6, anapophysis length (AL) and height (AH), and metapophysis angle (MA) increased starkly from values close to zero prior to T7 (Fig. 6A,B). Metapophysis height (MH) followed a similar trend, although it increased steadily rather than sharply from T6 caudad, decreasing slightly in L5 (Fig. 6C). Inter-zygapophyseal length (IZL) increased and post-zygapophyseal length (PozL) decreased at T7, reflecting the post-diaphragmatic, post-zygapophyseal morphology of the diaphragmatic vertebra, T7 (Fig. 6C).

The second component, which explained 21.7% of the variance, was most heavily weighted in diapophysis and transverse process measures (Table 4, Fig. 5), with transverse process/diapophysis craniocaudal angle (TPDAC) decreasing and transverse process/diapophysis dorsoventral angle (TPDAD) increasing between T10 and L1 (Fig. 6D). These changes mark the transition from diapophysis to transverse process, which signifies a morphological switch from thoracic to lumbar vertebrae.

The third component (PC3), which accounted for 8.0% of the variance, was most heavily weighted in centrum width (CW) and pre-zygapophyseal width (PrZW), with both measures being highest in T1 and decreasing to reach a somewhat steady level by T3 (Table 4, Figs 5, 6E). PC3 therefore appears to distinguish the most anterior thoracic vertebrae from the remaining thoracolumbar region.

A series of stepwise linear regressions using the first three principal components derived above as independent variables were performed on compliance constant and range of motion data in all three bending directions. Models composed of PC1 and PC2 were found to be predictive of compliance in all three directions, and of range of motion in dorsiflexion and lateroflexion (Table 5). A model composed of PC1 and PC3 was found to be predictive of range of motion in ventroflexion (Table 5). Several correlated measures were found to be predictive of compliance and range of motion in all three directions (see Tables S7, S8).

To illustrate the relationship between morphological regionalization and bending mechanics in the vertebral column of the nine-banded armadillo, we examined morphology, compliance and maximum range of motion during multi-axis bending in thoracolumbar vertebral joints. Although there are significant differences in mechanical properties between bending directions, compliance and range of motion follow similar patterns in all three directions: compliance decreases and range of motion increases caudally along the column in both dorsoventral (sagittal) and lateral flexion. These patterns in joint mechanics correspond to morphological regionalization prescribed by both the transition from pre-diaphragmatic, pre-xenarthrous vertebrae to post-diaphragmatic, xenarthrous vertebrae, and the transition from thoracic to lumbar vertebrae.

Of note is the first thoracic joint (T1–2), which is characterized by a very low range of motion and high compliance (Fig. 4). These traits define the intervertebral joint formed by T1 and T2 as one that will readily bend, but is restricted to a very small degree of bending, especially in lateroflexion and ventroflexion. Centrum and pre-zygapophyseal width are both substantially higher in the first two thoracic vertebrae than in all other thoracic vertebrae (Fig. 6E). Additionally, although not measured, visual observation revealed that the intervertebral disc formed between T1 and T2 is substantially thinner than in all caudad intervertebral joints. Such a wide and thin area of contact between centra has been suggested to restrict mobility of vertebral motion segments (Buchholtz and Schur, 2004; Long et al., 1997). This morphology and lack of mobility are also characteristic of the cervical region of D. novemcinctus. Although we did not measure compliance or range of motion in cervical vertebrae, qualitative manipulation revealed that it contains only two centers of flexion: between the atlas and axis, and between the sixth and seventh cervical vertebrae. Excluding the third and fourth cervical vertebrae, which are fused to the axis in the mesocervical bone characteristic of armadillos (Galliari et al., 2010), there is remarkable rigidity from the fourth to sixth cervical vertebrae, and between the seventh cervical and first thoracic vertebrae. With the exception of the long neural spines of the first two thoracic vertebrae, the morphology of the joint formed between them is remarkably similar to those in the cervical region. The cervical-like morphology and rigidity of T1–2 is a reflection of the transitional nature of this intervertebral joint, and especially of the first thoracic vertebra itself. In Mus musculus, the first thoracic vertebra is morphologically more similar to the cervical vertebrae than to the other thoracic vertebrae, and is also located in a position of overlap of Hox genes prescribing cervical and thoracic patterning (Head and Polly, 2015; reviewed in Wellik, 2007). The morphological and mechanical similarities between most of the cervical vertebrae and the first two thoracic vertebrae suggests a shared function for these vertebrae.

Within each bending direction, thoracic joint T3–4 is the most compliant of all motion segments, and has the highest range of motion of all pre-diaphragmatic motion segments (Fig. 4). Additionally, visual examination during experimentation revealed that T3–4 undergoes rotational motion quite readily, particularly during lateroflexion. Together, these traits define a highly flexible and mobile joint. The high compliance found in T3–4 fits in with the other pre-xenarthrous thoracic motion segments. Barring T1–2, which has a cervical nature, T3–4 is part of a series of compliant motion segments that decreases in compliance until reaching a post-diaphragmatic plateau. Its range of motion, however, deviates from the low values of the pre-xenarthrous thoracics in ventroflexion (Fig. 4B). While range of motion in ventroflexion is predicted by a model composed of PC1 and PC3, neither principal component can account for such a spike at T3–4. Indeed, while the regression model composed of PC1 and PC3 is significantly predictive of range of motion in ventroflexion, it had the least support of all directions (Table 5). Therefore, although our methods appear successful in developing predictive models of compliance in all three bending directions, and of range of motion in lateroflexion and dorsiflexion, they are not as successful at predicting range of motion in ventroflexion. This shortcoming agrees with findings from domestic horses that bony morphology is not as powerful at predicting range of motion in ventroflexion as in dorsiflexion (Jones, 2016). Instead, architecture and elasticity of the supraspinous ligament and ligamenta flava are likely more predictive of range of motion in ventroflexion (Gál, 1993b; Jones, 2016). Gál (1993b) found that the ligamenta flava, the ligament extending from axis to sacrum that connects successive vertebral laminae, was a primary component of resistance to flexion in several small mammals, whose size range includes D. novemcinctus. To identify the soft tissue drivers behind range of motion in ventroflexion in D. novemcinctus, it may be worthwhile to perform successive lesion experiments.

Upon xenarthry, compliance and range of motion reach a relative plateau in all three bending directions (Fig. 4), with post-diaphragmatic, xenarthrous motion segments being characterized by a similarly low compliance (from T7–8 caudad) and high range of motion (from T9–10 caudad). The decrease in compliance clearly aligns with the initial appearance of xenarthry and the diaphragmatic transition at the T7–8 motion segment, and is supported by the prominence of PC1 in regression models predicting compliance (Table 5, Figs 4, 6). These results support the theory that the xenarthrous region confers a lower compliance, or a higher stiffness, to the xenarthran vertebral column (Frechkop, 1949; Gaudin, 1999; Gaudin and Biewener, 1992; Jenkins, 1970). Xenarthrous articulations thus appear to function as a stiffening mechanism in concert with the post-diaphragmatic zygapophyseal articulations found on the same vertebrae.

In their investigation into vertebral bending mechanics in D. novemcinctus and the Virginia opossum, Didelphis virginiana, Gaudin and Biewener (1992) found that the xenarthrous region of D. novemcinctus is stiffer in lateroflexion and dorsiflexion than is the corresponding region of D. virginiana. In contrast, they found that bending in ventroflexion was significantly less stiff than in lateroflexion and dorsiflexion, and was comparable to that in D. virginiana. Additionally, they found no difference in the degree of axial rotation between D. novemcinctus and D. virginiana. Our results appear to both support and contradict those of Gaudin and Biewener, as we found that the higher stiffness in lateral and dorsal bending is statistically indistinguishable from the similarly high stiffness in ventral bending in the region spanning T9–L4 (Table 2). This difference may be accounted for by the methods employed in our studies. In their protocol, Gaudin and Biewener measured stiffness in vertebrae including and caudal to T6 as a single unit connected by ligaments, while we measured compliance across the thoracolumbar region in units of two vertebrae. As our results demonstrate, T6, and to some extent T7, are more similar morphologically to craniad pre-xenarthrous vertebrae than to caudad xenarthrous vertebrae (Fig. 6A–C), and the joint formed between T5 and T6 is also mechanically more similar to more cranial vertebrae (Table 3, Fig. 4). Although we did not measure compliance in the joint formed between T6 and T7, based on morphology, we assume that its compliance is likely intermediate between that of T5–6 and T7–8 (Fig. 6). In Gaudin and Biewener's results, the stiffness imparted to the column by the intervertebral joint connecting T6 and T7 cannot be distinguished from the resultant stiffness of all caudad intervertebral joints. It is therefore possible that the relative compliance of the cranialmost intervertebral joint measured by Gaudin and Biewener masks the stiffness of the remainder of the column, resulting in a low stiffness in ventroflexion that is comparable to that of D. virginiana.

In contrast to compliance, the link between high range of motion and post-diaphragmatic, xenarthrous morphology is less striking, because of the apparent misalignment between the diaphragmatic vertebra, T7, and the increase in range of motion at T9–10 (Figs 4B, 6A–C). Nonetheless, our data do implicate xenarthrous articulations and post-diaphragmatic zygapophyses with an increase in range of motion, as PC1, which was heavily weighted in xenarthrous and zygapophyseal measures, was found to be predictive of range of motion in all three bending directions (Tables 4, 5). It appears that, rather than acting as bony stops, xenarthrous articulations may help to amplify the observed effects of the curved, protruding post-diaphragmatic zygapophyses in mammals, which permit and secure flexion by increasing the surface area upon which vertebrae can slide without disarticulating, thereby enhancing range of motion (Filler, 2007; Jenkins, 1974; Pierce et al., 2011; Russo, 2010; Shapiro, 1995). In functioning in this way, xenarthrous articulations may extend beyond a purely stabilizing role and increase range of motion by providing greater surface area for movement. Hebrank et al. (1990) found a similar result in blue marlins. In a manner analogous to xenarthrans, marlin vertebrae have a unique interlocking morphology that allows for bending in lateroflexion, while restricting sagittal flexion and extension. Hebrank et al. (1990) suggest that these interlocking facets increase lateral bending while preventing disarticulation of the involved vertebrae. Further analysis into regional variation in bending mechanics in mammals without xenarthrous articulations would help clarify the contributions of xenarthrous articulations to range of motion, thereby elucidating whether xenarthrae enhance the suggested effects of post-diaphragmatic zygapophyses on mobility.

In addition to intervertebral articulations, we propose that enlarged metapophyses influence the high range of motion seen in the motion segments caudal to T7–8. Though metapophyses are present on T7 and T8, they are quite short (Fig. 1B, Fig. 6C, MH). Caudal to T8, metapophyses increase in height to become massive processes with the capacity for increased ligament and muscle insertion, thus enabling a higher range of motion (Buchholtz and Schur, 2004; Granatosky et al., 2014; Pierce et al., 2011; Shapiro, 1995). As evidenced by the involvement of PC2 in models predicting both compliance and range of motion, lumbar morphology is also correlated with higher range of motion (Table 5). Our results agree with those of Gaudin and Biewener (1992), who found that in D. novemcinctus, lumbar vertebrae reach higher levels of angular displacement than do thoracic vertebrae. Interestingly, Gaudin and Biewener found that thoracic and lumbar vertebrae in D. virginiana do not differ markedly in angular displacement, thus substantiating our claim that the xenarthrous and metapophyseal morphology of D. novemcinctus confer an increase in range of motion relative to a purported ancestral mammal.

Investigations into stiffness and compliance of the vertebral column, either as a whole unit or separated into distinct motion segments, are common (e.g. Gaudin and Biewener, 1992; Hebrank et al., 1990; Granatosky et al., 2014; Long et al., 1997; Molnar et al., 2014). Also common are data on maximum range of motion achievable by manual flexion or by simulation (e.g. Jones, 2016; Jeffcott and Dalin, 1980; Molnar et al., 2015; Townsend et al., 1983). However, the interplay between compliance or stiffness and range of motion is one that is often overlooked (Gál, 1993a; Molnar et al., 2015). Our study has shown that xenarthrous articulations do indeed appear to stiffen the armadillo vertebral column, but they may also confer an increase in range of motion. Consequently, the xenarthrous region of D. novemcinctus is very stiff, but also highly mobile under a substantial enough load. Gál (1993a) also found that stiff vertebral columns do not necessarily preclude mobility. Much like in post-diaphragmatic bending in D. novemcinctus, lumbosacral bending in monkeys and wallabies is characterized by relatively low compliance and high range of motion (Gál, 1993a). As stressed by these results, it is imperative that further investigations into vertebral bending address both compliance/stiffness and range of motion to arrive at comprehensive conclusions on bending mechanics.

The stiffness and range of motion conferred by xenarthrous articulations have significant implications for movement. Although Gaudin and Biewener (1992) concluded that the increase in lateral and dorsal stiffness they found in D. novemcinctus would aid in resisting the laterally and dorsally directed forces generated by digging (Gasc et al., 1986), they argued that the low stiffness they found in ventroflexion does not support Frechkop's claim that stiffness in ventroflexion would be necessary to support the weight of the forelimbs and trunk during digging (Frechkop, 1949). In contrast, our results suggest that stiffness in all three directions is fostered by xenarthrous vertebrae, allowing for the resistance of forces generated ventrally by the weight of the forelimbs and trunk, dorsally by the digging stroke against the substrate, and laterally by the alternating digging stroke (Frechkop, 1949; Gasc et al., 1986; Gaudin and Biewener, 1992). Passive stiffness of the xenarthrous region of the D. novemcinctus vertebral column thus precludes the need for substantial trunk muscle exertion when digging. An analogous role for passive vertebral stiffness in maintaining posture has been described by Smeathers and Gál in large mammals (Gál, 1993a; J. E. Smeathers, A mechanical analysis of the mammalian lumbar spine, PhD thesis, University of Reading, 1981).

Upon active muscle exertion, the xenarthrous region has the potential to become highly mobile because of an increase in joint range of motion. Such mobility is advantageous to an animal whose default vertebral condition is rigidity, as it facilitates the non-digging behaviours typical of mammals generally and of armadillos specifically, such as running, jumping, kicking and the rolling behaviour found in select species. Furthermore, the xenarthrous region roughly overlaps with the flexible banded portion of the nine-banded armadillo's carapace, suggesting that this potential for mobility is not limited by the carapace in D. novemcinctus. The extent to which xenarthry-associated range of motion can be realized in armadillos, however, likely varies according to species, reflecting the variation in carapace shape (Superina and Loughry, 2012). Also expected to be varied among armadillos is the prominence of the xenarthrous morphology itself. Armadillos demonstrate a wide range of digging behaviours, from fully subterranean to almost entirely cursorial (Vizcaíno and Milne, 2002). Although xenarthrous morphology has been described extensively in both orders of Xenarthra (Gaudin, 1999), a quantitative examination of relative size of xenarthrous articulations between species is lacking, as is a functional investigation into the relationship between extent of xenarthry and degree of fossoriality across armadillo species. Following the proposed relationship between xenarthry and fossoriality, we expect there to be a positive correlation between the relative prominence of xenarthrous articulations and the degree of fossoriality across armadillos. As such, we also expect the most fossorial armadillos (e.g. pink fairy armadillo, giant armadillo) to be characterized by xenarthrous regions of the highest stiffness and range of motion.

The dual role of the post-diaphragmatic region in D. novemcinctus is conferred by the unique shape of the xenarthrous articulation, and by the associated enlarged metapophyses. The xenarthrous articulation itself stiffens and imparts mobility to the region through increased articular surface area. Adding to this, the metapophysis provides more room for ligament and muscle attachment and increases mechanical advantage, thereby facilitating the powered movement necessary for flexion under xenarthry (Buchholtz and Schur, 2004; Granatosky et al., 2014; Pierce et al., 2011; Shapiro, 1995). Our results demonstrate that the post-diaphragmatic region of D. novemcinctus is characterized by decreased compliance and increased range of motion, two traits that appear to be especially conducive to a semi-fossorial lifestyle.

We thank Judy Chupasko, Mark Omura and Mark Renczkowski [Museum of Comparative Zoology, Harvard Univeristy (MCZ)] for all of their help during the course of this project, Frank M. Knight (University of the Ozarks) for donating armadillo carcasses, and several colleagues in the MCZ for helping with experimental setup and data collection (Robert Kambic, Brianna McHorse, Blake Dickson and Hanna Barnes). A version of this paper was submitted by J.D.O. in partial fulfillment of the Erasmus Mundus Master Programme in Evolutionary Biology.