Biological variation is often considered in a scalable hierarchy, e.g., within the individual, within the populations, above the species level. Morphological integration, the concept of covariation among constituent parts of an organism, is also hierarchical; the degree to which these ‘modules’ covary is a matter of the scale of the study as well as underlying processes driving the covariation. Multilevel analyses of trait covariation are a valuable tool to infer the origins and historical persistence of morphological diversity. Here we investigate concordance in patterns of integration and modularity across three biological levels of variation: within a species, within two genera-level radiations, and among species at the family level. We demonstrate this approach using the skull of mammalian family Leporidae (rabbits and hares), which is morphologically diverse and has a rare-among-mammals functional signal of locomotion adaptation. We tested three alternative hypotheses of modularity; from the most supported we investigated disparity and integration of each module to infer which is most responsible for patterns of cranial variation across these levels, and whether variation is partitioned consistently across levels. We found a common pattern of modularity underlies leporid cranial diversity, though there is inconsistency across levels in each module’s disparity and integration. The face module contributes the most to disparity at all levels, which we propose is facilitating evolutionary diversity in this clade. Therefore, the distinctive facial tilt of leporids is an adaptation to locomotory behavior facilitated by a modular system that allows lineages to respond differently to selection pressures.

All analyses were performed using the R Statistical Environment v. 4.2.1 (R Development Core Team 2022), using the geomorph package v.4.0.4 (Adams et al. 2022) unless otherwise stated, with statistical significance estimated using a permutation approach (1000 iterations) and evaluated at the significance level of 5%.

Samples

We sampled 317 specimens from 22 species across Leporidae representing the main lineages (Table S1), many of which featured in a previous study (Kraatz and Sherratt 2016), using predominantly museum collections (details in Supplementary file 1). For Oryctolagus cuniculus and Lepus europaeus, we could increase sampling with whole carcasses scavenged from pest control activities in Australia (ethics approved, details in acknowledgements), where these species were introduced on multiple occasions during the 1800s by European settlers from UK stocks (Peacock and Abbott 2013; Stott 2015). As the most widely distributed of all leporids, we restricted the O. cuniculus (herein Oryctolagus for brevity) samples in this study to those from Europe and Australia to limit population-level variation.

Specimens were scanned by X-ray computed tomography (CT) to obtain digital models of the cranium for measurement. Scans of museum specimens were made using a J. Morita Veraviewepocs 3D R100 system (College of Dental Medicine, Western University of Health Sciences), typically operated at 75 kV and 3 mAs with voxel size ranging from 125–160 µm. Whole specimens of Oryctolagus and L. europaeus were scanned using a Siemens SOMATOM Force CT scanner (Dr Jones & Partners, The South Australian Health and Medical Research Institute [SAHMRI], Adelaide). The scanner was operated at 120 kVp and 200 mA, with one second exposure and a slice thickness of 0.4 mm.

Scan acquisition and landmark placement

The CT tomograms were processed with Checkpoint (Stratovan Corporation, Davis, CA), thresholding by voxel grey value to obtain an isosurface representing bone. Then landmarks were manually digitized on the crania models by one author (BK); following Kraatz and Sherratt (2016), 44 landmarks were placed at homologous points on the cranium, over the left and right sides, and the curve tool was used to place eight equally spaced semilandmarks along the sagittal axis of the cranial roof to capture the arching curvature (Figure 1, Table S2). Zygomatic arches were then digitized with the landmark tool using as many points necessary to capture the complex curves, and these landmarks were converted into 11 equally spaced semilandmarks each using the ‘digit.curves’ R function. Missing landmarks, typically reflecting breaks in supraorbital process or regions of the basicranium, were an issue for 45 specimens (details in Figshare repository). The positions of these missing landmarks were estimated using a thin-plate spline approach (Gunz et al. 2009), where each missing landmark is predicted among all other homologous landmarks of complete specimens of the same species, implemented using the ‘estimate.missing’ function.

The landmark coordinate data were aligned using a generalized Procrustes superimposition (Rohlf and Slice 1990), taking into account object symmetry and allowing semilandmarks of the cranial roof and zygomatic arches to slide along their tangent directions in order to minimize bending energy (Gunz et al. 2005), implemented with ‘bilat.symmetry’ function. The resulting symmetric component of shape (sensu Klingenberg et al. 2002) was used in the following analyses.

Hypotheses of modularity

Preliminary analysis rejected the hypothesis of no modular structure in the leporid cranium, so we proceeded with evaluating three alternative modular hypotheses proposed for mammalian skulls (Cheverud 1995; Willmore et al. 2006; Hallgrímsson et al. 2007) (Table 1). The ‘tissue origin’ hypothesis (Figure 1b) use two modules based upon whether bones originate from neural crest or mesoderm cells (Willmore et al. 2006). The ‘developmental groups’ hypothesis (Figure 1c) identifies three modules based upon embryological development (Hallgrímsson et al. 2007), where the basicranium is derived from the chondrocranium, the neurocranium from dermatocranium bones of the cranial vault, and the face from the splanchnocranium initially with subsequent influence of dermatocranial elements. The ‘functional groups’ hypothesis (Figure 1d) has four modules based upon different functional regions of the cranium (Cheverud 1995; Willmore et al. 2006) and related to the six-module hypothesis found across many mammals (Cheverud 1982a; Goswami 2006). The six-module hypothesis was not investigated here due its similarity to the ‘functional groups’ hypothesis and because the derived morphology of the elongated rostrum of leporids, which is divided dorso-ventrally in the six-module hypothesis, has been shown in previous studies (Kraatz and Sherratt, 2016) to be a key region that varies coherently within the leporid cranium.

Datasets

The Procrustes landmark data were divided into three datasets that represent different biological levels. The first level is within-species variation, and the monospecific Oryctolagus was used because it has the largest sample size (N=60) and is an exemplary biological model, and often used to represent leporids as a whole. The second level is the within-genus level, where the hares (or jack-rabbits, Lepus) and cottontail rabbits (Sylvilagus), were used because they are the most speciose in Leporidae and both represent monophyletic clades. Here they are represented by six and seven species each (N=97 and N=110 respectively). We followed Cano-Sánchez et al. (2022) in their taxonomic recommendation that the pygmy rabbit Sylvilagus idahoensis (Merriam, 1891; previously Brachylagus) be included with the Sylvilagus radiation (e.g., Grinnell et al. 1930) because phylogenetic inference has consistently recovered this taxon as sister or nested within Sylvilagus (Matthee et al. 2004; Cano-Sánchez et al. 2022). To account for the group structure in each of these levels, an analysis of variance (ANOVA), implemented with ‘procD.lm’ function, was used to evaluate the proportion of variance attributed to population for Oryctolagus, and species for Sylvilagus and Lepus. Residuals of these ANOVAs were extracted to represent the pooled-within group shape for subsequent analyses of integration and disparity.

The third level is the among-species “evolutionary” level, represented by species-averaged cranial shape data of 22 species across the family (e.g., Klingenberg and Marugan-Lobon 2013). Other than Lepus and Sylvilagus, the family is represented by monotypic or low diversity genera. We used a published phylogenetic hypothesis based upon seven genes (five nuclear and 2 mitochondrial) (Matthee et al. 2004). The tree was pruned to 21 species included in this study using ‘drop.tip’ function in ape v.5.6-2 (Paradis et al. 2004). Lepus europaeus was not included in Matthee et al.’s original analysis. Based upon other molecular analyses and the propensity for L. europaeus to hybridise (Alves et al. 2008b; Ben Slimen et al. 2008; Suchentrunk et al. 2008; Melo-Ferreira et al. 2012; Ashrafzadeh et al. 2018), we applied the analyses to two alternate topologies: L. europaeus sister to L. capensis, and sister to L. timidus. We grafted L. europaeus onto the tree using ‘bind.tip’ function in phytools v.1.0-1 (Revell 2012); preliminary analysis tested different branch lengths with inconsequential changes to results, so we positioned the node halfway along the original branch in both topologies (Figure S2).

Statistical Analyses

Principal components analysis was used to identify over how many axes the shape variance was distributed at each level, and thus gives an indication of the magnitude of integration for the whole cranium, where higher integration results in more variance contained in relatively few dimensions. This was implemented using the ‘gm.prcomp’ function. Because integration is expected to be influenced by allometry, the covariation of shape and size, we investigated two types of allometry: static allometry, which is the relationship of shape and size among individuals of the same age class in a species/clade; and evolutionary allometry, which is the relationship of shape and size among species, and defined as the covariation of shape and size along branches of the phylogenetic tree. Multivariate regressions were used to assess the amount of shape variation of the whole cranium attributed to allometry (Monteiro 1999) at each level. The regression score approach (Drake and Klingenberg 2008) was used to visualise the allometric relationships with log-transformed centroid size, a measure of size derived from the landmark coordinates and calculated during Procrustes superimposition. The degree of static allometry was assessed in Oryctolagus, and Sylvilagus and Lepus genera by performing a regression on all specimens using ‘procD.lm’. The regression model included species (or population for Oryctolagus), allowing the interaction of size and species to be evaluated. Evolutionary allometry was examined among species using a phylogenetic ANOVA (PGLS) (Adams), with species mean values of shape and size, implemented with ‘procD.pgls’. The residuals of the statistically significant models were extracted and used to represent shape with the allometric component removed (“allometry-free”).

We tested the three a priori defined hypotheses of modularity and evaluated which was the most supported for the data, using the covariance ratio (CR) and post-hoc effect size approach (Adams and Collyer 2016; Adams and Collyer 2019). The CR coefficient was calculated for each hypothesis using the ‘modularity.test’ function. Effect sizes from each modularity analysis were evaluated using the ‘compare.CR’ function to infer which of the modular hypotheses is most supported by the data. This test was performed using the pooled-within species residuals and allometry-free datasets of Oryctolagus, Lepus and Sylvilagus genera. For the evolutionary level, we used the ‘phylo.modularity’ function to calculate CR in a phylogenetic context, which used the evolutionary covariance matrix among traits found under a Brownian motion model of evolution (Adams and Felice 2014).

Analyses of modularity and integration typically require large sample sizes (e.g., minimum of 30 specimens per species (Haber 2011)), which can be problematic for macroevolutionary studies. To understand the uncertainty that sample size may give in our results, we used a subsampling approach taking 10 specimens at random from the Oryctolagus dataset, calculated the test for modularity (details above), and repeated this 100 times to obtain an estimate of the margin of error. We then calculated the CR of each hypothesis for all species with a sample size of at least 10 to provide preliminary insights into modular variation across the clade.

The consensus, best supported modular hypothesis at each level was used to investigate which module contributes the most to the observed variation among individuals, and how cranial modules vary in terms of morphological disparity and strength of within and between module integration. Each module was subjected to a separate GPA to avoid correlations between modules due to the superimposition (Cardini 2019b), transformed to account for pooled-within group variation as above, and multiplied by module centroid size to ensure each module represents the original relative proportions (necessary for disparity). Since GPA influences the covariance structure among landmarks, whether to use a single or separate GPA is a matter of discussion in geometric morphometrics (Cardini 2019b; Goswami et al. 2019; Klingenberg 2021); our preliminary analysis showed a separate GPA changes the magnitude of results, but the overall among-module patterns are equivalent to a single GPA (Figure S1). We used Mantel’s test to compare a distance matrix of between-individual distances calculated from all landmarks to one calculated from landmarks within a module. This was implemented with the ‘mantel’ function in vegan R package v.2.6-2 (Oksanen et al. 2022). The highest correlation value indicates that module contributes the most to the overall morphospace. To measure morphological disparity, we used the Procrustes variance (the sum of the diagonal elements of the covariance matrix (Zelditch et al. 2012)) divided by the number of landmarks in the module. This was implemented using the ‘morphol.disparity’ function. Higher values represent greater variation in shape among observations (wider spread in morphospace). To measure the magnitude of integration within each module, we calculated the relative eigenvalue variance (Pavlicev et al. 2009) and from this the standardised effect score (Z-score) as per (Conaway and Adams 2022) implemented with ‘integration.Vrel’ function. This provides a directly comparable measure of the magnitude of integration across levels and modules, based upon the degree of eigenvalue dispersion, that is how variance is distributed across eigenvalues. Higher negative values of the effect score represent lower covariance between shape traits, which means weaker integration. The effect sizes were compared for significant differences using the ‘compare.ZVrel’ function. Finally to measure the magnitude of integration between modules at each level, two-block partial least squares analysis (PLS) was used (Rohlf and Corti 2000; Adams and Felice 2014), implemented with functions ‘two.b.pls’ and ‘phylo.integration’. This approach calculates the correlation coefficient (r-PLS) between two matrices (‘blocks’) of traits, where values closer to 1 mean stronger integration between modules.

To understand how much phylogenetic relatedness is driving morphological variation, phylogenetic signal was evaluated for the whole skull and per module using the K_mult approach (Adams 2014a). The metric, based upon K (Blomberg et al. 2003), measures the proportion of the tip variance that is explained by a Brownian motion (BM) process for the tree, thus testing whether relatives resemble each other less than expected under BM (K less than 1), which suggests a departure from Brownian motion evolution and indicates homoplasy. Where there is a K value greater than one, close relatives are more similar than expected under BM (strong phylogenetic structure). This was implemented with ‘physignal’ function.

R Statistical Environment, R 4.2.1

library(geomorph) # v.4.0.4

library(ape) # v.5.6-2

library(vegan) # v.2.6-2

library(phytools) # v.1.2-0