# Size, microhabitat, and loss of larval feeding drive cranial diversification in frogs

Carla Bardua, Anne-Claire Fabre, Julien Clavel, Margot Bon, Kalpana Das, Edward L Stanley, David C Blackburn, Anjali Goswami

► To cite this version:

Carla Bardua, Anne-Claire Fabre, Julien Clavel, Margot Bon, Kalpana Das, et al.. Size, microhabitat, and loss of larval feeding drive cranial diversification in frogs. Nature Communications, 2021, 12 (1), 10.1038/s41467-021-22792-y. hal-03231761

HAL Id: hal-03231761

https://univ-lyon1.hal.science/hal-03231761

Submitted on 26 Nov 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

ARTICLE

https://doi.org/10.1038/s41467-021-22792-y

OPEN

<figure><img src="image_3.png" /></figure>

Check for updates

# Size, microhabitat, and loss of larval feeding drive cranial diversification in frogs

Carla Bardua ID 1,2, Anne-Claire Fabre1,3, Julien Clavel1,4, Margot Bon1, Kalpana Das ID 5, Edward L. Stanley ID 6, David C. Blackburn6 & Anjali Goswami ID 1✉

Habitat is one of the most important factors shaping organismal morphology, but it may vary across life history stages. Ontogenetic shifts in ecology may introduce antagonistic selection that constrains adult phenotype, particularly with ecologically distinct developmental phases such as the free-living, feeding larval stage of many frogs (Lissamphibia: Anura). We test the relative influences of developmental and ecological factors on the diversification of adult skull morphology with a detailed analysis of 15 individual cranial regions across 173 anuran species, representing every extant family. Skull size, adult microhabitat, larval feeding, and ossification timing are all significant factors shaping aspects of cranial evolution in frogs, with late-ossifying elements showing the greatest disparity and fastest evolutionary rates. Size and microhabitat show the strongest effects on cranial shape, and we identify a “large size-wide skull” pattern of anuran, and possibly amphibian, evolutionary allometry. Fossorial and aquatic microhabitats occupy distinct regions of morphospace and display fast evolution and high disparity. Taxa with and without feeding larvae do not notably differ in cranial morphology. However, loss of an actively feeding larval stage is associated with higher evolutionary rates and disparity, suggesting that functional pressures experienced earlier in ontogeny significantly impact adult morphological evolution.

<footnote n="1">¹Department of Life Sciences, Natural History Museum, London, UK. ²Department of Genetics, Evolution & Environment, University College London, London, UK. ³Paläontologisches Institut und Museum, Universität Zürich, Zürich, Switzerland. ⁴Univ Lyon, Université Claude Bernard Lyon 1, CNRS, ENTPE, UMR 5023 LEHNA, F-69622, Villeurbanne, France. ⁵Museum für Naturkunde, Leibniz Institut für Evolutions und Biodiversitätsforschung, Berlin, Germany.</footnote>

<footnote n="6">⁶Department of Natural History, Florida Museum of Natural History, University of Florida, Gainesville, FL, USA. ✉email: a.goswami@nhm.ac.uk</footnote>
Adaptations to habitats are responsible for some of the most striking examples of phenotypic evolution, driving profound shifts in morphology¹⁻⁶, in shape diversity ('disparity')⁷ and in rate of shape evolution¹,⁸⁻¹¹. Adaptations to similar habitats also result in phenotypic convergence when responses to selection produce similar shifts in phenotypic optima across diverse clades²,¹²⁻¹⁵, as in the body forms of pelagic fishes, ichthyosaurs and whales. Habitat exploitation clearly plays a major role in shaping morphological evolution and macro-evolutionary dynamics across organisms. However, an organism's ability to specialise for an adult habitat may be constrained by adaptations to conditions experienced at earlier developmental stages. Species with morphologically and ecologically distinct developmental stages may experience confounding and competing selective pressures throughout ontogeny which may limit their ability to evolve specific morphologies¹⁶⁻¹⁸. The potential constraints of divergent selection at different developmental stages may drive selection for a simpler life history, i.e. one with no abrupt ontogenetic change in morphology or niche. Accordingly, the loss of a free-living, actively feeding larval stage has been credited with promoting morphological novelty¹⁹,²⁰ and driving increased rates of morphological evolution across salamanders and caecilians¹¹,²¹,²². Decoupling adult and embryonic form via an ecologically distinct larval stage may circumvent this potentially antagonistic selection by creating developmental modules that allow adult phenotype to evolve independently of the constraints of earlier life stages²³⁻²⁵. Investigating the relative influences of adult ecology and potential constraints of ontogenetic form on morphological evolution is complex given that the two are inherently linked. Developmental changes may act as exaptations that promote the exploitation of a habitat²⁶, or, alternatively, adaptation to environment may drive shifts in developmental strategy²⁷. Furthermore, both factors likely interact with other traits, such as body size or developmental timing²⁸,²⁹, and their influence may vary across morphological structures. Untangling the roles of ecology and development on morphological diversification thus necessitates a macroevolutionary and multifaceted approach.

With more than 7300 extant species, frogs (Lissamphibia: Anura) are one of the most diverse orders of vertebrates and show remarkable variation in both life history and ecology. Although frogs are predominantly generalist predators, they are ecologically diverse, including species specialised for fossorial, aquatic, arboreal and terrestrial microhabitats³⁰⁻³², which has facilitated their widespread and near-global distribution³³. Because anurans vary widely in body size, ranging from a snout–vent lengths of ~7 to 320 mm, there are ample opportunities for exploring the relationship between microhabitat specialisation and body size, which across organisms affects many ecological traits³⁴, and the evolution of shape³⁵. Life history strategies in frogs range from direct development (without a free-living, feeding larval phase) to fully biphasic, with free-living feeding larval phases (i.e. tadpoles) specialised to a wide range of aquatic microhabitats and diets that differ from adult frogs. Biphasic frogs undergo radical morphological changes during metamorphosis; the magnitude of ontogenetic change is less in direct-developing species³⁶,³⁷ though the extent of reorganisation of internal anatomy remains poorly explored. All frogs thus experience some degree of morphological transformations through ontogeny, but species vary extensively in the degree of interaction with their environment as larvae. As feeding is one of the most critical functions of the cranium, larval feeding likely imposes strong selection pressures on the developing cranium. In particular, taxa with feeding larvae require distinct functional adaptations for their larval and adult environments, given the distinct changes in feeding mechanisms across the two life stages³³,³⁸⁻⁴⁰. For this reason, the presence of an

actively feeding larval stage has been hypothesised to be an ontogenetic constraint on the adult morphology, whereas the loss of an actively feeding larval stage has been suggested to promote morphological novelty¹⁹,²⁰,⁴¹,⁴² and drive increased rates of morphological evolution¹¹,²¹,²². This ontogenetic transition in function may also differentially impact structures based on their relative timing of development, for example with greater constraints on early-ossifying bones⁴³.

The developmental and ecological complexity of frogs provide an ideal opportunity to test competing hypotheses on the relative influences of ecological and developmental factors on the diversification of frog morphology. Some previous studies have addressed how larval morphology impacts adult morphology in frogs, primarily with data on body size⁴⁴⁻⁴⁶, which is easily comparable across stages and taxa, or with studies focused on one or a few anuran families⁴⁴,⁴⁷,⁴⁸. While there is a correlation between larval and adult body size across all frogs⁴⁵, the strength of this relationship varies across families⁴⁴, and there is no significant correlation in rates of evolution of adult and larval size⁴⁵. Studies capturing aspects of morphology other than body size support decoupling of larval and adult evolution but are limited both in taxonomic breadth and morphological data⁴⁷,⁴⁸. Despite long-standing hypotheses of decoupled evolution between life history stages and its importance to understanding the prevalence of complex life cycles in animals, there are no studies that test this hypothesis with data that captures both the complexity of adult morphology and its diversity across frogs. Moreover, no studies have explicitly compared the relative influence of developmental effects to other factors, such as adult ecology. If the anuran larval stage acts as a developmental module that fully decouples larval and adult phenotypic evolution, then we would expect no differences in adult morphological evolution between taxa with different life history strategies. If so, then we would expect factors such as adult microhabitat to be more significant influences on adult morphological evolution. In contrast, if larval and adult evolution is correlated, then the compound effects of divergent selection pressures acting on larval and adult stages are expected to influence, and potentially constrain, phenotypic evolution of the adult form in taxa with feeding larvae. Furthermore, developmental and adult ecological factors may interact, meaning that both adult microhabitat and adult morphological evolution might be influenced by the degree to which larval forms interact with their environment.

Given the extensive shape variation⁴⁹,⁵⁰ and the potential heterogeneity of ecological and developmental influences across cranial regions, the adult skull provides an ideal system for testing these hypotheses. As in all vertebrates, the anuran cranium comprises multiple osteological units with different functions, embryonic origins and types and timings of ossification⁴³,⁵¹⁻⁵³. Anurans have the added complexity that many cranial bones are variably present between species⁵⁰,⁵⁴,⁵⁵, and their crania likely exhibit substantial mosaic evolution⁵⁶⁻⁶³. To tackle this challenge, we quantify anuran cranial complexity using a high-density surface morphometric approach that represents a significant advance over previous work which used a limited number of landmarks and excluded structures that are variably present, including a number of key functional elements such as the neopalatine, sphenethmoid, vomer, quadratojugal and stapes³,⁶⁴. With these data, we reconstruct the morphological evolution of the anuran skull, spanning the full range of frog ecological, developmental and morphological diversity (Fig. 1). We estimate the relative influences and interactions of adult microhabitat, skull size, larval feeding, and ossification timing on the morphology and diversification of the frog cranium and provide an empirical analysis of the larval-adult decoupling hypothesis using 3D shape data across the full phylogenetic breadth of frogs.
<figure><img src="image_4.png" /><figcaption>**Fig. 1** Cranial shape diversity across Anura. **A** Defined cranial regions displayed on one anuran skull from every family across the phylogeny (lateral aspect, skulls not to scale). Phylogeny for the 173 species in this study modified from Jetz and Pyron¹¹³ with the addition of *Thaumastosaurus gezei*. Branch lengths are scaled to time, with alternate shaded bands indicating periods of 50 million years (outer extreme at Recent). For details on specimens visualised here, see Supplementary Data 1. Landmarks and semilandmarks in **B** lateral, **C** dorsal and **D** ventral views, shown on *Adenomus kelaartii* (FMNH: Amphibians and Reptiles:1580). Points are coloured according to the 15 regions defined in this study. Regions are as follows: Occ (light purple) occipital region, FP (black) frontoparietal, Oj (red) quadratojugal, Max (orange) maxilla, Na (green) nasal, Neo (hot pink) Neopalatine, Otic (lime green) otic region, Pm (pale green) premaxilla, PS (purple) parasphenoid, Pt (light blue) pterygoid, SphD (brown) sphenethmoid (dorsal surface), SphV (light pink) sphenethmoid (ventral surface), Sq (blue) squamosal, St (yellow) stapes, Vo (grey) vomer.</figcaption></figure>
In this work, we show that skull size, adult microhabitat, larval feeding and ossification timing are all significant factors shaping aspects of cranial evolution in frogs, with late-ossifying elements showing the greatest disparity and fastest evolutionary rates. We identify a 'large size-wide skull' pattern of anuran, and possibly amphibian, cranial evolutionary allometry driven by structural constraints and feeding adaptations. Skull size and microhabitat consistently show the strongest effects on cranial shape, with fossorial, semi-fossorial and aquatic adult microhabitats occupying distinct regions of morphospace and displaying increased rates of evolution and high disparity, likely reflect specialisations of adult feeding and hearing functions. Taxa with and without feeding larvae do not notably differ in their cranial morphology, supporting partial adult-larval decoupling. However, loss of an actively feeding larval stage is associated with higher rates of evolution and disparity in most cranial regions, suggesting that functional pressures experienced earlier in ontogeny still significantly impact patterns of adult morphological evolution in frogs. Our analysis thus provides a comprehensive assessment of the relative effects of components of larval and adult ecology, and their interaction, on the morphological evolution of frogs, with relevance for understanding the evolution, and consequences, of complex life cycles across the tree of life.

## Results

**Cranial shape.** Principal components (PCs) analysis of the entire adult cranium (Fig. 2 and Supplementary Table 1) showed that the predominant shape changes along PC1 (28% of variance) were the height of the cranium, the anteroposterior translation of the quadratojugal, the lateral expansion of the frontoparietal, the reduction of the maxillary arcade and the shortening of the anterior process of the squamosal (Supplementary Note 1,

Supplementary Figs. 1, 2 and see Fig. 2 for the cranial regions). PC2 (14%) represented variation associated with the height of the maxilla, a larger nasal, and an anteroposteriorly compressed frontoparietal, as well as a taller, wider cranium (Supplementary Fig. 3).

The first two PCs separated fossorial (including semi-fossorial) and aquatic (including semi-aquatic) taxa, whereas terrestrial, arboreal and semi-arboreal species overlapped with both (Fig. 2). Taxa with and without an actively feeding larval stage overlap entirely in cranial morphospace (Supplementary Fig. 4). Phylogenetic regression of cranial shape on centroid size, microhabitat and larval feeding supported all three factors as significantly associated ($p < 0.01$) with shape of the entire skull and most cranial elements (Table 1 and Supplementary Table 2). Presence of an actively feeding larval stage consistently displayed the lowest effect size of the three factors, was not significantly associated with the shape of the nasal, stapes or dorsal sphenethmoid, and was only marginally significant ($0.05 > p > 0.01$) for the otic, maxilla and quadratojugal. Skull size and adult microhabitat showed similar effect sizes for most elements and were significantly associated with shape of all elements, with the exception of adult microhabitat showing no significant association with the neopalatine and only a marginally significant one with the dorsal sphenethmoid (Table 1).

There were significant interactions among size, adult microhabitat and larval feeding for the entire skull and most elements except for the nasal, neopalatine and dorsal sphenethmoid ($p > 0.05$) and marginally significant interactions ($0.05 > p > 0.01$) for the paraspheoid, premaxilla and quadratojugal (Supplementary Table 2). Size and adult microhabitat had highly significant interactions ($p < 0.01$) for most elements (Supplementary Fig. 4), as did adult microhabitat and larval feeding (Supplementary Fig. 5), whereas size and larval feeding showed a highly significant interaction in only six elements.

<figure><img src="image_5.png" /></figure>

**Fig. 2** *Phylomorphospace of anuran skull variation.* Each specimen's location is coloured by microhabitat and the symbols indicate larval feeding. The maximum morphospace occupation for fossorial (including semi-fossorial) and aquatic (including semi-aquatic) species is illustrated by shading in the respective regions brown and blue. Specimen meshes are included to demonstrate the general shape differences between fossorial (*Nasikabatrachus sahyadrensis*) and aquatic (*Pipa parva*) species, as well as the specimens with the largest (*Conraua goliath*) and smallest (*Stumpffia pygmaea*) crania. Extreme shapes along both axes are visualised (in lateral aspect) by deforming the shape data along each axis and are coloured by cranial region (see Supplementary Fig. 1). Estimated cranial morphology at minimum and maximum cranial size is also presented, shown in dorsal (top) and lateral (bottom) aspects. Landmark data were mirrored for visualisation purposes only.
<table><caption>Table 1 Region-specific analyses of shape.</caption><thead><tr><th>Cranial region</th><th>Disparity (Procrustes variance) (× 10<sup>-5</sup>)</th><th>Evolutionary rate, σ<sup>2</sup>mult, (× 10<sup>-9</sup>)</th><th>Ossification sequence rank</th><th>Evolutionary allometry (SES)</th><th>Habitat (SES)</th><th>Larval feeding (SES)</th></tr></thead><tbody><tr><td>Frontoparietal</td><td>1.04</td><td>0.349</td><td>2</td><td>7.85***</td><td>8.17***</td><td>4.06***</td></tr><tr><td>Maxilla</td><td>1.35</td><td>3.56</td><td>6</td><td>5.67***</td><td>7.22***</td><td>2.13*</td></tr><tr><td>Nasal</td><td>1.07</td><td>1.91</td><td>7</td><td>7.21***</td><td>7.77***</td><td>0.30</td></tr><tr><td>Neopalatine</td><td>1.91</td><td>34.8</td><td>12</td><td>5.47***</td><td>1.01</td><td>2.72**</td></tr><tr><td>Occipital</td><td>0.72</td><td>2.39</td><td>3</td><td>7.50***</td><td>8.34***</td><td>3.72***</td></tr><tr><td>Otic</td><td>1.13</td><td>3.87</td><td>4</td><td>8.18***</td><td>5.66***</td><td>1.67*</td></tr><tr><td>Parasphenoid</td><td>0.60</td><td>3.84</td><td>1</td><td>7.51***</td><td>7.01***</td><td>3.49***</td></tr><tr><td>Premaxilla</td><td>0.89</td><td>5.16</td><td>5</td><td>5.41***</td><td>7.57***</td><td>3.30**</td></tr><tr><td>Pterygoid</td><td>1.20</td><td>7.86</td><td>10</td><td>7.06***</td><td>6.27***</td><td>3.19***</td></tr><tr><td>Quadratojugal</td><td>2.92</td><td>27.8</td><td>11</td><td>9.21***</td><td>6.81***</td><td>2.31*</td></tr><tr><td>Sphenethmoid (dorsal)</td><td>1.29</td><td>12.7</td><td>13</td><td>5.66***</td><td>2.10*</td><td>0.08</td></tr><tr><td>Sphenethmoid (ventral)</td><td>1.96</td><td>2.33</td><td>13</td><td>8.38***</td><td>7.28***</td><td>3.67**</td></tr><tr><td>Squamosal</td><td>1.81</td><td>5.03</td><td>8</td><td>6.91***</td><td>9.92***</td><td>3.63***</td></tr><tr><td>Stapes</td><td>1.86</td><td>17.5</td><td>NA</td><td>10.84***</td><td>12.65***</td><td>1.35</td></tr><tr><td>Vomer</td><td>0.99</td><td>5.37</td><td>9</td><td>7.78***</td><td>4.53***</td><td>3.83***</td></tr></tbody></table>

Disparity is measured as Procrustes variance scaled by semi/landmark number. Evolutionary rate estimates include error. Median ossification sequence rank taken from Weisbecker$^{43}$. Effect sizes (SES) and significances for evolutionary allometry (size), microhabitat and larval feeding from Type II phylogenetic regressions of cranial shape with all three predictors and all interactions using the MCC tree. All tests are one-sided, with significances based on permutations (n = 999). Interactions, Pillai's test results, Pagel's lambda and results for sample of 100 trees are detailed in Supplementary Table 2. Significance of results for evolutionary allometry, influence of microhabitat, and influence of larval feeding is as follows: *p* values significant at the following alpha levels: *≤0.05, **≤0.01, ***≤0.001.

Replicating these analyses across a distribution of 100 phylogenetic trees demonstrates that these results are also robust to phylogenetic uncertainty (Supplementary Table 2).

Given the significant influence of size on cranial shape and the significant interaction of microhabitat with size for most cranial elements, we further investigated the variation associated with change in skull size across anurans. Larger species had wider, more triangularly-shaped skulls with relatively smaller cranial vaults (i.e. frontoparietal and parasphenoid) and occipital regions (Fig. 1) and a more horizontally-oriented suspensory apparatus (pterygoid, squamosal and quadratojugal). Aquatic species had larger crania than taxa from other microhabitats (Supplementary Fig. 4 and Supplementary Note 1) although size ranges overlapped across all microhabitats. Removing allometric size from shape data prior to morphospace analysis increased the separation of fossorial and aquatic taxa (Supplementary Fig. 6). Shape variation captured by the first four PC axes was similar for both the original and size-corrected shape data (Supplementary Fig. 7), with the main difference being increased shift in the lateral extension of the maxilla and otic regions on PC1 for the size-corrected data.

**Ancestral shape and trait estimation.** Ancestral state estimation of microhabitat and larval feeding mode supported the ancestral condition for frogs as terrestrial and with actively feeding larvae (Fig. 3). Transitions between all microhabitats are widespread across frogs, and a non-feeding larval stage has arisen at least 15 times among the taxa included in our analysis. A terrestrial ancestor with feeding larvae is consistent with the ancestral reconstruction of anuran cranial shape (Fig. 3), which we estimate as most similar to *Minervarya nilagirica* (which likely is terrestrial with an actively feeding larval stage).

**Evolutionary rate and disparity.** Cranial regions show considerable heterogeneity in evolutionary rate (σ²ₘₜₐₗ) and disparity (Procrustes variance). The highest evolutionary rates and disparity were concentrated in the neopalatine and quadratojugal, reflecting extensive variation in feeding mechanisms that impact

the relative position, orientation and extent of ossification of the jaw joint (Table 1 and Fig. 4). The frontoparietal, occipital and nasal evolved slowest and were least disparate, likely reflecting the conserved functions of these regions (e.g. protection of the brain and sensory organs and neck articulation).

When taxa were divided based on adult microhabitat, striking differences in evolutionary rate and disparity were observed for the whole cranium and most cranial regions (Table 1, Fig. 5 and Supplementary Fig. 8). The crania of fossorial, semi-fossorial and aquatic species evolved at the fastest rates (Fig. 5), with semi-aquatic, arboreal and terrestrial species exhibiting the slowest rates. Across individual cranial regions, fossorial, semi-fossorial and aquatic species also exhibited considerably, and, in many comparisons, significantly, faster rates of evolution than the remaining microhabitat specialists, especially for the quadratojugal and squamosal regions (and the stapes for fossorial and semi-fossorial species) (Fig. 5). The whole cranium and most cranial regions of fossorial, semi-fossorial and aquatic species were also significantly more disparate than those of taxa from the remaining microhabitats (Supplementary Fig. 8 and Supplementary Table 3).

No significant difference was found between the rate of cranial evolution for taxa with or without feeding larvae across the whole skull, but taxa with non-feeding larvae showed significantly faster evolutionary rates in eight of the 15 analysed regions: the nasal, otic, neopalatine, stapes, dorsal sphenethmoid, vomer and occipital (Supplementary Fig. 9). Rates of evolution for taxa with non-feeding larvae were also higher for most other cranial regions, although not significantly so (Supplementary Fig. 9 and Supplementary Table 4). Despite representing a minority of species in our dataset (39 non-feeding vs 124 feeding), taxa with non-feeding larvae had higher disparities for the whole cranium (non-feeding: 0.0120, feeding: 0.0095 *p* = 0.135) and most individual regions, although this difference was only significant for four cranial regions (premaxilla, otic, sphenethmoid (ventral) and vomer) and not for the entire cranium (Supplementary Fig. 9 and Supplementary Table 4).

Ossification sequence rank (relative timing) was significantly correlated with disparity (*ρ* = 0.76, *p* = 0.002) and evolutionary
<figure><img src="image_6.png" /><figcaption>**Fig. 3 Ancestral state and shape reconstructions.** A Ancestral state estimation for microhabitat (Total *N* for microhabitat = 170, three specimens with missing microhabitat data, as described in Methods). B Ancestral state estimation for larval feeding (Total *N* for larval feeding = 163, ten specimens with missing larval feeding mode data, as described in Methods). C Ancestral shape reconstruction for cranial shape in dorsal (left) and lateral (right) aspect. D Specimen cranium closest in shape to the ancestral shape reconstruction (Minervarya nilagirica, MNHN:RA:19842337) in dorsal (left) and lateral (right) aspect.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>**Fig. 4 Disparity across the frog cranium.** Landmarks and semilandmarks shown on Adenomus kelaartii (FMNH: Amphibians and Reptiles: 1580) colour graded by magnitude of disparity (Procrustes variance), from low (purple) to high (red), in **a** dorsal, **b** anterior, **c** ventral, **d** posterior and **e** lateral aspect.</figcaption></figure>

rate (mean $\rho$ = 0.50, $p < 0.01$ with error estimated, mean $\rho$ = 0.61, $p < 0.01$ without error), with later-ossifying bones exhibiting greater disparity and faster rates of evolution (Fig. 6 and Supplementary Figs. 10, 11).

**Discussion**

Our analyses support reconstructions of the ancestral frog as terrestrial and with an actively feeding larval stage⁶⁵, which was subsequently lost at least 15 times. Given our high-level sampling
<figure><img src="image_8.png" /><figcaption>**Fig. 5** Region-specific rates of evolution for different microhabitats. Evolutionary rate for each of the 15 cranial regions and the cranium across different microhabitats estimated across a random sample of 100 trees drawn in the posterior distribution of Jetz and Pyron¹¹³. Intraspecific variation and measurement error were jointly estimated along with evolutionary rates during model fit to account for uncertainty in traits values that may cause model departure. Rates are standardised by the mean, and box bounds represent the median with the first and third quartile, while whiskers indicate the minimum and maximum values within 1.5x the interquartile range.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**Fig. 6** Ossification sequence and evolutionary rate. Ossification sequence rank versus evolutionary rate for each cranial region. Evolutionary rates are estimated with error to account for intraspecific variance and departures from Brownian Motion (results estimated without error in Supplementary Fig. 11), with lines derived from a robust regression fit. Two-sided Spearman's rank correlation of evolutionary rate with ossification sequence rank (from Weisbecker & Mitgutsch⁴³) was significant (mean Spearman's rho = 0.50, *p* < 0.01). Colours for cranial regions are as indicated in Fig. 1.</figcaption></figure>

and focus on one aspect of life history, transitions to derived life history modes, including direct development, are likely even more complex and more widespread in anuran evolution. Transitions in adult microhabitat are also common across frogs, providing an ideal study system for comparing the relative influences of potential selective pressures in larval and adult stages. Our results support microhabitat and skull size as the dominant factors driving cranial shape evolution in frogs, with aquatic, fossorial and semi-fossorial taxa displaying the fastest rates of evolution (Fig. 5), the highest cranial morphological diversity (Supplementary Fig. 8) and distinct cranial morphologies (Fig. 2). Larval feeding mode is also a significant influence on skull shape across the whole skull and most cranial regions, with loss of feeding larvae consistently associated with greater disparity and faster evolutionary rates than those that maintain an actively feeding larval stage.

Adaptation to adult microhabitats has resulted in substantial variation in cranial shape, disparity and rate of evolution across frogs, particularly for aquatic, fossorial and semi-fossorial taxa. Subterranean habitats can minimise exposure to environmental and biotic fluctuations, but are also functionally challenging, creating considerable locomotory and sensory system constraints⁶⁶,⁶⁷. Similarly, the dense aquatic medium imparts significant functional pressures on locomotion and feeding⁶⁸. These challenging habitats can promote¹¹,²²,⁶² or limit⁶⁹ phenotypic evolution, and here we find evidence of the former, likely due to divergent specialised feeding modes within species considered the same microhabitat specialist (e.g. *Pyxicephalus* and *Nasikabatrachus* are both fossorial frogs). Interestingly, the crania of semi-aquatic taxa are the least disparate and slowest evolving, suggesting their more generalist approach does not impose the same functional pressures. Alternatively, semi-aquatic taxa may be limited by the divergent selection pressures of both aquatic and terrestrial microhabitats⁷⁰, an interpretation supported by our
findings that semi-aquatic and terrestrial frogs overlap in mor-
phospace and display similar patterns of cranial disparity and
evolutionary rate. For the shape of individual elements, those that
are most strongly associated with adult microhabitat also show
the largest differences in evolutionary rates across microhabitats.
Recent work has also demonstrated that adult microhabitat
interacts significantly with presence of hyperossification of the
skull, which is found in all microhabitat specialists, further sup-
porting adult microhabitat as a key factor influencing many
aspects of cranial shape in anurans⁶⁴.

Size is also a major influence on morphology in many clades (e.g.³⁵), and our analyses support allometry as a significant component of anuran cranial variation, with cranial width increasing with size. Our results for anurans, in combination with previous studies of caecilians and salamanders²²,⁶⁴,⁷¹, suggest an amphibian counterpart to the mammalian-based craniofacial evolutionary allometry rule (CREA) that large body size results in long faces³⁵,⁷². Instead, however, large size results in wide skulls²²,⁶⁴,⁷¹ in amphibian cranial evolutionary allometry (ACREA). This pattern may be driven by both structural constraints and feeding adaptations. Structurally, we find that larger crania have proportionally smaller braincases (congruent with mammals⁷²), and the occipital region shows the strongest relationship to size, as found across caecilians²². Smaller crania exhibit relatively larger occipital regions, otic regions and cranial vaults, and more vertically-oriented suspensory apparatuses in order to maintain the minimum space requirements for the brain and sensory organs⁷³. Thus, size constraints related to the nervous and sensory systems have a greater influence on the morphology of smaller crania. In terms of feeding, frogs are gape-limited predators⁷⁴, and the consumption of large prey requires wide skulls⁷⁵. Functional requirements related to feeding may therefore contribute to the ‘large size-wide skull’ pattern of skull allometry. For example, suction feeding, which is present in some aquatic species, requires a wider mouth⁷⁶⁻⁷⁸, supporting a link between cranial size, feeding function and microhabitat. Relatedly, we recovered a highly significant interaction between microhabitat and cranial size across the skull and in most cranial regions, with aquatic taxa exhibiting larger crania.

Feeding mechanics appears to be a considerable driving force behind the cranial evolution of frogs. The jaw articulation region of frogs (quadratojugal) evolves at the fastest rate, mirroring previous findings<sup>22,62</sup>, and it is the most morphologically diverse cranial region. Previous study has demonstrated that the suspensorium of anurans, comprising the quadratojugal, jaw joint, squamosal and pterygoid, is highly integrated, likely facilitating functional evolution of this region<sup>79</sup> (and not only in anurans<sup>80,81</sup>). While evolutionary rate and disparity vary widely within and across cranial regions, both are high in the quadratojugal and adjacent surfaces (Fig. 4). Moreover, differences in evolutionary rates across microhabitats are pronounced for the jaw suspensorium cranial regions, further supporting that different feeding modes across microhabitats drive divergences in evolutionary tempo in the anuran skull (Fig. 5). There is extensive convergence in feeding mechanisms across Anura<sup>82</sup>, with feeding specialisations (hydrostatic elongation and suction feeding) occurring in some microhabitat specialists such as fossorial and aquatic taxa (although information on feeding mode and diet remains poor across frogs). Suction feeding, restricted to aquatic species, requires wide skulls<sup>76,77</sup>, resulting in a more lateralised jaw suspensorium. In contrast, hydrostatic elongation, which is observed in some fossorial frogs<sup>83,84</sup> (and their close non-fossorial relatives, especially in the Microhylidae), requires minimal jaw opening, where only the lingual tip needs to leave the mouth<sup>83</sup>. This, and the small prey size (i.e. ants, termites) of some fossorial species<sup>83-86</sup>, allows the jaw suspensorium and the

pectoral girdle to move anteriorly, which effectively eliminates the
neck region⁸³,⁸⁵. Other fossorial species have wide, heavily ossi-
fied crania for consuming large prey including vertebrates (e.g.
*Pyxicephalus adspersus*), and a recent study also identified higher
rates of evolution, but not higher disparity, in hyperossified frogs,
with a significant interaction with microhabitat⁶⁴. Thus, the ele-
vated rate of evolution and disparity across the crania of aquatic,
fossorial and semi-fossorial frogs may be linked to the multiple
feeding modes found within each of these microhabitats.

Hearing ability also may vary across microhabitats. The stapes displays the highest disparity in fossorial, semi-fossorial and aquatic taxa and show significant associations with adult micro-habitat and size, but not larval feeding. Differences in anuran stapes evolution across microhabitats and its relatively high disparity may also be driven by the convergent losses, reappearances and re-losses of this structure⁵⁵, which has also complicated its inclusion in most previous morphometric analyses. The stapes may be lost in frogs with little loss of acoustic communication⁸⁷,⁸⁸, and middle ear development can be heavily protracted⁸⁹, suggesting the stapes may be 'evolutionarily dispensable'. The other tympanic middle ear structures (tympanic membrane and tympanic annulus) are also repeatedly lost throughout anuran evolution⁵⁵. Variation in extratympanic hearing sensitivity⁸⁷,⁹⁰ suggests species-specific extratympanic hearing mechanisms may exist⁹⁰, for example sound detection through the mouth⁸⁸. The different acoustic potentials of substrates⁹¹ may drive the evolution of alternative, extratympanic hearing mechanisms and promote lability in stapes morphological evolution across microhabitat specialists, meriting further investigation.

In addition to skull size and microhabitat, larval feeding mode
is a significant albeit relatively weaker influence on anuran cranial
evolution. The adult crania of frogs with and without a feeding
larval stage do not differ notably in morphology (Supplementary
Fig. 3). Of the factors studied here, there are the fewest significant
associations between larval feeding and morphology of cranial
regions. Previous study also demonstrates that taxa with or
without actively feeding larvae do not differ in patterns of evo-
lutionary integration for the cranium⁷⁹, in contrast to the
markedly different patterns observed among salamanders with
different reproductive strategies⁸¹. Our results thus indicate that
functional adaptations to the larval environment are less influ-
ential for adult morphology than are adaptations for the adult
habitat and support previous suggestions for the independence of
the larval and adult stages. These results further support previous
studies identifying similarities in adult cranial morphology
despite differences in early development in frogs and in marine
invertebrates³⁷,⁹²⁻⁹⁴. However, we found that taxa with non-
feeding larvae consistently showed faster rates of evolution, sig-
nificantly so for the majority of cranial regions, as well as higher
disparity, suggesting that loss of a feeding larval stage still pro-
motes morphological diversification in anurans.

As noted above, previous studies directly comparing larval and
adult traits have variably found them to be correlated<sup>45,95-97</sup> or
independent<sup>45,47,98,99</sup>. These conflicting results among studies
may reflect different effects across cranial elements, which were
not treated individually in previous studies and which our results
demonstrate display differential relationships with adult micro-
habitat and larval feeding. We found highly significant interac-
tions between adult microhabitat and larval feeding for
approximately half of the cranial elements (Supplementary
Table 2). Elements with a weaker or no significant interaction
between adult microhabitat and larval feeding, including the
nasal, neopalatine, dorsal sphenethmoid and vomer, show large
differences in evolutionary rate between taxa with or without a
feeding larval stage. This difference was observed even when there
is not a significant association between region shape and larval
feeding. In contrast, those same elements show little difference in evolutionary rate across microhabitat categories, even if there is a significant association between shape and microhabitat. These results suggest two unexpected patterns. First, differences in shape and evolutionary rates are decoupled in many elements of the anuran cranium. This mosaicism in both the drivers of cranial shape and pace of evolution supports recent work demonstrating extensive modularity in the frog skull⁷⁹. Second, there is a trade-off in the influence of adult microhabitat and larval feeding on evolutionary rates that varies across cranial elements. More specifically, differences in evolutionary rates among taxa with or without feeding larvae are most pronounced, and differences in evolutionary rates among microhabitats least pronounced, when there is not a significant interaction between adult microhabitat and larval feeding. In these cases, loss of an actively feeding larval stage consistently spurs significantly faster evolution of individual cranial regions.

It may be expected that frogs with non-feeding larvae have relaxed their constraints on certain aspects of larval cranial morphology, given the lack of selection pressures related to feeding at earlier stages of development. This may manifest as higher evolutionary rates, rather than a straightforward association with specific morphologies. Our finding that the mean evolutionary rate is higher for taxa with non-feeding larvae for nearly all cranial regions, and significantly so in most, provides strong support for the hypothesis of released constraints on cranial morphology. Direct development also has been hypothesised to promote morphological novelty across plethodontid salamanders¹⁹,²⁰, suggesting that presence of a larval stage may constrain adult morphology. However, we would expect to observe this effect most strongly in the jaw suspensorium region, as frogs with feeding larvae exhibiting dramatic ontogenetic changes in diet, transitioning generally from small-mouthed microphagous suspension feeders to large-mouthed carnivorous adults³³,³⁸⁻⁴⁰. The primary cranial shape change through metamorphosis results from a posterior translation of the quadrate bar (which supports the jaw joint)⁴⁰,¹⁰⁰, reflecting the change in feeding mode¹⁰⁰. In direct developers (all of which have non-feeding larvae), the quadrate bar varies widely, forming in positions that correspond to larval, mid-metamorphic or adult stages⁴⁰. A posterior translation of the quadratojugal (jaw joint) also contributes to the major axis of shape variation across adult frog crania (Fig. 2 and Supplementary Fig. 1), and the quadratojugal is the fastest-evolving cranial region in anurans. Yet we did not find a significant difference in evolutionary rate based on larval feeding for this structure or the adjacent and functionally related maxilla and squamosal. Morphological repatterning through ontogeny related to larval feeding mode may play a major role in adult morphological evolution. However, the wide range of ontogenetic trajectories among frogs and the retention of ancestral features in some taxa with non-feeding larvae⁹³ complicates identification of a more subtle relationship between larval feeding and adult morphological evolution. Quantifying ontogenetic trajectories of cranial shape across the breadth of frog diversity—and even across different larval feeding modes and microhabitats—is a challenging task for future studies but would bring critical data towards refining our understanding of this relationship.

Overall, we interpret our results to indicate that adult micro-habitat and body size shape large-scale differences in morphology for most structures, but both adult microhabitat and loss of actively feeding larvae promote diversification of cranial form. The relationship between larval feeding and adult morphology in anurans is likely to be more complex than that of size and microhabitat, with effects manifested within specific systems, whether clades, anatomical elements or microhabitats, but not at

a broader phylogenetic scale. This effect may extend to specific morphologies, as well as evolutionary rate and disparity. For example, larval feeding has been suggested to influence adult feeding mode in aquatic amphibians, specifically with suction feeding adults arising predominantly from suction feeding larvae⁹⁷. Our results similarly suggest that these effects of larval feeding mode are likely to manifest within specific structures and specific habitats rather than across the entire skull. Whether specific larval feeding modes and larval diets might differentially impact the disparity and rates of evolution of cranial structures remains unknown.

Another aspect of development that is more readily quantifiable is ossification sequence. Timing of bone formation varies across taxa, with, for example, earlier formation of jaw suspensorium bones in direct developers, related to their lack of larval feeding and earlier requirement of adult-mode feeding²⁸,¹⁰¹,¹⁰². Thus, ossification sequence rank may reflect functional requirements that differentially constrain the morphological evolution of some cranial elements, as has been demonstrated in marsupial mammals¹⁸. We find strong evidence that the early onset of ossification constrains the evolution of individual skull elements, and that later-ossifying bones are more evolutionarily labile, with a strong positive relationship between median ossification sequence rank (relative timing) and both evolutionary rate and disparity (Fig. 6 and Supplementary Figs. 10, 11). The early-ossifying parasphenoid, frontoparietal and occipital regions are involved in the protection of the brain and otic capsules³³,¹⁰³ and display the least restructuring during metamorphosis²⁸, suggesting they experience strong functional constraints. In contrast, many later-ossifying bones, such as the vomer, are variably present across Anura (possibly a result of progenesis⁴³ or miniaturisation⁷³) and undergo greater morphological changes during metamorphosis²⁸, suggesting their functions are less conserved and their morphologies thus less constrained. Previous studies have demonstrated that ossification in anurans occurs first in areas experiencing high stress¹⁰⁴. Ossification sequence rank is thus possibly linked to functional demands¹⁰⁵, with the later-ossifying bones being less functionally critical²⁸ and thus less constrained in morphology. Interestingly, recent study has suggested that later-ossifying elements are also more evolutionarily integrated⁷⁹, contrary to expectations that more functionally conservative regions may be more integrated. It may be that extensive variability in the late-ossifying elements, particularly with several elements that are variably present across anurans, requires tightly coordinated changes that drive high evolutionary integration within those elements, as well as greater modularity among elements. Overall, the results presented here support the long-standing hypothesis that early developmental stages are more evolutionarily conservative than late developmental stages¹⁰⁵,¹⁰⁶ and suggests that late ossification of cranial elements, and perhaps strong integration of those elements, may allow or even promote variability in form and function.

The coupling of dense phenotypic data with dense taxonomic sampling presented here has provided a framework with which a detailed, quantitative investigation into the drivers of phenotypic evolution across the anuran skull is achievable. Here we find microhabitat and size are the primary drivers of cranial morphological evolution across frogs, with aquatic and fossorial microhabitats associated with distinct cranial forms, faster evolutionary rates and greater morphological diversity. Evolutionary rate and disparity are particularly elevated in the jaw suspensorium regions, indicating that feeding mechanics related to micro-habitat is likely a significant driver of anuran cranial evolution. Feeding mechanics also contributes to a 'large size-wide skull' pattern of allometry across anurans, with large prey specialists requiring wider skulls. In combination with similar results for
caecilians²²,⁷¹, we suggest that this pattern may form a general rule of ACREA. Surprisingly, taxa with and without feeding larvae do not notably differ in their cranial morphology. Although larval feeding mode is significantly associated with shape of most cranial elements, it has a much smaller effect size than microhabitat or size on cranial morphology. Significant differences in evolutionary rate and disparity are observed in half of the cranial elements examined between taxa with and without actively feeding larvae. Where differences do occur, they consistently support higher rates of evolution in taxa that have lost an actively feeding larval stage. Our results thus support both that larval feeding constrains at least the rate of morphological evolution and that amphibian larvae serve as a distinct developmental module²³. Even if only a weak influence, functional pressures experienced earlier in ontogeny appear to still impact patterns of adult morphological evolution, and our results demonstrate a significant interaction between adult microhabitat and larval feeding across the skull and in most cranial regions. An interesting area for future research will be to conduct finer-scale comparisons of morphological evolution among anuran clades that have diversified into similar sets of microhabitat specialists but differ in their life histories. For example, do lineages with non-feeding larvae evolve more quickly to fill this morphospace than those with typical feeding tadpoles? Future research directly incorporating quantitative shape data through ontogeny and making explicit comparisons between clades with contrasting life histories or larval ecologies will be instrumental for further clarifying the complex interactions between adult and larval ecology and morphological evolution across >200 million years of frog evolution.

## Methods

**Specimens.** We reconstructed and processed meshes from microCT scans of crania for 173 anuran species, including representatives from all extant frog families and the mummified Eocene *Thaumastosaurus gezi*¹⁰⁷ (Fig. 2 and Supplementary Table 5). A recent phylogenetic analysis suggests *Thaumastosaurus gezi* is a ranoid with African affinities¹⁰⁷, nested within the Natatanura clade from Frost et al.¹⁰⁸. Although complete, undeformed 3D fossil material for anurans is limited, inclusion of any fossils increases our temporal sampling and has been demonstrated to improve estimates of evolutionary rate¹⁰⁹,¹¹⁰. Meshes from microCT scans were created in Avizo Lite 9 (FEI Visualisation Sciences Group, Burlington, MA, USA) and VG Studio MAX¹¹¹ and processed in Geomagic Wrap 2017 (3D Systems, Rock Hill, South Carolina, USA) to remove noise and small surface foramina, as previously described¹¹². Scan information can be found in Supplementary Data 1.

**Phylogeny.** We reconstructed the maximum clade credibility tree (MCC) from the posterior sample of the most recent, comprehensive, time-calibrated phylogeny of Anura¹¹³ using treeAnnotator. We pruned the MCC tree to match our species using the 'drop tip' function in the R package ape v.5.3¹¹⁴,¹¹⁵ and modified it by adding four specimens without species assignments to their respective genera within the phylogeny. Specifically, *Raorchestes* sp., *Dendrobates* sp., *Capensibubo* sp. and *Xenohina* sp. took the positions of *Raorchestes anili*, *Dendrobates auratus*, *Capensibubo rosei* and *Xenohina varia*, respectively. In addition, the fossil species *Thaumastosaurus gezi* was added as sister taxon to the clade comprising *Ceratobatrachus*, *Aubria* and *Pyxiscephalus*, as recovered from a recent phylogenetic analysis¹⁰⁷. We further randomly drew 100 trees from the posterior distribution¹¹³ to assess the effect of phylogenetic uncertainties in downstream comparative analyses (the fossil taxon was excluded from the analyses of the 100 trees).

**Morphometric data collection.** Fifteen cranial regions were defined (Fig. 2 and Supplementary Table 6) using 58 landmarks, 410 curve semilandmarks and 527 surface semilandmarks (Supplementary Tables 7–9). Landmarks and curve-semilandmarks were digitised in iDAV Landmark Editor v.3.6¹¹⁶, while surface semilandmarks were applied to each cranium semi-automatically¹¹² using the R package Morpho v.2.7¹¹⁷. Regions that were variably present across taxa were represented as ‘negligible regions’ when absent²²,¹¹² (Supplementary Table 5). The effect of negligible regions on trait integration was investigated previously for the dataset in the current study⁷⁹, where the pattern of trait integration across the anuran skull was found to be near identical for datasets including and excluding specimens with negligible regions. Negligible regions represent the loss of that region and are therefore important to retain when accurately reconstructing evolutionary rate. Exclusion of negligible regions would underestimate evolutionary rate as the evolutionary loss of structures, an extreme form of variation and evolutionary change, would not be captured. Semilandmarks were slid to minimise bending energy globally,

then all data were mirrored to improve Procrustes alignment¹¹⁷. Only right-side and midline data were retained for analysis to reduce redundant dimensionality⁶⁰–⁶²,¹¹⁸.

**Size.** Log centroid size of the cranium was extracted from shape data during Procrustes superimposition and was used as a proxy of overall size (Supplementary Table 9).

**Habitat and larval feeding.** We used functionally relevant, discrete categories for microhabitat and larval feeding, which were necessarily limited by data availability for most of the taxa in our dataset. Microhabitat was divided into the following seven categories: aquatic, semi-aquatic, arboreal, semi-arboreal, fossorial, semi-fossorial and terrestrial, based on adults observed outside of the breeding season. The habitat categorisation reflects the four main habitats identified across frogs³⁰, with the recognition that intermediate forms between terrestrial and the other three habitats also exist³⁰. In addition, these categories are in line with previous studies using habitat categories³¹,³²,¹¹⁹–¹²¹. Ecological categories were based on predominant (or partial for ‘semi-’ categories) use of a substrate (water for aquatic, vegetation for arboreal and soil for fossorial) compared with land (terrestrial). Whilst fossorial species can dig head-first, forelimb first or hind limb first, and whilst some species found in burrows may not actively burrow, this information is not always available, and we do not have sufficient numbers of representatives from each category. We therefore grouped these into ‘fossorial’ (and ‘semi-fossorial’) categories, as this would capture general selection pressures associated with endogeicity (living in soil), as previously discussed for caecilians²². Furthermore, studies have recognised that whilst some frogs may burrow backwards, they move forward underground and excavate termite colonies by digging forwards through tunnels e.g. *Rhinophrynus dorsalis*⁸³, so would likely experience similar selection pressures underground to forward fossorial. For a detailed investigation into burrowing influences on morphological evolution for myobatrachid frogs, we refer to Vidal-García and Keogh³. Data from the primary literature were used, as well as categories established in previous studies³¹,³² (Supplementary Data 1). Three taxa did not have microhabitat information available and were excluded from analyses of microhabitat (Habitat: *N* = 170).

Larval feeding categorisation used the following definitions³⁸: taxa with feeding (exotrophic) larvae (i.e. feed on external food sources) and taxa with non-feeding (endotrophic) larvae (i.e. solely provisioned with yolk). The larval feeding categorisation distinguishes frogs experiencing extrinsic selection pressures from two environments (larval and adult), to those experiencing only pressures from the adult environment. Whilst direct developers can vary in their ontogeny, all are assumed to be endotrophic¹²². Trophic egg feeding or feeding on skin secretions were both counted as ‘feeding’ larval forms. We use this categorisation rather than metamorphic vs direct developing because metamorphic forms include those with and without feeding larvae and thus that categorisation does not capture the differences in cranial function and selection pressures imposed by a feeding larval stage. Data were taken from the primary literature (Supplementary Data 1). Ten taxa did not have information on larval feeding and were excluded from those analyses (larval feeding *N* = 163; 162 taxa have data on both larval feeding and microhabitat). While finer categories of both developmental and ecological traits are possible, these broader categorisations of adult microhabitat and larval feeding are necessary for robust statistical analyses, which require a minimum of five species per group. Other traits would suffer from much more missing data, while finer categories would result in many bins with only few species within them.

**Ossification timing.** Data on relative timing of cranial bone ossification were taken from Weisbecker and Mitgutsch (2010)⁴³ using median rank position for all bones except the stapes (which is not included in that study). The occipital region defined in this study corresponds to the oto-occipital bone, and the otic region corresponds to the prootic bone (Supplementary Table 5).

## Data analyses

**Cranial shape.** Shape analyses were conducted using the R package geomorph v3.1.3¹²³ unless otherwise noted. PCs analysis was used to assess the main axes of shape variation for the whole cranium and individually for the 15 cranial regions. Phylomorphospaces were plotted in phytools v0.6-60¹²⁴, and the primary axes of shape variation were visualised by generating extreme shapes along the main PC axes. Size-related shape variation (allometry) was extracted for the whole cranium using the ‘procD.lm’ function, and we repeated the morphospace analyses on the residuals of shape variables regressed against log centroid size to visualise the distribution of microhabitat and larval feeding categories independent of allometry. Morphologies representing the extreme shapes along each PC axis that represented over 5% of the variation were visualised for the entire cranium, and along PC1 for each cranial region. For region-specific analyses of shape, only specimens with the relevant region present were retained. It is important to note that the directionality of the PC axes was arbitrary (Supplementary Fig. 2), as the PCs were conducted on individual cranial regions.

Type II phylogenetic MANOVAs (phylogenetic regressions) were performed on the landmark and semilandmark data for the whole cranium and individual cranial regions with log centroid size, microhabitat and larval feeding as predictors using the ‘mvglsl’ and ‘manova gls’ functions in the R package mvMORPH 1.1.4¹²⁵. Specifically, we fitted multivariate phylogenetic linear models with Pagel’s lambda
by penalised likelihood using 'mvgl's¹²⁶ and assessed the significance of each
predictors in this model by type II MANOVA tests using the 'manova.gls' function
with Pillai's statistic and 1000 permutations¹²⁷. Using Pagel's lambda corresponds
to fitting a phylogenetic mixed model which allows accounting for departure from
Brownian motion and usually provides increased flexibility in estimating the error
structure¹²⁷⁻¹³⁰. MANOVAs were performed on both the MCC tree and from the
100 trees randomly sampled in the posterior distribution¹¹³ to assess the effect of
phylogenetic uncertainties.

*Ancestral shape and trait estimation.* Ancestral cranial morphology for anurans was estimated by maximum likelihood using the ‘anc.recon’ function in the R package Rphylopars v0.2.11¹³¹⁻¹³³, and the ancestral states for microhabitat and larval feeding were estimated using the ‘ace’ function in *ape* v5.3.

*Disparity and evolutionary rate.* Evolutionary rate and morphological disparity (measured as Procrustes variance scaled by number of landmark/semilandmarks) were calculated for each cranial region across the entire dataset and for each microhabitat and larval feeding category after removing the effect of size (log cranial centroid size), using the state-specific Brownian motion (BMM) model in the ‘mvgl’s function in mvMORPH and ‘morphol.disparity’ functions in geomorph. The reconstructed history for microhabitat and larval feeding categories on which BMM models were fitted was obtained through stochastic character mapping across the sample of 100 trees using the ‘make.simmap’ function and an ‘ARD’ model in the R package phytools. Model fit in ‘mvgl’s was performed by jointly estimating the contribution of measurement error and intraspecific variation using the option ‘error=TRUE’ to mitigate the potential sources of BM departure that may lead to evolutionary rates inflation. Overall morphological disparity and evolutionary rate—taken as the average of rates estimated across coordinates and calculated across a sample of 100 trees—for each region were compared to ossification sequence rank (except for the stapes, which lacks data on ossification timing) using Spearman’s rank correlation analysis. Disparity and evolutionary rate were also calculated for each individual landmark and semilandmark in geomorph to visualise their distributions across the anuran skull on a finer scale.

**Reporting summary.** Further information on research design is available in the Nature Research Reporting Summary linked to this article.

**Data availability**

Meshes are available on the online repositories Morphosource.org and phenome10k.org for use by other researchers. Morphosource DOIs are detailed in Supplementary Data 1. Landmark and semilandmark data are available at https://github.com/anjigoswami/frogs_modularity¹³⁴. All other trait data and sources (e.g. ecological and developmental traits) are provided in Supplementary Data 1.

Received: 21 April 2020; Accepted: 25 March 2021;
Published online: 04 May 2021

**References**

1. Collar, D. C., Schulte, J. A., O'Meara, B. C. & Losos, J. B. Habitat use affects morphological diversification in dragon lizards. *J. Evol. Biol.* **23**, 1033–1049 (2010).

2. Da Silva, F. O. et al. The ecological origins of snakes as revealed by skull evolution. *Nat. Commun.* **9**, 1–11 (2018).

3. Vidal-García, M. & Keogh, J. S. Phylogenetic conservatism in skulls and evolutionary lability in limbs – morphological evolution across an ancient frog radiation is shaped by diet, locomotion and burrowing. *BMC Evol. Biol.* **17**, 1–15 (2017).

4. Fabre, A.-C., Cornette, R., Goswami, A. & Peigné, S. Do constraints associated with the locomotor habitat drive the evolution of forelimb shape? A case study in musteloid carnivorans. *J. Anat.* **226**, 596–610 (2015).

5. Dumont, M. et al. Do functional demands associated with locomotor habitat, diet, and activity pattern drive skull shape evolution in musteloid carnivorans? *Biol. J. Linn. Soc.* **117**, 858–878 (2015).

6. Baeckens, S., Goeyers, C. & Van Damme, R. Convergent evolution of claw shape in a transcontinental lizard radiation. *Integr. Comp. Biol.* https://doi.org/10.1093/icb/icz151 (2019).

7. Price, S. A., Holzman, R., Near, T. J. & Wainwright, P. C. Coral reefs promote the evolution of morphological diversity and ecological novelty in labrid fishes. *Ecol. Lett.* **14**, 462–469 (2011).

8. Price, S. A., Tavera, J. J., Near, T. J. & Wainwright, P. C. Elevated rates of morphological and functional diversification in reef-dwelling haemulid fishes. *Evolution* **67**, 417–428 (2012).

9. Millien, V. Morphological evolution is accelerated among island mammals. *PLoS Biol.* **4**, 1863–1868 (2006).

10. Salvidio, S., Crovetto, F. & Adams, D. C. Potential rapid evolution of foot morphology in Italian plethodontid salamanders (*Hydromantes strinatii*) following the colonization of an artificial cave. *J. Evol. Biol.* **28**, 1403–1409 (2015).

11. Ledbetter, N. M. & Bonett, R. M. Terrestriality constrains salamander limb diversification: implications for the evolution of pentadactyly. *J. Evol. Biol.* **32**, 642–652 (2019).

12. McGhee Jr, G. R. *Convergent Evolution: Limited Forms Most Beautiful* (MIT Press, 2011).

13. Vullo, R., Allain, R. & Cavin, L. Convergent evolution of jaws between spinosaurid dinosaurs and pike conger eels. *Acta Palaeontol. Pol.* **61**, 825–828 (2016).

14. Stayton, C. T. Testing hypotheses of convergence with multivariate data: morphological and functional convergence among herbivorous lizards. *Evolution* **60**, 824–841 (2006).

15. Mahler, D. L., Ingram, T., Revell, L. J. & Losos, J. B. Exceptional convergence on the macroevolutionary landscape in island lizard radiations. *Science* **341**, 292–5 (2013).

16. Sears, K. E. Constraints on the morphological evolution of marsupial shoulder girdles. *Evolution* **58**, 2353–2370 (2004).

17. Bennett, C. V. & Goswami, A. Statistical support for the hypothesis of developmental constraint in marsupial skull evolution. *BMC Biol.* **11**, 1–14 (2013).

18. Goswami, A. et al. Do developmental constraints and high integration limit the evolution of the marsupial oral apparatus? *Integr. Comp. Biol.* **56**, 404–415 (2016).

19. Wake, D. B. & Hanken, J. Direct development in the lungless salamanders: what are the consequences for developmental biology, evolution and phylogenesis? *Int. J. Dev. Biol.* **40**, 859–869 (1996).

20. Wake, D. B. & Larson, A. Multidimensional analysis of an evolving lineage. *Science* **238**, 42–48 (1987).

21. Bonett, R. M. & Blair, A. L. Evidence for complex life cycle constraints on salamander body form diversification. *Proc. Natl Acad. Sci. USA* **114**, 9936–9941 (2017).

22. Bardua, C., Wilkinson, M., Gower, D. J., Sherratt, E. & Goswami, A. Morphological evolution and modularity of the caecilian skull. *BMC Evol. Biol.* **19**, 1–23 (2019).

23. Schlosser, G. in *Modularity: Understanding the Development and Evolution of Natural Complex Systems* (eds. Callebaut, W. & Rasskin-Gutman, D.) (MIT Press, 2005).

24. Moran, N. A. Adaptation and constraint in the complex life cycles of animals. *Annu. Rev. Ecol. Syst.* **25**, 573–600 (1994).

25. Ebenman, B. Evolution in organisms that change their niches during the life cycle. *Am. Nat.* **139**, 990–1021 (1992).

26. Mallarino, R. et al. Two developmental modules establish 3D beak-shape variation in Darwin’s finches. *Proc. Natl Acad. Sci. USA* **108**, 4057–4062 (2011).

27. Liedtke, H. C. et al. Terrestrial reproduction as an adaptation to steep terrain in African toads. *Proc. R. Soc. B Biol. Sci.* **284**, 20162598 (2017).

28. Harrington, S. M., Harrison, L. B. & Sheil, C. A. Ossification sequence heterochrony among amphibians. *Evol. Dev.* **15**, 344–364 (2013).

29. Bonett, R. M., Phillips, J. G., Ledbetter, N. M., Martin, S. D. & Lehman, L. Rapid phenotypic evolution following shifts in life cycle complexity. *Proc. R. Soc. B Biol. Sci.* **285**, 20172304 (2018).

30. Laurent, R. F. Adaptive modifications in frogs of an isolated highland fauna in Central Africa. *Evolution* **18**, 458–467 (1964).

31. Moen, D. S., Morlon, H. & Wiens, J. J. Testing convergence versus history: convergence dominates phenotypic evolution for over 150 million years in frogs. *Syst. Biol.* **65**, 146–160 (2016).

32. Moen, D. S., Irschick, D. J. & Wiens, J. J. Evolutionary conservatism and convergence both lead to striking similarity in ecology, morphology and performance across continents in frogs. *Proc. R. Soc. B Biol. Sci.* **280**, 1–9 (2013).

33. Duellman, W. E. & Trueb, L. *Biology of the Amphibians* (McGraw-Hill publishing company, 1986).

34. LaBarbera, M. in *Patterns and Processes in the History of Life* (eds. Raup, D.M. & Jablonski, D.) (Springer, 1986).

35. Cardini, A. & Polly, P. D. Larger mammals have longer faces because of size-related constraints on skull form. *Nat. Commun.* **4**, 2458 (2013).

36. Callery, E. M. & Elinson, R. P. Thyroid hormone-dependent metamorphosis in a direct developing frog. *Proc. Natl Acad. Sci. USA* **97**, 2615–2620 (2000).

37. Ziermann, J. M. & Diogo, R. Cranial muscle development in frogs with different developmental modes: direct development versus biphasic development. *J. Morphol.* **275**, 398–413 (2013).

38. McDiarmid, R. W. & Altig, R. (eds) *Tadpoles: The Biology of Anuran Larvae* (University of Chicago Press, 1999).
39. Altig, R. & Johnston, G. F. Guilds of anuran larvae: relationships among developmental modes, morphologies, and habitats. *Herpetol. Monogr.* **3**, 81–109 (1989).

40. Rose, C. S. & Reiss, J. O. in *The Skull Volume 1: Development* (eds. Hanken, J. & Hall, B. K.) (The University of Chicago Press, 1993).

41. Callery, E. M., Fang, H. & Elinson, R. P. Frogs without polliwogs: evolution of anuran direct development. *BioEssays* **23**, 233–241 (2001).

42. Wake, D. B. & Roth, G. (eds). Complex Organismal Functions: Integration and Evolution in Vertebrates (Wiley, Chichester, UK, 1989).

43. Weisbecker, V. & Mitgutsch, C. A large-scale survey of heterochrony in anuran cranial ossification patterns. *J. Zool. Syst. Evol. Res.* **48**, 332–347 (2010).

44. Dehling, J. M. & Sinsch, U. Partitioning of morphospace in larval and adult reed frogs (Anura: Hyperoliidae: Hyperolius) of the Central African Albertine Rift. *Zool. Anz.* **280**, 65–77 (2019).

45. Phung, T. X., Nascimento, J. C. S., Novarro, A. J. & Wiens, J. J. Correlated and decoupled evolution of adult and larval body size in frogs: larval and adult size evolution. *Proc. R. Soc. B Biol. Sci.* **287**, 20201474 (2020).

46. Werner, E. E. Amphibian metamorphosis: growth rate, predation risk, and the optimal size at transformation. *Am. Nat.* **128**, 319–341 (1986).

47. Sherratt, E., Vidal-García, M., Anstis, M. & Keogh, J. S. Adult frogs and tadpoles have different macroevolutionary patterns across the Australian continent. *Nat. Ecol. Evol.* **1**, 1385–1391 (2017).

48. Wollenberg Valero, K. C. et al. Transcriptomic and macroevolutionary evidence for phenotypic uncoupling between frog life history phases. *Nat. Commun.* **8**, 15213 (2017).

49. Trueb, L. in *The Skull: Patterns of Structural and Systematic Diversity* (eds Hanken, J. & Hall, B. K.) (The University of Chicago Press, 1993).

50. Trueb, L. in *Evolutionary Biology of the Anurans: Contemporary Research on Major Problems* (ed. Vial, J. L.) (University of Missouri Press, 1973).

51. Reiss, J. O. The phylogeny of amphibian metamorphosis. *Zoology* **105**, 85–96 (2002).

52. Moore, M. K. & Townsend, V. R. Jr Intraspecific variation in cranial ossification in the tailed frog, *Ascaphus truei*. *J. Herpetol.* **37**, 714–717 (2003).

53. Yeh, J. The evolution of development: two portraits of skull ossification in pipoid frogs. *Evolution* **56**, 2484–2498 (2002).

54. Schoch, R. R. Amphibian skull evolution: the developmental and functional context of simplification, bone loss and heterotopy. *J. Exp. Zool. B Mol. Dev. Evol.* **322B**, 619–630 (2014).

55. Pereyra, M. O. et al. The complex evolutionary history of the tympanic middle ear in frogs and toads (Anura). *Sci. Rep.* **6**, 1–9 (2016).

56. Long, J. A., Young, G. C., Holland, T., Senden, T. J. & Fitzgerald, E. M. G. An exceptional Devonian fish from Australia sheds light on tetrapod origins. *Nature* **444**, 199–202 (2006).

57. Daeschler, E. B., Shubin, N. H. & Jenkins, F. A. Jr A Devonian tetrapod-like fish and the evolution of the tetrapod body plan. *Nature* **440**, 757–763 (2006).

58. Barton, R. A. & Harvey, P. H. Mosaic evolution of brain structure in mammals. *Nature* **405**, 1055–1058 (2000).

59. Schlosser, G. Mosaic evolution of neural development in anurans: acceleration of spinal cord development in the direct developing frog *Eleutherodactylus coqui*. *Anat. Embryol.* **206**, 215–227 (2003).

60. Felice, R. N. & Goswami, A. Developmental origins of mosaic evolution in the avian cranium. *Proc. Natl Acad. Sci. USA* **115**, 555–560 (2018).

61. Felice, R. N. et al. Evolutionary integration and modularity in the archosaur cranium. *Integr. Comp. Biol.* **59**, 371–382 (2019).

62. Watanabe, A. et al. Ecomorphological diversification in squamates from conserved pattern of cranial integration. *Proc. Natl Acad. Sci. USA* **116**, 14688–14697 (2019).

63. Owen, R. On the *Archaeopteryx* of Von Meyer, with a description of the fossil remains of a long-tailed species from the lithographic stone of Solnhofen. *Philos. Trans. R. Soc. Lond.* **153**, 33–47 (1863).

64. Paluh, D. J., Stanley, E. L. & Blackburn, D. C. Evolution of hyperossification expands skull diversity in frogs. *Proc. Natl Acad. Sci. USA* **117**, 8554–8562 (2020).

65. Gomez-Mestre, I., Pyron, R. A. & Wiens, J. J. Phylogenetic analyses reveal unexpected patterns in the evolution of reproductive modes in frogs. *Evolution* **66**, 3687–3700 (2012).

66. Nevo, E. Adaptive convergence and divergence of subterranean mammals. *Annu. Rev. Ecol. Syst.* **10**, 269–308 (1979).

67. Nevo, E. Mammalian evolution underground. The ecological-genetic-phenetic interfaces. *Acta Theriol.* **3**, 9–31 (1995).

68. Vogel, S. Life's Devices: The Physical World of Animals and Plants (Princeton Univ. Press, 1988).

69. Sansalone, G. et al. Impact of transition to a subterranean lifestyle on morphological disparity and integration in talpid moles (Mammalia, Talpidae). *BMC Evol. Biol.* **19**, 1–15 (2019).

70. Nauwelaerts, S., Ramsay, J. & Aerts, P. Morphological correlates of aquatic and terrestrial locomotion in a semi-aquatic frog, Rana esculenta: no evidence for a design conflict. *J. Anat.* **210**, 304–317 (2007).

71. Sherratt, E., Gower, D. J., Klingenberg, C. P. & Wilkinson, M. Evolution of cranial shape in caecilians (Amphibia: Gymnophiona). *Evol. Biol.* **41**, 528–545 (2014).

72. Cardini, A., Polly, P. D., Dawson, R. & Milne, N. Why the long face? Kangaroos and wallabies follow the same 'rule' of cranial evolutionary allometry (CREA) as placental. *Evol. Biol.* **42**, 169–176 (2015).

73. Yeh, J. The effect of miniaturized body size on skeletal morphology in frogs. *Evolution* **56**, 628–641 (2002).

74. Wells, K. D. *The Ecology and Behavior of Amphibians* (University of Chicago Press, 2010).

75. Emerson, S. B. Skull shape in frogs: correlations with diet. *Herpetologica* **41**, 177–188 (1985).

76. Carreño, C. A. & Nishikawa, K. C. Aquatic feeding in pipid frogs: the use of suction for prey capture. *J. Exp. Biol.* **213**, 2001–2008 (2010).

77. Fernandez, E., Irish, F. & Cundall, D. How a frog, *Pipa pipa*, succeeds or fails in catching fish. *Copeia* **105**, 108–119 (2017).

78. Herrel, A. et al. in *Feeding in Vertebrates: Evolution, Morphology, Behavior, Biomechanics* (eds. Bels, V. & Whishaw, I. Q.) (Springer, 2019).

79. Bardua, C. et al. Evolutionary integration of the frog cranium. *Evolution* **74**, 1200–1215 (2020).

80. Bon, M., Bardua, C., Goswami, A. & Fabre, A.-C. Cranial integration in the fire salamander, *Salamandra salamandra* (Caudata: Salamandridae). *Biol. J. Linn. Soc.* **130**, 178–194 (2020).

81. Fabre, A. et al. Metamorphosis and the evolution of morphological diversity in salamanders. *Nat. Ecol. Evol.* **4**, 1129–1140 (2020).

82. Nishikawa, K. C. in *Feeding: Form, Function and Evolution in Tetrapod Vertebrates* (ed. Schwenk, K.) (Academic Press, 2000).

83. Trueb, L. & Gans, C. Feeding specializations of the Mexican burrowing toad, *Rhinophrynus dorsalis* (Anura: Rhinophrynidae). *J. Zool.* **199**, 189–208 (1983).

84. Nishikawa, K. C., Kier, W. M. & Smith, K. K. Morphology and mechanics of tongue movement in the African pig-nosed frog *Hemisus marmoratum*: a muscular hydrostatic model. *J. Exp. Biol.* **202**, 771–80 (1999).

85. Henrici, A. C. Digging through the past: the evolutionary history of burrowing and underground feeding in rhinophrynid anurans. *Palaeobiodivers. Palaeoenviron.* **96**, 97–109 (2015).

86. Van Dijk, D. E. Osteology of the ranoid burrowing African anurans *Breviceps* and *Hemisu*. *Afr. Zool.* **36**, 137–141 (2001).

87. Womack, M. C., Christensen-Dalsgaard, J., Coloma, L. A. & Hoke, K. L. Sensitive high-frequency hearing in earless and partially eared harlequin frogs (*Atelopus*). *J. Exp. Biol.* **221**, 1–8 (2018).

88. Boistel, R. et al. How minute sooglossid frogs hear without a middle ear. *Proc. Natl Acad. Sci. USA* **110**, 15360–15364 (2013).

89. Womack, M. C., Stynoski, J. L., Voyles, M. K., Coloma, L. A. & Hoke, K. L. Prolonged middle ear development in *Rhinella horribilis*. *J. Morphol.* **279**, 1518–1523 (2018).

90. Womack, M. C., Christensen-Dalsgaard, J., Coloma, L. A., Chaparro, J. C. & Hoke, K. L. Earless toads sense low frequencies but miss the high notes. *Proc. R. Soc. B Biol. Sci.* **284**, 20171670 (2017).

91. Hetherington, T. E. in *The Evolutionary Biology of Hearing* (eds. Webster, D. B., Fay, R. R. & Popper, A. N.) (Springer, 1992).

92. Hanken, J., Klymkowsky, M. W., Summers, C. H., Seufert, D. W. & Ingebrigtsen, N. Cranial ontogeny in the direct-developing frog, *Eleutherodactylus coqui* (Anura: Leptodactylidae), analyzed using whole-mount Immunohistochemistry. *J. Morphol.* **211**, 95–118 (1992).

93. Hanken, J., Klymkowsky, M. W., Alley, K. E. & Jennings, D. H. Jaw muscle development as evidence for embryonic re patterning in direct-developing frogs. *Proc. R. Soc. B Biol. Sci.* **264**, 1349–1354 (1997).

94. Wray, G. A. & Raff, R. A. The evolution of developmental strategy in marine invertebrates. *Trends Ecol. Evol.* **6**, 45–50 (1991).

95. Watkins, T. B. A quantitative genetic test of adaptive decoupling across metamorphosis for locomotor and life-history traits in the Pacific tree frog, *Hyla regilla*. *Evolution* **55**, 1668–1677 (2001).

96. Wilson, A. D. M. & Krause, J. Personality and metamorphosis: is behavioral variation consistent across ontogenetic niche shifts? *Behav. Ecol.* **23**, 1316–1323 (2012).

97. O'Reilly, J. C., Deban, S. M. & Nishikawa, K. C. in *Topics in Functional and Ecological Vertebrate Morphology: A Tribute to Frits de Vree* (eds. Aerts, P., D'Août, K., Herrel, A. & van Damme, R.) (Shaker Publishing, 2002).

98. Philips, P. C. Genetic constraints at the metamorphic boundary: morphological development in the wood frog, *Rana sylvatica*. *J. Evol. Biol.* **11**, 453–463 (1998).

99. Johansson, F., Lederer, B. & Lind, M. I. Trait performance correlations across life stages under environmental stress conditions in the common frog, Rana temporaria. *PLoS ONE* **5**, e11680 (2010).
100. Wassersug, R. J. The adaptive significance of the tadpole stage with comments on the maintenance of complex life cycles in anurans. *Am. Zool.* **15**, 405–417 (1975).

101. Vassilieva, A. B. Heterochronies in the cranial development of Asian tree frogs (Amphibia: Anura: Rhacophoridae) with different life histories. *Dokl. Biol. Sci.* **473**, 110–113 (2017).

102. Kerney, R., Meegaskumbura, M., Manamendra-Arachchi, K. & Hanken, J. Cranial ontogeny in *Philautus silus* (Anura: Ranidae: Rhacophorinae) reveals few similarities with other direct-developing anurans. *J. Morphol.* **268**, 715–725 (2007).

103. Heatwole, H. & Davies, M. (eds.) *Amphibian biology (volume 5), osteology*. (Surrey Beatty & Sons, 2003).

104. Hanken, J. & Hall, B. K. Skull development during anuran metamorphosis: I. Early development of the first three bones to form—the exoccipital, the parasphenoid, and the frontoparietal. *J. Morphol.* **195**, 247–256 (1988).

105. Fink, W. L. The conceptual relationship between ontogeny and phylogeny. *Paleobiology* **8**, 254–264 (1982).

106. Strathmann, R. R. in *Echinoderm Phylogeny and Evolutionary Biology* (eds. Paul, C. R. C. & Smith, A. B.) (Clarendon Press, 1988).

107. Laloy, F. et al. A re-interpretation of the Eocene anuran *Thaumastosaurus* based on MicroCT examination of a “mummified” specimen. *PLoS ONE* **8**, e74874 (2013).

108. Frost, D. R. et al. The amphibian tree of life. *Bull. Am. Mus. Nat. Hist.* **297**, 1–370 (2006).

109. Quental, T. B. & Marshall, C. R. Diversity dynamics: molecular phylogenies need the fossil record. *Trends Ecol. Evol.* **25**, 435–441 (2010).

110. Slater, G. J. & Harmon, L. J. Unifying fossils and phylogenies for comparative analyses of diversification and trait evolution. *Methods Ecol. Evol.* **4**, 699–702 (2013).

111. Volume Graphics. VGStudio MAX v. 2.0 (Volume Graphics GmbH, 2001).

112. Bardua, C., Felice, R. N., Watanabe, A., Fabre, A.-C. & Goswami, A. A practical guide to sliding and surface semilandmarks in morphometric analyses. *Integr. Org. Biol.* **1**, 1–34 (2019).

113. Jetz, W. & Pyron, R. A. The interplay of past diversification and evolutionary isolation with present imperilment across the amphibian tree of life. *Nat. Ecol. Evol.* **2**, 850–858 (2018).

114. Paradis, E., Claude, J. & Strimmer, K. APE: analyses of phylogenetics and evolution in R language. *Bioinformatics* **20**, 289–290 (2004).

115. Paradis, E. & Schliep, K. ape 5.0: an environment for modern phylogenetics and evolutionary analyses in R. *Bioinformatics* **35**, 526–528 (2018).

116. Wiley, D. F. et al. Evolutionary morphing. In *Proc. Visualization Conference* (IEEE, 2005).

117. Schlager, S. in *Statistical Shape and Deformation Analysis* (eds. Zheng, G., Li, S. & Szekely, G.) (Academic Press, 2017).

118. Cardini, A. Left, right or both? Estimating and improving accuracy of one-side-only geometric morphometric analyses of cranial variation. *J. Zool. Syst. Evol.* **55**, 1–10 (2016).

119. Marshall, A. F. et al. High-density three-dimensional morphometric analyses support conserved static (intraspecific) modularity in caecilian (Amphibia: Gymnophiona) crania. *Biol. J. Linn. Soc.* **126**, 721–742 (2019).

120. Bossuyt, F. & Milinkovitch, M. C. Convergent adaptive radiations in Madagascan and Asian ranid frogs reveal covariation between larval and adult traits. *Proc. Natl Acad. Sci. USA* **97**, 6585–90 (2000).

121. Young, J. E., Christian, K. A., Donnellan, S. C., Tracy, C. R. & Parry, D. Comparative analysis of cutaneous evaporative water loss in frogs demonstrates correlation with ecological habits. *Physiol. Biochem. Zool.* **78**, 847–856 (2005).

122. Portik, D. M. & Blackburn, D. C. The evolution of reproductive diversity in Afrobatrachia: a phylogenetic comparative analysis of an extensive radiation of African frogs. *Evolution* **70**, 2017–2032 (2016).

123. Scott, E. A phylogeny of ranid frogs (Anura: Ranoidea: Ranidae), based on a simultaneous analysis of morphological and molecular data. *Cladistics* **21**, 507–574 (2005).

124. Adams, D. C. & Otárola-Castillo, E. geomorph: an R package for the collection and analysis of geometric morphometric shape data. *Methods Ecol. Evol.* **4**, 393–399 (2013).

125. Revell, L. J. phytools: an R package for phylogenetic comparative biology (and other things). *Methods Ecol. Evol.* **3**, 217–223 (2012).

126. Clavel, J., Escarguel, G. & Merceron, G. mvmorph: an r package for fitting multivariate evolutionary models to morphometric data. *Methods Ecol. Evol.* **6**, 1311–1319 (2015).

127. Clavel, J., Aristide, L. & Morlon, H. A penalized likelihood framework for high-dimensional phylogenetic comparative methods and an application to new-world monkeys brain evolution. *Syst. Biol.* **68**, 93–116 (2019).

128. Clavel, J. & Morlon, H. Reliable phylogenetic regressions for multivariate comparative data: illustration with the MANOVA and application to the effect of diet on mandible morphology in phyllosomid bats. *Syst. Biol.* **69**, 927–943 (2020).

129. Housworth, E. A., Martins, E. P. & Lynch, M. The phylogenetic mixed model. *Am. Nat.* **163**, 84–96 (2004).

130. Revell, L. J. Phylogenetic signal and linear regression on species data. *Methods Ecol. Evol.* **1**, 319–329 (2010).

131. Freckleton, R. P., Harvey, P. H. & Pagel, M. Phylogenetic analysis and comparative data: a test and review of evidence. *Am. Nat.* **160**, 712–726 (2002).

132. Goolsby, E. W., Bruggeman, J. & Ane, C. Rphylopars: phylogenetic comparative tools for missing data and within-species variation. *R package version* 0.2.11 https://CRAN.R-project.org/package=Rphylopars (2019).

133. Goolsby, E. W., Bruggeman, J. & Ane, C. Rphylopars: fast multivariate phylogenetic comparative methods for missing data and within-species variation. *Methods Ecol. Evol.* **8**, 22–27 (2017).

134. Bardua, C. & Goswami, A. Frog skull shape data for modularity and macroevolution. https://doi.org/10.5281/zenodo.4619880 (2020).

## Acknowledgements

We would like to thank Ryan Felice, Akinobu Watanabe and Benjamin Poole for help with analyses, and Katie Thomas, Jeff Streicher, David Gower, and members of the Goswami Lab for helpful discussions. We would like to thank Anne-Marie Ohler for help with CT scanning at the MNHN collection, Paris, and Kristin Mahlow, Nadia Fröbisch, Frank Tillack, Johannes Müller and Mark-Oliver Rödel for facilitating the CT scanning at the MfN, Berlin. This work was supported by ERC grant STG-2014-637171 to A.G. This research also received support from the SYNTHESYS Project http://www.SYNTHESYS.info/ (FR-TAF-5583) and Horizon2020 Marie Skłodowska-Curie actions (EVOTOOLS, H2020-MCSA-IF-2017-797373 to J.C.), which are financed by European Community Research Infrastructure Action under the FP7 Integrating Activities Programme. Data collection performed at the University of Florida was supported in part by NSF (DBI-1701714) to D.C.B. and E.L.S. Funding for scanning *Chirodamis rufescens* MCZ A 136738, *Rhacophorus reinwardtii* FMNH 252384 and *Theloderma stellatum* FMNH 253619 was provided by an NSF grant (DEB-0345885) to Dr. Hanken, and funding for image processing was provided by an NSF Digital Libraries Initiative grant to Dr. Timothy Rowe of the University of Texas at Austin. Funding for scanning *Incilius (Bufo) periglenes* TNHC 39447 was provided by Dr. Cossel and funding for image processing was provided by the High-Resolution X-ray CT Facility. We thank Renaud Boistel for use of the *Thaumastosaurus gezei* MNHN QU 17279 scan.

## Author contributions

C.B. primarily collected the data, A.G., C.B., A.-C.F., J.C. and D.C.B. designed the analyses. C.B. and J.C. performed the analyses. C.B., K.D., E.L.S., M.B. and D.C.B. generated and reconstructed the scans. A.G. C.B., J.C., A.-C.F. and D.C.B. wrote the manuscript.

## Competing interests

The authors declare no competing interests.

## Additional information

**Supplementary information** The online version contains supplementary material available at https://doi.org/10.1038/s41467-021-22792-y.

**Correspondence** and requests for materials should be addressed to A.G.

**Peer review information** *Nature Communications* thanks the anonymous reviewers for their contribution to the peer review of this work.

**Reprints and permission information** is available at http://www.nature.com/reprints

**Publisher's note** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

<figure><img src="image_10.png" /></figure>

**Open Access** This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/.

© Crown 2021