Vulnerability to climate change of a microendemic lizard species from the central Andes

Given the rapid loss of biodiversity as consequence of climate change, greater knowledge of ecophysiological and natural history traits are crucial to determine which environmental factors induce stress and drive the decline of threatened species. Liolaemus montanezi (Liolaemidae), a xeric-adapted lizard occurring only in a small geographic range in west-central Argentina, constitutes an excellent model for studies on the threats of climate change on such microendemic species. We describe field data on activity patterns, use of microhabitat, behavioral thermoregulation, and physiology to produce species distribution models (SDMs) based on climate and ecophysiological data. Liolaemus montanezi inhabits a thermally harsh environment which remarkably impacts their activity and thermoregulation. The species shows a daily bimodal pattern of activity and mostly occupies shaded microenvironments. Although the individuals thermoregulate at body temperatures below their thermal preference they avoid high-temperature microenvironments probably to avoid overheating. The population currently persists because of the important role of the habitat physiognomy and not because of niche tracking, seemingly prevented by major rivers that form boundaries of their geographic range. We found evidence of habitat opportunities in the current range and adjacent areas that will likely remain suitable to the year 2070, reinforcing the relevance of the river floodplain for the species’ avoidance of extinction.


Scientific Reports
| (2021) 11:11653 | https://doi.org/10.1038/s41598-021-91058-w www.nature.com/scientificreports/ accentuated thermal and hydric constraint imposed by the environment as they challenge the persistence of the current L. montanezi population. We hypothesize that the long-term persistence of L. montanezi could be threatened by climate change considering its microendemic condition and the present thermal, humidity, and dispersion restraints of its natural habitat. Herein, we analyze its current thermal niche, use of microhabitat, and activity patterns to predict L. montanezi warming tolerance, thermal safety margins, and the potential hours of restriction of activity. Using ecophysiological, geographic, and climatic parameters in Species Distribution Models (SDM) we predict probabilities of persistence, dispersion, or extirpation of this narrowly distributed population under global warming 2050 and 2070 scenarios 62 . Finally, we propose strategies for present and future conservation of L. montanezi.

Results
Relationships of T b , T pref, VT max , CT max , CT min with body size (SVL), body mass (BM) and body condition index (M i ). T b and T pref of L. montanezi were not related with SVL, BM or M i (Supplementary   Table S1). However, CT max was negatively related to BM and SVL, and VT max was positively related to M i (Supplementary Table S1).
Determination of the main thermal sources for thermoregulation, and relationship among T b and climatic factors: wind velocity (Wd) and relative humidity (RH%). T b s of L. montanezi were lower than T s , while T b and T s were significantly higher than T a (Supplementary Table S2). T b s increased with T a s, and did not depend on T s s (Stepwise Regression, F Tb-Ta (18) = 11.98, P < 0.001; F Tb-Ts (18) = 8.95, P > 0.3, Fig. 1).
Operative temperatures and thermal quality of the habitat. Overall, the T e s were different between models exposed to sun or shad conditions (Two-way ANOVA, T e -sun vs. T e -shade, F (1) = 786.4, P < 0.001), with the T e -sun higher than T e -shade (Supplementary Table S3, Fig. 2). Moreover, the hourly mean T e -sun was higher than the mean T e -shade (Two-way ANOVA, T e -sun vs. T e -shade * Hour F (8) = 12.5, P < 0.001, Supplementary Table S3, Fig. 2).
Only 10% (N = 2) of the T e s records were included within T set , whereas 65% of T e s (N = 13) were lower than T set , and 25% (N = 5) of T e exceeded the T set . However, the 38% (N = 8) of the operative temperatures in open habitat (T e -sun) exceeded the T set from 11:00 to 15:00 h, whereas the operative temperatures from shaded habitat (T e -shade) were lower than the T set of the species. Also, T e exceeded CT max during two hours in sunny habitats (Fig. 2).
Thermal quality (d e ) was highly variable throughout the L. montanezi activity period, showing a trimodal distribution in open habitat (d e -sun), and a bimodal pattern in shaded habitats (d e -shade; Fig. 2). The mean   Fig. 2). A majority (65%, N = 13) of the T b s of L. montanezi were lower than the minimum T set , whereas 10% (N = 2) of the T b s were higher than the maximum T set . The remaining 25% of the T b ' records (N = 5) were within the T set range. Males, females and juveniles of this species exhibit negative values of E, indicating there are other constraints that prevent them from thermoregulating within their set-points of T pref (Table 1). Daily activity and microhabitat use: bare soil exposed to sun (BS-sun), bare soil in shade (BS-shade), and rocks in shade (WR-shade). Liolaemus montanezi (N = 21) was active for 9 h (09:00 to 18:00 h). Mean activity, [considered as the direction (vector) of the "Ɵ" angle and dispersion (± circular standard deviation) of the activity which is stretching within the origin to the end of the last observation] was at 12:58 h (Circular ± s.d. = 2 h 33 min = 38.261°, Fig. 3). Overall, the observed frequencies of individuals depicted a bimodal pattern of activity, with peaks occurring between 10:00 to 12:00 h in which 33% (N = 7) occurred and between 15:00 to 16:00 h, when 14% (N = 3) of lizards were seen during these time intervals (Fig. 3). Nearly 62% (N = 13) of the observed lizards were recorded using BS-shade, 14% (N = 3) occupying BS-sun, and 24% (N = 5) associated to WR-shade microhabitats (Supplementary Table S4). Observed frequencies were different than expected among the microhabitats (Chi-square χ 2 Test 2; 21 = 8.0, P > 0.01), and the deviation of the expected frequencies indicated that BS-shade was more frequently used than BS-sun or WR-shade (Supplementary Fig. S1). However, when combined frequencies of temporal and spatial activity, the observed frequencies were marginally significantly higher than expected in the BS-sun at the 10:00 -11:00 h, and WR-shade at the 17:00 -18:00 h (Pearson Chi-Square χ 2 Test (16) = 29.47, P < 0.02; Multiple Comparison, z BS-sun; 10-11 h = 2.66, P < 0.005; z WR-shade;17-18 h , = 9.08, P < 0.002, Supplementary Table S4), suggesting that lizards observed at these hours preferably occupies BS-sun and WR-shade microhabitats, respectively.   (Fig. 2).

Pairwise comparisons between T b and T pref , and comparisons between T b , T pref , CT
The h a and h r as a function of T a, max of L. montanezi were explained by the logistic Richard's curve (Table 2;     Table S7) to predict the current known distribution of L. montanezi. However, SDMs based on GAM failed for ecophysiological predictors, as well as some based on ANN (Supplementary  Table S7). For projections, most SDMs indicated the locality where L. montanezi exists (Supplementary Figs. S2, S3). However, the ensemble derived from the RF algorithm was considered the most accurate model as it has fewer over-projections that describe the contemporary and future locations of habitat suitability for L. montanezi (Figs. 5,6). Overall, the occupancy likelihood and habitat suitability modeled based on RF shows that the probable potential area with suitable climatic properties for L. montanezi represents a narrow stretch located toward the south of the actual distribution in the junction of Blanco and La Palca rivers (San Juan, Argentina) by 2050 and 2070, under 4.5 and 8.5 RCP (Figs. 5, 6).

Discussion
Given the rapid decline and local extirpations of many lizards populations as consequence of climate change 13 , it has been emphasized that natural history traits and physiology are crucial to determine what critical environmental factors might be stressing and driving the decline in threatened species 63 . Liolaemus montanezi is a taxon that thermoregulates via a complex use of heat sources. Individuals were found basking directly exposed to the sun during the morning and retreat to shaded patches of Bulnesia retama and Larrea divaricata bushes to avoid overheating before midday during late spring. The mean T s was significantly higher than T b and T a and, T b s only increased according to the air temperature, suggesting the usage of air temperature on thermoregulation 64,65 . Cooling by convection is common in some psammophilous species of Liolaemus inhabiting sandy environments from Brazil (L. occipitalis 66 76 ). Indeed, downslope winds with strong convection are common phenomena observed along the eastern slopes at medium or high elevations in the Andes, especially from midday and thenafter 77 . In the case of L. montanezi, the population is located in alluvial fans in the western margin and riverbed of the Blanco river, an open gully with steep slopes that results in a wind corridor. The convection of heat, rather than radiation or conduction, appears to be the main mechanism of heat loss that explains the T b of L. montanezi (Fig. 1). Table 2. Values of the parameters of the first half of Richard's curve fitted for the hours of activity (h a ) and hours of restriction (h r ) on air maximum temperature (T a, max ) of L. montanezi. A (asymptote), k (rate), i (inflection point) and m (shape) using nonlinear least squares ('nls' function R-stats) with 'modpar' and 'SSposnegRichards' of the R-package 'FlexParamCurve' . The 95% c.i. were determined with the function 'confint2' of the R-package 'nlstools' 108 . The m parameter is estimated from the data and fixed to allow others to behave as free parameters.   www.nature.com/scientificreports/ In addition, body temperatures of L. montanezi varied according to the moisture level of the habitat. The inverse relationship between T b and the percentage of relative humidity (RH %) suggests that the loss of environmental moisture at midday (12:00 h) might intensify the restriction of the activity, particularly in the sun-exposed microhabitats. Thus, when humidity decayed to near zero, lizards were frequently sighted using shaded habitats even when they still have high T b s in the shade (Figs. 1, 2). Sheltering under the bushes probably helps them to avoid dehydration which could also affect the accuracy of thermoregulation 78 . Furthermore, Nicholson et al. 79 studying populations of Anolis stratulus and Ameiva exsul had pointed out that in water-restricted locations the interaction of temperature and humidity has strong influence on the daily activity patterns. Such environmental conditions have a noticeable demographic impact in many lizards since juveniles are more vulnerable to water loss than adults. Nevertheless, we found there were no differences in the ways in which males versus females, or juvenile versus adult lizards interact with their environment. The T b of L. montanezi were not different between sexes and age classes, nor did T b depend on body size or body condition. But, the CT max was lower in lizards of larger body size, leading to the assumption of a more accurate thermoregulation in adults. Thus, we hypothesize that the drastic decline of humidity at 12:00 h might entice lizards to refuge in more suitable habitats that provide hiding shelters in shaded sites. During these restriction hours, lizards will wait until milder heat and humidity conditions are given to resume their daily active behavior.
These activity patterns usually reflect the extent in which the behavioral thermoregulation buffers the environmental restrictions to maintain the T b to optimize physiological processes 80 . Thus, changes in ambient temperature during the day could cause the frequency of active lizards to describe a distribution that is unimodal or bimodal. Likewise, these patterns of behavior might also vary among the seasons of the year 80 . Hence, during spring L. montanezi might be active during 9 h from 09:00 to 18:00 h with activity peaks occurred between 10:00 to 12:00 h and between 15:00 to 16:00 h showing a bimodal pattern of activity (Figs. 2, 3) with a retreat during the time when the T e max was higher than 44 °C (Fig. 2). When considering the operative temperatures of the models deployed over open habitats [i.e., those with high temperatures (T e -sun) or those placed in shaded habitats (T e -shade)], we found that they also support a pattern of activity likely associated with the thermal quality (d e = d e -sun + d e -shade) of the environment (Fig. 2). Thus, we predict that the daily activity of these lizards starts when the environmental thermal quality is low (d e -sun and d e -shade exhibit the highest values) and increases while the thermal quality rises (d e -sun and d e -shade decrease). However, the number of active lizards outside shelters progressively decreases during the hotter hours and they are forced to retreat into shaded microhabitats where the thermal quality is optimal (d e -shade tending to zero). Based on our observations, we predict that a second peak of activity occurs in the afternoon when the thermal quality in shaded habitats decreases while it improves in open habitats allows the lizards to resume active behaviors (Fig. 2). Therefore, when considered as trimodal and bimodal patterns of d e -sun and d e -shade, respectively, both describe scenarios in which lizards emerge in the morning, retreat before the midday, and emerge again at the afternoon. This pattern of activity of L. montanezi during spring support the h a and h r observed considering the model for global warming (Fig. 4, Supplementary Table S5). Our operative temperature models likely provide more accurate representations of available temperatures at the capture than measurements of ambient air temperature, they do not represent "instantaneous" measures of equilibrium temperatures as recommended by Bakken et al 81 and therefore our analyses involving operative temperatures should be interpreted as such.
Microhabitat selection by individuals may allow populations to maximize the geographic range 82 . In this sense, shrubs play a significant role in thermal ecology of the species. Liolaemus montanezi maintain the T b near the operative temperatures recorded in the shaded habitats (Fig. 2), behaving as a thermoconformers during the hottest hours of the day. Thus, in a well-identified habitat, most of lizards were seen using sandy, bare soil in shaded places, sheltered beneath shrubs. However, when combining such microhabitat use with their diel activity pattern, we predicted that lizards would occupy strictly sunny, bare soil in the morning and would perch on rocks only in the afternoon (Supplementary Table S4). The fact that the individuals might use sunny bare soil at the beginning of the activity period, shuttling to cooler microhabitats during the warmest hours suggests that the species could exhibit spatial and temporal flexibility in microhabitat use 83 but their thermoconformity in the shady refuges point out their current vulnerability.
Preferred temperature of L. montanezi during spring was higher than T b as shown for other Liolaemus species 84 , even when the mean T e was closer than T set (Fig. 2). These results suggest that the species could obtain temperatures close to T pref in their natural environment with low costs of thermoregulation. However, the effectiveness of thermoregulation (E = − 0.29) shows that L. montanezi behaves almost as a thermoconformer (E = 0). But, the negative value shows that this species is avoiding thermally optimal microenvironments probably to minimize the risk of overheating 85 and dehydration (during the warmest hours of the day) selecting among the coolest and shaded microhabitats. At least a third (33%) of the recorded T e values were higher than the T pref during the warmest hours which increases the risk for these lizards of attaining near-lethal body temperatures close to their CT max . Accordingly, we found that their thermal safety margin (TSM = 3.57 °C) and warming tolerance index (WT = 8.8) were notably lower than reported in others liolaemid lizards from cold temperate environments such as Phymaturus tenebrosus (TSM = 15.58 and WT = 23.3 86 ) and P. calcogaster (TSM = 3.98 °C and WT = 10.67 °C 87 ). Therefore, these observations suggest that L. montanezi has a narrow TSM range to prevent it from reaching lethal temperatures 45 . However, the thermal environments in the sub-Andean region characterized by high thermal variation could have a dramatic influence on such thermal tolerance. Thus, the mean breadth of thermal tolerance (CT max − CT min , sensu Litmer & Murray 88 ) for L. montanezi was 31.68 °C showing that this species is likely eurythermic.
The east side of the Andes, inhabited by L. montanezi, is expected to experience thermal alteration during global warming. An increase of T a and a reduction of rainfall are expected to impose droughts and habitat modification 61,89 . Added thermal stress in the habitat is expected to cause a decline in the water availability and reduce plant growth and recruitment 61 13 , which might compromise population persistence and cause local extirpations and threaten extinction 13,50 in particular in microendemic taxa. The current known distribution of Liolaemus montanezi is explained by abiotic and ecophysiological factors (Supplementary Table S6). As we pointed out above, the small area occupied by this species is confined between rivers to the southeast and steep slopes to the west and northwest. Because the known population of L. montanezi occupies sandy patches interspersed between alluvial fans in the western margin and riverbed of Blanco river, the availability of moisture as indicated by evapotranspiration (one of the predictors in the ecophysiological models) may play a main role in L. montanezi habitat physiognomy. The role of evapotranspiration is more conspicuous in the model than h r and h a for L. montanezi (Supplementary Table S6). Contemporary bioclimatic and ecophysiological models have produced comparable results. Hence, the SDM projections show a narrow range of suitable habitat for L. montanezi located along their current riverine habitats (La Palca and Blanco Rivers) southeast of their current range (Figs. 5, 6). Future SDMs predictions based on RF (Supplementary Table S7), includes the current range of L. montanezi with a higher probability of persistence by 2070 under RCP 8.5 (Figs. 5, 6).
The SDM reinforces the hypothesis of the river floodplain south of their current range serving as the more suitable habitat for the species in the future. However, numerous surveys were done in the predicted area of occupation and confirm that L. montanezi is not currently occupying these southern locations. The predicted pattern of habitat suitability is probably influenced by a fluvial effect of moisture supply which promotes the recruitment and persistence of the dominant plants used by L. montanezi. One alternative interpretation is that the model forecasts increased rainfall as projected between 2081-2100 by Barros et al. 61 for that area. The MaxEnt, GAM, GLM, and ANN SDMs showed unreliable projections by overestimating the present and future scenarios localizing the species in currently unoccupied environments (Supplementary Figs. S2, S3). Instead, the RF approach projected a low likelihood of occupancy in neighboring areas in which the reduced range of this species will be intensified by adverse changes to its current habitats. Therefore, the microendemic character plus the strong dependence of this species on sandy patches covered with shrubs of B. retama and L. divaricata, together with the actual habitat fenced by the rivers suggest high constraints on niche tracking. This study also suggests the vulnerability of other species present in the area such as L. eleodori, L. parvus, and Phymaturus punae. The assessment of population vulnerability to climate change is a powerful tool to highlight how this and other species could become threatened within a relatively short time (~ 30 or 60 years). Our projections should be interpreted with caution given our low physiological sampling which were only held at the beginning of the activity season in spring. Despite of our modelled T e obtained for 4123 days (see Supplementary material for details) the habitat suitability projections should be seen as preliminary. Additional work, with deeper sampling across more of the year, is necessary to fully understand the potential vulnerability to climate change in this species. This approach, although employing low sample size, provides enough species-specific ecophysiology and vulnerability information to prompt immediate reassessment of the conservation status of L. montanezi. Moreover, the SDM predictions can play an important role in support of a habitat conservation initiative. Although we cannot discard the notion of natural dispersal across the rivers to the south resulting in colonization of areas projected to be suitable in the future, prudent conservation would target new surveys of these areas to search for yet-undetected populations and to assess habitat quality for active translocation plans.

Materials and methods
Field work. Lizard sampling, study area, and climate. We captured a total of 21 specimens (9 males, 7 females, 5 juveniles) during one day in late spring (November) by noose in sandy environments in the riverbed at the NW margin of the junction of the Blanco and La Palca rivers, Iglesia Department, San Juan Province, Argentina (− 29.55 S; − 69.19 W, 2185 m asl). The study area belongs to the Monte phytogeographic region dominated by xerophilous plants such as shrubs of Bulnesia retama, Larrea divaricata, Prosopis alpataco, and Atriplex crenatifolia 91 . The climate corresponds to arid cold desert (BWk) 92 with rainfall occurring mainly in the summer season and the annual mean temperature is < 18 °C.
Field data. Body temperatures (T b ) were taken in active lizards using a catheter probe TP-K01 (1.62 mm diameter) introduced ca. 3 mm into the cloaca. Individuals were handled by the head to avoid heat transfer and the temperature was recorded within 20 s of handling. In order to evaluate the main heating resources used by L. montanezi, microenvironmental temperatures were recorded at the exact site of capture for each lizard: substratum temperatures sand, rocks or beneath dwarf shrubs (T s, TP-K03 substrate probe), and air temperature at 1 cm above the ground (T a , TP-K02 gas probe). Probes were connected to a TES 1302 thermometer (TES, Electrical Electronic corp., Taipei, Taiwan, ± 0.01 °C). Thereafter, to determinate the relationship among environmental variables and lizard's activity, we measured the wind velocity using an anemometer (Proster, ± 0.1 m/s) taken at 1 cm above the substrate at the exact first sighting of the lizard and we measured relative humidity (RH %; HOBO, Pro V2, HR%/T°C ± 2%). Snout-vent length (SVL, Mitutoyo, type Vernier digital caliper ± 0.01 mm) and body mass (BW, 10 g Pesola, spring scale ± 0.5 g) were also registered.
To determine the available spatiotemporal heterogeneity of the microenvironmental temperatures for thermoregulation, operative temperatures (T e , sensu Bakken 93 ) were obtained using nine polyvinyl chloride (PVC) pipe models connected to external dataloggers (HOBO, U23-003, T°/T°C ± 1 °C) during the collection of lizards. The model with the greatest correspondence of thermal data to reflect the lizard T b was a PVC pipe (80 mm length × 2.15 mm thickness) sealed at both ends with silicone Fastix (Regression: Adjusted R 2 = 0.846, N = 2836, slope = 1.09, confidence interval = 1.05 -1.14) 94,95 . Subsequently, the models were deployed in the study area in the most representative microhabitats used by L. montanezi: (i) Bare soil with sunny, sandy substrate, (ii) Bare soil with shady, sandy substrate (beneath shrubs), and (iii) Weathered rocks in the shade beneath shrubs www.nature.com/scientificreports/ (3 microhabitats × 3 replicates). The T e was recorded every 5 min during the lizard activity from 8:00 to 19:00 h during one day of fieldwork.
Activity pattern and microhabitat use. Daily activity and microhabitat use were recorded during one day throughout randomly visual surveys in ~ 2 km 2 area. The survey was done walking from 08:00 to 19:00 h at a very slow pace (5-6 m/min) to provide enough time to scan all available habitat. Then, we registered the number of active lizards (frequencies of lizards) and time of day of each lizard sighting. All lizards seen were captured to avoid any possibility of temporal pseudoreplication. Due to the complexity of the habitat structure, the microhabitat used by the individuals were registered according to three main categories also used to deploy the T e models, unifying all the suitable records as follows: (i) BS-sun or bare soil with sandy substrate exposed to sun when the lizards can be seen moving outside shelters; (ii) BS-shade or bare soil with a sandy substrate in the shade when the lizards can be seen perching in bushy-edges or sheltering beneath shrubs, and (iii) WR-shade or weathered rocks in the shade when the lizards were seen on shaded rocky substrate or on rocks in the shade.

Estimation of preferred body temperatures (T pref ), and effectiveness of thermoregulation (E).
Lizards were provided with a thermal gradient produced by a 100-W incandescent lamp from 15° to 50 °C (taken 1 cm above the substrate) to thermoregulate during the experiments. Body temperatures were taken using ultra-thin (1 mm) catheter thermocouples located ca. 5 mm inside the cloaca and taped at the base of the lizard's tail to prevent the thermocouple from being dislodged during the experiment. The temperature of each lizard was obtained every 5 min for 3 consecutive hours by connecting the thermocouple to an 8-channel data-logger (Measurement Computing 1.2 kHz Data Acquisition Device, OMEGA, TC-08 ± 0.5 °C, Stamford, CT, USA).
The preferred body temperatures registered for each individual during 3 h trial (T pref-i ) were used to obtain the mean and the interquartile of T pref-i for each individual (T pref and T set-i , respectively). The accuracy of thermoregulation (d b ) of L. montanezi in their natural environment, was calculated as the mean of the absolute values obtained from the deviations of T b-i from T set-i (individual deviation; d b-i ). The index of the average thermal quality of the habitat from the organism's perspective (d e ) was calculated as the mean of the absolute values of the deviations of mean T e (obtained within the hour of the capture for each lizard) from the T set-i . The effectiveness of temperature regulation, E, was calculated as 1 − (mean d b-i /mean d e-i ) 38 . The values of E range from − 1 to 1, and E ~ 0 represents thermoconformers, E ~ 0.5 moderate thermoregulators, and E ~ 1 effective thermoregulators. Negative E values occur when lizards avoid thermally high-quality habitats with T e near or within the range of T pref 85,97 . We also calculated the d e in the shade (d e -shade) and d e in the sun (d e -sun) in order to determine their differences in thermal quality, and to examine their contribution to the lizard's activity.
Thermal tolerance: Determination of critical thermal minimum and maximum, and voluntary thermal maximum. Critical thermal minimum (CT min ) and critical thermal maximum (CT max ) were determined by means of cooling and heating trails, respectively. For CT min and CT max , lizards were placed individually in a plastic transparent box (20 cm × 20 cm × 20 cm) with a layer (5 mm) of high-density Styrofoam covering the bottom of the box to prevent thermal conductance. Lizards were connected to ultra-thin K-type thermocouples (OMEGA, 5SC-TT-K-30-36; diameter = 0.076 mm, Norwalk, Connecticut, USA) introduced ca. ~ 5 mm inside the cloaca, and taped at the base of tail to prevent the thermocouple from being dislodged during the experiments, and T b was recorded every 5 s using a data-logger (Measurement Computing 1.2 kHz Data Acquisition Device, OMEGA, TC-08 ± 0.5 °C, Stamford, CT, USA).
The lizards were placed individually in a 2 °C chest refrigerator with glass-top door cooled at constant rate (approximately, − 0.7 °C/30 s). During cooling, the lizards were turned onto their back (no more than four times per individual) until they reached CT min , considered as the T b at which the lizard was no longer able to right itself when placed on its back 86 .
After at least 48 h, the heating experiments were performed. Lizards were heated at a constant rate (0.7 °C/10 s) using a 150-W infrared lamp placed 30 cm overhead. In the same experimental system, the voluntary thermal maximum (VT max ) defined as the body temperature that induces a behavioral response seeking to cool down, was recorded 98 . Each lizard was monitored throughout the trials to record the critical thermal maximum (CT max ), defined as the temperature at the higher extreme of tolerance in which the lizard was no longer able to right itself when placed onto its back 27,86,99 . Finally, after the determination of CT min , VT max and CT max , each lizard was promptly removed from cooling or heating systems and steadily returned to room temperature to prevent the death of the individuals under study.
Vulnerability to global warming. To estimate population vulnerability to global warming we calculated the Warming Tolerance index (WT), defined as the difference between mean CT max and mean T e to show how close lizards are environmental temperatures to detrimental or lethal temperatures for physiological processes 86,100 . We also calculated the thermal safety margin (TSM) as the difference between mean T pref max (upper T pref ) and mean T e (sensu Deutsch et al. 45 ). Since most physiological processes are considered to occur at T pref , the upper T pref Scientific Reports | (2021) 11:11653 | https://doi.org/10.1038/s41598-021-91058-w www.nature.com/scientificreports/ threshold is a good proxy of the upper optimal temperature (T opt ) for most physiological functions (i.e., digestion, locomotion or embryo development) that could exhibit a large degree of plasticity within the population 89 .

Species distribution models (SDMs) using bioclimatic and ecophysiological predictors.
We estimated the present and future habitat suitability derived from the known distribution locations for L. montanezi and its physiological parameters described above. For this purpose, we implemented two approaches for SDMs, first using generic bioclimatic method and secondly ecophysiological (species specific) predictors. The latter uses georeferenced grids (rasters) derived from predicted behavioral response of our focal species to environmental parameters and primary productivity estimates associated with the species' habitats and community. We describe the SDM methods and their corresponding variable derivations in the Supplementary Material S1.
Statistical analyses. Variability in thermophysiological variables was described using descriptive statistics (mean ± standard deviation, minimum and maximum). Normality and variance homogeneity assumptions were tested using Kolmogorov-Smirnov's test and Levene's test, respectively. When normality or variance homogeneity assumptions were not met, non-parametric correlation, Mann-Whitney, and Kruskal-Wallis rank sum were used 101  Daily activity patterns were described using circular statistics. We estimate the mean vector (hour) and standard deviation among the frequencies of individuals with Oriana 4.02 102 . We also summarized all records plotting the frequencies of lizards sighted per hour using an angular histogram with Oriana 4.02 102 . The use of microhabitats was analyzed by contrasting the observed versus expected frequencies of lizards using a Chi-square test (χ 2 ) under the assumption of uniformity in the expected frequencies (H 0 = all microhabitat categories are equally used by lizards) 101 . Then, if the differences between observed versus expected frequencies were significant, a diagram of the deviation of the occurrence of individuals from its expected frequency (− 1 to 1 scale) was drawn to represent the relative changes of individuals among the microhabitats 103 .
A third approach that combined the periods of activity (by time of day) and microhabitat use (by categories of microhabitat) allowed the contrast of these two types of lizard frequencies. We then analyzed a 9 × 3 contingency table containing the frequencies of lizards by time intervals (N = 9, from 09:00 to 18:00 h) and microhabitat categories (N = 3, BS-sun, BS-shade, and WR-shade). Pearson's Chi-square test (χ 2 ) was performed to evaluate significant differences between observed versus expected frequencies under assumption of uniformity in the expected frequencies (H 0 = all microhabitat categories are equally used in any hour by the lizards). Post hoc analysis was performed from the adjusted residuals to obtain z-scores for each combined case. Subsequently, P-values were calculated from the transformation of Chi-square values derived from multiple tests 104,105 . Correction of Bonferroni was applied dividing the standard significance level (α = 0.05) by the number of multiple tests (N = 27, derived from 9 × 3 contingency table). Thus, the adjusted significance P-level obtained for the whole model was P < 0.0018.
Snout-vent length (SVL) and body mass (BW) were included in the scaled mass index of body condition in each individual (M i , sensu Peig and Green 106 ) to determine the scaled mass index of condition as an indicator of the health or quality assumed to be related to fitness. The scale mass index was calculated as: where M i and SVL i are the BW and SVL of each individual, SVL 0 is the arithmetic mean SVL of the population, and bSMA is the standardized major axis slope from the regression of ln(BW) on ln(SVL) for the population 106 . The scaling bSMA exponent was calculated directly using the software RMA v. 1.21 107 . We then tested for the influence of M i on thermal traits (T b , T pref , CT max , VT max , CT min ) according to sex and age classes (males, females, and juveniles).

Ethical statement.
Capturing and handling of the lizards were conducted in accordance with international standards on animal welfare (ASIH/HL/SSAR Guidelines for Use of Live Amphibians and Reptiles), being compliant with Argentinian regulations (Argentinean National Law #14.346). All individuals were collected under permits Exp. Number 1300-2643. Field and laboratory protocols were approved in the UNSJ-SEADS-2017-RBSG research engagement and in the CICITCA-UNSJ 21/E1101 plan (Universidad Nacional de San Juan).

Data availability
All data needed to produce the results and discussion in the paper are present in the paper and/or the Supplementary Material.