Ecophysiological constraints outperform environmental predictors in forecasting climate-driven extinction risk of a Cerrado endemic lizard

Forecasting species responses to climate change using only environmental predictors may underestimate extinction risk because it ignores physiological constraints. We compared species distribution models for Micrablepharus atticolus built with environmental predictors, thermal performance, hours of activity, and a hybrid combination of these predictors. We projected habitat suitability to 2050 and 2070 using 12 global circulation models, three representative concentration pathways, and three land-use scenarios. Model accuracy, quantified using the True Skill Statistic (TSS), was highest for the thermal performance model, followed by hours of activity, environment-only, and the hybrid model. The thermal performance model predicted larger current suitable areas and substantially higher future extinction risk than the environment-only model. For M. atticolus, and likely for species with similar ecological and physiological traits, ecophysiology-based models provide more reliable TSS-based predictions than models based only on environmental predictors and show that ignoring ecophysiological information may underestimate climate-driven extinction risk.

Cumulative anthropogenic impacts on natural ecosystems have reached such levels that global ecological processes are being disrupted (Ruddiman, 2013; Steffen et al., 2018), and a major biodiversity crisis is unfolding (Barnosky et al., 2011; Ceballos et al., 2015). When confronted with unfavorable changes in their environments, natural populations may persist through phenotypic plasticity or adaptation (Chevin et al., 2013; Valladares et al., 2014) and by tracking areas of environmental suitability through migration and range shifts (Parmesan, 2006; Thuiller et al., 2008; Tomiolo and Ward, 2018); otherwise, they undergo local extirpation or global extinction (Sinervo et al., 2010). Forecasting species responses to such changes is critical to mitigate anthropogenic impacts and ensure their persistence, often accomplished using predictive statistical models (Mouquet et al., 2015; Pereira et al., 2010). Models based on environmental predictors are the most popular approach to forecasting species range shifts under future change (Elith et al., 2010; Elith and Leathwick, 2009). These models correlate species occurrences or abundances with environmental characteristics within their ranges to predict areas with suitable environmental conditions. However, they have been criticized for assuming that species are at equilibrium with their environments (Araújo and Pearson, 2005; Munguia et al., 2012) and for having low transferability—in space or time—to novel environments (Davis et al., 1998; Menke et al., 2009). By explicitly incorporating spatial formulations of relevant biological processes, such as physiology, dispersal, demography, and species interactions, models based on ecophysiological constraints attempt to avoid the shortcomings of the correlative approach; yet, they require data that are often costly and of limited availability (Connolly et al., 2017; Kearney and Porter, 2009; Leroux et al., 2013; Urban et al., 2016).

Addressing the merits of these different approaches is often accomplished through accuracy metrics based on the fit of model predictions to observed species distributions (Allouche et al., 2006; Liu et al., 2011). However, the most accurate model may not necessarily produce the best predictions under novel conditions, and some models may be more sensitive to sampling bias, which might also reduce model transferability (Randin et al., 2006; Rodda et al., 2011). Further, simpler models produce less accurate predictions but have improved transferability (Guisan et al., 2017). Some have advocated the integration of ecophysiological and environmental predictors into hybrid models, with the goal of improved model transferability for novel conditions (Buckley et al., 2011; Ceia-Hasse et al., 2014; Martínez et al., 2015; Rodríguez et al., 2019). Alternatively, different model predictions can be averaged using model accuracies as weights, which presumably accounts for the many sources of model uncertainties and produces more robust predictions (Buisson et al., 2010; Diniz et al., 2009; Watling et al., 2015). Which of these several approaches is more accurate in forecasting global change impacts on species range shifts is still a matter of debate (Dormann et al., 2012; Fordham et al., 2018; Lurgi et al., 2015; Morin and Thuiller, 2009).

The need to accurately forecast responses to environmental changes is more urgent for species from highly threatened regions, such as the Cerrado hotspot (Oliveira and Marquis, 2002). Due to the combination of high levels of biodiversity/endemism and high rates of conversion of natural habitats into crops and pastures, the Cerrado is a global biodiversity hotspot (Mittermeier et al., 1998; Myers et al., 2000). Nevertheless, protected areas cover only 7.7% of the Cerrado, more than 50% of the native habitats have been cleared, and increased pressure from agriculture expansion will likely wipe out the remaining Cerrado by 2050 (Beuchle et al., 2015; Oliveira et al., 2017; Strassburg et al., 2017). Previous studies predicting biodiversity responses to environmental changes in the Cerrado produced conflicting results. Diniz-Filho et al. (2009) predicted gains in areas of suitability for 76% of 753 species of vertebrates in 2050, but gains were predicted to occur in areas of intense cattle ranching activities. Conversely, climatically suitable areas for the endemic snake Phalotris lativittatus were predicted to shrink by 65-70% in 2050 and 80-90% in 2080 (Vasconcelos, 2014). Moreover, Cerrado bats were predicted to lose, on average, from 54.3% (not allowing dispersion) to 18.1% (allowing dispersion) of their areas of climate suitability (Aguiar et al., 2016). None of these studies, however, used ecophysiological information.

Lizards are model organisms for ecological and evolutionary studies (Pianka and Vitt, 2003). As terrestrial ectotherms, which depend on external heat sources for thermoregulation, lizards are highly susceptible to changes in the environment, particularly those that affect their thermal biology (Nowakowski et al., 2018). Here we used the lizard Micrablepharus atticolus (Squamata, Gymnophthalmidae), an endemic species of the Cerrado, to compare predictions of range shifts based on models that differ in environmental inputs. We generated habitat suitability models using (1) only environmental predictors, (2) only ecophysiological predictors, and (3) hybrid models including both kinds of predictors. Then, we projected models considering scenarios of climate change and land use for 2050 and 2070. Micrablepharus atticolus is a small-bodied, blue-tailed lizard that lives primarily amid the leaf litter; it is active during the hottest hours of the day, feeds on diverse groups of arthropods, breeds during the dry season, and lays multiple clutches during the year (Sousa et al., 2016, 2015; Vieira et al., 2000).

We assessed the potential accuracy gains by incorporating ecophysiological data on hours of activity and thermal performance into models predicting climate-induced range shifts and extinction risk. We also account for model uncertainty resulting from different modeling algorithms and climate change and land-use scenarios. Rather than providing a universal test of mechanistic versus correlative models, our goal is to evaluate how ecophysiological constraints shape climate-driven extinction risk in this specific biological and ecological context. We hypothesized that models explicitly incorporating ecophysiological constraints (thermal performance and hours of activity) would outperform environment-only models and predict greater climate-driven extinction risk for this species.

Data

Environmental, land use, and species occurrence data

Environmental data

We obtained climate and elevation data for each occurrence record (below) from the WorldClim database (Hijmans et al., 2005), at a resolution of 2.5 minutes, using the RASTER package (Hijmans, 2019) in R (R Core Team, 2025). Climatic variables consisted of minimum, mean and maximum temperature, precipitation, and 19 derived bioclimatic variables (BIO1 through BIO19) for the present (1960-1990) and the years 2050 and 2070. Future climate data derived from the IPCC5 climate projections from global climate models (GCMs) for three representative concentration pathways (RCPs), which predict climate patterns under different amounts of greenhouse gases emissions (Hijmans et al., 2005). To incorporate uncertainties in estimating future climates, we used data from twelve GCMs (Supporting Information Table A.1) and three RCPs (2.6, 4.5, and 8.5). We cropped all WorldClim layers to South America between 10° and -40° latitude and -80° and -30° longitude.

To eliminate spatially correlated variables, which may lead to increased uncertainty in model parameters and reduced power of the model (de Marco and Nobrega, 2018), we first selected 1,000 random points from the study area using a mask (BIO1) to exclude regions with no data and extracted values of environmental variables for the present (1960-1990) at each point, using the DISMO package (Hijmans et al., 2017). Next, we excluded highly collinear variables regarded as those with variance inflation factors (VIFs) larger than 10 (Dormann et al., 2013), using the USDM package (Naimi et al., 2014). The analysis retained ten variables: altitude, mean diurnal range (BIO2), isothermality (BIO3), mean temperature of wettest quarter (BIO8), mean temperature of driest quarter (BIO9), precipitation of wettest month (BIO13), precipitation of driest month (BIO14), precipitation seasonality (BIO15), precipitation of warmest quarter (BIO18) and precipitation of coldest quarter (BIO19). We used the same set of variables for model projections in the future (2050–2070).

Land use data

We estimated areas of environmental suitability for Micrablepharus atticolus in the present, considering (1) the original extension of the Cerrado and (2) the remnants of Cerrado natural habitats. For future projections, we considered three scenarios for the year 2050: (1) no deforestation (original extension of the Cerrado), (2) if the Brazilian Forest Code (Soares-Filho et al., 2014) is put in practice, and (3) “business as usual,” i.e., without the Forest Code. We obtained the last two scenarios from GLOBIOM-Brazil, an adaptation of the global economic model (GLOBIOM) to project future land use and agriculture production in Brazil up to 2050, based on governmental data and maps on vegetation, remote sensing land cover, crops, livestock and planted forests (Câmara et al., 2015).

Species occurrence data

We obtained 48 non-overlapping geographic distribution records of Micrablepharus atticolus from scientific collections and the literature (Supporting Information Table A.2). To reduce possible sampling biases in the dataset, which might affect model predictions, we applied an environmental filter with the ENVSAMPLE function (Varela et al., 2014), based on the four environmental variables with the smallest VIF values (above): BIO2, BIO15, BIO18, BIO19. This procedure retained 35 occurrence points, which we used in further analyses.

Ecophysiological data

Thermal performance

We obtained the preferred body temperature (Tpref), thermal sensitivity of sprint speed (speed), maximum critical temperature (CTmax), and minimum critical temperature (CTmin) of Micrablepharus atticolus from 90 specimens collected with a combination of pitfall traps and drift fences in four Cerrado localities (Brasília, Distrito Federal; Gaúcha do Norte, Mato Grosso; Lagoa da Confusão, Tocantins; and Nova Xavantina, Mato Grosso). We checked traps daily (often twice a day), took captured lizards to a field laboratory, and placed them in individual plastic terrariums with water and food ad libitum. We obtained Tpref by placing each lizard in a 100 cm x 15 cm x 30 cm MDF (medium density fiberboard) thermal gradient, with an open top and a 2cm–deep substrate of mixed sand and vermiculite. We produced a temperature gradient (20–50 °C) by placing an incandescent light bulb (60 W) at one end, simulating the light and heat of the day, and an ice pack at the other, simulating the temperature of a shaded refuge. Each lizard placed on the thermal gradient received an insulated thermocouple (5SC-TT-T-36-72, Omega Engineering) attached to the ventral region and connected to an Eight Channel Thermocouple USB Data Acquisition Module (TC-08, Omega Engineering), programmed to record the temperature every minute during 65 min. We regarded the mean of these records as Tpref and the upper and lower limits as the maximum voluntary temperature (VTmax) and minimum voluntary temperature (VTmin) of each lizard, respectively (Kubisch et al., 2016).

We obtained locomotor performance data by inducing the animal to run along a 200 cm x 30 cm x 40 cm MDF racetrack at three body temperatures (20, 25, and 30 °C). We raced lizards two hours after completing the measurements of the Tpref experiment, twice at each temperature, with intervals of at least 30 min between each trial. We achieved the required body temperatures by exposing lizards to a heat source (60 W incandescent light bulb) or placing them in a cooler with ice packs and using a single-channel digital thermometer (HH-91, Omega Engineering) to record cloacal temperatures. As soon as lizards reached the desired temperature, we placed them on the track and stimulated them to run as fast as possible by gently touching their tails (Miles, 2004). We recorded each race at 420 fps with a digital camera (Casio HS EX-FH25) located 1.5 m above the center of the track and supported by an aluminum tripod. Using software Tracker 4.80 (https://physlets.org/tracker/), we obtained the maximum running speed of each lizard at each body temperature.

On the day after collecting Tpref and locomotor performance data, we obtained data on critical temperatures. To obtain CTmin, we placed each lizard on a plastic box along with an ice gel pack to lower its body temperature; to obtain CTmax, we used a heat source (60 W incandescent bulb) to increase its body temperature slowly. We observed its behavior and checked every minute its ability to return to the right position when placed supine. Upon losing this capacity, we immediately recorded the lizard’s body temperature with a quick-reading cloacal thermometer (T-6000, Miller & Weber). We collected CTmax data at least one hour after obtaining CTmin. An ethics committee approved all procedures on animal use (Process number: 33786/2016).

Hours of activity

To estimate microclimatic conditions available to lizards in their natural environment, we measured the operative environmental temperature (Te) using PVC (Polyvinyl chloride) models (1.5 cm diameter x 6.0 cm length), colored in gray with spray paint. We placed one probe of a two-channel automatic data logger (Onset U23–003 HOBO Pro v2 2x External Temperature Data Logger) within each model to record the temperature every 10 min. We placed 24 models of each size in different microhabitats used by lizards, with varying degrees of sunlight exposure (e.g., exposed to sun or shade on the ground, under logs, inside termite nests, within the leaf-litter), for at least 20 days. In this way, we recorded the temperature variations to which animals could be exposed in each microhabitat. At each locality, microclimatic data were recorded continuously for at least 20 days during warm-wet season, when environmental conditions are representative of the active season of M. atticolus.

Close to each pitfall trap, we placed another data logger (Onset U23–001 HOBO Pro v2 Temperature/Relative Humidity) that recorded air temperature and relative humidity every 10 min. We protected these data collectors from rain and direct sunlight using a PVC cover and attached them to iron stakes at 50 cm above the ground in a microhabitat that best represented a 6 m radius around each trap. The data obtained from these collectors allowed us to estimate the environmental temperatures (Tenv) experienced by lizards in their environments.

Modeling

Models based on environmental predictors

We built habitat suitability models for M. atticolus using the occurrence records as the response variable and the environmental data as predictor variables, with the BIOMOD2 package (Thuiller et al., 2019). Most statistical models we implemented require a binary response variable, i.e., both species presence and absence records. Because our data consisted of species presence records only, we sampled pseudo-absence (background) points across the study area. We fitted models with eleven different algorithms: three machine learning algorithms (Generalized Boosting Model – GBM, Classification Tree Analysis – CTA and Random Forest – RF) and seven regression algorithms (Generalized Linear Modeling – GLM, Generalized Additive Model – GAM, Artificial Neural Network – ANN, Surface Range Envelope – SRE, Flexible Discriminant Analysis – FDA, Multiple Adaptive Regression Splines – MARS and Maximum Entropy – MAXENT.Phillips and MAXENT.Tsuruoka) (Beaumont et al., 2016; Qiao et al., 2015). When fitting models based on machine learning algorithms, we generated as many pseudo-absence points as our species distribution records; for models based on regression algorithms, we used 10,000 pseudo-absence points (Barbet-Massin et al., 2012). We produced ten replicates of pseudo-absence data in both cases by randomly selecting points outside a 0.025 quantile surface range envelope model based on the M. atticolus presence data. We used 75% of the species presence data to fit models and 25% to assess model accuracy, repeating this process ten times. To evaluate model accuracy, we used the true skill statistic (TSS) (Allouche et al., 2006). Hereafter we refer to these as ENV models.

Models based on ecophysiological constraints

Thermal performance

We produced a thermal performance curve (TPC) (Huey and Stevenson, 1979) using a generalized additive mixed model (GAMM) (Wood, 2017; Zajitschek et al., 2012). The model had maximum running speed as the response variable, body temperature (smooth term based on cubic regression spline) and snout-vent length (covariate) as fixed effects, and individual (lizard’s identification) as a random effect to control for pseudoreplication resulting from multiple runs for each lizard. We used the model to produce thermal performance rasters from the environmental temperature rasters (minimum, mean, and maximum temperature; BIO1–BIO11). We implemented these procedures with the MAPINGUARI package (Caetano et al., 2020).

Hours of activity (Ha)

We used the Te data collected with PVC models to calculate the number of daytime hours animals were between VTmin and VTmax, averaged across all microhabitats sampled, which we regarded as Hact. Next, we built a nonlinear model based on the logistic function to relate Ha (response variable) to Te and used the model to produce a Ha raster from the environmental temperature rasters (minimum, mean, and maximum temperature; BIO1–BIO11).

Modeling

We analyzed the thermal performance and hours of activity data separately. First, we excluded highly collinear variables using VIFs as described above (section 2.1.1 Environmental data). The following were retained: perf_BIO2, perf_BIO3, perf_BIO5, perf_BIO6, perf_BIO8, and perf_BIO9 (thermal performance rasters); ha_Tmax, ha_BIO2, ha_BIO3, ha_BIO4, ha_BIO5, ha_BIO6, ha_BIO7, ha_BIO8, ha_BIO11 (hours of activity rasters). Then, we used the same steps outlined above (section 2.2.1 models based on environmental predictors), using the occurrence records of Micrablepharus atticolus as the response variable and the thermal performance and hours of activity rasters as predictor variables, to produce species distribution models based on ecophysiological constraints. Hereafter we refer to these as PER (thermal performance) and HOU (hours of activity) models.

Hybrid model

Here we grouped the three sets of predictor variables used to build the ENV, PER, and HOU models. The following were retained after removing highly collinear variables using VIFs (section 2.1.1 Environmental data): BIO2, BIO13, BIO14, BIO15, BIO18, BIO19, perf_BIO5, perf_BIO9, ha_Tmax, ha_BIO2, ha_BIO3, ha_BIO4, ha_BIO6, ha_BIO7, ha_BIO8, and ha_BIO11. Next, we used the steps outlined above (section 2.2.1 Environmental models) to produce hybrid models (hereafter HYB).

Model comparison and evaluation

Using the BIOMOD2 package (Thuiller et al., 2019), we obtained TSS values for each combination of four models (ENV, PER, HOU, and HYB), eleven algorithms (GBM, CTA, RF, GLM, GAM, ANN, SRE, FDA, MARS, MAXENT.Phillips, and MAXENT.Tsuruoka), ten pseudo-absence data replications, and ten runs of randomly partitioning the species occurrence data into training (75%) and testing (25%) datasets. Next, we compared model accuracy using a generalized linear mixed model (GLMM) with the LMER4 package (Bates et al., 2015). This model regarded TSS as the response variable and algorithm, pseudo-absence data replication, and run as nested random factors. We used the best linear unbiased predictions (BLUP) to compare the effects of the four modeling strategies on model accuracy (Henderson, 1975; Liu et al., 2008). These BLUPs provide model-based estimates of mean TSS for each modeling strategy that account for variation attributable to algorithms, pseudo-absence replicates, and runs, offering a transparent basis for comparing strategies while respecting the mixed-effects structure of the data.

Model predictions and projections

To decrease uncertainty associated with choices made for each model (Buisson et al., 2010), we built ensemble predictions by selecting models with TSS value ≥ 0.8 and weighing each model’s contribution to the ensemble by its TSS score, using the BIOMOD2 package (Thuiller et al., 2019). Using the original occurrence data and model predictions, we then calculated the optimal threshold and transformed each model’s prediction into a presence-absence raster, with the maxSSS method (Liu et al., 2016) of the BIOMOD2 package (Thuiller et al., 2019) using the PresenceAbsence package (Freeman and Moisen, 2008). We used these data and consequent projection evaluation metrics (confusion matrix, sensitivity, and specificity). Next, we calculated the total area of environmental suitability for M. atticolus for the present and each future scenario by multiplying the number of presence cells by the corresponding cell size. With this data, we calculated the difference between present and future projections, assessing changes of environmental suitability under 72 different scenarios—12 GCMs x 3 RCPs x 2 years (2050, 2070)—for each model. For future projections, we considered three scenarios for the year 2050: (1) no deforestation (original extension of the Cerrado), (2) if the Brazilian Forest Code (Soares-Filho et al., 2014) is put in practice, and (3) “business as usual,” i.e., without the Brazilian Forest Code.

Overall, models based on ecophysiological constraints were more accurate (higher TSS values) than those based on environmental predictors or hybrid models, and the GLMM indicated this difference was significant (Figure 1). The ensemble PER model, which had the highest accuracy, predicted a habitat suitability area of 1,657,000 km2 for Micrablepharus atticolus in the Cerrado, whereas the ensemble ENV model predicted an area of 1,048,866 km2 (Figure 2). Only considering Cerrado remnants in 2010 resulted in an overall decrease of approximately 50% in habitat suitability: 884,256 km2 for the ensemble PER model and 496,558 km2 for the ensemble ENV model. Predicted areas of habitat suitability were consistently more extensive for the PER model (36.7% larger for the entire Cerrado and 43.9% larger only for Cerrado remnants in 2010).

Combining GCMs, RCPs, years, and models resulted in 144 future projections of habitat suitability areas for Micrablepharus atticolus (Supporting Information Figures A.1-A.4), with nearly 100% of projections showing loss of habitat suitability (Supporting Information Tables A.3–A.6). Projections based on the ensemble PER model resulted in higher extinction risk for M. atticolus in all three scenarios (no deforestation, Brazilian Forest Code, business as usual) and, in the no-deforestation scenario, this difference reached approximately 60% (Figure 3). Our results indicate that, regardless of the scenario, habitat suitability areas for M. atticolus will likely undergo a severe reduction soon.

The projections for the twelve GCMs showed significant variation in extinction risk, with percent loss varying between 71.4% and 3.4% relative to the current area and the ENV model (For 2070 in the scenario “ND,” Table A.6). Regarding the RCPs, under a more optimistic prediction (RCP 2.6), the highest extinction risk was 97.0% for the PER model (2070 – BAU, Table A.4) and 89.5% (2070 – BAU, Table A.6) for the ENV model. For the intermediary prediction (RCP 4.5), the highest value was 97.8% for the PER model and 91.5% for ENV model (both 2070 – BAU, Tables A.4, A.6). Under more pessimistic modeling (RCP 8.5), extinction risk would reach 98.5% for the PER model and 91.1% for the ENV model (both 2070 – BAU, Tables A.4,A.6). Thus, it is clear the need to use the largest number of models available to minimize the bias inherent to each GCM and RCP.

The predicted loss of climate suitability areas based on the ensemble PER model was concentrated mainly in the northern, southern, and central regions of the Cerrado (Figure 4). For the ENV model, the loss was regularly distributed throughout the biome (Figure 4). In relation to the three deforestation scenarios (ND, BFC and BAU), the PER model presented mean area loss of ca. 1,522,000 km2, 1,604,000 km2 and 1,619,000 km2 (ND, BFC and BAU respectively), while for the ENV model the average loss was ca. 343,718 km2, 907,255 km2 and 952,422 km2 (ND, BFC and BAU respectively).

We found that incorporating ecophysiological constraints can significantly improve the predictive accuracy of habitat suitability models. Here, “improve” refers specifically to higher TSS-based performance; we did not directly evaluate model transferability or ecological realism under future climates, and complementary evaluation metrics could further strengthen this conclusion. Previous studies that used lizard body temperatures already indicated an improvement in predictive accuracy (Sinervo et al., 2010), and the refinement of these parameters show that this is an excellent way to increase model realism. Incorporating ecophysiological constraints into suitability models in other taxonomic groups has increased confidence in the results (Gamliel et al., 2020). Locomotor performance (sprint speed) is a widely used proxy for fitness in lizards because it is linked to survival and reproductive success via prey capture, predator avoidance, and social interactions (Miles, 2004; Zajitschek et al., 2012), but it cannot fully capture other components of fitness such as reproductive output, the effectiveness of behavioral thermoregulation, or performance in non-locomotor physiological functions.

Our findings do not support the notion that more complex models improve predicting habitat suitability and extinction risk (Fordham et al., 2018). However, this complexity can occur in several paths, and, in the case of this study, the improvement of predictive power was given by the use of ecophysiological predictors (performance and hours of activity). The combination of environmental and ecophysiological predictors was ineffective as in other studies (Ceia-Hasse et al., 2014) because the combined model had the lowest accuracy. This poor performance of the hybrid model is likely related to the structural dependence between the climatic BIO variables and the ecophysiological predictors, because thermal performance and hours of activity are themselves derived from climatic drivers; combining them in the same model may introduce redundancy and inflate the effective dimensionality of the predictor set without adding fully independent information. A possible way to improve the accuracy of the models would be the use of behavioral thermal data obtained in the field (Guo et al., 2020).

Our modeling approach avoids misleading results that omit essential sources of variation (Leach et al., 2016). The improvement in data and analysis has some effect, but it is vital to recognize the limitations of any model mainly because it is challenging to incorporate all the predictors and responses that influence the occurrence of the species (Thuiller et al., 2013). First, our models were calibrated with only 35 occurrence records, a relatively small sample size that may increase uncertainty in spatial predictions despite the use of multiple algorithms, pseudo-absence replicates, and ensemble approaches. Nevertheless, because all modeling frameworks were fitted to the same dataset and we focus on relative differences among them, we expect our comparative conclusions to be robust. Second, our estimates of hours of activity are based on microclimatic measurements collected over a relatively short 20-day window at each site, which may not fully capture seasonal variation in thermal environments and activity opportunities and therefore adds uncertainty to the performance of the HOU predictor relative to the PER and ENV models. Third, using multiple algorithms implemented in BIOMOD2 with a relatively small calibration dataset (35 occurrences) may increase variability among individual models and the potential for overfitting, but our TSS-based weighting and ensemble forecasting were designed to down-weight poorly performing algorithms and to provide more robust predictions. Fourth, because we evaluated accuracy exclusively with the threshold-dependent True Skill Statistic (TSS), which can be influenced by prevalence and by the lack of full independence among repeated data partitions and pseudo-absence replicates, our performance estimates should be interpreted with some caution, although the same evaluation framework was applied consistently across all models.

Another critical limitation of models using ecophysiological constraints is that species responses to thermal gradients can vary in time and space (Gamliel et al., 2020). Using better models, we can more accurately predict the future extinction risk of the species analyzed. To Micrablepharus atticolus, this risk is worrisome since the species is predicted to lose more than 90% of its area of habitat suitability in some scenarios. The results pointed to a more significant loss of suitability for models based on ecophysiological constraints. Work with mountain salamanders comparing correlative and mechanistic models designated larger areas suitable for the latter approach (Lyons and Kozak, 2020). This suggests that species responses may be idiosyncratic and that different modeling approaches should ideally be evaluated. Indeed, several comparative studies have shown that correlative BIO-based models can match or even outperform mechanistic approaches in some systems when the chosen physiological proxy does not fully capture the dominant constraints on species distributions (e.g., Robertson et al., 2003). Thus, our results should be interpreted as evidence that ecophysiology-based models can outperform environment-only models for M. atticolus, a small heliothermic lizard whose fitness is tightly constrained by thermal performance, and perhaps for species with similar ecological and physiological traits, rather than as a universal expectation that mechanistic models will always be superior across taxa and contexts. In our analysis, besides eliminating the leading causes of uncertainties associated with modeling by incorporating many GCMs (Buisson et al., 2010; Lyons and Kozak, 2020), we assessed which iteration of environmental, ecophysiological, or a combination of both achieved the highest TSS-based accuracy, which gave us greater confidence in our comparative results, while acknowledging that this does not directly address transferability under novel conditions.

The use of land use masks with the deforestation pattern for the future adds further realism to the predictions. We note that our estimate of current suitability in the hybrid models combines baseline climate data from 1960–1990 with land-use layers representing more recent Cerrado remnants, which introduces a temporal mismatch. This may slightly bias absolute estimates of present suitability and extinction risk, but because all modeling frameworks were calibrated against the same baseline climate and our main inferences rely on relative differences among models and scenarios, we do not expect this mismatch to qualitatively affect our conclusions. Models that do not incorporate this parameter can erroneously predict a lower risk of extinction considering the adequacy in deforested areas. Our predictions show a loss of suitable climatic conditions, mainly in the northern, southern, and central portions of the species’ current distribution in the Cerrado (Figure 4), reducing the area drastically by 2050 and even further by 2070, suggesting the extirpation of many of the current populations.

Despite the difficulties in obtaining ecophysiological data, a sampling of only four localities was enough to produce a more accurate model than using only environmental data. These four sites lie within the core region of the Cerrado and do not encompass the full geographic, climatic, and vegetational variability of the biome, particularly toward southern regions and peripheral transitions with neighboring biomes. As a result, extrapolating these physiological parameters across the entire Cerrado may overlook fine-scale intraspecific variation and adds uncertainty to our extinction risk estimates. So, an earlier cost-benefit analysis must be made before deciding which model to adopt. Models that consider the characteristics of the species have the disadvantage of not being able to project suitability (Willis et al., 2015). In this work, we can make these projections and thus increase the applicability of these models to conservation planning. However, applying these models to other animal species and communities is necessary to test how general these patterns are. In particular, we expect our conclusions to transfer most directly to small, diurnal, ground-dwelling lizards in open tropical savannas that share similar thermal ecology with M. atticolus, whereas extensions to other functional groups should be made with caution.

We assessed the impacts of climate change on the thermal physiology of Micrablepharus atticolus and predicted the location of climatically suitable habitats. The performance curve analysis indicated a drastic decrease in this parameter above Tbs of 34 °C (Figure A.5); this indicates that these lizards are very sensitive to high temperatures. However, lizards have strategies to avoid overheating. One of the primary mechanisms by which lizards can thermoregulate is modifying their activity period (Huey et al., 1977). An increase in hours of restriction (Hr) would mean individuals of M. atticolus spend more extended periods in refuges, and therefore, the net energy gain would become insufficient for growth, survival, or reproduction (Vicenzi et al., 2017). In this way, we conclude that M. atticolus is vulnerable to climate change but incorporating these parameters into climatic adequacy models is scarce, making it challenging to apply this logic. More broadly, our results highlight the potential for ecophysiological constraints to strongly modulate extinction risk estimates, but they also underscore that such inferences are intrinsically species-specific and contingent on life history, habitat use, and thermal strategies. Future work applying comparable mechanistic approaches to other taxa will be essential to evaluate how general these patterns are across ectothermic vertebrates.