Study area

The study area is located in Trento Province in the Central-Eastern Italian Alps (46°02’N, 10°38’E), across three chamois hunting districts: Adamello (area = 373 km2), Presanella (146 km2) and Brenta (263 km2). The area is forested up to the tree-line at about 2,000 m, above which it consists of Alpine meadows, rocky outcrops, scree fields and open rock faces. The average altitude varies among the districts, though with considerable overlap (mean altitude ± SD: Adamello, 1,901 ± 616 m; Presanella, 2,098 ± 540 m; Brenta, 1,594 ± 603 m). Adamello and Presanella are characterised by nutrient-poor siliceous vegetation whilst Brenta is characterised by nutrient-rich calcareous vegetation [37]. Typically, meadows in Adamello and Presanella are dominated by Festuca scabriculmis and Carex curvula, whilst those in Brenta are composed of Sesleria albicans and Carex firma. Throughout the study area, meadows are grazed by small herds of livestock (sheep, goats and cows) during summer, a practice that has been maintained at consistent levels throughout the study period. Several potential predators of chamois were present during the study, including a small, stable population of brown bear (Ursus arctos) in Brenta, a very small number of Eurasian lynx (Lynx lynx) and the golden eagle (Aquila chrysaetos). However, predation on chamois is very rare here (personal communication, Adamello Brenta Nature Park, Trento Province, Italy).

Body mass data

Chamois are hunted every year between mid-September and late-December. Data were collected on the eviscerated body mass and day of shooting of 10,455 yearling (≈1.5 year olds; hereafter juveniles) Alpine chamois (5,762 males and 4,693 females), hunted between 1979 and 2010 (see Additional file 1 for annual breakdowns of sample size). Hunting is heavily regulated and there is little potential for artificial selection by hunters, as chamois can easily detect hunters in the predominantly open habitat and will flee from hunters at particularly large distances [38],[39]. Moreover, there is no evidence of hunters preferentially harvesting larger bodied age-classes in these populations [38]. Hunting pressure on yearlings varies among sites (mean proportion of yearlings hunted in census years: Adamello males, 0.40 ± 0.01; Adamello females, 0.32 ± 0.01; Presanella males, 0.32 ± 0.01; Presanella females, 0.24 ± 0.02; Brenta males, 0.37 ± 0.01; Brenta females, 0.31 ± 0.02). In order to account for intra-seasonal variation in body mass, which is not the focus of this study, a published model of seasonal body mass change [38], which considers inter-annual mass variation, was used to estimate juvenile mass standardised to a specific day of the year. Annual estimates (n = 32) of mean juvenile body mass were produced for each sex, within each site, standardised to day 300 of the year (27th October) (see Figure 1). Body mass was estimated after the vegetation growing season (hereafter `growing season') because body condition at that time will have been influenced by the spring and summer environment, which is thought to have a strong influence on ungulate body mass [40].

Figure 1 Temporal juvenile body mass trends. Long-term temporal trends in body masses of a) male and b) female juvenile chamois in the three study sites between 1979 and 2010. Points are annual mass estimates standardised to day 300 and straight lines are fitted trends. Full size image

There were clear negative temporal body mass trends in all sexes and sites (Figure 1). In order to examine drivers of deviations from the long-term trends (i.e. years in which mean body mass was particularly high or low, even given the trend), the body mass time series were detrended by fitting linear models and calculating residuals. However, detrending can remove long-term fluctuations related to environmental trends [41],[42], which are of primary interest to us. As a result, we modelled body mass data, to examine drivers of long-term trends, and also modelled body mass residuals, to examine drivers of deviations from the trends.

Environmental and demographic data

A range of climatic and non-climatic factors might be expected to influence chamois body mass. Negative effects of population density on mass are common in ungulates [26],[27]. In the absence of natural predation, these effects generally operate through increased intra-specific competition at higher population densities, resulting in lower per-capita food intake, particularly during periods when food is scarce [43]-[45]. To investigate density-dependence in these chamois populations, site-specific population density estimates were used from total population censuses performed in September every year between 1981 and 2009 (with the exception of 1990 and 1991; data from these years were excluded from the analysis). Each year, a set of simultaneous censuses was performed from vantage points across different blocks of each hunting district. It was assumed that density estimates from this time of year would reflect the population density over the previous growing season.

To investigate a possible direct thermoregulatory link between climate and body mass, we calculated yearly site-specific estimates of mean daily growing season temperature between 1982 and 2007 from high-altitude meteorological stations located in each of the three study sites (Data provided by The Forecasts and Organization Office, Civil Protection and Infrastructures Department, Trento Province, Italy). Differences in the elevation of weather stations among sites contributed to inter-site temperature differences (see Figure 2a). However, this did not affect our analysis since the drivers of body mass trends were examined separately in each site and, additionally, temperatures were standardised within sites, along with all other environmental predictors (see `Modelling variation in mass and mass residuals'). The bounds of the growing season were estimated using snow cover data, also from meteorological stations located within each site. The growing season was defined as the period between the snow-melt in spring, when snow cover was reduced to 0% (which generally occurs between late March and early May), and the first significant snowfall in winter that results in new snow settling on the ground (which generally occurs between early November and late December).

Figure 2 Temporal trends in population density and mean growing season temperature. Long-term variation in a) mean growing season daily maximum temperature, between 1982 and 2007, and b) population density, between 1981 and 2009, in Adamello (black), Presanella (red) and Brenta (green). Gaps show years with missing data. Whilst the three study areas do differ in their climate, some of the observed inter-site differences in temperature in a) are due to variation in the elevation of weather stations among areas. Full size image

To investigate the effect of vegetation productivity and phenology on mass, NDVI (normalised difference vegetation index) data were used as a measure of vegetation productivity, processed by the Global Inventory Modelling and Mapping Studies group (GIMMS; [46],[47]). These data are global at a 0.07 degree resolution (approximately 8 km by 8 km) and are available at fortnightly intervals between 1982 and 2006 (thus slightly restricting the yearly data range for analyses). In order to focus on vegetation types utilised by chamois for foraging, such as alpine meadows and sparsely vegetated areas, only NDVI pixels dominated by such vegetation types were considered. To do this the Corine land-cover 2006 data-set at a 100 m resolution [48] was used to select only NDVI pixels within each site containing less than 25% coniferous woodland. As each of the three sites encompassed a number of these NDVI pixels, mean NDVI from these pixels was calculated for each fortnightly time period, within each site. Previous studies have implicated a number of metrics relating to annual NDVI variation as being important to ungulate body condition (e.g. [24],[25]). Here, we seek to derive NDVI metrics in a standardised fashion, despite inherent noise in NDVI estimates caused by factors such as cloud cover, water, snow or shadow [47]. As in previous studies [49], we used a smoothed function to characterise variation in NDVI with time in a given year. The following function was used (see Additional file 2, for an illustration of the functional form):

p - s , y , t = α s , y + β s , y - α s , y exp - t - t * σ s , y z s , y .

Here, p - s , y , t is predicted NDVI at time-period t in site s and year y, α s,y and β s,y are minimum and maximum NDVI respectively in site s and year y, σ s,y is a parameter related to the width of the function and z s,y is a parameter describing the shape of the function. Variation in NDVI data, p(s,y,t), about the predicted mean was beta distributed. Thus, the likelihood of the model parameters, θ s , y , given the data, parameterised by p - s , y , t and the dispersion coefficient ϕ s,y , is

L θ s , y = Π t Γ a s , y , t + b s , y , t Γ a s , y , t Γ b s , y , t p s , y , t a s , y , t − 1 1 − p s , y , t b s , y , t − 1 ,

where Γ(x) is the gamma function, a s , y , t = p - s , y , t / ϕ s , y and b s , y , t = 1 - p - s , y , t / ϕ s , y . The most parsimonious fit was identified using Akaike’s Information Criterion (AIC) [50],[51].

The most parsimonious fitted relationships between mean NDVI and time in each year and site were calculated. Using these relationships, four NDVI metrics, described below, were calculated relating to vegetation productivity and phenology. All four of the metrics selected have been highlighted as important either to juvenile chamois specifically or to other ungulate species [47],[52]. Previously, climate-induced changes in spring growing conditions [53] have been linked to higher juvenile body mass in ungulates, including chamois [52], due to longer growing-seasons and higher vegetation quality [24]. However, warmer springs have also been linked to negative impacts on body mass as higher temperatures lead to faster rates of vegetation `green-up' and, thus, a shorter period of access to nutrient-rich emergent vegetation associated with early spring [25],[54]. Here, the following four metrics were used: maximum rate of spring green-up, growing season duration, maximum NDVI and total growing season NDVI. Maximum rate of spring green-up (NDVI rate ) was calculated as the maximum first derivative of p - s , y t (i.e. the maximum rate of NDVI increase). The duration of the growing season (NDVI dur ) was calculated as the length of time between the maximum second derivative of p - s , y t (the start date of the growing season; when the rate of NDVI increase is increasing at its maximum rate) and the minimum second derivative of p - s , y t (the end date of the growing season; when NDVI is decreasing most rapidly). Maximum NDVI (NDVI max ) was calculated as the maximum value of p - s , y t and total growing-season NDVI (INDVI) as the integral of p - s , y t within the bounds of the growing-season. An illustration of the calculation of these metrics can be seen in Additional file 2.

Modelling variation in mass and mass residuals