The radiation continuum and the evolution of frog diversity – Nature.com

This study contained no experimental component or data collection on live animals and so no ethical oversight was necessary.

We measured 4628 adult museum specimens of 1234 species from around the world. Most of these data were novel, whereas 901 specimens from 194 species came from previously published datasets1,20,21,22. Our sample included 51 of 54 anuran families19. The three remaining families (Calyptocephalellidae, Ceuthomantidae, and Nasikabatrachidae) are scarce in museum collections. We chose species within families based on their availability in museum collections, with species sampling proportional to the described species diversity of each family (r=0.923). However, for eight families we were only able to sample a single species, which prevented calculating rates of morphological evolution. Thus, we excluded them (to total 1226 species from 43 families) from those analyses and all downstream analyses based on those rates. We also note that some studies of rates of morphological evolution have removed clades with low numbers of species (e.g., less than four8). In our dataset, 11 families had between 24 species sampled for morphological data. However, some of these families have four or fewer total extant species, and thus excluding these families would result in biasing our analyses to ignore clades with low species richness. Moreover, while lower sampling may increase the variance in estimates of a clades true rate of evolution, such estimates are unbiased1. Finally, to reduce potential effects of sexual-size dimorphism on our sampling76,77,78, we measured male specimens when possible (89% of all specimens sampled; 82% of our sampled species were represented only by males). Males tend to be better represented in collections than females, presumably because of their calling behavior. We include all raw intraspecific data as Supplementary Data1 and species means, sample sizes, standard deviations, and standard errors as Supplementary Data2.

We quantified body shape using linear, area, and volumetric measurements of traits that are ecologically and functionally relevant to locomotion and microhabitat use21,22,27,28. First, we measured snout-vent length, head length, head width, upper arm, forearm, hand, thigh, crus, tarsus, and foot lengths to the nearest 0.01mm using a Mitutoyo digital caliper (Kanagawa, Japan). We took each measurement only once, as our measurements were highly precise; preliminary repeated measurements showed a coefficient of variation of less than 0.03 for all measurements, with most <0.015. We summed the linear limb element measurements together (i.e., front limb length, hindlimb length). Second, we photographed the foot and hand of each specimen and measured the areas of digit tips on both the front and hind limb, interdigital webbing of the hind limb, and the inner metatarsal tubercle using ImageJ79. We summed the areas of the digit tips separately for the front and hind limbs and interdigital webbing across the foot. Detailed descriptions of all measurements are given in Supplementary Table4.

Finally, we quantified leg muscle volume using external linear measurements. We used thigh and crus muscle volume among the traits for characterizing anuran body shape. Muscle mass is strongly related to locomotor performance and microhabitat use in anurans21,22,26,55. However, we could not calculate mass by dissecting muscle tissue from museum specimens at this scale of sampling. Thus, we estimated leg muscle volume, which should scale 1:1 with mass80 and could be quantified using external linear measurements. We estimated muscle volume of the thigh and crus separately, considering each leg segment as two cones sharing an elliptical base (i.e., the approximate cross-sectional area of the underlying muscle). We measured the depth and width of the thigh and crus at their mid-points as the axes of the ellipse. To ensure our approximation of muscle volume adequately represented its mass, we took advantage of the previously published subset of our data (641 specimens from 132 species21,22) that included masses of dissected thigh and crus muscles. For these specimens, we natural-log transformed (ln) thigh and crus masses and volumes to linearize the relationship, then checked the correlation between thigh (or crus) muscle mass against estimated volume at the specimen level. We found that mass and volume were strongly correlated (r=0.974 and 0.965 for thigh and crus, respectively), which suggests that our volume approximation accurately represents muscle mass.

We lacked width and depth measurements but had muscle masses and lengths for thigh and crus for 238 specimens from 49 species. To include these 238 specimens in our analyses, we estimated the muscle cross-sectional area, which we could then use with observed leg segment length to estimate volume. We thus regressed the ln thigh (or crus) cross-sectional area on ln mass for the aforementioned 641 specimens with both data. We then used this model to predict cross-sectional areas for the 238 specimens that lacked width and depth measurements. These regressions showed that the mass of thigh and crus strongly predicted cross-sectional area (R2=0.949 and 0.931 for thigh and crus, respectively; Supplementary Fig.6).

We used previously published microhabitat data23 and additional natural history descriptions to classify most species to microhabitats; new classifications determined for this study are provided in Supplementary Data3. Most species can be categorized into eight different microhabitat states22,23. Four of these states are base microhabitat states that broadly categorize adult frog ecology: aquatic (found primarily in water), arboreal (found primarily in trees and brushes), burrowing (found primarily in self-dug burrows), and terrestrial (found primarily on the surface or under shallow leaf litter). Three additional categories combine terrestrial microhabitats with others, when ecological descriptions indicate that species spend time in both microhabitats: semi-aquatic, semi-arboreal, and semi-burrowing. The torrential microhabitat is characterized by occupying vegetation and rocks along high-gradient streams and rushing waters, thus combining aquatic and arboreal states.

We used the posterior distribution of time-calibrated, multi-locus trees generated by Jetz and Pyron29 for comparative analyses. We chose this phylogeny because it included all species in our morphological dataset. Whereas most more recent phylogenies81,82 may have more molecular data per species and potentially more accurate clade ages, they have far fewer taxa (i.e., they would leave out about 90% of our species). Moreover, recent comparative analyses of diversification rates in anuran families show similar results regardless of the tree used to calculate clade ages1.

We first pruned the posterior distribution to include only anuran species with genetic data (3449 species), because trees with taxa placed based on taxonomy alone may inflate rates of phenotypic evolution83. We used tools available at VertLife (www.vertlife.org/phylosubsets; date accessed: 25 January 2021) to download a random draw of 1000 trees. We then used TreeAnnotator84 to calculate the maximum-clade credibility (MCC) topology and summarize branch lengths in millions of years, doing so with the Common Ancestor heights option. This option generally produces more accurate estimates of clade age than mean branch lengths85.

Previous analyses have shown that adaptive morphological diversification in frogs is often unrelated to body size1,21,22,86. Thus, to focus on shape-based morphology, we size-corrected each trait using log-shape ratios30,31,32, wherein we divided variables by SVL and then ln-transformed the resulting ratios32. Traditional log-shape ratios consider size as the geometric mean of all morphological variables31. However, we only used SVL as a metric of size, given that we measured the other variables precisely because we expected them to differ based on ecology21,22. By contrast, SVL does not differ based on microhabitat22 and can differ greatly among species with similar body shape (e.g., refs. 57,63). For area and volume measurements, we took the square root or cube root of the raw values prior to size-correction to ensure equal scaling across variables80. We performed all size corrections on raw (i.e., intraspecific) data, then calculated species means from the size-corrected intraspecific values. For this and nearly all other analyses, we used the R computing environment87, version 4.1.0.

To ensure that size standardization did not affect pPC axis interpretation, we also performed interspecific size-correction using residuals33 of each trait regressed against SVL, using phytools in R. We then conducted a phylogenetic PCA on these residuals. We found high correlations between the eigenvectors of each PC axis resulting from this alternative method of size standardization and our preferred ratio method (rMantel=0.987; P<0.001). Thus, the method of size standardization is unlikely to change our interpretation of downstream analyses30. Furthermore, several papers have cautioned against treating residuals from linear regressions as data88,89,90. For these reasons and for brevity, we only present results obtained from the log-shape ratio method of size-correction.

We summarized diversity in body form using a phylogenetic principal components analysis (PCA) on species means, as implemented in the phytools package91, version 0.747. We included size-corrected measurements described above of head length and width, front and hindlimb lengths, volumes of the thigh and crus muscles, areas of foot webbing and theinner metatarsal tubercle, and area of the digit tips of the foot and hand. We assumed a Brownian motion model of evolution, and we conducted the PCA on the phenotypic covariance matrix, given our prior standardization of all variables to the same scale and units92. We also performed a non-phylogenetic PCA to ensure that the interpretation of body form was insensitive to analytical method92,93. We compared the results of these two types of PCA by conducting a Mantel test (10,000 permutations) on the PCA eigenvectors, as implemented in the package vegan94 version 2.5.7. This analysis showed a strong correlation (rMantel=0.885; P<0.001) between phylogenetic and non-phylogenetic PCAs. Thus, PCA method seemed unlikely to affect downstream analyses or interpretations, so we used the resulting phylogenetic PCA scores for later analyses of morphological diversity.

Our approach necessitated comparing many different clades. We chose families as the unit of analysis. Anuran families range from 1 to >1000 species and show substantial variation in diversification rates23. Families are also sufficient in number (54 total) to examine patterns with robust statistics. At shallower taxonomic levels (e.g., genera), we may see similar patterns as families57 but would generally have too few species per clade to robustly calculate rates of phenotypic evolution. In contrast, anurans have relatively few formally named clades above families81, which would leave a limited sample size for statistical analysis.

We recognize that using taxonomy to define clades may impact analyses95,96 (but also see respective responses97,98). To avoid possible biases from clade selection, we also used clades of the same age as alternative units for analysis96. We selected age-based clades by considering the most inclusive clade of a given age or younger. With the tree used here29, a threshold for clade selection much younger than 80 million years would return many groups with few species, limiting variation in net diversification rates. In contrast, a threshold much older than 120 would not return enough clades for robust statistical analysis (e.g., the 120million-year threshold produced 19 clades; Supplementary Fig.4). We therefore repeated the radiation-space analyses described below on clades defined by ages of 80, 100, and 120 million years old.

We estimated morphological diversity of all anurans, families, age-defined clades, and radiation-space quadrants (see below). We defined morphological diversity as the volume of n-dimensional morphospace occupied by a group of species. We used two approaches: a convex-hull volume35 and a hypervolume36. Convex hulls are effectively n-dimensional ranges35. They likely overestimate shape volume because they are sensitive to outliers and are unable to detect holesgaps between observationsin n-dimensional space36. Hypervolume methods use machine-learning algorithms to determine boundaries around points in n-dimensional space and are able to detect and exclude outliers and holes36,99,100. Hypervolumes likely underestimate shape and volume depending on the nature of the dataset. For these reasons, the convex hull and hypervolume approaches likely produce a maximal and minimal volume estimate (respectively) of morphological diversity. In consequence, correspondence of results from these two methods should indicate insensitivity to methods of quantifying morphological diversity.

Both methods assume that each axis considered is orthogonal to others, so we used scores from our phylogenetic PCA (pPCA) as data for morphospace calculations. Because both methods are computationally burdensome, we limited analyses to the first five pPC axes. We found in preliminary analyses that five was the best compromise between comprehensiveness and computation time. Moreover, a scree plot (Supplementary Fig.7) showed a considerable drop in variation explained after five axes101,102. These first five axes collectively explained 92.4% of the morphological variation in our dataset (Supplementary Table1). Most importantly, our results were similar for more (six) and fewer (four) dimensions (Supplementary Table5). To estimate the convex hull, we used the Quickhull algorithm implemented in the geometry package103, version 0.4.5. To estimate the hypervolume, we used the one-class support vector machine method as implemented in hypervolume100, version 2.0.12.

We estimated multivariate rates of morphologicalevolution for families and age-defined clades using the method of Adams37. This method calculates a single Brownian-motion rate of evolution that accounts for correlations among characters. Brownian motion is the simplest and most general model of continuous-trait macroevolution and is consistent with many different underlying evolutionary scenarios (e.g., stabilizing selection with randomly evolving optima)46,104,105. Moreover, previous work has shown that the evolution of these same traits is consistent with a Brownian-motion model in 217 species across many families1. Furthermore, given that our sampling of species within families averaged 25% of each familys extant species richness, we emphasize that incomplete clade sampling does not bias this metric. That same previous study1 (of anurans, with the same traits) used simulations to show that sampling as low as 2.3% of total species diversity has no effect on the accuracy of rate estimation.

We present our raw estimated rates as Supplementary Data4. However, comparing rates estimated here to previously published rates for other groups is incredibly challenging. While the method we used37 is increasingly employed for estimating multivariate rates of phenotypic evolution92, such rate estimates are influenced by different methods for size standardization (e.g., ratios, residuals, General Procrustes Analyses in geometric morphometrics106) and different numbers of traits107.

We followed the classification of AmphibiaWeb19 for defining families and counting their species diversity. For clades from 80, 100, and 120 million-year time slices, we established species richness using the full tree from Jetz and Pyron29, which included all known species at the time of their analysis. This tree provides an underestimate of current species richness19, but this step was necessary to calculate the species diversity of time-sliced clades when genera were separated into multiple clades. It also allowed us to include the species diversity of genera unsampled in the genetictree of Jetz and Pyron29, which we used for all other analyses.

We initially estimated net diversification rates using the method-of-moments estimator38. This method only requires species richness and clade ages, which are available for all anuran families. Moreover, recent simulation studies show that this method is accurate under many diversification scenarios, including faster rates in younger clades, rate variation over time within clades, and rate heterogeneity across subclades68,69,70. We recognize that many other methods of calculating diversification rates are available. However, the estimator we used allows as many different rates as families, far more than other methods typically find (e.g., see refs. 108,109). Moreover, this method allows one to estimate the potential effect of extinction on downstream analyses: we can compare how our results (potentially) differ based on low or high extinction fractions. This may be particularly important in anurans, whose oldest families may have low diversity due to high historical extinction rates110,111. Yet simple diversification metrics (like the method-of-moments estimator we use) may avoid problems associated with trying to extract too much information from phylogenies of extant taxa72. We also emphasize that adaptive radiation may be a temporal phenomenon (i.e., groups characterized as adaptive radiations now may not have been 100 million years ago), as are other macroevolutionary patterns. However, what we see in present-day groups is what we study here: we focus on what led to current species and phenotypic diversity, not how past adaptive radiation led to diversity we no longer see. Thus, using a diversification metric that integrates over the history of clades to the present day is what is most relevant to our study.

We also compared these rates (based on species diversity and ages) with birth-death rates (based on branch lengths) estimated by Moen et al.1. Because the birth-death rates could only be estimated for the 38 families with sufficient sampling (at least five species in Jetz and Pyron29), we added our originally estimated method-of-moments rates under stem ages and medium extinction fraction for the remaining five families to total 43 families, as in our other diversification-rate analyses. We found that the birth-death rates and method-of-moments rates were highly correlated (Supplementary Table2). Moreover, our radiation-space results were broadly similar using birth-death rates for diversification (Supplementary Fig.2). However, we prefer the method-of-moments estimates because we could include all 43 anuran families in this study under a single method of rate estimation.

To be consistent with our morphological analyses, we calculated the stem and crown ages for each family from our MCC consensus tree. Other phylogenies give younger ages for anuran families81,82. However, recent diversification analyses using ages from both Jetz and Pyron29 and Feng et al.81 showed high correlations in rates across families1. For example, rates based on stem ages and an extinction fraction of 0.5 showed a correlation of r=0.953 between the two trees. Here, we calculated rates using three extinction fractions (; 0, 0.5, and 0.9), following standard practice112,113,114. We present results based on rates estimated using moderate extinction fractions (=0.5). Low and high extinction fractions gave similar results in downstream analyses (Supplementary Fig.3). Moreover, we present results based on stem ages, which are estimated from the origination of the clade and are less sensitive to sampling density than crown ages115. Results for the latter were highly similar (Supplementary Fig.3).

Moen and Wiens23 showed a strong correlation between species diversity and net diversification rates of anuran families. Here, we re-evaluated this correlation for the 43 families examined in this paper, given updated species richness of families (i.e., >10% of anuran species have been described since 2016; ref. 19). We then tested the relationship between rates of multivariate morphological evolution and morphological diversity across families. We estimated morphological diversity using five-dimensional convex hulls and hypervolumes, as described above. Here, we only examined families with six or more species measured (n=27), because n+1 observations are required to define an n-dimensional volume. We then calculated the fifth root of the resulting values. For all variables, we ln-transformed, mean-centered, and scaled them to unit variance (Supplementary Data5). We then used phylogenetic generalized least-squares (pGLS) correlations to estimate correlations between morphological diversity and morphological rates of evolution, and net diversification rates and species richness. To be consistent with our calculation of rates, we again used the phylogeny of Jetz and Pyron29 for our pGLS analyses. However, we expect results to be highly similar with other recent phylogenomic trees81,82, given that pGLS is highly robust to tree misspecification116. We calculated pGLS correlations following Rohlf117 and using a custom R script from Moen et al.21.

We next tested the strength of the relationship between rates of net diversification and morphological evolution. This allowed us to examine whether rates were strongly correlated (producing a linear radiation continuum) versus weakly correlated or uncorrelated (yielding a two-dimensional radiation space). We calculated pGLS correlations on the mean-centered and scaled rates of net diversification and morphological evolution (n=43), as described above. We also visualized the relationship between rates by plotting them on the phylogeny with the ggtree R package118, version 3.0.2, with ancestral states estimated by maximum likelihood in phytools91. Given that we found a weak correlation (see Results), we next describe the continuum along its two dimensions.

To characterize an adaptive-radiation space, we separated clades into quadrants by rates of net diversification and morphological evolution, where the origin (0, 0) represented mean values among clades for both rates. Clades with rates of net diversification and morphological evolution >0 were assigned to the adaptive-radiation quadrant. Clades with rates of net diversification and morphological evolution <0 were considered non-adaptive non-radiations. Clades with net diversification rates >0, but rates of morphological evolution <0, were placed in the non-adaptive radiation quadrant. Clades with net diversification rates <0, but rates of morphological evolution >0, were considered adaptive non-radiations.

We also repeated clade assignments after redefining the quadrant boundaries as medians of rates. This alternative scheme allowed us to explicitly examine how robust our results were to quadrant limits. Because all analyses (i.e., families and clades extracted at 80, 100, and 120 million-year time slices) had an odd number of observations, the median clade always straddled at least two quadrants. To avoid omitting any clades, we split clades with median values for either net diversification or morphological rates equally between the quadrants these clades straddled. For morphospace volume calculations (see next section), this meant randomly assigning half (when straddling two quadrants) or a quarter (when straddling all four) of the median-clade species to each quadrant the clade straddled when estimating volumes. For species diversity, we simply divided the number of species in the clade by two (or four) and added them to the number of species observed in the quadrants they straddled.

We then characterized the phylogenetic distribution of evolutionary dynamics (i.e., our four quadrant types) by calculating the D statistic of phylogenetic signal for binary traits45 as implemented in the caper package119, version 1.0.4. We conducted four analyses, one for each radiation type, with each analysis estimating D for a binary trait consisting of one radiation type versus all others (e.g., for non-adaptive radiation, a trait with one state as non-adaptive radiation and the other state as all other types). A D of 0 or lower (i.e., negative) would indicate phylogenetic clustering, whereas a D of 1 or higher would indicate a random (D=1) or overdispersed distribution (D>1)45. Thus, we tested for a significant deviation from D=0.0 (which would suggest significant random distribution or overdispersion) and from D=1.0 (which would suggest significant clustering). We only conducted this analysis for quadrants delimited by mean evolutionary rates for families, given we found that no quadrant type showed a D significantly different from 0 or 1 (Supplementary Table3). We did not expect different results for other ways of characterizing clades or the radiation space.

Our primary goal was to determine the role of adaptive radiation in driving diversity in a major clade. For this goal we needed to first quantify total and quadrant-specific species and morphological diversity, then the proportion of diversity each quadrant of radiation space contained. For species diversity, we tallied the total species richness of the sampled families from AmphibiaWeb19 to represent total anuran species diversity. This diversity (7359 species) represents >99% of extant described anuran species (7426 species). Thus our results for these 43 anuran families should basically apply to all Anura. We then calculated species diversity for each quadrant by summing the currently described species richness of all families within that quadrant. We divided each quadrant total by the total diversity we analyzed (7359 species) to calculate the proportion of total diversity explained by each of the four types of radiations.

We quantified total morphological diversity as the morphospace volume occupied by all species in our morphological dataset (i.e., the 1226 species for which we could calculate rates of evolution). We then divided the pPC scores into four subsets of species, one subset for each quadrant of the adaptive radiation plane. Each subset included all the species from the clades that we categorized as belonging to that quadrant. We then estimated each quadrants morphological diversity using the methods described above.

We divided each quadrants volume by the total anuran volume to calculate the contribution of each radiation type to total anuran morphological diversity. Unlike species diversity, where each quadrants species contribute independently to total species diversity, morphospaces of different quadrants may overlap. When this occurs, the sum of quadrant percentages may total more than 100%. Alternatively, quadrant percentages may not sum to 100% if quadrant morphospaces occur in mutually exclusive regions of the total anuran morphospace (i.e., gaps between quadrants within the total anuran morphospace)99.

Both net diversification rates and rates of phenotypic evolution include time in their estimation. While time is directly used in the calculation of net diversification rate, it is involved in morphological rates through phylogenetic branch lengths. Such a shared dimension could, in principle, lead to similarity in these two types of rates (e.g., a family with a high net diversification rate could have a high rate of phenotypic evolution). Moreover, both rates often show a negative relationship with time across many groups of organisms120. Thus, we further explored the potential effect of shared time on our net diversification and multivariate morphological rate estimates. For brevity, we circumscribed these analyses to include only net diversification rates estimated using stem ages and moderate extinction fractions (=0.5). First, we assessed the relationship between age and rate by using phylogenetic generalized least squares (pGLS) regression under Brownian motion, as implemented in the R package phylolm, version 2.6.2121. We regressed net diversification rates on stem age (i.e., rather than crown age) because it was the age used to calculate the rates on which we focused here. In contrast, we regressed rates of phenotypic evolution on crown age, given that only the crown phylogeny of each family was used for estimating rates of evolution (using stem ages led to even weaker relationships). These regressions showed weak but statistically significant relationships between each rate and their respective family age estimates. Surprisingly, morphological rate of evolution had a significant positive slope (=0.0140.006; R2Adj=0.077; P=0.040), contrasting with the typically negative relationship122,123,124. Net diversification rate showed a significant negative relationship with time (=0.0200.005; R2Adj=0.231; P<0.001), as expected when regressing a ratio against its denominator125,126.

We then assessed whether time-independent net diversification rates and morphological rates of evolution were correlated. We did this by calculating residuals from each of the regression models; such residuals represent time-independent measures of net diversification rate and morphological rate of evolution. We examined the correlation with pGLS, as in our other correlation analyses. Similar to our main correlation analyses, which did not account for time explicitly, time-independent rates were uncorrelated (r=0.035; P=0.825).

Further information on research design is available in theNature Portfolio Reporting Summary linked to this article.

View original post here:

The radiation continuum and the evolution of frog diversity - Nature.com