Forecasting the response to global warming in a heat-sensitive species

Avoiding hyperthermia entails considerable metabolic costs for endotherms. Such costs increase in warm conditions, when endotherms may trade food intake for cooler areas to avoid heat stress and maximize their energy balance. The need to reduce heat stress may involve the adoption of tactics affecting space use and foraging behaviour, which are important to understand and predict the effects of climate change and inform conservation. We used resource selection models to examine the behavioural response to heat stress in the Alpine ibex (Capra ibex), a cold-adapted endotherm particularly prone to overheating. Ibex avoided heat stress by selecting the space based on the maximum daily temperature rather than moving hourly to ‘surf the heat wave’, which minimised movement costs but prevented optimal foraging. By integrating these findings with new climate forecasts, we predict that rising temperatures will force mountain ungulates to move upward and overcrowd thermal refugia with reduced carrying capacity. Our approach helps in identifying priority areas for the conservation of mountain species.


-Small scale resource selection analysis
.1 -GLMM beta estimates for small-scale resource selection by male ibex observed from 2010 to 2011 in the Gran Paradiso National Park, Italy. Beta coefficients were plugged in the exponential resource selection function RSF after dropping the intercept, resulting in the resource selection patterns depicted in Figs. S1. 1    Figure S1.4 -Relative probability of selection for aspect (cos-transformed, x-axes) interacted with wind direction (different colours represent different scenarios for wind direction) and wind speed (the two plots represent two different scenarios) as predicted by the smallscale resource selection function, which was built using male ibex observations collected from May to October (2010-2011) in the Gran Paradiso National Park, Italy. Figure S2.1 -Large-scale resource selection function (RSF) evaluation: area-adjusted frequency of categories (bins) of RSF scores. The evaluation implied calculating the correlation between RSF ranks and area-adjusted frequencies for a withheld sub-sample of data, e.g. 1/5 of the data in a 5-fold crossvalidation scheme. We investigated the pattern of predicted RSF scores for partitioned testing data (presence-only) against categories of RSF scores (10 bins). A Spearman rank correlation between area-adjusted frequency of cross-validation points within individual bins and the bin rank was calculated for each cross-validated model. A model with good predictive performance would be expected to be one with a strong positive correlation, as more use locations (area-adjusted) would progressively fall into higher RSF bins. In this case, the 5-fold cross-validation showed that the largescale resource selection model (Table 1 of the main manuscript) with daily maximum temperature as predictor had outstanding predictive ability on withheld data (Spearman correlation coefficients: ρ fold1 = 0.988, ρ fold2 = 0.988, ρ fold3 = 0.976, ρ fold4 =0.976, ρ fold5 = 0.988).

Figure S2.2 -
Small-scale resource selection function (RSF) evaluation: area-adjusted frequency of categories (bins) of RSF scores. The evaluation implied calculating the correlation between RSF ranks and area-adjusted frequencies for a withheld sub-sample of data, e.g. 1/5 of the data in a 5-fold crossvalidation scheme. We investigated the pattern of predicted RSF scores for partitioned testing data (presence-only) against categories of RSF scores (10 bins). A Spearman rank correlation between area-adjusted frequency of cross-validation points within individual bins and the bin rank was calculated for each cross-validated model. A model with good predictive performance would be expected to be one with a strong positive correlation, as more use locations (area-adjusted) would progressively fall into higher RSF bins. Compared to the large-scale resource selection model ( Fig.  S2.1), the small scale RSF had a weaker -but still very good -predictive ability on withheld data (ρ fold1 =0.881, ρ fold2 =0.912, ρ fold3 = 0.952, ρ fold4 =0.967, ρ fold5 =0.939). c)

-Temperature interpolation models Background
We collected high-resolution temperature data in our study site and built interpolation models to predict fine-scale temperature variation over space and time during the monitoring period. We combined weather data collected at the nearest local weather station with temperature logger data (iButton DS1922L, Maxim Integrated) that we located throughout the study site over the study period.
We had two different expectations about the effect of temperatures on ibex habitat selection: we hypothesised that ibex either select habitat based on the actual temperature (i.e., hourly temperature)meaning that they are more likely to be located where the conditions are optimal in that moment -or based on the overall daily temperature (i.e., daily maximum temperature) -meaning that they are located at the elevation where the conditions will be optimal when the maximum temperature will be recorded. As a consequence, we built two temperature interpolation models predicting hourly temperature and daily maximum temperature, respectively.

Temperature loggers
We air temperature every hour. We described below the two steps that were necessary to build the temperature interpolation models.  Step 1: filling in missing data in the dataset of the temperature loggers Some of the loggers stopped recording temperature data, providing measurements for about 70-80% of the study period. For each data logger with missing data, we fitted linear models with air temperature as a response variable and the following covariates as predictors: time of the day (including quadratic and cubic term to account for nonlinearity), Julian day (including a quadratic term), year of study, their interactions (Julian date × time of the day, Julian date × year of study), and the air temperature recorded by the Pont weather station or by any other temperature logger.
Regarding the latter predictor, we selected the temperature series that had the strongest correlation with the air temperature data collected by the target logger for which we wanted to replace missing data. Because of difference in altitude between the weather station or temperature loggers and the target logger, in some linear models we temporally shifted (1-3 hours) the response variable (temperature of the target logger) and the predictor (temperature recorded by the weather station or another logger selected based on the strength of the correlation) in order to achieve the best predictive model.
Once we built the starting linear model structure for each target logger, we then ran a forward and backward stepwise algorithm (step function of the stats package in R) and selected the best model structure with the lowest AIC values. We reported the adjusted R 2 as a measure of the predictive ability of each model, which we used to predict missing data and complete hourly temperature series of iButton loggers. After filling in all gaps, we calculated the daily maximum air temperatures.
Step 2: predicting hourly temperature and daily maximum temperature over space and time

(interpolation models)
We used the logger temperature datasets to build spatial and temporal interpolation models We repeated the same procedure and built two alternative GAMs to model daily maximum temperatures. We used the same set of predictors as for the hourly temperature models, this time removing the time of the day and its interactions because meaningless when modelling the daily maximum temperature.
We random-cross-validated our alternative models to verify their ability to predict on new data. We trained the best model on 80% of the data, predicted on the remaining 20%, and computed the R 2 of the relationship between predicted and observed temperature data. We extracted the average R 2 after repeating the procedure 50 times.

Results
Step 1: filling in missing data in the dataset of the temperature loggers Best linear models used to fill in missing data of temperature loggers were reported in  Step 2: predicting hourly temperature and daily maximum temperature over space and time

(interpolation models)
We reported the two alternative models explaining the variability of hourly temperatures in Table S6.2. The model with elevation fitted as smooth term in the GAM outperformed the alternative model (with elevation fitted as linear and quadratic term) and was used to predict hourly temperature over space and time during the study period. Likewise, we reported the two alternative models explaining the variability in daily maximum temperatures in Table S6. 3. Also in this case, the model with elevation fitted as smooth term in the GAM outperformed the alternative model (with elevation fitted as linear and quadratic term) and was used to predict daily maximum temperature over space and time during the study period.

-Calculation of buffer sizes needed to depict random availability in small-scale resource selection analysis.
To define the buffer size for small-scale resource selection analysis (i.e., the area around each ibex presence location where to depict random availability), we made use of satellite telemetry data.
We screened telemetry data available for male ibex in the Levionaz valley (year 2013), and we estimated ibex monthly mobility. Such movement behaviour can give us a clue on how far an ibex can move monthly, and thus define where to set the limit for sampling available resources in small-scale resource selection analysis. Because we used telemetry data collected in 2013 to define small-scale availability for ibex observed in 2010 and 2011, we decided to keep a more conservative approach and use a monthly temporal scale rather than depicting buffer sizes daily. The latter approach might have been biased by different environmental conditions occurring in the two different study periods.

-Sensitivity analysis aimed at defining the minimum number of random available points to be associated with ibex presence data in resource selection analyses.
Aim of this supplementary information is to run a sensitivity analysis and define the minimum number of available points that need to be associated to ibex presence data in resource selection functions.

Methods
Following recommendations by Ciuti et al. 3 and Roberts et al. 4 (see box 2 therein), we fit a generalized linear mixed-effect model with Bernoulli distribution of errors, presence (1) and availability (0) as response variable, and air temperature predicted by interpolation models as predictor. We used air temperature because it was the environmental covariate collected at the finest spatial resolution (10 x 10 m). We fitted the stratum-ID (identifying each pair of used location with its associated random available locations) nested within the individual-ID (identifying the individually recognizable ibex) as random intercept in the model; we also fitted the group-ID (identifying the identity of the group where the marked ibex was observed) as (crossed) random intercept.
We started the sensitivity analysis with an available:used location ratio 1:1 by running the GLMM 10 times, each time drawing a new random spatial sample of available locations. We repeated this procedure several times, stepwise increasing the number of random available points associated with each used point (i.e., available:used 2:1, 3:1, …, 50:1), and extracted the estimated model parameters (beta estimates). We thus screened the variation of beta estimates as a function of the number of random available locations, and fit a generalized additive model (response: beta estimates; predictor: number of random available locations per used location) to describe the relationship and detect model parameter stabilisation. Based on the visual inspection of GAMs, we selected the thresholds for the minimum number of random points needed to get stable parameter estimates.      Data from temperature loggers were combined with those collected by the weather station and used to build interpolation models predicting hourly and maximum daily temperature (°C) for each 10 x 10 m pixel of the study area at any given day within the study period RADIATION Solar radiation (W/m 2 ) recorded at the weather station WIND_SPEED Speed of the wind (m/s) recorded at the weather station.

WIND_DIRECTION
Direction of the wind recorded at the weather station, cosine-transformed to range between -1 with wind blowing from the South and +1 with wind blowing from the North.