Environmental correlates of phenotypic evolution in ecologically diverse Liolaemus lizards

Abstract Evolutionary correlations between phenotypic and environmental traits characterize adaptive radiations. However, the lizard genus Liolaemus, one of the most ecologically diverse terrestrial vertebrate radiations on earth, has so far shown limited or mixed evidence of adaptive diversification in phenotype. Restricted use of comprehensive environmental data, incomplete taxonomic representation and not considering phylogenetic uncertainty may have led to contradictory evidence. We compiled a 26‐taxon dataset for the Liolaemus gracilis species group, representing much of the ecological diversity represented within Liolaemus and used environmental data to characterize how environments occupied by species' relate to phenotypic evolution. Our analyses, explicitly accounting for phylogenetic uncertainty, suggest diversification in phenotypic traits toward the present, with body shape evolution rapidly evolving in this group. Body shape evolution correlates with the occupation of different structural habitats indicated by vegetation axes suggesting species have adapted for maximal locomotory performance in these habitats. Our results also imply that the effects of phylogenetic uncertainty and model misspecification may be more extensive on univariate, relative to multivariate analyses of evolutionary correlations, which is an important consideration in analyzing data from rapidly radiating adaptive radiations.

the association of adaptive variation to speciation rates (Marques et al., 2019). Phenotype-environment associations, or covariation between species or population phenotypes and environmental variables, often result from the evolution of divergent functional capabilities for optimal performance in specific ecological contexts (i.e., Arnold, 1983). Therefore, revealing the evolutionary association between environmental and phenotypic traits (i.e., Irschick et al., 2008) is important for understanding the process of adaptive radiation (Gavrilets, 2009;Glor, 2010;Harder & Schluter, 2001).
Nevertheless, to date, because of the ease of obtaining data, most studies have largely been limited to studying climatic variables (namely average temperatures and precipitation) and their impacts on body size evolution. However, the climate is only one aspect of an organism's functional environment which may have limited impacts on organism performance, and body size represents a limited perspective on phenotypic evolution (Wainwright, 2007). Furthermore, various aspects of an organism's environment may impose selection in different ways on alternate elements of the phenotype (Irschick et al., 2008). Such variation likely has led to inconsistency in the results of studies that have assessed phenotypic-environment correlations in various groups. On top of largely consisting of tests of evolution restricted to univariate dependent variables (namely body size) relative to climate, many macroevolutionary studies have also had limited exploration of the impacts of phylogenetic uncertainty. This is problematic given the amount of research indicating that rapidly evolving clades also result in a large amount of discordance in phylogenetic estimation and hybridization-leading to uncertainty in reconstructing topological relationships and timing of divergences (Ortego & Knowles, 2022;Rangel et al., 2015;Talavera et al., 2013).
Squamate reptiles display a startling array of phenotypic and ecological diversity, but with inconsistent results with respect to the impact of the environment on phenotypic evolution depending on the phenotypic dependent variable and the environmental correlate tested. Multiple squamate families display strong evolutionary relationships between body shape (i.e., limb and body dimensions) and structural habitat characteristics (Losos, 1990;Melville & Swain, 2000;Tulli et al., 2009;Wiens et al., 2006), providing evidence for the evolution of optimal locomotor performance for different microhabitat characteristics and not climate-driven evolution. Alternatively, body size variation, an important ecomorphological axis of variation in snakes, is not correlated with climate in Lampropeltine snakes (Pyron & Burbrink, 2009) and in general squamates show inconsistent patterns of body size evolution in response to the occupation of different climatic regimes (Slavenko et al., 2019). It may be that body size evolution is not that prominent in many squamate clades, and that shape variation is a more important axis of phenotypic evolution that is predicted by structural habitat characteristics, rather than climate. However, shape variation in invasive Anolis lizards is not consistent with structural habitat differences, and rapid adaptation, between native and invasive populations (Kolbe et al., 2007) indicating that environment does not always lead to rapid morphological change even in adaptive radiations. Thus, the exploration of evolutionary relationships between multiple phenotypic and environmental axes is essential to our understanding of how adaptive evolution progresses and what factors are the most important in aiding the diversification of rapidly radiating groups.
The lizard genus Liolaemus occupies one of the widest climate ranges observed in a vertebrate clade, extending from sea level on the southern Pacific and Atlantic coasts, through the Atacama Desert, Validivian Forest, cold Patagonian steppes, to the high Andes (~5000 m elevation; (Cei, 1986;Donoso-Barros, 1966;Espinoza et al., 2004;Schulte et al., 2000). Given that the 260 currently described species of Liolaemus (Uetz, 2019) cover such ecological breadth, hypotheses abound regarding the evolutionary mechanisms by which its species diversity arises. Speciation within Liolaemus has been characterized by a slow, but consistent, accumulation of lineages with significant diversification in body size (Harmon et al., 2003; but see Olave et al., 2020). Body size has also been shown to covary with ecological, life history, and other morphological traits, like ovipary/viviparity and altitude in Liolaemus (Meiri, 2008;Pincheira-Donoso et al., 2011). Harmon et al. (2003) suggested that convergent body size evolution among sub-clades had occurred in response to the occupation of different ecophysiological environments across elevational gradients. Espinoza et al. (2004) further showed that smaller animals resided in cooler climates, and hypothesized that the evolution of herbivory and thermal constraints limited the evolution of body size. Body size was shown to increase with latitude, relating to thermal amplitude and temperature, in the Liolaemus geotschi species complex, supporting the notion that variation in body size relates to the thermal biology of Liolaemus species in different climatic regimes (Azócar et al., 2013).
The occupation of cooler climates has also been associated with the evolution of viviparity in Liolaemus (Esquerré et al., 2019;Pincheira-Donoso et al., 2013).
In contrast, traits like body shape and limb proportions may not evolve in response to habitat diversification within Liolaemus (Jaksić et al., 1980;Schulte et al., 2004;Tulli et al., 2012). However, this is not consistent and studies have also largely focused on broad evolutionary patterns across a sparsely sampled Liolaemus phylogeny. Such limited taxon sampling may impede meaningful analyses of evolutionary patterns (e.g., Heath et al., 2008), particularly the detection of evolutionary correlations between ecology and morphology (e.g., Ackerly, 2000). Analyses of macroevolutionary trends in Liolaemus have largely been undertaken without assessing the impacts of phylogenetic uncertainty, which is problematic in a group that has experienced extensive hybridization and rapid radiation (Esquerré et al., 2021;Olave et al., 2018). In support of environmentally driven morphological shape change, Tulli et al. (2009)) suggested that limb morphology in Liolaemus evolved for maximal performance in different habitats (Tulli et al., 2011), especially in response to substrate type (Tulli et al., 2016). Foot morphology may also be evolutionary labile in Liolaemus (Tulli et al., 2012). Narrower pelvises and longer tibias and metatarsals (i.e., foot length) have been connected with greater sprint speeds in lizards (Irschick & Jayne, 2000), especially in phrynosomatids which have similar morphologies to Liolaemus (Bergmann & Irschick, 2010). Furthermore, foot and limb morphology covaries with escape behavior in Liolaemus (Schulte et al., 2004).
We use multivariate comparative phylogenetic analyses to test the relationship between body size and shape evolution and multiple classes of environmental data in the well-sampled and ecologically diverse L. gracilis group and allied taxa. This clade includes up to 33 candidate species, including described and undescribed species. The distribution crosses the cold Patagonian steppe deserts (i.e., Liolaemus bibronii sensu stricto), mixed deciduous forests (i.e., Liolaemus vhagar), and highland steppe biomes (i.e., Liolaemus puna), and thus represents much of the environmental diversity encompassed by the genus Morando et al., 2007;Olave et al., 2011). This is the first analysis of multivariate morphological evolution patterns at this phylogenetic scale within Liolaemus, with extensive taxon sampling. We conduct all our analyses by explicitly incorporating phylogenetic uncertainty into our interpretation of results. Specifically, we assess the diversification of body size and shape, along with vegetation, temperature, precipitation, and elevation data, test if environmental and morphological evolution are evolutionarily correlated, and discuss how patterns of trait evolution within this Liolaemus clade compare with those inferred across Liolaemus as a whole, and more broadly in squamates. We also discuss the impact of model choice and phylogenetic uncertainty where evolutionary correlation assesses phenotype-environment correlations with dependent multivariate and univariate data.

| Taxon sampling
There are currently 33 species considered within the L. gracilis species complex (Quinteros et al., 2020), including 25 described and 8 candidate species. Recent taxonomic papers have addressed the phylogenetic history and taxonomic composition of this species complex, but there are caveats. Portelli and Quinteros (2018) suggested some supported groups were not congruent with those Quinteros et al. (2020) propose. Also, cases of hybridization with phylogenetically distantly related species confuse taxonomy (e.g. Liolaemus abdalai may not be part of the L. gracilis complex; Morando and Avila (2020)). Within the L. gracilis complex, there are two clades (Portelli & Quinteros, 2018), one mostly central-northern clade, the Liolaemus robertmertensi clade and a mostly southern clade, the L. bibronii clade. Our taxon sampling includes representatives of both these clades, with additional candidate species that Esquerré et al. (2021) show is monophyletic.
We sampled tissues from 47 individuals across 26 terminal taxa of the L. gracilis species complex (Table S1; Figure 1). Several candidate species, assigned to L. bibronii, including clades 1-2, 3c, 5, 9-12, and 14-16, are novel to this paper. These novel candidate species were delimited using methods outlined in Wiens and Penkrot (2002) using tree-based methods with morphological data (Martinez, 2012 Several candidate species were included in our molecular dataset for reconstructing phylogenetic relationships but were dropped from F I G U R E 1 Distribution map of the L. gracilis species complex shown overlain over a digital elevation model of Patagonia, position relative to South America is shown inset. Different species are indicated by the various colored and shaped symbols as outlined in the inset key analyses due to a lack of morphological data (Table S1), this included L. bibronii clades 6, 17, and 19 (Morando et al., 2007). Outgroup taxa, Liolaemus walkeri and Liolaemus punmahuida, were included in phylogenetic reconstruction, but not included in comparative analyses.

| Environmental data
Vegetation data are from satellite data taken from the NASA MODIS/Terra database (https://modis.gsfc.nasa.gov/data/). For measures of vegetation structure, we used mean normalized difference vegetation index (NDVI; a measure of vegetative biomass) and standard error of NDVI (describing seasonality in NDVI), both data from MODIS product MODIS13A1 (see https://lpdaac.usgs. gov/produ cts/mod13 a1v00 6/), and proportional tree cover, using MODIS product MODIS44b (see https://lpdaac.usgs.gov/produ cts/mod44 bv006/). MODIS data from the years 2002-2011 were downloaded using the ModisDownload.R source code v3.2, now part of the rts package v1.1.3 (Naimi, 2021). Datasets (HDF files) were converted to TIFF files using the MRT (Modis Reprojection Tool) Batch Reprojection java tool, MRTBatch.jar. Tree cover data are collected annually, so we were not able to calculate anything other than an average across the 10-year period using the calc function of the raster package v3.5.11 (Hijmans, 2021). For NDVI, data are collected at 16-day intervals. First, we calculated mean annual variables, then averaged these annual subsets of data across the 10-year period using the calc function, as above for tree cover, for mean NDVI. Then we calculated yearly standard deviations from the 16-day data and averaged these across the 10-year period using the calc function, to represent average seasonal fluctuation in vegetative biomass.
Vegetation variables were chosen to reflect structural vegetation features of the habitat. NDVI measures the greenness of vegetation and provides an index of vegetative productivity and biomass (Fensholt et al., 2004) and for our purposes provides a measure of the productivity of habitats-which is possibly linked to phenotypic evolution. Climate data were obtained from publicly available GIS layers from the WorldClim Database (http://www.world clim.org; Hijmans et al., 2005). Temperature variables included WorldClim layers 1-11, while precipitation variables included WorldClim layers 12-19. Climate and vegetation data were extracted for each individual sample using the spatial coordinates associated with morphological specimens. All variables were log-transformed prior to analysis.

| Phenotypic data
Phenotypic data were collected by one of the authors from specimens utilized in genetic analyses where possible (135 specimens; 5.7 ± 4.3 [mean ± standard deviation] individuals per species; Table S1). Specimens are located in the LJAMM Herpetological Only males were used for this data set to avoid complications of sexual dimorphism in characters. Measurements of snout-vent length (SVL; mm) were collected as a proxy for body size. Morphological measurements collected were axilla-groin length (AG), arm length (AL), tibia length (TbL), foot length (FL), and pelvic girdle width (PW; all in mm). Measurements were missing for L. bibronii candidate species 5 for tibia length only, therefore for completeness values were imputed using a regression approach using the relationship between body size and tibia length to predict the missing tibia value from body size information available for that candidate species.
Regression analyses were performed using the lm function of the stats package v4.0.2 (R Core Team, 2020).

| Ecological and phenotypic trait data reduction for comparative analyses
Prior to analysis, all morphological variables were standardized via log-transformation. Body shape in lizards is connected to locomotory performance (Bergmann & Irschick, 2010;Irschick & Jayne, 2000), is labile in Liolaemus (Tulli et al., 2012) and connected to escape behavior (Schulte et al., 2004). To allow for more functional analyses of the evolution of shape in relation to performance, by allowing for the inclusion of allometry (Bergmann & Irschick, 2010), we corrected for the effects of body size on shape variables by calculating the ratios of AG, AL, TbL, FL, and PW with SVL (e.g., AG:SVL).
All morphological variables displayed significant regressions against SVL validating the use of ratios (Table S3). Given the potential for differences in patterns of evolution, and in evolutionary correlations between shape, size, and environmental variables (vegetation, precipitation, temperature, and elevation), several separate datasets were analyzed. For size, shape, temperature, precipitation, and vegetation data variables were averaged for each species. To simplify environmental datasets containing potentially correlated variables, each dataset underwent variable reduction prior to further analyses using principle components analyses (PCA) implemented in the dudi.
Temperature, precipitation, and vegetation axes comprising 99% of the variation were retained for further analysis. Output from these analyses, including variable loadings, are in Table S4. Variations in datasets across the phylogeny can be found in Figures S1-S5.
Univariate datasets were size (SVL) and elevation, multivariate datasets were shape, vegetation, precipitation, and temperature.

| Time-calibrated species tree analyses
We used *BEAST v1.8.2 to estimate the species tree and lineage divergence times simultaneously (Heled & Drummond, 2010). Lineage divergence times were estimated using calibration methods described in (McCormack et al., 2011), as no internal fossil calibrations exist for the L. gracilis species group. Species tree analyses were time-calibrated using a mean rate of evolution for the cytochromeb locus (0.0193555 mutation rate/site/million years) calculated by Olave et al. (2015) using a fossil that places the divergence between Eulaemus and Liolaemus at 20 million years ago (Albino, 2008). This rate was fixed for the cytochrome-b locus only in our species tree analyses under a strict clock model (McCormack et al., 2011) using multilocus genetic data in *BEAST (Heled & Drummond, 2010) for 500 million generations, sampling every 5000 steps with a 10% burnin. Tree files for each run were combined and summarized after confirming singular run convergence (i.e., high ESS values), topological convergence across runs, and removing burnin using Tracer v1.7 (Rambaut et al., 2018). The median time-calibrated species tree was then used in all further comparative analyses, along with 1000 trees randomly sampled from the posterior distribution where appropriate.

| Patterns of ecological and phenotypic diversification
Disparity-through-time (DTT) analyses (Harmon et al., 2003) were used to assess changes in subclade disparity in phenotypic and environmental traits throughout the evolutionary history of the L. gracilis group. The disparity is calculated as the average pairwise Euclidean distance between species (MDI). To provide a temporal context to trait change, the observed disparity is calculated as the average relative disparity at each node, including all subclades present at that particular time (Harmon et al., 2003). Values of MDI closer to 1 indicate that species within subclades overlap substantially, having independently evolved to occupy similar convergent regions of trait space. MDI values that are closer to 0 suggest, alternatively, that variation is distributed among subclades (Harmon et al., 2003).
Some have suggested that this latter case, when lower than the null model of Brownian Motion (BM) evolution, represents constrained trait evolution (Smith et al., 2011). These analyses rely on comparing observed data to a BM null model on slow, constant change through time.
DTT analyses were undertaken using means of all variables for multivariate datasets and for univariate datasets. Analyses were performed using the dtt function of the geiger R package v2.0.7 (Pennell et al., 2014) with the effects of phylogenetic uncertainty mapped using the methods outlined in Edwards et al. (2015), for iterating DTT analyses over 1000 trees randomly drawn from the posterior distribution of species trees. We incorporate phylogenetic uncertainty because posterior trees can vary in not only the topological relationships inferred but also the timing of divergencies among species. Therefore, testing many trees from the posterior allows us to incorporate uncertainty in topological reconstruction and divergence timing in interpreting our results. This analysis allows for an interpretation of observed subclade disparity in light of phylogenetic uncertainty, both regarding the magnitude of observed subclade disparity relative to the null and the timing of shifts throughout history.
The magnitude of difference between observed disparity and the null model can be used to determine if trait evolution constrained (below the BM null), or if traits are have diversified at any point in history (above the null) among subclades (e.g., Smith et al., 2011). For instance, traits typically showing an early burst model of evolution will show increased disparity early in evolution followed by declining to present relative to the BM null model.

| Mode of environmental and phenotypic trait evolution
Because there could be differences in evolutionary trajectories between clades (Uyeda et al., 2018) driving evolutionary correlation, we sought to test for these differences in a phylogenetic framework using MANOVA and ANOVA analyses for multivariate and univariate data respectively. These analyses rely on deciphering the appropriate model of trait evolution first which DTT analyses suggest might not always be BM. To determine the model of trait evolution that best fit each dataset, we compared BM, Ornstein-Uhlenbeck (OU), and Early Burst (EB) models of trait evolution for both multivariate and univariate trait datasets. We assessed model fit using the mvBM, mvOU, and mvEB functions of the mvMORPH package in R v1.1.4 (Clavel et al., 2015). For these analyses, we used the "rpf" method to impose a computationally intensive generalized least squares approach; iterating over 1000 trees from the posterior distribution. For each tree, we used the AICc score weights to assess which model of trait evolution best fit the data using the aicw function of the mvMORPH package v1.1.4 (Clavel et al., 2015).
To test if there were differences in environmental and phenotypic traits that were largely driven by a northern and southern lineage split, we assessed if there were significant differences between these groups in trait values. For multivariate trait datasets, we used a phylogenetic MANOVA implemented in the mvgls function setting model to the model of best fit chosen (see above), error to "true" to allow for estimation of measurement error from the data, a method to "LOOCV" for leave one out cross-validation of the penalized loglikelihood, penalty to "RidgeArch" to linearly shrink the covariance structure to the target of "unitVariance" which is proportional diagonal variance relative to identity. Posthoc tests were undertaken using the manova.gls function with 999 permutations and the Wilks test statistic. Both mvgls and manova.gls are part of the mvMORPH package v1.1.4 (Clavel et al., 2015). Wilks statistic is equivalent to an F-statistic, which we calculated for univariate analyses, making comparisons possible. For univariate trait analyses (i.e., elevation), we undertook a phylogenetic ANOVA using gls from the package nlme v3.1-152 (Pinheiro et al., 2013) and undertook an ANOVA using the anova function of the base R package v4.0.2. All ANOVA/ MANOVA analyses were iterated over 1000 trees from the posterior distribution to account for phylogenetic uncertainty. Analyses were implemented using the model of trait evolution that best fit the data.

| The evolutionary relationships between environment and phenotype
Significant evolutionary relationships between environmental and phenotypic traits typify ecologically driven phenotypic evolution (Garland Jr. et al., 1992;Glor, 2010). To investigate how environmental variables explain variance in phenotypic trait values, we investigated evolutionary correlations using phylogenetic linear models to understand the relationships between phenotype and environment.
Evolutionary correlation tests were undertaken using the mvgls function of the mvMORPH package v.1.1.4 (Clavel et al., 2015). The evolutionary model of best fit, between the BM, OU, and EB models, was determined using mvgls function setting error to "TRUE" and using the "LL" method and iterated across 1000 randomly chosen posterior trees. The AIC values of respective models were then compared for each individual tree and the model that fit the majority of trees was selected as the best fit using the aicw function. The evolutionary correlation was then estimated using the same parameter settings model testing analyses. Univariate data were analyzed in a similar manner using the gls function of the nlme package v3.1.153 (Pinheiro et al., 2013) and using the corBrownian, corMartins, and corBlomberg correlation structures from the ape r package v5.6-1 (Paradis & Schliep, 2019) for the BM, OU, and EB models respectively. To ensure we were testing an early burst, and not a decelerating model for the corBlomberg structure, values of g were fixed at 0.5.
Both univariate and multivariate analyses were undertaken using the best-selected model of trait evolution identified for the correlation structure between compared datasets. Analyses were iterated over 1000 randomly chosen trees, randomly chosen from the posterior distribution of trees to account for phylogenetic uncertainty.
Dependent variables consisted of the PC scores of multivariate datasets (temperature, precipitaiton, and vegetation) or univariate data (elevation), analyzed separately. From this analysis of phylogenetic uncertainty, we calculated the proportion of trees supporting relationships between dependent and predictor variables, and the mean and 95% CI of test statistics and model coefficients.  (Figures 1 & 2), including L bibronii sensu stricto and many other lineages assigned to L. bibronii, termed the L. bibronii group of Quinteros et al. (2020), has a more recent common ancestor arising 1-2 million years ago, and shows a subsequent recent and rapid divergence of lineages and species throughout southern Patagonia. We refer to these clades as the northern and southern lineages respectively throughout.

| Patterns of ecological and phenotypic diversification
The impact of phylogenetic uncertainty varies among phenotypic and environmental traits (Figure 3), with phylogenetic uncertainty impacting the interpretation of if traits have higher relative disparity than the null (i.e., diversification in trait values) or if traits are neutrally evolving for all environmental data (Figure 3c-f). This is especially exaggerated for vegetation (Figure 3d) which the median tree shows fits a more classic pattern seen in adaptive radiation, with increasing disparity early in the tree followed by a decline toward the present-yet posterior tree distributions indicate that broad patterns that cannot be differentiated from the null distribution of trees (gray shading). The same can be seen for precipitation (Figure 3e)

| Mode of environmental and phenotypic trait evolution
All individual trait datasets were shown to evolve via an OU model of trait evolution, with the exception of elevation (Table 1). For the most part, model weights were clear (with values close to or at 1) among competing models suggesting distinguishability of the model of best fit. The exception was elevation where most trees from the posterior supported a BM model (average AICw 0.65; Table 1) over the OU and EB models, but AICc differences were low (~2.57) indicating that models were not as distinguishable. Phylogenetic uncertainty seemed to have little impact on our results, given the tight confidence intervals observed on the AICw scores, parameter values (EB and OU models), and the differences in AICc scores between models tested across trees (Table 1). Body size (100% of trees), body shape (100% of trees), vegetation (95.5% of trees), precipitation (84.9%), temperature (83.4%), and elevation (100% of trees) data showed little to no evidence of significant differences between the northern and southern clades (Table 2). Again, tight confidence intervals among trees on the inferred test statistic and p-value estimates indicate that phylogenetic uncertainty has little impact on results.

| The evolutionary relationships between environment and phenotype
Model of evolution tests showed that all multivariate phenotypeenvironment correlation structures best fit the BM model (Table   S5) for environmental correlations with body shape. However, the model fits for a BM model were generally marginally favored (i.e., AICc < or = 2) but consistently favored across 1000 posterior trees (i.e., >99% of trees supported a BM model). Model weights ranged from 0.42 to 0.58 here with tight confidence intervals across trees on all estimated parameters, indicating a limited effect of phylogenetic uncertainty on the interpretation of results or favored model.
For body size, the model of evolution tests showed that phenotypeenvironment correlation structures best fit an OU model (Table S5), with greater power to differentiate among models (i.e, AICc <8 and > 99% of posterior trees supporting model of best fit). Again, according to the spread of estimated parameters being narrow, there seemed to be a limited effect of phylogenetic uncertainty on inferring the model of best fit for correlations with body size.
Greater power to distinguish models in univariate tests, relative to F I G U R E 2 Median time-calibrated species tree showing the evolutionary relationships among species within the Liolaemus gracilis species complex, including a predominantly southern and northern clade. The 95% confidence intervals on divergence dates are shown in the gray node bars, and as the text above the major nodes for the northern, southern, and the most recent common ancestor of the group. A time axis is also displayed below the tree. Numbers around nodes represent the posterior probabilities for species tree support for each node multivariate analyses, may relate to the relative power needed to distinguish the best-fit model of trait evolution. Table 3 outlines the results from phenotype-environment correlations. Body shape is significantly (p < .01) correlated with vegetation PC1 which corresponds to vegetative biomass (NDVI) and seasonality in biomass. The results for body shape appear robust (i.e., low confidence limits) to phylogenetic uncertainty when estimated across 1000 posterior trees. We do not find significant, and robust, estimates of a phenotype-environment correlation between body size and any variable, although ~25% of trees favor a correlation between vegetation PC1 (mean and seasonality of NDVI), and ~42% favor a correlation with temperature (wet and cold temperature; Table 3). Confidence limits are extremely wide for the test statistic and p-value estimates with most comparisons with body size yielding a few trees that supported a particular correlation.
However, the same comparison may yield the directly opposing pattern in a set of trees. For example, p-values range from <.001 to .95 for the correlation between elevation and body size. To test if this was due to the difference in the model, between an OU and BM, we used the OU model for body shape analyses and the BM for body size analyses (Table S5). These results showed no impact of both phylogenetic uncertainty and model choice on interpretation for multivariate analyses of correlations between body shape and environmental data. However, the univariate analyses were both susceptible to phylogenetic uncertainty, with often large confidence intervals around the test statistic and p-values, but also sensitive to model choice, showing drastically more correlation between body size and environmental data.

| DISCUSS ION
We focus on evolutionary patterns at shallow phylogenetic scales in the L. gracilis species complex with extensive taxon sampling.
Our results suggest that two major lineages (Figure 2) whose body shape and size variation and possible variation in the vegetative structure of the habitats they experience have diversified recently, with phenotypic diversification increasingly prominent toward the  ; Table 1). However, these two lineages do not just represent two significantly different lineages each with its own evolutionary trajectory (Table 2). We show that body shape is the most important aspect of phenotypic variation in this clade, an often underappreciated axis of phenotypic evolution in macroevolutionary studies of animals but especially in Liolaemus (Figure 3). Shape evolution is evolutionarily correlated with the occupation of habitats with various vegetation structures and degrees of biomass (Table 3). Lastly, we explore important analytical caveats when undertaking evolutionary correlation analyses with respect to phylogenetic uncertainty in rapidly radiating groups and the importance of model choice between dependent multivariate and univariate tests of phenotype-environment correlation.

| Phylogenetic relationships and timing of diversification
The L. gracilis species complex, which includes much of the environmental diversity and geographic diversity represented by the Liolaemus clade, comprises northern and southern lineages (Figures 1 & 2) that our results show diversified over the last ~4-6 million years. The evolutionary history of these two lineages 107.90 ± 1.26 2.57-2.57 β = −1.72 ± 1.14 Note: Model fit was assessed using AICc calculated from across 1000 posterior trees. Shown are the mean and 95% confidence intervals on model weight (brackets), the mean ± standard deviation of AICc, the 95% confidence intervals on AICc, and the mean ± standard deviation of the ⍺ and β parameters for the OU and EB models, respectively. differs, with an older phylogenetic history ( However, southern Patagonian habitats were either under permafrost or shifted off the current coastline by up to four degrees during glaciation (Ponce et al., 2011;Rabassa, 2008;Rabassa et al., 2005;Rabassa et al., 2011), which could have resulted in the shallower phylogenetic history of the southern L. gracilis clade. Nevertheless, deeper phylogenetic histories were found in the Liolaemus lineomaculatus (Breitman et al., 2012) and Liolaemus elongatus (Medina et al., 2018) groups suggest that this shallow phylogenetic history may not be entirely due to a lack of refugia in the south. Shallow histories may simply be a consequence of a natural southern expansion of the southern L. gracilis clade or competitive exclusion by other diverse Liolaemus clades, further investigation is warranted here.

| Patterns of ecological and phenotypic diversification
Despite long being recognized as a source of error in comparative analyses (Huelsenbeck et al., 2000), methods accounting for phylogenetic uncertainty have only recently become available.
Phylogenetic uncertainty within Liolaemus is not entirely unexpected, given the well-documented cases of interspecific introgression observed among species (e.g., Esquerré et al., 2021;Olave et al., 2011;Olave et al., 2018), the rapidity with which species have accumulated, and the small number of loci used for phylogenetic studies for most species within Liolaemus species groups (Heled & Drummond, 2010). We found that phylogenetic uncertainty impacts been observed in other lizards (Edwards et al., 2015), but such studies are rare in other animal groups.
For body size and shape, the impacts of phylogenetic uncertainty do not actively impede the interpretation of patterns of diversification, where both shows increasing diversification toward the present (Figure 3). Our results indicate that body shape is the most labile aspect of phenotypic diversification in the L. gracilis clade, much more so than body size. There are many studies that suggest body shape diversification has been important for multiple vertebrate groups (Allen et al., 2013;Bonett & Blair, 2017;Friedman et al., 2019Friedman et al., , 2020Goodman et al., 2009;Gray et al., 2019;Vanhooydonck & Damme, 1999). Therefore, contrary to numerous studies and higher phylogenetic scales (i.e., Espinoza et al., 2004;Harmon et al., 2003), our study at a shallower phylogenetic scale and with detailed sampling, does not support the hypothesis that body size is a major axis of variation in Liolaemus. This result concurs with more recent studies with more extensive taxon sampling, showing that body size evolution does not correspond to diversification in Liolaemus (Olave et al., 2020). Note: We provide the mean (large type) and 95% confidence intervals (brackets) for the test statistic, and mean and 95% confidence intervals (brackets) for slope significance p-values. Parameter estimates are taken from 1000 randomly drawn trees from the posterior distribution. Trait correlations that are significant are bolded. p-value ≤ .01 **. The proportions of variance described by multivariate axes (i.e., vegetation, temperature, and precipitation) are shown (large type proportions) along with the main loaded variables on those axes (see Table S4). and invertebrates in terms of patterns of body size diversification (Cooper & Purvis, 2010;Itescu et al., 2018;Laurin, 2004;Pallarés et al., 2019;Rainford et al., 2016). The finding that body shape evolution may be more extensive and possibly independent of body size evolution has been found in fishes (Colombo et al., 2015;Friedman et al., 2019), but see (Hendry et al., 2006), and lizards (Goodman et al., 2009). Although such studies that compare body shape diversification directly to body size diversification are rare, and our results imply these may warrant further study.

| Mode of phenotypic trait evolution
We show that the evolution of environmental and phenotypic traits, individually, clearly follow an OU model (Table 1) and we find no evidence that the northern and southern lineages significantly differ in the environmental and phenotypic axes tested (Table 2). Our results indicate that the traits we tested seem to be evolving according to adaptive processes, which is consistent with analyses of body size data across Liolaemus (Olave et al., 2020). This study is the first to show that an OU process fits the evolution of a variety of environmental and phenotypic features, which has not been tested broadly across Liolaemus. This has largely been due to shape either not being considered variable (Jaksić et al., 1980) or irrelevant to adaptive variation due to a lack of phenotype-environment correlations (Schulte et al., 2004). This warrants further investigation across liolaemids. OU models generally model stabilizing selection toward an adaptive optimum (Beaulieu et al., 2012), but nevertheless have their limitations, including high Type-I error rates with limited sample sizes and difficulty interpreting parameter values (Cooper et al., 2016).
Therefore, we must view our results within these caveats.

| The evolutionary relationships between environment and phenotype
Our results show that body size has limited or no evolutionary relationship with any environmental variable tested in L. gracilis species (Table 3; but see discussion below of potential methodological issues). This contrasts with studies in other Liolaemus species where correlated evolution between body size and temperature have been reported (Azócar et al., 2015;Espinoza et al., 2004), but see (Jaksić et al., 1980;Pincheira-Donoso et al., 2007). Such variation in evolutionary relationships between body size and environmental variables among ectotherm clades is common because habitat preferences may vary and evolutionary lability in size is not necessarily consistent across clades (Pallarés et al., 2019). Our results show a clear and robust correlation between body shape and vegetation, specifically biomass average and seasonality (Table 3). Yet, studies across Liolaemus have failed to find consistent evidence that the evolution of environmental and body shape traits are related (Jaksić et al., 1980;Pincheira-Donoso et al., 2009;Schulte et al., 2004;Tulli et al., 2012). Body shape variation is particularly important for functional locomotory performance in lizards (Bergmann & Irschick, 2010;Garland & Losos, 1990;Melville & Swain, 2000) that are also important in other animals (Irschick et al., 2008). High biomass and low seasonality versus low biomass with high seasonality typify the Argentinian humid northwest Andes, and the southern Patagonian cold and arid steppe environments, respectively (Fabricante et al., 2009;León et al., 1998). Therefore, such differences in habitat structure likely impact escape behavior, which has been associated with body shape in Liolaemus (Schulte et al., 2004), and selection of body shape in these vastly different environments.

| Phylogenetic uncertainty and model choice in studies of evolutionary correlation
One important limitation of our study appears to be the robustness  (Table S6). This contrasts with studies indicating that the univariate phylogenetic regression is robust to tree misspecification, but can be impacted by extreme variation in topological and branch length reconstruction (Stone, 2011). Likely important here given hybridization in Liolaemus (Esquerré et al., 2021;Olave et al., 2018). We are not aware of any such explorations of the robustness of multivariate regression approaches considering tree misspecification, nevertheless, our results suggest they may be more robust to phylogenetic uncertainty.
Model misspecification has been noted with increasing trait dimensions in multivariate analyses (Adams & Collyer, 2018, 2019, and it may be that our dataset here is small enough to avoid such artifacts. Nevertheless, model misspecification is clearly noted as impacting the interpretation of univariate dependent trait phylogenetic regression analyses and an OU model is often chosen for such datasets despite not being able to adequately explain the correlation structure (Pennell et al., 2015). Early burst models are rare in comparative data , therefore, it may be that a BM model does explain the correlation structure of phenotypeenvironment correlations here, and that body size does correlate with environmental features. Regardless, our results indicate that further explorations of the relative robustness of multivariate and univariate phylogenetic regression approach with respect to model and tree misspecification may be warranted.

| CON CLUS IONS
Body shape is an important predictor of performance, escape behavior, and an animals function in its environment. We show that body shape has diversified more than body size in the L. gracilis clade and that the evolution of shape in this group is correlated with vegetation, likely reflecting a relationship with locomotory performance and escape tactics in different structural environments. While our results for relationships between environmental variables and body size are less clear, body size has diversified toward the present in this clade. It will be important for future studies to consider the relative importance of size versus shape in adaptive radiations. Not incorporating phylogenetic uncertainty has likely impacted the interpretation of adaptive evolution in Liolaemus but also many other adaptive radiations. It is critical that in analyses of adaptive radiations and rapidly radiating clades that the impact of phylogenetic uncertainty be included, particularly because topological and branch length variation is highest in these clades. Interestingly, fundamental differences between methods testing multivariate versus univariate dependent variables may fundamentally differ in their sensitivity to model and tree misspecification and should be explored through simulation studies, ideally also with a unified statistical approach that can incorporate both these data types.

ACK N OWLED G M ENTS
We thank fauna authorities from Santa Cruz, Chubut, Río Negro, Neuquén, Mendoza, and La Pampa provinces for assistance with permitting. Funding for this project came from the National Science

CO N FLI C T O F I NTE R E S T
The authors declare no conflict of interest.

DATA AVA I L A B I L I T Y S TAT E M E N T
The data that support the findings of this study are available in Dryad (https://doi.org/10.6071/M36X1M). Some data (Genbank numbers, morphological, and spatial data) are also in the Supplementary Material, which can also be found in Dryad.