|Home | About | Journals | Submit | Contact Us | Français|
A problem in climate studies has been on how to treat causal chains of explanations and both direct and indirect effects. Mammals in strongly seasonal environments of the boreal forest typically lose condition during winter and gain mass (and reproduce) during the summer season when biomass and plant quality peak. Mass decay of large herbivores during winter is due to direct effects of winter weather, such as increased costs of movement, thermoregulation and reduced access to food when snow is deep. Deer condition during summer is thought to be affected mainly indirectly by weather through plants. High spring temperature speeds up plant development, and deep snow can delay phenology in early summer. Current statistical modelling does not take into account these mechanistic pathways. We used hierarchical Bayes modelling to more mechanistically link global climate, local weather and plant phenology to autumn body mass of red deer in Norway. Red deer were much more affected indirectly through trophic interactions. No solid evidence of direct effects of snow depth was found on autumn body mass. We discuss the implications of our results relative to our ability to predict effects of global change on large mammalian herbivores in the boreal forest.
There is increasing consensus among ecologists that climate change is affecting ecosystems (Stenseth et al. 2002; Walther et al. 2002, 2005; Parmesan 2006). Climate change and variability are not only affecting many species directly, but also through trophic interactions among species (Schmitz et al. 2003; Krebs & Berteaux 2006). There is a huge variety in actual mechanisms by which animals can be affected by climatic variability. Even the same species may be limited by different weather variables depending on year or area. Increased winter mortality in the Soay sheep (Ovis aries) population on St Kilda, Scotland was due to years with strong winds, heavy rainfall or extreme cold when the population density was high (Hallett et al. 2004). While the red deer (Cervus elaphus) population on Rum, Scotland was most strongly affected by heavy winter rainfall (Post & Stenseth 1999), the red deer further north along the west coast of Norway were more strongly influenced by snow depth (Mysterud et al. 2000, 2001b; Pettorelli et al. 2005a). What unites all of these examples is the fact that they can be related to variation in a single variable; the state of the North Atlantic Oscillation (NAO) during winter. However, the sign of the correlation of individual performance and the NAO may be opposite for the same species in different regions (Mysterud et al. 2003), suggesting local responses are complex.
There is currently somewhat of a gap between detailed studies of behavioural responses to local weather patterns and the literature linking global climate systems such as the NAO to animal performance. It has been argued that if you are mainly interested in good predictions, you should use the NAO and other coarse-scale indices (Stenseth et al. 2003). If you are interested in learning about mechanisms, you should use local weather variables (Krebs & Berteaux 2006). We believe most ecologists are interested in both, and, clearly, to predict well it would be good to understand the mechanism. A well-studied ungulate system is the red deer along the west coast of Norway, having revealed both direct (Loison et al. 1999) and indirect (Mysterud et al. 2001b; Pettorelli et al. 2005a) effects of winter climate and direct effect of spring climate on population ecology and life history (figure 1). During winter, snow depth affects body mass loss (at least of yearlings, Loison & Langvatn (1998)), possibly due to extra energetic costs of movement in snow (Parker et al. 1984), thermoregulation (Parker & Robbins 1985) and reduced access to forage (Mysterud et al. 1997). However, since snow accumulation during winter affects early plant growth, plant phenology during spring is delayed at high elevation, leading to a more variable phenology after winters with a high NAO index and thus positively affecting body mass (Pettorelli et al. 2005a). Also, a warm April linked to a high phase of the NAO speeds up plant phenology, benefiting deer (Pettorelli et al. 2005a). Based on all of these studies, one can build a tentative mechanistic model on how the NAO probably impacts individual performance for this particular system (figure 1).
The main approach to investigate direct versus indirect effects is based on path analysis (Shipley 2002). Such methods are similar to linear models in the sense that they assume linear relationships among predictor and response variables, but they differ in the sense that they can assume that specific relationships do not exist: for example, one can assume that the NAO has an effect on red deer body mass only through its effect on snow depth and plant phenology. In our case, classical path models cannot be used because we have two levels of variability in the data: individual red deer body mass which depends on individual characteristics such as age or date of harvesting, as well as repetitive variables such as plant phenology, red deer density and climatic variables that are measured only on an annual basis. Our situation could therefore be seen as an example of hierarchical models (Draper 1995), but with specific restrictions concerning linear relationships among the variables. Such approaches are now implemented within the framework of hierarchical Bayes modelling (e.g. Clark & Gelfand 2006). Since data are recorded at two scales, these approaches are therefore suitable to more mechanistically link global climate (NAO), local weather (temperature and snow depth) and plant phenology (Normalized Difference Vegetation Index, NDVI) to autumn body mass of red deer (C. elaphus) along the west coast of Norway.
Red deer are at their northern limit in Norway, mainly living in the hemiboreal (or boreonemoral) region (Abrahamsen et al. 1977). Exceptions are a small area around the Hardangerfjorden, Hordaland county, which is in the nemoral zone, and some areas not far from the Trondheimsfjorden which is in the southern boreal zone. In Kvinnherad municipality (figure 2), Hordaland county, forests are dominated by deciduous and pine forest. However, there has been a large-scale programme for planting Norway spruce (Picea abies L.), which is an important winter habitat for red deer (Ahlén 1965). Snillfjord and Orkdal municipalities (figure 2), Sør-Trøndelag county, are dominated by birch (Betula spp. L.) and pine (Pinus sylvestris L.).
The data derive from municipalities Kvinnherad, Hordaland county in the south and the adjacent municipalities Snillfjord and Orkdal, Sør-Trøndelag county in the north of the red deer distribution range along the southwest coast of Norway. The data used are a subset of the much larger dataset from the whole of the southwest coast (e.g. Mysterud et al. 2001b; Pettorelli et al. 2005a). The reason for selecting these areas is that data from most years between 1982 and 2004 are available (when data on NDVI have been gathered, see below), and as long a time frame as possible is important for considering climate effects.
Data on red deer derive from hunted animals. The hunting period lasts from 10 September to 15 November. Hunting is controlled through licenses issued by local wildlife boards in each municipality. For each deer, the sex, date, location (municipality) and dressed body mass, together with the mandibles, were provided by the hunters. Recorded body mass is dressed mass (i.e. live mass minus head, skin, viscera, bleedable blood and metapodials). Dressed mass in hinds is approximately 58% of live mass (Langvatn 1977) and highly correlated (r=0.994) with total mass in Cervids (Wallin et al. 1996). From the mandibles, we could age calves and yearlings based on tooth eruption patterns (Loe et al. 2004), whereas older animals were aged using annuli in the cementum of the first incisor (Hamlin et al. 2000).
It is well known that harvesting data may give bias in some cases (e.g. Martínez et al. 2005; Mysterud et al. 2006). However, for this time period in Norway, it is unlikely that harvesting bias is strong. Certainly, quotas (given as calf, yearling and adults of each sex) ensure that harvesting is not random. However, within these categories, and based on our detailed experience with the system, there is not likely to be strong selection for phenotypic traits such as body mass. Part of this is due to the hunting tradition, which is for meat and not for trophies (Milner et al. 2006). Part of this is also due to the large quotas and limited time for hunting, which for moose (Alces alces) populations have been shown to reduce selectivity (Solberg et al. 2000). It is simply a high cost in terms of lost chances to be selective in the forested habitat along the west coast of Norway, where group size of red deer is small (average of 1.5 female per male in mixed groups; Bonenfant et al. 2004). The only probable bias we may have (though not documented) is that females with offspring are less likely to be shot than females without (as reported for Norwegian moose populations; Nilsen & Solberg 2006). However, since most females from 3 to 13 years of age ovulate and are likely to calve (Langvatn et al. 2004), this is probably not an important factor. Further, any small bias this may yield has probably been constant over the time period, so that it is unlikely to interfere with the main factor of interest here, i.e. the climate effects.
As a proxy for vegetative phenology (greenness), we used the NDVI which has been used extensively as a proxy of vegetation dynamics (review in Pettorelli et al. 2005b). Data are collected by the National Oceanic and Atmospheric Administration satellites (http://eosdata.gsfc.nasa.gov/). From these, NDVI values (ranging from −1 to 1) have been produced from visible and near-infrared reflectance measurements (NDVI=(NIR−RED)/(NIR+RED), where NIR is the near infrared light reflected by the vegetation and RED is the red visible light reflected by the vegetation). The spatial scale of resolution (pixel size) is 64km2 and an NDVI value is available on a 10-day basis, from 13 July 1981 to present. Based on an earlier study (Pettorelli et al. 2005a), we chose to use NDVI on 1 May as a measure of the greenness in spring, which is the average date of vegetation start in the study area along the west coast of Norway. There are 25 pixels for Kvinnherad, 11 pixels for Orkdal and 10 pixels for Snillfjord.
The NDVI data can be processed in different ways to remove errors (review in Pettorelli et al. 2005b). In a previous analysis (Pettorelli et al. 2005a), we used the so-called PAL data that are no longer updated, while we here use the GIMMS data which are considered the best time series of the NDVI currently available. At the scale of each specific municipality, the mean and standard deviation of NDVI–GIMMS on 1 May (per year for all pixels) are highly correlated (Kvinnherad: r=0.806 [r=0.939]; Orkdal: r=0.813 [r=0.901]; Snillfjord: r=−0.093 [r=0.686]—values in brackets are correlations obtained after removing the extreme outlying value for 1987). As this is not the case with the PAL data (figure 3), it is probable that this is related to the different error processing. It is likely, but remains to be determined, whether this is related to the use of a spatial smoother (krieging). We therefore only use the mean of NDVI in further modelling.
We used the principal component analysis (PCA)-based winter (December–March; hereafter termed NAOw) and April index of the NAO (Hurrell 1995; Hurrell et al. 2003) as our global-scale climate index. The NAOw has a strong impact on the population ecology of red deer along the southwest coast of Norway (Post et al. 1997; Mysterud et al. 2000, 2001b, 2005; Pettorelli et al. 2005a). In previous studies we have shown that there is an interaction between global climate (NAOw) and topography that determines local weather conditions (Mysterud et al. 2000; Pettorelli et al. 2005a). We therefore know that a high NAOw index is correlated with less snow at low elevation, which corresponds to the main wintering areas of red deer (Albon & Langvatn 1992). We therefore use snow depth at low elevation when it comes to the possible direct effect of winter weather on red deer body mass. In contrast, a high NAOw index is correlated with more snow at high elevation, which extends into the summer ranges of red deer (some deer do remain at low elevation during summer as well). We therefore used data deriving from high elevation to measure any indirect effect of snow depth working through plant phenological development. The NAO in April is positively correlated with temperature and an early start of the growing season (Pettorelli et al. 2005a).
There are a large number of meteorological stations along the west coast of Norway, but not all municipalities have weather stations and not all of them cover the years 1982–2004. A further problem is that few of them measure temperature, while many measure precipitation owing to the importance of hydropower in Norway. It is also experienced that temperature is more spatially correlated than precipitation (including snow depth). We therefore searched for stations within or as close to the municipality of interest as possible. However, we also paid particular attention to the elevation and distance from the coast, so as to find as representative weather stations as possible for both high and low elevation for each municipality.
For Kvinnherad, we used the weather station at Eidfjord-Bu (no. 49580) at 165m a.s.l. for temperature in April, Sauda (no. 46610) at 5m a.s.l. for snow depth (average of max) at low elevation during December–March, while we used the station at Liset (no. 49750; Eidfjord) at 748m a.s.l. for snow depth at high elevation. For Snillfjord and Orkdal, we used stations at Tingvoll-Hanem (no. 64550; Møre og Romsdal) at 69m a.s.l. for temperature in April, Sandstad (no. 65450 Hitra; Sør-Trøndelag) at 20m a.s.l. for snow depth at low elevation and Hafsås, Sunndal and Møre og Romsdal (no. 63530) at 698m a.s.l. for snow depth at high elevation.
As data for each analysis below derive from a single or two adjacent municipalities, we did not need to enter the spatial covariates we frequently use to control for imbalance in the sampling design (e.g. Mysterud et al. 2001a). For Snillfjord and Orkdal, we entered municipality as a factor. We also entered an index for population size, which has affected body mass more than climate over this time period (Mysterud et al. 2001c). When analysing data from a single municipality, we only entered the harvest size. When analysing the two adjacent municipalities Snillfjord and Orkdal, we used instead the harvest size relative to the amount of red deer habitat (as in previous analysis; Mysterud et al. 2001c). Over the study period, the number of harvested red deer increased from 314 in 1982 to 805 in 2004 in Kvinnherad, from 57 to 280 in Orkdal and from 100 in 1982 to 323 in 2004 in Snillfjord. Based on reconstruction of cohorts, it is safe to assume that the annual harvest is fairly well correlated with the actual population size (Mysterud et al. 2007). We also entered date of harvest to control for possible growth or decay during the hunting season.
We used hierarchical Bayes modelling with sampling-based methods for fitting and analysis (see Clark & Gelfand (2006) for an introduction). Models included a specification of relationships among environmental variables (climate and NDVI) on one hand (with years as units) and individual body mass and other covariates (date of harvest and density index) on the other (with deer as units). Linear models were used throughout (except splines for modelling nonlinear relationship between body mass and age; see below). These linear models were therefore defined at two levels of variability: years for phenology and climatic variables, and individual deer for body mass. Yearly values for variables such as phenology that are used to predict body mass are assumed to be functions of other climatic variables such as spring temperature or snow depth, but uncertainties in the relationships is taken into account in such models (something one cannot do using, for example, residuals, such as body mass standardized by age and shooting date). graphical models with two levels of variation can be fitted in a Bayesian framework. Bayesian methods assume prior distributions for model parameters (regression coefficients and residual variances) and use Bayes' theorem to derive the posterior distributions of parameters. Numerical methods based on Markov Chain Monte Carlo simulations are used to derive posterior distributions. Model fitting and model diagnostics (in particular, the convergence of numerical simulations) were done using the library BRugs and coda in the R environment (R Development Core Team 2006; Thomas et al. 2006), using R v. 2.4.0. Prior distributions are usually chosen so as to be uninformative, for example, a normal distribution with mean 0 and a very large variance for a regression coefficient. The following prior distributions were used: normal for regression intercepts and slopes with a mean of 0 and a variance of 106, and Gamma for the inverse of residual variances with shape and scale parameters 0.001 and 0.001, respectively. An initial burn-in of 20000 iterations was used, and posterior distributions of parameters were based on 20000 more iterations. We used three chains to check for the stability of posterior distribution estimates (as recommended, for example, in Clark & Gelfand (2006)), and assessed convergence and autocorrelation using the diagnostics tools available in BRugs and coda (e.g. convergence test based on Gelman & Rubin (1992)). We fitted different models with varying levels of complexity regarding the relationships among climatic variables and phenology. Model selection was partly based on our a priori expectations regarding direct and indirect effects, credible intervals and the deviance information criterion (DIC), which trades off the fit of the model as measured by the posterior expected deviance and the complexity of the model as measured by the effective number of parameters (PD). DIC is a criterion similar to AIC, in the sense of finding a compromise between model fit and model complexity . As sample size is hard to define in models with two levels of variability (e.g. Vaida & Blanchard 2005), DIC is based on an estimate of effective sample size directly based on data and model used.
We only present the final results, focusing on the median of posterior distributions and 95% credible intervals (i.e. intervals that cover 95% of the posterior distribution) to interpret effect sizes and their associated uncertainty. More complex models gave similar estimates for regression coefficients, but did always include parameter estimates with credible intervals including 0.
Linearity and residuals were assessed initially using additive models (AM; Hastie & Tibshirani 1990; Wood 2006). Based on previous work on the system, we have a good idea on what to expect. We focused on females only, as males lose much body mass during the rutting season, making it much more complex to model annual mass (Yoccoz et al. 2002). As age explains most of body mass variation, it is important to model this correctly (Mysterud et al. 2001c; Yoccoz & Gaillard 2006). We ended up pooling all individuals above 12 years of age (few individuals are older than 12 years of age and estimates were therefore not reliable). Based on the AM plotting, we then fitted a third-order polynomial from 0 to 3 years of age, and thereafter a linear segment. This is thus a kind of threshold model for age but using the same first-order derivative for 3 years of age as for a regular spline. At the scale of municipality, the mean and standard deviation of the NDVI was highly positively correlated (see above and figure 3), and we therefore retained the mean of NDVI for further modelling. All variables except age were standardized (mean=0, s.d.=1) to make estimates more stable and easier to compare.
An overview of predicted relations as well as estimated effects is presented in figure 1 (all parameters in table 1). The behaviour of numerical chains was found to be satisfactory (e.g. all scale reduction factors were estimated as 1; Gelman & Rubin 1992) and the sample size corrected for autocorrelation used to calculate posterior distributions and credible intervals was at least 15000 for all parameters. All posterior distributions were approximately normal (figure 4). As expected from previous work (Mysterud et al. 2000), the NAO during winter had a negative effect on snow depth at low elevation (winter range), while the effect is opposite at high elevation (summer range; figure 1; table 1). This effect was slightly more marked in the south (Kvinnherad) than in the north (Snillfjord/Orkdal), probably due to the fact that with increasing latitude the overall winter temperature is colder, making the difference in snow fall due to elevation smaller (snowing further down for a given NAOw value). Also as expected, the NAO in April has a positive effect on temperature in April (Pettorelli et al. 2005a); however, the 95% credible intervals included 0 and the evidence was therefore not strong. Temperature in April in turn strongly influenced the NDVI on 1 May. The effect of snow depth on the NDVI on 1 May (associated with increased variability in NDVI) was weakly negative, but the evidence for the effect was weak. There was no evidence for a direct effect of the NAOw on NDVI, after controlling for the weak effect of snow depth on NDVI.
There was a fairly strong effect of the NDVI on body mass, especially in the south (the evidence for such an effect being weaker in the north), while we found no clear evidence of a direct effect of winter snow depth on red deer (figure 1). Mass decay of deer during winter in these coastal areas was thus much less important for autumn mass than the increase in deer condition during summer which is mainly affected indirectly by climate through plants (i.e. trophic interactions). As known, the effects of age and population density on body mass of red deer were marked (figure 4).
Boreal forests are situated at northern latitudes and there is naturally a strong seasonality in the way climate affects deer populations. Our analysis suggests that summer fattening linked to plant phenology development during summer is a more important climatic factor for body condition in autumn than winter decay related to snow depth. Though body mass development of yearlings during winter is linked to winter harshness (snow and temperature pooled in PCA; Loison & Langvatn 1998), it is evident from our analysis that the magnitude of these effects are much smaller, being overshadowed by the indirect effect operating through plants. The winter NAO markedly affects pattern of snow accumulation (Mysterud et al. 2000) and thus increases variation in plant phenology as measured by the NDVI at broad population scales (Pettorelli et al. 2005a). It has been shown that the NAO operates through increased soil moisture, which affects plant growth during summer (Kettlewell et al. 2006). Our study lends some further support to such indirect effects, but evidence pointed mainly to an effect of spring climatic conditions on plant phenology (see also Stewart et al. 2005). This does not necessarily imply that summer is more important than winter for population dynamics. The ‘Klein hypothesis’ (1965) states that while summer season determines the quality of the animals (body mass), winter is the more critical period for limiting deer population sizes. Summer foraging conditions are decisive for condition during rut and thus most important for breeding (Langvatn et al. 2004). Body mass in autumn also has an effect on winter survival of calves (Loison & Langvatn 1998).
Current work on climate effects on mammals has been criticized: ‘correlations abound but mechanistic understanding is limited owing to long causal chains and indirect effects’ (Krebs & Berteaux 2006). Path analysis is intuitively appealing to biologists, as it offers a good solution to issues related to causal chains of explanations. It is indeed very common that there are clear mechanistic hypotheses for how the environmental factors of interest are related to each other (in our case reviewed in figure 1). There are a number of reported mechanistic ways by which weather is known to affect different deer populations in northern areas in general (e.g. Weladji et al. 2002), and explicitly being able to quantify the relative role of each of these factors for a specific population is exciting. However, a limitation with ordinary path analysis is that it requires data on a similar scale. Data on seasonal climate indices typically have one value per year, while data on animals are typically numerous. Using a yearly average is problematic because other covariates such as age vary among individuals and years. Moreover, this will greatly decrease the precision of estimated effects. While it is clearly important to have long time series in order to estimate precisely temporal changes, methods that explicitly take into account different levels of variation as well as direct and indirect effects represent an important development (Clark & Gelfand 2006). Our use of hierarchical Bayes modelling is novel in this respect, enabling us to fully use the potential of the data at a subpopulation scale.
Some of the relationships we have reported in previous analyses at broader scales (Mysterud et al. 2000, 2001b; Pettorelli et al. 2005a) are nevertheless not significant or weaker here (figure 1). The relationships are much stronger in the southern (Kvinnherad) compared with the northern area (Orkdal/Snillfjord); this result was expected due to the stronger impact of the NAO further south (figure 1). The relationship between the NAO and temperature in April was as predicted, though not significant, probably due to inclusion of data from only a single weather station. A more complicated issue is the one between the NAOw, snow depth and NDVI. More snow is predicted to delay the NDVI, but since increases in the NAOw give more snow at high elevation and less snow at low elevation, we expect mainly increased variance in NDVI with increasing NAOw. However, mean and variance of the NDVI were highly correlated and not possible to separate. It remains to be determined with certainty whether this is related to the spatial smoothing of the NDVI–GIMMS data; in any case, it limits the value of this part of the analysis. Also, the scale of the NDVI is 8×8km; thus, the NDVI is measured over many elevations, making it difficult to capture year-related variation in phenology at the relevant scale within a municipality. In the previous analysis (Pettorelli et al. 2005a) we found that the effect of the NAOw interacted with the proportion of high-elevation habitat within a municipality (i.e. a very coarse scale), thus increasing variance in NDVI at high NAO. Clearly, more fine-scale data on the NDVI would be rewarding to determine with more precision how plant phenology is linked to snow accumulation at the scale of deer movements. Newer satellites such as MODIS might provide insights in that respect.
Forecasting and predicting effects of global climate change on mammal populations are challenging (Berteaux et al. 2006). To progress, we need statistical tools to link the various datasets in better ways so as to solve frequent problems of data that are measured at different scales. Hierarchical path analysis as outlined here solves some of these issues. We can now link the NAO to environmental variables in the same model as the deer data, even though datasets are not replicated at the same scale. Whether such approaches will enable better predictions remains to be seen. Integrating possible interactions with age and/or density is also important in future studies (Coulson et al. 2001). There is somewhat of a limitation that the correlation between the NAO and local weather variables may change between periods (Mysterud et al. 2003). Indeed, the strength of correlation between the NAO and performance of animals and birds therefore similarly varies over time (Durant et al. 2004), and a stronger focus on regime shifts in ecology linked to climate may be necessary (Hörnfeldt 2004; Hörnfeldt et al. 2005; Smith et al. 2006).
We acknowledge the financial support of the Research Council of Norway to A.M. (YFF). We are grateful to John Fryxell, Anne Loison, Kiyoko Miyanishi and one anonymous referee for their very useful comments on a previous draft.
One contribution of 12 to a Theme Issue ‘The boreal forest and global change’.