Several processes and mechanisms have previously been suggested to explain the cyclic appearance of many grouse populations (Angelstam, Lindström & Widén 1984, Lindström 1996, Hudson, Dobson & Newborn 1998, Watson, Moss & Rae 1998). For instance, Finnish populations of capercaillie Tetrao urogallus, black grouse Tetrao tetrix and hazel grouse Bonasa bonasia display 6–7 year periodicities (Lindén 1989, Lindström, Ranta, Kaitala & Lindén 1995), which may be the result of delayed density dependence (Lindström, Kokko, Ranta & Lindén 1999) and external perturbations (Kaitala, Ranta & Lindström 1996a,b). Spatial correlation of external perturbations could further explain why the three species fluctuate in synchrony both in time and space (Moran 1953, Royama 1992, Ranta, Lindström & Lindén 1995, Ranta, Kaitala & Lindström 1997, Lande, Engen & Sæther 1999, Ranta, Kaitala & Lindström 1999). Stochastic perturbations of breeding success have been suggested to be a likely disturbance factor in Finnish grouse, with the effect of altering the dynamics from damped to quasi-periodic fluctuations (Kaitala et al. 1996a,b) and synchronising the spatial dynamics (Ranta et al. 1995). However, as pointed out by Kaitala et al. (1996b), there is no reason to believe that the perturbation factor should be the same every time, nor does it have to be the same factor for all species within a year.

An alternative, or complementary perturbation factor, is harvesting. Some of the most well-known time series in ecology, such as the fur trade records of the Canada lynx Lynx canadensis, are found in exploited populations (Kendall, Prendergast & Bjømstad 1998). In fact, the time series we analyse are often bag records or catch per effort plotted over time. To study the harvesting process and its dynamic consequences are important if we are to understand to what extent different processes (including harvesting) contribute to the population dynamics in exploited populations (Jonzén, Ripa & Lundberg 2002).

The purpose of this paper is twofold. First, we will demonstrate that harvesting per se may cause a population with damped internal dynamics to show quasicyclic fluctuations if the harvest fraction is time-variant. Second, we will examine the population dynamic effects of two non-linear harvesting strategies; an ‘adaptive’ harvesting strategy, where the harvest fraction is set higher when the population is increasing than when it is decreasing, and a threshold harvesting strategy where a proportion of the population is harvested above a threshold density. The purpose is not to investigate all possible harvesting strategies, but rather to highlight some of the potential population dynamic consequences of harvesting in terms of autocorrelation structure and frequency properties of the population dynamics. We will achieve this by analysing a model shown to represent the population dynamics of black grouse in Finland reasonably well (Lindström et al. 1999), thus extending the linear analysis recently presented by Jonzén et al. (2002).

Methods

Following Kaitala et al. (1996b), we used a non-linear autoregressive model structurally modified from Royama (1992) and added a harvesting term:

where N_t+1 denotes population density at time t+1, a₁ and a₂ give the strength of density-dependence at lag 1 and 2, respectively. This model (without harvesting) has previously been found to perform well, compared to alternative time series models, in describing the temporal dynamics of Finnish grouse populations, including black grouse (Lindström et al. 1999). The model is scaled such that the equilibrium without harvest is 1.

We investigate two different harvesting strategies. First, we let the fraction of the population harvested, H_t be beta distributed, hence, restricted to [0, 1], which must be the case regardless of the distribution of annual kills (Lauck, Clark, Mangel & Munro 1998). The mean and variance of a beta distribution (see e.g. Gelman, Carlin, Stem & Rubin (1995)) are determined by two parameters, α and [β. If α and [β are equal, the probability density function is symmetric around 0.5. Harvesting is taking place after reproduction and all individuals are equally vulnerable. There is no stochasticity except in the harvest term.

Second, we introduce environmental stochasticity according to

where e, is a random normal deviate with mean zero and unit variance. The value of o sets the magnitude of the environmental stochasticity. Without harvesting, this model generates fluctuations with a period of 5–7 years for certain values of a₁ and a₂. Now, instead of using a beta distributed harvest fraction, we let the harvest fraction depend on whether the population size before harvesting has increased or decreased compared to the previous year. This seems to be the pattern of grouse shooting in Finland (Lindén 1991). We call this an ‘adaptive’ harvesting strategy and set H_t = 0.15, 0.25 or 0.35 after an increase, and set H_t = 0.05 after a decline. For comparison, we also evaluate a threshold strategy where a constant fraction is harvested when the density is above 75% of the unexploited equilibrium density.

The effect of a₁ and a₂ has recently been treated in detail (Kaitala et al. 1996b), and we will illustrate our results using a₁ = -0.12 and a₂ = -0.71, estimated for black grouse in the province of Turku-Pori (Lindström 1996). The black grouse in this province demonstrates very obvious cyclic dynamics, and we use the parameter estimates from that population to show how cyclic populations in general can be affected by different harvesting strategies. We ran model 1 for different values of a and (3, varying the mean harvest fraction and the coefficient of variation. In model 2, the value of a was set to 0.1 or 0.3. We simulated the two models for 800 generations. Each parameter combination was repeated 1,000 times, and we analysed the last 500 generations. Periodicity was judged from the autocorrelation function (ACF). When analysing the dynamic effect of the adaptive and threshold harvesting strategies, we also plotted the periodogram (see e.g. Chatfield 1999).

Figure 1.

Population fluctuations of black grouse in the province of Turku-Pori. Finland during 1964–1984 (A) with data detrended and standardised to mean zero and unit variance. The probability density function (Pdf; in panel B) of the beta distributed harvest fraction is used in panel C. The parameter values in the beta distribution in panel B are α = 21.1 and β = 119.6, corresponding to a mean harvest fraction of 0.15 and coefficient of variation of 0.2. Panel C gives the trajectory of standardised population density generated by model 1, when the harvest fraction is beta distributed according to panel B. The other parameter values are a₁ = -0.12 and a₂ = -0.71, estimated from the time series on black grouse in Turku-Pori (Lindström 1996). We ran 1,000 iterations for 800 years and analysed the last 500 years. For illustrative purposes, only the first 40 years of this period are shown in panel C. Panel D gives h autocorrelation function (ACF) for the simulated population trajectory (for the last 500 years) when harvest fraction is beta distributed, and the dotted lines indicate the approximate 95% confidence limit (Chatfield 1999).

Results

Whether model 1 shows periodic fluctuations with significant negative time lags at year 3–4 and positive lags at year 6–7 (indicating a period length of 6–7 years) depends on the value of a₁, a₂ and the coefficient of variation of the harvesting. For a moderate mean harvest fraction, say 0.15, with a coefficient of variation around 20%, model 1 generates fluctuations with a period of 6–7 years, in accordance with the dynamics in Finnish grouse (Fig. 1). However, for a high target harvest fraction combined with very low coefficient of variation, the periodicity corresponding to grouse dynamics can not be found by the autocorrelation function. The same is true if the parameters in the beta distribution (α and β) are both equal to one, hence, the harvest fraction is uniformly distributed with a mean of 0.5. These two cases seem to be very unlikely in grouse hunting or any other exploited population.

The periodicity is, however, dependent on the deterministic part of the model. We may also be interested in the amplitude of the fluctuations, which increases with increasing variation in the harvest fraction. In real grouse populations, however, we will find the dynamics to be affected by not only a mixture of intrinsic dynamics and a variable harvest fraction, but also by environmental noise. If variable harvesting is independent of environmental noise, the two components of the variance are additive. Therefore, variable harvesting and environmental noise will work in concert, potentially reinforcing the amplitude of population fluctuations.

In Figure 2 we show the effect of implementing the ‘adaptive’ harvesting strategy. The overall effect is stabilising, changing the periodicity from 5–6 year cycles to long-term fluctuations without any clear cyclic pattern. This is demonstrated by the autocorrelation function and further substantiated by inspection of the power spectrum. The effect is more pronounced when the harvest fraction is 0.25 (see Fig. 2E–F) compared to 0.15 (see Fig. 2C–D) after an increase in population size. If we increase the harvest fraction after a population increase even more (h = 0.35), the resulting population dynamics are a mixture of short and long-term variation. However, this effect of harvesting disappears if the magnitude of environmental stochasticity is changed from Σ = 0.1 to Σ = 0.3 (see Fig. 2G–H).

Also the threshold harvest strategy has the potential to dampen the cyclic dynamics to long-term fluctuations (Fig. 3), but when a very high fraction of the population is harvested (h = 0.6 above the threshold) and the magnitude of the environmental stochasticity is moderate (Σ = 0.1), the population jumps between low and high density in a regular fashion and high frequency variation dominates the dynamics (see Fig. 3G–H).

Figure 2.

Estimated autocorrelation functions (ACF) and spectral powers expressed in periodograms for model 2, assuming either no harvest (A–B), or an ‘adaptive’ harvest strategy (C–H). In the adaptive harvest strategy, the harvest fraction is low (h = 0.05) after a decrease and higher after an increase in population density: h = 0.15 (C–D), h = 0.25 (E–F), h = 0.35 (G–H). We ran 1.000 iterations for 800 years and analysed the last 500 years. The dotted lines in the left-hand panels (A,C,E,G) indicate the approximate 95% confidence limit (Chatfield 1999). Parameter values are: a₁ = -0.12, a₂ = -0.71 and Σ = 0.1 (solid line) or 0.3 (dashed line).

Figure 3.

Estimated autocorrelation functions (ACF) and spectral powers expressed in periodograms for model 2. assuming either no harvest (A–B), or a threshold harvest strategy (C–H) in which harvest was abandoned when the population density dropped below 75% of the unharvested equilibrium size. Above that threshold, the harvest fraction, h. was 0.2 (C–D), 0.4 (E–F), or 0.6 (G–H). We ran 1,000 iterations for 800 years and analysed the last 500 years. The dotted lines in the left-hand panels (A,C,E,G) indicate the approximate 95% confidence limit (Chatfield 1999). Parameter values are: a₁ = -0.12, a₂ = -0.71 and Σ = 0.1 (solid lines) or 0.3 (dashed lines).

Discussion

Considering the economic importance of many exploited populations, it is important to identify sources of variation and how they interact with internal dynamics to produce the time series we analyse. There is nothing new in suggesting that damped internal dynamics topped with noise can produce quasi-cyclic oscillations (Nisbet & Gurney 1982, Potts, Tapper & Hudson 1984, Stenseth 1999), but harvesting has rarely been treated as a stochastic variable (but see Lauck et al. 1998, Patterson 1999, Mangel 2000, Jonzén et al. 2002), and the potential role of harvesting as an external noise factor keeping periodic fluctuations going has never been treated in detail (e.g. Kendall et al. 1998). It has, however, been noted that if there is a lagged response of harvest fraction to population size, the effect may be periodic fluctuations of the exploited resource (Botsford, Method & Johnston 1983).

Historically, fixed exploitation rate (constant fraction) strategies have often been implemented in population management (Hilbom & Walters 1992, Walters & Parma 1996). This policy has recently been questioned by Lande, Sæther & Engen (1997), who have shown that strategies involving a threshold density below which no harvest takes place are preferable when the risk of extinction is taken into account. We have shown here that one such threshold strategy (where a constant fraction is harvested above the threshold) may affect the autocovariance structure of the population dynamics and, hence, the relative contribution of short and long-term variability. Also the ‘adaptive’ strategy had such effects. Interestingly, we saw some effects of harvesting that were not predicted by the linear analysis by Jonzén et al. (2002). In a linear model, harvesting decreases the magnitude of autocorrelation in the population dynamics. However, the non-linear analysis performed here show that harvesting can even shift the sign of the first-order autocorrelation coefficient (see Fig. 3G) or give rise to a mixture of short and long-term variability (see Fig. 2G–H). Hence, we feel that the topic deserves further interest. Changing the periodicity and/or the relative contribution of high and low frequent variation could of course have management implications if not only the average yield and population size is of concern, but also the variation.

Much of the knowledge about the population dynamics in Finnish grouse species is based on hunting statistics and an intensive monitoring programme (Rajala 1974, Lindén & Rajala 1981). However, there are just a few studies on the impact of hunting on Finnish grouse (e.g. Lindén 1981, Lindén & Raijas 1986, Baines & Lindén 1991, Lindén 1991, Lindström 1994). Recently, Lindström (1994) investigated the possible effects of hunting on the dynamic behaviour of a population. He compared different hunting strategies, using a logistic growth model, and concluded that, generally, the hunting pressure has to be high in order to change the dynamic behaviour. His study differs from ours in that he compared different hunting strategies and analysed a deterministic model, rather than letting the harvest fraction per se be a stochastic variable, as assumed in our first model. An analysis of bags vs density over time (see Figure 1 in Lindén 1981) suggests that the actual fraction taken varies substantially over time, hence, supporting our idea about time-variant harvest fraction. However, the harvest fraction is higher during the increase phase and considerably lower during the decrease phase, suggesting that our ‘adaptive’ harvest fraction is a better description of grouse hunting in Finland. If so, one would not expect hunting to be the perturbation factor maintaining the grouse population cycles going, since the ‘adaptive’ harvest strategy may change the periodicity and even remove the cycles (as seen in Fig. 2). On the other hand, the time series of grouse densities are of course the result of various deterministic and stochastic elements, including harvest. Hence, it is difficult to judge the effect of harvesting in practice, if we do not have a control group which is not exploited. The exact cause of the Finnish grouse cycles therefore remains elusive.

In some exploited populations, harvesting is the most important source of mortality and potentially an equally important source of variation (Jonzén et al. 2002). We therefore feel that it is time for the harvest process to receive more attention from ecologists studying the dynamics of exploited population. Often, harvesting has been considered as nothing but a disturbance that lowers the population equilibrium unless completely compensated for. We have shown here that harvesting has the potential to modulate the temporal dynamics of an exploited population in many ways. Studying the effect of harvesting on resource dynamics may therefore teach us something about population dynamics in general.