INTRODUCTION
Nematode infections are a major threat to sustainable animal farming, and their control has been complicated further by the emergence of resistance to widely used anthelmintic treatments such as benzimidazoles (BZ) in herbivores (for review see Kaplan, 2004). The molecular basis of BZ resistance in nematodes implicates one recessive locus, with no apparent fitness disadvantage in the absence of treatment (Elard and Humbert, 1999). Hence, the best way of managing resistance is to delay the early increase in the frequency of the resistance allele. However, the practical consequences in terms of controlling the early phase of resistance selection have not been investigated. The first cases of anthelmintic resistance in parasitic nematodes were recorded 5 years after the introduction of BZ (for review see Mehlhorn, 2001), and so we investigated changes in low resistance allele frequency over 5 years.
Modelling anthelmintic resistance calls for both population genetics and epidemiological approaches (Paterson and Viney, 2000). Several mathematical models have been devised to explain how drug resistance develops, and how it can be slowed or even avoided. However, models of resistance selection often include simplifications that may not be without effect on their predictions. Firstly, compartment models average the effects of intra-host dynamics on both parasitic infection and resistance selection (Smith, 1990). Nematode infections are characterized by an over-dispersed distribution between hosts (Gaba et al. 2005). Interaction between density dependence and over-dispersion should influence the spread of anthelmintic resistant nematodes throughout the population. A relatively low nematode burden may have a substantial representation in the genetic composition of the transmission stages after treatment. Secondly, previous models of resistance evolution have generally assumed an invariable nematode transmission rate within the host population (e.g. Smith, 1990; Leathwick et al. 1995). However, infective stages of nematodes are heterogeneously distributed in clumps on the pasture and Tallis and Leyton (1969) have shown that the stochastic transmission of infective clumps exaggerates the aggregated distribution of nematodes among their hosts. Aggregation transmission might also increase the reproduction success of some infra-populations relative to others (Criscione and Blouin, 2005) and, stochastic transmission of infective clumps was shown to be a significant force in promoting the spread of resistant genotypes (Cornell et al. 2003). However, the model of Cornell et al. (2003) did not describe the dynamics of each developmental stage in the parasites' life-cycle including weather-dependent survival of free-living stages on pasture. The relationship between the development of the parasites and the timing of treatment was not taken into account, which may explain why they observed a rapid process during which the resistance allele either was fixed within the nematode population or became extinct within 3 or 4 grazing seasons, which is never observed in real conditions (for chronological data about first appearance of anthelmintic resistance see Mehlhorn, 2001). Finally, several modelling studies of anthelmintic resistance only set out to delay the spread of resistance alleles, without attempting to find out whether this decrease in frequency of the resistant allele would be accompanied by any reduction in the nematode burden of the hosts, e.g. (Dobson et al. 1987; Smith, 1990; Leathwick et al. 1995). However, since morbidity and mortality of hosts increase with the parasite burden (Poulin, 1998), sustainable control requires both a decrease in the nematode burden of the hosts and the least possible selection for drug resistance.
In the present paper, we report the development of an individual-based model, in which the infection process of each sheep consists of the random sampling of clumps of infective larvae. Individual nematode infestations in sheep display considerable variability (Shaw and Dobson, 1995; Gaba et al. 2005) and so average values are not truly representative. In contrast to mathematical models (using first-order differential equations) where the biological parameters of the sheep population are averaged, individual-based models explicitly model an individual sheep, and seem to be more appropriate for nematode infestation. All parasitic and free-living stages were explicitly modelled to take into account the interaction between the treatment effects and the parasite dynamics on the selection of resistance and on the control of infection. We studied both the genetic and demographic processes of a nematode infesting sheep, investigating the influence of the number of treatments and their timing on the early phase of resistance allele selection.
MATERIALS AND METHODS
We modelled naive lambs that grow old during the simulation, and acquire protection against nematode infection. For simplicity's sake, the hosts will be referred to as ‘sheep’ throughout this paper. A formal description of the individual-based model is appended (Appendix 1).
The grazing management scenario of our model was as follows. In mid-March (Julian week 10), a group of 30 naive sheep is turned out onto a 2·5 ha permanent pasture, contaminated with over-wintered infective larvae (106 per ha), and kept on the pasture until the end of November (Julian week 48); they were kept indoors during winter. The larvae were of 3 genotypes found in benzimidazole resistance (rr, resistant to anthelmintics, or rS and SS, susceptible to anthelmintics). The initial frequency of the anthelmintic resistance allele was 0·5% (P=0·005) as selected in other models (Gettinby, 1989; Barnes et al. 1995; Saul, 1995). The sheep were treated every year while leaving the sheepfold (week 10). The biological data (Table 1) used as the basis for the parameters in our model were derived from field data obtained in a temperate area for the parasitic nematode Teladorsagia circumcincta in sheep (Cabaret et al. 2002).
At the beginning of the simulation, 2·5×106 infective larvae were randomly distributed in clumps on the pasture. In the following years, the pasture infection originated from over-wintered infective larvae. A ‘clump’ consisted of a defecation produced by an individual sheep (Cornell et al. 2004). A clump harbours eggs that develop into larvae that are the offspring of the mature nematodes in a sheep (Cornell et al. 2004). As the number of resistant infective larvae was low at the beginning of the simulation (65 rr larvae out of a total of 2·5×106 larvae), only a few clumps contained resistant individuals. The nematode life-cycle was explicitly modelled in 6 stages: the first 3 free-living stages (egg-L1 confounded, L2 and L3) were developing on the pasture, and the parasitic next 3 stages (ingested L3, L4 and adults) developed in sheep. The established larvae were developing into sexually mature nematodes. Each stage lasted approximately 1 week, and so the time-step of the model is a week.
Each week, sheep on pasture randomly ingested a proportion (ρ, ingestion rate) of available clumps of infective larvae of mixed genotypes. Therefore the number of ingested clumps directly depended on the total number of clumps on the pasture. Each ingested clump was randomly selected from among the defecation clumps available on the pasture, leading to a variable transmission rate between the hosts (Fig. 1). In extensive grazing systems, with a moderately gregarious sheep flock, sheep may be infected by infective larvae derived from several defecation clumps produced by a given sheep, which is consistent with our scenario. The fraction of the ingested infective larvae that became established depended on the host's protective immunity. The immune response reduced nematode establishment rate and depended on the age of the host (Dobson et al. 1990a,b). We modelled the decrease in the establishment rate with the age of the sheep in a deterministic way, decreasing from 0·6 week−1 when the sheep were naive, to 0·47 week−1 for a 1-year-old sheep and 0·35 week−1 for sheep from 2 years of age.
Fig. 1. Diagram of the nematode-sheep model. The size of the excreted clumps varies with nematode burden in the hosts. The transmission of nematodes and the mortality rate of adult nematodes are both density dependent.
The probability that established infective larvae will develop into 4th-stage larvae (L4) was 0·9 week−1 and as was the probability that they would go on to form polygamous adults (0·9 week−1). We assumed a balanced sex ratio and random mating between adult nematodes within the host. Egg production (350 per parasite and per week, estimated from Cabaret and Ouhelli, (1984) depended on the total nematode burden (Bishop and Stear, 2000), so each additional nematode decreased the reproductive capacity of those already established. This crowding effect is coupled with a mating probability (as described by May (1977)). Adult survival is also affected by the crowding effect, and so we used a lower adult survival rate during the grazing period (0·88 week−1) than during the sheepfold period (0·94 week−1). We assumed that the administration of an anthelmintic would result in the death of all mature susceptible nematodes (SS and Sr). The frequencies of the rr, rS and SS genotypes in eggs produced by random mating were pi2, 2piqi and qi2, where pi is the frequency of the allele r of the mature nematodes in sheep i and qi=1−pi. Eggs are deposited on the pasture in defecation clumps in the form of ‘packets of eggs’. The size of the defecation clumps varies with the nematode burden in the sheep, and so the offspring of mature nematodes are assumed to be randomly distributed in ‘a’ defecation clumps, where ‘a’ is the number of defecation clumps produced in 1 week by a given sheep. All the eggs in the defecation clumps hatched into 1st-stage larvae (L1), and the probability that they would mature into 2nd-stage larvae (L2) was 0·9 week−1, and the probability that L2 would go on to develop into infective larvae (L3) was 0·9 week−1. The survival rate of infective larvae depended on environmental conditions and was lower during the summer (0·5 week−1 from weeks 28 to 40) than during the rest of the year (0·97 week−1). We used 3 criteria to identify the optimal nematode control strategy: the resistant allele frequency (to estimate the selection of resistance), the abundance of free-living stages and the nematode burden of the sheep (to assess the nematode infection rate).
A sensitivity analysis was performed to test the impact of changing the values of various parameters on the spread of resistance and on the infection level: the initial infection level of the pasture, the initial resistance allele frequency, the number of ingested clumps, the stocking rate, the number of defecation clumps produced/sheep/week (a), and the survival rate of infective larvae L3 during the summer. Except from flock size (we tested 10–50 sheep for constraints of time-consuming calculations instead of an average of 250–500 sheep in farms in France and Europe (Guesdon, 2003)), the values tested were all reasonable and matched those observed in French meat sheep farms. All the sheep were treated once a year, at week 10. We then investigated the effect of treating between zero and 14 times per grazing season on the spread of resistance, and on infection levels in the sheep and pasture. When more than 1 treatment was applied, the treatments were regularly spaced over the grazing season (between week 10 and week 48). We also investigated the influence of treatment timing by treating the sheep once per grazing period at different weeks (from week 10 to week 46 by steps of 4 weeks).
In total, 1000 simulations were performed for each scenario. Each simulation covered 260 weeks, equivalent to 5 grazing seasons in temperate areas. Results are divided in classes of final frequency of resistance allele, chosen in order to show rare events (genetic drift and loss of resistant allele, or fixation of resistant allele).
RESULTS
Analysis of sensitivity of the model
Initial resistance, allele frequency and initial level of infection of the pasture
For a given initial frequency of the resistant allele, a large variability of the final frequency of the resistant allele is observed, illustrating the random sampling of clumps by sheep. Genetic drift may result in the loss of resistance allele when the initial frequency of resistant allele does not exceed 10−3 (Table 2). For higher initial frequency of resistance allele, 1 treatment a year is sufficient to fix resistance in some cases (89 simulations out of 1000). The initial resistance allele frequency has no significant effect on the mean infection level of pasture or on the mean nematode burden.
The higher the initial infection of the pasture, the higher the mean density of larvae on pasture (Table 3). Conversely, the relationship between the initial level of infection of the pasture and the final frequency of the resistant allele or the mean nematode burden is a non-monotonic relationship (Table 3). This relationship results from the interaction between selection events and the density-dependent fecundity of parasites. A low initial level of infection of pasture leads to a low initial level of infection of sheep and, when treatment is applied, the mean nematode burden is too low to favour an efficient selection of the resistant allele. When time goes on, more clumps are produced leading to an increase of both the level of infection of the pasture and the mean nematode burden, with a moderate selection of the resistance. When the level of infection of pasture increases (or when the initial infection level is high), sheep harbour a larger mean nematode burden and hence more resistant parasites which survive to treatment. However, their contribution to the next generation is modulated by the density dependence constraints on fecundity, leading to fewer offspring per resistant nematodes.
Number of defecation clumps (a) and the selection of resistance
The number of defecation clumps per sheep and per week, ‘a’, has no significant effect on the average infection level of pastures or on the mean nematode burden of sheep (Table 4). Conversely, larger numbers of defecation clumps increase the probability of exclusive rare events (loss or fixation of resistance allele) and, on average, resulting in an apparent stable final frequency of the resistant allele.
Effect of the stocking rate
The greater the stocking rate, the greater the effect of selection (Fig. 2). As the area of pastures remained unchanged, when flock size increased, the stocking rate increased. For a given density of infective larvae, the probability that a sheep ingested clumps with rare resistant larvae increased. Moreover, for a greater number of sheep, the annual treatment will be performed on a larger number of nematodes, increasing the size of nematode population submitted to selection, and hence increasing the final frequency of the resistant allele.
Fig. 2. Effect of the stocking rate on the final value of the resistance allele frequency (r) after drug application for 5 years (1 per year). The results of 1000 runs for each value of stocking rate were divided into 4 classes: class 1 (0–0·005), in white, is characterized by having a frequency lower than or equal to the initial value (0·005); class 2 (0·005–0·1) in light grey, represents an undetectable resistance level with phenotypic diagnosis; class 3 (0·1–0·9) in dark grey, represents a treatment failure identified by phenotypic diagnosis, and class 4 (0·9–1) in black, corresponds to a very high level of resistance. For example, when the stoking rate is 4 sheep/ha, the final frequency of the resistance allele ranges between 0 and 0·005 in 40% of simulations, and between 0·005 and 0·1 in 60% of simulations.
Intensity of selection by treatments
The intensity of selection was investigated using the following parameter values: a flock of 30 sheep producing 50 defecations/sheep/week that resulted in clump of infective larvae.
Number of treatments
The number of treatments had a strong effect on resistance selection: the higher the number of treatments, the higher the final frequency of the resistance allele (Fig. 3A). However, this effect was highly modulated by the timing of the treatment (as shown below).
Fig. 3 (A) Final value of the resistance allele frequency (r) after 5 years of different regimens of drug application. The results of 1000 runs for each number of clumps were divided into 4 classes as described in Fig. 2. (B) Average nematode burden (solid black line) and average density of infective larvae on the pasture (solid grey line) over 5 years of different drug application regimens. The average has been computed with the results obtained for 1000 runs. The dotted lines represent the 95% interval of confidence.
Without any treatment, random sampling of infective larvae leads to genetic drift by which the resistant allele could be lost. Surprisingly, the absence of treatment did not fail to control sheep or pasture infection, but led to infection levels close to those recorded when numerous treatments per year were undertaken. Conversely, 1 to 2 treatments per year (applied on week 10 or weeks 10 and 48) strongly reduced sheep and pasture infection (Fig. 3B). More than 3 treatments per year exerted a strong selection, systematically leading to the fixation of resistant allele, impeding an efficient control of sheep infection. Consequently, infection levels on pasture and in sheep reach their highest values (4·5ċ106 larvae/ha and 7ċ103 nematodes) (Fig. 3B).
Timing of treatment
The timing of treatment markedly affected the early stages of resistance selection (Fig. 4A and B): the treatment effect was modulated by the weather (lower survival rate of infective larvae in summer). With a constant survival rate of infective larvae during the summer, the effect of the treatment was limited as a result of the re-infection of the sheep during the pasture period (Fig. 4B). Conversely, with a lower survival rate of infective larvae during the summer, treatments applied during autumn (week 42 and week 46) favoured the spread of resistance (Fig. 4A): sheep were kept indoors after the treatment and adult nematodes that survived treatment made a large contribution to the contamination of pasture at the beginning of the following spring.
Fig. 4. Effect of the timing of treatment with a reduced survival rate during summer (A) and with a constant survival rate (B) on the resistance allele frequency (r) after 5 years. (A) The survival rate of the infective larvae is lower (s3=0·5) during the summer (from week 28 to week 40). (B) The survival rate of infective larvae remains constant (s3=0·97) throughout the grazing season. The results of 1000 runs for each time of treatment were divided into 4 classes, as described in Fig. 2.
The lower survival rate of infective larvae during the summer shows that the size of the refuge (provided by the density of susceptible larvae on the pasture) was a significant factor in the development of anthelmintic resistance (Fig. 4A). Treatments applied when survival of infective larvae was low (from week 28 to week 40, during the summer) lead to greater selection for resistance than treatments applied during the period of favourable weather. Lower survival rate of larvae during the summer reduced the density of infective larvae on pasture before the sheepfold period. As a consequence, the re-infection of sheep treated a few weeks before the sheepfold period (week 42 and 46, Fig. 4A) was reduced and resistant nematodes made a high contribution to the contamination of pasture at the next grazing period. Taking into account the strong effect of differential survival rate of infective larvae, we tested the reliability of our results to several survival rates (data not shown). A marked effect of treatments applied during the summer on the resistance selection was observed for survival rate from 0·7 to 0·8 week−1. However, when the survival rate of infective larvae decreased dramatically (0·4 week−1), few infective larvae (susceptible and resistant) survived, hence reducing both the probability of infection and of ingestion of resistant larvae.
DISCUSSION
We developed an individual-based model to investigate the impact of genetic and demographic processes in a sheep/nematode life-cycle on the early spread of a rare resistance allele. Our model allowed quantification of the contribution of those factors in the early evolution of the resistance gene: we confirmed the importance of the number and the time of treatments, in the selection of resistance gene. We found that the stocking rate that was known to favour high nematode intensity in sheep or cattle (Hansen et al. 1980; Grenfell et al. 1987; Aumont et al. 1992) was also a force in promoting resistance.
Unlike Cornell et al. (2003), we did not observe a rapid extinction of the nematode population within 3 or 4 grazing seasons. The selection of the resistance allele is slower in our model in comparison with Cornell et al. (2003). Certain differences between these two models may account for this discrepancy. First, in the study of Cornell et al. (2003), the number of infective larvae ingested resulted either in the nematode population becoming extinct (never observed) or in extremely high nematode burdens (106 nematodes/lamb after 1 grazing season), whereas, in our model, nematode burdens were less variable and consistent with field data for a temperate area (Euzeby et al. 1964; Mangeon, 1986; Chartier and Reche, 1992; Hong et al. 1996). The differences in the assumption of the modelling approaches may explain this discrepancy. Modelling the parasite life-cycle with 6 stages (3 free-living and 3 parasitic, each transition between stages during 1 week) versus 2 stages (time step of the day) in their model entails a delay in the spread of resistance and in the host infection because of the longer developmental time required for offspring to infectious stages and to sexual adults. Secondly, in the study of Cornell et al. (2003), each year a new cohort of naive lambs was turned out to pastures. In contrast, we monitored the same sheep for 5 years, which allowed us to introduce immunity by making older sheep less susceptible to infection as they get older. Even if our approach is somewhat approximated (immunity is a complex phenomenon), this factor cannot be ignored when investigating nematode population dynamics and the spread of resistance: Saulai et al. (2001) showed that immune protection may strongly reduce parasite population. Our scenario, in which the sheep graze the same area for 5 years, attempts to simulate nematode infection on natural extensive pastures.
As expected, frequent administration of anthelmintic treatments reduced the time taken to reach a high level of resistance. Nevertheless, the stochastic ingestion of infective larvae highlights the impact of genetic drift at the beginning of anthelmintic selection when resistance selection is low (1–2 treatments per year). Frequent treatments are necessary but not sufficient in selecting for resistance: without differential survival rate of infective larvae during the summer, the probability that selection for resistance occurs is strongly reduced (less than 100 simulations out of 1000). We used 0·5 and 0·97 per week as the survival rate of infective larvae during the summer and the rest of the year, respectively. This corresponds to a mean life-span of 1 week and 25 weeks, respectively, which is in agreement with field data from temperate areas (Smith and Grenfell, 1985). One may expect that with higher rainfall, the survival rate of infective larvae may be greater and the size of the population in refuge would be larger than in our results show. As a consequence, we should observe a still higher effect of the refuge on the delay of drug resistance selection.
The absence of treatment leads to a moderate infection of sheep (mean nematode burden of 6000 nematodes/sheep) and pasture (4×10e6 larvae/ha), together without any drug resistance spread. This raises the question of treatment necessity under certain circumstances. Treatments eliminate parasitic stages and favour the establishment of newly ingested larvae. The migration of larvae throughout the mucosa is the most pathogenic aspect of T. circumcincta infection. As a consequence, prophylactic treatments should be avoided when feasible (medium to low infection) in order to maintain the balance between the intensity of the infection and its control by sheep. It has already been observed on experimental paddocks (unpublished data on T. circumcincta infection in 1st-year grazing lambs) that the multiplication of treatments did not increase the bodyweight gain and did not control for nematode intensity in comparison with an untreated group of sheep. These considerations are only reliable under medium to low sheep nematode infection. In our model, the sheep acquire a protective response to infection during the first two grazing seasons, which is commonly observed in ewes that may often remain untreated in farms. The model showed that the contribution of resistant parasites to the next generation is highly dependent on the size of the population in refuge, e.g. the free-living stages not submitted to treatment selection. A small refuge favours the selection of the resistance allele, but a larger one is not sufficient to completely avoid the resistance selection, due to the random sampling of clumps of infective larvae. The size of the refuge varies with season, being smaller during the summer when the mortality of infective larvae is high. Therefore, a sustainable management strategy involves a trade-off between preserving a population of susceptible larvae in refuge on pastures to slow down selection of the resistance allele, and keeping the number of infective larvae low enough on pasture to reduce sheep nematode infection. The size of the population in refuge is an important factor, as has already been demonstrated in insecticide resistance management (Lenormand and Raymond, 1998) and in pest resistance management in crops (Vacher et al. 2003). Modelling may help us to identify an optimum usage of refuge: rotation of pastures may be explored to optimize the size of populations in refuge to control both resistance selection and nematode infection, especially since the numbers of possibilities make it impracticable to test these features experimentally.
This model was intended for use in investigating the spread of the resistance alleles during the early phase of selection. In this study, we did not integrate the ewes' reproduction, but we monitored a cohort of lambs over 5 years in order to include the acquisition of protection against nematode infection. Nearly all the data were based on those established in lambs, which display higher nematode burdens and have the largest contribution to the production of infective larvae (Donald and Waller, 1973). This simplification should not affect the qualitative findings reported here. Our model is flexible and could possibly incorporate other aspects as, for instance, an age-structured host population including breeding management (ewes with differing physiological status, pregnant or lactating) as the immune protection of ewes is known to relax at the pre and post-lambing period (Barger, 1993). Both periods are high risk periods for nematode infection. We should expect that accelerated lambing (3 lambing periods per 2 years or 5 lambing periods per 3 years) may further increase pasture infectivity and sheep infection. The model was developed to correspond with the benzimidazole resistance mechanism but we may expect that it would be different for polygenic resistance mechanisms (e.g. levamisole or macrocyclic lactones).
S. Gaba is the grateful recipient of a Ph.D. grant from the INRA Epidemiology Modelling Program. The model is implemented with a computer program written with R. Copies of the program and instructions for use are available from A. Silvestre, the corresponding author. The English text has been checked by Monika Ghosh.
APPENDIX
The model
The dynamics of parasites may be described by the following set of equations. The equations 1 to 3 describe the dynamics of the free-living stages, and the equations 6 to 8 the dynamics of the parasitic stages. The effect of the clumping process is described by the equation 4 and 5.
with:
Parasite fecundity
Parasite fecundity depends on:
- Density-dependence
- Egg output from one sheep decreases with worm burden. The crowding effect can be expressed as:with α the parameter characterizing the severity of density-dependent constraints.
- Effective mated females
- Let denote the distribution of the number of mated females in host i at time t. The number of mated females depends on parasite sexual habits (here: dioecious polygamous) and on parasite distribution. We supposed that this distribution is aggregated (for more details see May, 1977). We assumed that the sex ratio is balanced in hosts.
Therefore, parasite fecundity at date t can be expressed as:
Egg genotypes are calculated from Hardy-Weinberg mating rates. The probability of hatching eggs of genotype j is Fi,j[ratio ]pi2 eggs of genotype j=1(rr), 2pi(1−pi) eggs of genotype j=2(rS) and (1−pi)2 eggs of genotype j=3(SS). pi is the frequency of the resistant allele r in host i.
