Unravelling intravertebral integration, modularity and disparity in Felidae (Mammalia)

Morphological integration and modularity, which describe the relationships among morphological attributes and reflect genetic, developmental, and functional interactions, have been hypothesized to be major influences on trait responses to selection and thus morphological evolution. The mammalian presacral vertebral column shows little variation in vertebral count and therefore specialization for function occurs primarily through modification of vertebral shape. However, vertebral shape has been suggested to be under strong control from developmental canalization, although this has never been explicitly tested. Here, we assess hypotheses of developmental modules in the vertebrae of felids to determine whether developmental interactions are a primary influence on vertebral modularity. Additionally, we analyze the magnitudes of both intravertebral integration and disparity to evaluate if level of integration varies along the vertebral column and, if so, whether integration and disparity are associated. Our results confirm the hypothesis of vertebral developmental modularity, with most presacral vertebrae displaying two modules. Exceptions are concentrated in the boundaries among traditional and functional regions, suggesting that intravertebral modularity may reflect larger‐scale modularity of the felid vertebral column. We further demonstrate that overall integration and disparity are highest in posterior vertebrae, thus providing an empirical example of integration potentially promoting greater morphological responses to selection.


INTRODUCTION
The dichotomy between maximum individual trait adaptation and cohesion between functioning parts is one that directly affects phenotypic response to selection Badyaev et al. 2005;Hansen and Houle 2008;Porto et al. 2009;Goswami and Polly 2010a;Goswami et al. 2014). The basis for understanding how organisms are organized was laid by the seminal work by Olson and Miller (1958) in which they described the fundamental concepts of phenotypic integration and modularity as can be ascertained through quantification of patterns of trait covariation. In line with these, modules are a set of traits that show higher covariation among them than with other parts of the organism due to shared genetic or developmental origins or function, while integration is the overall pattern of intercorrelation (e.g., Hansen and Houle 2008;Klingenberg 2008;Goswami and Polly 2010b; Marug an-Lob on 2013). Interestingly, however, those two definitions are not contradictory and complex traits may present overall high integration and still be modular (Bookstein 2015), such as the mammalian skull (Goswami 2006a,b;Goswami and Polly 2010c). Specifically, trait units might present significant covariation among the whole structure (i.e., integration), while still showing higher organization into smaller sets which present consistently higher within-set covariation than across the whole phenotype. Moreover, trait integration has been shown to reflect shared developmental pathways in early ontogeny, postnatal function, and heterochronic shifts (Zelditch and Carmichael 1989;Goswami et al. 2009;Zelditch et al. 2009;Bennett and Goswami 2011;Goswami et al. 2012Goswami et al. , 2014, and to be susceptible to reorganization by extreme changes in selection (Drake and Klingenberg 2010).

Intravertebral developmental modularity
Morphological traits in fully grown organisms may present correlations due to developmental modularity, by which variation in a set of traits is dependent on a common embryonic origin or other shared developmental history (Cheverud 1996;Arthur 1997Arthur , 2002Raff and Sly 2000;Klingenberg 2003;Buchholtz et al. 2012). However, trying to infer developmental modularity through the organization of adult morphology can be problematic due to repatterning of integration through ontogeny (Hallgrimsson et al. 2009). Nevertheless, knowledge of trait developmental origin can be used in confirmatory analyses to test hypotheses of developmental modularity West-Eberhard 2003;Klingenberg 2013). Here, we test the hypothesis that developmental origins of vertebral components (i.e., centrum vs. neural spine attributes) dictate adult vertebral morphology in cats (i.e., Felidae, Mammalia).
Mammalian vertebral column development has been suggested to be under strong canalization and developmental stability (Galis 1999;Narita and Kuratani 2005;Wellik 2007;Buchholtz and Stepien 2009;Hautier et al. 2010;M€ uller et al. 2010;Asher et al. 2011;Varela-Lasheras et al. 2011;Buchholtz et al. 2012;Fleming et al. 2015), and the derivation of somitic segments into the tissues involved into limb and vertebral column formation has been described in great detail (Christ et al. 2007).
For mammals, in which presacral vertebral count shows very little variation when compared to other vertebrate clades (Narita and Kuratani 2005;M€ uller et al. 2010;Buchholtz et al. 2012;Buchholtz 2014), changes in the axial skeleton are typically manifested in changes in vertebral shape. Buchholtz (2007) summarized the types of evolutionary change that have been observed in vertebral column morphology; those concerning changes in the mammalian axial skeletal morphology may reflect "diversifying" or "skeletogenetic" changes (caused by effects of Hox genes and growth factors) or be due to changes in "module association" of these vertebrae (Raff 1996;Polly et al. 2001;Buchholtz 2007).
Additionally, Christ et al. (2007) have described how vertebral components are derived from distinct somitic origins through segmentation of the sclerotome. The vertebral body (centrum) originates from the ventral and, to a lesser degree, central regions of the sclerotome, while the neural arch, spinous process, pedicles, and transverse processes originate from the dorsal and posterior central regions of the sclerotome and integrated somitocoel cells (Fig. 1). The condensation of these two vertebral parts has also been shown to be distinct, with the centrum-related sclerotome condensing around the notochord, while the same is not true and not yet fully understood for the development of the other vertebral elements (Hall 1977;Christ et al. 2000Christ et al. , 2007. Additionally, Boyd (1976) has confirmed that all presacral vertebrae in cats originate from two ossification centres, the only exception being C2 (axis) with a third ossification centre for the dens.

Trait integration can direct responses to selection
Research in the last few decades has built on the work on integration and modularity by demonstrating how trait relationships can both shape responses to selection and be affected by extrinsic perturbations such as environmental stress (West-Eberhard 1989;Badyaev and Foresman 2004;Hansen and Houle 2008;Badyaev 2010;Goswami and Polly 2010a,b;Buchholtz et al. 2012 Goswami et al. 2014Goswami et al. , 2015. Some of the direct ways integration and modularity have been suggested to affect trait evolution are by either constraining or promoting the spectrum of responses to selection (Cheverud 1996;Hansen and Houle 2008;Marroig et al. 2009;Klingenberg and Marug an-Lob on 2013;Sears et al. 2013). Integration has been traditionally hypothesized to constrain these responses to a smaller portion of the morphospace because high correlation among traits means that any change in the trait directly affected by selection can be hindered by stabilising selection on other covarying traits. Similarly, modularity has been hypothesized to counter this effect, by breaking larger sets of correlated traits into smaller modules, allowing newly independent modules to respond more freely (i.e., potentially promoting larger phenotypic variation). However, Goswami et al. (2014) demonstrated through the use of simulation analyses that integration may promote both lower and higher degrees of morphological disparity, and that range in disparity can be considerably larger in correlated traits than in uncorrelated ones, confirming previous hypotheses on the possible effects of integration (Schluter 1996;Klingenberg 2005). By directing variation along particular axes of the total possible morphospace, the maximum range of variation can be increased (Schluter 1996;Goswami et al. 2014).
Here, we first test the hypothesis that developmental origin drives intravertebral modularity, resulting in two intravertebral modules in adult morphology: the centrum and the neural spine. We subsequently quantify the magnitude of overall integration in individual vertebrae of felids by measuring relative eigenvalue standard deviation (Pavlicev et al. 2009) and compare these results to vertebral morphological disparity to determine whether higher integration is associated with higher or lower disparity. We conduct these analyses in the presacral vertebral column of felids and discuss our results in relation to previous analyses of ecological specialization in felid vertebral morphology (Randau et al. 2016b) and previous studies of the evolutionary significance of phenotypic integration and modularity.
In order to capture the most detail in vertebral morphology, and due to morphological differences throughout the vertebral column, different sets of landmarks were collected in some vertebral regions: 12 landmarks were gathered on C1 (atlas), 14 on C2 (axis), 18 on C4, 20 on C6, 16 on C7-T10, 16 on T11, 17 on T12-T13, 19 on L1-L4, and 17 on L6-L7 (see Table S2 for landmarks identity). Additionally, in order to facilitate direct comparisons across as many vertebrae as possible, 16 landmarks are homologous in C4-T10 and L1-L7, and thus only these landmarks and vertebrae were used in analyses of disparity and an additional analysis of integration for direct comparison to the disparity results (described in the Data analysis section). Vertebrae C1, C2, and T11-T13 were not included in the disparity analysis due to their unique morphology (e.g., vertebrae T11-T13 lack transverse processes but present accessory processes).

Intravertebral modularity
Vertebra-specific landmark coordinates for C1-L7 were assigned to modules based on models of developmental origins and ossification centers (Table 1). All vertebrae were hypothesized to be composed of two developmental modules: "centrum" and "neural spine" (Christ et al. 2007) as depicted in Figure 1. Additionally, a three-module hypothesis was also tested for C2 (axis) as the dens has been shown to originate from an additional ossification center that fuses with the centrum early in vertebral development (Boyd 1976).
The degree of modularity and the significance of these models were evaluated by using two alternative methods: RV coefficient analysis (Escoufier 1973;Klingenberg 2009) and Covariance Ratio analysis (CR; Adams 2016). Both methods are similar in their outputs, but differ in that CR disregards within-  Table S2 for landmark definitions. trait variation and uses only the covariation between and among traits for its calculations, while RV accounts for both measures.
We have chosen to present both results because, while RV has been one of the most used confirmatory analyses of modularity in recent years (Goswami and Polly 2010a;Klingenberg 2009), it has recently been shown to be sensitive to sample size and landmark number (Fruciano et al. 2013;Adams 2016).
Significance of the hypothesis of modularity in both methods is obtained by randomly assigning landmarks to 10,000 alternative models of modularity to generate a distribution of values. Significant results are indicated if the observed signal is small (here, P < 0.05) relative to the randomly generated distribution.

Accounting for phylogenetic relationships
Modularity results prior to any phylogenetic correction to vertebral shape or analyses are displayed due to the following reasons: (i) The mammalian vertebral column has been suggested to be under strong developmental control and, especially with the Felidae family being very constrained in count, there is no reason to assume that individual felid species should present distinct developmental pathways to vertebral formation (Narita and Kuratani 2005;M€ uller et al. 2010;Buchholtz et al. 2012); (ii) Removal of any potential phylogenetic signal on shape may conceal real patterns of morphological modularity or integration driven by genetic or developmental origins ); (iii) Tests for phylogenetic signal in shape were significant for only two anterior vertebrae in felids, the atlas, and the axis (Randau et al. 2016a), while tests for phylogenetic signal in both shape and centroid size of all other studied vertebrae were not significant. Instead, we corrected for grouping multiple species into a single analysis by first calculating a pooled within-species variance-covariance matrix (VCV) for each vertebrae and then used this VCV matrix in CR analysis of vertebral modularity. This pooled within-species VCV matrix was calculated using the "covW" function in the "Morpho" package (Schlager 2016) in R. It is important to raise the caveat that this is a new implementation of the CR method, one which has not yet been tested through the use of simulations, and therefore caution should be kept in mind when applying this methodology to other studies. Nevertheless, modularity results both with and without using the pooled within-species VCV matrix were similar and therefore there is no obvious reason to think that the properties of the CR method would not hold in this case.

Overall vertebral integration and disparity
Vertebrae C4-L7 (excluding T11-T13) containing the 16 homologous landmarks were individually subjected to a General Procrustes Superimposition for extraction of shape coordinates (i.e., excluding information on size, rotation, and translation). The correlation matrix was obtained from these shape coordinates, and this was subsequently used to calculate the singular-value decomposition to generate matrix eigenvalues.
The overall morphological integration per vertebra was calculated using relative eigenvalue standard deviation (i.e., eigenvalue dispersion) as detailed by Pavlicev et al. (2009). High numbers of eigenvalue dispersion indicate strong integration, as variance is concentrated on fewer eigenvectors due to high covariance of traits, at the cost of low variance explained by higher eigenvectors. This measure of integration has been shown to be highly correlated with r 2 (mean squared correlation coefficient; Marroig et al. 2009;Goswami et al. 2014; not to be confused with the coefficient of determination R 2 ), and to be independent of trait number, and thus can be readily used for comparison across datasets. Therefore, we also calculated this measure using the specific vertebral landmark datasets for C1, C2, C4, C6, T11, T12-T13, L1-L4, and L6-L7 for maximum shape information, after subjecting individual vertebrae to General Procrustes Superimposition.
Morphological disparity per vertebra (e.g., T1) was calculated on the C4-L7 shape coordinates (homologous landmarks) both as Procrustes variances and as maximum Procrustes distance between specimens (Zelditch et al. 2012). The Procrustes variance analysis was performed both with and without centroid size as a covariate, as vertebral size has been shown to correlate with shape throughout the spine (Randau et al. 2016a,b). Both measures of disparity were calculated first per individual species per vertebra, and then across taxa per vertebra, using the species mean shapes.

Intravertebral modularity
Results from both RV and CR analyses of modularity were consistent in all but one case, and strongly supported the twomodule model (P < 0.01) for all but six (C2, C7, T1, T8, L6, and L7) of the 19 analyzed vertebrae. They differed only with regards to T13, which was marginally significant for the tested modules with RV analysis, but significant when analysed with CR (P values ¼ 0.051 and 0.011, respectively; Table 2). The threemodule model tested for C2 was not supported (P > 0.05). When testing the modularity model using the pooled within-species VCV matrix, three vertebrae presented different results: the threemodule model was supported for the axis, and the two-module model was significantly supported for C7 but not for C4. As these are the most conservative results, and similar to the raw RV and CR results, our discussion focuses on them.

Overall vertebral integration and disparity
Results from the eigenvalue dispersion analysis using the homologous-only landmarks for C4-L7 (and therefore not sampling C1, C2, and T10-T13) or the vertebra-specific landmark coordinates were extremely similar for the vertebrae analyzed with both datasets (Table 3). Values for eigenvalue dispersion ranged from 0.226 to 0.307 in the C4-L7 homologous dataset (mean 0.267; median 0.263), and from 0.215 to 0.300 in the vertebra-specific C1-L7 dataset (mean 0.261; median 0.253). Although these values can be considered moderate in the integration spectrum (Pavlicev et al. 2009), in both datasets, vertebrae T10 and L1-L7 presented the highest values of eigenvalue dispersion (>0.27), with the addition of C2 and T11 for the vertebra-specific landmarks analysis.
Procrustes variances across species for the C4-L7 homologous coordinates were the same both before and after accounting for centroid size, and ranged from 0.002 to 0.012, with a mean and a median of 0.005 (Table 3). However, only six vertebrae displayed values of Procrustes variance higher than the mean: C4, T10, L1, L2, L4, and L6, with Procrustes variances of 0.006 for all of these, except T10 with a variance value of 0.012 (Table 3 and Fig. 2). Similarly, the maximum Procrustes distance across specimens per vertebra ranged from 0.109 to 0.296, and only five vertebrae (T10, L1, L2, L4, and L6) presented values higher than the mean and median (0.181 and 0.159, respectively). With both measures of disparity, L7 showed values very close to the mean and higher than the disparity values observed for the anterior vertebrae (with the exception of variance in the atlas).
Regarding the disparity results per vertebra and per individual species, more species presented disparity values that were higher than the mean and median for all vertebrae, both as Procrustes variance and as maximum Procrustes distance, in the general region of T10-L7, and consistently on vertebrae T10, L6, and L7 (Fig. 3, Tables 4 and 5).

DISCUSSION
Here, we have analyzed the patterns of intravertebral modularity and integration throughout the presacral vertebral column and evaluated these patterns in the combination with data on function and morphological disparity. Combined, this work provides a novel view of the evolutionary and developmental forces contributing to vertebral shape differentiation.
The results from our modularity analyses are consistent with the hypothesis that distinct somitic contributions and separate ossification centers in vertebral development result in similar modules in adult vertebral morphology throughout the presacral vertebral column. Only five out of the 19 analyzed vertebrae failed to show support for the two hypothesized modules (centrum and neural spine) based on somitic origins. However, these five vertebrae (C4, T1, T8, L6, and L7) all either form, or are adjacent to, boundaries of traditional vertebral morphological regions, as discussed in detail below.
C4 is part of a previously suggested mammalian developmental module composed of mid-cervicals C3-C5 (Buchholtz et al. 2012). Buchholtz et al. (2012) argued that the commitment of migratory muscle precursor cells from the C3-C5 somites to the formation of the muscularized diaphragm resulted in modular organization of this cervical region, which secondarily contributed to the fixation of cervical number in mammals. Additionally, the cervical region has previously been shown to present its own regionalization, with vertebrae divided into "upper" cervicals (i.e., atlas and axis) responsible for skull articulation, intermediate cervicals (i.e., C3-C5), and "lower" cervicals (i.e., C6 and C7) responsible for neck movement and with morphologies more similar to those seen in the anterior thoracic region (Vidal et al. 1986;Graf et al. 1995;De Iuliis and Puler a 2006;Buchholtz 2014;Arnold et al. 2016). Interestingly, the modular results for felid neck vertebrae expand on the conclusions of a recent study of integration in dog vertebrae (Arnold et al. 2016), in which they found high integration in the cervical (C3-C7) morphology of domestic dogs and suggested that this result can be expanded to a general mammalian pattern.
Our results of eigenvalue dispersion support a hypothesis of moderate integration in the cervical region, although we emphasize that this does not contradict support for developmental modularity within cervical vertebrae. Modularity and integration should not be interpreted as the opposing ends of a spectrum, as modules are typically highly integrated within themselves Klingenberg 2013;Bookstein 2015). As the method used by Arnold et al. (2016) (i.e., PLS, Partial Least Squares; Rohlf and Corti 2000) is mostly suited for testing hypotheses of integration, rather than providing an output value that can discriminate between whole integration and modularity, we suggest that the pattern observed here of developmental modularity for 13 out of 19 vertebrae, including C1, C2, C6, and C7, may also represent a broader mammalian pattern.
Continuing with vertebrae that failed to support the developmental two-module model, the first thoracic vertebra (T1) is at the boundary between the cervical and thoracic regions, with the highly conserved number of seven vertebrae in the mammalian neck (Buchholtz et al. 2012;Buchholtz 2014), and appearance of ribs and consequent reduced mobility in T1. T8 may also be involved in the boundary between two large and more inclusive vertebral regions, although the lack of T9 in our dataset hinders further testing of this hypothesis. Nevertheless, T8 is only two vertebrae away from another previously defined boundary which divides the vertebral column into pre-and postdiaphragmatic regions, T10. This boundary marks the transition between rib-bearing vertebrae, which are restricted by the diaphragm and surround vital organs such as heart and lungs, and the end-thoracics and lumbar region, which undergo more pronounced sagittal bending (Polly et al. 2001  The "eigenvalue dispersion-16" column shows the results for the C4-L7 16 homologous landmarks, while the "eigenvalue dispersion" values are regarding the C1-L7 vertebra-specific landmarks. Bold results mark results higher than the mean and median for the eigenvalue dispersion analyses. analyses also suggest that the posterior region, especially the T10-L7 region, may be more evolutionarily responsive as it shows stronger ecological signal than the anterior column and greater distinction in shape between species showing distinct locomotory specialization (Randau et al. 2016a,b). This vertebral boundary hypothesis can also be adopted toward interpreting the results from L6 and L7, which are the last vertebrae of the presacral region of felids and display higher overall integration. Given the identities and locations of these five vertebrae that do not show a modular structure related to somite origin, we therefore suggest that a functional overprinting of developmental vertebral patterning may occur in these structures in order to maintain larger modular organization of the whole vertebral column (Polly et al. 2001;Buchholtz 2007). However, further analyses in other datasets are needed to confirm this hypothesis.
One unexpected exception to this pattern is T10, which forms one of the most significant boundaries in the vertebral column, but shows significant support for the model of developmental modularity. Based on our hypothesis of vertebral regional boundaries, as well as the results for vertebral disparity and integration, we expected this vertebra to also be an exception to the developmental signal pattern, but it is instead a good example of a structure presenting both a modular organization and an elevated overall integration index. T11 (i.e., the anticlinal vertebra, which marks the change in anteroposterior orientation of the neural spine from caudally to cranially oriented processes) also exhibited high overall integration, and most importantly, the lowest RV and CR values, suggesting the strongest modular organization of the vertebral shape. The consecutive T12 also displayed similarly low and significant RV and CR values. These results suggest that T10-T12 are under strong developmental control (West-Eberhard 1989Arthur 2002;Badyaev et al. 2005), which is maintained even when subjected to varied selection pressures that likely drive the high disparity observed for T10 (and presumably for T11 and T12, although they were not directly compared in disparity due to the lack of homologous landmarks).
The results presented here support the hypothesis that phenotypic integration may promote morphological disparity (Goswami et al. 2014), as observed in the association between higher vertebral overall integration and higher values of morphological disparity (Figs. 2 and 3,. Posterior vertebrae (T10-L7) exhibited the highest degree of overall vertebral integration, as demonstrated by eigenvalue dispersion values higher than both the observed mean and median throughout the vertebral column. These results are particularly interesting when considered with the observation that those vertebrae presented markedly higher values of morphological disparity (both as Procrustes variance and maximum distance) than other vertebrae. We have previously shown that the posterior region is the vertebral section that presented the highest shape differentiation and correlation with ecological specialization in felids, in terms of locomotory mode and prey size, and also allometry (Randau et al. 2016a,b), and this may suggest that this region might display the greatest evolutionary respondability (i.e., raw magnitude of response in any direction to selection, Hansen and Houle 2008) across felids, or even more broadly. Goswami et al. (2014) demonstrated through the use of simulations that integration might increase disparity by coordinating the evolution of traits within functional units and directing this response through paths of higher trait covariance (Klingenberg 2010), although this association has only rarely been supported by empirical data. They have additionally shown that eigenvalue dispersion was highly and significantly positively correlated with respondability. By concentrating variance within determined evolutionary paths the range of morphological diversity is increased, meaning more disparate morphologies may occur than if traits are uncorrelated (Goswami et al. 2014).
Here, we have conducted analyses of vertebral morphological integration and disparity throughout the presacral column of felids, and demonstrated that both measures present their highest values in the posterior axial skeleton, linking these measures to previously demonstrated high levels of ecological diversification (Randau et al. 2016a,b). With this, we add an empirical example of positive association between high integration and disparity to the existing discussion of the role of covariation in promoting versus constraining evolution (Klingenberg 2010;Goswami et al. 2014). Finally, we provided confirmation for the hypothesis that a two-module intravertebral organization driven by somatic origins dominates in the presacral vertebral column in felids, but that this pattern is disrupted, or overprinted, at the boundaries of vertebral regions.  C4  C6  C7  T1  T2  T4  T6  T8  T10  L1  L2  L4  L6 C7  T1  T2  T4  T6  T8  T10  L1  L2  L4  L6