Alpine viper in changing climate: thermal ecology and prospects of a cold-adapted reptile in the warming Mediterranean

Abstract

In a rapidly changing thermal environment, reptiles are primarily dependent on in situ adaptation because of their limited ability to disperse and the restricted opportunity to shift their ranges. However, the rapid pace of climate change may surpass these adaptation capabilities or elevate energy expenditures. Therefore, understanding the variability in thermal traits at both individual and population scales is crucial, offering insights into reptiles' vulnerability to climate change. We studied the thermal ecology of the endangered Greek meadow viper (Vipera graeca), an endemic venomous snake of fragmented alpine-subalpine meadows above 1600 m of the Pindos mountain range in Greece and Albania, to assess its susceptibility to anticipated changes in the alpine thermal environment. We measured preferred body temperature in artificial thermal gradient, field body temperatures of 74 individuals in five populations encompassing the entire geographic range of the species, and collected data on the available of temperatures for thermoregulation. We found that the preferred body temperature (T_p) differed only between the northernmost and the southernmost populations and increased with female body size but did not depend on sex or the gravidity status of females. T_p increased with latitude but was unaffected by the phylogenetic position of the populations. We also found high accuracy of thermoregulation in V. graeca populations and variation in the thermal quality of habitats throughout the range. The overall effectiveness of thermoregulation was high, indicating that V. graeca successfully achieves its target temperatures and exploits the thermal landscape. Current climatic conditions limit the activity period by an estimated 1278 h per year, which is expected to increase considerably under future climate scenarios. Restricted time available for thermoregulation, foraging and reproduction will represent a serious threat to the fitness of individuals and the persistence of populations in addition to habitat loss due to mining, tourism or skiing and habitat degradation due to overgrazing in the shrinking mountaintop habitats of V. graeca.

Introduction

Reptiles are exposed to increasing threats worldwide, with about 20% of species facing the risk of extinction due to the interaction or synergy of factors including habitat loss, fragmentation and degradation, environmental pollution and spread of diseases, among others^1,2. These impacts are further intensified by the rapidly changing thermal environment, which present significant ecological and physiological challenges to all ectotherms³. Despite their widely acknowledged thermal constraints, increasing evidence underscores the diverse coping strategies that reptiles employ to adapt to variation in climate. These strategies may include physiological and behavioral plasticity, shifting their geographic ranges, or ultimately undergoing evolutionary adaptation^4,5,6. Central to these responses is the essential role of thermal trait variation. The evolution of all physiological and behavioral traits hinges on the existence of inherent phenotypic and genetic variation, which is presumed to be linked to fitness⁷. Therefore, comprehending the variability in thermal traits at the individual and population level is essential. This understanding can serve as a predictive tool for evaluating reptiles' vulnerability to climate change^8,9.

Temperature is a critical environmental factor which determines the geographical distribution and abundance of populations of any species, particularly reptiles^10,11. Ectotherms depend on effective thermoregulation, balancing physiological and behavioural mechanisms to maintain their body temperature (T_b) within an optimal range¹². This temperature regulation is crucial for essential life processes, including foraging, reproduction, and predator avoidance¹³. Reptiles exhibit a range of thermoregulatory strategies, influenced by the specific environmental conditions and the relative costs and benefits of various thermoregulatory behaviours^14,15. These strategies span from thermoconformity, where there is minimal active temperature regulation, to active thermoregulation, where reptiles actively seek or avoid heat sources to control their body temperature¹⁶. The variability in thermoregulatory behaviours among reptile species reflects their adaptation to local environmental conditions^17,18. For instance, the timing of daily and seasonal activities plays a significant role in determining body temperature, furthermore reptiles often select specific thermal microhabitats that facilitate optimal thermoregulation during active periods⁵. The rapid changes in temperature, particularly warming, are pushing many species towards the limits of their thermal tolerance¹⁹. This shift leads to altered activity patterns, which can result in reduced fitness, and ultimately in changes in population size and distribution³. Cold-adapted species and those in higher latitudes or altitudes are, particularly, exposed to acute challenges due to significant temperature increases, as their narrow thermal tolerances could amplify their vulnerability²⁰.

Forecasting the impacts of environmental change on reptiles depends on evaluating reptiles' ability to regulate their T_b ²¹. However, even differentiating whether variations in T_b arise from environmental fluctuations (i.e., thermal constraints) or from changes in thermal behaviour or physiology, such as plasticity in preferred or optimal T_b, represents a challenge. Fortunately, a detailed protocol for defining the accuracy and effectiveness of thermoregulation has been developed and refined over time^16,22,23. Furthermore, the suite of thermal ecology traits defined as measurable proxies of thermoregulatory behaviors (Table 1)^21,24, directly or indirectly influence an individual's habitat choice, temporal and spatial activity patterns, and behavior within the thermal environment^25,26,27,28. Consequently, these traits also shape and can serve as indicators of the distribution and abundance patterns of populations or species within ecosystems^29,30.

Table 1 Variables used in this study to assess the thermal ecology of Vipera graeca.

In a rapidly changing thermal environment, local adaptation could be essential for the survival of reptile populations, given their limited ability to disperse and the narrow window of time available for shifting their distribution range or undergoing adaptive evolutionary changes^5,6,31. Physiological and behavioral plasticity can help ectotherms cope with warming by reducing the thermal sensitivity of vital processes and enhancing physiological tolerances^32,33. Additionally, these plastic responses can minimize exposure to harmful or lethal temperatures^11,34. However, the pace of climate warming might outstrip the capacity for plastic responses or increase the costs associated with behavioral strategies^26,35,36.

This study aims to comprehensively assess the thermoregulatory behaviour of Vipera graeca, a cold adapted endangered alpine snake particularly threatened by climate change and to evaluate its degree of vulnerability to its effects. Perviously, the thermoregulation of the species was almost unknown. In a recent paper, we suggested that vipers do not fully exploit the thermally optimal time window available to them, likely because they shift their activity to periods with fewer avian predators¹³, and we provided the measurement of voluntary thermal maxima (VTM), which revealed high upper thermal tolerance of grassland vipers in alpine environments, which was 36.6 °C in the case of V. graeca³⁷.

We use data from on-site measurements, field tests and the thermal landscape to address five pivotal questions. (Q1) Does their preferred body temperature depend on the sex and body size of individuals? (Q2) Do the populations of V. graeca differ in preferred body temperature? (Q3) Are the differences in thermal characteristics among populations, if any, influenced by latitudinal or phylogenetic position? (Q4) To what extent can individuals accurately and effectively thermoregulate to achieve the preferred range of body temperatures and exploit the thermal landscape or quality of their habitats? (Q5) Do current climatic conditions represent thermal limitations for the activity of V. graeca and, if so, how will these limitations change concerning future climate change? By addressing these questions, the study aims to provide a detailed understanding of the thermoregulatory strategies and adaptations in V. graeca, considering a range of biological and environmental influences at both the individual and the population levels. Additionally, mitigating the potential impacts of climate change on these vipers is crucial for effective conservation efforts^38,39.

Results

In total, we collected temperature data on n = 74 Vipera graeca individuals from five populations (Fig. 1; Table 2). The overall preferred body temperature (T_p) was 28.77 ± 0.43 °C (mean ± SE, Table 2). The sex of individuals did not have an effect on T_p (Wilcoxon W = 425, P = 0.824; Linear model coefficient slope b = 0.146, P = 0.873). We found that in males, snout-vent length (SVL) did not influence T_p (b = − 0.009, P = 0.717), however, it had a significant positive effect in females (b = 0.02, P = 0.003; Fig. 2A). Gravid reproductive state in adult females did not influence T_p (W = 46, P = 0.23; b = 1.63, P = 0.22).

Figure 1

Panel (A) The distribution of Vipera graeca at 50 × 50 km spatial resolution Mizsei et al.⁵², and the area of panel (B) (black rectangle). Panel (B) Buffered polygons of suitable habitats of V. graeca based on a distribution model Mizsei et al.⁴⁶ and sampling sites of this study (red dots), and a picture of a representative female V. graeca individual sampled during the data collection. Maps were generated in QGIS 3.14 using the QuickMapServices plugin (QGIS.org (2024). QGIS Geographic Information System. Open Source Geospatial Foundation Project. http://qgis.org).

Table 2 Sample size of temperature measurements, number of individuals and snout-vent length (SVL) means recorded in this study.

Figure 2

Influence of snout-vent length (SVL) on T_p in Vipera graeca females. Solid lines show the prediction of LM, and dashed lines delimit 95% CI (panel A). T_set (solid interval lines) and T_p (dots) of measured populations of Vipera graeca (panel B). Grey lines show individual T_set values. For populations acronyms see Fig. 1.

We found that most populations did not differ from each other in T_p, except for a significant difference in T_p among the southernmost and northernmost localities, the Vardoussia and the Tomorr populations (b = 3.3, P = 0.004, Fig. 2B). Latitude positively influenced T_p (b = 1.557, P = 0.006), lower set-point range of body temperature (T_set) (b = 2.375, P = 0.002) and upper T_set (b = 0.872, P = 0.059; Fig. 3). The phylogenetic position in the sampled V. graeca populations did not constrain T_p as we found no evidence of a significant phylogenetic signal (κ = 1.13, P = 0.185; λ = 1.276, P = 0.601).

Figure 3

Influence of latitude on T_set and T_p of Vipera graeca. Dots show individual values, solid lines show prediction of LM, and dashed lines delimit 95% CI.

We found that the overall mean accuracy of body temperature (d_b) was 0.45 ± 0.26 °C (mean ± SE), which indicated high accuracy of thermoregulation (Table 3). We found no difference among the populations regarding d_b based on the overlapping 95% CI of estimates. The mean thermal quality of the habitat (d_e) was 4.97 ± 0.014 (mean ± SE, Table 3). At the habitats of Vardoussia and Tymfi, we found significantly lower d_e compared to the other habitats, indicating higher thermal quality at these sites. The cross-population mean of the effectiveness of thermoregulation sensu Hertz (E) was 0.92 ± 0.055 (mean ± SE, Table 3) and the effectiveness of thermoregulation sensu Boulin-Demers and Weatherhead (I) was 4.9 ± 0.3 (mean ± SE, Table 3). Spot-on T_b measurements never exceeded the upper limit of T_set, and overall, they overlapped with the range of T_set, except for some early morning T_b observations (Fig. 4). The overall exploitation of thermal landscape (E_x) was 2.57.

Table 3 Thermal ecology indices of Vipera graeca.

Figure 4

Spot-on T_b measurements of Vipera graeca individuals (dots) observed basking (red dots) or in shade of vegetation (green dots), and the limits of T_set (target body temperature, dashed lines) of each population. Solid lines indicate 95% CI range of T_e at full sun (red lines) and in the shade of vegetation (green lines) measured at four populations of Vipera graeca in July and August.

At current climatic conditions, the mean annual restriction time (h_r) of V. graeca was h_r = 1278.4 ± 9.7 h (mean ± SE) across all localities (Fig. 5). Latitude significantly decreased h_r (b = − 210.9, P < 0.0001). For future conditions, regarding the SSP1-2.6 scenario, it is expected to increase to h_r = 1606 ± 9.0 h (mean ± SE), or based on the SSP5-8.5 scenario, to h_r = 2086 ± 8.4 h (mean ± SE). Thus, the time when T_e is above the upper limit of T_set is expected to increase by 328 or 806 h for the SSP1-2.6 and SSP5-8.5 scenarios, respectively. h_r was significantly higher in SSP5-8.5 compared to SSP1-2.6 (W = 2519, P < 0.0001).

Figure 5

Expected change in h_r at Vipera graeca habitats from current climate to two future scenarios for all known localities. Inset plots show the change in mean h_r for each population. Colour of the arrows in the inset plot correspond to the expected change in h_r depicted in the maps. For populations acronyms see Fig. 1.

Discussion

Our results indicate that T_p of Vipera graeca is unaffected by sex or gravid state in females, although female size showed a positive influence on T_p. We observed a significant difference in T_p only between the northernmost and the southernmost studied populations of V. graeca, however, latitude positively influenced T_p, while we found no evidence that the phylogenetic position has a constraint. The study revealed a high accuracy of thermoregulation in V. graeca populations and some variation in the thermal quality of habitats among the studied locations. The overall effectiveness of thermoregulation was high, indicating that V. graeca is successful in achieving its target temperatures and we also found that this species effectively exploits the thermal landscape of its habitats. Under current climatic conditions, V. graeca faces a mean annual activity restriction of approximately 1278 h across all localities, and this restriction time expected to increase significantly under future climate change scenarios.

The thermal ecology of grassland vipers, including species like Vipera ursinii, Vipera renardi, and their relatives, remains largely unexplored due to the scarcity of published studies. This gap in research limits our ability to interpret our findings on V. graeca with existing literature. However, available data on V. ursinii rakosiensis suggest an optimal body temperature "somewhat below 35 °C", as observed in radiotagged individuals⁴⁰. This optimal temperature is notably higher than the T_p and T_set we recorded for V. graeca. For V. u. rakosiensis, this temperature might align more closely with their voluntary thermal maxima (VTM) rather than their upper T_set³⁷. Further comparison with V. ursinii moldavica, based on spot-on T_b recordings, shows some overlap with our observations for V. graeca⁴¹. Beyond grassland vipers, additional insights can be drawn from studies on Vipera berus, which can also be categorized as a cold-adapted species. Our findings indicate that V. graeca exhibits a lower T_p and a much wider range of T_set than V. berus¹⁴, suggesting that V. graeca may possess a broader thermal niche.

Consistent with our findings, similar patterns have been observed in other viper species regarding the lack of significant difference in T_b between reproductive and non-reproductive females. Specifically, in V. berus, a study found no significant difference in T_p among reproductive states⁴². A parallel observation was made in V. ammodytes, where no distinct differences in T_b were noted concerning the reproductive state; however, gravid females tended to select warmer microhabitats, possibly to achieve more precise thermoregulation⁴³. In contrast, V. aspis presents a differing scenario, where gravid females, particularly in the latter half of gestation, tend to maintain a higher T_b⁴⁴. This variability across species underscores the complexity of thermoregulatory behaviors in vipers and suggests potential species-specific strategies related to reproductive state and thermoregulation.

Our study revealed that V. graeca in higher latitude habitats, which are characteristically colder, exhibit a preference for higher temperatures, both in terms of T_p and T_set. This finding aligns with studies indicating that reptiles in higher elevations and latitudes tend to have increased upper thermal tolerance limits, such as thermal-safety margin, critical thermal maximum (CT_max), and VTM^11,37. This pattern may be attributed to reptiles becoming more effective thermoregulators in challenging thermal environments with limited heat resources⁴⁵. The absence of a significant phylogenetic influence on T_p in our study suggests that the observed variations in T_p across different populations of V. graeca could be the result of local adaptation and/or phenotypic plasticity rather than of evolutionary descent. Future research in this field should, therefore, include phenotypic plasticity and the capacity for reversible phenotypic plasticity as a fundamental component of experimental designs, providing a comprehensive understanding of how reptiles adapt to their thermal environments⁴⁶.

Our data shows that V. graeca generally maintains its T_b within the preferred T_set, except for three “cold” individuals, which were likely in the initial stages of morning basking. This accuracy of thermoregulation potentially leads individuals to maintain high physiological performance by active thermoregulation in a challenging environment. Interestingly, these snakes do not utilize the entire available thermoregulatory time window despite suitable environmental temperatures to keep T_b in T_set. This behaviour might be influenced by factors beyond thermal conditions alone. For instance, a study suggest a behavioral adaptation in V. graeca, where shifts in their daily bimodal activity patterns—earlier mornings and later afternoons—could be attributed to predator avoidance strategies, particularly evasion of visually-searching raptor birds such as the short-toed snake eagle (Circaetus gallicus)¹³. Additionally, variations in the exploitation of thermally suitable times are likely to be observed across seasons, which may reflect changes in environmental conditions and the physiological needs of the individuals. Factors such as size (or age), foraging requirements, and reproductive objectives could also play significant roles in shaping these thermoregulatory activities. While most existing studies in thermal ecology, including this study, focus on short-term, intensive data collection periods, future research should aim to extend these observational windows. Covering longer periods would provide a more comprehensive understanding of how seasonal variations, life-history traits, and ecological interactions, such as predator–prey dynamics, influence the thermal ecology of V. graeca and other reptilian species.

In a previous study employing correlative species distribution modelling, we projected that around 90% of V. graeca’s current habitats might be lost by the end of the twenty-first century due to climate change⁴⁷. These correlative models, praised for their simplicity and applicability across various species with available distribution and environmental data⁴⁸, contrast with more accurate process-based mechanistic models⁴⁹. The latter includes detailed physiological mechanistic models that solve coupled energy and mass balance equations to establish an explicit link between the organism’s requirements and the environmental availability of resources⁵⁰. Our current analysis, while not forecasting changes in distribution, focuses on how thermoregulatory time budgets may shift. We observed a significant increase in h_r, the duration when environmental temperatures surpass the upper limit of the species’ temperature preference. This rise in h_r implies reduced opportunities for thermoregulation, foraging, and reproduction, potentially leading to direct impacts on individual fitness and population growth rates. For a species such as V. graeca, which is already endangered and confined to shrinking and deteriorating mountaintop habitats, this reduction in key demographic factors could exacerbate the risk of local extinctions.

In conclusion, our study on V. graeca revealed critical insights into the species’ thermal ecology. We found that preferred body temperature varies significantly with latitude and female body size, but is not influenced by sex or reproductive state. The species shows high accuracy in thermoregulation and level of exploitation of the thermal landscape. However, our findings also indicate an impending challenge: as climate change progresses, V. graeca faces increasing restrictions in thermoregulatory opportunities, which could significantly impact its foraging and reproductive behaviours, thereby affecting overall fitness and population growth. There is a crucial need for considering both local environmental conditions and potential impacts of climate change in the conservation strategies for V. graeca, to mitigate the risks of habitat loss and population decline in this endangered species.

Methods

Study species

The Greek meadow viper (Vipera graeca Nilson & Andrén, 1988) is a small-sized grassland viper living exclusively in the subalpine and alpine grasslands of the Pindos mountain range in Southern Albania and Central Greece, usually between 1600 and 2100 m above sea level^47,51,52. It has a severely fragmented distribution forming 17 known isolated populations on mountaintops above the tree line (Fig. 1)^47,53. The habitats of the species are threatened by climate change and anthropogenic degradation such as overgrazing, and most likely 90% of the habitats will probably disappear by the end of this century⁴⁷. The species is currently listed as endangered by the IUCN Red List⁵³. The populations genetically form two major lineages, a southern (Vardoussia among the sampled populations in this study) and a northern one (all other populations sampled in this study, Fig. 1), without signs of inbreeding or genetic drift. This species is a dietary specialist on Orthopterans⁵⁴.

Data collection

We collected data on five populations of Vipera graeca along a latitudinal gradient from north to south: Tomorr, Kulmak (Albania), Tymfi, Lakmos and Vardoussia mountains (Greece), including the northernmost and the southernmost populations (Fig. 1). All populations were sampled in the activity peak in summer: 29 July–14 August 2017 (Tymfi), 14–22 August 2017 (Lakmos), 23 July–12 August 2018 (Vardoussia), July 25–August 15 2019 (Tomorr and Kulmak).

The evaluation of the thermoregulatory characteristics of ectotherm animals such as V. graeca, requires information on the availability and distribution of body temperatures a non-thermoregulating individual could or would achieve based purely on energy flux from radiation, conduction, and convection¹⁶. Accordingly, we collected three types of temperature data for this study: operative temperature (T_e), field body temperature (T_b) and target body temperature in a thermal gradient (i.e. set-point temperature range T_set). Operative or environmental temperature (T_e) is the equilibrium of T_b, as it gives the null distribution of T_b expected from non-thermoregulating individuals⁵⁵. We measured T_e by using the physical model (operative temperature model, OTM) made of a copper tube that mimicked the size (18 × 350 mm), shape, and heat absorption of the study species. To measure T_e, we equipped the models with temperature data loggers (iButton DS1921G-F5#, Thermochron Ltd., Castle Hill, NSW, Australia; 0.5 °C resolution, ± 1.0 °C accuracy) pre-set to record data in 5-min intervals. We placed the models in two different micro-environments, with one exposed to the sun and the other under the shade of vegetation cover. The models were placed in the micro-environments closest to the exact location of capture of V. graeca individuals.

The distribution of body temperature of an actively thermoregulating animal is expected to differ from T_e. To measure field T_b we sampled V. graeca by spot sampling often called “grab and jab” sampling¹⁶. We intensively searched for vipers in their habitats, and when a viper was spotted, we carefully captured it using protective gloves and immediately measured its cloaca temperature with a thermometer (Testo 103, Testo SE & Co. KGaA, Baden-Württemberg, Germany; 0.1 °C resolution, ± 0.5 °C accuracy). We recorded the GPS coordinates and the snake was subsequently transported to our field camp, where other measurements were conducted. Each individual was housed in a linen bag in a shaded area to ensure their well-being during the temporary keeping.

To better understand the thermoregulation of ectotherms it is central to identify the target body temperature that an individual would achieve⁵⁶. Target (or preferred, or selected) body temperature can be measured in artificial thermal gradients, in an environment that is independent from the ecological costs and constraints that can influence thermoregulation in the field¹⁶. To measure T_set we kept the individuals in a thermal gradient set up in the field, from 6:00 AM to 6:00 PM and measured cloaca temperature hourly with a thermometer (Testo 826-T4, Testo SE & Co. KGaA, Baden-Württemberg, Germany; 0.1 °C resolution, ± 0.5 °C accuracy). The thermal gradient was set up in 100 × 30 × 30 cm polycarbonate terrariums (dimensions: 100 × 30 × 30 cm; n = 4) positioned under permanent shade. To establish the desired thermal gradient (max. 20 °C on the cool and min. 40 °C on the hot end), the hot side of the terrarium was heated using a 200 W ceramic heat wave source (Exo Terra PT2046; Rolf C. Hagen, Inc. Montreal, QC H9X 0A2, Canada). Additionally, to maintain an optimal temperature range, evaporative cooling was utilised whenever the air temperature at the cooler end of the gradient reached 20 °C. This was achieved by covering the initial 30 cm of the terrarium's top surface with a water-soaked textile sheet, placed at the cooler end. Thermal measurements were followed by sexing, determination of the gravidity status of females and measuring the snout-vent length (SVL) of the individuals. After measurements, snakes were released at their exact capture location.

To assess the influence of change in the thermal landscape on the activity time of vipers we compiled a database of all V. graeca occurrence records available from our previous studies and the observations recorded in this study (n = 378 records from all known populations)^{37,47,51,52,53}. We downloaded data on monthly air temperatures at 30″ spatial resolution (~ 720 × 930 m grid) from the Worldclim 2.1. database⁵⁷ for the entire geographic range of V. graeca. To characterize future climatic conditions (years 2081–2100), we selected three Global Climate Models (GCMs): HadGEM3-GC31-LL, IPSL-CM6A-LR, MIROC6⁵⁸ and two Shared Socioeconomic Pathway (SSP) scenarios, an optimistic one (SSP1-2.6) and a pessimistic one (SSP5-8.5). The optimistic SSP1-2.6 scenario represents a future world characterized by low population growth, strong sustainability efforts, and significant reductions in greenhouse gas emissions, i.e., representing a future in which ambitious actions are taken to mitigate the worst effects of climate change⁵⁹. The pessimistic SSP5-8.5 represents a high-emission scenario in which greenhouse gas concentrations continue to rise throughout the twenty-first century, leading to a world with high levels of carbon dioxide and other greenhouse gases⁵⁹.

Data analysis

To describe the thermoregulation of Vipera graeca we calculated the indices and variables commonly used to describe the thermoregulation of free-living animals (Table 1)²⁴. All data operations and statistical analyses were conducted in the R 4.1.3 statistical environment in a fully reproducible way⁶⁰.

Because most of the ectotherms appear to thermoregulate between an upper and lower temperature limit rather than around a single temperature value (T_p), it is straightforward to determine T_set. The upper and lower limits of preferred T_b could be estimated by calculating the central 50% of the distribution of all T_b values selected in a thermal gradient¹⁶. To calculate T_p and T_set, we used an approximation function on the density distribution of body temperatures measured in thermal gradient for each individual. First, we fitted the density function to the data, second, we transformed the probability values (y-axis) to range between 0 and 1. Third, using the approxfun function we estimated the probability values for a sequence of temperature values from 15 to 45 °C by 0.05 °C steps. Finally, we read the T_p, and T_set values: at probability = 1 as T_p, min. value at probability = 0.5 as T_set lower value and max. value at probability = 0.5 as T_set upper value.

To assess the influence of SVL on T_p (Q1) we fitted linear models (LM) using the lm function for female and male individuals separately. To assess the effect of gravidity on females’ T_p we used Wilcoxon rank test using the wilcox.test function and fitted an LM. We investigated the differences in T_p across populations (Q2) by fitting a LM, and then we applied multiple comparisons on the LM using the glht function of the multcomp package to identity significantly (P < 0.05) different groups⁶¹. To assess the influence of latitude on T_p and T_set, we also applied LMs at the level of individuals (Q3). We fitted three LMs separately for each dependent variable (T_p, T_set lower, T_set upper) using the same independent variable (latitude).

To test for the presence of a phylogenetic signal in T_p (Q3), we used the phylogenetic tree reconstruction based on RADseq data (see supplementary material). To import the tree file of that phylogenetic reconstruction and for the subsequent analysis, we applied the ape and phytools packages^62,63. To calculate branch lengths we used the compute.brlen function of ape. We tested for the presence of phylogenetic signal by computing Pagel’s lambda⁶⁴ and Bloomberg’s K⁶⁵ on T_p using the phylosig function of phytools package.

We applied the index of accuracy of body temperature (d_b) by calculating the mean of the deviations of T_b from T_set (Table 1)¹⁶. To assess the thermal quality of the habitats (d_e), we used the mean of the deviations of T_e from T_set (Table 1)¹⁶. For this calculation, we used a subset of T_e measurements to encompass the diurnal activity time (6:00 AM to 8:00 PM) of the study species using the lubridate package⁶⁶. d_b and d_e values with 95% confidence intervals were calculated for each measured population (n = 4 and n = 5 respectively). We assessed the effectiveness of thermoregulation (Q4) by calculating the E index (E = 1 − (d_b/d_e))¹⁶ and the I index (I = d_e − d_b)²². To calculate the thermal exploitation of the habitats E_x²³, we counted cases when spot-on T_b measurements were inside the range of T_set and divided by the number of all T_b measurements when the corresponding T_e was overlapping with T_set.

To estimate the expected change in activity time (Q5), we calculated the restriction time of activity as the number of hours during a full year (24 × 365) when T_e exceeded the upper limit of T_set²⁴. To get T_e for an entire year and all occurrence areas (n = 102 at the spatial resolution of 30″) of V. graeca, we used a microclimate model to establish a simulated environment and an ectotherm model to predict T_e using the NicheMapR package⁶⁷. As the default temperature data in the micro_global modelling function has a low spatial resolution (10 × 10 km grid), we used the current (1981–2010) monthly mean air temperature data extracted for each viper occurrence grid cell from the Worldclim 2.1. database. All raster operations were done using the terra, raster and dismo packages^68,69,70. We set the warm parameter to the difference of means at the fine resolution of the default monthly temperature values to establish a microclimate model for current climatic conditions. In this simulated environment, we fitted an ectotherm model by default parameters, except for the weight of the animal model (33.5 g, mean of all adult V. graeca weighted, n = 192), the shape of the animal (cylinder), and diurnality. We fitted the ectotherm model as “dead” (live parameter set to 0), meaning non-thermoregulating, as we were interested in getting T_e measurements. We extracted the T_e estimates from the ectotherm model output, and counted the number of hours when T_e exceeded the T_set (= restriction time, h_r). To get h_r for future climate, we used the warm parameter similarly as above using the extracted monthly air temperature of the n = 6 future data sets (3 GCM × 2 SSP). We repeated this procedure for each location to map h_r for the whole distribution of V. graeca. Finally, we calculated the difference between future and climate restriction times to predict the change in h_r of activity, and we mapped the h_r results for presentation purposes to the same raster resolution as the Worldclim data we used for microclimate modelling.