Modelling the impact of targeted anthelmintic treatment of cattle on dung fauna

The insecticidal properties of many anthelmintics pose a risk to dung fauna through the e ﬀ ects of drug residues in dung on the activity, oviposition and development of dung-dwelling invertebrates. Reductions in dung fauna numbers can inhibit dung degradation, which may impact biodiversity and nutrient cycling on farms. A simu- lation model was created to predict the impact of antiparasitic drugs on cattle dung fauna, and calibrated using published data on the dung-breeding ﬂ y Scathophaga stercoraria . This model was then tested under di ﬀ erent e ﬀ ective dung drug concentrations (EC) and proportions of treated cattle (PT) to determine the impact under di ﬀ erent application regimens. EC accounted for 12.9% of the observed variation in S. stercoraria population size, whilst PT accounted for 54.9%. The model outputs indicate that the tendency within veterinary medicine for targeted selective treatments (TST), in order to attenuate selection for drug resistance in parasite populations, will decrease the negative impacts of treatments on dung fauna populations by providing population refugia. This provides novel evidence for the bene ﬁ ts of TST regimens on local food webs, relative to whole-herd treatments. The model outputs were used to create a risk graph for stakeholders to use to estimate risk of anthelminthic toxicity to dung fauna.


Introduction
Anthelmintic drugs are widely and routinely administered to grazing livestock to control gastrointestinal nematodes and other parasites. Anthelmintics are typically not fully metabolized within the host animal and residues of the drugs are often excreted in dung (McKellar et al., 1993) (and urine (McKellar, 1997)) and can therefore exert non-target effects on invertebrate fauna which spend part, or all, of their life cycle in dung (Floate, 1998a(Floate, , 1998bGover and Strong, 1995;Madsen et al., 1990;Sommer et al., 1992;Sutton et al., 2014) and also on soil invertebrates (Scheffczyk et al., 2016). Such effects include inhibited motility, oviposition, emergence, and reduced dung pat colonisation (Floate, 1998a(Floate, , 1998bGover and Strong, 1995;Suarez et al., 2003). Invertebrate dung fauna significantly contribute to the degradation of dung through physical processes and therefore reductions in the activity and populations of degradative fauna can to slow dung degradation (Madsen et al., 1990;Wall and Strong, 1987) with potential knock-on effects on important local processes., including local ecology (Beynon, 2012;Strong, 1993;Wall and Beynon, 2012) and epidemiology. In recent years, the mounting resistance of gastrointestinal parasites of domestic livestock to anthelmintic drugs has led to a shift away from whole-herd treatments, and recommendations for targeted selected treatment (TST) (Charlier et al., 2014) of only part of the herd. This strategy aims to generate refugia from drug exposure among parasite populations, slowing the development of resistance. In principle, refugia from drug residues ought also to be generated for dung fauna, supporting their populations; however, to date no systematic attempts have been made to evaluate this possibility.
The ability to assess and predict the impact of anthelmintics and other routine veterinary medicines on the wider environment is essential for informed drug development and policy in agriculture. In particular, parasite control practices that slow the development of resistance to commonly administered anthelmintics are essential to sustainable livestock production systems. However, the scale and complexity of the drug-dung-fauna system is challenging to observe and quantify in vivo and is difficult to fully represent under controlled laboratory conditions. Modelling techniques are the best alternatives to address these issues by allowing for the manipulation of a wide range of variables specific to individual field scenarios, and rapid assessments of the potential impacts of new parasite control and other management practices on dung fauna. Boxall et al. (2007) developed a screening index for assessing the impact of veterinary medicines on dung flies. The index was simple and allowed for estimates to be calculated with relatively small amounts of data, allowing for rapid screening of multiple drugs. The index assessed impact by multiplying three variables: proportion of cattle treated, proportion of time of faunal contact with dung, and dung toxicity. A central assumption was that the three variables are equally weighted, but this assumption inadvertently creates a potential mathematical ceiling to drug toxicity. Vale and Grant (2002) took a different approach in their development of a model to assess the impact of insecticide-contaminated dung on dung fauna. The model considered a broad and novel range of variables including the response to distinct adverse ecological events on insect life cycle stages and dung-insect interactions which aided the understanding of the importance of refugia for the ecology different species of invertebrates.
Here, we test the hypothesis that the proportion of cattle treated (PT) with anthelmintics has a greater influence on Scathophaga stercoraria populations than the strength of drug residue in dung (EC). We build on previous theoretical and modelling approaches to create a new modelling approach to simulate the drug-dung-fauna system and evaluate the potential impacts of antiparasitic drug use in grazed cattle production systems. We use the model to consider how varying treatment regimens administered by veterinarians for the purpose of livestock health and welfare have non-target influences on dung invertebrates, and to provide a risk graph to inform stakeholders in sustainable livestock production systems.

Model description
A simulation model was created using NetLogo 5.0.4 (Wilensky, 1999) to estimate the impact of a hypothetical anthelmintic that expressed insecticidal properties when excreted in dung by cattle in a grazed field, upon a model dung invertebrate. A 2-dimensional virtual pasture system was created, occupied by a herd of cattle and a population of the model invertebrate. All actions and interactions presented were simulated hourly time-steps for each individual cattle or invertebrate, as appropriate.

Model components
The model simulated the interaction between a model dung invertebrate and cattle defecation behaviour, and the potential for invertebrate survival to be changed by different concentrations of anthelmintic residues in the dung.
The model invertebrate was the yellow dung fly Scathophaga stercoraria. The model utilized published data (Table 1) to simulate the life cycle of S. stercoraria in a temperate cattle grazing system. Scathophaga stercoraria is a well-studied dung fauna species, for which detailed information on life cycle parameters is widely available. The species is highly abundant across the northern hemisphere, and some of its life cycle stages are dependent on dung.
The model cattle were based on published data on temperate grazing commercial beef and dairy herds (Table 1). There were two components to cattle behaviour: (1) defaecation frequency, and (2) randomized movement across a field. The cattle were treated or untreated with a hypothetical anthelmintic, producing toxic or non-toxic dung, respectively. The proportion of cattle treated (PT) ranged from 0 to 1 in increments of 0.1 and was specific as an independent variable in each simulation, The rate of defecation of model dung by the model cattle and its mean carrying capacity for S. stercoraria was based on published data for temperate commercial beef and dairy systems ( Table 1). The model dung were toxic or non-toxic. The strength of the toxicity, i.e. effective concentration (EC) ranged from 0 to 1 in increments of 0.1 and was specific as an independent variable in each simulation. The dung became unattractive for S. stercoraria regardless of toxicity after a simulated 120 h.
A starting population of 100 individuals of S. stercoraria, covering a random distribution of ages within typical life expectancy for S. stercoraria, were simultaneously introduced to the system. They actively sought out cattle dung in order to produce off-spring with no preference for toxic or non-toxic dung. Population fitness responses of the S. stercoraria to contact with toxic dung was based on the interaction between PT and the specific EC.
Primary assumptions were: (i) the model dung toxicity retained a constant toxicity for 120 h (ii) there were no sub-lethal effects of the anthelmintics upon S. stercoraria (iii) there were no other sources of mortality exist for S. stercoraria other than toxicosis or exceedance of life span (iv) the population of S. stercoraria is isolated.
No values or weightings of variables within the model were assumed or given arbitrary values.

Application of modeling approach
The model was run 605 times. Each run simulated 4380 h (6 months) using all combinations of 11 PT values and 11 EC values, totaling 121 unique sets of parameter values. There were five repeats of each set, with variable outcomes depending on values simulated from normal distributions: the mean of each set of repeats was used for statistical analyses. The Anderson-Darling normality test was conducted on residuals for the dependent variable of final population size at the end of the simulated period to ensure appropriateness for parametric testing. This was followed by Pearson's correlation analyses of final population size versus PT and EC. Multiple regression analyses were then conducted to attribute how much of the variation in final population size was due to PT and EC, respectively.
A number of individual paired simulations were run to evaluate the index created by Boxall et al. (2007). These simulations were performed in pairs in which the product of PT and EC were equal, but the individual values of PT and EC in each pair were not equal. To achieve this the values for PT and EC of pair 1 were switched to form pair 2 ( Table 2). For the Boxall et al. (2007) model to agree with the presented model, there should be no significant different between pairs that meet the aforementioned assumptions. Final population numbers from  Martin et al. (2004), 5. Aland et al. (2002), 6. Gary et al. (1970), 7. Oudshoorn et al. (2008), 8. Sahara et al. (1990), 9. Villettaz Robichaud et al. (2011), 10. Floate (1998, 11. Vale and Grant (2002), 12. Geiger (2010), 13. Parker (1970 simulations were then subject to the Paired T-test.

Results
The distribution of final population sizes across all simulations was non-normal (Anderson-Darling, p = < 0.005). The data shows two distinct groupings based on final population size, one at 0 and the other in the region of 3100-4300 (Fig. 1). This latter group, the 'maximum fitness' group, had a normal distribution (p = 0.383). Quartiles for the maximum fitness group were measured as Q0 = 3259, Q1 = 3597, Q2 = 3703, Q3 = 3798, Q4 = 4197. PT and EC combinations that resulted in final populations of < Q0, and therefore outside of this group, were considered as high risk. Combinations that fell between Q0 and Q1 were considered medium risk, and all over combinations resulting in final populations > Q1 were considered low risk (Fig. 2).
In general, incremental increases in PT and EC at low levels had little effect on final population size (no. of individuals), but a tipping point was reached beyond which the population decreased exponentially (Table 3). Rising EC values from 0.0 to 0.5 brought about gradual decreases in final population size; however, as EC exceeded 0.5 its effect on population size reduced. In contrast, rising PT values of 0.0-0.5 had little impact upon population sizes, but as PT exceeded 0.5 there was a rapid drop in population size.
The residuals of the complete data for all experimental runs were normally distributed (as tested by Anderson-Darling test, p = 0.281) and thus no transformation was required for parametric analyses. A Pearson's correlation analysis showed that final population size was significantly correlated with PT (−0.694, p = < 0.001) and EC (−0.336, p = < 0.001). A subsequent multiple regression calculated the total variance of final population explained by PT and EC together, R 2 , to be 67.8% (p < 0.001). Further individual regressions showed that PT explained 54.9% (p = < 0.001) and EC explained 12.9% (p = < 0.001) of total variance in final population size.
The paired t-tests, for the purpose of evaluating Boxall et al., showed a statistically significant intra-pair difference (t = 2.43, p = 0.023) and therefore the H 0 was rejected in favour of the H a , that there is an intrapair difference. That is: simulations of which the sum of PT and EC are equal do not yield equal results.

Discussion
In this study, we used a novel simulation to test the hypothesis that PT had a greater impact on the population size of S. stercoraria than EC. The outcomes of 605 simulations of 121 parings of PT and EC indeed indicate that this hypothesis can be accepted. The distribution of data predicted that populations of S. stercoraria were generally resilient and can maintain stable numbers up until a tipping point at which mortality becomes probable. As such, our model develops the concept of the screening level index (Boxall et al., 2007) through simulation modelling using published data about key life cycle parameters that could strongly influence drug-insect interactions. We propose that this new approach provides a better justified mechanistic framework for impact assessment, which will improve recommendations of use of veterinary medicines with consideration for livestock dung ecology and wider impacts on the environment.
Cow pats in grazed systems without drug residues may provide an important reservoir of biodiversity, allowing maintenance populations of coprophagic fauna that are important for ecosystem services including nutrient cycling, carbon cycling and soil quality, e.g. dung beetles and insect larvae. Therefore, TST, as opposed to a whole-herd treatment, is recommended to reduce the impacts of drug treatment on local ecosystems, with additional economic benefits through reduced inputs on-farm (Charlier et al., 2012), ecosystem service delivery, and ensuring sustainable parasite control options in the longer term through slowing of anthelmintic resistance. However, in our model, S. stercoraria populations were assumed to be isolated, but wider consideration of the spatial variation in the local food web including the availability of drug-free dung in the wider environment would more closely represent the complexity of farming systems.
Our model provides a framework that is adaptable to dung-breeding insect species other than S. stercoraria. Its application to other target Table 2 Values of PT and EC for paired simulations in order to evaluate Boxall et al. (2007)  A.S. Cooke et al. Environmental Toxicology and Pharmacology 55 (2017) 94-98 species, however, would require further empirical information on the toxicity of various drugs, as faecal residues, on specific fauna. Moreover, life cycle parameters specific to other species would be required, although the model could also be used to explore parameter space and identify broad characteristics of species that are likely to be vulnerable to anthelmintic residues in dung, and the extent to which these might be attenuated by TST. Since the model framework was developed using a bottom-up approach, it lends itself to constructive adaptation and expansion. With sufficient observational data there is scope for future models, within such a framework, to increase in complexity and realism. Expansion of the model to represent multiple invertebrates at farm level would enable holistic landscape-scale impact assessments and attenuation strategies. The use of veterinary medicines, with non-target insecticidal properties, is ubiquitous and therefore the applicability of observed results may be equally wide, and this model framework can be adapted to any system, anywhere, given workable parameter estimates. The data from the literature that provided the foundation for the model was predominantly derived from studies in temperate regions, and we recognize that climatic differences may have a significant impact on the ecotoxicity of such medicines (Kryger et al., 2005). Moreover, the model was based upon a set-stocked system but could be modified to represent more extensive systems, including ranch/range rearing of cattle in the USA, South America, and Australia. The model provides a framework for the development of future similar work and could be applied to scenario-testing using the specific characteristics of different cattle production systems across the world.
Despite the high profile, global threat of drug resistance, the longterm impacts of drugs, especially antiparasitics with non-target insecticidal properties, are largely unknown. The topic is a key area for future work to enable effective assessment and regulation of the use of veterinary medicines, with regards to their impact on all aspects of biodiversity (Adler et al., 2016). Future work should also include economic analysis, in order to balance short-term production gains with longer term environmental impacts. There is likely to be a utilitarian argument to use veterinary medicines in a more sustainable manner, including the utilization of preventative and non-pharmaceutical methods (Kaplan and Vidyashankar, 2012;Papadopoulos, 2008;Fig. 2. Risk thresholds for the impact of anthelmintics on S. stercoraria. "Low Risk" refers to PT (proportion of treated cattle) and EC (effective concentration) combinations that result in final populations exceeding Q1 of the maximum fitness group, "Medium Risk" to those falling between Q0 and Q1 and "High Risk" to those falling below Q0.  Wolstenholme et al., 2004). The emergence of part-herd anti-parasitic treatments, or TST, is an example of a more efficiently targeted approach to chemical utilization in agricultural systems, which has potential long-term economic benefits, as well as reduced environmental impacts. The current model shows this synergy in quantitative terms for a model insect species, and provides a framework for impact assessment and optimization of TST strategies across a wider range of dung fauna, including those of conservation relevance.

Author contributions
All authors contributed to drafts of the paper. EM devised the study and advised on model assumptions and parameter values. AC compiled and ran the model, analyzed results and led the writing of the paper. JD supervised AC and provided editorial support.