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 (AIC_{c}, 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 (N_{1977}/N_{1976} = 0.64, N_{t} being the autumn population density of year t) and in 1977 (N_{1978}/N_{1977} = 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

where N_{t}is the total population density, S

_{t}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(N

_{t+1}/N

_{t}) = f(N

_{t},S

_{t}) 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 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

_{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 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.

_{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 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 β_{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 (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^{2}.

### 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_{1981} and use the snow depth series S_{1981}, S_{1982}, ..., S_{2004} 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 Y_{t}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 and 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 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 N_{t}:

## Results

### 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_{t}N_{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_{t}S_{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 β_{2} multiplying S_{t} were positive while the coefficients β_{4} 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^{2}, 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} = X_{1981} 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 β_{0}.

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_{t}D_{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.

## Acknowledgements

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.

## References

*Capra ibex*). - International Journal for Parasitology 40(11): 1285–1293. Google Scholar

*Capra ibex*populations? - Wildlife Biology 14(4): 489–499. Google Scholar

*Capra i. ibex*). - Zoo Biology 6(4): 331–339. Google Scholar

*Capra ibex ibex*). - Canadian Journal of Zoology 75(1): 75–79. Google Scholar

*Capra ibex*) in Bargy, French Alps. - Ethology 101(4): 291–300. Google Scholar

*Capra ibex*population in Gran Paradiso National Park (North-Western Italian Alps). - V World Conference on Mountain ungulates. Granada, Spain. Google Scholar

## Appendices

**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.

## Table 1.

The best models according to selection criteria SRM, BIC and AIC_{c} 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. D_{t} is the compartment considered for density dependence: ‘all’ indicates total population density (i.e. D_{t} = N_{t}), ‘females’ indicates D_{t} = F_{t} + K_{t} + Y_{t}/2 while ‘males’ indicates D_{t} = M_{t} + Y_{t}/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 R^{2} 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.

*Capra ibex ibex*population dynamics," Wildlife Biology 18(3), 318-332, (1 September 2012). https://doi.org/10.2981/11-084