Climate drives global functional trait variation in lizards

A major challenge in ecology and evolution is to disentangle the mechanisms driving broad-scale variation in biological traits such as body size, colour, thermal physiology traits and behaviour. Climate has long been thought to drive trait evolution and abiotic filtering of trait variation in ectotherms because their thermal performance and fitness are closely related to environmental conditions. However, previous studies investigating climatic variables associated with trait variation have lacked a mechanistic description of the underpinning processes. Here, we use a mechanistic model to predict how climate affects thermal performance of ectotherms and thereby the direction and strength of the effect of selection on different functional traits. We show that climate drives macro-evolutionary patterns in body size, cold tolerance and preferred body temperatures among lizards, and that trait variation is more constrained in regions where selection is predicted to be stronger. These findings provide a mechanistic explanation for observations on how climate drives trait variation in ectotherms through its effect on thermal performance. By connecting physical, physiological and macro-evolutionary principles, the model and results provide an integrative, mechanistic framework for predicting organismal responses to present climates and climate change. The authors incorporate mechanistic information about lizard physiological responses to heat into predictions of trait variation across time and space, finding that the range of functional traits is more constrained in locations where the local climate strongly selects for thermal performance.

A major challenge in ecology and evolution is to disentangle the mechanisms driving broad-scale variation in biological traits such as body size, colour, thermal physiology traits and behaviour. Climate has long been thought to drive trait evolution and abiotic filtering of trait variation in ectotherms because their thermal performance and fitness are closely related to environmental conditions. However, previous studies investigating climatic variables associated with trait variation have lacked a mechanistic description of the underpinning processes. Here, we use a mechanistic model to predict how climate affects thermal performance of ectotherms and thereby the direction and strength of the effect of selection on different functional traits. We show that climate drives macro-evolutionary patterns in body size, cold tolerance and preferred body temperatures among lizards, and that trait variation is more constrained in regions where selection is predicted to be stronger. These findings provide a mechanistic explanation for observations on how climate drives trait variation in ectotherms through its effect on thermal performance. By connecting physical, physiological and macro-evolutionary principles, the model and results provide an integrative, mechanistic framework for predicting organismal responses to present climates and climate change.
Climate imposes constraints on organismal trait diversity and evolution through its effects on physiological performance and fitness 1,2 . As climatic conditions vary across space and time, these constraints also change, leading to sometimes predictable variation in functional traits such as body size and shape, skin colour, physiological thermal tolerance and behaviour (ecogeographical rules 3,4 ). Understanding the mechanisms underpinning these patterns lies at the core of macroecology and macrophysiology 5 and represents a milestone towards explaining and predicting organismal responses to past, present and future climates. Yet, predicting how climate modulates functional trait diversity and evolution is challenging. This task requires a thorough understanding of the mechanisms that govern organism-environment interactions 2 , environmental information at the proper scale of these interactions 6,7 and a framework that enables evaluation of the joint influence of multiple traits on physiological performance and fitness across environmental gradients.
In ectotherms, environmental temperature plays a fundamental role in regulating physiological performance and fitness 8 . Temperature is thought to impose selective pressures on functional traits involved in body temperature regulation that give rise to broad-scale spatial gradients in body size 9 , body shape 10,11 , skin colour 12 , thermal tolerance limits 13,14 and behaviour 15,16 . However, these traits interact with each other, obfuscating how we might predict them to be altered by selection when allowed to mutually vary. For example, behavioural adjustments of body temperature can buffer the evolution of thermal tolerance traits, a phenomenon known as the Bogert effect [16][17][18] . Previous Article https://doi.org/10.1038/s41559-023-02007-x observed in places where the slope of performance sensitivity is positive for that trait and smaller trait values should be observed in places where the slope of performance sensitivity is negative for that trait. Second, under the hypothesis that selection is stronger in places where performance sensitivity to trait variation is higher, we predict that trait variance will be greater in places where performance sensitivity to trait variation is closer to zero (weak selection) and less variable where performance sensitivity to trait variation is high in absolute value (strong selection) 20 . We can apply these predictions to macro-evolutionary processes considering variation among species within a phylogeny or to environmental filtering processes by considering variation and means among co-existing species. Hence, we evaluated the first prediction, that sensitivity drives trait values across species, by analysing the relationship between species' mean trait values and their mean sensitivity (averaged across each species' geographical range). We also evaluated this prediction at the assemblage level by computing the mean trait value of all the species that co-exist in a grid cell and analysing its relationship with the sensitivity estimated in that cell. We evaluated the second prediction, that performance sensitivity constrains trait variance, only at the assemblage level by analysing the relationship between the variance of traits within-species assemblages in each grid cell and the performance sensitivity to that trait at that location.

Prediction I for sensitivity drives trait values
The species-level analyses comparing the functional traits of lizards with the sensitivity values averaged across their geographical range revealed positive relationships for body mass, preferred temperature, field body temperature and critical thermal minimum with their respective sensitivity values ( Fig. 3 and Table 1). The relationship was non-significant for the critical thermal maximum. Preferred and field body temperatures, both used as indicators of parameter T pref in the model, showed similar relationships with T pref sensitivity. The phylogeny generally explained a large amount of variance suggesting that historical factors and the phylogenetic relatedness among species played an important role in determining functional trait variation (Table 1). Overall, these results support the prediction of a positive relationship between trait values and sensitivity-a possible proxy for the strength and direction of selection-at the species level (Fig. 3). Interestingly, the analyses using sensitivity as predictor outperformed correlative models based on climatic variables by improving the explanatory capacity of the linear models ( Supplementary Information 4).
At the assemblage level, the relationship between mean trait values and sensitivity was weaker overall than for species-level analyses. We found no significant relationship for mean body mass, critical thermal maximum or critical thermal minimum, although the relationship between field body temperature and T pref sensitivity was significant and positive (Supplementary Information 5).

Prediction II for sensitivity constrains trait variance
Species assemblages facing higher sensitivity (in absolute value) to changes in body mass displayed lower variance in body mass, whereas assemblages facing lower sensitivity values showed either low or high variances. The relationship between trait variance and sensitivity was also triangular for preferred temperature, field body temperature and critical thermal minimum. That is, the variance in these traits was delineated by quantile regression across the upper 0.9 quantile (Fig. 4), suggesting that sensitivity imposes an upper limit on trait variance. In contrast and against expectations, the relationship was positive for critical thermal maximum (Table 2).

Discussion
Our global analysis of performance sensitivity to trait variation predicted geographical patterns in the direction and strength of selection on different functional traits involved in thermoregulation and thermal performance in lizards. The sensitivity analysis provided new studies have focused on the analysis of statistical associations between functional traits and climatic variables, which has allowed us to identify the main environmental factors hypothesized to drive these patterns. However, these correlative methods lack a mechanistic representation of the physical and biological processes underpinning these patterns, which limits their capacity to predict the responses of organisms to new environmental conditions.
Here, we propose a mechanistic approach combining a biophysical model and a sensitivity analysis to numerically predict the effects of different functional traits (body size, skin colour, thermal tolerance limits, preferred body temperature and thermoregulatory capacity) on thermal performance of ectotherms (Fig. 1). The biophysical model calculates body temperature and thermal performance as a function of both climatic conditions and multiple independently varying functional traits ( Fig. 1; Methods; Supplementary Information 1-3). Then, the sensitivity analysis quantifies the effect of each functional trait on thermal performance by assessing the relationship between changes in each trait value and the resulting changes in thermal performance. By combining this approach with spatially explicit microclimatic data 19 , we compute a putative direction and strength of selection on each trait across the globe. We then test two hypotheses about the evolution and assembly of temperature-related phenotypes. The first is that selection has acted on trait variation to maximize thermal performance and the second is that selection strength on traits varies in proportion to their predicted impact on performance. We test these hypotheses by comparing model-generated predictions on the direction and strength of selection with empirical trait and assemblage data of global lizards.

Results
We analysed the sensitivity of cumulative thermal performance to changes in body mass, skin colour (absorbance to short-wave radiation), preferred temperature (T pref ), thermoregulatory ability (λ), critical thermal maximum (CT max ) and critical thermal minimum (CT min ; Fig. 1). Sensitivity values, that is, the slope of the predicted relationship between change in performance and change in each trait value, displayed geographical patterns that recapitulate patterns expected with climate. For body mass, sensitivity values were negative in tropical regions meaning that increasing body mass reduced thermal performance in these regions. These negative sensitivity values vanished towards higher latitudes suggesting that changes in body mass did not affect thermal performance in colder environments (Fig. 2a). Conversely, skin absorbance showed stronger, positive sensitivity values in colder regions, negative in subtropical regions and deserts and weak, negative relationships in tropical regions (Fig. 2b). These results suggest that a higher skin absorbance improves thermal performance in colder environments since it increases body temperature by raising absorption of solar radiation, whereas a lower skin absorbance would be favoured in warmer environments. Sensitivity values close to zero in the tropics suggest that the effect of skin absorbance on thermal performance is weak in these regions.
Sensitivity values of behavioural thermoregulation were higher in the tropics, suggesting that tropical ectotherms are more pressed to increase T pref and thermoregulatory ability (λ) than are temperate organisms (Fig. 2c,d). Similarly, CT max exhibited higher, positive sensitivity values at lower latitudes, suggesting that increasing heat tolerance improves performance in these regions (Fig. 2e). Finally, CT min sensitivity was higher and negative towards colder regions, as expected to reflect greater performance with lower CT min values (Fig. 2f).

Model predictions for direction and strength of selection
Our model allows us to make predictions about trait variation according to two hypotheses. First, under the hypothesis that selection on trait variation has acted to maximize thermal performance, we predict that trait values among species will vary in proportion to their predicted impact on performance. In other words, larger trait values should be to simulate body temperature and thermal performance of a lizard-like ectotherm that thermoregulates behaviourally by moving between sunexposed and shaded conditions. The model uses information on body mass, skin absorbance, preferred temperature (T pref ), thermoregulatory ability (λ) and critical thermal limits (CT max and CT min ). b, To perform the sensitivity analysis of the biophysical model, we changed the value of one trait at a time drawing random values from their input distributions. Then, the effect of each trait on cumulative performance is estimated by computing the slope of the relationship between the change in performance and the change in each trait value. This process is repeated across locations (grid cells) on a global map to investigate emerging geographical patterns of trait sensitivity. Article https://doi.org/10.1038/s41559-023-02007-x hypotheses for where relevant phenotypes are expected to be under strong versus weak selection and the direction of selection. Although the observed patterns of trait variation can still preserve the imprint of historical processes acting on species' distributions and trait evolution 21 , our findings provide evidence that selective regimes on body mass, skin absorbance, thermoregulatory ability and thermal tolerance limits vary predictably with climate across space. Our model predicts the direction and strength of selection of multiple traits at once, providing numerical predictions that can be tested against empirical data. This work thus provides a first step in building an integrated view of how climate modulates functional trait evolution and composition of species assemblages through its effect on thermal performance.
The sensitivity of thermal performance to each evaluated trait exhibited strong geographical gradients. For instance, sensitivity to changes in body mass was negative in tropical latitudes, meaning smaller sizes had higher cumulative performance than larger sizes and approached zero (a flat, invariant relationship) towards colder latitudes, meaning body size variation had no effect on performance. This pattern suggests that selective pressures favouring smaller body sizes in warm latitudes weaken at higher latitudes. Hence, the selective pressure towards maximizing thermal performance seems to drive body size evolution among tropical species. The effect of climate on large-scale patterns in body size might arise through the direct influence of temperature accelerating growth rates and reducing maximum size (a proximate cause, for example, ref. 22 ) or via selection on maximum body size to avoid overheating in warm environments (an ultimate cause). Our findings build on previous work in macroecology focused on the ultimate causes of body size evolution suggesting, for example, that larger sizes evolve to facilitate heat conservation in cold environments 9 . Yet, latitudinal body size patterns in lizards have been previously disputed, with empirical analyses reporting negative 9 , positive 23 or even no relationships between body size and environmental temperature 24 (Supplementary Table 4). Our mechanistic analysis resolves this lack of consensus, by predicting different relationships to arise from spatial variation in the relative importance of body size evolution upon thermal performance. Mechanistically, large bodies reduce convective heat dissipation relative to the absorption of solar radiation, increasing body temperature and overheating risk 25,26 , especially in the tropics, where solar radiation levels are higher 27 . Among temperate species, however, the many factors that might drive body size evolution (for example, historical, ecological and constraints on life-history traits 9,28,29 ) might supersede the constraints linked to overheating risk. Yet, the idea that larger species experience higher overheating risk due to lower effectiveness of convective heat dissipation has been questioned 30 , since larger animals might also experience lower temperatures due to greater height above the ground. Although we used size-specific microclimates in our modelling, our results might be sensitive to our specific assumptions about allometric relationships between body size and height above the ground 30 and this could lead to an overestimate of size sensitivity in the tropics in our models.
Under the perspective of our mechanistic model, the empirical analyses supported the predicted association between body size and climate at the species level (Fig. 3a). At the assemblage level, we found Article https://doi.org/10.1038/s41559-023-02007-x no significant relationship between mean body mass and sensitivity (Supplementary Table 5), whereas the relationship between the assemblage variance and absolute sensitivity was negative as predicted (Fig. 4a). Because the accuracy of the estimated mean body mass at the assemblage level depends on the assemblage variance, this change in variance in relation to sensitivity might affect the observed relationship between mean assemblage body mass and sensitivity. Our predictions could be further tested with intraspecific trait variation data as they become available. Assemblages at environments with higher sensitivity to body mass variation displayed lower variance in body mass, depicting a triangular relationship between variance and sensitivity (Fig. 4a). That is, the variances at the assemblage level could be either high or low in regions with lower sensitivity to body mass but only low in regions of higher predicted climatic constraints. These results suggest that thermally related selection on body size influences body mass evolution across lizard species in warm environments and constrains assemblage composition constituting an environmental filter. The variation in different functional traits may generate different adaptive pathways to maximize thermal performance across climatic gradients. For example, light skin colours (low skin absorbance)  Table 1). The predicted geographical patterns of skin absorbance supported the thermal melanism hypothesis, which proposes that a darker colouration (low reflectance or high absorbance) is beneficial in cold environments because it accelerates heating rates and increases body temperature 12 . In addition, the thermal melanism hypothesis predicts that-once potential confounding factors such as crypsis, fecundity and life history are controlled for-skin absorbance might be negatively related to body mass, with larger species showing lighter colours 12 . Our model and analyses support this prediction and suggest that the expected relationship between body size and skin colour may vary geographically. Thus, although reducing body size lowers the risk of overheating in the tropics, changes in body size may not significantly affect thermal performance in temperate and subtropical latitudes. Instead, animals at these latitudes may effectively increase body temperature via increasing skin absorbance. Thus, tropical species may experience selective pressures to reduce skin absorbance and body size (potentially leading such variables to covary), whereas the skin colour-body size correlation is expected to vanish among temperate and subtropical species because body size variation in these regions is less influenced by thermal constraints. In the end, the actual interspecific relationship between these traits likely depends on the evolutionary lability of each trait, time and the relative importance of other selective forces such as crypsis or fecundity 31,32 .  The sensitivity analysis predicts that animals living in warmer latitudes select higher body temperatures (T pref ) than those living in colder regions. Similarly, thermal tolerance limits are predicted to display latitudinal gradients, with CT max increasing towards warmer latitudes and CT min decreasing towards colder environments. Following the expected effect for T pref , lizards showed higher preferred temperatures in warm environments (Fig. 3b) and the assemblage-level variance in preferred temperature was higher in regions where T pref sensitivity values were lower (Fig. 4b). The empirical patterns of heat tolerance, however, did not followed the prediction of the model. Thus, the mean heat tolerance of the species was unrelated to CT max sensitivity and the assemblage-level variance in this trait increased in relation to sensitivity. This result is consistent with previous observations that found invariance in heat tolerance across latitudes 15,[33][34][35][36][37] . The mismatch between the observed and predicted patterns for CT max indicates that other factors might be driving or constraining the evolution and geography of heat tolerance. One explanation is that thermal tolerance limits-generally measured as the temperature at which locomotor activity declines-are not relevant to fitness. However, the observed patterns for the cold tolerance limit were consistent with the model's prediction for CT min suggesting that it is relevant to fitness and modulated by climate. Therefore, an alternative explanation for the mismatch between observed and predicted patterns for CT max is that heat tolerance is relevant to fitness but evolutionarily constrained because elevated temperatures compromise biochemical stability of membranes and proteins, imposing a hard boundary to heat tolerance evolution 14 . Another explanation involves the role of thermoregulatory behaviour through the Bogert effect, under which greater behavioural capacity to control body temperature may inhibit selection on CT max 16,18 . Indeed, our model suggests that both CT max and thermoregulatory ability are more strongly selected in the tropics which might open two alternative evolutionary pathways: either increasing CT max or improving thermoregulatory ability to avoid overheating in exposed conditions. In tropical regions, the greater availability of cool and stable microclimates under the forest canopy may facilitate lizards to behaviourally avoid high body temperatures in sun-exposed environments inhibiting selection on CT max 36 . This work offers a mechanistic framework to predict the direction and strength of the effect of climate on traits involved in thermoregulation and thermal performance in ectotherms. Our approach predicts broad-scale patterns of trait variation, providing numerical descriptions of patterns that might be observed as 'ecogeographical rules' (such as geographical clines in body size or skin colour across latitudes) from their fundamental biophysical and physiological mechanisms. By deepening how we consider the mechanisms underpinning lizard-microclimate interactions, the sensitivity analysis offers insights into the climatic drivers of trait variation that would not be possible to test using traditional, correlative models. Mechanistic predictions demonstrate that not only the direction but the strength of selection on different functional traits varies with climate, explaining previous contradictory results on how functional traits like body mass should covary with latitude in terrestrial ectotherms. Finally, the sensitivity analysis raised new hypotheses and provided a way to numerically evaluate the environmental filtering effect of climate on trait variance in species assemblages. Our framework connects concepts of interspecific phenotypic trait evolution with abiotic filtering of assemblages and patterns of functional diversity. We propose that this framework, complemented with further modelling work describing, for example, water exchange and energy assimilation rates, provides a key step forward towards unifying frameworks in ecology and evolution.

Modelling body temperature and thermal performance
Body temperature in ectotherms depends on the heat exchange between the body and its environment and thus on morphological traits such as body mass, shape or skin colour and on microclimatic conditions such as solar radiation or air temperature at sun-exposed or shaded conditions. The variation in experienced body temperature over a time interval determines cumulative thermal performance or fitness 33,38,39 . To develop a mechanistic model that captures body temperature variation and its effect on thermal performance, we needed to arrange (1) a transient-heat model that computes heat exchange between the animal and its environment as a function of microclimatic conditions and the animal's morphology and colour, (2) a behavioural model to simulate the capacity to select different micro-environments behaviourally (for example, sun or shade) to maintain body temperature close to the preferred temperature and (3) a physiological model to compute cumulative thermal performance as a function of experienced body temperature.
We developed a transient-heat model to compute changes in body temperature (T b ) as a function of the absorbed solar radiation, the exchange of thermal radiation with the environment, convective cooling and conductive heat exchange with the soil surface 25,26,40 (Supplementary Information 1). To compute heat flow and body temperature, the model uses information on the microclimatic conditions (solar radiation, air and soil temperature and wind velocity 19 ), together with morphological and radiative properties of the body (body mass, skin surface area, body length, skin absorbance and emissivity of thermal radiation). We used this model to calculate how body temperature varies over time in one micro-environment or when the animal moves between sites.
To model behavioural thermoregulation, we considered that animals prefer the micro-environment (sun-exposed or shade) closer to their preferred temperature (T pref ) 41 . Animal's decisions on what micro-environment is selected at each instant of time balance multiple costs and benefits, not only on thermal performance but also on nutritional needs or predation risk 38,42 . To address this complexity, behavioural thermoregulation can be modelled as a probabilistic process using the maximum entropy framework where the thermoregulatory constraint weights-but does not entirely determine-the probability of selection of each micro-environment in the repertoire 26 . Thus, the probability that the animal selects sun-exposed conditions depends on the difference between T pref and body temperature in the sun and a non-dimensional parameter, λ, that characterizes thermoregulatory ability: 1) and the probability of selecting the shade is 1 − P(T b,sun , λ, T pref ) 26 . The denominator, that is, the summation across all possible micro-environments (here, sun-exposed and shade conditions; subscript j is the index of summation across micro-environments), normalizes probabilities so that the sum of all probabilities is equal to 1. Both the preferred body temperature (T pref , °C) and parameter λ weight the importance of thermoregulatory behaviour in determining shifts between micro-environments. When λ is 0, the probabilities of sun-exposed and shaded conditions are the same: P(T b,sun , λ = 0, T pref ) = P(T b,shade , λ = 0, T pref ) = 0.5, meaning that the animal shifts between micro-environments irrespectively of what its body temperature is at each site. By contrast, high λ values led to an increase in the probability of selection of the micro-environment that is closer to the animal's preferred body temperature 26 . Together, the behavioural thermoregulation and transient-heat models generate a series of body temperatures over a time interval, which may be expressed as a body temperature distribution n(T b ) representing the probability density of body temperatures as a function of microclimatic conditions, animal's morphology, radiative properties, T pref and thermoregulatory ability, λ. This distribution can be obtained by numerically integrating the master equation of the system to find Article https://doi.org/10.1038/s41559-023-02007-x the stable T b distribution across every hour of the day (details in ref. 26 ). Then, daily T b distributions are generated by adding hourly distributions for the period of interest, for example, during the daytime. This approach reduces computational time and provides accurate predictions of daytime T b distributions of lizards 26 . Validation analyses suggest that the model accurately predicts body temperatures observed in the field (within-species validation, R 2 values mostly ranging between 60% and 90%; among-species validation, R 2 = 0.86; Supplementary  Information 3).
Finally, we were interested in modelling how experienced body temperatures determine cumulative performance or fitness over time.
To do so, we modelled physiological performance using thermal performance curves (TPCs), which describe thermal dependence of physiological processes such as locomotion or digestion and thus assume that fitness varies with temperature following a TPC 33,38 . TPCs are typically asymmetric curves with their maximum shifted towards high body temperatures. Thus, performance declines rapidly for body temperatures exceeding the optimum and more gradually towards temperatures below the optimum (Fig. 1). As such, TPCs can be characterized by three parameters: the body temperature at which maximum performance is attained (T opt ) and the critical thermal minimum (CT min ) and maximum (CT max ), representing temperatures at which thermal performance becomes zero. The cumulative performance over a time period is the integral of thermal performance weighted by the distribution of body temperatures over that period 33,38,39 . Here, we used the equation by ref. 33 and computed relative performance (ranging between 0 and 1) as a function of CT max and CT min . We assumed that T opt corresponds to T pref (see ref. 43 for a discussion) to reduce the complexity of the model and because empirical T opt data are scarce in the literature. We then calculated cumulative performance using the simulated T b distribution, n(T b ), by computing the integral ∫TPC(T b )n(T b )dT b .
Microclimatic conditions in the sun and the shade were modelled using the NicheMapR microclimate model 19 in R v.3.5.3 (ref. 44 ). To simulate conditions in the shade, NicheMapR requires input information on the percentage of shade in that location. To estimate the level of shade, we modelled the percentage of solar radiation intercepted by the canopy using remotely sensed data on the leaf area index (LAI) 45 , obtained from a monthly averaged global gridded database of LAI worldwide 46 . Specifically, the proportion of solar radiation that makes it through the canopy was modelled as an exponential function of LAI as ω p = exp (− √ 0.8K b,c LAI) where K b,c is the extinction coefficient for direct solar radiation at zenith angle z, K b,c = √1+tan 2 z 2.00132 (ref. 45 ). Then, we used the value ω p to compute the maximum shade level in the Niche-MapR microclimate model and correct the radiation level under the canopy. To parametrize height above ground in the microclimate model, we estimated the height of a lizard as a linear function of body mass (Supplementary Table 1) thus considering size-specific microclimates 30 . With this approach, we computed solar radiation (W m −2 ), air temperature (°C), soil surface temperature (°C), wind speed (m s −1 ) and relative humidity (%) every hour for the middle day of each month. To reduce computational times, we focused the simulations on the daytime and on the period of activity, which was determined to be months in which maximum temperature (operative temperature in the sun) was higher or equal to T pref 47 .

Sensitivity analysis
We are interested in quantifying the effects of different functional traits on cumulative performance and how these effects change across climatic gradients. That is, we aim to assess the sensitivity of cumulative performance to changes in six traits: body mass, skin absorbance to short-wave radiation, preferred temperature (T pref ), thermoregulatory ability (λ), critical thermal maximum (CT max ) and critical thermal minimum (CT min ). Sensitivity analysis is a well-stablished approach in ecological modelling for analysing the effect of input parameters (here, functional traits) on the output of the model (here, cumulative thermal performance) 48,49 . For numerical models, this analysis often requires running simulations using different combinations of parameter values to evaluate their effect on the output of the model. Because simulating each combination of parameter values requires substantial computational effort, exploring the total factorial space is effectively impossible. Therefore, sampling methods are required to navigate effectively throughout the multidimensional space of parameter values. Here, we used the Latin hypercube sampling method in a one-at-a-time sensitivity analysis design [48][49][50] . The method starts by sampling a set of start values from their input distributions and computing cumulative performance for this set of trait values. The second step changes the value of only one variable and calculates the resulting change in cumulative performance as the difference in performance between the second and the first run. The process continues changing one variable each time and calculating the change in cumulative performance with the previous time step until all variables have been changed. We ran this process 20 times, with a different set of start values each time, totalling 140 runs (20 × (k + 1) where k = 6 is the number of input traits). To obtain an efficient sampling of these traits, we used the Latin hypercube sampling method (using the R package lhs) 50 , which provides a near-random, well-distributed walk through the hypercube of parameter values. We tested the repeatability of sensitivity estimations by running the model 20 times, obtaining a different set of starting values for each run, across ten randomly selected cells. The repeatability of the slopes was high for all traits: 85.6% ± 0.081 (mean ± s.d.) for body mass; 99.1% ± 0.008 for skin absorbance; 98.3% ± 0.016 for T pref ; 93.7% ± 0.047 for thermoregulatory ability; 95.9% ± 0.032 for CT max ; and 99.2% ± 0.007 for CT min . We parameterized trait distributions using values that are representative for lizards: log-normal distributions for log body mass (values used for parameterization are given in Supplementary Table 1) and thermoregulatory ability; and normal distributions for skin absorbance, T pref , CT max and CT min (Supplementary Table 1).

Data collection from lizards
We gathered mean species' trait data of diurnal, heliothermic lizards, including body mass 51 (estimated using family-specific allometric equations based on snout-vent lengths), preferred temperatures 52 and thermal tolerance limits 53 . We also used field body temperature 51 as an indicator of the model parameter T pref (relationship between field and preferred temperatures: r = 0.70, t = 5.9, d.f. = 36, P < 0.01) because field temperature was available for more species in our database. Predictions for skin absorbance and thermoregulatory ability (λ) were not analysed due to the lack of empirical data for these traits among heliothermic lizards. Finally, we used the geographical ranges of lizards in ref. 54 and the phylogeny by ref. 55 . Our database included 1,219 lizard species, with body mass being available for 1,215 species, preferred temperature for 49 species, field body temperature for 351 species, critical thermal maximum for 103 species and critical thermal minimum for 68 species.

Data analyses
To perform the species-level analyses, we used phylogenetically controlled linear models (phylogenetic generalized least squares, PGLS) to test the prediction of a positive relationship between each trait with its corresponding sensitivity value averaged across the geographical range of the species. All statistical analyses were evaluated using two-sided tests. Thus, we compared mean species body mass versus body mass sensitivity across the species' range; preferred temperature and field body temperature versus T pref sensitivity; critical thermal maximum versus CT max sensitivity; and critical thermal minimum versus CT min sensitivity. We then computed pseudo-R 2 values to estimate the variance explained by either the ecological factor (trait sensitivity) or the phylogeny using variance partitioning methods for models with autocorrelated residuals 56 . In an additional analysis, we Article https://doi.org/10.1038/s41559-023-02007-x compared the explanatory capacity of linear models using either trait sensitivity or climatic variables as predictors to evaluate the potential of mechanistic variables to describe functional trait variation when compared with traditional, correlative approaches ( Supplementary  Information 4).
For the assemblage-level analyses, we projected the observed average trait values of species across their geographical ranges. Then, for each trait, we extracted both the mean and variance of the species assemblage at each location. Because the estimation of variance of the assemblage is sensitive to sample size (number of species with known trait data co-existing in a grid cell), we performed additional analyses controlling for missing data in species assemblages to evaluate whether our results are robust to incomplete trait information 57,58 (Supplementary Information 6). To analyse the relationship between trait variance and sensitivity, we used both ordinary linear regression and quantile regression (regression at the upper 0.9 quantile) 59 . Quantile regression estimates specified quantiles of the response variable (trait variance) and computed a linear relationship between estimated quantiles and the predictor variable (trait sensitivity) to detect a possible triangular relationship between these variables.
To control for spatial autocorrelation in assemblage analyses, we obtained Moran eigenvectors and included them as covariates in linear regressions. This method effectively reduces the spatial autocorrelation and type-I error in the assemblage-level analyses by accounting for the spatial dependence in the residuals 60,61 . We implemented Moran eigenvectors using the R package spdep, which selects the number of eigenvectors required to reduce the Moran's I index for regression residuals up to a tolerance set to 0.1 (ref. 61 ). To represent the relationships between spatially corrected variances and sensitivity values (Fig. 4), we first computed the residuals of the trait variance in relation to the Moran eigenvectors and then represented these residuals in relation to the sensitivity values. For the cross-species analyses, species' trait sensitivities were obtained by averaging sensitivity values across their geographical ranges. Therefore, the spatial component is not explicit in these analyses and spatial eigenvectors may not be accurate. However, to confirm that our results were not affected by spatial autocorrelation, we performed additional analyses including Moran eigenvectors computed from the centroids of species' geographical ranges. The results of these spatially corrected models were similar to those excluding spatial eigenvectors (Supplementary Information  7 and Supplementary Table 7).

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

Data availability
The data that support the findings of this study 62