Royal Society Publishing

Can parasites drive population cycles in mountain hares?

Sunny E Townsend, Scott Newey, Simon J Thirgood, Louise Matthews, Daniel T Haydon

Abstract

Understanding the drivers of population fluctuations is a central goal of ecology. Although well-established theory suggests that parasites can drive cyclic population fluctuations in their hosts, field evidence is lacking. Theory predicts that a parasite that loosely aggregates in the host population and has stronger impact on host fecundity than survival should induce cycling. The helminth Trichostrongylus retortaeformis in the UK's only native lagomorph, the mountain hare, has exactly these properties, and the hares exhibit strong population fluctuations. Here we use a host–parasite model parametrized using the available empirical data to test this superficial concordance between theory and observation. In fact, through an innovative combination of sensitivity and stability analyses, we show that hare population cycles do not seem to be driven by the parasite. Potential limitations in our parametrization and model formulation, together with the possible secondary roles for parasites in determining hare demography, are discussed. Improving our knowledge of leveret biology and the quantification of harvesting emerge as future research priorities. With the growing concern over the present management of mountain hares for disease control in Scotland, understanding their population drivers is an important prerequisite for the effective management of this species.

Keywords:

1. Introduction

Understanding what drives population cycles is a central goal of ecology; yet despite more than 75 years of debate, there is no clear consensus on their causation (Turchin 2003). There is however a growing view that trophic interactions play an important role (Berryman 2002). While predator–prey and herbivore–plant systems have been well studied, the role of parasites has received less attention. Despite a strong theoretical basis that parasites can drive host cyclic dynamics (Anderson and May 1978), empirical support is limited. While parasite-mediated effects are thought to contribute to unstable dynamics in Soay sheep Ovis aries (Gulland 1992) and Svalbard reindeer Rangifer tarandus (Albon et al. 2002), empirical evidence that parasites can drive cyclic dynamics in their wild host is presently limited to the red grouse–Trichostrongylus tenuis system (Hudson et al. 1998, but see Lambin et al. 1999).

The mountain hare is the only lagomorph species native to the UK with Scotland containing 99 per cent of the UK population (McGradySteed et al. 1997). Mountain hares are believed to be under threat from habitat loss and fragmentation, local overexploitation, hybridization and competition with the introduced brown hare and a growing concern over large-scale culls of mountain hares to control ticks and louping-ill (McGradySteed et al. 1997; Macdonald et al. 1998; Battersby 2005; Kinrade et al. 2008). Mountain hares are listed in Annex V of the EC Habitats Directive (1992) requiring the UK to ensure their conservation status and sustainable management. In response to growing concerns over the long-term conservation status and present management of the species in 2007, the mountain hare was made a UK biodiversity action plan species. The factors causing fluctuations and long-term changes in the numbers and distribution of mountain hares remain unknown and complicate attempts to inform management through analysis of patterns in abundance. A greater understanding of the species population dynamics is essential for their sound management.

Scottish populations of mountain hares on grouse moorland are characterized by large amplitude fluctuations of variable regularity with a mean periodicity of 9.2 years (Newey et al. 2007b). The reasons for cyclic dynamics remain unclear (Newey et al. 2007a). Mountain hares are non-territorial and social interactions are not thought to be important (Flux 1970; Hewson 1976), and there is no evidence of food limitation (Keith 1983). Mammalian and avian predators are controlled on moorland managed for red grouse in Scotland and therefore, unlike the situation in Scandinavia, predators are not thought to be important in driving mountain hare populations (Newey et al. 2007a). Hares are, however, susceptible to parasite infections, in particular the helminth Trichostrongylus retortaeformis and recent field studies have demonstrated that T. retortaeformis is loosely aggregated in the mountain hare population and that parasite-mediated effects on survival are small compared with parasite-induced reductions in host fecundity (Newey & Thirgood 2004; Newey et al. 2005). These features of the mountain hare–T. retortaeformis system are consistent with the characteristics that analytical host–parasite models suggest can lead to instability and population cycles (May & Anderson 1978).

Major advances in understanding the causes of population dynamics have come from synthesizing modelling and empirical work (Kendall et al. 1999; Turchin 2003). Here we combine empirical field experiments, time-series analysis and modelling to assess whether parasites can drive mountain hare population cycles. A range of observed dynamical patterns have been quantified from time-series analysis, cross-sectional studies and field experiments to generate a list of characteristic properties with which modelled population dynamics can be compared. Hare population densities fluctuate from 20 to 200 hares km−2 (Watson et al. 1973; Hewson 1976), with a range of periods between 4 and 15 years (Newey et al. 2007b). T. retortaeformis burdens average approximately 2000 worms per individual (Newey et al. 2005). Our approach was to contrast these listed properties with equivalent characteristics in modelled mountain hare populations in order to: (i) test whether our present empirical understanding supports parasite driven hare dynamics, (ii) in the case that it does not, identify plausible parameter changes that would lead to population dynamics with the observed properties, (iii) determine whether small changes in parameters can account for the wide diversity of observed dynamics across Scottish populations, and (iv) improve our understanding of the system and prioritize future empirical research activities.

2. Material and methods

We used a variant of the classic Anderson & May macroparasite model (Anderson & May 1978; May & Anderson 1978) introduced by Diekmann & Kretzschmar (1991), which describes continuous growth equations for a host population of density, H, which interacts with a parasite population, PEmbedded Image(2.1)Embedded Image(2.2)Parameters are defined in table 1. The structure of the model encapsulates important elements of the system that include: (i) the negative binomial distribution of parasites among hosts (Newey et al. 2005) described by the mean parasite load P/H and aggregation parameter k, (ii) a transmission rate dependent on host density, and (iii) host fecundity modelled through the use of a multiplicative term to avoid biologically meaningless negative host birth rates (Diekmann & Kretzschmar 1991).

View this table:
Table 1

Empirically determined parameter point estimates and plausible limits. Values in italics generate dynamics closest to those observed (modified parameter set).

Point estimates and plausibility envelopes for parametrizing equations (2.1) and (2.2) are given in table 1. The data sources and methods of estimation are described in the electronic supplementary material, appendix S1. Rather than strict confidence envelopes, plausible ranges of parameters were most practically based on the best available empirical information.

The dynamical properties of the parametrized model were derived using standard analytical techniques and numerical simulations (refer to the electronic supplementary material, appendix S2 for further details). Elasticity analyses were performed to compare the proportional effects of changing each parameter in table 1 on dynamical properties of the model populations.

3. Results

(a) Model parametrization within empirically defined plausibility envelope

Parametrizing the model with the point estimates presented in table 1 resulted in rapidly damped oscillations to a stable equilibrium point where parasite burdens were far greater than those found in mountain hare populations (figure 1a). Elasticity analysis identified that: an increase in hare intrinsic mortality (b) or parasite-reduced hare fecundity (δ) or a decrease in hare intrinsic fecundity (a) would bring about a simultaneous reduction in both stability and parasite burdens. Increasing parasite-induced hare mortality (α) reduced parasite loads but was stabilizing while the parasite parameters (fecundity (λ), adult mortality (μ) and transmission inefficiency (H0)) had little effect on equilibrium parasite load or stability. A new modified parameter set was identified by increasing the values of hare intrinsic mortality (b) and parasite-reduced hare fecundity (δ) and decreasing hare intrinsic fecundity (a) to empirically plausible limits (table 1). The simulated population dynamics maintained a weakly stable equilibrium hare density characterized by weakly damped oscillations with a period within the observed range (figure 1b). However, these changes could not bring parasite loads down sufficiently to be consistent with those found in mountain hares.

Figure 1

Simulated population dynamics based on the empirical information available on the hare–parasite interaction. (ac) Hare density, (df) worms per hare. The model was parametrized with (a) point estimates and (b) the modified parameter set chosen to be the best-fitting combination within the identified plausibility envelope. (c) Realistic population dynamics generated using a=1.8; b=0.61; δ=0.001; α=0.00004; and λ=600 (solid line, hares; dotted line, parasites). For the limit cycle in state space, the velocity within the limit cycle is indicated by the length of the dashes, one per year.

(b) Parameter changes that generate dynamics with the observed properties

We reverse engineered changes to the modified parameter set that would reduce parasite loads while maintaining all other dynamical properties in the vicinity of those observed. Using elasticity analysis, the key parameters in determining equilibrium parasite load were identified as hare intrinsic fecundity (a), parasite-reduced hare fecundity (δ), hare intrinsic mortality (b) and parasite-induced hare mortality (α) with some interactions also being important. Figure 2 shows the four parameters plotted pairwise revealing that to generate stable limit cycles with the observed properties requires either one of two possible parameter set modifications, both of which require increasing a parameter outside its plausibility envelope set by empirical data. To generate observed dynamics, the effect of the parasite on hare fecundity (δ) can be increased by approximately 10-fold. Alternatively, hare intrinsic mortality (b) can be increased by approximately 0.8 adult hares per year (reducing mean hare lifespan by approx. 0.8 years) combined with a small increase in parasite-induced mortality (α) within the plausible envelope.

Figure 2

Parameter changes required to obtain dynamical properties observed in wild hare populations are given by the distance between the asterisk and the polygon. The asterisk is the position of the modified parameter set, the closest we get to observed dynamics within the empirically determined plausibility envelope, while the polygon represents the observed range of dynamics specified by the observed equilibrium hare densities (20–200 km−2), equilibrium parasite load (1000–3000) and period of 4–15 years (period contours indicated). Stability contours are shown (dashed lines: value of real part of dominant eigenvalue −0.1 (stable), 0, and 0.1 (unstable)) with stable limit cycles occurring at low positive values. Other parameters were held constant at the modified values. (a) α versus δ, (b) α versus b, (c) α versus a, (d) δ versus a, (e) b versus a, (f) δ versus b (α, parasite-induced hare mortality; δ, parasite-reduced hare fecundity; b, hare mortality rate; a, hare intrinsic fecundity).

As we will discuss, we believe that parasite-reduced fecundity (δ) and hare intrinsic mortality (b) may have been empirically underestimated. Increasing parasite-reduced fecundity (δ) from 0.0001 to 0.001 hare parasite−1 resulted in a qualitative change from a stable point to a stable limit cycle with a 15-year period (figure S1a). Increasing hare intrinsic mortality (b) from 0.61 to 1.40 year−1 (annual survival of 0.25–0.54) resulted in a stable limit cycle with a period of 18 years. Subsequently increasing parasite-induced mortality (α) to 0.000014 reduced the period of the limit cycle to 15 years (figure S1b). Increasing parasite-induced mortality (α) alone generated rapidly damped oscillations with a small period. It was not possible to obtain the observed population dynamics by changing hare intrinsic fecundity (a) alone.

For both sets of dynamics in figure S1a and b, the peak parasite loads were unrealistically high (105) which, if we assume that the parasite load at the peak of the cycle corresponds to maximum parasite loads counted in the field, should be approximately 16 000 worms per hare. Parasite loads of a more realistic amplitude were obtained by increasing parasite-induced mortality (α) above 0.00004 (figure 3a), which lies well within the plausibility envelope. Additionally, the simulated hare populations shown in figure S1a and b spend most years at numbers much below the lower observed limit for hare density. Changes in parasite fecundity (λ) and transmission inefficiency (H0) affected the amplitude of hare oscillations but not of parasite burdens (figure 3b,c). Thus, a set of parameters was identified, which produced realistic dynamics in both hare and parasite populations (figure 1c).

Figure 3

The effect of small, plausible changes in parameters on the amplitude of hare–parasite limit cycles. (a) Increasing parasite-induced hare mortality, α, from our point estimate of 0.000008 (solid line) to 0.00004 (dashed line) year−1 reduces the parasite load oscillation to below 16 000 worms per hare. The cycle shrinks further as α is increased to its upper plausible limit of 0.0001 (dotted line) year−1. Variation in (b) parasite fecundity, λ, and (c) transmission inefficiency, H0, controls the amplitude of the hare oscillation. λ and H0 were varied between lower and upper plausible limits (dotted lines). Solid lines are the limit cycles using point estimate values of λ and H0. Other parameters were kept constant: a=1.8, δ=0.001, b=0.61, α=0.00004, μ=0, k=1/2.

(c) Can small changes in parameters account for variability in dynamical properties?

Scottish populations of mountain hares exhibit a wide diversity of observed dynamics. We used the model to look for parameters that may vary across Scotland and affect period of cycles and amplitude of limit cycles within plausibly small changes in their value. Figure 4 shows the sensitivity analysis of stability and period to small changes in individually varied parameters around the system, which generated realistic dynamics (figure 1c). Stable limit cycles occurred where the system crossed the boundary from stable to unstable, and amplitude increased with increasing instability. Variation in hare intrinsic fecundity (a), parasite-reduced hare fecundity (δ), parasite-induced hare mortality (α) and adult parasite mortality (μ) could account for the range of periods observed in natural populations. Variation in all parameters in figure 4 except adult parasite mortality (μ) could account for variability in stability and amplitude, which occurs across the species range in Scotland. Finally, although the parasite transmission parameters (λ and H0) were not found to influence stability or period, the amplitude of the hare density limit cycle was sensitive to small changes in their value (figure 3).

Figure 4

The effect of small, plausible changes in parameters on stability and period. One parameter was varied at a time, while others were held at values for the system that generated realistic dynamics (figure 1c): parasite-induced hare mortality (α) was varied from its point estimate to upper plausible limit; hare intrinsic fecundity (a), hare intrinsic mortality (b) and adult parasite mortality (μ) were varied from lower to upper plausible limits; degree of overdispersion (k) values 0.5, 1, 2; δ was varied from its point estimate to 0.001 hare parasite−1 and follows the opposite path to α. The dashed line marks the boundary between a stable (negative) and an unstable (positive) equilibrium.

4. Discussion

This model of the mountain hare–T. retortaeformis interaction cannot predict observed population dynamics of mountain hares with realistic parasite burdens within the broad range of parameter space we judge to be plausible. We now discuss three possible interpretations of this observation. (i) Parasites are the main drivers of hare cycles, but the model, while including the key elements of the interaction, represents them insufficiently realistically. (ii) Parasites are the main drivers of hare cycles and the model has altogether omitted important ways in which the parasites influence hare demography. (iii) Parasites are indeed not the main drivers of hare cycles.

(a) Are key elements of the interaction represented sufficiently realistically?

To represent the hare–parasite system sufficiently realistically requires both adequate model parametrization and formulation. As several of our plausible parameter ranges were based on small sample sizes or indirect data sources, it is possible that our estimated parameter ranges are wrong. The key difficulty is to find a model where the dynamics are unstable and parasite loads realistic. Parasite burdens were particularly sensitive to the level of parasite-induced hare mortality (α), and our estimate was based on a single study (Newey & Thirgood 2004). However, the study found almost no difference in survival between parasite-reduced and untreated hares (Newey & Thirgood 2004) and increasing parasite-induced mortality (α) while lowering parasite burdens towards more realistic levels has a strong stabilizing influence on the resulting dynamics.

Parasite burdens were also sensitive to hare fecundity (a), parasite-reduced hare fecundity (δ) and hare mortality (b). Here, fecundity was used as the measure of recruitment and leveret pre- and post-natal mortalities were not included because the effect on parasites at these stages is unknown (but see next section). Mountain hares are killed for sport, pest and disease control, but this mortality is not included in the present analysis. The relationship between density and harvesting has not yet been studied in mountain hares (Newey et al. 2007a). Game harvesting bags are a good proxy for population abundance in red grouse (Cattadori et al. 2003), but there is likely to be more inconsistency across years in mountain hares. If, as suspected, the relationship between mountain hare density and harvesting is not density dependent, then the present model formulation holds and hare mortality rate (b) should be increased. Decreasing hare recruitment, strengthening parasite suppression of hare recruitment or higher hare mortality is destabilizing and it reduces parasite burdens; but to attain realistic dynamics, a large parasite-induced mortality rate (α) is still required.

If our estimate of parasite-induced hare mortality (α) is reasonable, then it seems unlikely that a key element of the parametrization has been omitted, and now we query the formulation of our model. Hare recruitment and parasite development were represented as purely continuous processes. Time delays and seasonality are well known to be destabilizing the population dynamics of infectious diseases (Greenman et al. 2004; Altizer et al. 2006), and both occur in mountain hares and the parasite T. retortaeformis. Mountain hares do not mature in their year of birth but in the following year, and the breeding season is restricted to approximately nine months of the year (Flux 1970). T. retortaeformis is a direct life cycle parasite; eggs voided in the host's faeces develop to an infective stage outside the host over a period of time which depends on climatic conditions (Crofton 1948). Although it is well established that the developmental time lag has a destabilizing influence on model host dynamics (May & Anderson 1978) the present model does not incorporate a time delay in parasite recruitment. In favourable conditions, development time is short and the assumption of negligible time delay in relation to changes in hare densities is reasonable. However, development may last several months over winter. We have explored discrete-time formulations of our model, which incorporated: a step function that restricted hare reproduction (at an accordingly increased rate) to a nine-month breeding season; a delay in the maturation of leverets until the start of their first breeding season; and a simple delay (that ranged between 1 and 12 weeks) in parasite maturation that was constant across the year. With these alternative formulations, damping times increase but we still do not recover sustained limit cycles within the plausible parameter ranges. However, we do not have to go as far outside these ranges as we do with the purely continuous-time formulation.

(b) Are important ways in which the parasites influence hare demography omitted?

The red grouse–T. tenuis system in Scotland has similar characteristic features as the mountain hare–T. retortaeformis system, such as low predation, greater parasite-reduced fecundity than survival and range of parasite burdens. Yet in applying a similar approach as here, Dobson & Hudson (1992) were able to reproduce grouse cycles. A striking difference is that parasite effect sizes (α,δ) were estimated at approximately 30 times greater for grouse than those estimated here. As parasite effect on fecundity (δ) was calculated from seven-week-old brood sizes, we now discuss the possibility that parasites may affect the aspects of hare recruitment other than the number of ova shed in females.

Parasites may have a larger impact on recruitment through influences on leveret survival, growth rate or timing of breeding. The timing of breeding is important for reproductive success in mammals (e.g. Clutton-Brock et al. 1982) and is influenced by parasite infection in a range of species (e.g. Feore et al. 1997). Time of first breeding in mountain hares is influenced by temperature, female age, size and weight with older, larger and heavier females attempting to breed earlier than those younger and smaller (Flux 1970). Young born earlier in the year have a longer growing season, enter the winter heavier and larger than late born young and have higher overwinter survival and greater future fecundity (Iason 1989, 1990). Thus, females may seek to breed earlier to produce more, higher quality young and we suggest future studies could profitably investigate a maternal effect of parasite load on the timing of breeding, survival and growth of leverets.

Maternal effects may destabilize population dynamics and promote cycles (Beckerman et al. 2002). To model this would require a new hare–parasite model formulation that could encapsulate: a maternal body size effect on the birthdate of young; a maternal parasite load effect on the birthdate of young; adult body sizes determined by birthdate; and adult hare mortality related to body size. This additional parasite-mediated effect may reduce the extent to which parameters need to deviate from our point estimates to generate realistic dynamics. However, they require a move from simple ordinary differential equation formulations of host–parasite dynamics to a partial differential equation model or individual-based approach, which is beyond the scope of the present paper.

Other forms of environmental variation that would lead to sufficiently large stochastic variation in the parameters of the host–parasite model at realistic frequencies could result in the generation of sustained limit cycles of a realistic magnitude from the damped oscillations predicted by the deterministic model. However, there is as yet no empirical data to inform the magnitude, covariation or frequencies of these stochastic processes.

(c) A secondary role for parasites?

If parasites are not the main driver of mountain hare cycles, could they still have a role in hare population dynamics? Parasitic nematodes of Soay sheep increase the depth of population crashes initiated by winter food shortage (Gulland et al. 1993). Similarly, reduction of parasitic nematodes from a red grouse population (Hudson et al. 1998) arguably does not remove a tendency to cycle (Lambin et al. 1999; Tompkins & Begon 1999), but parasites might deepen the extent of grouse crashes rather than determine their frequency. This notion is supported by our analyses, which showed that hare cycle amplitude was very sensitive to parasite transmission parameters whereas period was relatively insensitive.

(d) Conclusion

Despite the observation of large parasite burdens in mountain hares, and the perceived absence of predation and food limitation, we have found limited support for parasite-driven hare cycles. The results of our sensitivity analysis suggest that lower recruitment rates, stronger parasite suppression of recruitment or raised adult hare mortality than presently realized would allow a closer fit between model predictions and observed dynamics. Therefore, we identify leveret biology and the quantification of harvesting of hares for sport, pest and disease control as research priorities. If parasites do drive hare cycles, the model presented here suggests that our understanding of the full effects of parasites on hare demography is importantly incomplete.

Acknowledgements

Empirical data to parametrize the models was collected during a long-term study. This study has been approved by the Local Ethical Review of the Macaulay Institute. Capture and handling of mountain hares is conducted under licence from Scottish Natural Heritage. Parasite reduction experiments were conducted under licence from the UK Home Office.

S.E.T. was funded by a NERC PhD studentship. S.N. and S.J.T. were funded by the Scottish Government Rural Environment Research and Advisory Department. We are grateful to Russ Hobbs, Christina Cobbold, Margaret Kerlin, Douglas Kerlin, Isabella Cattadori, Kyrre Kausrud and an anonymous reviewer for their advice and comments on the manuscript.

Footnotes

    • Received November 14, 2008.
    • Accepted December 19, 2008.

References

View Abstract