## Abstract

Interactions among traits that build a complex structure may be represented as genetic covariation and correlation. Genetic correlations may act as constraints, deflecting the evolutionary response from the direction of natural selection. We investigated the relative importance of drift, selection, and constraints in driving skull divergence in a group of related toad species. The distributional range of these species encompasses very distinct habitats with important climatic differences and the species are primarily distinguished by differences in their skulls. Some parts of the toad skull, such as the snout, may have functional relevance in reproductive ecology, detecting water cues. Thus, we hypothesized that the species skull divergence was driven by natural selection associated with climatic variation. However, given that all species present high correlations among skull traits, our second prediction was of high constraints deflecting the response to selection. We first extracted the main morphological direction that is expected to be subjected to selection by using within- and between-species covariance matrices. We then used evolutionary regressions to investigate whether divergence along this direction is explained by climatic variation between species. We also used quantitative genetics models to test for a role of random drift versus natural selection in skull divergence and to reconstruct selection gradients along species phylogeny. Climatic variables explained high proportions of between-species variation in the most selected axis. However, most evolutionary responses were not in the direction of selection, but aligned with the direction of allometric size, the dimension of highest phenotypic variance in the ancestral population. We conclude that toad species have responded to selection related to climate in their skulls, yet high evolutionary constraints dominated species divergence and may limit species responses to future climate change.

## 1. Background

All organisms are complex systems composed of many traits that interact with each other. In such systems, investigation of whether genetic constraints impose limits to evolution is of fundamental importance [1–4]. By constraints, we mean any restrictions or limitations on the course or outcome of evolution (*sensu*; [1]). The heritable associations among traits within a complex structure are represented by the additive genetic variance–covariance matrix (the G-matrix; [5]), in which the magnitudes of genetic covariation between traits results from the linkage disequilibrium and pleiotropy of the underlying genes [6,7]. On the one hand, correlations among traits can constrain the evolutionary trajectories of population means along the adaptive landscape if the G-matrix has multivariate directions (i.e. linear combinations of traits) with low genetic heritable variance [2,8]. On the other hand, multivariate dimensions that accumulate most genetic variance embedded in the G-matrix may act as lines of least evolutionary resistance, attracting population divergence in their direction, even if selection is not pulling in the same direction [9–13].

Quantitative genetics theory has been used to model the response of populations to environmental fluctuations [14–16]. These models posit that environmental alterations may produce shifts in the adaptive landscape, as indeed has been seen in nature [17,18]. Thus, genetic constraints may facilitate, hamper, or even prevent populations from adapting to these new selective pressures. Which scenario will occur depends on the degree of alignment between the direction of selection and directions of highest genetic variance [9,14–16]. Climate change is expected to act as a new and probably strong selective pressure, given that it influences several aspects of species biology, on ecological as well as on evolutionary timescales [19–21]. The interaction between selection, climate, and constraints has mainly been investigated in an ecological framework [22–24], overlooking the potential interplay between climate and constraints on phenotypic evolution at a wider phylogenetic scale [25].

A potentially good model for testing an effect of genetic constraints on multivariate evolution is a group of species that are distributed along environmental gradients and that also have differences in several quantitative traits. Populations or related species that present a wide distributional range are predicted to have gone through adaptive processes connected to environmental variation [26], in special climatic differences [27–30]. We therefore chose to study a group of closely related toad species, the *Rhinella granulosa* group [31,32], which are distributed in distinct open-habitats in Central and South America, ranging from the warmer and more arid Caatinga and Cerrado habitats to the colder and more humid Pampas and Oriental Chaco habitats (figure 1; [32]). These tropical toad species differentiated around 12–35 Ma and are distinguished mainly by the shape of their heads and by specific skull traits [32].

Specifically, we aimed to investigate the relative contributions of natural selection and constraints in shaping skull divergence in the studied toad species and the potential interplay of those factors with climate. Given that several species of the *R. granulosa* complex are to some extent restricted to specific habitats that present diverse climatic regimes, Narvaes [33] has speculated that species differentiation was linked to quaternary climatic cycles [35] in which populations became isolated from each other and adapted to their specific habitats (morphoclimatic domains; [36]). Thus, we hypothesize that some of the skull divergence across the *R. granulosa* species complex was produced by selection associated with climate. Yet, considering that the toad species present very high skull trait correlations [37], we expect evolutionary constraints to be high, potentially biasing most of the skull divergence towards the directions of high within-species genetic variance.

## 2. Methods

### (a) Analysis workflow

In order to test our prediction of climate-associated selection acting on skull divergence, we performed a complex workflow of analysis. We first present a summary of our analysis to make the connections among all the steps clear:

(i) three-dimensional landmarks, (ii) 21 linear distances, (iii) within-group species phenotypic variance covariance matrices (P-matrices), (iv) ancestral matrices as weighted averages of species P-matrices (W-matrices), (v) covariance matrix of species means (B-matrix), (vi) C-matrix: within- and between-species difference in variation (*C* = *W*^{−1} *B*), (vii) PC1 of C-matrix: axis describing the linear combination of skull traits subjected to selection (‘most selected axis’), (viii) species means projected on the most selected axis, (ix) evolutionary regressions (ERs): how much of between-species variation in the most selected axis is explained by climatic variables, and (x) formal test that species divergence was due to selection.

### (b) Sample and linear distances

We analysed 1 034 specimens of 11 toad species belonging to the *Rhinella granulosa* group (electronic supplementary material, table S1). We scanned all the specimens with an X-ray microtomography system (SkyScan 1176, Konitch, Belgium) from the Instituto de Biociências at the Universidade de São Paulo (IB-USP, Brazil), with methods and parameters described in Simon & Marroig [38]. We placed 22 three-dimensional landmarks in the toad skulls (electronic supplementary material, table S2) using TINA manual landmarking tool software [39] in order to extract 21 linear distances that represent individual bone dimensions, thus capturing local developmental and/or functional processes. The distances were allocated into three distinct functional sets: neurocranium, snout, and suspensorium ([40]; electronic supplementary material, figure S1). We opted to use linear distances to represent skull variation instead of landmark configurations, normally used in the geometric morphometrics (GM) framework [41] because of potential problems in estimating species P-matrices with the later method and issues related to constructing selection gradients. When the variation across homologous landmarks is not isotropic (i.e. homogeneous), the superimposition procedure used in the generalized Procrustes analysis (GPA) spreads the variation of more variable landmarks to the rest of the landmarks, therefore confounding variation across the structure [42,43]. Even though Procrustes distances has been shown to be a better method to estimate shape explicitly and to detect shape differences between samples than linear distances [44,45], the adequacy of using GPA to estimate covariance matrices is still not clear and solutions proposed have not been widely used [42,46]. In addition, because the interpretation of selection gradients obtained with GM is disputed [47,48] and because the evaluation of constraints in a quantitative genetics framework demands the interpretation of selection gradients [49], our choice of method is better suited to the objectives that we delineated here.

### (c) Evolutionary regression of morphology on climate

We performed ERs [50] to test our first prediction that divergence in the mean skull morphology of species is caused by selection associated with variation in mean climatic variables drawn from species geographical distributions. ER is an extension of Hansen's [51] model of evolution by stabilizing selection to continuous predictor variables that are evolving randomly through time (i.e. Ornstein–Uhlenbeck model; [52,53]). The analysis uses two parameters: adaptation half-life (*t*_{1/2}), the time it takes for an adaptation to a new selective regime to overcome the influence of ancestral states (i.e. the strength of phylogenetic inertia), and *v _{y}*, the expected variance between species that evolved for a period under the same selective regime [54]. Given that not all variation in the skulls of species was necessarily produced by selection, we first calculated the linear combination of skull traits that is expected to be subjected to selection, using the following equation [55,56]
2.1where

**W**^{−1}is the inverse of the pooled within-group covariance matrix of the basal most node in the toad phylogeny and

*is the covariance matrix of the species multivariate means. Both the W-matrix and the B-matrix were divided by their trace (total variation). The first eigenvector of the C-matrix indicates a multivariate direction of most difference between the W-matrix and the B-matrix [55,56]; that is, the combination of skull traits that differ the most when comparing within-species variation and variation among species means. We refer in the text to the first eigenvector of the C-matrix as the ‘most selected axis’.*

**B**We estimated the ancestral W-matrix as an average of all species phenotypic matrices (P-matrices) deriving from the basal node, weighted by species sample sizes and phylogeny structure ([57–60]; function PhyloW from the ‘evolqg’ R package, [61]). We controlled for skull trait mean differences owing to geography and sex within each species to construct P-matrices (using residuals of linear models; see electronic supplementary material, table S1 for factors controlled on each species). The inversion of the W-matrix as done in equation (2.1) is problematic because it is estimated with error owing to sampling, accumulating noise in the smallest eigenvalues [62]. Given that the smallest eigenvalues dominate the inverted matrix, causing biases in the estimation of *C*, we performed a noise control approach delineated in Marroig *et al*. [62]. By analysing the distribution of the W-matrix eigenvalues, we determined the last reliable eigenvalue (minimum value of the second derivative distribution) as the seventh one and subsequent eigenvalues (8–21) were replaced by this value (using the function ExtendMatrix of the ‘evolqg’ R package [61]).

We used the species multivariate means projected on this most selected axis as the dependent variable in the ER analysis and we predict that variation among species in this axis will be explained by variation in climatic variables. We extracted species climatic data from the WorldClim database that refers to the period between the years 1950 and 2000 [63]. We used all the localities in which each species were reported to occur [33], and we extracted bioclimatic variables using DIVA-GIS software. Given that several bioclimatic variables are highly correlated (see electronic supplementary material, table S4), we chose to use as predictors in the ER only variables that presented the highest contrast in the first two PCs extracted from a climatic correlation matrix (electronic supplementary material, table S5) and that did not have a correlation higher than 0.7. The bioclimatic variables tested were mean diurnal range (°C), temperature seasonality (°C), mean temperature of warmest quarter (°C), annual precipitation (mm), and precipitation seasonality (mm). Because we only have 11 species, we calculated the regression parameters for weak and strong phylogenetic inertia (*t*_{1/2} reaching 10% or 100% of the full phylogeny length, respectively) by using a recently published molecular Bayesian phylogeny of the toad species [34]. The weak and strong models were compared by AICc values only if the ER were significant for both models (i.e. 95% CI for the slope did not contain the zero value). We did not interpret all the parameter values drawn from the preferred models because we have a small species sample size [60]. Thus, we only used the intercept and slope estimates of the ER and the percentage of variance explained. ER was performed using the function oubm.fit from the ‘slouch’ R package.

### (d) Random drift tests

Even though we extracted an axis of most difference between the ancestral W-matrix and the B-matrix, we cannot be sure that these two matrices are different owing to selection or to other evolutionary processes, such as random drift. Thus, to formally test that selection is underlying species divergence, we used the analyses grounded on quantitative genetics theory to investigate a role of genetic drift versus natural selection in skull divergence [5,57–59,64,65]. The disparity in phenotypic multivariate means, as a result of random drift, is a function of the ancestral G-matrix, the time since divergence and effective population sizes of descendant populations [5]. The key implication is that variation between descendant populations (i.e. variation between species multivariate means) is expected to be proportional to the within-group variation of the ancestral population. Once again, we used weighted averages of species P-matrices to estimate ancestral matrices along the toad phylogeny. We considered our procedure adequate to estimate the ancestral matrices, because we previously compared all species P-matrices and they are very similar (ranging from 0.88 to 0.97 using Random Skewers; [37]). Furthermore, the regression drift test has nominal type I error rates when matrix correlation between the G-matrix and the P-matrix are higher than 0.7 [65]. The W-matrices are used as proxies of the ancestral G-matrices [57–59].

To compare the within- and between-population variation, we used the regression test developed by Ackermann & Cheverud [59] which is based on principal components (PCs; function TreeDriftTest; [61]). Within-population variation is represented by the variance explained by each of W-matrices' PCs (eigenvalues), whereas between-population variation (B-matrix) corresponds to the variance of the phylogenetic independent contrasts (PIC; [66,67]), calculated using species/ancestors multivariate means, projected on W-matrices' PCs (yielding species scores). By using PIC, we acknowledge that not all the descendant species originated from the same ancestor [67] and that skull trait means of closely related species may be more similar to each other than to distantly related species [66]. We estimated ancestral trait means using the maximum-likelihood (ML) approach developed by Schluter *et al*. [68] that assumes a Brownian motion (BW) model of evolution in which changes in the means do not follow any tendency (function fastAnc of ‘phytools’ R package; [69]). We have also estimated the ancestral means using linear parsimony (using Mesquite; [70,71]) and the results are similar to ML (mean vector correlation among ancestral means estimated by both methods equals 0.99). The regression slope is expected to be one if W-matrix eigenvalues and B-matrix (in log scale) are proportional (see §1.1 in the electronic supplementary material for mathematical details). To investigate deviations from proportionality, we compared the empirical regression slopes with 95% confidence intervals obtained by simulating random drift, while accounting for the phylogeny (details in §1.1 of the electronic supplementary material). We rejected drift whenever the empirical slope was below or above the lower and upper limits of the drift intervals, respectively. We performed the regression test only for five nodes in the phylogeny (black nodes in figure 3), because we did not test nodes that had less than four descendant populations (grey nodes in figure 3) owing to potentially low statistical power. We also checked that differences in P-matrices across species did not bias the drift regressions tests (see §1.2 in electronic supplementary material, SI and figure S2).

In addition to the regression test, we also computed Pearson product–moment correlations between species scores on W-matrix eigenvectors [59,60,64], also using PIC. These PCs are by definition orthogonal to each other at the within-group level, and therefore independent. Under random drift, we expect the species divergence in the PCs scores to remain independent. However, under the alternative hypothesis of natural selection producing population divergence, species PC scores could be correlated if selective pressures on these dimensions were also correlated (selective covariance; *sensu* Felsenstein [64]). We correlated *n* − 1 number of PCs, with *n* corresponding to the number of taxa in each of the five nodes tested. We performed univariate tests of PC pairwise correlations, also computing type I error rates for these tests (see electronic supplementary material, table S3). We rejected drift when at least one PC correlation was significantly different from zero at *α* = 0.05.

### (e) Selection and evolutionary constraints

We tested our second prediction of an effect of constraints on skull divergence by comparing the direction accumulating most phenotypic variance within the ancestral W-matrices (PC1 of each W-matrix) with both the direction of selection (** β** vectors) and the direction of the evolutionary response (

**Δ**vectors). The direction accumulating most variance (PC1) is called

*z**p*

_{max}[5], and it is expected to present the highest constraining effect. Any selection gradient non-orthogonal to

*p*

_{max}will produce a response biased in this direction [72]. A good estimation of constraints, therefore, depends on the quality of the ancestral matrix estimation. Considering that ancestral matrices were estimated using the species P-matrices (average number of individuals = 94) and these present very high repeatabilities (ranging from 0.96 to 0.99; [37]), we argue that constraints are reasonably well estimated in our study. All the toad species present very high trait correlations resulting in a high percentage of the total variance explained by the PC1 (ranging from 50% to 80%; [37]). If directions accumulating phenotypic variance do indeed influence evolution, then we would expect high vector correlations between PC1 and

**Δz**, indicating that divergence is aligned with constraints.

The evolutionary response is a vector of differences in skull traits between descendant – ancestral pairs ([species means directly derived from the node] − [ancestral trait means from the specific node]). Selection gradient reconstruction was done by rearranging the multivariate breeder's equation [5,11]
2.2where *W*^{−1} is the inverse of the pooled within-group matrix for each node in the phylogeny (figure 3*a*). We once again used the noise control approach on the W-matrices to avoid biasing the estimates of ** β**. Marroig

*et al*. [62] showed that selection gradient reconstruction is improved by this procedure using simulations. We reconstructed 20

**vectors in total for the whole phylogeny, using each node ancestral W-matrix (nodes 12–21 in figure 3**

*β**a*) and its corresponding Δ

*zs*. To test if the evolutionary response is aligned with the direction of selection, we computed the vectors correlations between each

**Δz**with its corresponding

**vector for each branch in the phylogeny [11]. To test if selection is aligned with constraints, we computed vector correlations between PC1 and**

*β***vectors. Finally, we evaluated the alignment of**

*β***vectors with the most selected axis to inspect the contribution of individual reconstructed selection vectors to the divergence between W-matrix and B-matrix variation. All vector correlations were done with normalized vectors (each vector element divided by the vector length), so that only vector directions are compared. The significance of the vector correlations was determined by constructing a random distribution of 1 000 21-element vector correlations drawn from a random multivariate distribution of zero mean and s.d. = 1.0. 95% of all values ranged between −0.45 and 0.45, and any value outside this range can be considered significant at**

*β**α*= 0.05.

## 3. Results

### (a) Between-species climatic variation explains skull morphological variation

Between-species variation in temperature and precipitation seasonality (figure 3*a*,*b*), as well as in mean temperature of the warmest quarter, explains 90%, 79.5%, and 73.5% of variation on species means projected on the most selected axis, respectively (electronic supplementary material, table S7). Although these climatic variables do not have the highest correlations in the climatic correlation matrix (see electronic supplementary material, table S5), they are moderately correlated. Thus, part of the morphological variance explained by each climatic variable is shared with the other two variables. The most selected axis represents a contrast between prenasal, nasal, and frontoparietal bones (distances in red in figure 2) against distances from maxillary, squamosal, and pterygoid bones (distances in blue in figure 2; electronic supplementary material, table S5), which belong to the suspensorium functional unit. Species that have longer prenasal, nasal, and frontoparietal bones and shorter maxillary, squamosal, and pterygoid bones are subjected to less temperature seasonality and more precipitation seasonality than species with the opposite morphological pattern.

### (b) Natural selection signature on toad skulls

We rejected drift in all nodes of the phylogeny that were tested, three of them with the regression test and four of them with the PC correlation test (figure 3*a* and electronic supplementary material, table S3). The highest regression slope was found for the basal most node (figure 3*b*). PIC scores on morphological PC5 are correlated with PIC scores on PC1 and PC2 (figure 3*c*) for the basal most node. PC5 is also correlated with PC3 for nodes 13 and 16, being the highest correlation for this last node (electronic supplementary material, table S3). Type I error rates are acceptable for both the regression and the PC correlation tests (see electronic supplementary material, table S3). Interestingly, the most selected axis, when extracted using PIC (i.e. B-matrix is a covariance matrix of empirical PIC), has a vector correlation of −0.65 (*p* = 0.03) with PC5, the morphological axis with significant PC correlations.

Coefficients for the first five morphological PCs are shown in electronic supplementary material, table S6. PC1 vectors from the ancestral matrices may be considered allometric size vectors (all of its coefficients have the same sign, indicating that skull traits all increase or decrease together along its direction), although presenting very high correlations with an isometric vector in log scale (0.96–0.99). PC2 is a contrast between distances from the parasphenoid bones and all other distances. PC3 opposes distances between the nasal bones and frontoparietal, maxillary, pterygoid, and mandible bones. PC5 is a contrast between frontoparietal and squamosal bones against nasal and maxillary bones.

### (c) Interplay between selection and constraints on skull divergence

Several of the evolutionary responses are aligned with the direction of the reconstructed selection gradients along the toad phylogeny (12 out of 20 phylogeny branches, **Δz** × ** β** mean vector correlation = 0.58 ± 0.21), yet most of them are also aligned with PC1 (

**Δz**× PC1 mean vector correlation = 0.75 ± 0.27; table 1). However, no selection gradient is aligned with the PC1 direction (

**× PC1 mean vector correlation = 0.06 ± 0.05). Interestingly, only when**

*β***is orthogonal to PC1 (vector correlation less than 0.02) is the evolutionary response highly aligned with**

*β***(**

*β***Δz**×

**vector correlation ≥ 0.88; table 1). Several selection gradients are aligned with the most selected axis (11 out of 20 branches;**

*β***× most selected axis mean vector correlation = 0.54 ± 0.32; table 1).**

*β*## 4. Discussion

### (a) Climatic variables as selective agents acting on toad skull variation

In this study, we integrated quantitative genetics and phylogenetic methods to investigate divergence in skull morphology and its potential relationship with distinct evolutionary processes as well as with climatic variation in a group of related toad species. The ER results suggest that climatic pressures if not directly were at least related to the selective agents that pulled a set of skull traits towards different optima in the different species. The rejection of random drift in all nodes tested supports that natural selection has acted directly on specific skull traits or as a co-selection of different sets of traits, mostly on the snout region. Narvaes [33] hypothesized that these toad species had an ancestral population widely distributed in South American open-habitats that got split owing to the expansions and retractions of forest habitats during the Quaternary climatic cycles. This hypothesis implies that species differentiation was facilitated by gene flow interruption and subsequent divergence among lineages by accumulation of genetic differences owing to drift/mutation and selection. This scenario may seem incompatible with our results favouring selection as a driver of skull divergence among species. However, rejecting drift at a particular node in the phylogeny indicates the action of selection, but does not give any insight of whether speciation was linked to selection. Thus, initial divergence on the history of the *R. granulosa* species complex could have occurred by neutral processes and afterwards each lineage could start to respond to divergent selective pressures specific to their habitats [73]. Even though we cannot be certain that climate-related selection was relevant to speciation events of the *R. granulosa* complex, given that the climatic data used in the ERs refer to very recent time periods (1950–2000), high vector correlations between the selection gradients with the most selected axis suggest that climatic selective pressures were relevant in the deep divergence of these toad species.

An association between skull and climatic variation has also been found in other natural systems, such as rodents [74–76], lizards [30], and monkeys [27,77]. In general, the explanation for skull–climate associations are based on the effect that divergent climate has on food type and availability. Thus, climate is interpreted as having an indirect effect on skull morphology through feeding mechanics. However, most anuran amphibians present very similar diets, especially Bufonidae species, which are considered ant feeders [78]. There is little information on diet preferences for the species composing the *R. granulosa* complex, yet both *R. dorbignyi* and *R. granulosa* eat ants and beetles [79,80], despite being subjected to very distinct precipitation and temperature patterns. Thus, we argue that diet probably does not account for the differences found in skull morphology of the toad species. Instead, we suggest that these skull differences are related to the reproductive ecology of the species, especially, because precipitation seasonality is known to be very relevant for anuran life-history traits [81–83]. The species from the *R. granulosa* group are explosive breeders dependent on ephemeral pools to reproduce [32], and the ability to detect sites for reproduction is probably fundamental to population persistence, especially in more seasonal habitats. The part of the toad skull more relevant to perceive water cues is the snout region that protects and supports the olfactory capsules [40]. In *Rhinella arenarum*, Jungblut *et al*. [84] showed that a new epithelium emerges after metamorphosis and it can detect water-borne cues. Thus, we may speculate that this epithelium has a functional importance to the ability of the adult toads to reach water sources for reproduction or even for maintaining water balance in the dry season. In more seasonal and dry habitats, longer snouts may provide better protection for this olfactory epithelium or it may be more developed. Further studies with anatomical comparisons of the snout region of different toad species may support this speculation.

### (b) Size variation as the major constraint underlying skull divergence

Even though the toad species responded to some extent to climate-associated selection on their skulls (as variation in the most selected axis and on morphological PC5 shows), practically, all the divergence across species was governed by evolutionary constraints as we have anticipated (table 1). Given that the reconstructed selection gradients are not aligned with PC1, the size direction that concentrates most phenotypic variance, we may infer that almost all evolutionary responses were severely deflected to the allometric direction. Interestingly, even a very low alignment of selection with PC1 (vector correlations between 0.05 and 0.1) was enough to bias the evolutionary responses against following the direction of selection (see §1.3 in electronic supplementary material and figure S3). Variation in the most selected axis and in PC5 indicate local variation in the skull that reflect shape differences across the toad species. Hence, we may interpret that climate-associated selection was acting on skull shape (*sensu* [85,86]), but the evolutionary responses were strongly biased to the allometric size direction.

The evolutionary implications of multivariate genetic constraints expressed in the G-matrix are still underappreciated in natural systems [2], though studies that estimated G-matrices show that few dimensions accumulating most genetic variance is a common phenomenon [3,25,87]. Although there are theoretical models of genetic constraints hampering the response of populations to selection associated with environmental variation [14–16], few empirical studies have investigated the interplay between constraints and climatic variables as the agents of selection [22–24], especially for multivariate systems as we have done here. An understanding of the interaction between genetic, ecological, and evolutionary processes in the past is critical to the evaluation of the evolutionary potential of species to respond to future climatic changes. Yet, we must be aware that we used phenotypic matrices as substitutes for genetic matrices, ignoring not only the possibility of phenotypic plasticity, but also the possibility of other ecological responses to habitat change, such as niche tracking. Thus, toad species may be able to compensate for future climate change with physiological, behavioural, or even morphological plasticity. Distinguishing between genetic and plastic responses in relation to climatic variation is an important issue that can only be dealt with G-matrix estimation and common garden experiments. Still, divergence in size is a common pattern of several amphibian species [88–90], suggesting that variation within species might have also acted as a constraint in other systems. Even though body size has been shown to be an important fitness component of amphibians [91,92], the study of within-species variation is essential to test whether evolution of size was indeed adaptive or the consequence of genetic constraints [72].

## Data accessibility

The datasets supporting this article were uploaded in Dryad: http://dx.doi.org/10.5061/dryad.28b06 [93].

## Authors' contributions

M.N.S. and G.M. designed the study. M.N.S. collected the data. All authors analysed and wrote the manuscript.

## Competing interests

We declare we have no competing interests.

## Funding

This study was supported by Fundação de Amparo à Pesquisa do Estado de São Paulo (FAPESP), PhD scholarship 2011/07584-2 to MNS, PhD scholarship 2011/21674-4 to F.A.M., multi-user project 2009/54203-4 and thematic project 2011/14295-7 to G.M.

## Acknowledgements

We thank all museum curators that loaned specimens for this study: David Kizirian and Robert Pascocello (AMNH); Célio Haddad and Nadya Pupin (UNESP-Rio Claro); Gustavo H.C. Vieira (CHUFPB); Julián Faivovitch and Santiago J. Nenda (MACN); Márcia F. Renner (MCN-FZB-RS); Glaucia M.F. Pontes (MCT-PUCRS); José P. Pombal Jr and Pedro Henrique Pinna (MNRJ); Taran Grant and Carolina Mello (MZUSP); Jeremy F. Jacobs and Robert Wilson (USNM); and David Cannatella and Travis LaDuc (TNHC). M.N.S. also thanks Patricia Narvaes for helping with specimen identification and Miguel Trefaut Rodrigues for helping contacting museum curators from Brazil. M.N.S. also thanks David Houle for suggesting the eigenanalysis used to extract the most selected axis. We also thank all colleagues from our laboratory for discussing this work and the reviewers for their comments and suggestions to improve the manuscript.

## Footnotes

Electronic supplementary material is available online at https://dx.doi.org/10.6084/m9.figshare.c.3497571.

- Received August 12, 2016.
- Accepted September 29, 2016.

- © 2016 The Author(s)

Published by the Royal Society. All rights reserved.