Phenotypic integration in the carnivoran backbone and the evolution of functional differentiation in metameric structures

Explaining the origin and evolution of a vertebral column with anatomically distinct regions that characterizes the tetrapod body plan provides understanding of how metameric structures become repeated and how they acquire the ability to perform different functions. However, despite many decades of inquiry, the advantages and costs of vertebral column regionalization in anatomically distinct blocks, their functional specialization, and how they channel new evolutionary outcomes are poorly understood. Here, we investigate morphological integration (and how this integration is structured [modularity]) between all the presacral vertebrae of mammalian carnivorans to provide a better understanding of how regionalization in metameric structures evolves. Our results demonstrate that the subunits of the presacral column are highly integrated. However, underlying to this general pattern, three sets of vertebrae are recognized as presacral modules—the cervical module, the anterodorsal module, and the posterodorsal module—as well as one weakly integrated vertebra (diaphragmatic) that forms a transition between both dorsal modules. We hypothesize that the strength of integration organizing the axial system into modules may be associated with motion capability. The highly integrated anterior dorsal module coincides with a region with motion constraints to avoid compromising ventilation, whereas for the posterior dorsal region motion constraints avoid exceeding extension of the posterior back. On the other hand, the weakly integrated diaphragmatic vertebra belongs to the “Diaphragmatic joint complex”—a key region of the mammalian column of exceedingly permissive motion. Our results also demonstrate that these modules do not match with the traditional morphological regions, and we propose natural selection as the main factor shaping this pattern to stabilize some regions and to allow coordinate movements in others.

Nevertheless, although spine evolution has been recently studied in nonmammalian synapsids (e.g., Jones et al. 2018aJones et al. ,b, 2019Jones et al. , 2020, much is yet unknown on the advantages and costs of vertebral column regionalization in anatomically distinct blocks, its association with functional specialization, and how it channel new evolutionary outcomes. To this respect, studies on phenotypic integration (Olson and Miller 1958) between vertebrae (covariation) and how this integration is structured (modularity) may provide crucial insights to understand how regionalization in metameric structures evolves. Phenotypic integration originates from trait interactions at genetic, developmental, or functional levels (e.g., Goswami et al. 2014;Klingenberg 2014;Martín-Serra et al. 2015), and it provides important clues as to how these aspects influence phenotypic macroevolution (Goswami et al. 2014). Accordingly, phenotypic integration in the column of mammals and other amniotes has been investigated to understand how pleiotropy, developmental constraints, and functional adaptations influence the evolution of regionalization-and functional differentiation (e.g., Buchholtz 2007Buchholtz , 2012Randau et al. 2017;Randau and Goswami 2017a,b;Jones et al. 2018a,b;Polly et al. 2001;Ward and Brainerd 2007;Head and Polly 2015).
Among mammals, both Randau and Goswami (2017a) and Jones et al. (2018b) quantified the integration among vertebrae in the spine of felids and in the thoracolumbar region of mammals, respectively. Randau and Goswami (2017a) identified a five-modules pattern (plus a no-adjacent cervical-lumbar module) that reflects shared developmental pathways, ossification timing, and ecological diversification. On the other hand, Randau and Goswami (2017a) demonstrate a pattern of increasing disparity in the postdiaphragmatic region, and Randau et al. (2017) found stronger ecological signal in the lumbar of felids. Jones et al. (2018b) demonstrated that lumbar and diaphragmatic vertebrae present a higher ecological signal coupled with higher disparities and evolutionary rates, suggesting a link between adaptation and evolvability in the lumbar region. The authors concluded that thoracic and lumbar vertebrae form two semi-independent modules.
Despite both studies (Randau and Goswami 2017a;Jones et al. 2018b) found different modules, the same pattern of increasing disparity, evolvability, and adaptation in postdiaphragmatic vertebrae is evident. Although Randau and Goswami (2017a) performed the study at family level with 19 vertebrae of a total of 27 presacrals, Jones et al. (2018b) conducted the study at a wider phylogenetic scale (at class level) but restricted to five presacral vertebrae.
Here, we investigate phenotypic integration among all presacral vertebrae using a wide dataset with representatives of most families that belong to the mammalian order Carnivora. Our main goal is to investigate the modular organization of the vertebral column in carnivores and to discuss on the possible factor(s) shaping the pattern of integration in metameric structures. Specifically, we answer the following questions: (i) Is there a pattern of modular organization of specific subsets of vertebrae? If so, (ii) does it match with the classical morphological regions (cervical, thoracic, lumbar) or does it match with the more recently proposed five-region (cervical, pectoral, thoracic, diaphragmatic, lumbar) model of Jones et al. (2018b) and possibly with the extra cervical-lumbar region proposed by Randau et al. (2017)? And finally (iii) to hypothesize on those possible factors driving macroevolutionary patterns in metameric structures.

Methods
All presacral vertebrae of 41 species of mammalian carnivores ( Fig. 1; Table S1) were scanned in three dimensions with a medical CT or using a NextEngine ® surface scanner (http://www. nextengine.com). This resulted in 1107 models of vertebrae that were digitally processed with Meshlab (Cignoni et al. 2008).
We digitized a set of homologous landmarks (Table S2; Fig. 2) with the software Landmark from IDAV (Wiley et al. 2005) and the raw landmark data were incorporated into geomorph (Adams et al. 2017).

VERTEBRAE
Exploring integration analyses across the carnivoran vertebral column in a pairwise fashion is challenging because thoracics and lumbars are subject to homeotic variation (differences in thoracic and lumbar counts among species but keeping constant the total thoracolumbar count at 20 [i.e., 13T+7L, 14T+6L, 15T+5L, or 16T+4L]). We solve this problem using different procedures for grouping vertebrae (Fig. 3).
A joined thoracolumbar count procedure Counting cervicals and directly grouping the thoracolumbar vertebrae as dorsal vertebrae. This approach allows including all presacrals, allowing to keep adjacency among all the subunits (Fig. 3A). However, it mixes thoracic with lumbar vertebrae at the thoracolumbar boundary. The diaphragmatic vertebra is also variable in position, and hence, it mixes pre-and postdiaphragmatics. Therefore, this procedure leads to group together vertebrae that are not equivalent from a functional point of view across species.
Seven datasets were created for the C01-C07 with a fixed count and 20 datasets for TL01-TL20. Therefore, the TL14, TL15, or TL16 can be the first lumbar for some species but the last thoracic for others (Fig. 3A).
A thoracolumbar boundary count procedure To do not group nonhomologous vertebrae at the thoraciclumbar boundary that could bias integration, we counted from the thoracic-lumbar transition 13 thoracics in cranial direction and four lumbars in caudal direction (the maximum number of vertebrae shared by all species) (Fig. 3B). The last thoracics were joined together regardless of their position (TL13, TL14, TL15, or TL16), as well as the first lumbars. This procedure warrants homology among all grouped vertebrae in terms of regional morphology (thoracic vs. lumbar), and on position with respect to this regional transition. Moreover, it keeps the adjacency among vertebrae. However, it is not possible to include all the vertebrae because many species have more than 13 thoracics and/or more than four lumbars, and it still mixes together pre-and postdiaphragmatic vertebrae.   A diaphragmatic start count procedure

EVOLUTION LETTERS
We counted from the diaphragmatic vertebra defined as the last thoracic with vertebrosternal ribs and with pre-and postzygapophyses oriented in horizontal and vertical planes, respectively (Randau and Goswami 2017a) (Fig. 3C). All diaphragmatic vertebrae are joined together regardless of their position and the remaining vertebrae are established as pre-and postdiaphragmatic in cranial and caudal direction, respectively. The homology criterion is based on the diaphragmatic vertebra and on vertebral position. Therefore, it does not mix diaphragmatic vertebrae with pre-and postdiaphragmatic vertebrae, and it keeps adjacency in most vertebrae. However, the procedure still mixes thoracics with lumbars and it does not include all the presacral vertebrae, as the number of pre-and postdiaphragmatics vary among species ( Fig 3C).
A selected vertebrae procedure We selected the most important vertebrae according to their position and/or morphology regardless of the number of vertebrae between them: first thoracic, intermediate thoracic (TL06 if the species has 13 or 14 thoracics and TL07 if the species has 15 or 16 thoracics), diaphragmatic, last thoracic, first lumbar, and last lumbar (Fig. 3D). Therefore, following this procedure, homology is based on relative position and/or key anatomical markers, and EVOLUTION LETTERS JUNE 2021 2 5 5 it does not mix thoracics with lumbars or pre-with postdiaphragmatic vertebrae. However, this method only allows including a very small fraction of anteroposterior variation across the column. In addition, the comparisons are not completely equivalent due to adjacency differences (e.g., the last thoracic and the first lumbar are adjacent but the first and last lumbars are not, which could bias integration patterns), and it does not take into account the number of vertebrae between them (e.g., the number of vertebrae between the first and the last lumbar varies among species).

PRELIMINARY ANALYSES
For each vertebral position, we performed a Procrustes alignment accounting for bilateral symmetry (Adams et al. 2017). Afterward, two species (Mustela putorius and Suricatta suricatta) that were represented by two individuals were averaged. We computed a phylogenetic Procrustes ANOVAs for each vertebra using log-transformed centroid size as the independent variable. The phylogeny included for this analysis was taken from Nyakatura and Bininda-Emonds (2012).

INTEGRATION ANALYSES
The covariation analyses were performed from the PGLS residuals for each vertebra. We used two block Partial Least Squares (2B-PLS) analyses to explore morphological covariation between vertebrae and we compared each vertebra with the others. These comparisons were also performed to account for possible phylogenetic effects (2B-PLS phy ). Given the high number of comparisons, we used the Benjamini-Hochberg method to avoid false rejections of null hypothesis for each set of PGLSs (Benjamini and Hochberg 1995). We used the Z-score as a proxy for the strength of integration between pairs of vertebrae (Adams and Collyer 2016). To test the hypothesized boundaries between modules, we selected one standard comparison for each module (i.e., the maximum Z-score between adjacent vertebrae) and its Z-score was compared with the Z-score obtained from the comparison between pairs of vertebrae at hypothesized boundaries. The location of the boundary will be confirmed with changes from nonsignificant to significant in the P-value obtained from such comparisons.
To visualize the covariation between vertebrae, we modeled the shape variation accounted by the first axis of variation for some selected PLS comparisons only for the joined thoracolumbar method.
All the statistical analyses were performed with geomorph package (Adams et al. 2017) in R (R Core Team 2017).

Results
The phylogenetic Procrustes ANOVAs performed for each vertebra yielded varying significant results depending on the count-ing approach (see Tables S3-S6), and we used the allometric-free residuals for all the vertebrae to keep consistency among each set of analysis.

THE JOINED THORACOLUM BAR COUNT APPROACH
Both 2B-PLS and 2B-PLS phy analyses yielded significant results for most of the comparisons (Table S7; Fig. 4A). The nonphylogenetic Z-scores indicate three main regions of high integration (Table 1): (i) the cervical region, from C01 to C07 (being C02 relatively independent); (ii) the anterior thoracolumbar vertebrae, from C06-C07 to TL09; and (iii) the posterior thoracolumbar vertebrae, from TL12 to TL20. There is also a region of weak integration composed by TL10-TL11. The Z-scores obtained from the 2B-PLS phy show similar results. However, there is a general decrease of integration in the 2B-PLS phy that also reduces the significance of Z-scores differences, although the cervicothoracic and anterior-posterior thoracic boundaries are confirmed (Table 1).
To provide a morphological interpretation of this pattern of integration, we selected some pairwise comparisons to visualize those shape changes associated with vertebrae covariation. Three main types of pairwise comparisons were selected: (i) the C01-C02 due to their unique morphology among cervicals; (ii) representative between-vertebrae comparisons of integrated regions (cervical, anterodorsal, and posterodorsal modules); and (iii) comparisons between transitional vertebrae located in regional boundaries ( Fig. 5; Supporting Information).

APPROACH
Most of comparisons were significant (Table S8). Again, the Z-scores indicate the presence of two integrated regions for thoracics and lumbars, as well as a weakly integrated region composed by T04-T03 (starting from the last thoracic; Fig. 4B; Table 2). The P-values confirm the hypothesized boundaries at T04 or T03 (Table 2).

THE DIAPHRAGM ATIC START COUNT APPROACH
Most of the 2B-PLS and 2B-PLS phy were significant (Table S9; Fig. 4C). They indicate the presence of two integrated regions, anterior thoracic and posterior thoracolumbar, separated by the diaphragmatic vertebra, which is weakly integrated. The P-values of the 2B-PLSs confirm that the diaphragmatic vertebra does not belong to any integrated region, but this pattern is not as evident with the 2B-PLS phy (Fig. 4C; Table 3).

THE SELECTED VERTEBRAE APPROACH
Again, most of the 2B-PLS and 2B-PLS phy were significant (Table S10). Moreover, there is an integrated posterior thoracolumbar region, composed by the last thoracic and the first 2 5 6 EVOLUTION LETTERS JUNE 2021 and last lumbars (Fig. 4D). The anterior thoracic region is less clear than with the previous counting methods. The diaphragmatic vertebra has the weakest integration with its closest vertebrae. The P-values indicate, at least for the 2B-PLSs, that it does not belong to the posterior thoracolumbar region (Table 4).

Discussion
Almost all pairwise comparisons between vertebrae were significant (Tables S7-S10), suggesting that the mammalian vertebral column is a highly integrated structure despite its high degree of regionalization. However, this integration is not uniformly distributed along the vertebral column.
EVOLUTION LETTERS JUNE 2021 2 5 7 The results of the joined thoracolumbar count support three highly integrated regions (hereafter modules; e.g., Klingenberg and Marugán-Lobón 2013) in the presacral column of carnivores (Fig. 4A): C01-C07 (cervical), TL01-TL09 (anterodorsal), and TL13-TL19 (posterodorsal) and a region of weak integration composed by TL10-TL11. This indicates that the anterodorsal module corresponds with the prediaphragmatic vertebrae and the posterodorsal module corresponds with the postdiaphragmatic vertebrae. The break between the anterodorsal and posterodorsal regions likely occurs at/near the diaphragmatic vertebra for all taxa. This pattern is also present using the alternative counting methods. These modules match a pre-post zygapophyseal pattern, which is a widely recognized regional division in the mammalian column, but they do not match the classic morphological regions-that is, cervicals, thoracics, and lumbar.
Using the thoracolumbar boundary count, both dorsal modules are separated by a set of weakly integrated vertebrae (T3-T04) forming a transitional zone between both modules. However, using the diaphragmatic start count, it is evident that this transitional zone is mainly composed by the diaphragmatic vertebra with varying location across taxa (TL10, TL11, or TL12). For all counting methods, this vertebra shows a weak integration with contiguous vertebrae. This pattern is not solely explained by a phylogenetic patterning, as the three highly integrated modules as well as the weakly integrated diaphragmatic vertebra are still significant with 2B-PLS phy (Fig. 4). Despite this, the strength of integration among pairs of vertebrae decreases when taking phylogeny into account (e.g., Martín-Serra and Benson 2020).
The cross-taxonomic analysis of integration performed here may suggest that the "regions" found by Jones et al. (2018b)cervical, pectoral, anterior dorsal, diaphragmatic, and lumbar regions-based on morphological patterns of individual taxa may be formed in several ways. For example, the "pectoral" region seems to be the result of a module overlap between the cervical and the anterodorsal regions, whereas the diaphragmatic "region" seems to be a gap between the anterodorsal and 2 5 8 EVOLUTION LETTERS JUNE 2021 posterodorsal modules. Therefore, different hierarchical combinations of modules/integration patterns may produce more complex regionalization patterns in terms of realized morphological and functional variation.

THE CERVICAL MODULE
The first two cervicals-that is, especially the atlas but also the axis in a lower degree-are weakly integrated with other cervical vertebrae (Fig. 4A). The PLS models exhibit that most of the EVOLUTION LETTERS JUNE 2021 2 5 9 coordinate shape changes are concentrated in the craniocaudal length of the body and the length of the spinous process (Fig. 5).
Our results evidence a single module for almost all the cervicals, from C02-C03 to C07 (Fig. 4A). However, the C06-C07 also belong to the "anterodorsal" module. Our data do not support the hypothesized C03-C05 mammalian developmental module of Buchholtz et al. (2012). In contrast, our data confirm the results obtained by Randau and Goswami (2017a) with a significant covariation between C02, C04, and C06. The C03-C05 developmental module could be basal for mammals and posteriorly disrupted in carnivorans. The main shape changes associated with this module are the same throughout the entire cervical region, suggesting that the entire C03-C07 represents probably a single module.
The high integration found among cervicals could be caused by developmental constraints, otherwise it would be more labile and the phylogenetic patterning would be higher. The common developmental "regionalization" of the mammalian neck results from the common Hox code for living placentals (Buchholtz et al. 2012;Böhmer 2017). Moreover, this high integration from C02-C07 could not be explained by the functional roles in which these vertebrae are involved, because although the mid-cervical vertebrae (C03-C05) are involved in axial rotation, the anterior and posterior cervicals do contribute to sagittal movements (Graf et al. 1995;Arnold 2020).
Although some degree of integration is observed between C02 and other dorsal vertebrae, especially the lumbars (e.g., TL16-TL20 in Fig. 4A), we have not found evidence in the presacral column of integration between the cervicals and all the thoracics plus lumbars, as demonstrated by Randau and Goswami (2017a) for felids and explained for the coincidence in ossification timing among these vertebrae. The strength of integration of cervicals-thoracics-lumbars is high, but this integration is very weak with the TL10-TL12 (Fig. 4A).
Both C06 and C07 are integrated with the anterodorsal module, but C06 is more integrated with the cervicals and C07 is more integrated with the thoracics (Fig. 4A). Moreover, the shape covariation between C07 and TL01 is different from the shape covariation between C07 and other cervicals, suggesting that the morphology of this vertebra is the result of a trade-off between variation in the cervical module and variation in the anterodorsal one. In any case, our data do not support a C06-T02 "pectoral" module, as proposed by Randau and Goswami (2017a) for felids because these two vertebrae are not more integrated within them than with C05 and TL01.

THE ANTERODORSAL MODULE
This module is composed by the prediaphragmatic vertebrae (Fig. 6), which are articulated with vertebrosternal ribs. Although our data do not provide a direct link between vertebral motion and integration, we hypothesize that this module may be related to motion constraints of the thorax. The shape changes associated with this module are very uniform and are related to the length of the spinous process (Fig. 5).
Contrary to other amniotes, mammals were able to separate breathing from locomotion early in their evolutionary history (Buchholtz 2012). This was greatly facilitated by the appearance of the diaphragm (Carrier 1987;Ruben et al. 1987;Perry et al. 2010). As a consequence, mammals were able to abandon costal breathing, which depends on exceedingly lateral movements of the axial system and the attached ribs (Carrier 1987). This probably happens as a response of enhancing the aerobic capacity (Ruben et al. 1987) due to the greater levels of sustainable activity (Carrier 1987). Ventilation in many nonmammalian amniotes occurs very occasionally at rest, and at a very low frequency at maximum exercise (Perry et al. 2010). The ventilatory demands of mammals are high, and they should keep high rates of ventilation when locomoting rapidly (Filler 1986;Carrier 1987). This higher breathing efficiency is acquired by strict excursion movements of the ribs coupled with the movements of the diaphragm (Filler 1986). The rigid array of the anterior thorax is related to prevent exceedingly bending movements of the axial system as it collapses rib pattern excursion and compromises ventilation (Filler 1986). This highly integrated anterodorsal module coincides with the rigid array of the anterior thorax with vertebrosternal ribs, and we hypothesize that this highly integrated region could be advantageous for impeding large increments of bending at each vertebra to secure respiration during locomotion and counteract loading forces transmitted from the limbs (Schilling and Hackert 2006). Moreover, Jones et al. (2020) using data from cat cadavers showed that the degree of ventroflexion, lateroflexion, and dorsiflexion ( = extension) is consistently reduced from the cervicals until the diaphragmatic vertebra.

THAT MARKS THE TRANSITION BETWEEN MODULES
Using the joined thoracolumbar and the thoracolumbar boundary counts, there is a set of three weakly integrated vertebrae (TL9-TL12 or T03-T05, respectively). However, this should be the result of the varying position of the diaphragmatic vertebra across taxa. The diaphragmatic start count and the selected verte-brae approaches confirm that the diaphragmatic vertebra ( Fig. 6) is the less integrated vertebra of the entire carnivoran axial system (Fig. 4C), separating both dorsal modules. The diaphragmatic vertebra is contained within the "Diaphragmatic joint complex (DJC)" of mammals (Slijper 1946;Filler 1986), which marks the transition from thoracic to lumbar forms across seriation. This region connects the fairly rigid array of anterodorsal vertebrae that are stiffened by elements of the costal articulation with the rib-less lumbar region. Filler (1986) claimed that the motion of the axial column is concentrated from this region of the posterior thorax to the lumbars, representing a region with more relaxed constraints of permissive motion compared to the anterior thorax. This has been recently demonstrated in a study based on maximum range of vertebral motion of cat cadavers (Jones et al. 2020), suggesting that ventroflexion and dorsiflexion of the column begin to increase substantially from the diaphragmatic vertebra (Jones et al. 2020).
The diaphragmatic vertebra is characterized by a change in zygapophyses orientation: prezygapophyses are located on the dorsal plane and postzygapophyses are located on the sagittal plane (Filler 1986). This region has reduced rib articulations, and is less constrained by the costal system, allowing ventroflexion and lateral flexion (Filler 1986). The diaphragmatic vertebra also reduces the neural spines, and it is followed by the presence of anticlinality in most cases, which also increases extension (Slijper 1946). However, Schilling and Hackert (2006) stated that in small therians during asymmetrical gaits sagittal bending neither corresponds to the trunk subdivisions (i.e., possession or lack of ribs) nor to the number of floating ribs or the position of the diaphragmatic vertebra. Despite this, their in vivo study found that "sagittal spine movements clearly occurred in the postdiaphragmatic region, but in none of the species was there a simple relation to the zygapophyseal facetation" (Schilling and Hackert 2006). Moreover, Oliver et al. (2016) using cadavers of the EVOLUTION LETTERS JUNE 2021 2 6 1 nine-banded armadillo (Dasypus novemcinctus) demonstrated that the postdiaphragmatic region has a high range of motion and low compliance.

THE POSTERODORSAL MODULE
The posterodorsal module comprises a set of highly integrated vertebrae being the last vertebra weakly integrated with other vertebrae of this module (Figs. 4A and 4D; Fig. 6). Filler (1986) proposes that the transverse processes of the lumbar vertebrae limit extension of the column, because the ligaments integrate the processes and provide passive limitation against their separation from each other by engaging the horizontal septum. Therefore, Filler (1986) concluded that the role of the lumbar vertebrae in locomotion is to prevent the exceeding extension propagated caudally from the vertebrae of the DJC. Jones et al. (2020) found that the degree of dorsiflexion and ventroflexion substantially decreases from the thoracolumbar boundary until the L06-L07 intervertebral joint, where a substantial increase in both parameters is noticed until the lumbosacral articulation where both reach their maximum values. Wachs et al. (2016) quantified intervertebral motions using X-ray reconstruction of moving morphology and scientific rotoscoping on dogs at walking and trotting gaits, and concluded that the greatest intervertebral motions occurred at the lumbosacral articulation or at L06/L07 intervertebral joints. Moreover, Wachs et al. (2016) stated that only a fraction of the maximum mobility of the lumbosacral joint is used during locomotion, which could be supported with the hypothesis of Filler (1986) that axial musculature and the transverse processes prevent exceeding extension.
Other authors have demonstrated that there is a strong correlation between the ecology of taxa and lumbar morphology (Randau et al. 2016;Jones et al. 2018b). Traditionally, although this sagittal mobility has been associated with the characteristic lumbar region without ribs (e.g., Boszczyk et al. 2001), Schilling and Hackert (2006) claimed that they are not perfectly correlated with the thoracolumbar transition. In any case, this is not contradictory with the high ecological signal found by Randau et al. (2016) for felids and by Jones et al. (2018b) for crown mammals. Indeed, different degrees of flexion-extension could explain the strong ecological signal across the lumbars. This variation in motion movements may be subject to natural selection toward specific locomotory demands, which could increase morphological disparity and evolutionary rates (Jones et al. 2018b).

ACKNOWLEDGMENTS
This study has been funded by the Spanish Ministry of Economy and Competitiveness (MINECO) (grant numbers CGL2015-68300P and CGL2017-92166-EXP to BF). We are grateful to F. de Paz (University of Valladolid, Spain), J. Pearson (University of Málaga, Spain), Z. Timons (Natural History Museum of Edinburgh, UK), and M. Bastir and P. Osborne (Museo Nacional de Ciencias Naturales, Madrid, Spain). We also thank the editor and two anonymous reviewers for their insightful comments.

Supporting Information
Additional supporting information may be found online in the Supporting Information section at the end of the article. Table S1. Species sampled in this study, their vertebral count and the position of the diaphragmatic vertebra (Diaph). Table S2. Number of landmarks digitized and the anatomical criteria per region used in this study. See also Figure 2. Table S3. Results of the Procrustes ANOVAs performed for each vertebra, according to the joined thoracic-lumbar count, with log-transformed centroid size as independent factor. Table S4. Results of the Procrustes ANOVAs performed for each vertebra, according to the thoracic-lumbar boundary count, with log-transformed centroid size as independent factor. Table S5. Results of the Procrustes ANOVAs performed for each vertebra, according to the diaphragmatic start count, with log-transformed centroid size as independent factor. Table S6. Results of the Procrustes ANOVAs performed for the selected vertebrae method, with log-transformed centroid size as independent factor. Table S7. Results of the significance tests of the phylogenetic and non-phylogenetic two-block PLS analyses according to the joined thoracic-lumbar count. Table S8. Results of the significance tests of the phylogenetic and non-phylogenetic two-block PLS analyses according to the count starting at the thoracic-lumbar boundary. Table S9. Results of the significance tests of the phylogenetic and non-phylogenetic two-block PLS analyses according to the count starting at the diaphragmatic vertebra. Table S10. Results of the significance tests of the phylogenetic and non-phylogenetic two-block PLS analyses with the selected relevant vertebrae.