Habitat use affects morphological diversification in dragon lizards

Habitat use may lead to variation in diversity among evolutionary lineages because habitats differ in the variety of ways they allow for species to make a living. Here, we show that structural habitats contribute to differential diversification of limb and body form in dragon lizards (Agamidae). Based on phylogenetic analysis and ancestral state reconstructions for 90 species, we find that multiple lineages have independently adopted each of four habitat use types: rock-dwelling, terrestriality, semi-arboreality and arboreality. Given these reconstructions, we fit models of evolution to species’ morphological trait values and find that rock-dwelling and arboreality limit diversification relative to terrestriality and semi-arboreality. Models preferred by Akaike information criterion infer slower rates of size and shape evolution in lineages inferred to occupy rocks and trees, and model-averaged rate estimates are slowest for these habitat types. These results suggest that ground-dwelling facilitates ecomorphological differentiation and that use of trees or rocks impedes diversification.


Introduction
One of the great questions in evolutionary biology concerns the causes of differences in diversity among clades. Ecological factors are often implicated to explain this pattern because the ecological circumstances available to the members of a lineage contribute to the mode of natural selection they experience and thus shape ecological divergence, morphological adaptation and the evolution of new species. Although much work has focused on the role of biotic interactions within communities (e.g. competition-driven divergent selection in Anolis lizards [Williams, 1972;Losos, 2009], Hawaiian silverswords [Carlquist et al., 2003] and Darwin's finches [Schluter, 1988;Grant & Grant, 2008]), other aspects of a lineage's ecology may also be important for diversifica-tion. In this study, we test the hypothesis that diversity varies as a function of habitat.
For many reasons, some habitats may foster greater diversity than others. Some habitat types may be readily subdivided, perhaps because of spatial complexity (e.g. coral reefs [Alfaro et al., 2007]) or geographical area (e.g. arid habitats in Australia [Rabosky et al., 2007]) and may thereby present evolutionary lineages with many alternative means for microhabitat specialization or local adaptation. Other habitats may impose stringent functional constraints that lead to strong selection resisting ecological and phenotypic divergence away from an adaptive peak (Butler & King, 2004;Collar et al., 2009). Habitat types may also contribute differently to diversification because they vary in the number and type of species interactions they present, such as the presence or absence of predators (McPeek & Brown, 2000). In addition, some habitats may provide opportunities if they are variable across space in the strength of species interactions (McPeek, 1996) or in their functional demands.
The consequences of habitat use for diversification have been investigated primarily in the context of explaining variation in species richness among clades. Evolutionary transitions between habitat types that differ in the opportunities they provide for ecological divergence are often implicated to explain shifts in rates of lineage diversification. For example, several invasions of coral reefs in tetraodontiform fish lineages are temporally coincident with increases in rates of diversification (Alfaro et al., 2007), the transition into arid habitats is associated with elevated rates of diversification in Australian skinks (Rabosky et al., 2007), and diversification rates in damselflies vary across gradients of pond permanence (McPeek & Brown, 2000).
One explanation for these associations between habitat use and species richness is that the process of lineage splitting is mechanistically linked to niche differentiation. However, species richness and ecological diversity need not be correlated during evolution (Foote, 1993;Losos & Miles, 2002;Adams et al., 2009), and elevated rates of lineage diversification within a habitat type do not require increases in rates of ecological evolution. Indeed, the neutral theory of biodiversity emphasizes the extent to which species diversification may occur in the absence of ecological differentiation (Hubbell, 2001).
Using morphological variation in ecologically relevant characters as a surrogate for ecological variation, we asked whether a relationship exists between habitat use and ecological diversity. We focused on dragon lizards (Agamidae), an ecologically and morphologically diverse radiation of iguanian lizards comprising roughly 400 species distributed throughout the Old World. Agamid lizards vary in their structural habitat use, including species that primarily use rocks, trees or terrestrial surfaces as well as some semi-arboreal species that frequently use both trees and terrestrial surfaces. Because the ability to move about and hold position in the environment is partly a consequence of structural habitat use and because movement is important to the performance of ecological tasks, such as foraging, evading predation and defending territory (Losos, 1990;Irschick & Garland, 2001), these four types of habitat use may contribute differently to ecomorphological diversification in agamid lineages.
To evaluate this hypothesis, we applied a phylogenetic approach that tests for associations between habitat and rates of morphological evolution in agamid lineages. We inferred phylogenetic relationships for 90 agamid species based on mitochondrial DNA sequences, reconstructed ancestral habitat use and used these reconstructions as the basis for fitting models of evolution to species values for morphological traits. We then compared fit and parameter estimates for models that differ in the number of evolutionary rates (based on the Brownian motion model), where rates are allowed to vary among lineages inferred to use different habitat types. We interpreted habitat types associated with high rates of morphological evolution to be those that facilitate ecological divergence.

Reconstructing phylogeny
We reconstructed phylogenetic relationships for 90 agamid species-representing nearly one-quarter of the group's recognized species diversity-as well as four outgroup species. Our molecular data set included 1.2 kb of mitochondrial DNA including partial sequences for the protein-coding genes, NADH dehydrogenase subunit 1 (ND1) and cytochrome c oxidase subunit I (COI), and the complete sequence for NADH dehydrogenase subunit 2 (ND2). This analysis excluded intervening tRNA-coding regions because they are highly variable among the sampled taxa, making unambiguous alignment of these regions difficult and potentially unreliable (Schulte et al., 2004a;Schulte & de Queiroz, 2008). All sequences were obtained from GenBank (accession numbers are in Table S1) and aligned by eye. Base positions inferred to have ambiguous homology at the ends of ND1 and ND2 were excluded from phylogenetic analyses (198 of 1281 aligned positions). Alignment is available in TREE REEBASE (Study accession number S2669, Matrix accession numbers M5148; to be added upon acceptance of manuscript).
We used these sequences to simultaneously infer phylogenetic relationships among agamid species and estimate branch lengths in relative time using Bayesian phylogenetic analysis and a relaxed molecular clock approach implemented in the program BEAST BEAST (Drummond et al., 2006;Drummond & Rambaut, 2007). We partitioned mtDNA sequences by codon position and, for each partition, separately fit a general time reversible model of nucleotide substitution that allows for gammadistributed substitution rate variation among sites and invariant sites (Yang, 1994) because previous analysis of these sequences for a subset of the agamid species included in this study showed that this model provided the best fit relative to simpler substitution models (Schulte et al., 2004a). Variation in substitution rates among lineages was modelled by a lognormal distribution in which the mean rate was set to 1.0 (i.e. no external calibration was used to estimate divergence times), and no correlation was assumed between ancestor and descendant branches (Drummond et al., 2006;Drummond & Rambaut, 2007). Uninformative priors were applied for all parameter estimates.
We used BEAST BEAST to sample the posterior probability distribution of phylogenetic trees and substitution model parameters given species' sequence data according to a Markov chain Monte Carlo (MCMC) algorithm (Drummond & Rambaut, 2007), which we ran twice for 25 · 10 6 generations per run. For each run, a random starting tree was generated under a Yule (pure-birth) process (Drummond & Rambaut, 2007), and the first 2.5 · 10 6 generations were discarded as burn-in. We verified the adequacy of sampling from the posterior probability distribution using the program TRACER TRACER (Drummond & Rambaut, 2007) to determine that effective sample sizes for model parameter estimates were greater than 200 (Drummond et al., 2006) and confirmed convergence of the two MCMC runs using the program AWTY AWTY (Nylander et al., 2008) by inspecting the correlation between split frequencies.
The central question of this study focuses on the role of structural habitat use in morphological evolution. Although reconstruction of agamid phylogeny is necessary to address this question, strong inference about a single, best phylogeny is not. Rather than base our analyses on a single consensus tree, we retained a set of 1000 phylogenies sampled from the posterior distribution for use in ancestral state reconstructions (see next section). To do this, we randomly sub-sampled 500 trees from each of the MCMC runs. Because MCMC algorithms of BEAST BEAST sample trees in proportion to their posterior probability, performing subsequent analyses on this set of trees incorporates uncertainty in agamid phylogeny into our analyses in a manner similar to the method of Huelsenbeck & Rannala (2003).

Reconstructing ancestral habitat use
To reconstruct ancestral habitat use in agamid lineages, we used stochastic character mapping, which is a Bayesian method that implements MCMC to sample character reconstructions in proportion to their posterior probability under a Markov process of state transitions given species' character states and a phylogeny (Nielsen, 2002;Bollback, 2006). We assigned structural habitat types-rock-dwelling, arboreal, terrestrial and semi-arboreal-to the 90 agamid species included in our phylogenetic analysis based on Stuart-Fox & Owens's (2003, their online appendix 3) and Stuart-Fox & Ord's (2004, their online appendix A) syntheses of ecological data in Agamidae. Our data set excluded species that these studies categorized as generalists-those that are reported to occur on rocks, trees and terrestrial surfaces-because this category appeared to contain a heterogeneous set of species, including species comprised of generalized individuals, polymorphic populations or populations that vary in habitat use. We used the program SIMMAP SIMMAP (Bollback, 2006) to generate 10 stochastic maps of habitat use for each of the 1000 phylogenies sampled from our phylogenetic analysis (described earlier). From the resulting 10 000 habitat use reconstructions, we randomly sampled 500 for use in model-fitting analyses. The set of 500 reconstructions thus represents sampling from the posterior probability distributions of both trees and ancestral reconstructions. Because the MCMC algorithms of BEAST BEAST (Drummond & Rambaut, 2007) and SIMMAP SIMMAP (Bollback, 2006) sample trees and character state histories, respectively, in proportion to their posterior probabilities, use of this set of structural habitat reconstructions in subsequent model-fitting analyses allowed us to integrate over uncertainty in both phylogeny and ancestral states.

Quantifying species' morphology
We quantified species values for morphological features of the body and limbs that have functional consequences for iguanian lizard locomotor performance (e.g. Reilly & Delancey, 1997;Spezzano & Jayne, 2004), including snout-vent length, tail length, body width, pectoral width, pelvic width, humerus length, ulna length, carpal length, IV metacarpal length, femur length, tibia length, tarsal length and IV metatarsal length. Descriptions of landmarks used to delimit these morphological variables are detailed in Schulte et al. (2004b). Species values are means of measurements made on preserved specimens of adult males and females (species data are available from the authors by request). Sample sizes within species ranged between one and 36 individuals (median = four individuals; Table S1).
To account for correlations between variables and to reduce the dimensionality of the morphological data set, we performed principal components analysis (PCA) on the correlation matrix of agamid species values. Species scores on principal component (PC) axes were then used as species character values in subsequent model-fitting analysis. Additionally, we quantified sampling error for species PC scores. We used the eigenvalues and eigenvectors from PCA on species values to transform morphological values for all individuals into PC scores and estimated the pooled within-species variance for each PC. Each species' sampling variance was then taken as the pooled within-species variance divided by the number of individuals sampled for that species.

Fitting models of morphological evolution
To assess the effects of structural habitats on morphological diversification in Agamidae, we fit several models of evolution to species' PC scores and reconstructions of agamid phylogeny and ancestral habitat states. We examined Brownian motion models of character evolution that differed in the number of rates of evolution, defined as the time-independent variance parameter, r 2 , of the Brownian motion model of character evolution (Felsenstein, 1985(Felsenstein, , 1988Garland, 1992;Martins, 1994;Collar et al., 2005;O'Meara et al., 2006;Thomas et al., 2006). These models specified separate evolutionary rates for lineages inferred to use different habitat types following O' Meara et al. (2006) and Collar et al. (2009). Inferred habitat states in agamid lineages were based on ancestral reconstructions from SIMMAP SIMMAP (Bollback, 2006). The most complex model specifies separate rates for rockdwelling, arboreal, semi-arboreal and terrestrial lineages (4-rate full: r 2 rock , r 2 arbor , r 2 semi , r 2 terr ), whereas the simplest model specifies a single rate for all agamid lineages, regardless of inferred habitat state (1-rate no effect: r 2 rock¼arbor¼semi¼terr ). We also explored the fit of eight additional models in which the effects of some habitat types were assumed to be equal. We note that the following models are not an exhaustive set of all possible models given these four habitat categories, but rather a subset of models that we deemed most plausible based on hypothesized shared and unique properties of the different surface types. We fit three three-rate models that set rate categories to be equivalent according to possible shared effects of habitats: the effects of arboreality and semi-arboreality may be the same because these habitat types both require that species occur in forests and move along and cling to trees, which may result in similar ecological and functional demands (3-rate shared tree effect: r 2 rock , r 2 arbor¼semi , r 2 terr ); alternatively, the effects of terrestrial and semi-arboreal habitats may be equivalent if the demands of ground surfaces prevail over those of trees in semi-arboreal species (3-rate shared ground effect: r 2 rock , r 2 arbor , r 2 semi¼terr ); additionally, because both rocks and trees present steeply inclined surfaces to the species that use them, rock-dwelling and arboreality may have similar effects in agamid lineages (3-rate shared incline effect: r 2 rock¼arbor , r 2 semi ,r 2 terr ). We also fit four two-rate models to test for unique effects of each habitat type: the unique effect of rock-dwelling (2-rate rock effect: r 2 rock , r 2 arbor¼semi¼terr ); of arboreality (2-rate arboreal effect: r 2 arbor , r 2 rock¼semi¼terr ); of semi-arboreality (2-rate semiarboreal effect: r 2 semi , r 2 rock¼arbor¼terr ) and of terrestriality (2-rate terrestrial effect: r 2 terr , r 2 rock¼arbor¼semi ). Finally, because species that predominantly use steeply inclined surfaces (rock-dwelling and arboreal species) may experience demands that differ from species that use the ground frequently (terrestrial and semi-arboreal), we fit a fifth two-rate model that separates the effects of steeply inclined surfaces from those of ground-dwelling (2-rate incline-ground: r 2 rock¼arbor , r 2 semi¼terr ). To evaluate the effects of different habitat states on morphological evolution, we fit multiple-rate Brownian motion models rather than multiple-peak Ornstein-Uhlenbeck (OU) models, which describe evolution under selection towards fixed phenotypic optima (Felsenstein, 1988;Hansen & Martins, 1996;Hansen, 1997;Butler & King, 2004). This decision was based on two considerations. First, multiple-rate Brownian models allowed us to assess whether the process of diversification differed across lineages characterized by different evolutionary regimes (O'Meara et al., 2006;Thomas et al., 2006;Collar et al., 2009). In contrast, current multiple-peak OU models allow inferences to be made about the positions of phenotypic optima corresponding to different selective regimes, but they assume that the process of evolving towards those values (i.e. the strength of selection and the magnitude of the rate of stochastic change) is the same for all selective regimes (Hansen, 1997;Butler & King, 2004;Scales et al., 2009). Because our goal was to test the hypothesis that habitat types contribute differently to the process of evolution, we were more interested in evaluating habitat-associated variation in the rate of morphological change than in finding fixed adaptive morphologies. Second, Brownian motion has been shown to adequately describe adaptive evolution under a variety of evolutionary conditions (e.g. when environmental change causes shifts in the position of adaptive peaks, or when unconsidered lineage-specific differences in selection, environment or genetics have large effect relative to the strength of selection because of the considered selective regime [Hansen & Martins, 1996;Hansen, 1997]). Nevertheless, we note that both Brownian motion and OU processes are relatively simple models used to describe complex reality, and the multiple-rate models we chose to fit were those that we deemed most appropriate for evaluating our hypothesis given the aforementioned considerations.
We used the computer program Brownie 2.1 (O'Meara et al., 2006;O'Meara, 2008) to fit each of the 10 models to species scores for each PC across the 500 habitat reconstructions. The method uses maximum likelihood and extends the noncensored approach of O' Meara et al. (2006) to accommodate models that specify multiple evolutionary rates for phylogenetic branches associated with states of a categorical variable (in this case, habitat use). The method also incorporates sampling error for species values using the approach of Martins & Hansen (1997). For each of the four PCs, we fit each model to the set of species' scores iterating over the 500 habitat reconstructions. This process resulted in distributions of model parameter estimates and fit scores for each combination of model and PC, and the spread of the distributions represent variation that is because of uncertainty in phylogeny and ancestral habitat reconstructions.
To assess model fit, we used the Akaike information criterion (AIC), which is a function of the likelihood, L, of the data given the model and the number of parameters, k: AIC = 2k ) 2ln (L). More specifically, we used AICc, which is AIC with a correction for small sample size-appropriate when the number of observations, in this case species, is < 40 times the number of estimated parameters (Burnham & Anderson, 2002). Lower AICc scores indicate better fit. To select the best fitting model for each PC given uncertainty in phylogeny and ancestral habitats, we evaluated the average model fit over the 500 reconstructions and compared mean AICc across models. We note that averaging AICc in this way is valid because the data for species are the same for all iterations of model fitting. The 500 reconstructions should not be considered as different, independent data sets but as alternative estimates of phylogeny and ancestral habitat states sampled in proportion to their posterior probabilities given the same data, and averaging AICc over this sample allowed us to quantify model fit in a way that integrates over uncertainty in these estimates. We also calculated Akaike weights from mean AICc as an additional descriptor of each model's fit to each PC. Akaike weights describe the proportion of support a model receives relative to the total support available for all models (Burnham & Anderson, 2002).
In addition to comparing mean AICc, we also explored the sensitivity of model selection to alternative phylogeny and habitat reconstructions. We compared AICc among models on each reconstruction and generated distributions for DAICc-the difference between each model's AICc and the best fitting model's AICc.
We were unable to unambiguously select a single best fitting model for any of the PCs because several models received substantial support; DAICc was < 2 (Burnham & Anderson, 2002) and several models provided the best fit for a similar proportion of the reconstructions. To assess the effect of habitat on morphological evolution in light of this uncertainty in model selection, we compared model-averaged estimates of the rate of evolution associated with each habitat use type, where the modelaveraged rate is the average rate across all models weighted by each model's Akaike weight (Burnham & Anderson, 2002). This weighting strategy ensures that parameter estimates from better supported models count more towards the overall model-averaged rate estimates. The resulting model-averaged rates of PC evolution for each habitat state thus average across uncertainty in the model of character evolution as well as uncertainty in the reconstruction of agamid phylogeny and ancestral habitat states.

Results
Bayesian phylogenetic analysis yielded a set of 1000 ultrametric trees that are generally consistent with previous phylogenetic hypotheses involving agamid lizards (Macey et al., 2000;McGuire & Kiew, 2001;Melville et al., 2001;Schulte et al., 2004a;Hugall et al., 2008). Figure 1 shows the topology and branch lengths in relative time for the maximum clade credibility tree (i.e., the tree with the greatest posterior probability summed over all nodes) from our sample of the posterior distribution; we summarize node support on this tree by identifying nodes whose posterior probabilities are less than 0.95. We found strong support for the monophyly of three recognized biogeographical groups, a clade of southwest Asian and African species, a clade of southeast Asian agamids, and a clade of species from Australia and New Guinea which is the sister group to the southeast Asian species, Physignathus concincinus.
Stochastic mapping of habitat use across the set of agamid phylogenies provides strong support for multiple, independent transitions to each of the four structural habitat types. All sampled reconstructions infer multiple origins of semi-arboreality (median = 7, minimum = 5; maximum = 10), more than 99% of reconstructions infer multiple origins of rock-dwelling (median = 4, maximum = 7 origins) and terrestriality (median = 5, maximum = 10 origins), and 91.3% reconstruct more than one origin of arboreality (median = 3, maximum = 7 origins). Figure 1 depicts one stochastic habitat reconstruction from SIMMAP SIMMAP (Bollback, 2006) on the maximum clade credibility tree and shows the median number of origins of each habitat type.
PCA on species values for morphological traits provides four axes that together account for 96.3% of the total variation between species. PC 1 accounts for 82.9% of the variation and loads strongly and similarly across all variables (Table 1); we interpreted PC 1 to represent variation among species that is because of differences in size. The three subsequent PCs collectively explain 78.0% of the variation in shape (i.e., the variation not explained by PC 1). Loadings for PCs on the original morphological variables are reported in Table 1 and the distribution of species on PCs 2 and 3 is shown in Fig. 2. Notably, PC 2 separates arboreal species from species that use the other habitat types; all arboreal species have negative PC 2 scores, reflecting generally long, narrow bodies and relatively long forelimbs, whereas nearly all other species have positive scores on this axis, indicative of shorter, wider bodies and relatively short forelimbs (Fig. 2). Habitat groups do not appear to separate as clearly on the other shape axes (PC 3 and 4) or on the size axis (PC 1).
Morphological PCs are generally best fit by multiplerate models that infer evolutionary rates to be slower in rock-dwelling and arboreal lineages than in terrestrial and semi-arboreal lineages. Table 2 presents parameter estimates and fit scores (-ln likelihood and AICc) for each model as means and standard errors taken across 500 habitat reconstructions. Also shown in Table 2 are DAICc and Akaike weights (calculated from the mean AICc scores), which served as the basis for choosing the preferred model for each PC. We also compared model fit on each of the 500 habitat reconstructions, and Table 3 reports each model's mid-95% density interval for DAICc as well as the percent of reconstructions for which each model is preferred (DAICc = 0.0) and the percent for which each is relatively unsupported (DAICc > 2.0).
In general, models that allow habitat-associated rate variation are more strongly supported than the singlerate model. The best fitting multiple-rate model is much more strongly supported on average than the single-rate model for PCs 1, 2 and 3; for PC 4, several multiple-rate models provide better fit than the single-rate model, but the single-rate model receives substantial support (Table 2). Looking across reconstructions, we found that for all PCs the single-rate model is preferred in 5% or fewer of the reconstructions. In addition, the single-rate model is unsupported in the vast majority of reconstructions for PCs 1, 2 and 3 (100%, 70% and 96%, respectively; also see Fig. S1 for histograms of the single-rate model's DAICc for each PC). The superior fit of the multiple-rate models over the single-rate model supports the hypothesis that rates of morphological evolution vary in association with habitat.
Size evolution in agamids is best fit by the four-rate model, which infers an exceptionally high rate associated with semi-arboreality (r 2 semi = 146.07 ± 46.02), a slower but intermediate rate for terrestriality (r 2 terr = 20.33 ± 4.67), and yet slower rates for arboreality (r 2 arbor = 13.22 ± 2.15) and rock-dwelling (r 2 rock = 2.24 ± 2.06).  1 Maximum clade credibility phylogeny for 90 agamid species illustrating a single reconstruction of structural habitat use. Nodes are supported by > 0.95 Bayesian posterior probabilities unless otherwise noted, and branch lengths are proportional to time (i.e. root node depth is 1.0). Colour/shade of branches indicates inferred habitat use based on stochastic character mapping (see key). Habitat states for species are given by colour/shade of terminal nodes. Transitions between habitat states are highlighted by vertical, black bars. We used this reconstruction of phylogeny and habitat state and 499 others as the basis for fitting models of evolution in which rates of morphological diversification are allowed to differ in lineages that use different habitat types.
This model is strongly preferred on average over seven of the other models (DAICc > 4; Table 2) and provides the best fit for 77% of the habitat reconstructions (also see Fig. S2 for histograms of DAICc for the models preferred for each PC). However, the four-rate model is only somewhat preferred over the two-rate semi-arboreal effect model (r 2 semi = 166.73 ± 51.75, r 2 rock¼arbor¼terr = 14.13 ± 0.86, DAICc = 1.42), which provides the best fit for 18% of reconstructions. In addition, there is some support for the three-rate shared incline effect model (r 2 rock¼arbor = 11.59 ± 1.96, r 2 semi = 148.40 ± 48.17, r 2 terr = 20.29 ± 4.65, DAICc = 1.94), but this model provides the best fit for only 2% of reconstructions. Each of the three models receiving support infers an elevated rate of evolution in semi-arboreal lineages, and the two best fitting models estimate a separate, intermediate rate for terrestriality. Evolution of PC 2 is best described by the two-rate terrestrial effect model, which infers a four-fold higher evolutionary rate associated with terrestriality (r 2 terr = 1.63 ± 0.39) compared to the rate shared by other habitat types (r 2 rock¼arbor¼semi = 0.40 ± 0.12). This model provides the best fit in 73% of habitat reconstructions (also see Fig. S2); however, on average it is only somewhat preferred over the three-rate shared tree effect model (r 2 rock = 0.08 ± 0.14, r 2 arbor¼semi = 0.45 ± 0.16, r 2 terr = 1.58 ± 0.40, DAICc = 0.99), which is preferred in 20% of the reconstructions. Other models receiving support are the two-rate incline-ground model (r 2 rock¼arbor = 0.40 ± 0.16, r 2 semi¼terr = 1.51 ± 0.42, DAICc = 1.23), and the three-rate shared incline effect model (r 2 rock¼arbor = 0.42 ± 0.15, r 2 semi = 0.32 ± 0.44, r 2 terr = 1.61 ± 0.42, DAICc = 1.71), though these models provide the best fit for only 3% of the reconstructions. Nevertheless, models that receive at least moderate support are similar in that they infer rates of PC 2 evolution to be higher in terrestrial lineages than in rock-dwelling and arboreal lineages.
The two-rate incline-ground model best fits the evolution of PC 3 and estimates a shared evolutionary rate for rock-dwelling and arboreal lineages (r 2 rock¼arbor = 0.33 ± 0.16) that is nearly six times slower than the rate shared between semi-arboreal and terrestrial lineages (r 2 semi¼terr = 1.84 ± 0.43). On average this model is only weakly or moderately supported over the three-rate models (shared tree effect: r 2 rock = 0.03 ± 0.09, r 2 arbor¼semi = 0.41 ± 0.16, r 2 terr = 1.87 ± 0.41, DAICc = 0.39; shared ground effect: r 2 rock = 0.03 ± 0.07, r 2 arbor = 0.37 ± 0.20, r 2 semi¼terr = 1.79 ± 0.43, DAICc = 0.11; shared incline effect: r 2 rock¼arbor = 0.34 ± 0.17, r 2 semi = 0.88 ± 0.88, r 2 terr = 1.96 ± 0.48, DAICc = 0.99), the two-rate terrestrial effect model (r 2 terr = 1.91 ± 0.42, r 2 rock¼arbor¼semi = 0.37 ± 0.13, DAICc = 0.55), and the four-rate model (r 2 rock = 0.03 ± 0.08, r 2 arbor = 0.38 ± 0.20, r 2 semi = 0.82 ± 0.87, r 2 terr = 1.91 ± 0.48, DAICc = 1.12). The two-rate incline-ground model was the most commonly preferred model but is the best fit for only 47% of habitat reconstructions (see Fig. S2). Other models that provide the best fit for a substantial proportion of reconstructions include the three-rate shared tree effect model, the tworate terrestrial effect model and the two-rate rock effect model, whose fit varied widely among reconstructions (Table 3). Despite the variability in model fit across reconstructions, the models that consistently receive support share similarities in their parameter estimates; the evolutionary rate associated with terrestriality is greater than the rates associated with arboreality and rock-dwelling.  Fig. 2 Scatterplots of agamid species in a morphospace defined by principal components 2 and 3. Color-coding for species' habitat states is the same as in Fig. 1. Loadings of original variables on PCs are described for each axis. For brevity, we use 'hindlimb length' and 'forelimb length' to describe loadings on PCs when more than one element of that limb loads strongly on that axis. See Table 2 for details about PCA.

Habitat use and diversification in agamids 1039
On average, PC 4 is best fit by the three-rate shared ground effect model, which infers the shared evolutionary rate for semi-arboreality and terrestriality (r 2 semi¼terr = 0.52 ± 0.06) to be much higher than the rate for arboreality (r 2 arbor = 0.23 ± 0.02) and rock-dwelling (r 2 rock = 0.01 ± 0.01). However, support is quite evenly distributed among the 10 models examined (Tables 2 and 3). In fact, five models receive only somewhat less support than the preferred model, including the three-rate shared tree effect (r 2 rock = 0.01 ± 0.01, r 2 arbor¼semi = 0.25 ± 0.02, r 2 terr = 0.52 ± 0.05, DAICc = 0.26), two-rate rock effect (r 2 rock = 0.01 ± 0.01, r 2 arbor¼semi¼terr = 0.33 ± 0.02, DAICc = 0.23), two-rate terrestrial effect (r 2 terr = 0.50 ± 0.05, r 2 rock¼arbor¼semi = 0.22 ± 0.02, DAICc = 0.66), two-rate incline-ground (r 2 rock¼arbor = 0.21 ± 0.02, r 2 semi¼terr = 0.50 ± 0.06, DAICc = 0.23) and single-rate (r 2 rock¼arbor¼semi¼terr = 0.29 ± 0.01, DAICc = 1.02) models. Moreover, although it is preferred on average, the three-rate shared ground effect model is not the most commonly preferred among reconstructions; it is best fit for 26% whereas the tworate rock effect model is best fit for 32%. Two additional models-the three-rate shared tree effect and two-rate incline-ground models-also provide the best fit to a substantial proportion of reconstructions (Table 3). Although the single-rate model receives at least moderate support for a large proportion of reconstructions, multiple-rate models consistently provide better fit to the evolution of PC 4 (Table 3).
In spite of the ambiguity in selecting the absolute best fitting multiple-rate model, the better fitting models are those that allow rates of evolution to be faster in lineages inferred to be terrestrial or semi-arboreal. This pattern is captured in comparisons of the model-averaged estimates of the rates of PC evolution for each habitat type, which account for uncertainty in model selection. Comparisons of these rate estimates reveal a consistent pattern across PCs that describe limb and body shape variation (PCs 2, 3 and 4); terrestrial lineages have experienced the fastest rates of evolution, semi-arboreality is associated with intermediate evolutionary rates, and arboreal and rockdwelling lineages evolve at similarly slow rates, though rates are slower for rock-dwelling (Fig. 3). The pattern is somewhat different for model-averaged rates of size (PC 1) evolution; semi-arboreal lineages exhibit a much higher rate than the other three habitat types, though terrestriality is associated with a somewhat higher rate than arboreality and rock-dwelling (Fig. 3).

Discussion
Habitat use has had strong effects on the evolution of limb and body form in agamid lizards, suggesting that habitats contribute differently to ecological diversification. Models that allow the rate of morphological evolution to vary between lineages that use different habitats provide better fit to the distribution of species' trait values than the single-rate model, which constrains the rate to be the same in all agamid lineages (Tables 2  and 3). Although we found ambiguity in selection of the preferred multiple-rate model, the better fitting models are generally similar in that they infer slower rates of morphological evolution in rock-dwelling and arboreal lineages than in terrestrial or semi-arboreal lineages (Table 2). Moreover, model-averaged estimates of the rates of evolution for PCs that describe morphological shape (PCs 2-4) reveal a clear and consistent pattern of rate variation associated with habitat states: terrestrial lineages evolve fastest, semi-arboreal lineages evolve at an intermediate rate, and arboreal and rock-dwelling lineages experience similarly slow rates (Fig. 3). For PC 1, which describes variation in size, semi-arboreality is associated with the highest rate, though arboreality and rock-dwelling again exhibit the slowest rates of size evolution (Fig. 3). These results suggest that terrestrial habitats facilitate microhabitat differentiation or evolution along additional ecological axes, whereas the use of trees or rocky surfaces impedes such diversification.

Diversification within habitat categories
The effects of structural habitat use on diversification occur across multiple lineages that have independently derived these habitat use types. Nearly all stochastic character maps sampled from the posterior distribution infer multiple gains of each habitat use type, and given this set of habitat reconstructions, the best fitting models are those in which rates of morphological evolution vary with habitat (Table 3). Moreover, transitions to each habitat type are inferred to have occurred independently in clades that represent agamid radiations in different geographical regions (see Fig. 1). Multiple transitions across phylogenetically and geographically distant lineages allow the opportunity to evaluate the generality of the effects of habitat use on diversification and to detect the possible influence of lineage-specific effects unrelated to habitat (Read & Nee, 1995). To this end, in the following paragraphs, we review the major groups included in each of the habitat categories and qualitatively compare habitat-associated rate estimates in the major agamid clades (see Fig. 4).
Terrestrial species are spread across three phylogenetically distant agamid groups. The Australian agamid radiation is renowned for its ecomorphological diversity (Pianka, 1986;Melville et al., 2001Melville et al., , 2006Harmon et al., 2003), and most of the 70 recognized species within the clade are terrestrial, including forms as diverse as Moloch, Pogona (bearded dragons), the spindly legged Diporiphora, Tympanocryptis and a variety of Ctenophorus. Indeed, the rate of shape evolution is inferred to be higher in the Australian radiation than in the other major continental radiations (Fig. 4). In addition, southeast Asian terrestrial agamids include not only the species rich Calotes, but also some highly disparate taxa from the Indian subcontinent Table 2 Parameter estimates and model fitting for multiple-rate Brownian models in Agamidae, where evolutionary rates are associated with inferred habitat use states.
Values are means ± standard errors from model fitting across 500 reconstructions of habitat use on agamid phylogenies sampled from posterior probability distribution. Bold denotes the best fitting model.   , 1994]). The terrestrial lineages of the southeast Asian clade exhibit a high rate of size evolution relative to other terrestrial agamid lineages (Fig. 4). The rate of shape evolution in this clade, however, is somewhat lower than the estimate across all terrestrial agamids, though shape evolves more rapidly in terrestrial than in arboreal southeast Asian agamids (Fig. 4). By contrast, the southwest Asian ⁄ African group including Trapelus and Pseudotrapelus contains lesser variation and has experienced substantially lower rates of shape evolution than terrestrial species in the other continental radiations, though the rate of shape evolution in these lineages is somewhat elevated relative to southwest Asian ⁄ African rock-dwellers (by about a factor of two; see Fig. 4). We note, however, that Trapelus and another terrestrial lineage within this clade, Phrynocephalus, are not well sampled in our study. Also, by one account Pseudotrapelus siniatus is considered to be rockdwelling, rather than terrestrial (El Din, 2006), which would lessen our sample for estimating the rate of evolution associated with terrestriality in this clade.
Relatively few semi-arboreal species are included in the data set, none forming large clades. However, many of these taxa are quite distinct, including the two extremely large and semi-aquatic Physignathus species, which turn out not to be closely related ( Fig. 1; Macey et al., 2000;Hugall et al., 2008). In addition, Chlamydosaurus (the frilled lizard) and the long-headed, longlegged Lophognathus are closely related, but highly morphologically disparate. The high rate of size evolution in semi-arboreal agamids seems to be largely driven by a high rate in several semi-arboreal Australian lineages, though the elevated rate of shape evolution for semiarboreality relative to arboreal and rock-dwelling  Fig. 3 Model-averaged estimates for the rates of PC evolution in rock-dwelling (grey), arboreal (green), semi-arboreal (blue) and terrestrial (yellow) lineages. Point estimates are means of rate estimates from the 10 multiple-rate models and have been weighted by Akaike weights. Error bars are standard errors, representing uncertainty in habitat and phylogenetic reconstructions. Note that the y-axis for PC 1, an axis of size variation, is different from the y-axis for PCs 2, 3 and 4, which describe shape.
agamids is largely because of a high rate inferred for multiple semi-arboreal southeast Asian lineages (Fig. 4). Within southeast Asia, shape seems to evolve at a faster rate in semi-arboreal lineages than in terrestrial lineages, though the rate associated with semi-arboreality is highly variable across habitat reconstructions. In addition, within Australia, the rate of shape evolution is estimated to be lower in semi-arboreal lineages than in the arboreal clade, Hypsilurus, though again there is substantial uncertainty in the rate estimate for semi-arboreal lineages.
Arboreal agamids have generally experienced relatively slow morphological diversification. Low rates are inferred for two large clades, Draco and Hypsilurus as well as a third paraphyletic group of arboreal species from southeast Asia, including Gonocephalus, Lyriocephalus and Ceratophora (Fig. 4). The estimates for the rates of size and shape evolution in the latter group are somewhat higher than the rates in Draco or Hypsilurus, but are still lower than the rates in southeast Asian lineages that use other habitat types.
Rock-dwelling agamids have experienced the lowest rates of morphological evolution. In contrast to the high rates in lineages of the Australian radiation that use other habitats, two clades of rock-dwelling Ctenophorus species within this larger group have experienced very slow rates of size and shape evolution (Fig. 4). Rock-dwelling lineages of the southwest Asian ⁄ African clade are comprised mostly of Laudakia species, and the rates of evolution of size and shape in these lineages are nearly as low as in the rock-dwelling Ctenophorus lineages (Fig. 4). African rock-dwelling Agama are represented by only one species in this study; however, most species in this genus appear morphologically homogeneous, in agreement with trends exhibited by the other rock-dwellers.
For the most part, habitat has consistent effects on diversification in the three major continental radiations of agamids. Although the magnitude of influence of some habitat types varies in the three major continental radiations, these clade and regional effects do not confound our general conclusions about the effects of habitat use. Within the major clades, rates of size and shape evolution are slower in rock-dwelling and arboreal lineages than in lineages that use ground surfaces. Terrestrial lineages of the Australian radiation, which occur primarily in deserts, diversified in shape more rapidly than desert rock-dwellers or forest species. Among southeast Asian taxa-all forest-dwelling-two groups of arboreal species have diversified more slowly than related lineages that occupy other habitat types. Rock-dwelling lineages from southwest Asia ⁄ Africa have diversified slowly, but we also note that under-sampling of terrestrial southwest Asian ⁄ African species prevents us from ruling out a generally slow rate of morphological evolution in this clade. Finally, we note one exception to the general trend is the somewhat slower rate of shape evolution in semi-arboreal Australian lineages relative to the Australian arboreal clade, Hypsilurus. Although rates of size evolution are very rapid in semi-arboreal Australian agamids, rates of shape change in these lineages may be just as slow or even slower than in Australian rockdwelling and arboreal lineages. However, the rate estimate for semi-arboreality in this clade varies substantially across phylogeny and ancestral reconstructions (see Fig. 4), and thus the rank of this rate estimate relative to arboreal and rock-dwelling lineages is uncertain.

Mechanisms by which habitat affects diversification
A variety of factors could result in differences in rates of evolution in different structural habitats. We discuss several of these factors in the following paragraphs. We Fig. 4 Rates of size and shape evolution for each habitat category estimated within the major agamid clades. Point estimates for the rates of size evolution are means of estimates for the rates of PC 1 evolution across the 500 habitat reconstructions, and estimates for the rates of shape evolution are sums of the mean rate estimates for PCs 2, 3 and 4. Error bars are standard errors for the rate estimates across the 500 reconstructions. Shapes correspond to clade identity (see Fig. 1): diamonds are southwest Asian ⁄ African lineages, triangles are southeast Asian lineages, and squares are Australian lineages. Horizontal lines represent rate estimates for all Agamidae (based on the full, four-rate Brownian model), and grey boxes represent ± one standard error taken across reconstructions.
Habitat use and diversification in agamids 1045 note, however, that at this point we have insufficient data to distinguish among possibilities.

Functional constraints
Some structural habitats may impose stricter functional constraints than others, leading to stronger selection resisting diversification away from morphologies that confer effective use of that surface (Simpson, 1953;Butler & King, 2004). In particular, movement and position holding on the steeply inclined surfaces that rocks and trees present to the species that use them may require morphological features and performance abilities that prevent falling (Sinervo & Losos, 1991;Revell et al., 2007;Goodman et al., 2008). In contrast, terrestrial habitats are generally broad and flat, and clinging and climbing are functional considerations that generally apply to a much lesser extent to ground-dwellers. Terrestrial surfaces may therefore impose weaker functional constraints on form and allow for morphology to diverge by neutral evolution or adaptation for other functions (Hansen, 1997;Alfaro et al., 2005). Indeed, ground-dwelling seems to permit a broader variety of forms; an extreme example of this is the thorny devil, Moloch horridus, which moves slowly over sandy deserts  and possesses short limbs and a wide body that is unique among terrestrial agamid species (Fig. 2a, M. horridus has the highest score on PC 2 and the most negative score on PC 3).

Habitat complexity
The converse of functional constraints, some structural habitats may provide more ways of making a living. Terrestrial habitats, for example, provide opportunities for burrowers, species that live in leaf litter, in grass, that run quickly in open areas and that remain cryptic against the soil. Terrestrial members of the Australian Ctenophorus radiation, for example, have evolved different refugeseeking strategies, including digging burrows in sand or loose soils or hiding in areas covered by shrubs or grasses (Melville et al., 2001). Burrowing and manoeuvering through structurally complex habitats may impose additional selective demands on locomotor performance that contribute to morphological diversification among terrestrial lineages (Thompson & Withers, 2005). In this way, terrestrial habitats may provide for finer microhabitat differentiation than rocky or arboreal habitats.
In contrast, it is conceivable that fewer ways exist to adapt to rocky or arboreal habitats. Recent work demonstrating convergent and parallel evolution of morphology and performance in rock-dwellers from several different radiations support the hypothesis that there are few ways to make use of rocky habitat (Revell et al., 2007;Goodman et al., 2008). However, ecomorphological diversification in arboreal habitats is well documented in some lizard groups, such as anoles (Williams, 1972;Losos, 2009), geckoes Vitt et al., 2003) and chameleons (Bickel & Losos, 2002), and these groups challenge the generality of our finding that arboreal habitats are diversity limiting. Unlike these groups, agamids lack specialized toe-pads or foot structures, which confer exceptional clinging abilities in the species that possess them (Irschick et al., 1996;Zani, 2000). Relatively limited clinging capabilities may prevent agamids from diversifying to make use of the different structural niches utilized by these other lizard clades.

Ecological interactions
If the number of co-occurring, related species is greater in some habitat types than in others, then selection for resource partitioning may lead to adaptive divergence. On the other hand, some habitats may have more competing species of other taxa-such as insectivorous birds or mammals-which may limit the ability of agamids to diversify. Certainly, many terrestrial Australian agamid species occur in sympatry in some areas of the Australian desert, perhaps accounting for diversification in this clade. On the other hand, as many as seven species of the arboreal southeast Asian clade, Draco, occur sympatrically (Inger, 1983;McGuire & Dudley, 2005), yet rates of evolution in this clade are low ( Fig. 4; though we note that sympatric Draco differ in wing size (McGuire, 2003), a morphological attribute which we did not measure). Asian rainforests may have more competing species of other taxa (e.g. birds, small mammals) than Australian deserts-which are known to be dominated by squamates (Pianka, 1986)-and this difference could explain these discrepant patterns.

Geographical distribution of species
The converse of large numbers of sympatric species, some habitat types may not facilitate co-occurrence of ecologically similar species. Rather, species may replace each other across the geographical landscape, and thus may occupy the same or similar niches, only in different places. The southwest Asian ⁄ African clade is one example in which sympatry of clade members is generally quite low, and the slow rate of morphological evolution in rock-dwelling members of this clade may be a consequence of similar selection pressures acting on members of this clade.

Caveats
Comparisons of model fit and model-averaged rate estimates provide evidence that morphological diversification varies as a function of habitat use in agamid lineages. We note, however, that the model-fitting approach employed in this study is limited to detecting the best of the evolutionary models we specified to evaluate our hypothesis (Burnham & Anderson, 2002). Therefore, the scope of our conclusions is limited to the relative fit of these models. Although the superior fit of the multiple-rate models over the single-rate models supports the hypothesis of habitat-associated variation in morphological diversification, we cannot exclude the possibility that an alternative, unspecified model provides a better fit to agamid species values and phylogeny. Likewise, our analysis does not rule out roles for other factors that may have influenced morphological evolution in agamid lineages. For example, diversification may also vary with differences in intrinsic factors, such as genetic constraints or origins of novel structures.
A related point concerns the susceptibility of our approach to lineage-specific factors unrelated to habitat that might speed or slow morphological evolution within clades (Read & Nee, 1995). Our conclusions are somewhat protected against such confounding factors because multiple transitions into each of the four habitat use categories are likely to have occurred; however, large clades characterized by a single habitat type (e.g. Draco) could exert undue influence on the estimated rate of evolution associated with that habitat. In such a case, we would not be able to disentangle the role of habitat from any other derived factor shared within that clade. Indeed, for this reason we qualitatively assessed possible lineagespecific effects in Fig. 4, but we note that a more general test of habitat's effects on morphological evolution would involve fitting separate habitat-associated rates within major clades or comparing evolutionary rates between samples of sister clades that differ in habitat use.
In addition, our sampling of agamid species likely influenced the pattern of evolutionary rate variation we document. Our data set represents approximately onequarter of recognized agamid species diversity, though this proportion is not uniform across the major clades. The under-sampling of species from the southwest Asian ⁄ African clade may be the most problematic with respect to our conclusions because at least two primarily terrestrial genera, Trapelus and Phrynocephalus, may exhibit relatively little morphological disparity. If this is indeed true, their under-representation may have resulted in an inflated estimate for the rate of morphological evolution across terrestrial agamids. However, based on a more extensive morphological data set (J. Schulte, unpublished), disparity in Trapelus and Phrynocephalus is similar to that of other terrestrial clades that we sampled more extensively in this analysis (e.g. Tympanocryptis [size and shape], Diporiphora [size]; see Fig. S3). Consequently, rates of evolution in these clades may be comparable to those in other terrestrial clades, and thus their under-representation in our data set may not have biased our results.

Conclusions
Our results provide compelling evidence that habitat use shapes diversification of limb and body form in Agamidae. The pattern of variation in rates of morphological evolution suggests that terrestrial habitats promote ecological differentiation whereas diversification is slower in rocky and arboreal habitats; however, the precise mechanism by which these habitats contribute differently to diversification remains speculative. We recommend that future research investigate the extent to which functional constraints, habitat complexity or biotic interactions within habitats influence the pattern we document. For example, locomotor performance tests could be applied across a range of agamid forms to determine whether arboreal and rocky surfaces impose more stringent functional demands than terrestrial surfaces. Also, comparisons of the number of co-occurring agamid species within each habitat type could assess whether these habitats present different numbers and types of species interactions. Application of such studies to Agamidae or other animal taxa has the potential to provide further insights into how habitat contributes to differential diversification among evolutionary lineages.

Supporting information
Additional Supporting Information may be found in the online version of this article: Table S1 Agamid species' GenBank numbers and morphological specimen data. Figure S1 Distributions of DAICc for the single-rate model for each PC. Figure S2 Distributions of DAICc for the model that on average best fits each PC. Figure S3 Disparity in several groups of primarily terrestrial agamids.
As a service to our authors and readers, this journal provides supporting information supplied by the authors. Such materials are peer-reviewed and may be reorganized for online delivery, but are not copy-edited or typeset. Technical support issues arising from supporting information (other than missing files) should be addressed to the authors.