Population and data

The GPNP is located in the northwestern Italian Alps (45°25‘N, 7°34‘E). The wardens of the GPNP perform two ibex counts every year; one in late spring and one in the autumn. They count the population by walking over established routes. The numbers of ibex in the two counts are highly correlated (von Hardenberg et al. 2000); also, the autumn counts include all summer newborns (discounted for neonatal mortality). Therefore, we used the autumn counts for our analysis. Details on the count techniques used to obtain the data are provided in Appendix A of Jacobson et al. (2004). In the same Appendix, on the basis of the correlation between the two count series, the authors suggest that counts are reliable.

The replication of counts within the same year, together with identification of outliers, has been suggested by Largo et al. (2008) as a methodology to make the ibex counts more reliable and to avoid huge underestimates of the population size. Largo et al. (2008) define the data of years with a growth rate > 1.35 as ‘obvious outliers’. In the GPNP case, this value was never exceeded. Nevertheless, the minimum and maximum growth rate occurred in two consecutive years; in 1976 (N₁₉₇₇/N₁₉₇₆ = 0.64, N_t being the autumn population density of year t) and in 1977 (N₁₉₇₈/N₁₉₇₇ = 1.28). Counts of 1977 are then candidate to be a notable underestimate of the actual population size. We therefore compared these growth rate values with those obtained using spring counts of year t () which are / = 0.61 in 1976 and / = 1.25 in 1977. Thus, the growth rate values are comparable using the two different counts, and we therefore decided to include 1977's data in the analysis. Moreover, 1976 was the year characterised by the largest snow depth in the time series, so that a very low growth rate was not surprising.

Data of the ibex population in the GPNP are available from 1956 until today and constitute one of the longest continuous existing time series of mountain ungulate counts worldwide. In addition, the data are structured in age and sex classes, which, for clarity, hereafter, we will call population compartments. Figure 1A shows the total number of ibex (autumn counts) while Figure 1B reports the counts of each of four population compartments: kids (0-1 years), yearlings (1-2 years), adult (> 2 years) females and males.

Figure 1A reports the mean winter snow depth (S_t) measured at the IREN ENERGIA meteorological station at Lake Serrù (45°16‘N, 7°8‘E, 2,240 m a.s.l.), averaged from November of year t to May of year t + 1. In the threshold models by Jacobson et al. (2004), the critical value of S̄ = 154 cm is chosen as the average plus half a standard deviation of the snow time series in the period of 1961-2000.

From 1961 to 1982, the total number of ibex in GPNP fluctuated around a mean value of about 3,300 individuals. Then, an interesting temporal pattern clearly emerged: the total abundance displayed a rapid increase (from ca 3,250 individuals in 1982 to almost 5,000 in 1993) and then a sharp decrease, with only ca 2,700 ibex surviving in 2008. This unimodal variation of the population size occurred under a monotonically decreasing trend of snow depth, somehow suggesting that the interaction between climatic conditions and population numbers is not at all as simple as one might expect.

Unstructured and sex- and age-structured models

The use of thresholds (see Jacobson et al. 2004, Corani & Gatto 2007) represents a form of non-linearity in the demographic growth rate model. We wished to extend this approach to better capture the non-linearities of the system. To this end, we considered continuous models that incorporate all the quadratic terms in the explanatory variables. In its complete form, a model of this type has the polynomial expression

The second extension formulated here is the introduction of sex- and a simplified age-structure. Since adult counts were not partitioned into yearly age classes, the models we propose are not fully age-structured; however, they provide a first attempt toward the bottom-up approach to ibex demography invoked by Jacobson et al. (2006) and Yoccoz & Gaillard (2006). According to the counts data reported in Figure 1, in each year t, the population was partitioned into four compartments: kids K_t, yearlings Y_t, adult males M_t and adult females F_t. Juvenile classes K_t and Y_t included animals of both sexes that were not yet reproductive. At year t, the rate of demographic variation of kids (σ_K,t), males (yearlings and adults,σ_M,t) and females (yearlings and adults, σ_F,t), as well as the weaning success w_t were defined as

By using the same approach introduced above for the unstructured counterparts, structured models were defined through non-linear relationships between the logarithms of rates (2) and two covariates: the average snow depth S_t and the density D_C,t of animals that was relevant for intraspecific competition within each compartment C (see below). As in the unstructured case, such non-linear relationships can be continuous, i.e.

Apart from the short breeding season, Alpine ibex live in groups spatially segregated by age and/or sex (Villaret et al. 1997). Therefore, in addition to accounting for dependence on total density D_C,t = D_W,t = N_t, we also considered the case of sexually-segregated density dependence. Although not all male yearlings depart from their mothers, we considered the simple hypothesis of two separate groups, one with ‘males’ (adult males and male yearlings) and the other with ‘females’ (mothers with all kids and female yearlings). For equations 3 and 4 above, this translates into considering D_W,t = D_F,t = D_K,t = F_t + K_t + Y_t/2 and D_M,t = M_t + Y_t/2. Note that equation 4 is dynamic, in the sense that weaning success in the spring of year t + 1 is dependent on the state of the population in the autumn of year t and average snow depth in the winter between year t and year t + 1. In fact, meteorological conditions and population density of the year preceding births can describe the physiological state of mothers during the reproduction period well (see the ‘tap-hypothesis’ in Grøtan et al. 2008).

Model selection criteria

The number of candidate models that emerged from the scheme outlined in the previous section was huge. In fact, we defined two families of models: unstructured and sex- and age-structured. Every unstructured model was fully characterised by the function relating the total population growth rate (log(N_t+1/N_t,)) to the different covariates. Every structured model instead required the mathematical description of four demographic quantities (log(σ_K,t), log(σ_M,t), log(σ_F,t) and log(w_t)) that can be either related to the density of the entire population (no spatial segregation) or to the compartment-specific density of ‘males’ or ‘females’, as defined in the previous section (spatial segregation). Each of the above rates could be related to the covariates by either a continuous formulation (like in equations 1, 3 and 4) or a threshold formulation. Independent of the chosen formulation, candidate models included all possible combinations obtained from the most complex formulation in which one or more terms were dropped, with the only exception of the constant term β₀. Simple combinatorics revealed that there were 32 continuous formulations and eight threshold formulations for the suite of unstructured models. For the structured case, because of the two different hypotheses on density dependence (with or without spatial segregation) and considering that some of the models depend on snow only, we had a total of 74 potential model candidates for each of the four demographic parameters.

To orient ourselves inside such a large dimensional space of candidates, we used standard model selection techniques: the second-order Akaike's Information Criterion (AIC_c; see Sugiura 1978), the Bayesian Information Criterion (BIC; see Schwarz 1978) and a criterion based on Structural Risk Minimisation (SRM; see Cherkassky et al. 1999). Explaining in detail the logic and the methodology that underlie each of the above criteria goes beyond the purpose of our study; the interested reader might refer to Burnham & Anderson (2002) for AIC_c and BIC, and to Corani & Gatto (2007) for SRM. Here, it suffices to mention that all used criteria are apt at analysing short temporal data series, as is the case for most demographic time series like that of ibex counts in GPNP. For each selection method, the ‘best’ model is the one that minimises the value of the criterion. We used multiple criteria rather than a unique selection method, because their different ways of accounting for model parsimony can produce a difference in the results which is worth investigation.

In general, model selection does not provide one winner model, but a hierarchical set of optimal models. The number of these optimal models vary between criteria and is not known a priori. For AIC_c, Richards (2005) suggests for example to include in the set of optimal models those candidates whose AIC_c's differ from the minimum value by < 4. The most parsimonious among these models (i.e. the one with the smallest number of parameters) might be chosen as the very best option, but other choices within that set are also possible. For BIC, a similar selection rule is available (ΔBIC < 2; Raftery 1995), while to date no rule has been proposed for the SRM criterion. Here we used a rather arbitrary cut-off threshold and selected as optimal those models where SRMs differed by < 6% from the minimum (SRM_best). In our study, we thus qualified as best models those with the lowest number of estimated parameters among the models that met all these conditions: ΔAIC_c < 4, ΔBIC < 2 and SRM < 1.06SRM_best. In principle, it might be possible that no model satisfies the three conditions simultaneously, but this was not the case in our study. For the best models, we also calculated the adjusted R².

Performances of models

To evaluate errors of parameter estimates, we used bootstrap (Efron 1979). The bootstrap method provides an unbiased estimation of each parameter (but see caveats in Efron et al. 1993) and an estimation of their variances. The statistics reported in Appendix I were obtained by calibrating the parameter values over 1,000 bootstrap samples, each consisting of n data values drawn with replacement from the original n-sized data set.

In order to assess the predictive ability of the best models over longer time scales, we repeated the same parameter tuning and simulation experiments performed in Jacobson et al. (2004). To compare our results with theirs, we recalibrated all parameters of the selected models using only the first 20 years of data (1961-1980). Then, based on the recalibration, we simulated the ibex population trends until 2005 with both unstructured and structured models.

Simulating population trends after 1980 with the unstructured models was quite simple, because it was sufficient to initialise the systems with N₁₉₈₁ and use the snow depth series S₁₉₈₁, S₁₉₈₂, ..., S₂₀₀₄ as inputs. Predicting the total population abundance with structured models was more cumbersome, because the different rates of demographic variations and the weaning success of the structured models must be aggregated into one global model.

The main advantage of the structured models is to provide information on the role played by covariates in each specific demographic process (i.e. adult rates of demographic variation, kid rate of demographic variation, weaning success). Also, structured models can be used to understand whether simulations of any particular population compartment deviated from data significantly more than others, thus pointing out the weakest links of the chain in the global model. To specifically investigate this point, we performed additional simulations with global models in which part of the state variables were computed via the model equations, while others were directly equated to data. More precisely, we predicted the number of adult males and females from yearling counts using the following equations

The estimated total population abundance N̂_t can be derived by simply summing the abundances of all compartments. A direct comparison of the long-term prediction obtained from structured vs unstructured models was then possible. As an index of predictive ability, we used the root mean square error between N̂_t and N_t:

Predictions

By using the selection criteria described in the previous section, we obtained the best models reported in Table 1. Within the family of unstructured models, the two selected systems are both discontinuous, which is in agreement with the original findings of Jacobson et al. (2004). NT1 includes pure density dependence and the interaction term S_tN_t, while NT2 includes pure snow-dependence and the interaction term. This latter model has the same structure as the model selected by Corani & Gatto (2007). Interestingly enough, the inclusion of higher-order terms in the continuous systems (equation 1) for the total growth rate did not result in better performance of the unstructured models.

The picture emerging from the analysis of structured models was more complex. First, we noticed that the set of best models among the myriad of potential candidates was indeed very small (see Table 1). Most interesting was the fact that the best density-dependent models selected by our procedure tended to be those with spatial segregation, i.e. the dynamics of ‘males’ was influenced more by ‘males’ than by the entire population and the same was true for ‘females’. The only exception was the rate of demographic variation of kids, for which both a model incorporating the segregation hypothesis (KC1, in which D_t = ‘females’) and one excluding segregation (KC2, in which D_t = total population density) passed the model selection. As for the rate of demographic variation of males, three models that did not incorporate the spatial segregation hypothesis would also satisfy the Δ's criterion described above, but were excluded because they were not as parsimonious as the others. While for females and kids two alternative models were selected, the rate of demographic variation of males and the weaning success had a unique best functional formulation.

For structured models, there was no systematic prevalence of threshold over continuous formulations. More precisely, while threshold models are best for adult compartments (male and female rates of demographic variation), continuous models were selected for kid demography and weaning success. These included the linear and quadratic terms in the snow depth S_t and the interaction term D_tS_t. A closer look at the signs of the best estimates of parameter values (see Appendix I: Tables 6, 7 and 8) revealed that the coefficients β₂ multiplying S_t were positive while the coefficients β₄ multiplying were negative, an indication that kid demography and weaning success had a non-monotonic dependence on snow depth (Fig. 2). This suggests that years characterised by particularly low snow depth can be detrimental to the juvenile compartments of the Alpine ibex, a result that is in keeping with the recently observed drop of the relevant rates (see von Hardenberg et al. 2009).

In terms of adjusted R², the best performances were provided by models for the total population and for the rate of demographic variation of females, while the models for male and kid demography and weaning success displayed a poorer fit to the data.

Most coefficients of variation of the estimated parameters, as shown in Appendix I (Tables 1-7) were of the order of 10⁻¹ (see also the Discussion), thus showing that best fits were rather robust. Using the bootstrapped parameter distributions, we assessed the predictive ability of the best models under parameter uncertainty. First, we performed one-step-ahead predictions, whose distribution (5th-95th percentiles) was obtained from the 1,000 parameters of the bootstrap analysis. The result is shown in Figure 3. Despite all models being calibrated on data over the entire period 1961-2004, the predictive ability deteriorated at the end of the 1970s. In fact, while data fell within the prediction range during the first part of the time series, deviations of predictions from data were more frequent after the beginning of the 1980s. Particularly evident was the mismatch in the case of the rate of demographic variation of kids.

Long-term simulations

Having selected two best models for the rates of demographic variation of females and kids and one best model for the other rates, we obtained four structured global models, namely STR1 (consisting of models MT, FT1, KC1 and WC), STR2 (MT, FT2, KC1 and WC), STR3 (MT, FT1, KC2 and WC) and STR4 (MT, FT2, KC2 and WC). Similarly to the case of unstructured models, the initial condition was considered as known X̂₁₉₈₁ = X₁₉₈₁ for all X's) as well as the snow depth time series S_t over the entire simulation horizon (from 1981 to 2004). Figure 4 shows the distributions of values (in terms 5th-95th interpercentiles and 25th-75th interquartiles) simulated with the best models under parameter uncertainty (evaluated via bootstrap). The plots revealed that the unstructured models NT1 and NT2 reproduced the recent trends in quite a similar way. The main difference was that model NT2 exhibited more oscillations than NT1. In terms of RMSE_N, the performances of the two unstructured models were comparable (RMSE_N = 506 for NT1 and 531 for NT2), and they were both significantly more effective than the model by Jacobson et al. (2004; RMSE_N = 852) which did not include the double estimate of β₀.

The structured global models did not perform qualitatively better than the unstructured model, or provide quantitatively significant improvements (RMSE_N = 577 for STR1, 570 for STR2, 524 for STR3 and 520 for STR4). As in the case of unstructured models, trends simulated with structured models underestimated the peak and overestimated the abundances of recent years (see Fig. 4).

The temporal evolution of RMSE_N(k) for the best structured and unstructured models is shown in Figure 5. Structured models performed better than unstructured ones in terms of RMSE_N(k) along almost the entire simulation period, with the exception of the first and last years. On the other hand, recent data were systematically included between the 5th and 95th percentiles of all simulations obtained with unstructured models (see Fig. 4), while this was not true for simulations obtained with structured models. However, unstructured models displayed very high variability of long-term predictions, simulating unrealistic abundances as high as 10,000 individuals or more.

Even though structured models did not significantly improve the prediction of total population numbers, it was useful to explore the contribution of the different compartments to the simulated dynamics. The simulated numbers of adult males and females obtained from the models in equation 5, while considering the time series of yearlings Y_t as a known input, are shown in Figure 6A-C. The two best models FT1 and FT2 performed rather similarly and both underestimated the population peak of the 1990s but reproduced fairly well the subsequent decreasing trend. The model for males was more precise than the two models for females. Observations were almost always included between the 5th and 95th percentiles of bootstrapped simulations for both males and females. Simulations for the adult male compartment were more variable than those for females. This was expected because the coefficients of variation of the estimated parameters of model MT were higher than unity when snow depth was below threshold (see Appendix I: Table 5). Data in Figure 6D-F show the juvenile compartments of kids and yearlings and were obtained with equations 6. While the counts data were often included in the prediction range before and during the population peak, they were frequently below the 5th percentile starting from the mid-1990s, especially in the case of yearlings.

Discussion and conclusions

The models proposed and analysed in our study relate the total annual growth rate (unstructured models) or the rates of demographic variation and weaning success (structured models) to population density and snow depth. For the unstructured models, the idea is not new, but here we accounted for richer functional forms than those available in the literature. As for the novel compartment-structured approach, we proposed models aimed at including the likeliest factors impacting on the juvenile and adult ibex rates of demographic variation and on mother weaning success.

The sexually-segregated density-dependence hypothesis almost systematically improved the model performance and was specifically selected for the female and male rates of demographic variation and for the weaning success. For kids, one of the two selected models was based on sexual segregation, while the other was not. With the only exception of one model for adult female demography, all the structured models selected among the many potential candidates included the interaction term between snow depth and density (S_tD_t). This result reinforced the evidence that climate strongly intensifies density dependence of ibex in GPNP, in accordance with the results of Jacobson et al. (2004). In fact, not only the total annual growth rate, but every single vital rate appeared to be crucially dependent on the joint effect of snow depth and population density.

In contrast to what happened for the rate of demographic variation of adult males and females, the best models for kids and weaning success were continuous rather than threshold-like. Also, they revealed that an intermediate value of snow depth was optimal for both demographic parameters. This non-monotonic snow dependence is presumably related to different biological mechanisms. Winters with high snow cover are detrimental to all ibex, and especially adults, because food is scarcer than in low-snow winters and more energy is required to move in the high snowpack and dig out the dry grasses from below the snow. High ibex densities amplify the snow effect as in this case intraspecific competition becomes more severe. In addition, winters with large snow cover have a correspondingly higher probability of avalanches and the associated ibex casualties. This mechanism, where high winter snow has negative effects on ibex dynamics, has been thoroughly discussed in Jacobson et al. (2004).

Here, however, another effect was discovered. Winters with very low snow cover (such as the winters in the last 20 years) were also detrimental to ibex, but this time through their effects on kid survival and weaning success. This seems to reflect a major sensitivity of juveniles to a lack of snow during winter. Pettorelli et al. (2007) have shown that, possibly because of climate change, the green-up of GPNP vegetation has become faster. Indeed, the annual maximal increase in the normalised difference vegetation index (NDVI), a satellite-based measure that is strongly correlated with the net primary productivity, appeared to increase over time (Pettorelli et al. 2007). This may lead to a shorter period of availability of high-quality forage over a large spatial scale, decreasing the opportunity for mountain ungulates to exploit high-quality forage. The rate of demographic variation in kids might thus be influenced either directly or possibly via the state of mothers during lactation.

Although simulations obtained with the best structured and unstructured models showed quite high variability within the range of bootstrapped parameter values, the actual time series of animal counts was not always included in that range. Reference simulations made using unbiased parameters (see Fig. 4) showed that all models underestimated the population peak (occurring from the mid-1980s to the mid-1990s) and overestimated recent counts.

The long-term simulations of each single compartment showed an underestimation of the adult compartments during the growing phase of the peak and an overestimation of the juvenile compartments (kids and yearlings) during the population decline which started from the mid-1990s. These results indicate that weaknesses in modelling adult rates of demographic variation are responsible for the underestimation of the population peak in the 1990s, while the overestimation of the recent declining trend can be mainly ascribed to inadequate modelling of recruitment and rate of demographic variation of kids.

Our findings suggest that mechanisms other than direct climate effects and population density could influence the dynamics of ibex in GPNP. As mentioned above, the recent drop in kid survival might also be related to the state of pastures. Other important, yet poorly explored factors, are parasite infections and interspecific competition. In fact, emergence of parasitic infections that critically affect the demography of ungulates have been recorded in different Arctic populations (Kutz et al. 2004). Such epidemics appear to be favoured by climate warming. Modelling their effects on ibex population dynamics is thus a promising avenue for investigation (see for example Ferrari et al. 2010). As for interspecific competition, it is known that ibex share their habitat with chamois Rupicapra rupicapra in the GPNP. Actually, the populations of chamois and ibex showed similar trends until 1993 (Picollo 2002). After that year, however, chamois stayed approximately constant, while ibex started to decrease. Therefore, testing whether competition between the two species (documented in Pfeffer & Settimo 1973) might be responsible for the decrease of ibex is worth exploring.

Finally, a big step ahead and a possible remedy to the shortcomings of currently available models would consist in formulating age-structured models that consider a more realistic subdivision into yearly classes of individuals of the same sex. Models at this finer scale would include the effects of senescence on male and female rates of demographic variation, which has been clearly documented, for instance in a French population of ibex (Toïgo et al. 2007). As reported for other species of large herbivores (Gaillard et al. 2000), it is plausible that old ibex suffer much more than young individuals in years characterised by particularly unfavourable environmental conditions. Compartmented models, like those proposed in our present study, cannot account for animal senescence, because all adults are included in the same class independent of their age. To develop truly age-structured models, one could integrate the count data with capture-mark-resight data collected on cohorts of individually marked animals (currently under way at GPNP) and use inversion techniques such as those based on the Kalman filter approach (for an example on Soay sheep Ovis aries see Tavecchia et al. 2009). In any case, the long and unique time series of population counts at GPNP will again provide an excellent empirical database against which to test different hypotheses on the mechanisms at work.