Portal de Periódicos da CAPES

A modified matrix model to describe the seasonal population ecology of the European tick Ixodes ricinus

2011; Wiley; Volume: 48; Issue: 4 Linguagem: Inglês

10.1111/j.1365-2664.2011.02003.x

1365-2664

Andrew D. M. Dobson, Thomas Finnie, Sarah Randolph,

Vector-Borne Animal Diseases

Journal of Applied EcologyVolume 48, Issue 4 p. 1017-1028 Free Access A modified matrix model to describe the seasonal population ecology of the European tick Ixodes ricinus Andrew D. M. Dobson, Andrew D. M. Dobson Department of Zoology, University of Oxford, Oxford, OX1 3PS, UKSearch for more papers by this authorThomas J. R. Finnie, Thomas J. R. Finnie Health Protection Agency, Porton Down, Salisbury, SP4 0JG, UKSearch for more papers by this authorSarah E. Randolph, Corresponding Author Sarah E. Randolph Department of Zoology, University of Oxford, Oxford, OX1 3PS, UK Correspondence author. E-mail: [email protected]Search for more papers by this author Andrew D. M. Dobson, Andrew D. M. Dobson Department of Zoology, University of Oxford, Oxford, OX1 3PS, UKSearch for more papers by this authorThomas J. R. Finnie, Thomas J. R. Finnie Health Protection Agency, Porton Down, Salisbury, SP4 0JG, UKSearch for more papers by this authorSarah E. Randolph, Corresponding Author Sarah E. Randolph Department of Zoology, University of Oxford, Oxford, OX1 3PS, UK Correspondence author. E-mail: [email protected]Search for more papers by this author First published: 18 May 2011 https://doi.org/10.1111/j.1365-2664.2011.02003.xCitations: 77AboutSectionsPDF ToolsRequest permissionExport citationAdd to favoritesTrack citation ShareShare Give accessShare full text accessShare full-text accessPlease review our Terms and Conditions of Use and check box below to share full-text version of article.I have read and accept the Wiley Online Library Terms and Conditions of UseShareable LinkUse the link below to share a full-text version of this article with your friends and colleagues. Learn more.Copy URL Summary 1. The sheep tick Ixodes ricinus is the most multicompetent vector in Europe, which is responsible for significant diseases of humans and livestock throughout the northern hemisphere. Modelling the tick’s complex seasonal dynamics, upon which pathogen transmission potential depends, underpins the analysis of tick-borne disease risk and potential tick control. 2. We use laboratory- and field-derived empirical data to construct a population model for I. ricinus. The model is a substantially modified stage-classified Leslie matrix and includes functions for temperature-dependent development, density-dependent mortality and saturation deficit–meditated probability of questing. 3. The model was fitted to field data from three UK sites and successfully simulated seasonal patterns at a fourth site. After modification of a single parameter, the model also replicated divergent seasonal patterns in central Spain, but any biological factors underlying this geographical heterogeneity have not yet been identified. The model’s applicability to wide geographical areas is thus constrained, but in ways that highlight gaps in our knowledge of tick biology. 4. Sensitivity analysis indicated that the model was generally robust, particularly to changes in density-independent mortality values, but was most sensitive to changes in parameters related to density-dependent mortality. 5. Synthesis and applications. Vector population models allow investigation into the effects of individual environmental factors on population dynamics in ways not easily possible by experimental manipulation of in situ populations. Our model can be used to evaluate public health risk, tick management strategies and potential effects of future environmental change. Introduction The sheep tick Ixodes ricinus (Linnaeus) is a haematophagous ectoparasite of considerable medical and veterinary importance. Its ubiquity within Europe and Eurasia, abundance and catholic use of vertebrate hosts make it an effective vector for a large number of bacterial, viral and protozoan infections. In the United Kingdom alone, it is responsible for transmitting the agents of Lyme borreliosis, human granulocytic anaplasmosis (HGA), babesiosis and Louping ill. Lyme borreliosis is the most common vector-borne disease of humans in the northern hemisphere and may also affect domesticated animals. A quantitative understanding of the complex life cycle of vector ticks must underpin any serious attempt to predict disease risk and investigate potential control strategies. Marked environmental seasonality and the consequent adaptive diapause commonly seen make modelling far more complex for temperate species than for tropical species. In the United Kingdom, peak abundance of I. ricinus nymphs and adults is seen in April/May, followed approximately 6 weeks later by the peak in larvae. All three stages display a lesser resurgence in numbers in late summer. These basic patterns vary geographically according to local environmental conditions, with significant consequences for the tick’s vectorial capacity, especially for pathogens with short duration of infectivity in the vertebrate (Randolph et al. 2000). The same is true of I. scapularis in North America (Ogden et al. 2007; Gatewood et al. 2009). Each instar of I. ricinus (larva, nymph, adult female) takes a single blood meal from a separate vertebrate host. Larvae and nymphs typically feed for 2–4 days, whereas adult females require c.10 days and increase their body weight 100 times. The processes of development and moulting between instars takes several months, during which time ticks must survive on nutrients from this single meal. Microclimatic factors impose strong seasonal constraints on tick population dynamics, not least because most of the life span is spent off-host, either actively questing for hosts on vegetation or sheltering in moist microhabitats on the ground to avoid desiccation during development and seasonal diapause (Randolph 2004). I. ricinus may undergo both behavioural and developmental (or morphogenetic) diapause (Belozerov 1982). Behavioural diapause is typically defined in the acarological literature as a preprogrammed period of functional inactivity in unfed ticks, whereas developmental diapause is a physiological state in which the ticks’ bodies undergo biochemical changes, particularly related to a slowing of metabolic processes (Belozerov 1982; Randolph et al. 2002). Both forms of diapause allow ticks to minimise energy use during times when it would be climatically unfavourable to quest for hosts. In I. ricinus in temperate Europe, behavioural diapause occurs to variable degrees through autumn and winter and developmental diapause occurs as a precursor to development when ticks feed in late summer or autumn (reviewed in Randolph 2004). In this paper, we present the first biological, process-based population model for I. ricinus. It incorporates all available empirical information on the biological processes that control I. ricinus populations, including environmental determinants of diapause initiation and cessation (Randolph 2004), and the probability of questing (Alekseev et al. 2000; Perret, Rais & Gern 2004), day-degree summation modelling of development (Randolph et al. 2002) and estimates of density-dependent mortality (Randolph 1997). The model is constructed using data from three separate field sites in the United Kingdom where detailed tick abundance and microclimatic variables were measured over a period of several years (Randolph et al. 2002), and tested against data from independent sites in the United Kingdom and in central Spain. A sensitivity analysis is conducted to assess the model’s robustness to deviations in individual parameter values. One of the most recent models of Ixodes spp ticks (Ogden et al. 2005) captured seasonality well and offered insights into potential range expansion of I. scapularis in Canada (Ogden et al. 2006), but none has been applied to I. ricinus. Our model is characterised by its distinct underlying mathematical approach, protracted egg-laying phenology, use of saturation deficit (as opposed to temperature) to mediate questing probability and simulation of development periods by incremental day-degree summation. Materials and methods Model description General The tick model, which simulates relative numbers of ticks of each instar and physiological stage, is written in R language (R Core Development Team, 2005) and based around a stage-classified Leslie matrix (Leslie 1945; Caswell 2001), with some significant modifications. It moves by daily increments; ticks in each stage experience a daily probability of moving to the next stage (Fig. 1), although certain transitions incur fixed time delays to simulate such processes as feeding and post-eclosion cuticular hardening. Active processes such as temperature-dependent development and egg-laying, which take more than a single day to complete, are also simulated in multielement vectors outside the matrix. Figure 1Open in figure viewerPowerPoint Schematic outline of part of the model. The full model is repeated three times, once for each instar. Autumn questing and behavioural diapause (with subsequent spring questing) are mutually exclusive alternatives. Ticks enter developmental diapause if they have fed after a certain date in late summer (see text); thus, the proportions vary depending upon the development rate of the previous stage and the host contact rate. Labelled processes: FL, Larvae feeding (combines host encounter rate with probability of actively questing, mediated by saturation deficit); LD, Engorged larvae entering developmental diapause; LE, Larvae exiting developmental diapause; DevL, Engorged larvae beginning development; AQN, Developed nymphs questing in autumn (following a 20-day cuticular hardening period); BDN, Developed nymphs going into behavioural diapause in autumn (following a 20-day cuticular hardening period); QN, Nymphs exiting behavioural diapause to begin questing; FN, Nymph feeding (combines host encounter rate with probability of actively questing, mediated by saturation deficit); ND, Engorged nymphs entering developmental diapause; NE, Nymphs exiting developmental diapause. Developmental diapause Developmental diapause delays the onset of development if ticks experience shortening day length (Belozerov 1998) by feeding after a certain date in late summer. Empirically, at three latitudinally graded sites across Europe, the observed date of diapause onset (see Randolph 2004 for details) coincided with the period (e.g. end of July, mid August, etc.) at which day length shortens at a maximal rate, measured according to when the centre of the sun dips 9° beneath the horizon (‘tick twilight’). This is applied to the model. Because brief exposure to cold (at least 0 °C) breaks diapause (Campbell 1948) and development is determined by temperature above a threshold of c.7·5–9 °C (depending on the stage), as long as ticks have emerged from diapause before spring warming starts, the timing of emergence (onset of development) is not critical and is controlled by a Gaussian probability function that gradually releases ticks from 6 February onwards, a date chosen to accommodate these criteria adequately. Development Daily increments of development are modelled as a function of temperature with day-degree summation, using the following stage-specific equations derived by Randolph et al. (2002) from empirical observations of Campbell (1948), where T is daily mean temperature: From egg to emerged larva: (eqn 1) From engorged larva to emerged nymph: (eqn 2) From engorged nymph to emerged adult: (eqn 3) Development is complete when the cumulative development reaches one, and eclosion at each instar is followed by a 20-day period of cuticular hardening. Questing and behavioural diapause Ticks with the potential to quest (i.e. not in diapause or development) are considered to be ‘questing’, although, importantly, this is a nominal state in the model and does not imply that they are necessarily actively questing for hosts. Probability of active questing is incorporated into the probability of contacting a host. The fat content of field ticks demonstrated that, whether or not they have undergone developmental diapause, ticks emerge for questing in a largely synchronised fashion each year in the late summer/early autumn (Randolph et al. 2002). These ticks then either quest at this juncture or undergo a period of behavioural diapause. It is not known what determines which individuals start to quest during the autumn rather than entering behavioural diapause, and here the choices are modelled as mutually exclusive. The proportions that go into either autumn questing (0·1) or behavioural diapause (0·9) were estimated by reference to the relative heights of spring and autumn peaks of abundance and are constant for all UK field sites described in this paper. Ticks that undergo behavioural diapause over winter emerge to quest the following spring (Randolph et al. 2002), thought to be triggered by increasing temperatures; the threshold for onset of questing is set at a smoothed weekly average daily maximum temperature of 7 °C for nymphs and adults and 10 °C for larvae, in line with field observations (Perret et al. 2000; Randolph 2004). In years when this temperature index does not drop below 7 °C, the date of activity onset is set at 22 December, the shortest day; increasing day length is thought to stimulate questing where temperatures are sufficiently high (Randolph 2004). In reality, an increasing proportion of the tick population may start to quest as more of its microhabitat reaches the threshold temperature, so a Gaussian function is applied after the threshold date. The parameter values of this function were calibrated against field data so as to capture the timing and rate of spring appearance of questing ticks, but are nevertheless constant across the three UK field sites. Host density and feeding probability Effective host density is defined by the stage-specific host contact rate, represented as the probability that a tick in a questing position will make contact with a host (Table 1). This is separated into the relative contributions from large and small hosts: small hosts are programmed to feed 0% of adults, 55% of nymphs and 95% of larvae, according to limited empirical data from woodland habitats (Talleklint & Jaenson 1994), but on upland moorlands in the United Kingdom, large hosts (predominantly sheep) feed >90% of immature stages (Steele & Randolph 1985; Ogden, Nuttall & Randolph 1997). Here, we consider an individual’s probability of finding a host, not the proportional rate of change in the questing population. Host contact rates were expected to vary between field sites, and values were scaled according to observed differences in tick density at the three sites used in model fitting (the only such site-specific parameter values). Table 1. Model parameter values, specific for each field site only where shown. Values for the two validation sites are as per Exmoor A tick’s feeding probability is calculated as the probability of its active questing (i.e. being positioned on the vegetation ready to respond to a passing host) multiplied by the host contact rate. The probability of active questing is modelled as a function of saturation deficit (SD) (eqns 4 and 5). Data on the effects of SD on the probability of questing by nymphs (also applied to larvae in lieu of anything more suitable) and adults are based on the proportion of ticks seen questing in field arenas over several months (Perret, Rais & Gern 2004) but do not account for tick mortality, so the data were first corrected, modelling mortality as an exponential decay curve. The SD values of Perret, Rais & Gern (2004) were calculated from relative humidity and temperature measured at a nearby meteorological station rather than by a data logger at the site. These values were adjusted using the relationships between relative humidity and temperature at the one of our field sites (Dorset) and a nearby meteorological station. The relationship between the corrected proportions of questing ticks and adjusted SD was used to generate the functions: (eqn 4) (eqn 5) The feeding probability (F) for each stage can therefore be described by the generic equation: (eqn 6) where pq = the probability of questing shc = small host contact rate lhc = large host contact rate Host contact rates for each stage, combined for both host sizes, are listed in Table 1. The feeding period (3 days for larvae and nymphs, 10 days for adults) is incorporated into the model as a fixed delay. Density-dependent mortality is applied during this period. Egg-laying A fully engorged adult female I. ricinus is assumed to lay 2000 eggs (Gray 1981). Egg laying is modelled as this fixed number of eggs distributed over a 60-day period according to the egg-laying curve for Rhipicephalus appendiculatus (Randolph & Rogers 1997), shown in eqns 7 and 8. Adult females are removed from the model after laying, when they die. (eqn 7) (eqn 8) LV is a 60-element vector that contains the factors by which to multiply the 60-element laying adult vector to generate the correct number of eggs laid per female per day. The term x(1,60) describes the vector elements 1–60; the number of eggs laid by each female varies as she moves sequentially down the vector. The constant k, which depends upon the values of r1 and r2, ensures that the maximum number of eggs laid over 60 days does not exceed 2000; the function is applied to adults of both sexes, so the maximum number of eggs is halved to account for a 50:50 sex ratio (eqn 8). The constants r1 and r2 determine the steepness of the climb (r1) and decline (r2) of the egg-laying curve; we use the same values as Randolph & Rogers (1997), who adjusted the values to fit the curve to that of Branagan (1973). Mortality Daily density-independent mortality rates are applied to all stages. Background rates are based upon field observations on ticks held in permeable containers in the Czech Republic (Daniel et al. 1976) and Ireland (Gray 1981), but were modified post hoc to ensure that they were low enough to allow population equilibrium at the climatically least favourable of the field sites used for validation (Powys, Wales; see below for description of sites). To this end, all the published mortality values were multiplied by a factor of 0·68. Mortality rates were not available for unfed ticks in behavioural diapause, which may suffer lower predation from arthropod predators than engorged ticks (Gigon 1985; Samish & Alekseev 2001) and be less susceptible to freezing stress (Dautel & Knulle 1997), but may be more vulnerable to desiccation owing to their higher surface-area-to-volume ratio and the relatively low levels of deposition of waxy lipids in their cuticles (Yoder, Selim & Needham 1997). In the absence of empirical data, the most parsimonious solution is to assume equal background mortality and apply the same rates as for the engorged stages. Mortality rates vary with the degree of desiccation stress determined by atmospheric SD but increase only when SD rises beyond approximately 3 mm Hg (Mount, Haile & Daniels 1997). This was virtually never recorded in the tick habitats monitored in the United Kingdom, so this term was omitted rather than relying on the only available empirical parameters and relationships derived for a different species (I. scapularis) on a different continent (N America). In addition to the background daily mortality rate, questing duration is limited to 60 days in larvae and 120 days in nymphs and adults, whether or not preceded by behavioural diapause, according to the approximate maximum life span of unfed ticks imposed by the metabolic costs of questing (Steele & Randolph 1985). Besides mortality through desiccation, starvation and predation, all ticks are likely to experience density-dependent mortality, potentially mediated via host immunity (Sutherst, Wagland & Roberts 1978; Randolph 1994; but see Ogden et al. 2002), host grooming (Levin & Fish 1998) and host avoidance of infested areas (Sutherst et al. 1986). In East and South Africa, relatively simple demographic processes and seasonal recruitment in tropical climates allowed density-dependent effects to be estimated for R. appendiculatus from the relationships between interstadial mortality indices (k-values) and tick density (Randolph 1997). For example, in South Africa, the overall density-dependent relationship showed a slope of 0·387 and intercept of −0·814 during the larva-to-nymph interstadial period and a slope of 1·24 and intercept of −2·579 during the nymph-to-adult period. These data came mainly from farms where cattle act as the dominant host. A more varied host assemblage might produce weaker density-dependent effects if some hosts are less likely to display tick avoidance, grooming and/or acquired resistance (Randolph 1979). Where rodents feed a greater proportion of ticks, less density-dependent tick population regulation might be expected than in a cattle-only system, owing to the rapid turnover of hosts relative to the timing of acquired immune response. Because the I. ricinus model was designed to simulate environments with diverse wildlife host assemblages (and because ‘mortality’ in the African study also included density-independent mortality, which is accounted for separately in our model), the slope of density dependence was set below one and the intercept below zero (i.e. a positive threshold tick density). Within these constraints, the slope and intercept for the larva/nymph transition were iteratively adjusted to produce matches with the observed data; the nymph/adult values were fixed to retain the same proportional relationship to the larva/nymph values as in Randolph’s (1997) study, with the same values applied to the adult/egg transition in the absence of any other empirical data. Density-dependent mortality is applied during the engorgement period according to the generic eqns 9 and 10 for the probability of survival (S); individual values for slope and intercept for the three stages are given in Table 1: (eqn 9) (eqn 10) where kv = ‘k-value’; i = intercept of line of best fit in relationship between kv and log of population size (see Randolph 1997); s = slope of line of best fit in relationship between kv and log of population size; qp = active questing population; y = a scaling value inversely proportional to the host contact rate. Higher host contact rates (i.e. host density) should cause concomitant reductions in density-dependent mortality, because the on-host density of a fixed number of ticks is reduced when the host density is increased. Field data Primary data (density of each tick stage) for model construction came from standard sampling for I. ricinus by dragging a blanket over the ground, carried out at approximately fortnightly intervals between 1997 and 2000 at three sites in the United Kingdom: Powys, Wales, and Dorset and Somerset (Exmoor), south-west England (Randolph et al. 2002 for sites and methods). Data for independent model validation were provided by blanket sampling on Exmoor in a similar habitat to above at 3-weekly intervals from March 2008 to December 2009 (Dobson, Taylor & Randolph 2011) and for a site in central Spain where ticks were sampled monthly for 9 years (Estrada-Pena et al. 2004). Climate data From April 1996 to May 2001 in Dorset, April 1997 to September 2001 on Exmoor and April 1997 to September 2000 in Wales, hourly temperature (ground level and 30 cm above) and relative humidity (30 cm above ground level) at the three field sites were recorded using automatic weather stations (based on Squirrel data loggers; Grant Instruments, Cambridge, UK). Hourly data were converted to daily minimum, maximum and mean figures. Saturation deficit at 30 cm above ground level was calculated from temperature and relative humidity (equation in Randolph & Storey 1999). One gap in the Dorset data owing to logger malfunction was filled by calibrated data from local meteorological stations. For model simulations, actual climate data for each site were preceded by a 10-year series in which each day was given a value equal to the mean for that date across the full observation period for that site; this allowed population establishment. Comparisons with observed data were based on model simulations using the actual data after these 10 years of mean values. Exits from behavioural diapause are triggered by threshold values in the smoothed weekly average of the daily maximum temperatures taken from meteorological stations closest to the field sites and smoothed by loess interpolation. Climate data (temperature and relative humidity at questing height) for the Exmoor test site were taken during the field visits. Comparisons between these values and those from a nearby meteorological station were used to calibrate the latter, which were used as the input in model testing. Weekly climate data recorded for the field data from central Spain (Estrada-Pena et al. 2004) were converted to daily values by simple linear interpolation. Model evaluation Empirical validation Model outputs were generated for the three field sites according to their specific climate data. Other parameter values (Table 1) were held constant across all sites apart from host contact rates, which reflect estimated indices of varying host abundance. The goodness-of-fit between the predicted seasonal patterns of numbers of ticks of each stage and those observed on each date of field sampling over 3 years was assessed according to the R2 value and slope of the regression line (Randolph & Rogers 1997). The same procedure was used to test simulations for the two independent field sites (Exmoor and central Spain). In both cases, host contact rates were specified as per the original Exmoor site. The two sorts of diapause prevent the estimation of interstadial development periods from observations on field populations. Therefore, we cannot compare the observed and predicted relationships between interstadial mortality rates and population density, as did Randolph & Rogers (1997). Simulations for the independent Exmoor site are shown for the years of data collection (2008–2009), but those for Spain (a 9-year series) are presented as mean values per sample/model day across all years, because of an observed strong upwards trend over the 9 years that was demonstrably unrelated to climate (Estrada-Pena et al. 2004). Sensitivity The sensitivity of model predictions to deviations in parameter values was investigated following a previous method (Randolph & Rogers 1997). Models were run with each parameter value in turn alternately increased and decreased by 50%, except where parameters consisted of dates (e.g. the centre date of the curve controlling exit from behavioural diapause), which were moved forwards and backwards by 14 days. The intercepts of the density-dependent mortality functions were not altered, because they merely pertain to the hypothetical unit of density of the model simulations, and are thus effectively arbitrary. Models were run with each altered parameter in turn for a 13-year simulation using climate data from the Dorset data logger, with each day’s values comprising the mean values from that date in each of the years for which data were recorded. Using only the final 5 years of simulations, the effect on R2 values and slopes in regressions of each changed model against the final model was recorded. Changed models that returned mean population sizes lower than the final model give regression slopes between zero and one and are presented in 7, 8 as negative reciprocals to give equivalent magnitude of change as in models where means are higher (when slopes exceed one). A slope of 0·2 therefore becomes −5, indicating that the changed model’s mean was five times lower than that of the final model. Figure 7Open in figure viewerPowerPoint Sensitivity analysis: reduction in regression R2 value relative to final model after 50% increase (filled bars) or decrease (open bars) in parameter values. Parameter definitions: Controls on Gaussian function controlling exit from developmental diapause: EDbA, Adult curve centre date; EDaA, Adult curve height; EDbN, Nymph curve centre date; EDaN, Nymph curve height; EDbL, Larva curve centre date; EDaL, Larva curve height. Controls on Gaussian function controlling exit from behavioural diapause: EBbA, Adult curve centre point; EBaA, Adult curve height; EBbN, Nymph curve centre point; EBaN, Nymph curve height; EBbL, Larva curve centre point; EBaL, Larva curve height. Fraction questing in autumn: FractA, Adults; FractN, Nymphs; FractL, Larvae. Mortality values: EgS, Eggs; EngA, Engorged adults; DevDiaA, Adults in developmental diapause; BehDiaA, Adults in behavioural diapause; QuesA, Questing adults; EngN, Engorged nymphs; DevDiaN, Nymphs in developmental diapause; BehDiaN, Nymphs in behavioural diapause; QuesN, Questing nymphs; EngL, Engorged larvae; DevDiaL, Larvae in developmental diapause; BehDiaL, Larvae in behavioural diapause; QuesL, Questing larvae. Host contact rates: HCA, Adults; HCN, Nymphs; HCL, Larvae. Density-dependent mortality regression values: sA, Adult slope; sN, Nymph slope; sL, Larva slope. Figure 8Open in figure viewerPowerPoint Sensitivity analysis: change in regression slope relative to final model after 50% increase (filled bars) or decrease (open bars) in parameter values. Positive numbers are equal to the slope; negative numbers are equal to the negative reciprocal of the slope; thus, the axis reflects the factor by which slope is either multiplied or reduced. Parameter labels as in Fig. 7 legend. Population growth rate calculations involving matrix eigenvectors (Caswell 2001) are not suitable for this model owing to the number of life-history events (notably density-dependent mortality and