Regional differentiation of felid vertebral column evolution: a study of 3D shape trajectories

Recent advances in geometric morphometrics provide improved techniques for extraction of biological information from shape and have greatly contributed to the study of ecomorphology and morphological evolution. However, the vertebral column remains an under-studied structure due in part to a concentration on skull and limb research, but most importantly because of the difficulties in analysing the shape of a structure composed of multiple articulating discrete units (i.e. vertebrae). Here, we have applied a variety of geometric morphometric analyses to three-dimensional landmarks collected on 19 presacral vertebrae to investigate the influence of potential ecological and functional drivers, such as size, locomotion and prey size specialisation, on regional morphology of the vertebral column in the mammalian family Felidae. In particular, we have here provided a novel application of a method—phenotypic trajectory analysis (PTA)—that allows for shape analysis of a contiguous sequence of vertebrae as functionally linked osteological structures. Our results showed that ecological factors influence the shape of the vertebral column heterogeneously and that distinct vertebral sections may be under different selection pressures. While anterior presacral vertebrae may either have evolved under stronger phylogenetic constraints or are ecologically conservative, posterior presacral vertebrae, specifically in the post-T10 region, show significant differentiation among ecomorphs. Additionally, our PTA results demonstrated that functional vertebral regions differ among felid ecomorphs mainly in the relative covariation of vertebral shape variables (i.e. direction of trajectories, rather than in trajectory size) and, therefore, that ecological divergence among felid species is reflected by morphological changes in vertebral column shape.


Introduction
From species description to detailed studies of ecomorphology, analyses of form have long been used by researchers examining ecological and evolutionary trends in both living and fossil organisms (e.g. Dumont et al. 2015;Lauder 1995;Rudwick 2005;Davies et al. 2007;Gonyea 1978;Gould 1966;Benoit 2010;Boszczyk et al. 2001;Goswami et al. 2012Goswami et al. , 2014. The geometric morphometrics revolution has greatly improved the scientific capacity to extract detailed information from biological structures. Yet it has also been hindered by computation issues with statistical tests used and the constraints involved in analysing data that are dense (e.g. large numbers of landmarks) and multidimensional, with specimen/landmark ratios decreasing as a result of these new advances Collyer et al. 2014;Adams 2014b;Cardini and Electronic supplementary material The online version of this article (doi:10.1007/s13127-016-0304-4) contains supplementary material, which is available to authorized users. Loy 2013). Newly developed software and methods are rapidly tackling these analytical power issues, with a plethora of recent papers describing and applying these approaches to diverse morphometric datasets (e.g. Adams and Collyer 2009;Adams 2014a, b;Collyer et al. 2014;Sheets and Zelditch 2013;Mitteroecker and Gunz 2009;Monteiro 2013;Polly et al. 2013;Mitteroecker et al. 2013;Klingenberg and Marugán-Lobón 2013).
Among morphological studies in the vertebrate literature, both those using geometric morphometrics (GMM) and studies using linear or cross-sectional measurements, there is a clear bias towards the morphology of the skull (e.g. Meachen-Samuels and Van Valkenburgh 2009a;Slater and Van Valkenburgh 2008;Fabre et al. 2014;Stayton 2005;Figueirido et al. 2010;Goswami and Polly 2010;Goswami 2006;Pierce et al. 2008Pierce et al. , 2009Piras et al. 2013;Drake and Klingenberg 2010;Foth et al. 2012;Meachen et al. 2014), followed by studies of the limbs (e.g. Bennett and Goswami 2011;Fabre et al. 2013;Bell et al. 2011;Alvarez et al. 2013;Martin-Serra et al. 2014;Adams and Nistri 2010;Walmsley et al. 2012;Zhang et al. 2012;Andersson and Werdelin 2003;Ercoli et al. 2012;Sears et al. 2013;Meachen-Samuels and Van Valkenburgh 2009b;Doube et al. 2009). The axial skeleton, in contrast, is comparatively underrepresented in the morphological literature, with the majority of work on this structure taking a biomechanical or developmental perspective (e.g. Macpherson and Fung 1998;Boszczyk et al. 2001;Long et al. 1997;Molnar et al. 2015;Smeathers 1981;Wellik 2007;Gál 1993;Müller et al. 2010;Buchholtz et al. 2012Buchholtz et al. , 2014Galis et al. 2014;Schilling and Long 2014;Narita and Kuratani 2005;Chen et al. 2005;Breit and Künzel 2004;Chatzigianni and Halazonetis 2009). Additionally, due to the difficulties in studying a structure that is composed of discrete units, research on axial skeletal morphology has frequently focused on separate analyses of individual vertebrae, with a few studies presenting intervertebral comparisons of individual measurements or differential morphospace occupation of vertebral types, rather than combined analysis of the full column (e.g. Alvarez et al. 2013;Jones 2015;Arnold et al. 2016;Manfreda et al. 2006;Buchholtz et al. 2014). Nevertheless, the limited morphometric studies of vertebral form have demonstrated that ecological specialisations and developmental patterning are reflected in the morphology of individual vertebrae, as well as along the entire spine (e.g. Jones and German 2014;Pierce et al. 2011;Shapiro 2007;Ward and Mehta 2014;Head and Polly 2015;Randau et al. 2016;Werneburg et al. 2015;Jones and Pierce 2015;Böhmer et al. 2015;Johnson et al. 1999;Chen et al. 2005). Indeed, many large clades, including the vast majority of placental mammals, do not display meristic changes (i.e. variation in number) in the presacral axial skeleton; therefore, adaptation of this structure must happen through modifications of its shape (Müller et al. 2010;Narita and Kuratani 2005;Buchholtz 2014;Buchholtz et al. 2012).
Recently, we conducted a large-scale linear morphometric analysis of the felid (cats) presacral vertebral column and found that this method was unable to strongly differentiate taxa based on either prey size specialisation or locomotor mode (Randau et al. 2016). For instance, there were few statistical differences in vertebral profile plots (i.e. variation in linear measures along the column), and a principal components analysis found a locomotory signal only in the lumbar region. These results were surprising considering felid prey size specialisation has been shown to correlate with osteological measures of the skull and appendicular skeleton (Meachen-Samuels and Van Valkenburgh 2009a, b;Slater and Van Valkenburgh 2008), and similar linear morphometric studies on other mammalian groups (e.g. pinnipeds, whales) have found the vertebral column to hold a strong ecological signal (e.g. Pierce et al. 2011;Buchholtz 2001a, b;Hua 2003;Finch and Freedman 1986). As felids are a morphologically conservative group, with little variation in musculoskeletal anatomy across the clade (Doube et al. 2009;Cuff et al. 2016a, b;Day and Jayne 2007), it remains uncertain whether the felid vertebral column holds little ecological signal or if linear morphometric techniques are not powerful enough to discriminate more subtle variation in vertebral form. To investigate this further, we extend our work by quantifying vertebral morphology in felids using three-dimensional landmarksbased GMM, and we include a novel application of phenotypic trajectory analysis (Adams and Collyer 2009;Collyer and Adams 2013) to identify ecological signal in serial structures. Three-dimensional (3D) landmarks are expected to provide greater detail and biological information than linear data (e.g. Fabre et al. 2014;Cardini and Loy 2013), and thus this work expands and improves upon existing linear studies considering this clade (Randau et al. 2016;Jones 2015). To our knowledge, two previous uses of 3D GMM to study the shape of a complete vertebral region have been reported in the literature (e.g. the cervical region, Werneburg 2015;Böhmer et al. 2015). While Böhmer et al. (2015) analysed individually landmarked cervical vertebrae by plotting them together with a principal component analyses, which described main shape variation among those and allows for qualitative analyses of shape change across taxa, Werneburg (2015) described a complex methodology that may not be broadly applicable. Specifically, that method relied on finding landmarks on three-dimensional reconstructions which had been matched to photographs of either manually articulated cervical vertebrae to approximate in vivo orientations or on model reconstructions of CT scans obtained from living animals. Those conditions are not readily available for many taxa, and thus we believe that the approach described here will be useful for a broader range of future studies. Additionally, Head and Polly (2015) used two-dimensional landmarks to characterise the precoaclal axial skeleton of squamates; however, the methodology described was applied to investigate patterns of regionalisation in the axial skeleton instead of testing correlations between shape and ecology.
We first analyse the individual shape of selected vertebrae and test for the influence of factors known to affect the shape of skull and limbs, including size, locomotion and prey size specialisation (Carbone et al. 1999;Meachen-Samuels and Van Valkenburgh 2009a, b). We then conduct separate analyses of each region of the vertebral column (cervical, thoracic and lumbar regions, and hypothesised functional regions composed of different combinations of these regions) and assess shape differences and differential allometry associated with ecological groupings. Finally, we apply phenotypic trajectory analysis to the main dataset, a combined analysis of cervical, thoracic and lumbar vertebrae, and also to individual regions with significant ecological signal, to analyse the shape of the vertebral column as a succession of contiguous units, thus overcoming the long-standing issue of analysing vertebrae as independent objects in geometric morphometric studies. We use these approaches to test the following hypotheses: (1) ecology is a significant influence on the morphology of felid vertebral column, and (2) vertebral regions display different levels of ecological and phylogenetic signal due to the regionalisation of shape in the mammalian vertebral column.

Data collection
In order to compose our 3D dataset, landmarks were collected from 19 presacral vertebrae from nine species of extant cats using an Immersion Microscribe G2X (Solution Technologies, Inc., Oella). This dataset included the following vertebrae: atlas, axis, C4, C6, C7, T1, T2, T4, T6, T8, T10, T11, T12, T13, L1, L2, L4, L6 and L7. As time constraints hindered the ability to collect dense data for every vertebra, but sufficient data were needed to describe the full presacral vertebral column morphology, the selection of these vertebrae was based on the following criteria: vertebrae with measurements that accounted for the highest principal component loadings in a previous linear study (Randau et al. 2016), vertebrae comprising the boundaries between vertebral regions and immediately preceding and succeeding vertebrae (e.g. C7 and T1, and C6 and T2, respectively), and vertebrae which are thought to be of particular biomechanical importance (e.g. T11, the anticlinal vertebra). Landmarks were collected from 109 specimens, ranging from seven to 17 specimens per species, with the final dataset including a total of 1712 individual vertebrae (see Table S1 for specimen numbers). Analyses grouped these dataset in various ways, ranging from treating all vertebrae individually to pooling vertebrae in the most inclusive grouping (C4-L7, excluding T11-T13), as described further below. Vertebrae were also grouped into the following five regions for some analyses, including C4-T10, T1-T10, T1-L7, T10-L7 and L1-L7. These regions were selected because they correspond to or group clear anatomical regions (e.g., T1-T10, L1-L7, and T1-L7) or more inclusive regions demarked by anatomical transitions (i.e. anterior or posterior vertebral column defined by the dorsal limit of the diaphragm, e.g. C4-T10 and T10-L7, respectively; Gray et al. 2005;Buchholtz et al. 2012;Jones 2015).
Sixteen homologous landmarks were identified on 14 of these vertebrae (i.e. the post-atlanto-axial and pre-sacral C4-L7 except for the T11-T13). Twelve landmarks were gathered on C1 (atlas), and 14 on C2 (axis), due to their unique morphologies ( Fig. 1 and Table S2 of landmarks). Vertebrae T11 to T13 lack transverse processes and thus two out of the 16 selected landmarks (i.e. the right and left transverse process tips) could not be identified on those elements. Comparative analyses across all sampled vertebrae require all observations to have the same landmarks. For this reason, the majority of the following analyses, unless otherwise stated, only used the 14 vertebral types that contained the same 16 landmarks ( Fig. 1d-i, i.e. not including the axis and atlas, shown on Fig. 1a-b and j-k, respectively, due to their unique shape, or vertebrae T11 to T13).
In order to still include the T11-T13 vertebrae in our tests of ecological correlates of axial skeleton morphology, we conducted a second analysis using two alternative landmarks that represent the locations of the right and left accessory processes of these vertebrae (Fig. S1, landmarks 7 and 8). Accessory processes are slender processes that originate on the pedicle and extend posteriorly, laterally to each postzygapophyses, and reinforce the interzygapophyseal joint (De Iuliis and Pulerà 2007). Additionally, accessory processes were also present on vertebrae L1, L2 and L4 of all species analysed here. Therefore, the second analysis used the two accessory process landmarks instead of transverse process landmarks for the vertebrae T11-L4, while the remaining vertebrae (C4-T10 and L6-L7) continued to use the transverse processes landmarks. In this manner, a dataset of 16 landmarks was constructed for 17 vertebrae, although two of these landmarks are not homologous in all of the vertebrae.
As only the 14-vertebrae dataset (excluding C1-C2 and T11-T13) was composed of homologous landmarks, we focus on the 'multi-vertebrae' analyses of that dataset, hereafter referred to as the 'homologous dataset' (or C4-L7 for shortening, although not containing T11-T13 as stated). The results from the alternative dataset that includes T11-T13 by using two non-homologous landmarks (accessory processes landmarks instead of transverse process landmarks for T11-L4), hereafter referred to as the 'alternative dataset', were remarkably consistent and are presented in the Supplementary information.
Ecological data for all analyses were collated from the literature (Meachen-Samuels and Van Valkenburgh 2009a, b;Sunquist and Sunquist 2002). Prey size groupings include small, mixed and large prey specialists. Locomotory groupings include arboreal, cursorial, scansorial and terrestrial. Phylogenetic comparative analyses used the composite tree of Piras et al. (2013) pruned to the species sampled here.
Prior to all subsequent analyses, missing landmarks due to broken specimens were imputed using the multivariate regression ('Reg') method in the 'estimate.missing' function of 'geomorph'. This approach predicts the missing landmarks by using a multivariate regression of the specimen with missing values on all other landmarks in the set of complete specimens ). A total of 126 out of 30,695 (0.41 %) landmarks were imputed. All vertebrae were then subjected to Procrustes Superimposition within the relevant sample (i.e. either within same vertebral type sample or specific vertebral region analysed depending on the analysis level) to remove any effects due to scale, rotation and translation.
Phylogenetic and ecological signal of individual and regional vertebral shape Preliminary analysis of vertebral column shape was performed with a combined principal component analysis (PCA) of all of the vertebrae in the homologous landmark dataset (C4-L7, excluding T11-T13). A second PCA was performed on the region encompassing vertebrae T10-L7 in the homologous landmark dataset. Scans of individual cheetah (Acinonyx jubatus, USNM 520539) vertebrae were used to create an average reference mesh with the 'warpRefMesh' function in geomorph, and this mesh was used to warp the PC1 and PC2 minimum and maximum shapes in order to display vertebral shape changes across the main eigenvectors.
The effects of centroid size and ecological specialisation (both in terms of locomotion and prey size categories) on vertebral shape were evaluated with factorial MANOVAs of the vertebral Procrustes coordinates (i.e. shape − centroid size × ecology). Factorial MANOVAs with this size-ecology interaction accounts for the effect of 'size' while examining the other factors that describe shape and define the groups. Additionally, these non-parametric MANOVAs with 'RRPP' (residual randomisation permutation procedure) allowed for significance tests with multidimensional data that have fewer observations than dimensions (Collyer et al. 2014). These analyses were performed separately on each vertebra from  Table S2 C1 to L7, with each set composed of an across-species pool (i.e. C1 dataset contained all C1 vertebrae measured, across all nine species) as well as on the complete homologous dataset (see Supplementary information for further details on analyses of the alternative dataset). Additionally, factorial MANOVAs were applied to the five vertebral regions of described above, using the homologous dataset. Each described region contained all vertebrae of the named types, including all species listed here.
In order to assess the influence of phylogenetic relatedness on vertebral shape and centroid size (i.e. whether more closely related species were more phenotypically similar ;Felsenstein 1985), we first constructed the mean shape for each individual vertebra (C1 to L7) per species and calculated the phylogenetic signal with the 'Kmult' method (i.e. a multivariate version of the K-statistic; Adams 2014a) with the 'physignal' function in 'geomorph'. As L1-L4 have both transverse processes and accessory processes and thus are the only elements with different landmarks in the homologous and alternative datasets, this analysis was performed for both datasets for those elements. For individual vertebrae that presented a significant phylogenetic signal in their shape across the studied species, we also performed phylogenetic MANOVAs to assess the relationship between shape, centroid size and ecological factors. Phylogenetic MANOVAs use a phylogeny-informed context under a Brownian motion model of evolution to calculate a phylogenetic transformation matrix and the Gowercentred distance matrix from predicted variable values, which are then used to asses significance from comparisons between the values of statistical attributes obtained from those and the observed values (Adams 2014b;Garland et al. 1993). Phylogenetic MANOVAs were done using the 'procD.pgls' function in 'geomorph'.

The interaction of allometry and ecology in vertebral regions
Considering that previous studies of felid vertebral morphology have demonstrated the widespread influence of allometry in vertebral linear dimensions (see below; Randau et al. 2016;Jones 2015;Jones and Pierce 2015), we investigated whether prey size or locomotory ecomorphs presented different allometries in their vertebral shape. Based on the MANOVA results (see below and Table 5), the vertebral region with the highest absolute variance explained by the two ecological variables (i.e. T10-L7) was selected to examine differences in vertebral allometry with respect to ecological specialisation.
Using the 'PredLine' method of the 'plotAllometry' function in 'geomorph', the predicted allometric scores for these regions were calculated for each ecological group from the shape against centroid size regression. The method used produced allometric trajectories (i.e. plotted PC1 of the predicted values against size) which clearly exhibited allometric differences between ecological groups (Adams and Nistri 2010).
The significance of the differences in the log centroid sizeshape relationship between groups could be quantified by both the P value of the comparisons between slope distances, which itself measures differences in amount of shape change per unit of centroid size change, and the slope angle's P value, which indicates if the directions of these vectors point at different regions of the morphospace (Collyer et al. 2014;Collyer and Adams 2013). This last step was performed using the 'advanced.procD.lm' function in 'geomorph'.

Ecological signal across the vertebral column
Shape for the proxy of an entire vertebral column (i.e. C4-L7, excluding T11-T13), as well as for individual regions, was quantified using a novel application of phenotypic trajectory analysis (PTA). PTA identifies a shape trajectory among associated data points (vertebrae, in this case) and then compares this trajectory among vertebra within each predetermined group (e.g. mean shape of C7 for all arboreal taxa), and then traces the trajectory between these means (e.g. C6 to C7, C7 to T1, etc.) Collyer 2007, 2009;Collyer and Adams 2013). The trajectories can then be visualised in morphospace for a qualitative comparison between groupings, and differences in size, direction and shape of the trajectories for each group can also be quantitatively compared. As above, taxa were grouped by prey size and locomotory categories for analysis of ecological signal in phenotypic trajectories.
As noted in 'Materials and methods' section, all of the following results refer to the homologous dataset unless otherwise indicated. The PC1 minimum shape was generally mediolaterally and anteroposteriorly compressed and dorsoventrally elongated, with smaller centrum width and centrum length, smaller distances between transverse processes, prezygapophyses and post-zygapophyses, and larger heights for the centrum, neural canal and neural spine. The PC1 maximum shape showed larger centrum width and centrum length, larger distances between transverse processes and intra-zygapophyses, but shorter heights for the centrum, neural canal and neural spine. PC2, which separated the C4 cluster from the other two vertebral clusters, presented similar shape differences, with the PC2 minimum shape displaying even more exaggerated features related to mediolateral compression, but, in contrast, also exhibiting some anteroposterior elongation.
The main feature of PC2's maximum shape was the relative augmentation of the distances in the mediolateral dimension, with larger centrum width and intra-zygapophyseal distances. Results from the PCA applied to the 'T10-L7' region (Table 2  and Table S5, see below) showed that the majority of the variation (>90 %) was explained by the first five PCs, with PC1 explaining >60 % of total variance. When individual vertebral datasets were subjected to factorial MANOVAs of shape against centroid size, locomotion and prey size groups (Table 3), all vertebrae displayed significant correlations of shape with all three factors (P < 0.001-0.05), with the exception of the T8 × prey size (P > 0.05). After Bonferroni correction, only three correlations ceased from being significant (i.e. P > 0.003): C6 and T10 versus prey size, and L7 versus centroid size. The three examined factors explained a range between 3 and 23.77 % of vertebral shape (highlighted on Table 3). Further, estimating the influence of evolutionary relatedness on vertebral shape recovered a significant (i.e. P < 0.05) phylogenetic signal for the mean shape (i.e. Procrustes coordinates) of only five vertebrae: atlas, axis, C6, T1 and T2 (Table 4); however, after Bonferroni correction, this signal was only significant for the atlas and axis (i.e. P < 0.003). Conservatively, all of these five vertebrae were further subjected to a second round of MANOVAs using the same factors as above, while controlling for this phylogenetic signal. After this correction, none of ecological correlations were significant (P >> 0.05, Table 5). No phylogenetic signal was recovered for centroid size of any of the analysed vertebrae.
Factorial MANOVAs were also applied to five regions composed of multiple vertebrae for quantification of the influence of ecological factors on vertebral regions. The highest ecological signal in vertebral shape was observed in the region from T10 to L7, with ∼17.55 and~12.2 % of overall shape explained by prey size and locomotory categories, respectively (see MANOVAs in Table 6 for all results). This region also displayed the second highest values for the influence of centroid size on shape (∼7.8 % Table 6). No significant correlation with locomotory categories was found for the complete homologous dataset (C4-L7) or for the C4-T10 region, while significant (i.e. both prior and after Bonferroni correction) correlations with both locomotory and prey size groups were found for the other regions but those ranged between 2.0 and 11.9 % for locomotion and 1.6-12.6 % for prey size ( Table 6).

The interaction of allometry and ecology in vertebral regions
As stated above, the interaction factor between ecological groups and centroid size was significant and exhibited its highest values (Table 6) for the T10-L7 region, demonstrating that species belonging to different ecological groups displayed For each factor (i.e. centroid size, locomotion and prey size), the highest coefficient of determination (R 2 ) value is shown in bold, and the lowest value is displayed in italics. The sole test which was not statistically significant (i.e. P > 0.05) is underlined. Tests which are not significant after Bonferroni correction (i.e. >0.003) are marked with an asterisk distinct shape versus size relationships in the posterior presacral vertebrae. Plots of the predicted allometric trajectories for each ecological factor on both datasets are presented in Fig. 3a and b. The analysis using prey size groups for categorisation showed that, while 'small' and 'big' prey size groups possessed allometric trajectories that were very similar in slope distance (P > 0.1, Table 7), the 'mixed' prey size group's trajectory exhibited a slope distance that was significantly different from both the large and small prey size groups (P << 0.05). However, differences in the slope distance of the allometric trajectories between 'large' and 'mixed' prey size groups were not significant after Bonferroni correction (i.e. P > 0.006). Slope angles were significantly different between the 'large' and 'small prey' categories, but not after Bonferroni correction. Grouping species by their locomotory modes resulted in allometric trajectories that were similar in slope distance between 'arboreal' and 'cursorial' groups (P >> 0.05), but both differed in all other pairwise comparisons between locomotory groups (P << 0.05). Slope angles were only significantly different between the 'terrestrial' and 'scansorial' subsets (P << 0.05).

Ecological signal across the vertebral column
Phenotypic trajectory analysis was first performed using the most inclusive homologous dataset (i.e. C4-L7) to quantify the shape of the post-atlanto-axial presacral vertebral column (Table 8 and Fig. 4), followed by analysis of the T10-L7 region. When species were grouped by prey size specialisation, phenotypic trajectories for the full dataset were significantly different in shape. The 'small' prey size trajectory was also different from both the 'mixed' and 'big' prey size groups in terms of trajectory size. Grouping species by locomotory mode with the complete dataset was not performed because the MANOVA results for this region exhibited a nonsignificant correlation with locomotory groups (P >> 0.05, Table 6). Analysis of the T10-L7 vertebrae resulted in significant differences in phenotypic trajectories for both ecological factors (Table 9 and Fig. 5a, b). With prey size categorisation, the phenotypic trajectories were all significantly different in direction. The 'small' prey size trajectory was also different from both the 'mixed' and 'big' prey size groups in terms of shape. Locomotory group trajectories were different in direction for all pairwise comparisons, except between the 'scansorial' and 'terrestrial' groups. In terms of shape, the 'cursorial' phenotypic trajectory was statistically different from the 'arboreal' and 'scansorial' trajectories, but only before Bonferroni correction and not after (P <0.05 but >0.006, respectively).  The highest coefficient of determination (R 2 ) values for both prey size and locomotion were found in the T10-L7 region and are shown in bold.
The tests which were not statistically significant (i.e. P > 0.05) are underlined. All significant tests were still significant after Bonferroni correction (i.e. P < 0.008)

Discussion
When combined, analyses of the relationship among 3D vertebral shape, size, ecology and phylogeny provide a more complete understanding of the forces shaping the evolution of the felid vertebral column evolution. The results reported here have confirmed our initial hypotheses on ecological drivers in the vertebral column shape differentiation in felids, and we have detailed how specialisation towards the observed ecologies correlates with regionalisation of the presacral axial skeleton. While vertebrae in the anterior-most region of the felids' vertebral columns (i.e. atlas and axis, but also C6, T1 and T2) were more phylogenetically conservative in shape, the posterior regions of the vertebral column showed a stronger influence of ecological specialisations. That the strongest size and ecology correlations are observed in this more caudal region of the presacral vertebral column (i.e. T10-L7; see Supplementary information for similar results on the dataset using the accessory processes landmarks) supports the inference that this region may be subjected to stronger selection, or equally to weaker evolutionary constraints, and might present greater evolutionary respondability across felids, or even more broadly. This observation agrees with the work by Jones and German (2014), in which they found that in mammals, centrum length varied the most in the lumbar region both through ontogeny and interspecifically. As an osteological measurement that is informative towards the degree of passive robustness at intervertebral joints (Pierce et al. 2011;Shapiro 1995Shapiro , 2007Koob and Long 2000), centrum length can be used to Fig. 3 Allometric trajectories displaying the differences in the predicted shape-size relationship between ecological groups. a Species groups by their prey size, b species grouped by locomotory category  make inferential comparisons of resistance to intervertebral bending and general biomechanical properties between species or ecological groups. An additional PCA limited to the T10-L7 vertebrae (post-diaphragmatic homologous dataset) (Fig. 2c) shows that the anteroposterior vertebral axis, which primarily represents centrum length, is one of the main contributors to variation in this dataset.
When compared to our previous work on the linear morphological change in the felid axial skeleton (Randau et al. 2016), our present study supports our general conclusions of regionalisation of ecological signal in the vertebral column, with stronger locomotory signal present in the posterior region. However, contrary to results from linear data (Randau et al. 2016), the 3D analyses described here also found a significant correlation between vertebral morphology and prey size specialisation. Previous studies of individual vertebral attributes (e.g. centrum length) and different proxies for body size (e.g. total vertebral length, body mass) using length measurements have also identified significant allometry across felids (Randau et al. 2016;Jones 2015). Here, we were interested in investigating whether the influence of size (i.e. centroid size) on vertebral multidimensional shape was also regionalised and, most importantly, whether such scaling relationships differed with ecology. Our results reinforce the conclusion that size influences vertebral shape throughout the axial skeleton (i.e. C4 and post-T2 vertebrae), but that these size effects are strongest in T10 and the lumbars (Tables 3 and 6, and in the last thoracics in Table S6). Additionally, we have demonstrated that ecological specialists, especially in terms of locomotory specialisation, indeed exhibit a distinct scaling relationship between shape and centroid size (Table 7). Observed differences between prey size subsets were very consistent with both measures of differentiation (slope angle and distance). 'Small' and 'mixed' prey size groups were shown to have distinct allometric vertebral shapes. Although 'large' and 'small' prey groups were not significantly different in terms of the intensity of their allometries (i.e. the Procrustes distances between slopes), they displayed distinct angles in their slope vector, showing that the covariances between the variables are different in these ecological categories (Collyer and Adams 2013;Adams and Collyer 2009). However, these differences between 'large' and 'small' categories, or regarding the intensity of the allometry between 'large' and 'mixed' categories, were not significant after correction, suggesting differences in allometry between prey size specialist groups might be subtle. This could therefore be one of the factors which caused linear measurements were not to be successful in finding correlations between felid vertebral morphology and specialisation towards prey size (Randau et al. 2016). With regards to locomotory specialisation, the two statistical attributes presented different patterns. A better separation between the groups was found in terms of the intensity of their allometries than in their directions. Additionally, it is clear from the observation of regression slopes (Fig. 3b) that allometric shape changes are much greater in 'arboreal' and 'cursorial' species and, although significant, size-related changes in the posterior vertebral morphology are less demarked in 'scansorial' and 'terrestrial' felids. Although all but one pairwise comparisons were significantly different with regards to slope distance, the only  Statistically significant values (i.e. P < 0.05) are shown in bold. Pairwise comparisons which were not significant after Bonferroni correction (i.e. P > 0.006) are marked with an asterisk significant difference in the direction of the allometric trajectories was found between the 'terrestrial' and 'scansorial' categories. Hence, although these two more generalist locomotory groups show a comparatively smaller degree of vertebral allometric scaling, they are still distinct in the relative way size influence vertebral shape variables. As nearly all individual vertebrae showed some significant correlation between shape and ecology (i.e. Table 3), individual analyses alone provide little clarity in terms of regionalisation of ecological and phylogenetic signals. Such differentiation was only possible when sets of vertebrae were analysed together through PTA. With this method, we were able to quantitatively differentiate the vertebral shape gradient changes between locomotor and prey size specialist felid species, therefore extracting the subtle morphological changes between the recognised ecomorphs in this phenotypically conserved clade.
Of the two ecological factors examined in this study, only prey size specialisation as an isolated factor exhibited a significant correlation with total vertebral column shape, contrary to the results of linear analyses (Randau et al. 2016). This result once again supports the regionalisation of locomotory specialisation in the vertebral column, which was instead found to significantly correlate only to more posterior regions, while also highlighting the increased resolution provided by 3D data. However, because prey size specialisation is directly correlated to the species' body mass (Carbone et al. 1999(Carbone et al. , 2007, a significant correlation between this factor and vertebral shape is possibly an indirect reflection of overall body size influence on vertebral three-dimensional shape. When we focused our analyses on the vertebral regions with highest correlations between shape and the factors examined, the T10-L7 trajectories were best able to separate among ecological groups, both for the locomotion and prey size categories (Fig. 5a, b). All significant differences between trajectories were found in comparisons of the shape and direction of those trajectories (Table 9). This result suggests that no differences in the amount of shape variation (i.e. trajectory size) were found in the species of felids studied here. Additionally, this differentiation in trajectory direction implies that the differences found were primarily based on the distinct relative covariations of vertebral shape variables between ecological groups throughout the vertebral column (Collyer and Adams 2013;Adams and Collyer 2009). More interestingly put, these differences in trajectory direction between groups are evidence of ecological divergence between those groups Stayton 2006). As it follows, the only two groups that did not differ significantly in trajectory direction (the 'scansorial' and 'terrestrial' groups) show ecological convergence in the shape of the posterior vertebral column.
Combining the PTA and posterior region PCA results (Fig. 2c) provides additional information on the changes in vertebral morphology correlated with cursoriality in felids. Cheetahs (Acinonyx jubatus), as the species represented by the 'cursorial' locomotory group, presented an average lumbar morphology that exhibited longer centra, and overall less shortening of the centrum from L1 to L7, which could be visualised by the trajectory lumbar points presenting lower values on PC1 and higher values on PC2 (Fig. 5b). The relative length of centra has been shown to be associated with the degree of flexibility between two consecutive vertebrae (Koob and Long 2000;Long et al. 1997;Pierce et al. 2011), and results from a study by Jones (2015) on linear vertebral dimensions revealed allometric shortening of the lumbar region in felids (but see Randau et al. 2016 for alternative results showing isometric scaling of the lumbar region in this family, albeit with a different sample). Ergo, having lumbar vertebrae that are relatively longer might indeed contribute to greater sagittal bending and contribute to having the longer stride lengths observed in this highly specialised felid (Hildebrand 1959).

Conclusion
The vertebral column has been underrepresented in the functional morphology and morphometric literature, but recent studies have shown that the vertebral form carries rich developmental and ecomorphological signals. Here, through multivariate statistical analyses, we have demonstrated that the use of geometric morphometrics to study the axial skeleton can offer even more detailed ecomorphological information than what has been reported by linear studies. Additionally, we have here provided the first application of a method that allows for the shape analysis of a contiguous sequence of vertebrae as functionally linked osteological structures.
We have shown that ecological correlates influence the shape of the vertebral column heterogeneously, specifically with discrete regions such as the posterior axial skeleton presenting higher correlation with both locomotory and prey size specialisation. Furthermore, we suggest that the post-T10 vertebrae may be the most ecologically adaptable region among felid species. While anterior vertebrae may either have evolved under stronger phylogenetic constraints or are more ecologically conservative, posterior vertebrae show clearer differentiation between ecomorphs in Felidae.
Future studies, which may benefit from focusing on a more restricted species range or on smaller vertebral regions, would gain from including vertebrae that were not analysed here in order to compare the general patterns found to specific complete regional trends.