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-2100 meters above sea level (Mizsei et al. 2016, 2018a, 2021). It has a severely fragmented distribution forming 17 known isolated populations on mountaintops above the tree line (figure 1, Mizsei et al. 2018b, 2021). 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 (Mizsei et al. 2021). The species is currently listed as endangered by the IUCN Red List (Mizsei et al. 2018b). The populations genetically form two major lineages, a southern and a northern one, without signs of inbreeding or genetic drift (Mizsei et al. in prep.). This species proved to be a dietary specialist on arthropods, and preys almost exclusively on Orthopterans (Mizsei et al. 2019). The thermoregulation of the species was almost unknown. In a recent paper, Radovics et al. (2023a) 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 in Radovics et al. 2023b they provided the measurement of voluntary thermal maxima (VTM), which revealed high upper thermal tolerance of grassland vipers in alpine environments, which was 36.6±0.24 °C (mean±SE) in the case of V. graeca.
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 (figure 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 (Hertz et al. 1993). Accordingly, we collected three types of temperature data for this study: operative temperature (Te), field body temperature (Tb) and target body temperature in a thermal gradient (i.e. set-point temperature range Tset). Operative or environmental temperature (Te) is the equilibrium of Tb, as it gives the null distribution of Tb expected from non-thermoregulating individuals (Bakken and Gates, 1975). We measured Te by using the physical model (operative temperature model, OTM) made of a copper tube and a small thermometer that mimicked the size (18×350 mm), shape, and heat absorption of the study species. To measure Te, 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-minute 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 Te. To measure field Tb we sampled V. graeca by spot sampling often called “grab and jab” sampling (Hertz et al. 1993). 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 (Huey 1982). 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 (Hertz et al. 1993 and references therein). To measure Tset 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 previous studies (n=378 records from all known populations; Mizsei et al. 2016, 2018a, 2018b, 2021, Radovics et al. 2023b). We downloaded data on monthly air temperatures at 30” spatial resolution (~ 720×930 m grid) from the Worldclim 2.1. database (Fick & Hijmans 2017) 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 (Meehl et al. 2020) 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 (Tebaldi et al. 2020). The pessimistic SSP5-8.5 represents a high-emission scenario in which greenhouse gas concentrations continue to rise throughout the 21st century, leading to a world with high levels of carbon dioxide and other greenhouse gases (Tebaldi et al. 2020).
All the experimental protocols were performed according to relevant guidelines and regulations. Data collection protocols were approved and were then carried out according to the permits provided by the Ministry for Environment and Energy of Greece (number 158977/1757) and the Ministry of Environment of Albania (number 6584). All methods are reported in accordance with ARRIVE guidelines (https://arriveguidelines.org).
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, Taylor et al. 2020). All data operations and statistical analyses were conducted in the R 4.1.3 statistical environment in a fully reproducible way (R Core Team 2022).
Because most of the ectotherms appear to thermoregulate between an upper and lower temperature limit rather than around a single temperature value (Tp), it is straightforward to determine the range of target body temperatures called the set-point range of body temperature (Tset). The upper and lower limits of preferred Tb could be estimated by calculating the central 50% of the distribution of all Tb values selected in a thermal gradient (Hertz et al. 1993). To calculate Tp and Tset, 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-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 Tp, and Tset values: at probability = 1 as Tp, min. value at probability = 0.5 as Tset lower value and max. value at probability = 0.5 as Tset upper value.
To assess the influence of SVL on Tp (Q1) we fitted linear models (LM) using the lm function for female and male individuals separately. To assess the effect of gravidity on females’ Tp we used Wilcoxon rank test using the wilcox.test function and fitted an LM. We investigated the differences in Tp 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 (Hothorn et al. 2008). To assess the influence of latitude on Tp and Tset, we also applied LMs at the level of individuals (Q3). We fitted three LMs separately for each dependent variable (Tp, Tset lower, Tset upper) using the same independent variable (latitude).
To test for the presence of a phylogenetic signal in Tp (Q3), we used the phylogenetic tree reconstruction based on RADseq data (see supplementary material; Mizsei et al. in prep). To import the tree file of that phylogenetic reconstruction and for the subsequent analysis, we applied the ape and phytools packages (Paradis & Schliep 2019, Revell 2012). 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 (Pagel 1999, Freckleton et al. 2002) on Tp using the phylosig function of phytools package.
We applied the index of accuracy of body temperature (db; Hertz et al. 1993) by calculating the mean of the deviations of Tb from Tset (table 1). To assess the thermal quality of the habitats (de; Hertz et al. 1993), we used the mean of the deviations of Te from Tset (table 1). For this calculation, we used a subset of Te measurements to encompass the diurnal activity time (6:00 AM to 8:00 PM) of the study species using the lubridate package (Grolemund & Wickham 2011). db and de 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 (Hertz et al. 1993, E = 1 – (db/de)) and the I index (Boulin-Demers & Weatherhead 2001, I = de – db). To calculate the thermal exploitation of the habitats Ex (Christian & Weavers 1996), we counted cases when spot-on Tb measurements were inside the range of Tset and divided by the number of all Tb measurements when the corresponding Te was overlapping with Tset.
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 Te exceeded the upper limit of Tset (Taylor et al. 2020). To get Te 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 Te using the NicheMapR package (Kearney 2017). 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 (Hijmans 2023a), raster (Hijmans 2023b) and dismo (Hijmans et al. 2022) packages. We set the warm parameter to the difference of means of the fine resolution and 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 Te measurements. We extracted the Te estimates from the ectotherm model output, and counted the number of hours when Te exceeded the Tset (= restriction time, hr). To get hr 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 hr for the whole distribution of V. graeca. Finally, we calculated the difference between future and climate restriction times to predict the change in hr of activity. Finally, we mapped the hr results for presentation purposes to the same raster resolution as the Worldclim data we used for microclimate modelling.