DOI: 10.5061/dryad.4mw6m90f3
title: Data for: Multilevel analysis of integration and disparity in the mammalian skull

description: 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.
description: 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%.
 <strong>Samples</strong> 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 <em>Oryctolagus
 cuniculus</em> and <em>Lepus europaeus</em>, 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
 <em>O. cuniculus</em> (herein <em>Oryctolagus</em>
 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
 <em>Oryctolagus</em> and <em>L. europaeus</em>
 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.
 <strong>Scan acquisition and landmark
 placement</strong> 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.
 <strong>Hypotheses of modularity</strong>
 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.
 <strong>Datasets</strong> The
 Procrustes landmark data were divided into three datasets that represent
 different biological levels. The first level is within-species variation,
 and the monospecific <em>Oryctolagus</em> 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,
 <em>Lepus</em>) and cottontail rabbits
 (<em>Sylvilagus</em>), 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 <em>Sylvilagus
 idahoensis</em> (Merriam, 1891; previously
 <em>Brachylagus</em>) be included with the
 <em>Sylvilagus</em> 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
 <em>Oryctolagus</em>, and species for
 <em>Sylvilagus</em> and <em>Lepus</em>. 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
 <em>Lepus</em> and <em>Sylvilagus</em>, 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
 <em>ape</em> v.5.6-2 (Paradis et al. 2004). <em>Lepus
 europaeus</em> was not included in Matthee et al.’s original
 analysis. Based upon other molecular analyses and the propensity
 for<em> L. europaeus</em> 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: <em>L. europaeus</em> sister to <em>L.
 capensis</em>, and sister to <em>L. timidus</em>. We
 grafted <em>L. europaeus</em> onto the tree using ‘bind.tip’
 function in <em>phytools</em> 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).
 <strong>Statistical Analyses</strong>
 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 <em>Oryctolagus</em>, and
 <em>Sylvilagus</em> and <em>Lepus</em> genera by
 performing a regression on all specimens using ‘procD.lm’. The regression
 model included species (or population for
 <em>Oryctolagus</em>), 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<em> a priori</em> 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
 <em>Oryctolagus</em>, <em>Lepus</em> and
 <em>Sylvilagus</em> 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
 <em>Oryctolagus</em> 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 <em>vegan</em> 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 <em>within</em> 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 <em>between</em> 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<sub>mult</sub> 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.
description: 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