Study area and species
The study was conducted in Sierra Nevada, mostly within the National Park (37°05′N, 3°28′W; SE Spain; Additional file 1: Figure S1). This park extends over 85,883 ha and is composed of mountains that rise over 3000 m a.s.l., ranging from 1700 to 3500 m. It is dominated by a continental Mediterranean climate with altitudinal gradients of temperature and rainfall. Rainfall is more frequent in spring and autumn, whereas summers are hot and dry and winters are cold with snowfall from November until April. The park has a highly diverse vegetation, with 2100 plant species (some endemic), structured in forests, shrubland and grassland along altitudinal gradients.
The Iberian ibex (Capra pyrenaica) is an endemic species of the Iberian Peninsula that inhabits mountainous systems [19]. This species live in social groups, but show spatial sexual segregation for most of the year, only coming together during the courtship (rutting) season, usually from October to December [25]. Kids are born in late spring, usually May. The Iberian ibex is a generalist herbivore, foraging as both a browser and grazer on a wide and varied diet that includes grass, shrubs and sometimes trees [19, 26]. Diet and foraging mode depend on resource availability [19]. Ibex can also perform altitudinal movements as to track seasonal resources that become available depending on climatic factors such as temperature and snow cover [19].
We conducted the study in two different geographical areas inhabited by different ibex population nuclei to control for possible effects of geographical idiosyncrasy, for example in the composition and density of resources. The areas differed mainly in altitude and vegetation cover, with the eastern nucleus being at lower altitudes and having a denser forest cover.
Movement data
We equipped 22 Iberian ibex with GPS collars (Microsensory, Córdoba, Spain, and Vectronic Aerospace, Berlin, Germany) after capturing them by darting using an anaesthetizing mixture of xylazine (3 mg/Kg) and ketamine (3 mg/Kg). Their movement was monitored over a maximum of two years during 2005–2007 by obtaining positions every one to every four hours depending on the animal. Because some ibex died and some stopped transmitting data before completing at least a full season, we had a final sample of 18 animals (10 males and 8 females) living in two separated geographical areas of the mountains (9 in each of the two nuclei; Additional file 1: Figure S1). All the included ibex were tracked for multiple seasons. After removing the first five fixes and obvious relocation errors, we had a total of 2085–4639 fixes per animal. However, to homogenize the time lags between successive relocations across tracked animals, we subsampled the movement data to obtain relocations every four hours, rendering a total of 700–3230 fixes per animal over a temporal range of 206–576 days.
Data analysis
We defined four different seasons according to the period of the year and the biology and life history of the Iberian ibex: spring (kidding season; April–June), summer (July–September), autumn (mating season; October–December) and winter (January–March). For each of the seasons, we first estimated for each of the 18 ibex habitat selection models that accounted for both movement and resource availability. These models allowed us to estimate selection and movement coefficients. Then, we explored to what extent selection-independent movement, selection strength, resource availability, and intrinsic factors (sex) determined seasonal home range size. All the specific analyses are described below.
Habitat selection models
Movement allows organisms to track the environment, and when estimating habitat selection, failure to account for the movement process may produce biased selection estimates [3]. A recent approach, termed integrated step selection analysis (iSSA) [4], builds on resource and step selection functions [27, 28] and can be used to model habitat selection while accounting for individual differences in movement behaviour. As such, this model can also be used to obtain estimates of selection-independent movement coefficients. We performed iSSA for each animal in each season to obtain individual, rather than population-level, estimates (as recommended in [4, 29]). iSSA simultaneously estimates movement and habitat selection parameters by comparing each used movement step with a set of conditioned available steps randomly sampled from an analytical distribution parameterised based on observed steps (N = 10 in this study). Movement steps were characterised by their length, i.e. the distance between the start-point and end-point of a given step, and direction, defined as the angular deviation (or turn angle) between successive steps. In particular, available step lengths were randomly sampled from a Gamma distribution fitted to observed step lengths of each animal in each season by maximum likelihood, and directions were randomly sampled from a uniform distribution of turn angles between successive steps.
Habitat covariate values were extracted for the end-point of each step and consisted of environmental variables related to foraging and roosting habitat: terrain slope (derived from a digital terrain model; www.juntadeandalucia.es/medioambiente/site/rediam), heat load (derived from the same digital terrain model) [30, 31], both with a spatial resolution of 10 m, and a primary productivity index (the Normalized Difference Vegetation Index, NDVI) obtained from satellite imagery (NASA product MOD13Q1; spatial resolution = 250 m, temporal resolution = 16 days). Because NDVI varies over time, each observed and control location was associated to the specific NDVI value corresponding to that location at the closest date. Heat load is an index of incident solar radiation that takes into account the orientation of the terrain slope, thus, depending on the latitude and season, it is associated to vegetation cover and snow depth. All habitat variables were centred and standardized and the respective quadratic effects included also as explanatory variables, as habitat selection might also show non-monotonic responses. Habitat variables were checked for collinearity by performing pair-wise correlation tests (correlation coefficients were all below 0.70; only in one case out of 76, we found a correlation lower than −0.70; Additional file 1: Table S1).
Each iSSA model included movement covariates, including step-length, its natural logarithm, and the angular deviation (cosine of the turn angle), as well as all habitat covariates mentioned above. Models were estimated using conditional logistic regression in the R package “survival” [32, 33]. The importance of habitat selection to explain space use was determined by comparing iSSA models containing only the movement covariates with models containing both movement and habitat covariates by means of the Akaike Information Criterion (AIC). In order to estimate mean step length (lmean) while accounting for habitat selection, we combined the estimated coefficients of the step-length with the parameter estimates of the Gamma distributions (used for sampling available step-lengths) as follows [4]:
$$ {l}_{mean}=\frac{k+\beta \ln (l)}{\theta^{-1}-{\beta}_l}, $$
where k and θ are the shape and scale of the observed Gamma distribution, respectively, and β
l
and β
ln(l)
are the iSSA coefficients for the observed step length and natural-log step length, respectively.
Because selection coefficients are not explicit about the range and values of used habitat, and might be dependent on other habitat covariates, we used the relative selection strength (RSS) [34] to show and interpret habitat selection results. RSS is defined as the probability of habitat use in one location over other locations. To obtain “population” (as defined by a group of interest) rather than individual RSS, we bootstrapped the mean RSS and calculated population 95% confidence intervals (see a similar approach in [35]).
Home range estimation
Home range size was estimated by calculating the area used by the different tracked animals through a bivariate normal utilization-kernel using the R package “adehabitatHR” [36]. We used the reference smoothing parameter (h_ref) [37] to estimate home ranges and the respective 90% and 50% contours (percentages chosen based on [38]), the latter representing an estimate of core areas. Home ranges were estimated for each animal in each season.
Determinants of home range size
Linear mixed-effects models (LMM) were used to test how home range size varies within and across seasons and what drives this variation. We hypothesised that season, selection-independent movement, selection strength, resource availability, and sex drive home range size. Selection-independent movement was characterised by the mean step length estimated from the iSSAs and the iSSA movement coefficient for the turn angle, which represents a measure of directional persistence (i.e. the concentration parameter of a von Mises distribution; [4]). Home range size was expected to increase with both step length (i.e. displacement rate) and directional persistence. iSSA selection coefficients were used as selection strength covariates. For resource availability, we used several proxies that included altitude and population nucleus as well as the mean and coefficient of variation (CV) of habitat variables related to forage availability (heat load and NDVI). It is worth noting that ibex were found to select a defined range of resource values, as consistently significant unimodal relationships were predicted by the iSSA models (see Results); therefore, higher CV (i.e. lower resource density) meant less resource availability.
The log-transformed home range size was used as the response variable, and the predictors included season (four-level factor), population nucleus (two-level factor), altitude, mean step length, directional persistence, slope selection (the linear and quadratic selection coefficients), heat load selection (linear and quadratic), NDVI selection (linear and quadratic), heat load availability (mean and CV), NDVI availability (mean and CV), and sex (two-level factor). Ibex identification was used as a random intercept effect. We performed model selection by comparing all models that included all possible predictor combinations by means of corrected AIC (AICc). Model estimation and selection were performed with the R packages “lme4” [39] and MuMIn [40], respectively. The relative importance of the different home range drivers was assessed by the difference in AICc when removing the target group of predictors.