Translator Disclaimer
1 September 2012 Sex- and age-structured models for Alpine ibex Capra ibex ibex population dynamics
Andrea Mignatti, Renato Casagrandi, Antonello Provenzale, Achaz von Hardenberg, Marino Gatto
Author Affiliations +

Past analysis has shown that the population dynamics of Alpine ibex Capra ibex ibex are regulated by both population density and winter snow accumulation. However, recent time series of the ibex counts in the Gran Paradiso National Park, Italy, show interesting trends in comparison with historical snow data: while the winter snow depth has steadily decreased since the beginning of the 1980s, the ibex population experienced rapid growth during the 1980s and the early 1990s, followed by a strong decrease. To explain these dynamics, we built novel age-structured population models in which demographic parameters depended on density and snow depth. They included a non-monotonic effect of snow depth and density on the vital rates, the age and sex structure of the population, and spatial segregation between females and males. Using information criteria (AICc, BIC and SRM), we selected the best models and found that: 1) snow and density interacted in determining the demography of all population sex and age classes, thus confirming that unfavourable climatic conditions intensified the density dependence of the population, 2) the effect of snow was non-monotonic on weaning success and rate of demographic variation of kids, which were maximal for intermediate snow depths, and 3) accounting for spatial segregation between sexes improved the fitting of the models, which suggests that the different use of space made by males and females influenced intraspecific competition. When the selected models were recalibrated using past data and used to simulate recent trends, they underestimated both the rapid growth of the 1980s-1990s and the recent decline of the population. Using the novel sex- and age-structured models, we found that the underestimation of the peak was mainly due to deficiencies of adult demography models, while the overestimation of the recent population abundance was caused by shortcomings in the models of recruitment.

The population dynamics of ungulates is strongly determined by both density dependence and environmental drivers (Lande 1993, Post et al. 1997, Sæther 1997, Forchhammer et al. 2002), which may operate in an interacting way (Gaillard & Yoccoz 2003). Jacobson et al. (2004) were first in showing that the population dynamics of Alpine ibex Capra ibex ibex is significantly affected by a specific en-vironmental variable, i.e. the average snow depth during winter. Since Alpine ecosystems are extremely sensitive to climate change (Fischlin et al. 2007), and snow depth has been decreasing in the Alps in recent years (e.g. Terzago et al. 2010), studying how population dynamics of high-altitude species is influenced by abiotic disturbances and trends is of paramount importance.

The exceptionally long time series of Alpine ibex counts in the Gran Paradiso National Park (GPNP), Italy, provides a unique opportunity to study the complex interplay between population density and climatic conditions operating in the species dynamics. To reproduce the dynamics of the population and to possibly make predictions into its near future, different models have been proposed. The simple but powerful approach put forth by Jacobson et al. (2004) explored various possible relationships of total ibex population increase with both total animal density and snow depth. The analysis performed by Jacobson et al. (2004) revealed that the population rate of increase is significantly different in years characterised by low vs high snow depth. This observation led the authors to include ‘threshold-models’ in their model suite, in which different functional relationships are used for years of high and low snow depth, respectively. The process of parameter calibration and model selection actually picked up two threshold models as the best candidates. Both models acceptably reproduce the increase of the 1980s and the population peak occurring at the beginning of the 1990s, while the subsequent decline is only partially captured (see Fig. 3 in Jacobson et al. 2004). Unfortunately, if used to simulate more recent trends emerging from newly available data, such models largely overestimate the total abundance of ibex (see also the last three simulated years in Jacobson et al. 2004).

The models proposed by Jacobson et al. (2004) have been expanded in various ways. For example, Corani & Gatto (2007) found that a slightly modified version of Jacobson et al.'s threshold models, obtained by introducing two different values of the intrinsic rate of increase at low and high snow depth, had better performances in terms of model selection criteria. Other authors tried to avoid the use of thresholds and proposed smoother functional forms for the rate of increase. Bianchi et al. (2006), for example, used local linear models instead of a piecewise linear system, while Lima & Berryman (2006) proposed non-parametric, non-linear versions of the dynamic model using the so-called Generalised Additive Modeling (GAM) approach of Hastie & Tibshirani (1990). Although these models provided technical insights into the potential limitations of the approach followed in Jacobson et al. (2004), none of them altered the main biological assumptions of the original study significantly or was able to reproduce the population dynamics during the low-density phase in the last 10 years.

In our paper, we follow the lines suggested by Yoccoz & Gaillard (2006) and improve over existing models by incorporating peculiarities of the ibex life cycle that have never been included in previous modelling attempts, despite their ecological importance (Jacobson et al. 2006). The most evident among these characteristics are sexual dimorphism and age structure. Also, we contrast the threshold models by Jacobson et al. (2004) and Corani & Gatto (2007) with differentiable models (hereafter called continuous models) in which we use second-order polynomial approximations of non-linear unknown functional forms. The effects of sex and age are quite strong in Alpine ibex and must be taken into account to understand how density and environmental drivers might affect the various fitness components of ibex populations (Gaillard et al. 1998, 2000). According to sex and age, animals live in spatially segregated groups, using different habitat types. Males and females usually only join during the short breeding season from mid-November to mid-January (Nievergelt 1974). Kids live with their mothers throughout their first year of life while male yearlings gradually depart from female groups and form bachelor groups (2-3 years old) that join adult males (Villaret & Bon 1995). Recruitment and juvenile survival are usually considered more sensitive to both density and environmental variability than adult survival. Among adults, female survival is larger and more buffered against limiting factors than male survival (Gaillard et al. 1998, Toïgo & Gaillard 2003, Toïgo et al. 2007), although the difference is perhaps smaller for ibex than for other large mammals (Toïgo et al. 1997). Since males and females (with kids) are almost always spatially segregated, it is reasonable to imagine that intraspecific competition can involve members of the same sex only. In our present work, this hypothesis is contrasted against the alternative that intraspecific competition occurs among all individuals in the population.

Data for kids, yearlings, adult males and adult females are actually available, although they have never been used to formulate structured population models (Appendix E of Jacobson et al. 2004). Here, we fill the gap and use this detailed information to build models for the rate of demographic variation in kids and adults (both males and females) and for the weaning success. These sex- and age-structured models not only allow a demographic projection in time of specific population compartments, but can also reveal which demographic processes are most influenced by biotic (e.g. population density) or abi-otic (snow) factors.

Material and methods

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 (N1977/N1976 = 0.64, Nt being the autumn population density of year t) and in 1977 (N1978/N1977 = 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 (i0909-6396-18-3-318-ilm01.gif) which are i0909-6396-18-3-318-ilm02.gif/i0909-6396-18-3-318-ilm03.gif = 0.61 in 1976 and i0909-6396-18-3-318-ilm04.gif/i0909-6396-18-3-318-ilm05.gif = 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 (St) 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

where Nt is the total population density, St is the mean winter snow depth and ρt is a stochastic factor representing environmental noise and unmodelled processes. In previous models, the parameters β4 and β5 were set to zero. Continuous models of this kind can be viewed as a second-order Taylor expansion of a more general model log(Nt+1/Nt) = f(Nt,St) around the average ibex density and snow depth.

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 Kt, yearlings Yt, adult males Mt and adult females Ft. Juvenile classes Kt and Yt 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 wt were defined as

Definitions of the rates of demographic variation σM,t and σF,t implicitly entail the assumptions of balanced sex ratio at birth and no ratio distortion in the juveniles (kids and yearlings). For the realism of these assumptions, see Stüwe & Grodinsky (1987) and Toïgo et al. (1997).

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 St and the density DC,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.

with C equal to K, F or M, and
or, alternatively, they may be discontinuous. In this latter case, they contain fewer terms (β4 = β5 = 0), but the parameters β0, β1, β2 and β3 can take different values in years characterised by snow depths above or below the critical threshold S̄

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 DC,t = DW,t = Nt, 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 DW,t = DF,t = DK,t = Ft + Kt + Yt/2 and DM,t = Mt + Yt/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(Nt+1/Nt,)) 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(wt)) 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 β0. 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 (AICc; 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 AICc 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 AICc, Richards (2005) suggests for example to include in the set of optimal models those candidates whose AICc'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 (SRMbest). 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: ΔAICc < 4, ΔBIC < 2 and SRM < 1.06SRMbest. 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 R2.

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 N1981 and use the snow depth series S1981, S1982, ..., S2004 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

where Yt is the measured number of yearlings in year t, while M̂ and F̂ are the model-predicted numbers of adult males and females, respectively, and i0909-6396-18-3-318-ilm06.gif and i0909-6396-18-3-318-ilm07.gif are the best-estimated rates of demographic variations. Similarly, we predicted the number of kids from the mother counts and the number of yearlings from the kid counts using
where Ŵt and i0909-6396-18-3-318-ilm08.gif are the best estimates of weaning success and rate of demographic variation of kids, respectively.

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 Nt:

To better evaluate the performance of the different models throughout the simulation period, it was also convenient to define the root square error averaged from the beginning of simulation (year 1981) to year k, i.e.



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 StNt, 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 Dt = ‘females’) and one excluding segregation (KC2, in which Dt = 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 St and the interaction term DtSt. 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 β2 multiplying St were positive while the coefficients β4 multiplying i0909-6396-18-3-318-ilm09.gif 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 R2, 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−1 (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̂1981 = X1981 for all X's) as well as the snow depth time series St 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 RMSEN, the performances of the two unstructured models were comparable (RMSEN = 506 for NT1 and 531 for NT2), and they were both significantly more effective than the model by Jacobson et al. (2004; RMSEN = 852) which did not include the double estimate of β0.

The structured global models did not perform qualitatively better than the unstructured model, or provide quantitatively significant improvements (RMSEN = 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 RMSEN(k) for the best structured and unstructured models is shown in Figure 5. Structured models performed better than unstructured ones in terms of RMSEN(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 Yt 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 (StDt). 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.


this work could not have been performed without the dedication of the GPNP park wardens who have been collecting the ibex data for > 50 years. To them, we express our gratitude. Part of this work has been performed during Andrea Mignatti's PhD, supported by Fondazione Lombardia per l'Ambiente (project SHARE-Stelvio). The meteorological data were kindly provided by the IREN ENERGIA SPA, Torino, Italy. We also thank two anonymous reviewers and technical editor Jan Drachmann for the usefull comments and suggestions, which significantly contributed to improving the quality of the publication.



M Bianchi G Corani G Guariso& C Pinto 2006: Prediction of ungulates abundance through local linear algorithms. - Environmental Modelling and Software 21: 2004–2007. Google Scholar


K Burnham& D Anderson 2002: Model selection and multimodel inference: a practical information-theoretic approach. Springer Verlag, New York, New York, USA, pp. Google Scholar


V Cherkassky X Shao F Mulier& V Vapnik 1999: Model complexity control for regression using VC generalization bounds. - IEEE Transactions on Neural Networks 10(5): 1075–1089. Google Scholar


G Corani& M Gatto 2007: Structural risk minimization: a robust method for density-dependence detection and model selection. - Ecography 30(3): 400–416. Google Scholar


B Efron 1979: Bootstrap methods: another look at the jackknife. - The Annals of Statistics 7(1): 1–26. Google Scholar


B Efron& R Tibshirani 1993: An introduction to the bootstrap. Chapman & Hall/CRC, New York, New York, USA, pp. Google Scholar


N Ferrari R Rosà P Lanfranchi& K.E Ruckstuhl 2010: Effect of sexual segregation on host-parasite interaction: model simulation for abomasal parasite dynamics in alpine ibex (Capra ibex). - International Journal for Parasitology 40(11): 1285–1293. Google Scholar


A Fischlin G.F Midgley J.J Price R Leemans B Gopal C Turley M.D.A Rounsevell O.P Dube J Tarazona& A.A Velichko 2007: Ecosystems, their properties, goods and services. In: M.L Parry O.F Canziani J.P Palutikof P.J van der Linden& C.E Hanson (Eds.); Climate change 2007: Impacts, Adaptation and Vulnerability. Contribution of Working Group II to the 4th Assessment Report of the Intergovernmental Panel on Climate Change. Cambridge University Press, Cambridge, UK, pp. 211–272. Google Scholar


M Forchhammer E Post N Stenseth& D Boertmann 2002: Global climate change and phenotypic variation among red deer cohorts. - Population Ecology 44: 113–120. Google Scholar


J Gaillard M Festa-Bianchet& N Yoccoz 1998: Population dynamics of large herbivores: variable recruitment with constant adult survival. - Trends in Ecology & Evolution 13(2): 58–63. Google Scholar


J Gaillard M Festa-Bianchet N Yoccoz A Loison& C Toïgo 2000: Temporal variation in fitness components and population dynamics of large herbivores. - Annual Review of Ecology and Systematics 31: 367–393. Google Scholar


J Gaillard& N Yoccoz 2003: Temporal variation in survival of mammals: a case of environmental canalization? - Ecology 84(12): 3294–3306. Google Scholar


V Grøtan B-E Sæther F Filli& S Engen 2008: Effects of climate on population fluctuations of ibex. - Global Change Biology 14(2): 218–228. Google Scholar


T Hastie& R Tibshirani 1990: Generalized additive models. Chapman & Hall/CRC, New York, New York, USA, pp. Google Scholar


A Jacobson M Festa-Bianchet A Provenzale A von Hardenberg& B Bassano 2006: Comment on Lima & Berryman (2006): the Alpine ibex revisited. - Climate Research 32(2): 137. Google Scholar


A Jacobson A Provenzale A von Hardenberg B Bassano& M Festa-Bianchet 2004: Climate forcing and density dependence in a mountain ungulate population. - Ecology 85(6): 1598–1610. Google Scholar


S Kutz E Hoberg J Nagy L Polley& B Elkin 2004: “Emerging” Parasitic Infections in Arctic Ungulates. - Integrative and Comparative Biology 44(2): 109–118. Google Scholar


R Lande 1993: Risks of Population Extinction from Demographic and Environmental Stochasticity and Random Catastrophes. - Global Change Biology 142(6): 911–927. Google Scholar


E Largo J Gaillard& M Festa-Bianchet 2008: Can ground counts reliably monitor ibex Capra ibex populations? - Wildlife Biology 14(4): 489–499. Google Scholar


M Lima& A Berryman 2006: Predicting nonlinear and non-additive effects of climate: the Alpine ibex revisited. - Climate Research 32(2): 129. Google Scholar


B Nievergelt 1974: A comparison of rutting behaviour in the Ethiopian and Alpine ibex. In: V Geist& F Walther (Eds.); The Behaviour of Ungulates and its relations to Management. IUCN Publications new series no 24, Morges, Schwitzerland, pp. 324–340. Google Scholar


N Pettorelli F Pelletier A von Hardenberg M Festa-Bianchet& S.D Côté 2007: Early onset of vegetation growth vs. rapid green-up: impacts on juvenile mountain ungulates. - Ecology 88(2): 381–390. Google Scholar


P Pfeffer& R Settimo 1973: Déplacements saisonniers et compétition vitale entre mouflons, chamois et bouquetins dans la Réserve du Mercantour (Alpes Maritimes). - Mammalia 37(2): 203–219. (In French). Google Scholar


M Picollo 2002: Effetti della variabilitá climatica su ecosistemi alpini. - M.Sc. thesis, Facoltá di Scienze Matematiche, Fisiche e Naturali, Universitá degli studi di Torino, Torino, Italy, pp. (In Italian). Google Scholar


E Post N Stenseth R Langvatn& J Fromentin 1997: Global climate change and phenotypic variation among red deer cohorts. - Proceedings of the Royal Society, Series B 264: 1317–1324. Google Scholar


A Raftery 1995: Bayesian model selection in social research. - Sociological Methodology 25(1): 111–163. Google Scholar


S Richards 2005: Testing ecological theory using the information-theoretic approach: Examples and cautionary results. - Ecology 86(10): 2805–2814. Google Scholar


B-E Sæther 1997: Environmental stochasticity and population dynamics of large herbivores: a search for mechanisms. - Trends in Ecology & Evolution 12(4): 143–149. Google Scholar


G Schwarz 1978: Estimating the dimension of a model. - The Annals of Statistics 6(2): 461–464. Google Scholar


M Stüwe& C Grodinsky 1987: Reproductive biology of captive Alpine ibex (Capra i. ibex). - Zoo Biology 6(4): 331–339. Google Scholar


N Sugiura 1978: Further analysis of the data by Akaike's information criterion and the finite correction. - Communications in Statistics: Theory & Methods 7(1): 13–26. Google Scholar


G Tavecchia P Besbeas T Coulson B Morgan& T Clutton-Brock 2009: Estimating Population Size and Hidden Demographic Parameters with State-Space Modeling. - American Naturalist 173(6): 722–733. Google Scholar


S Terzago C Cassardo R Cremonini& S Fratianni 2010: Snow Precipitation and Snow Cover Climatic Variability for the Period 1971-2009 in the Southwestern Italian Alps: The 2008-2009 Snow Season Case Study. - Water 2(4): 773–787. Google Scholar


C Toïgo& J-M Gaillard 2003: Causes of sex-biased adult survival in ungulates: sexual size dimorphism, mating tactic or environment harshness? - Oikos 101(2): 376–384. Google Scholar


C Toïgo J-M Gaillard& J Michallet 1997: Adult survival pattern of the sexually dimorphic Alpine ibex (Capra ibex ibex). - Canadian Journal of Zoology 75(1): 75–79. Google Scholar


C Toïgo J-M Gaillard M Festa-Bianchet E Largo J Michallet& D Maillard 2007: Sex- and age-specific survival of the highly dimorphic Alpine ibex: evidence for a conservative life-history tactic. - Journal of Animal Ecology 76(4): 679–686. Google Scholar


J.C Villaret& R Bon 1995: Social and spatial segregation in alpine ibex (Capra ibex) in Bargy, French Alps. - Ethology 101(4): 291–300. Google Scholar


J Villaret R Bon& A Rivet 1997: Sexual segregation of habitat by the alpine ibex in the French Alps. - Journal of Mammalogy 78(4): 1273–1281. Google Scholar


A von Hardenberg B Bassano V Percacino A Jacobson J von Hardenberg& A Provenzale 2000: Preliminary analysis of the temporal variability of the Alpine ibex population in the Gran Paradiso National Park. - Ibex Journal of Mountain Ecology 5: 201–210. Google Scholar


A von Hardenberg B Bassano& A Provenzale 2009: Temporal variability in juvenile survival explains the strong population decline in the Alpine ibex Capra ibex population in Gran Paradiso National Park (North-Western Italian Alps). - V World Conference on Mountain ungulates. Granada, Spain. Google Scholar


N Yoccoz& J Gaillard 2006: Age structure matters for Alpine ibex population dynamics: comment on Lima & Berryman (2006). - Climate Research 32(2): 139. Google Scholar


[1] Edited by Associate Editor: Stefano Focardi


Appendix I. Parameter values of the best selected models

The following tables contain parameter values of the selected best models, together with their statistics. The labelling scheme of models and parameters is explained in the main text (see as references equations 3 and 4). Since there is no risk of confusing βi,C with βi,w or vice versa, we have omitted the second subscript in the tables. For threshold models, the superscripts refer to parameters calibrated using data of years with snow depth lower (L) or higher (H) than the threshold of S̄ = 154 cm used in Jacobson et al. (2004). The column ‘Best fit’ indicates the best fitted value of the parameter obtained by minimising the square errors (data until 1980), while ‘Unbiased value’ indicates the bias-corrected parameter as suggested in Efron et al. (1993). The column μ ± SD represents the mean ± the standard deviation of the parameter values obtained by using the bootstrapped 1,000 samples (see main text for details). CV is the coefficient of variation.

Appendix 1. Table 1.

Parameters for the logarithm of the total annual growth rate qualifying the unstructured threshold model NT1.


Appendix 1. Table 2.

Parameters for the logarithm of the total annual growth rate qualifying the unstructured threshold model NT2.


Appendix 1. Table 3.

Parameters for the logarithm of the adult female rate of demographic variation of the age-structured model FT1.


Appendix 1. Table 4.

Parameters for the logarithm of the adult female rate of demographic variation of the age-structured model FT2.


Appendix 1. Table 5.

Parameters for the logarithm of the adult male rate of demographic variation of the age-structured model MT.


Appendix 1. Table 6.

Parameters for the logarithm of the kid rate of demographic variation of the age-structured model KC1.


Appendix 1. Table 7.

Parameters for the logarithm of the kid rate of demographic variation of the age-structured model KC2.


Appendix 1. Table 8.

Parameters for the logarithm of the kid rate of demographic variation of the age-structured model WC.


Figure 1.

Available data for Alpine ibex in the GPNP. A) Fall counts of the total population abundance (solid line, left axis) and mean winter snow depth St (dashed line, right axis) estimated at Serrù station as the average from November of year t to May of year t + 1. B) Number of ibex in the four population compartments: adult males (solid line), adult females (dashed line), yearlings (dotted line) and kids (dash-dotted line).


Figure 2.

Models for kid rate of demographic variation (A) and weaning success (B) vs snow depth, where asterisks indicate the data and the dotted vertical line corresponds to the snow threshold value of 154 cm. The continuous line shows the fit from models KC1 (A, the model IDs are as shown in Table 1 ) and WC (B) while dashed lines correspond to optimal threshold models. As all the models included both the snow depth and the interaction between snow depth and density, curves were computed using the mean value of the population density over the entire study period.


Figure 3.

One-step-ahead prediction ranges of the best models (5th to 95th percentiles and mean; grey bars) vs the logarithm of observed vital rates (connected with black solid or grey dashed lines). Parameters, whose values and statistics are in tables of Appendix I, were calibrated using 1,000 bootstrap samples from data along the entire time span. The figure shows the total annual growth rates predicted by the unstructured models NT1 (A) and NT2 (B) and the results of the structured models: C) female rate of demographic variation predicted by model FT1 (left axis, black solid line) and FT2 (right axis, grey dashed lined), D) male rate of demographic variation (model MT), E) kid rate of demographic variation predicted by model KC1 (left axis, black solid line) and KC2 (right axis, grey dashed line) and F) weaning success (WC). The structure of each model is summarised in Table 1.


Figure 4.

Simulation ranges for the dynamics after 1980 of the total number of ibex obtained with unstructured models NT1 (A) and NT2 (B), and with global models STR1 (C), STR2 (D), STR3 (E) and STR4 (F). Parameter values were calibrated using 1,000 bootstrap samples over the first 20 years of data. The dark gray areas include the 25th to 75th percentiles, while light gray areas are for the range 5th to 95th percentiles. The thick continuous line connects data while the thin dashed line corresponds to the reference simulation obtained with the unbiased parameters reported in the tables of Appendix I.


Figure 5.

Temporal evolution of RMSEN(k) as defined in the main text. All parameter values as in the reference simulations of Figure 4.


Figure 6.

Simulations for adult males and adult females obtained with the models in equation 5 using parametrisation as in FT1 (A), FT2 (B) and MT (C), respectively. Simulated kids (D) and yearlings (E and F), respectively, were obtained by applying models WC, KC1 and KC2 in equation 6. Parameter values and curve types are as in Figure 4.


Table 1.

The best models according to selection criteria SRM, BIC and AICc as explained in the main text. The first letter of the Model ID indicates unstructured modelling (N) or structured modelling: females' (F), males' (M) and kids' rate of demographic variation (K) or weaning success (W). The second letter in Model ID is T for threshold formulations or C for continuous formulations. Dt is the compartment considered for density dependence: ‘all’ indicates total population density (i.e. Dt = Nt), ‘females’ indicates Dt = Ft + Kt + Yt/2 while ‘males’ indicates Dt = Mt + Yt/2. p is the number of free model parameters (not including variance), while symbol ‘x’ denotes inclusion of the corresponding polynomial term into the best models. The column R̂2 contains the adjusted R2 considering the degrees of freedom of the model according to the formula R̂2 =1 – (SSE(n – 1)/SST(n – p – 1)), where n is the number of calibration data, SSE is the sum of squared prediction errors and SST is the total sum of squares. All indicators of model performances are calculated over the entire data set. For all selection criteria reported in the relevant column, the corresponding indicator should be minimum.

Andrea Mignatti, Renato Casagrandi, Antonello Provenzale, Achaz von Hardenberg, and Marino Gatto "Sex- and age-structured models for Alpine ibex Capra ibex ibex population dynamics," Wildlife Biology 18(3), 318-332, (1 September 2012).
Received: 6 September 2011; Accepted: 1 April 2012; Published: 1 September 2012

Back to Top