Evolutionary integration of the frog cranium

Evolutionary integration (covariation) of traits has long fascinated biologists because of its potential to elucidate factors that have shaped morphological evolution. Studies of tetrapod crania have identified patterns of evolutionary integration that reflect functional or developmental interactions among traits, but no studies to date have sampled widely across the species‐rich lissamphibian order Anura (frogs). Frogs exhibit a vast range of cranial morphologies, life history strategies, and ecologies. Here, using high‐density morphometrics we capture cranial morphology for 172 anuran species, sampling every extant family. We quantify the pattern of evolutionary modularity in the frog skull and compare patterns in taxa with different life history modes. Evolutionary changes across the anuran cranium are highly modular, with a well‐integrated “suspensorium” involved in feeding. This pattern is strikingly similar to that identified for caecilian and salamander crania, suggesting replication of patterns of evolutionary integration across Lissamphibia. Surprisingly, possession of a feeding larval stage has no notable influence on cranial integration across frogs. However, late‐ossifying bones exhibit higher integration than early‐ossifying bones. Finally, anuran cranial modules show diverse morphological disparities, supporting the hypothesis that modular variation allows mosaic evolution of the cranium, but we find no consistent relationship between degree of within‐module integration and disparity.

Trait integration is an inherent property of biological systems and manifests at many scales of analysis, from populations to large clades. Correlations among traits, and their organization into highly integrated, semi-autonomous modules, can result from genetic, developmental, or functional interactions, and it is thought that patterns of genetic and developmental integration may evolve adaptively to promote functional integration (Riedl 1978;Cheverud 1984;Wagner and Altenberg 1996). These trait associations may be replicated in macroevolutionary patterns observed across larger clades, with traits evolving in either a coordinated or a modular fashion. The modular organization of phenotypic evolution is central to many fundamental concepts in evolutionary biology, including mosaicism (De Beer 1954), which have been suggested to promote the diversification of form. Several recent studies have explored this hypothesized relationship by quantifying phenotypic (or variational) and evolutionary integration and modularity and determining how they relate to disparity and rates of evolution in various species or clades. These studies provide extensive evidence that the presence of evolutionary modules is associated with higher disparity or evolutionary rate (Goswami and Polly 2010;Claverie and Patek 2013;Randau and Goswami 2017;Larouche et al. 2018;Dellinger et al. 2019;Watanabe et al. 2019; Bardua et al. 2019b). Many studies have also assessed whether the magnitude of evolutionary integration among traits (across an entire phenotype or within a module) correlates with higher or lower disparity or evolutionary rate. For this question, the answer is less clear than for the one above. Some studies find that high evolutionary integration among traits is correlated with low disparity and slow evolutionary rate (Claverie and Patek 2013;Martín-Serra et al. 2019). Conversely, others have found that highly integrated traits show greater disparity than less integrated ones. In addition, some have found that there is no relationship between evolutionary modularity and disparity or rate of evolution at this scale (Goswami et al. 2014;Watanabe et al. 2019; Bardua et al. 2019b). However, our understanding of the factors shaping the evolutionary integration among traits and its relevance for morphological evolution is incomplete, with most studies focused on a few clades, such as mammals and birds, both of which have relatively little developmental diversity. Here, we have expanded the study of evolutionary modularity to one of the most taxonomically, developmentally, ecologically, and morphologically diverse clades of tetrapods: the lissamphibian order Anura (frogs).
Lissamphibians constitute over 25% of extant tetrapod species (totalling 8,160 species;AmphibiaWeb 2020), and yet patterns of trait integration at any scale across this diverse clade remain largely unexplored compared to amniotes. Previous studies of both phenotypic and evolutionary integration in lissamphibian crania have found limited support for modular structure, with studies recovering either few modules or no modular structure at all (Ivanović and Kalezić 2010;Sherratt 2011;Ivanović and Arntzen 2014;Simon and Marroig 2017;Vidal-García and Keogh 2017). With the application of high-dimensional geometric morphometric data, strong modular structure has been identified at the intraspecific (phenotypic) and evolutionary levels across caecilian (Marshall et al. 2019; Bardua et al. 2019b) and salamander Fabre et al. 2020) crania, in analyses exploring much wider ranges of models. However, the most diverse lissamphibians, frogs, have not yet been incorporated into these studies. The few studies of cranial modularity in frogs to date vary widely in taxonomic and morphological coverage as well as data type, hindering our understanding of anuran cranial modularity and preventing direct comparison of modular structures across Lissamphibia (e.g., Simon and Marroig 2017;Vidal-García and Keogh 2017). Beyond their taxonomic diversity, anurans provide a unique opportunity to investigate patterns of cranial integration across lineages varying markedly in life history, including the repeated loss of a free-living and feeding larval stage. The presence of such a larval stage has been hypothesized to drive the decoupling of genetic and developmental traits with functional traits, as the distinct, divergent selection pressures associated with larval and adult ecological niches may drive low genetic correlations between these two life history stages ("adaptive decoupling hypothesis"; Ebenman 1992; Moran 1994). Indeed, tadpoles and frogs can evolve independently, responding to antagonistic selection pressures , and different tadpole morphologies can converge to the same adult morphology (Bragg and Bragg 1958;Pfennig 1990). Frogs with a freeliving, feeding larval stage may therefore experience fewer developmental or genetic constraints, which may lead to a decrease in strength of phenotypic integration and a greater partitioning of traits into modules, or to greater variation in patterns of phenotypic integration across taxa, both of which would be reflected in greater evolutionary modularity. Some support for this hypothesis is offered by recent analysis of evolutionary modularity in salamanders, with species undergoing complete metamorphosis exhibiting elevated cranial modularity compared with pedomorphic taxa . Life history may therefore have a significant and persistent influence on the strength and pattern of cranial integration, and this effect is expected to be particularly strong in lineages with disparate free-living developmental stages, such as frogs.
Other aspects of development may also influence evolutionary integration across frogs. For example, derivation from different cranial neural crest (CNC) streams may generate developmental modules that may be expected to evolve in a coordinated manner . The pattern of the contribution of CNC streams to the osteocranium is considerably different for the frog Xenopus than for the axolotl (a salamander) or for amniotes (Hanken and Gross 2005;Gross and Hanken 2008;Piekarski et al. 2014), suggesting a possible deep divergence in the pattern of cranial development between frogs and other tetrapod clades. Furthermore, frogs differ from other lissamphibian clades in cranial ossification sequence timing, with a heterochronic shift whereby frogs exhibit relatively later ossification of bones associated with adult feeding (Harrington et al. 2013). Larval and adult feeding modes differ more for frogs than for salamanders or caecilians, so larval frog mouthparts would likely be more impeded by early ossification of bones involved in adult feeding. However, direct-developing frogs, which feed like adults immediately following hatching, partially reverse this trend (Hanken et al. 1992;Kerney et al. 2007;Harrington et al. 2013). Cranial ossification sequences therefore vary across Anura (Weisbecker and Mitgutsch 2010), and can even vary intraspecifically (Reiss 2002;Moore and Townsend Jr 2003), which may result in greater developmental modularity between bones (Weisbecker and Mitgutsch 2010). Moreover, because early-ossifying bones are generally conservative in timing and function (protecting the brain and otic capsules, Duellman and Trueb 1986;Heatwole and Davies 2003), these elements may be expected to display higher evolutionary integration than later ossifying elements. Ossification sequence timing and derivation from CNC streams may therefore both influence the pattern and strength of evolutionary integration in frog crania and contribute to differences in evolutionary integration between frogs and other tetrapod clades.
Here, we quantify cranial morphology across 172 anuran species, sampling every family and an extensive range of ecologies and developmental strategies. We implement a high-density data approach (detailed in Bardua et al. 2019a), which has proved successful in capturing morphology across a range of structures and clades (Dumont et al. 2015;Parr et al. 2016;Watanabe et al. 2019;Bardua et al. 2019b). This work represents the most comprehensive analysis of amphibian cranial integration to date, in terms of both taxonomic sampling and density of shape data. With these shape data, we investigate patterns of evolutionary integration and modularity in the frog skull, testing a range of developmental and functional models. We further investigate developmental influences on cranial integration, first by comparing pattern and magnitude of integration in taxa with and without a feeding larval stage, and second by assessing the relationship between magnitude of integration and relative timing of ossification for each cranial bone. Finally, we address the macroevolutionary significance of anuran cranial integration by quantifying the relationship between magnitude of integration and morphological diversity (disparity).

SPECIMENS
Our dataset consists of the crania of 172 extant anuran species (one specimen per species), encompassing representatives from all extant families of frogs (Table S1). All specimens were spirit preserved. Determining juvenile from adult specimens can be difficult as adult crania are variably ossified across Anura (Trueb 1973;Nishikawa 2000), so the largest specimen available in the collections visited for each species was selected. Due to lack of data and/or the availability of specimens we did not control for sex in data collection, but it is unlikely that sex-related shape dimorphism would significantly affect results at this scale of macroevolutionary analysis (see Sherratt et al. 2014;Bardua et al. 2019b). Specimens were micro-CT scanned and threedimensional isosurface models ("meshes") were created from segmenting the volumes using Avizo version 9.3 (FEI, Hillsboro, OR, USA) and VG Studio MAX version 2.0 (Volume Graphics 2001). Meshes were prepared and cleaned in Geomagic Wrap (3D Systems) by removing the mandible from each specimen and removing small foramina texturing the surface, because these can hinder the collection of surface semilandmarks (Bardua et al. 2019a). Data were taken from the right side of each cranium, so 11 specimens were mirrored when this side was incomplete or damaged (Table S1).

PHYLOGENY
Anuran phylogenetic relationships have been the subject of extensive recent study. To incorporate evolutionary relationships into analyses, we used the most recent, most comprehensive phylogeny of Anura (Jetz and Pyron 2018) (Fig. 1). One hundred and sixty-eight anuran species that we sampled are present in this phylogeny, and 140 of our specimens have a phylogenetic position based on direct analysis of molecular data. The consensus phylogeny from Jetz and Pyron (2018) was pruned using the "drop.tip" function in the R package ape (Paradis et al. 2004). The four remaining specimens did not have species assignment, so these were added to the phylogeny at the appropriate position for their genus, to allow inclusion in the phylogenetically informed analyses (Raorchestes sp., Dendrobates sp., Capensibufo sp., and Xenorhina sp. were added as Raorchestes anili, Dendrobates auratus, Capensibufo rosei, and Xenorhina varia, respectively).

LIFE HISTORY
There is a high degree of variation in life history across Anura (Duellman and Trueb 1986;Wells 2010), with nearly 40 reproductive modes defined based on the site of egg development (Haddad and Prado 2005). We coded specimens by life history in terms of possession of feeding larval stage (Table S1). Directdeveloping taxa, and some metamorphic species, have nonfeeding larvae provisioned solely with yolk ("endotrophic"), whereas most metamorphic taxa have feeding larvae that feed on external food sources ("exotrophic") (McDiarmid and Altig 1999). This categorization distinguishes larvae that interact with their environment (and presumably experience strong environmental pressures as larvae) from those that do not or do so in a limited way.

MORPHOMETRIC DATA COLLECTION
Among living amphibians, anurans have posed challenges for collecting high-dimensional morphological data. With currently 7,204 species (AmphibiaWeb 2020), frogs are nearly 10 times more species rich than salamanders (742 species), and over 30 times more speciose than caecilians (214 species). Anurans are also the most diverse in larval and adult cranial structure of the extant amphibian orders (Trueb 1993). Many cranial bones are lost as discernible elements repeatedly across this clade, including through fusion between elements or failure to ossify, as well as novel bones evolving in specific subclades (e.g., see Trueb, 1970Trueb, , 1973Hall and Larsen 1998;Campos et al. 2010;Schoch 2014;Pereyra et al. 2016;Lambert et al. 2017). Traditional morphometric approaches would fail to capture the morphology of these variably present bones. A high-dimensional approach allows inclusion of this extreme variation in morphology and in bone presence, as this variation may have a crucial influence on patterns of trait integration and modularity across Anura.

FROG CRANIAL INTEGRATION
iii. iv. v.

Regions
We divided the right side of each frog cranium into 19 regions, which represented the maximum reasonable partitioning of the cranial morphology ( Fig. 1). In most cases, regions were simply whole bone surfaces, but some bones could be further divided into two distinct regions each (premaxilla, maxilla, otooccipital, quadratojugal, and sphenethmoid; see Table S2 for details). All regions were defined by clear anatomical structures and mostly represented homologous regions, although functionally analogous regions had to be combined in rare cases, for example, when novel bones were present (Table S2). In addition, the otooccipital (exoccipital and opisthotic) can in some cases indistinguishably fuse with the prootic, so for these specimens the "occipital" (∼ otooccipital) and "otic" (∼ prootic) regions were divided along the posterior epiotic ridge. Approximately half of our specimens had at least one absent region, and over one third of the regions were variably present. Removing regions that were variably present (or specimens lacking these regions) would greatly impact either our quantification of morphology or our sample size (Table S1). To incorporate the full range of skull variation in our analyses, we represented absent regions as "negligible regions," that is, a position of zero size (or near-zero following Procrustes Alignment) described below and previously (Bardua et al. 2019b).

Landmarks and curve semilandmarks
Each cranial region was defined by Type I and Type II landmarks (Bookstein 1991) linked by "curves" (consisting of sliding semilandmarks; Gunz et al. 2005) (Fig. S1). Landmarks and curves were digitized in IDAV Landmark Editor version 3.6 (Wiley et al. 2005). A total of 58 landmarks and 59 curves (ranging from two to 12 semilandmarks, Tables S3-S4) were placed on the right side of each cranium to define all regions present across all specimens. These curves were resampled (for a description and code, see Supporting Information in Botton-Divet et al. 2016), resulting in a total of 410 curve semilandmarks. An additional 24 landmarks and 24 curves were used to define variably present regions in specimens that had those regions represented.

Surface semilandmarks
A semi-automated procedure in the R package Morpho version 2.5.1 (Schlager 2017) was used to project surface semilandmarks from a template onto each present cranial region, as described (Schlager 2017) and previously implemented Marshall et al. 2019; Bardua et al. 2019b). Variably present regions were patched as normal when present, with the additional landmarks and curves defining these regions. These landmarks and curves were then removed prior to analyses, so that these regions were represented only by surface semilandmarks. When variably present regions were absent, these regions were represented by one coordinate position (that best represented the location of each missing region), replicated to achieve an array of dimensions matching that of present regions. This method has been previously implemented (Bardua et al. 2019b), and a similar method was suggested for incorporating novel structures (see fig 1B from Klingenberg 2008). Variably present regions were thus represented only by surface points. A total of 527 surface semilandmarks were applied evenly across each cranium (see Table S5), so that each cranium was represented by a total of 995 landmarks and semilandmarks.

Procrustes alignment
Non-shape aspects of our data (translation, rotation, and scale) were removed through Procrustes Alignment, using the "gpagen" function in the R package geomorph version 3.1.3 (Adams and Otárola-Castillo 2013; Adams et al. 2017). Our data were mirrored prior to alignment, as a bilaterally symmetrical structure aligns more successfully (Cardini 2016). Anuran crania only have two midline landmarks (anterior and posterior extremes of the parasphenoid), so additional midline positions were created by finding the midpoint of two bilaterally symmetrical landmarks (anteromedial and posteromedial extremes of the frontoparietal). Data were mirrored using the "mirrorfill" function in paleomorph version 0.1.4 (Lucas and Goswami 2017). Following alignment, the mirrored data were removed, as these data were redundant.

Phylogenetic signal and correction
We calculated K mult , the phylogenetic signal in our shape data under the assumption of Brownian motion, using the "physignal" function in the R package geomorph (Adams 2014). A phylogenetic correction was applied to our shape data to account for shared evolutionary history and presumed increased similarity between more closely related species. We computed phylogenetic independent contrasts (Felsenstein 1985) for our shape data and used these in further analyses. We also generated a morphospace for the first two phylogenetic principal components to visualize the spread of cranial shapes within our dataset (Fig. S2).

Allometric signal and correction
Size-related shape changes were investigated by quantifying evolutionary allometry in our shape data. We used the "procD.pgls" function in geomorph, which conducts a phylogenetic generalized least squares analysis using log centroid size. Centroid size is the square root of the sum of squared distances of all landmarks from the center (centroid) of a structure. For the centroid size of each cranium, see Table S6. We also visualized size-related shape changes using the "procD.allometry" function in geomorph v3.0.5.  To account for size-related shape changes, we corrected our shape data for allometry. We performed a Procrustes ANOVA using the "procD.lm" function in the R package geomorph, with log centroid size as a factor, and used the residuals from this analysis in modularity analyses.

Modularity
We hypothesized 27 different model structures, ranging from a fully integrated cranium (one "module"), to every bone or every region as its own module (14 or 19 modules) (see Table 1). We compared a range of models based on function, development, and ossification sequence rank, including models modified from previous studies (Simon and Marroig 2017;Vidal-García and Keogh 2017). These include models based on the contribution of CNC streams to cranial bones, and various divisions of the cranium based on hypothesized functional units. The ossification sequence rank model uses the median rank positions of cranial bones from Weisbecker and Mitgutsch (2010). We also included a model analogous to the 10-module model recovered across cae-cilian crania (Bardua et al. 2019b). For details on the models, see Text S1.
We investigated patterns of trait integration and modularity using two methods: Evaluating Modularity with Maximum Likelihood (EMMLi) (Goswami and Finarelli 2016) and Covariance Ratio (CR) analysis (Adams 2016), both of which have been implemented in previous analyses of modularity using highdimensional shape data Marshall et al. 2019;Watanabe et al. 2019;Bardua et al. 2019b). EMMLi is a maximum likelihood approach that allows the testing of multiple hypotheses of modularity, each of which can vary in its number of modules. This is implemented using the "EMMLi" function in the EMMLi R package (Goswami and Finarelli 2016). We tested our 27 different model structures with the shape data, as well as with the phylogenetically-and allometry-corrected shape data. To assess the robustness of our results, we subsampled our shape data down to 10% (using random jackknife resampling), and ran EMMLi iteratively 100 times with these subsampled data using the "subSampleEMMLi" function in EMMLiv2 (https://github.com/hferg/EMMLiv2/). We compared the average results from these 100 runs to the results from the original analysis. We also ran EMMLi with the landmark-only dataset to compare results using different data types. The landmark-only dataset excludes variably present regions, as these are represented only by surface points. To test whether the creation of "negligible regions" imposed artificially high integration on variably present regions (because surface points occupy˜identical positions in negligible regions), we ran phylogenetically-corrected EMMLi analysis for just the specimens that have every region present (N = 83). We compared the pattern of trait integration from this analysis to the original analysis, as well as determining whether within-region trait correlations were different.
We observed the pattern of trait correlations between each module for the best-supported model. However, EMMLi is not exhaustive in its comparison of models, and recent analyses suggest that EMMLi tends to favor more highly parameterized models (Adams and Collyer 2019;Felice et al. 2019;Marshall et al. 2019;Watanabe et al. 2019;Bardua et al. 2019b;Fabre et al. 2020). As such, we used the estimated within-and between-module correlations (ρ) for the best-supported model and grouped highly integrated regions into larger modules as previously described Marshall et al. 2019; Bardua et al. 2019b), to construct an alternative model of modularity. We grouped regions when the between-module correlation for two modules was within 0.2 of the lowest withinmodule correlation (Bardua et al. 2019b).
Model selection approaches such as EMMLi only select the best fit model, but they do not explicitly test hypotheses of modularity. For this reason, we also used Covariance Ratio analysis. Covariance Ratio analysis compares the overall covariation between hypothesized modules to the covariation within those modules (Adams 2016). We conducted CR analysis (Adams 2016) for the model identified from EMMLi analysis, using the "modularity.test" function in geomorph, to observe whether both methods support similar patterns of modularity. Specifically, we ran this analysis six times: for the uncorrected, allometry-corrected, and phylogenetically-corrected complete datasets, as well as for the uncorrected, allometry-corrected, and phylogenetically-corrected landmark-only datasets.
We further ran EMMLi and CR analyses on subsets of specimens, divided based on presence or absence of a feeding larval stage. We split the dataset into species possessing "feeding" and "nonfeeding" larvae prior to Procrustes alignment, so that both datasets were aligned separately. We then performed phylogenetically informed EMMLi and CR analysis for each subset. However, our imbalance in the number of species with (N = 124) and without (N = 39) a feeding larval stage may affect strength of integration. We therefore also ran EMMLi and CR analyses iteratively 100 times each on subsets of the larger "feeding larvae" group (N = 124), taking 39 specimens at random each time.
We then calculated the mean estimated correlations/covariations from these 100 runs and compared these results to the results from the original complete "feeding larvae" dataset (N = 124), and to the results from the smaller "nonfeeding larvae" group (N = 39). We then compared the effect sizes for strength of modularity (Z-scores) of the three datasets using the "compare.cr" function in geomorph.

Disparity
Disparity (morphological diversity) of each cranial module was defined as Procrustes variance (calculated using the "morphol.disparity" function in geomorph), divided by module landmark/semilandmark number, to correct for landmark/ semilandmark number (as this affects variance). A regression of disparity on the estimated magnitude of integration (withinmodule ρ) from the EMMLi analysis was conducted to understand the influence of the latter on the former.

Ossification sequence rank
The median rank position within the frog cranial ossification sequence was taken for each cranial bone from Weisbecker and Mitgutsch (2010). The stapes was excluded as ossification sequence information was absent for this bone, and the "suspensorium" module was split into its individual bones (quadratojugal, squamosal, and pterygoid). The Spearman's rank correlation coefficient was then calculated between the magnitude of integration (within-module ρ) and the median ossification position for each cranial bone to determine the relationship between these two metrics.

MODULARITY
EMMLi analyses for the uncorrected, phylogeneticallycorrected, and allometry-corrected data all recovered the most parametrized model (the 19-region model) as the best supported. In all three analyses, very similar patterns of trait integration were identified between cranial regions (Figs. 2 and S3; Tables S7-S9). For the reasons noted above, comparisons of withinand between-region trait correlations were conducted, resulting in a novel 13-module model recovered from all three analyses. This new model resulted in the construction of the following multi-region modules: (1) a suspensorium module (quadratojugal (jaw joint and lateral process), pterygoid, and squamosal); (2) a maxilla module (both dorsal and ventral regions); (3) a   Table S2.
Resampling our data down to 10% and taking the average of 100 runs returned very similar results to the full run. The pattern of trait integration was the same, and the identical 13-module model was recovered following the steps outlined above ( Fig. S3; Table S10).
The landmark-only EMMLi analyses (uncorrected, phylogenetically-corrected, and allometry-corrected) showed very similar results to one another with only subtle differences ( Fig. S4; Tables S11-S13). Two highly integrated modules were recovered for all three landmark-only EMMLi analyses: (1) an occipital module (the occipital and occipital condyle) and (2) a "facial module," consisting of both regions of both the maxilla and premaxilla. Besides these two modules, the within-region trait correlations observed for all landmark-only analyses were considerably lower than for the complete shape data analyses, and extremely low for the nasal, frontoparietal, and maxilla (dorsal) (0.14, 0.23, and 0.27, respectively). Between-region trait correlations were also lower, but only slightly, so that in some cases the between-region correlations were actually higher than the within-region correlations (e.g., occipital within: 0.57, frontoparietal within: 0.23, between: 0.32). Many within-region and between-region trait correlations were below 0.3, and the effect of this was that most regions did not stand out as independent modules, nor could they be grouped with other regions as being strongly integrated. Most of the cranium therefore appeared unintegrated, with the exception of the occipital and facial modules. The suspensorium module identified in the complete dataset analyses could not be investigated here, as two of these regions (the quadratojugal regions) were not present in the landmark-only data.
Rerunning EMMLi analysis excluding all specimens with any "negligible regions" revealed a near-identical pattern of trait integration to that recovered from the entire dataset for phylogenetically-corrected data (Fig. S5). For this analysis, the within-region trait correlations for the variably present regions were either unchanged, or only marginally different (Table S14) to analysis using the complete dataset. The quadratojugal regions exhibited no change in within-region trait correlations, and the largest difference was a 0.05 decrease in within-region trait integration (stapes and sphenethmoid (ventral) regions). We recovered a modular structure from this analysis that was identical to the model recovered using the complete dataset. The use of the "negligible region" method therefore did not artificially exaggerate within-region trait integration for regions that are variably present.
All three CR analyses of the complete shape data (using uncorrected, phylogenetically-, and allometry-corrected datasets) recovered significant modular signal (CR = 0.52, 0.48, and 0.49, respectively, P = 0.001 for all). When pairwise CR values were investigated, the pattern of trait relationships was extremely similar to the results from EMMLi for all three analyses of the complete dataset (Fig. S3). For all three analyses, the strongest covariances between regions were those within the identified occipital module (CR = 0.97-0.99), premaxilla module (CR = 0.95-0.97), and suspensorium module (CR = 0.78-0.96) (Tables S15-S17). The maxilla module also showed strong covariation (CR = 0.73-0.8), although the maxilla (dorsal) covaried more strongly with the occipital than the maxilla (ventral) using phylogeneticallycorrected data and also more strongly with the occipital condyle using uncorrected data.
Conducting CR for the landmark-only datasets (uncorrected, phylogenetically-, and allometry-corrected) revealed weaker, but still significant, modular signal compared with the full datasets (CR = 0.82, 0.77, and 0.78, respectively, P = 0.001 for all). The pattern of integration observed for these datasets was similar to the landmark-only EMMLi analyses; there was strong covariation within the occipital module (CR = 1.14, phylogenetically-corrected CR = 1.12, allometry-corrected CR = 1.17, respectively) and strong covariations between the four regions comprising the "facial" module (CR = 1.04-1.34, phylogenetically-corrected CR = 0.95-1.31, allometry-corrected CR = 1.03-1.17), although the two maxilla regions also covaried strongly with some other regions ( Fig. S4; Tables S18-S20). Species with and without a feeding larval stage exhibited very similar patterns of trait integration, although specimens without a feeding larval stage had very slightly more strongly integrated crania, as recovered from both phylogenetically informed EMMLi and CR analyses ( Fig. S6; Tables S21-S26). However, subsampling the "feeding larvae" group down to 39 specimens and finding the average of running phylogenetically informed EMMLi and CR analyses iteratively 100 times eliminated the very slight difference between the two groups. Z-scores for the three datasets (taxa with nonfeeding larvae, taxa with feeding larvae [N = 124] and subsampled taxa with feeding larvae) were not significantly different (pairwise P-values = 0.72-1.00).

OSSIFICATION SEQUENCE
There was a significant relationship observed between ossification sequence rank and magnitude of integration for individual cranial bones (Spearman's rho = 0.59, P = 0.027) (Fig. 4). Median Ossification Sequence Position Within Module Correlation (ρ) Figure 4. Influence of ossification sequence timing on integration, using the median rank position within the frog cranial ossification sequence for each cranial bone (Weisbecker and Mitgutsch 2010). The "suspensorium" module was split into its constituent bones (quadratojugal, squamosal, and pterygoid) and the stapes was excluded. Spearman's rank correlation between ossification sequence timing and integration (within-module correlation) revealed the relationship was significant (Spearman's rank correlation rho = 0.59, P = 0.027).

Discussion
Trait integration is of fundamental interest in reconstructing the evolution of phenotype and of phenotypic diversity. When integrated traits experience divergent selection pressures, their covariation can hinder the ability of each trait to evolve toward its respective optimum. This constraint may drive the fragmentation of sets of integrated traits into smaller, functionally linked units comprising traits with aligned selection pressures and with reduced linkages to other sets of traits, thus freeing modules from the constraints of integration with functionally unrelated traits. This modular organization may thus be expected to facilitate evo-lutionary diversification. Although it is the variational, rather than the evolutionary, relationships among traits that shape evolvability and influence response to selection, evolutionary modularity informs on the outcome of those processes.
Variational and evolutionary modularity often correspond, for example, in caecilians (Marshall et al. 2019;Bardua et al. 2019b) and in salamanders Fabre et al. 2020), suggesting that the genetic, developmental, and functional interactions of traits that shape phenotypic integration are replicated in patterns of trait evolution. As discussed above, many studies have suggested that evolutionary modularity promotes diversification of form. Studies have also found that high levels of phenotypic and evolutionary integration within modules are associated with limited morphological diversity (disparity) and/or evolutionary rate of those modules (Goswami and Polly 2010;, or conversely, that strong integration of traits is associated with higher disparity (Claverie and Patek 2013;Parr et al. 2016;Randau and Goswami 2017;Navalón et al. 2020). Many other studies suggest there may be no simple relationship between integration and either morphological disparity or rate at either level (Drake and Klingenberg 2010;Goswami et al. 2014;Marshall et al. 2019;Watanabe et al. 2019;Bardua et al. 2019b). Improving our understanding of the relationship between modularity in trait evolution and morphological diversification ultimately requires more (and more diverse) empirical data on modularity and morphological evolution within species and across clades. Here, we have expanded the study of evolutionary modularity to a taxonomically, developmentally, ecologically, and morphologically diverse tetrapod clade: frogs.

HIGHLY MODULAR ANURAN CRANIA
We find anuran crania evolve in a highly modular manner, exhibiting 13 distinct evolutionary modules. Our results contrast with previous analyses recovering only weak support for phenotypic or evolutionary modularity across frog crania (Simon and Marroig 2017;Vidal-García and Keogh 2017), and we find no support for the functional models investigated in these studies. However, those studies investigated modularity at a finer taxonomic scale, sampling at subgenus or family level, respectively, with likely greater similarity (and less variation) in function and ecology than is sampled here. Furthermore, differences in data type and study design hinder direct comparisons, with our study the first to investigate highly modular hypotheses with dense morphometric data. Our analyses suggest a complex model, where osteological units were mainly recovered as either distinct modules (premaxilla and maxilla) or split into multiple modules (e.g., dorsal and ventral surfaces of the sphenethmoid). One multi-region module, the "suspensorium," comprises three bones (quadratojugal, squamosal, and pterygoid) and corresponds to a similar multi-bone module in salamanders and caecilians (Bardua et al. 2019b;Fabre et al. 2020). Function (adult feeding) thus appears to be a more proximal driver of evolutionary integration than development (i.e., tissue origin) is for this region.
The suspensory apparatus braces and suspends the jaws against the neurocranium for autosystylic jaw suspension and has been long recognized as an important functional unit in frogs (Trueb 1973). The squamosal and pterygoid act to protect and strengthen the cartilaginous system (including the quadrate), and the quadrate (which ossifies in some species) is often invaded by the ossified quadratojugal (Trueb 1973). The suspensorium likely acts as a functional unit of selection with studies noting coordinated change in this unit. For example, the suspensorium has been found to be more vertically oriented in miniaturized anuran species, in which a more rostrally positioned jaw joint articulation leads to a smaller gape size (Yeh 2002). Also, direct-developing species have been shown to have accelerated appearance of the jaw and suspensorium bones through their cranial bone formation (Hanken et al. 1992;Kerney et al. 2007;Vassilieva 2017), demonstrating coordinated changes of the suspensorium module due to heterochronic repatterning. The range of feeding mechanisms across frogs associated with different environments, for example, aquatic suction feeding (Carreño and Nishikawa 2010;Fernandez et al. 2017) and hydrostatic elongation in many fossorial (and some terrestrial) species (e.g., Trueb and Gans 1983;Nishikawa et al. 1999), may have been facilitated by an integrated jaw suspensorium module, allowing these mechanisms to arise convergently multiple times (Nishikawa 2000).
Life history strategies are extremely varied across frogs, and differences in strategy (in particular, the presence or absence of a larval stage that interacts with its environment) may be expected to exert an influence on the structure of evolutionary modules. Frogs undergo major morphological restructuring of the cranium through metamorphosis (Rose and Reiss 1993), either inside (direct-developing species; Callery and Elinson 2000;Ziermann and Diogo 2013) or outside (biphasic species) the egg. Unlike frogs with nonfeeding larvae, frogs with a feeding larval stage need to adapt to two environments, each with unique functional pressures (Reiss 2002). Species with feeding larvae thus may have greater autonomy across stages and consequently fewer developmental constraints, which may be expected to manifest in greater variation in patterns of phenotypic integration and thus greater evolutionary modularity. In contrast, taxa with nonfeeding larvae may be more developmentally constrained, resulting in greater similarity in patterns of phenotypic integration across species and thus stronger evolutionary integration. Surprisingly, despite these expectations, we do not find evolutionary integration to differ substantially between taxa with and without a feeding larval stage. A recent study on salamander crania  found that taxa undergoing complete metamorphosis (biphasic and direct-developing) displayed more modular crania than did pedomorphic taxa. Thus, metamorphosis in salamanders may represent an adaptive mechanism through which developmental processes shaping morphological variation are decoupled across different life history stages (Moran 1994;Fabre et al. 2020). It may therefore be the process of metamorphosis, and not the interaction with a larval environment, that results in a decoupling of life history stages. Furthermore, life history may be a bigger influence on cranial morphological evolution of salamanders than frogs. Ecology (habitat) has been found to be correlated with skull shape in frogs (Paluh et al. 2020) and may be a primary factor shaping cranial morphological evolution, in contrast to salamanders, where ecology is a relatively minor factor compared to life history . Thus, the influence of life history on evolutionary modularity may be confounded by ecological influences, and the process of metamorphosis may be a greater driver of trait fragmentation than the possession of life stages with distinct ecological pressures.
Developmental associations of traits are often hypothesized to drive phenotypic and evolutionary integration (e.g., Simon and Marroig 2017). We found no evidence of strong relationships among regions sharing contributions from CNC streams or with similar ossification sequence ranks. However, frogs (evidenced by Xenopus) exhibit a considerably different pattern of CNC contributions to cranial elements compared with salamanders and other vertebrates (Hanken and Gross 2005;Gross and Hanken 2008;Piekarski et al. 2014). The pattern of CNC contributions therefore appears evolutionarily labile (Hanken and Gross 2005) and may have been changing throughout anuran evolution, which is not captured in the CNC model of modularity.
An additional developmental influence on the strength of cranial integration is the timing of bone ossification. Cranial ossification sequences are highly evolvable across Anura (Weisbecker and Mitgutsch 2010), with large variation in the relative timing of cranial bone ossification, although early and late ossifying bones are generally the most conservative. Contrary to our expectation, ossification sequence rank correlated positively with the strength of evolutionary integration within cranial bones. The latest-ossifying bones showed the highest evolutionary integration and are those that are lost in many taxa, possibly due to pedomorphosis (Weisbecker and Mitgutsch 2010). The intermediateossifying bones, which are generally the least conserved in ossification sequence rank (Weisbecker and Mitgutsch 2010;Vassilieva 2017), displayed the broadest range of magnitudes of integration. Direct and biphasic developers differ most in the relative timing of intermediate-ossifying bones (Vassilieva 2017), so the wide range of timings of intermediate-ossifying bones across taxa may therefore have weakened the relationship between ossification and integration across these bones. Regardless, the significant relationship between ossification sequence rank and cranial integration across frogs, compared with the lack of influence of life history on integration, suggests at least some aspects of development do significantly influence the coordinated evolution of traits in the anuran cranium.

EVOLUTIONARY MODULARITY
A large, multi-region module associated with feeding and prey capture appears to be common across lissamphibian lineages (Bardua et al. 2019b;Fabre et al. 2020), comprising the lateral and articular surfaces of the quadrate/quadratojugal, squamosal, and, in frogs and biphasic salamanders, pterygoid. This tight functional feeding unit is thus replicated across lissamphibians and selection on this module is likely to be a primary and persistent factor shaping morphological evolution across lissamphibians. The model of evolutionary modules in the anuran cranium supported here is remarkably similar in number and pattern to that recovered recently from analyses of both phenotypic and evolutionary modularity in caecilian (Marshall et al. 2019;Bardua et al. 2019b) and salamander Fabre et al. 2020) crania, despite substantial morphological differences among these clades. Replication of patterns of cranial modularity from static to evolutionary level across caecilians and salamanders has already been demonstrated (Marshall et al. 2019;Bardua et al. 2019b;Bon et al. 2020;Fabre et al. 2020), and now our results suggest that these evolutionary modules are shared across Lissamphibia. Lissamphibian crania also display more modular evolution than amniote crania, which may be partly attributed to their more frequent losses and regains of cranial elements as well as their greater number of cranial elements. Studies of amniote crania using similar morphometric data have identified numbers of evolutionary modules ranging from seven (Aves; , to nine (snakes; Watanabe et al. 2019) or 10 (dingos and dogs; Parr et al. 2016;lizards;Watanabe et al. 2019). Comparison of patterns of trait evolutionary integration recovered from these studies reveals some striking similarities. For example, the occipital region and the jaw joint region are both strongly integrated across caecilians (Bardua et al. 2019b), salamanders , birds , squamates , crocodylomorphs and nonavian dinosaurs , and, as demonstrated here, frogs. However, caecilians, salamanders, and frogs exhibit the strongest concordance in their patterns of cranial integration, suggesting a possible amniote-amphibian divergence in the structure of evolutionary modules of the cranium.
Other studies of tetrapod cranial modularity have recovered weaker or less complex modular structures, including six-module models (mammals, Cheverud 1982;Goswami 2006;Goswami and Polly 2010;Goswami and Finarelli 2016), twomodule models (Alpine newt, Ivanović Klingenberg and Marugán-Lobón 2013). However, many of these studies used many fewer two-dimensional or three-dimensional landmarks and thus necessarily tested less-complex patterns of modularity, hindering comparison with studies incorporating surface semilandmarks.

MORPHOLOGICAL EVOLUTION
With these new data, we return to the question of whether evolutionary modularity has promoted the diversification of the frog cranium. Our results demonstrate that anuran cranial modules vary widely in terms of morphological disparity, suggesting that modular variation has indeed promoted the mosaic evolution of the anuran cranium. However, the magnitude of evolutionary integration within each cranial module does not show a strong or consistent relationship with disparity, similar to findings from caecilian (Bardua et al. 2019b), salamander , and squamate  crania. Other studies have found high evolutionary integration to correspond with either higher or lower disparity or rates of evolution in different clades (Goswami and Polly 2010;Claverie and Patek 2013;Parr et al. 2016;Randau and Goswami 2017;Navalón et al. 2020). The observed relationship between evolutionary integration and morphological evolution depends on many factors, including the alignment of the direction of selection with the path that phenotypic integration facilitates in morphospace and the evolutionary stability of both of those attributes, as well as biases in fully capturing morphological diversity due to extinction. The relationship evolutionary integration has with morphological disparity is therefore likely complex, with patterns for each module reflecting their diverse selection histories, rather than following a single simple relationship across a larger structure.
Finally, on a more practical note, comparing landmark-only analyses to our analyses of more complete descriptions of cranial shape highlights the benefits of implementing a surfacebased approach for investigating modularity across frog crania. The potential for landmarks to amplify between-region trait correlations for adjacent regions and understate within-region correlations Marshall et al. 2019;Bardua et al. 2019b) may result in the cranium appearing less modular, and more integrated, than when incorporating curve and surface information. In addition, crucially, landmarks fail to capture seven cranial regions in anurans, because these are variably present across the clade. This exclusion of nearly one third of the regions from modularity analyses would prevent patterns of integration from being comprehensively explored across the whole cranium. We assigned "negligible regions" in taxa with absent regions (Bardua et al. 2019a,b), and our comparison of analyses including and excluding these taxa demonstrates that these "negligible regions" do not affect the observed pattern of integration, as results are near-identical with and without taxa lacking those regions. Thus, our high-dimensional surface-based approach allows a more representative quantification of trait integration and cranial morphological evolution across frogs.

Conclusion
We implemented a high-density approach to quantify anuran crania across the entire clade, with extensive sampling of taxonomy and morphology. We identified 13 evolutionary modules, which replicate patterns recovered across caecilians and salamanders, but are less similar to those observed in amniotes, suggesting a possible lissamphibian-amniote divergence in cranial modularity. Surprisingly, contrasting life histories have no impact on evolutionary integration of the anuran skull, in contradiction to the hypothesis that a feeding larval stage may have arisen as a mechanism to fragment evolutionary associations among regions (although it may well apply to different stages). Timing of ossification shows a significant correlation with evolutionary integration, with late-ossifying bones displaying higher magnitudes of integration than earlier-ossifying bones. Cranial modules display a wide range of morphological disparities, but magnitude of within-module integration did not correspond significantly with module disparity, supporting recent studies suggesting a complex relationship between evolutionary integration and morphological diversification. Finally, we have illustrated the utility of this highdimensional approach, for quantifying cranial morphology across every frog family, allowing for incorporation of regions that were absent in some specimens. Applying a similar high-dimensional approach to other extremely diverse clades, many of which are currently undersampled, will facilitate the study of phenotypic and evolutionary integration in these clades and allow broader comparisons across different scales of analysis and across the tree of life.

CONFLICT OF INTEREST
The authors declare no conflict of interest.

AUTHOR CONTRIBUTIONS
CB and AG designed the study. CB, ELS, KD, MB, and DCB generated and reconstructed the scans. CB undertook data collection, performed analyses, and primarily wrote the manuscript. All authors provided ideas and discussion, contributed to the writing of the manuscript, and read and approved the final version.

Supporting Information
Additional supporting information may be found online in the Supporting Information section at the end of the article. Figure S1. Landmark and semilandmark data displayed on specimen. Figure S2. Morphospace of the first two phylogenetic principal components for anuran crania. Figure S3. Network graphs displaying results from EMMLi and CR analyses. Figure S4. Network graphs displaying results from landmark-only EMMLi and CR analyses. Figure S5. Network graph for EMMLi analysis excluding specimens with absent regions. Figure S6. Effect of life history on trait integration.  Table S7. Results of EMMLi analysis using uncorrected data. Table S8. Results of EMMLi analysis using phylogenetically-corrected data. Table S9. Results of EMMLi analysis using allometry-corrected data. Table S10. Results of EMMLi analysis using subsampled data. Table S11. Results of EMMLi analysis using uncorrected landmark-only data. Table S12. Results of EMMLi analysis using phylogenetically-corrected landmark-only data. Table S13. Results of EMMLi analysis using allometry-corrected landmark-only data. Table S14. Results of phylogenetic EMMLi analysis excluding specimens with absent regions. Table S15. Covariance Ratio results using uncorrected data. Table S16. Covariance Ratio results using phylogenetically-corrected data. Table S17. Covariance Ratio results using allometry-corrected data. Table S18. Covariance Ratio results using uncorrected landmark-only data. Table S19. Covariance Ratio results using phylogenetically-corrected landmark-only data. Table S20. Covariance Ratio results using allometry-corrected landmark-only data. Table S21. Results of phylogenetic EMMLi analysis for species possessing a feeding larval stage. Table S22. Results of subsampled phylogenetic EMMLi analysis for species possessing a feeding larval stage. Table S23. Results of phylogenetic EMMLi analysis for species lacking a feeding larval stage. Table S24. Results of phylogenetic CR analysis for species possessing a feeding larval stage. Table S25. Results of subsampled phylogenetic CR analysis for species possessing a feeding larval stage. Table S26. Results of phylogenetic CR analysis for species lacking a feeding larval stage.