An analysis of the spatio-temporal occurrence of anthelmintic veterinary drug residues in groundwater

• Anthelmintic veterinary drugs are emerging groundwater contaminants


H I G H L I G H T S
• Anthelmintic veterinary drugs are emerging groundwater contaminants of concern. • In a spatial survey, anthelmintic residues were detected at 18% of groundwater sites. • Thirteen of seventeen detected residues belonged to the broad spectrum benzimidazoles. • Transformation products were more prevalent in groundwater than in surface waters. • Temporal variations at karst springs are linked with drug usage and meteorology.

G R A P H I C A L A B S T R A C T a b s t r a c t a r t i c l e i n f o
Anthelmintics are antiparasitic drugs used to control helminthic parasites such as nematodes and trematodes in animals, particularly those exposed through pasture-based production systems. Even though anthelmintics have been shown to be excreted into the environment in relatively high amounts as unmetabolized drug or transformation products (TPs), there is still only limited information available on their environmental occurrence, particularly in groundwater, which has resulted in them being considered as potential emerging contaminants of concern. A comprehensive study was carried out to investigate the occurrence of 40 anthelmintic residues (including 13 TPs) in groundwaters (and associated surface waters) throughout the Republic of Ireland. The study focused on investigating the occurrence of these contaminants in karst and fractured bedrock aquifers, with a total of 106 sites (88 groundwaters and 18 surface waters) samples during spring 2017. Seventeen anthelmintic compounds consisting of eight parent drugs and nine TPs were detected at 22% of sites at concentrations up to 41 ng L −1 . Albendazole and its TPs were most frequently detected residues, found at 8% of groundwater sites and 28% of surface water sites. Multivariate statistical analysis identified several source and pathway factors as being significantly related to the occurrence of anthelmintics in groundwater, however there was an evident localised effect which requires further investigation. An investigation of the temporal variations in occurrence over a 13 month period indicated a higher frequency and concentration of anthelmintics during February/ March and again later during August/September 2018, which coincided with periods of increased usage and intensive meteorological events. This work presents the first detections of these contaminants in Irish groundwater Science of the Total Environment 769 (2021) 144804

Introduction
Pasture based animal production systems play an important role for sustainable farming in many temperate regions throughout Europe, including Ireland, with grass grazing providing a low cost, on-farm feed source for industries such as milk production (Hennessy et al., 2020). However, such systems are not without their disadvantages in terms of animal health, with pasture grazed animals (primarily cattle and sheep) highly susceptible to infection by helminthic parasites, which can be an economic burden on the food production system. According to Animal Health Ireland (2016), gastrointestinal nematodes (GIN) (roundworms), and trematodes (flukes) are the most economically detrimental groups of helminths infecting cattle and sheep livestock in Ireland. Typical issues associated with infection of livestock with roundworms or flukes includes poor appetite and feed intake resulting in poor digestion and absorption of nutrients, ultimately causing negative impacts on growth rate, fertility and production yield (e.g. milk or wool quality) (Miller et al., 2012;Animal Health Ireland, 2013). Consequently, anthelmintic agents have become a critical component in animal husbandry and pasture-based production systems (Patten et al., 2011;Bloemhoff et al., 2014).
There are several classes of anthelmintic drugs that can be used to treat these helminthic parasites, each class differing by their mode of action, with some classes more effective toward particular species of helminths than the others. The 40 anthelmintic residues (including 27 parent drugs and 13 transformation products) included in this study, subdivided into the main anthelmintic classes, are listed in Table 1, with a broad classification of the primary usage of each drug also provided. In cattle and sheep production, the most important of these classes are the benzimidazole group (BZs), macrocyclic lactone group (MLs) and the imidazothiazole levamisole (LEV), given their broad-spectrum efficacy toward many helminths during most of their life-cycle (i.e. early immature, immature and adult parasites). Flukicides, a collective of drugs from across the different classes, used specifically to target liver flukes in intensive livestock production, are also of importance.
Although the classes of drugs used are similar, anthelmintic drug usage patterns are complex and can differ not only between production systems, but also from farm to farm. Recommended best practices are based around management systems which adopt the approach of dosing according to the onset and severity of infection, as informed by regular testing for faecal egg counts. However, despite best efforts, traditional chemoprophylactic approaches involving a regular dosing programme, regardless of the infection status of the animals, are still prevalent on many farms (Patten et al., 2011). Consequently, resistance of certain helminths to anthelmintics such as the BZs, has become an emerging issue, particularly in sheep, with reported treatment failure attributed to resistance (Keegan et al., 2017). This in turn adds burden on treatment regimes, often with higher doses necessary to achieve the required efficacy. Anthelmintic usage and application patterns are largely centred around protecting calves and lambs (due to their lack of immunity) during periods of highest parasite burden (Bloemhoff et al., 2014), which are largely driven by climate and ground conditions. A summary of the main seasonal events of cattle and sheep production (i.e. housing, turn-out to pastures and grazing), around which anthelmintic dosing occurs in Ireland, is provided in Fig. S1 (Supplementary file 1). Boxall (2010) provides a comprehensive overview of the factors controlling both the entry and movement of veterinary drugs such as the anthelmintics in the environment. Once administered, typically as an oral drench, topical pour-on, or parenteral injection, the drugs can undergo transformation by metabolism in the animal, the extent of which depends on both the drug and route of administration. For some (e.g. ivermectin (IVER) and abamectin (ABA)), extensive metabolism does not normally occur with a large proportion (often >90%) of parent drug shown to be excreted in urine or faeces (Liebig et al., 2010;Beynon, 2012;Horvat et al., 2012). For others (e.g. fenbendazole (FBZ), morantel (MOR) and Closantel (CLOS)), more extensive metabolism has been reported (Wardhaugh, 2005). In some instances when the drugs are metabolized, the excreted metabolite(s) are often more active and toxic than the parent drug themselves (Danaher et al., 2007). Entry of anthelmintic residues into the environment is therefore primarily as a result of the direct excretion of the parent drugs and/or metabolites onto pastures during grazing, or as a result of the land-spreading of manure and slurry containing anthelmintic residues, typically accumulated in large amounts during the housing period (Boxall et al., 2004). The latter of these exposure routes is likely to be of most importance; several studies have highlighted the potential for anthelmintic drugs to persist in the manure due to minimal or slow degradation during storage (Floate et al., 2005;Kreuzig et al., 2007). Application of stored manure following the housing period can therefore provide a concentrated source of anthelmintic residues in the environment with initial landspreading. Once in the environment, anthelmintics (both parent and excreted metabolites) can further undergo degradation by processes such as oxidation, hydrolysis and photolysis (Horvat et al., 2012). As a result, breakdown products present in the environment (as a result of both excretion and environmental degradation) are collectively referred to as transformation products (TPs).
On entry to the environment, anthelmintics have the potential to be transported between different environmental compartments, including soil, surface water and groundwater. The extent to which anthelmintics are mobile in the environment is controlled by a combination of factors such as the sorption behaviour of the substances (dictated by their physicochemical properties) and environmental factors such as soil characteristics, climate conditions and hydrogeology (e.g. aquifer type and properties) (Boxall, 2018). Based on their physicochemical properties (Table 1), with low water solubility and high sorption coefficients, the anthelmintics are likely to be found more associated with soil and sediment than water. Reflecting this, the majority of analytical methodologies and previous occurrence studies have focused on soil or sediment (Moreno-Gonzalez et al., 2015), with only few on manure leachate (Raich-Montiu et al., 2008) and surface waters (Bartelt-Hunt et al., 2009;Sim et al., 2013).
To date, there is an apparent lack of research focusing specifically on groundwater occurrence, which may be as a result of the expected immobility of anthelmintics in the environment due to their affinity for soil. However, several studies have demonstrated the transport of anthelmintics via particle associated surface transport (Kreuzig et al., 2007) and preferential macropore flow (Weiss et al., 2008). Such pathways are of particular importance in hydrogeological settings dominated by secondary permeability, with flow primarily via fissures and fractures (Swartz et al., 2003;EPA Ireland, 2010). Furthermore, karstic aquifers possess solutionally widened fractures and openings, and unique features such as swallow holes and sinking streams, which allow for rapid point recharge to groundwater with minimal attenuation (Coxon, 2011). Groundwater bodies within such settings are therefore potentially vulnerable to contamination by anthelmintics.
Groundwater is not only important as a natural drinking water reservoir, with up to 75% of European Union habitants relying on groundwater for their drinking water supply (European Commission, 2019), it also plays a vital role in supporting and maintaining many groundwater dependent ecosystems, on which anthelmintic drugs can have a nontargeted effect.
While some anthelmintics have been shown to exhibit toxicological properties such as teratogenicity, embryotoxicity, neurotoxicity, goitrogenicity and mutagenicity (Kinsella et al., 2009), the expected levels within environmental waters are not likely to be of major concern for human health (ACVM, 2019), certainly lower risk than those expected in foods of animal origin, for which safe maximum residue limits have been set. However, the main concern over the environmental occurrence of anthelmintics is due to their ecotoxicological effects on nontarget organisms (Floate et al., 2005). The ecotoxicological effects of some classes of anthelmintics, such as the macrocyclic lactones, have been well documented, with negative effects observed for dung beetle populations (O'Hea et al., 2010;Jacobs and Scholtz, 2015) and different aquatic organisms (Sanderson et al., 2007;Liebig et al., 2010). Information on other classes is lacking (Wagil et al., 2015a), with ecotoxicological information on TPs even more scant, as highlighted by Horvat et al. (2012).  Danaher et al., 2007;Krogh et al., 2008a,b;Horvat et al., 2012;Santaladchaiyakit and Srijaranai, 2012;Popova et al., 2013;van der Velde-Koerts, 2014;Zrncic et al., 2014) where Sw = water solubility, logKow = logarithm of octanol-water partition coefficient, pKa = dissociation constant and logKoc = logarithm of soil organic carbon-water partitioning coefficient. Early work by the Boxall research group highlighted gaps in current research regarding the environmental fate and ecotoxicity of veterinary medicines, with several anthelmintic drugs deemed of potential high priority in terms of risk to the environment (Boxall et al., 2003;Boxall et al., 2004). On assessing more recent reviews on this topic, it is still evident that these gaps have yet to be filled (Horvat et al., 2012;Kaczala and Blum, 2016;Snow et al., 2016). A European Medicines Agency (EMA) reflection paper (CVMP, 2016) noted that a lack of information prevented a conclusive decision on 20 substances identified as potentially persistent bioaccumulative toxic substances (PBTs), the majority of which were parasiticides used in food-producing species. The anthelmintic moxidectin was deemed suitable for classification as a PBT. More recently, as part of the concluding remarks of a report reviewing the method of supply of anthelmintic veterinary medicinal products carried out by the Health Product Regulatory Authority (HPRA), it was also noted that there is a deficit in the monitoring of antiparasitic drugs in the environment (ACVM, 2019). Monitoring of these drugs in the environment is also important in order to assess the potential impacts of climate change, which is likely set to increase the usage of anthelmintic drugs due to increased parasite burden associated with extreme weather events (Morgan and Wall, 2009;van Dijk et al., 2010;Bloschl et al., 2019).
Considering the above, the main aim of this study was to apply a newly developed comprehensive analytical procedure (Mooney et al., 2019) to investigate the spatial occurrence of anthelmintic drugs in environmental waters, particularly groundwaters, throughout Ireland. This study particularly focused on groundwater within karstic and fractured bedrock aquifers, given their prevalence throughout Ireland, and their potential to accommodate inadequately attenuated transport of contaminants to groundwater, through flow paths that are likely to be of most importance for the anthelmintics. The work also aimed to provide an insight into any seasonal controls on anthelmintic occurrence, the understanding of which is critical for assessing potential future variations in occurrence caused by the effects of climate change.

Materials and methods
2.1. Sampling sites 2.1.1. Spatial study site selection This study was carried out as part of a project which contributed to a larger research challenge focused on securing and protecting groundwater within karst and fractured bedrock aquifers. In Ireland, karst aquifers are subdivided into two categories depending on the extent of conduit and diffuse flow within the system (DELG/EPA/GSI, 1999). As a result, groundwater sampling sites were mainly selected to investigate the relevance of the following three flow regimes to anthelmintic occurrence, with a representative number of sites selected broadly across these three regimes: (a) karst aquifers dominated by conduit flow (Rkc), (b) karst aquifers dominated by fracture diffuse flow (Rkd) and (c) non-karst bedrock aquifers with fracture flow. Samples of several associated surface waters were also taken, as explained below. Besides this main selection criterion, other factors that were taken into account were pre-existing research on sites (notably sites from the Teagasc Agricultural Catchments Programme (ACP) and from joint Geological Survey Ireland (GSI) and National Federation of Group Water Schemes (NFGWS) work), stakeholder interest and site ownership/access considerations.
All groundwater sites were characterised based on the different land-uses (source factors) and hydrogeological properties (pathway factors) within their zone of contribution (ZOC), with the predominant class of each property recorded and input for statistical analysis (as described in Section 2.3). ZOCs for most of the sampling sites have been previously delineated using hydrogeological mapping and/or numerical modelling incorporating (but not limited to) data on topography, groundwater tracing experiments, recharge, extraction and previously determined boundaries, as described by Kelly (2010). ZOC boundary shapefiles for all Environmental Protection Agency Ireland (EPA Ireland) groundwater monitoring sites were downloaded from publicly available datasets (EPA Ireland, 2018), while ZOCs for private group water schemes (GWSs) were provided confidentially by the NFGWS in coordination with the GSI. ZOCs for several research observation piezometer boreholes were estimated in coordination with Teagasc ACP researchers.
Overall, 106 sites comprising 88 groundwaters (which included 41 boreholes (BHs) and 47 springs (SPs)) and 18 surface waters (SW) (which included 11 streams/rivers, 2 lakes, 3 swallow holes (SHs) and 2 drainage ditches), were sampled across the Republic of Ireland during March and April 2017 (Fig. 1). This sampling was selected to coincide with a period of active groundwater recharge and following the return of animals to pasture and the recommencement of manure spreading ( Supplementary Fig. S1). A large proportion (42%) of sites were sampled in coordination with the EPA Ireland in tandem with the national groundwater quality monitoring programme for the EU Water Framework Directive (WFD) (EPA Ireland, 2019). The remaining sites comprised private and/or semi-private group water schemes (GWSs) (28%) sampled in coordination with the NFGWS or observation/research sites (31%) sampled in coordination with Teagasc ACP. In total, 66% of groundwater sites were used for public drinking water supplies (including springs and BHs), with the remainder used for agricultural or observational purposes. Overall, 35% of groundwater sites selected fell within karstic conduit flow aquifers, with 24% and 41% of sites falling under karstic diffuse flow aquifers and non-karst fractured bedrock aquifers, respectively.
While this sampling network represents a variety of land use and hydrogeology, it is not necessarily representative of all groundwaters in Ireland, given the focus of the study was on karst and fractured bedrock aquifers. Some areas were also sampled at a catchment level, resulting in several small clusters of sampling sites. The surface waters included are not intended to be representative of surface water quality in Ireland, rather these surface waters were sampled to complement some groundwater sampling in areas of karst groundwater-surface water interactions and in ACP catchments.

Temporal study site selection
Following the spatial investigation in March/April, a temporal study was undertaken starting in November 2017. Temporal sampling was focused on karst groundwater with rapid flow and short travel times, carried out at a catchment level in two regions of Ireland and comprised 11 sampling sites, which included four karst springs and three swallow holes (sinking streams) in County Roscommon, and four karst springs in County Clare (Fig. 2). These sites were selected in consultation with the Irish hydrogeological community and based on pre-existing knowledge of, and research at, the study sites. Consideration was also given to such sites which had detections of anthelmintics during the spatial study. Based on previous dye tracing experiments, the three sinking streams sampled were individually linked to three of the springs sampled. Sampling of these sinking streams was therefore carried out to complement the sampling of the spring associated with the sink. Each site is summarised in Table S1 (of Supplementary information file 1), in terms of the predominant land-use and hydrogeological properties within their ZOC (characterised using the approach described in Sections 2.1 and 2.3). All 11 sites are associate with conduit-flow dominated karst aquifers (Rkc) which provides the potential for rapid underground travel, with short residence times (Drew, 2008). Groundwater velocities across both Roscommon and Clare sites typically range from 19 to 224 m h −1 , with travel times ranging from a few hours to a few days e.g. positive dye tracing experiments to Clare Site A showed groundwater velocities of 163-224 m h −1 , resulting in a travel time of 6-8 h, while similar studies in Roscommon showed a velocity of >157 m h −1 from a swallow hole (Roscommon site F) linked to Historical daily meteorological data were obtained for the entire sampling period (and preceding months) from The Irish Meteorological Service (Met Éireann, 2020), for the synoptic weather stations nearest to each sampling location (Mount Dillon station for Roscommon and Shannon Airport station for Clare sites). This dataset included daily rainfall and soil moisture deficits (SMD) (calculated for well drained soils, according to the Schulte model (Schulte et al., 2005)), which were used to determine the overall daily effective rainfall (ER). ER was used as a proxy for groundwater recharge given it is a measure of the amount of incident rainfall that has the potential to percolate through the soil matrix and contribute to recharge. Rainfall data were also examined to account for the potential for contaminant transport to a karst spring via shallow surface pathways, and/or rapid transport via solutionally widened conduits and fractures, which is plausible given the nature of the conduit flow karstic aquifers underlying each site. Overall, each of the 11 sites was sampled on a monthly basis for 13 months, from November 2017 to November 2018, with each sample analysed for the 40 anthelmintic drugs. Due to logistical issues with site access, there is a gap in sampling for all seven Roscommon sites and one Clare site for December 2017. In addition, the summer of 2018 was atypical in terms of rainfall, with drought conditions throughout the country recorded for several months. July had just 42.2 mm of rainfall compared to the 30-year mean of 73.1 mm, with soil moisture deficits in the range of 40-90 mm for Clare and Roscommon with no effective rainfall occurring (Met Éireann, 2020). As a result, there are some gaps in sampling for the three sinking streams for the month of July and extending into August in some instances, because of insufficient water for sampling.

Sample collection
Surface and groundwater samples were collected as once-off grab samples in accordance with relevant ISO EN 5667 standards (NSAI, 2006(NSAI, , 2009(NSAI, , 2016(NSAI, , 2018 and represented the water quality at that given moment in time. Surface water samples were collected directly into the sampling bottle, from a minimum depth of 0.5 m from the surface where possible. Groundwater samples were sampled from either a raw water sampling tap at pre-existing distribution pump houses, or directly from the supply (as for several springs). Observation wells were sampled with a discrete depth sampler (closed bailer). In all cases, sample containers were rinsed with the source water prior to sample collection. Samples were received at the laboratory in a chilled condition and stored at <4°C until analysis which was normally carried out within 3 days of collection, and no longer than 7 days. Quality control field blanks were prepared and incorporated into sample collection where appropriate, as described by Mooney et al. (2019).

Anthelmintic analysis
The extraction and instrumental detection of the anthelmintic was carried out according to the method described by Mooney et al. (2019). Samples (500 mL) were modified with methanol (100 mL methanol), adjusted to pH 7, and extracted using a UCT Enviro-Clean highly linked divinylbenzene (HL-DVB) solid phase extraction cartridge. Extracts were analysed for 40 anthelmintic compounds (consisting of 27 parent drugs and 13 transformation products) using ultra-highperformance liquid chromatography tandem mass spectrometry (UHPLC-MS/MS) detection. This comprehensive method allows the analysis of unfiltered raw water sample, thus the measurement of "whole water" concentrations (i.e. analyte fractions in solution and fractions adsorbed to suspended solids). Details of this analytical procedure, including the chemicals, standards, equipment and instrumental parameters, are as described by Mooney et al. (2019). The method was extensively validated and is fit for purpose for the quantitative confirmatory analysis of all anthelmintic residues, besides triclabendazole sulphoxide (TCB-SO) and triclabendazole sulphone (TB-SO 2 ), which are included for screening only. Details of the method validation procedure and results are available open access, as presented by Mooney et al. (2019). The main method detection and quantification limits for all 40 compounds are summarised in Table 1.

. Land-use properties (source factors)
Two different sources of data were used for classification of land-use within the ZOC of the sampling sites, both offering a different level of detail. The first dataset used was the publicly available Corine (Co-ORdinated INformation on the Environment) Land Cover (CLC) datasets, which consists of geo-spatial information on natural and built environments across Ireland (Lyndon and Smith, 2014). Data were obtained (EPA Ireland, 2012) for 44 different land-use classes, and for this study, these were amalgamated into four classes (arable, non-arable pasture, forest and other), which allowed discrimination of the two main agricultural activities in Ireland.
To allow for a more accurate characterisation of agricultural activity with the ZOC of each site, higher resolution data on land-use were obtained from different Department of Agriculture, Food and Marine (DAFM) datasets, namely the Land-Parcel Information System (LPIS), Animal Identification and Movement (AIM) system and the agricultural sheep census (data obtained under license to the Teagasc ACP through a formal data sharing agreement). LPIS is a database maintained by DAFM that is used to identify the location and shape of agricultural land units, known as parcels. The land-use recorded on the LPIS is the crop description (from a pre-defined category) provided by farmers/ landowners annually. Information was obtained for all parcels that fell within (fully or partially) the ZOC of each site, with information on the 95 original crop descriptions recorded for each parcel. For the purpose of statistical analysis, these 95 crop descriptions were amalgamated into four broad groups (grass, tillage, farmyard and other) as described in Supplementary information file 1 Table S2.
The AIM system is a database which contains various levels of details on the identification and movement (herd activity) of cattle, sheep, goats, pigs and horses through a tag number system. Integration of the AIM data with LPIS allows the number of animals to be determined for each land parcel for a calendar year. Given the primary source of anthelmintics in Ireland is due to administration to cattle and sheep, data were obtained for these species only, with data broken down depending on sex or age. Using these data, the number of livestock units (LU) were calculated by application of different LU coefficients as described in Annex I of Regulation 1200/2009/EC (European Parliament, 2009). LU values for individual parcels were then adjusted to account for the percentage of the parcel that fell within the ZOC of each site, with the overall total LU within the ZOC calculated by summing all adjusted values for the given ZOC. Cattle and sheep stocking density was determined by expressing the number of animals per unit area (hectare, ha), accounting for the entire area of each ZOC. Finally, the animal stocking data were also used to calculate the nitrogen load due to animal excretion, expressed as nitrogen per hectare (N ha −1 ). Nitrogen load was calculated for each parcel by applying the different excretion rates for each animal, as described in Table 6 of the GAP (Good Agricultural Practice) Regulations (Government of Ireland, 2017), and then related to the agricultural area within each ZOC.

Physical and hydrogeological properties (pathway factors)
For statistical analysis, groundwater sampling points were classified as either springs or BHs, with BHs further subclassified as abstraction boreholes or monitoring BHs. BHs pumped daily and used as drinking water supplies or for agricultural purposes, were classified as abstraction BHs, while BHs bailed or pumped for short periods less frequently (e.g. weekly or monthly) and used for research and observational purposes, were classified as monitoring BHs. The hydrogeological (pathway) properties used for sites characterisation are summarised in Supplementary file 1 Table S3, and included: (a) bedrock geology presented as 27 hydrostratigraphic units (GSI, 2016a) amalgamated into six lithological groups as described by Tedd et al. (2017), (b) GSI aquifer category comprising 11 aquifer classes (see Supplementary file 1 for a description of the 11 classes) (GSI, 2015b,c), (c) Water Framework Directive (WFD) flow regime (karstic, productive fractured and poorly productive fractured) as described by the Working Group on Groundwater (2001), (d) flow regime further categorised into three classes (conduit-dominated karstic, diffuse-dominated karstic and fractured) as previously described for site selection in Section 2.1.1, (e) groundwater vulnerability consisting of five vulnerability classes (X-Extreme, E-Extreme, H-High, M-Moderate and L-Low) (GSI, 2015d), (f) Irish Forests Soils (IFS) (Bulfin et al., 2002;Teagasc-EPA-GSI, 2006) consisting of 25 soil classes simplified and dichotomised based on the four main principal components (mineral/peat, acidic/ basic, well drained/poorly drained, shallow/deep) used in the classification system described by Fealy et al. (2009), (g) Quaternary sediments, presented as 53 sediment classes (GSI, 2016b) amalgamated into seven geneses and (h) subsoil permeability, extracted from the GSI Subsoil Permeability dataset (GSI, 2015a), and presented as four permeability classes (high, moderate, low and depth to bedrock less than 3 m (DTB < 3 m) where permeability was not determined).
Data on the relative percentage of each property class was extracted using the ARCGIS intersect tabulate tool, and the class occupying the largest percentage of the ZOC was assigned as the predominant class, as summarised in Supplementary file 2 for all groundwater sites.

Statistical analysis of relationships with site characteristics
The non-random selection of the data sets for analysis introduces a risk of bias in the outcomes but the authors believe that there is sufficient value in the data to warrant an analysis. The interpretation should be used to guide future research rather than as the basis for definitive conclusions.
The association between detections and the recorded characteristics (land-use and hydrogeological properties) of the groundwater sites was analysed using SAS 9.4 (SAS Institute Inc., 2014). The main test to be used was logistic regression with compound concentration as the response and site characteristics as the regressor(s). Due to the high degree of censoring (compound concentration below the limit of detection (LOD)), regression based on the Normal distribution was not feasible. A way forward was to recode the data in a binary format so that any samples with a detection greater than the LOD were recorded as a detect, with results <LOD recorded as a non-detect. The Logistic procedure was used to fit regression models with detection/ non-detection as a binary response. Marginal tests were conducted for each characteristic and a variable selection procedure was used to build a conditional model. The former established the usefulness of each variable on its own as a predictor of detections while the latter checked for any combination of the characteristics that was useful for such prediction. The main focus in evaluating the results was on the interpretation of odds ratios and their 95% Wald confidence intervals. Characteristics with continuous variables were also analysed using marginal tests with detection/non-detection as the regressors, with the means of sites with detections compared to those recorded as nondetection. Significance was assigned as p values ≤ 0.05, with values slightly higher than 0.05 noted as being of interest. Quasi-complete separation, a failure of the regression procedure for binary data, was encountered frequently in the analyses and therefore the Firth option in Logistic procedure, which adds a penalty term to the log likelihood, was used to achieve a solution to the equations. The overall analysis was complete for detections defined for several different groups of compounds as follows; (a) all anthelmintics, (b) BZs, (c) non-BZs, (d) parent compounds and (e) TPs.

Results
3.1. Analysis of anthelmintic spatial occurrence 3.1.1. Anthelmintic spatial occurrence summary The geographical spread of the 106 sampling sites, classified by sampling point type (surface water vs. groundwater) and detection vs. non-detection, is shown in Fig. 1. One or more anthelmintic drugs were detected at 20% of the total sites (23 of 106 sites) at levels greater than the LOQ, while a further 2% of sites had detections at levels less than the LOQ but greater than the LOD (i.e. present but not quantified). In total, 17 out of the 40 different compounds were detected throughout the campaign, with eight different anthelmintic compounds detected across 39% of SW sites at levels >LOQ (1.0-40.6 ng L −1 ) ( Table 2) and up to five compounds detected at any one SW site (Table A1). Seventeen compounds were detected across 16% of GW sites at levels >LOQ (0.9-22.2 ng L −1 ), with a further 2% of GW sites with detections <LOQ but >LOD (Table 2). Up to six different anthelmintics were detected at any one GW site (Table A1). The concentration range of each compound detected in SW and GW is also summarised in Table 2. Detection of an anthelmintic was shown to be 2.9 times more likely (p = 0.057) at SW sites than at GW sites.
ABZ and its sulphoxide (ABZ-SO), sulphone (ABZ-SO 2 ) and aminosulphone (ABZ-NH 2 -SO 2 ) transformation products (TPs) were the most detected compounds, each detected in 3.8%, 6.6%, 3.8% and 7.5% of total sites, respectively. ABZ and/or its TPs were detected at 8% of GW and 28% of SW sites. ABZ was also the compound detected at the highest concentration during the campaign, detected at 40.6 ng L −1 in a field drain, with all three transformation products also detected in this same sample (Site 089, see Appendix A Table A1). Converting the concentration of each of these TPs back to the parent equivalent, the combined concentration of ABZ in this sample was 73 ng L −1 . The ABZ-SO 2 TP was the compound detected at the highest concentration in GW (22.2 ng L −1 ) (Site 079). The next most detected compounds in SW were TCB (16.7% of SW) and OXYCLOZ (11.1% of SW), while OXYCLOZ was the second most detected in GW (4.5%), followed by FBZ, OXF, LEV and CLOS, all of which were detected at 3.4% of GW sites. Notably, TPs were more predominant in groundwater compared to surface waters, with twenty-two detections of eight different TPs in eight groundwater samples, compared to seven detections of three different TPs in five surface water samples.

Occurrence and ZOC land-use characteristics
Analysing categorical and continuous variables with detection/ non-detection as a binary response showed that the detections of (a) all anthelmintics, (b) BZs, (c) non-BZs, (d) parent drugs and (e) TPs, were significantly associated (p < 0.05) with one or more of the land-use properties used for site characterisation. The land-use characteristics shown to be associated with each detect grouping are summarised in Table 3, along with a description of odds ratios and their interpretation. Similar results were observed for the analysis of LPIS land-use and for the analysis of the Corine dataset, and this is believed to be due to aspects of the Corine data being derived from LPIS. To avoid duplication, results are therefore only presented for LPIS land-use in Table 3.
Detections of all anthelmintics, collectively, were related to sites with higher proportions of agricultural land in their ZOCs, in addition to higher sheep density (specifically ewe density and "other" nonbreeding sheep density). This same trend was observed for detection defined for parent compounds only and non-BZ compounds. Ram density also produced results with p < 0.05, however interpretation of these results was not practical due to the outcome being dependent on a limited number of sites only, producing a questionable model fit. Occurrence of BZ compounds was statistically related to higher proportions of agricultural land, a higher proportion of tillage crop and a lower proportion of grass crop, within the ZOCs. A similar trend was observed for TPs, with detections associated with a higher proportion of tillage and a lower proportion of grassland in the ZOC. However, the overall proportion of agricultural land within the ZOC was not statistically related to TP occurrence, neither were the LPIS land-use or Corine land cover classes. All other land-use characteristics not listed in Table 3 were not statistically related to anthelmintic occurrence, with a complete summary of p-values shown in Supplementary file 1 Table S5.
For several characteristics, analysis of the continuous variables using detect/non-detect as the regressors indicated a significant difference between the mean of sites that had a detection compared to the mean of sites with a non-detection. The mean percentage of grassland in the ZOC of sites with a BZ detection (52.4%) was significantly different from the mean of sites with no detection of BZs (74.9%) (p = 0.0356). The mean percentage of tillage within the ZOCs of sites with a BZ detection (46.5%) was also significantly different from the mean of sites with no BZ detections (9.6%) (p = 0.0102). Like BZs, the mean percentage of tillage within the ZOC of sites with detections of TPs (38.1%) was also significantly different from the mean percentage of tillage at sites with no detection (10.8%) (p = 0.0345). The mean density of ewes at sites Table 2 Summary statistics for the anthelmintic compounds detected above the limit of detection (LOD) and limit of quantification (LOQ), at surface water (SW) and groundwater (GW) sites. with a detection of non-BZ compounds (1.44 ewes per Ha) was significantly different to the mean density of ewes at sites with non-detection of non-BZ compounds (0.25 ewes per ha) (p = 0.008). Similarly, the mean density of "other sheep" at sites with a non-BZ detections (0.2 ewes per ha) was significantly different from the mean of sites with non-detection (0.11 ewes per ha) (p = 0.03).

Occurrence and ZOC hydrogeological characteristics
Analysis of the different groupings of anthelmintic detections indicated that occurrence of one or more groups was significantly related to several pathway factors including sampling point type, aquifer category, bedrock group, Quaternary sediments and IFS soil types. Detections of all anthelmintics collectively, were shown to be more likely in (a) monitoring BHs compared to abstraction BHs, or springs, (b) noncalcareous bedrock compared to pure limestone bedrock, and (c) poorly drained mineral soils (minPD), as opposed to deep well drained soils (minDW) ( Table 4). While statistically significant, it must be noted that this relationship with bedrock is biased as a result of several detections in a cluster of sampling points within the same catchment.
Assessing parent compounds only, detections were shown to be more likely in non-calcareous bedrock as opposed to pure limestone, however the same issue remains with these results, due to the localised cluster of detections in one catchment. While not fully significant in terms of p-value, IFS soil type odds ratios analysis indicated a higher likelihood of a detection of a parent compound in poorly drained soils, as opposed to deep well drained soils. Non-BZ compounds (i.e. CLOS, LEV, OXYCLOZ and RAFOX, all of which are parent compounds) show similar relationships with bedrock and IFS soil type. Aquifer class was also shown to be significant, with detections more associated with poorly productive Pl aquifers. This relationship is also questionable, due to the same localised cluster of detections.
TPs on the other-hand were shown to be statistically associated with locally important (Ll) aquifers, with a detection 7.7 times more likely in (Ll) aquifers compared to regionally important conduit dominated (Rkc) aquifers, and 1.021 times more likely, for every 1% increase in the percent of Ll aquifer within the ZOC. TPs were more associated with impure limestone bedrock than pure limestone bedrock, and with Quaternary sediments comprised of tills derived from Namurian sandstone and shales (TNSSs) as opposed to tills derived from limestones (TLs). Similar associations (to TPs) for aquifer category, bedrock group and Quaternary sediments were shown for BZ compounds, which is not unexpected given that all detected TP compounds also belong to the BZ class. However, additional pathway factors were shown to be related to BZ occurrences: detections were more likely in monitoring BHs than in springs and more likely in shallower soils than deeper soils (based on IFS soil data).

Analysis of anthelmintic temporal variations at catchment level
There were 11 anthelmintic compounds detected across the 11 sites at some point throughout the temporal study. These compounds included ABZ and its three TPs (SO, SO 2 and NH 2 SO 2 ), FBZ and its TP (OXF), TCB, LEV, CLOS, OXY and IVER. The anthelmintics detected at each of the 11 sites across the 13 months, with their respective concentrations, are shown in Tables S7 and S8 (Supplementary file 1). The frequency of detection of anthelmintic compounds, in addition to the average number of anthelmintic detections per site (DPS), is shown in Fig. 3(a), while Fig. 3(b) depicts the mean concentration of anthelmintics detected at both Clare and Roscommon sites, on each sampling occasion. Considering these data, there is a trend evident at both locations, with two spikes in anthelmintic detections, firstly during the months of February and March, followed by a gradual decline, and a second spike during the month of August. These spikes occur not only for the frequency of detection and mean DPS, but also the mean concentration of anthelmintics. During February and March, an average of 1.6 anthelmintic compounds were detected at 86% of Roscommon sites, with an average of 1.5 compounds detected at 75% of Clare sites. The mean concentration detected in February was 10.3 ng L −1 at Roscommon sites and 14.6 ng L −1 at Clare sites, which further increased to 23.7 and 18.9 ng L −1 in March, respectively. In August, the spike in detection is Table 3 Summary of land-use (source) characteristics showing a significant relationship with the occurrence of (a) all anthelmintics, (b) benzimidazoles, (c) non-benzimidazoles, (d) parent drugs and (e) transformation products (TPs), with corresponding p-values, confidence intervals and odds ratio likelihood interpretations. most evident in terms of the mean concentration, with the frequency and mean DPS slightly lower in comparison to the Feb/March period. Comparison of anthelmintic detections across all 11 sites using statistical analysis indicated a significant association with month (p = 0.0021). Detections of anthelmintic compounds were 26.6 times more likely in February compared to January (95% intervals: 2.66-266.5), 26.3 times more likely in March compared to January (95% intervals: 2.65-250), 10.1 times more likely in August compared to January (95% intervals: 1.14-90.1) and 11.6 times more likely in September compared to January (95% intervals: 1.34-100). Furthermore, detections were 9.9 and 11.4 times more likely in both February and March compared to July (95% intervals: 1.18-82.63) and November respectively (95% intervals: 1.96-66.18). These initial observed trends appear to coincide with periods of expected anthelmintics application as part of seasonal events for the different production systems e.g. animal turn-out to pastures and spread of manure around February/March, as well as summer dosing of spring lambs and calves. These trends are discussed in more detail in Section 4.3.
In order to further assess and explain these trends, daily rainfall and the calculated daily ER were examined for the period during and prior to each sampling event. Fig. 4 presents the concentration of total anthelmintics detected at each of the seven Roscommon sites, with the data overlaid by the daily rainfall (black bars) and daily effective rainfall (blue bars) (bottom chart). The breakdown of these concentrations for parent compounds compared to transformation products is also shown in Fig. 4 (top chart). Similar data are presented for the Clare sites in Fig. 5. It is evident from these graphs that the spikes in anthelmintics detections and concentrations in Feb/March, and August, also coincide with period of increased rainfall and/or ER, during the days immediately prior to the sampling date. The importance of rainfall is more evident for Clare sites, which show no apparent ER occurring prior to the spike in detections in August. 4. Discussion

Anthelmintic spatial occurrence
The timing of the spatial investigation (March/April 2017) coincided with a period which falls during, or just after, two of the main seasonal dosing occasions: animal housing and turnout to pasture (Supplementary Fig. S1). A general understanding of these usage patterns is therefore important in interpreting the anthelmintic occurrences. Housing represents a period of indoor confinement of the animals following the transition from pasture, typically during winter months (e.g. September/October onwards for cattle (Bloemhoff et al., 2014)), or at defined stages of the production calendar (e.g. drying-off of dairy cows or during lambing). Anthelmintic dosing during housing is important to treat worms and flukes acquired by the animal on pastures, prior to housing. The housing period also typically occurs during what is known as the "closed-period" (typically mid-October to mid-late January) whereby spreading of manure, slurry and chemical fertilisers is prohibited under the GAP (Good Agricultural Practice) Regulations (Government of Ireland, 2017). During housing, contaminated animal waste is stored for land-spreading during the growing season. Therefore, while anthelmintic drugs administered at housing are not likely to be at immediate risk of entry into the environment, there is the potential of their entry from mid-January onwards, as a result of the re-commencement of land-spreading. According to a National Farm Survey Report for 2017 (Buckley et al., 2020), 46% of farmyard waste produced during housing was spread during the period of January to April, thus contributing a potential source of anthelmintics. The turnout stage represents the period where housed animals are returned to grazing. Depending on ground conditions, this can occur as early as February for cattle (Bloemhoff et al., 2014) and from mid-March through to May for sheep (Patten et al., 2011). Treatment of cattle at turn-out focuses on treating roundworms in calves, with frequent dosing (every 3-6 weeks) until the end of July, while treatment of sheep focuses on preventing infection of lambs by Nematodirus spp., with the highest risk during April and May (Animal Health Ireland, 2020a). During these periods, administered anthelmintics can also be excreted directly onto the soil while grazing. Overall, it is certainly plausible that the detections during March and April relate to such usage patterns; however, this is not definitive, given the potential of varying time lags for the transport of contaminants to water bodies.
The BZs were the most detected class of anthelmintics, with 14 of 17 compounds detected belonging to this class. The prevalence of BZ detections is not surprising given they are one of the most commonly used classes of anthelmintics, due to their broad-spectrum efficacy toward nematodes at all stages of their life-cycle (i.e. immature and mature worms). Keegan et al. (2017) reported that BZ drugs were the most popular anthelmintic class used by sheep farmers in Ireland, with usage in 42% of the cases. ABZ is the most common BZ administered at housing (typically at the latter stage), owing to its dual efficacy toward both nematodes and mature liver flukes (Bloemhoff et al., 2014). This could account for the higher frequency of detection of ABZ and/or its TPs in this study (occurring in 52% of sites with detections). TCB and its TPs (TCB-SO and TCB-SO 2 ) were the second most detected BZs, the frequency of which (22% of sites with a detection) was much lower than the ABZ compounds. TCB is the only drug effective toward earlyimmature flukes; however due to resistance issues, it is only used for treating acute infections of liver fluke (Animal Health Ireland, 2013Ireland, , 2020b, and this is likely to account for lower frequency of detection. All other BZs that were detected are authorised for used in both cattle and sheep; however, their usage trends are less predictable and depend on factors such as availability, cost and the farmer's preference and past experience (Patten et al., 2011). FBZ and OXF are sometimes favoured for treating nematodes in lactating cows, due to the short milk withholding times, which could account for their detection at 19% of sites with detections. In order to further assess these trends in BZ detections, it is also important to consider other factors, such as the environmental mobility, that may influence the environmental occurrence. A comprehensive review of such factors was carried out by Horvat et al. (2012), who provide a critical overview of the current knowledge on the fate and the ecotoxicology of anthelmintics and their TPs. Considering their physiochemical properties (Table 1), on application to soil, it is expected that BZs will be relatively immobile and more associated with soil and sediment compartments than aqueous phase. Such effect is demonstrated by Kreuzig et al. (2007) who investigated the fate of FLU and FBZ in manure and manure amended soil. Both BZs demonstrated slow degradation in manure, and on application to soil, both were extractable from the near surface soil. Notably, this study also demonstrated the degradation of FBZ to its metabolites FBZ-SO (OXF) and FBZ-SO 2 in the manure amended soils, with OXF increasing to 45% within 14 days, after which is was further converted to the SO 2 derivative (12% of the amount applied). This is consistent with the findings of our study, with the FBZ parent and OXF TP detected more frequently than the FBZ-SO 2 ( Table 2).
Considering the pKa of ABZ, and other BZs, and the typical pH of soils in Ireland (typically between 5.5 and 6.5 (Teagasc, 2018)), there is the potential for the BZs to be ionized to a certain extent that may increase their mobility. However, Mutavdzic Pavlovic et al. (2018) indicated the medium to strong sorption of ABZ on sediment and soil and noted that the sorption of ABZ is particularly dependent on pH, with the neutral form of ABZ predominant between both pKa's. Therefore, this neutral form may restrict the mobility of ABZ (at this typical soil pH range), which could account for the lower occurrences of ABZ in groundwaters compared to its SO, SO 2 , and NH 2 SO 2 TPs. These TPs are more polar and likely to have greater water solubility, which could enhance their transfer to surface and groundwater (Horvat et al., 2012). The prevalence of BZ TPs in this study can also be accounted for by degradation of the parent drugs during transport through environmental compartments. As highlighted by Horvat et al. (2012), BZs including ABZ, FBZ, FLU, MBZ and TCB, can potentially undergo extensive transformation (by oxidation, reduction or hydrolysis) while in the environment, with photolysis indicated as one of the main degradation pathway, particularly for ABZ. However, if transport is facilitated by underground pathways, the contaminants will be less exposed and may persist for longer.
Besides the BZ compounds, the other detected anthelmintics are commonly used as flukicides (CLOS, OXYCLOZ and RAFOX), except for LEV. Spring calves are typically dosed with a flukicide, such as CLOS and OXYCLOZ, initially at turn-out (if not during housing) and again during late summer and early autumn when infection increases with the onset of wetter conditions. Sheep are much more susceptible to fluke than cattle since they graze closer to the ground, with treatments involving flukicides such as CLOS or RAFOX used for spring dosing to treat chronic and subacute fluke infection. Such usage can therefore account for the occurrence of these compounds. The environmental fate and degradation pathways of these compounds are not widely reported; however, CLOS and RAFOX are highly hydrophobic and may also be subject to similar transport mechanisms as the BZs. LEV is a drug used as a wormer, and is often used in combination with TCB, with products containing LEV often favoured later in the season, if BZ efficacy fails. While LEV has been shown to have poor mobility in the environment due to strong sorption to soil (Ma et al., 2019), its high solubility and relatively low LogK ow indicates it is susceptible to preferential and conduit flow pathways in which the bulk soil matrix is bypassed. The lack of detections of the ML class (endectocides) of anthelmintics is notable, and somewhat surprising, given that these drugs are used just as heavily as the BZs. IVER can often be favoured at housing due to its efficacy both as a wormer, but also as an ectoparasiticide for external parasites such as ticks and lice (Animal Health Ireland, 2013). The lack of groundwater occurrence of IVER may be explained by its relative instability (photodegradation DT 50 of 6-39 h in water) (Prasse et al., 2009) and lack of mobility in the environment (Krogh et al., 2008b). It must be noted however, given the typical concentrations of BZs detected (majority <10 ng L −1 ), the lack of detection of the MLs relative to the BZs may be of a consequence of the lower sensitivity of the analytical method for MLs, with method limits for the BZs 10 times lower than for the MLs.
To the best of our knowledge, there have been no reported occurrences of most anthelmintics in groundwaters, except for thiabendazole detected in a single farm well in Norway (Haarstad and Ludvigsen, 2007). The overall findings of this current work however are consistent with other previously reported environmental occurrences in surface waters, both in terms of the compounds detected, and the concentrations detected as part of the spatial study (1-41 ng L −1 ). Van De Steene et al. (2010) detected FLU in surface waters at concentrations ranging from 0.3-20.2 ng L −1 , while Wagil et al. (2015b) reported the detection of FLU and FBZ in surface waters in Poland in the range of 5.4-87.5 ng L −1 . Prior to this current study, the most comprehensive investigation of anthelmintics in environmental waters was carried out by Zrncic et al. (2014) who reported the detection of eight anthelmintics, including ABZ, FBZ, FLU, LEV, TCB and MOXI, typically at concentrations ranging from 1 to 5 ng L −1 , except for LEV (up to 39 ng L −1 ). One noticeable difference, however, is the prevalence of TP detections in the current study, which most other studies did not include in their analytical methodology. Detection of additional compounds in the current work, such as the flukicides OXYCLOZ and CLOS, may be more reflective of usage patterns, with the wetter, damper climate in Ireland having a higher fluke burden, thus need for flukicides, compared to Mediterranean countries such as Spain.

Relationship of anthelmintic groundwater spatial occurrence with site characteristics
Statistical analysis of the different source factors indicated that occurrence of non-BZ compounds was significantly related to the density of sheep, with detections associated with higher densities of ewes and "other sheep". This association can be accounted for by a combination of both source and pathway factors. The two frequently detected non-BZ compounds (OXYCLOZ and CLOS) at the sites with higher sheep density were not surprising given both are used for treating liver flukes, to which sheep are particularly susceptible. Geography may also have an effect, with most of the detections of these flukicides occurring in the west of Ireland, which has an expected higher fluke burden as a result of wetter climate and poorer ground conditions (Parr and Gray, 2000). Bloemhoff et al. (2014) reported a higher usage of flukicides on dairy farms in western regions and attributed this to the higher fluke burden, as well as the higher proportion of mixed cattle and sheep farms. Western regions of Ireland are also where the majority of sheep production occurs (DAFM, 2017). Therefore, it is likely that this relationship with sheep intensity is localised to such regions and may be dependent on environmental conditions. Statistical analysis also showed that occurrence of BZs and their TPs was more likely at sites with ZOCs dominated by tillage land compared to grassland. This may be initially explained as a result of the method of application of animal waste to both systems. A national study by Hennessy et al. (2011) reports that the majority of manure applied to tillage was solid manure incorporated by ploughing (compared to liquid manure on grassland by spray plate), which may potentially facilitate the direct entry of the contaminants to lower depths below the surface, where they may be less attenuated due to the lower soil organic carbon and microbial activity (Alletto et al., 2010). However, the overall interpretation of the relationship of detections with tillage is difficult given that the sites dominated by tillage land represented only a small sample size (n = 7) compared to grassland (n = 81), with the overall relationship based on BZ detections at just 3 of 7 tillage dominated sites, two of which were in the same surface catchment. Therefore, while significant, this relationship cannot be regarded as conclusive. It is also notable that although there was a negative association between detections and the percentage of grassland within the ZOC, there was still an appreciable proportion of grassland within the ZOCs of sites that had a detection (mean of 52%) compared to sites with non-detections (mean of 75%). Therefore, a much more focused and detailed study is required to properly disentangle these relationships.
Several pathway factors (including sampling point type, aquifer type, bedrock group, Quaternary sediment and IFS soil type) were shown to be statistically related to the occurrence of anthelmintic compounds (of one or more classes of these). However, meaningful interpretation of the relationship with a number of these hydrogeological factors (namely aquifer class, bedrock group and Quaternary sediment) is complicated as a result of clusters of sites at catchment level and/or a relatively low number of detections from a small sample size, causing a bias in the statistical analysis. It must therefore be stressed that these observations should only be considered in an exploratory manner and may be used for informing future studies.
The occurrence of anthelmintics (collectively) was shown to be more associated with monitoring BHs than abstraction BHs, and this may be reflective of the larger ZOCs associated with abstraction BHs, compared to the monitoring BHs (Median ZOC area of monitoring BHs was 0.09 km 2 compared to 0.50 km 2 for ZOCs of abstraction BHs). Monitoring BHs, having smaller ZOCs, are likely to be more sensitive to contamination by nearby localised activities and sources of anthelmintics. For abstraction BHs, a localised area of contamination within a larger ZOC is more likely to undergo dilution due to the larger volume of aquifer from which the water is drawn, with longer distances and travel times to the BH providing more opportunity for attenuation of the contaminants. Detections were also shown to be more likely in poorly drained soils (minPD) compared to deep well drained soils (minDW). This relationship appears to be driven more by the non-BZ and parent drug subgroups of detections and there is evidence that this association may be confounded by the relationship with sheep density. Transport of contaminants to the groundwater below the poorly drained soils may be facilitated by preferential flow through the soil via macropores, as demonstrated by Weiss et al. (2008) who report loss rates of up to 16% of BZs via preferential flow beneath permanent pastures. Transport may also occur via surface or near surface pathways (particularly with heavy rainfall) (Stamm et al., 2002;Kreuzig et al., 2007) to features such as karstic sinking streams, which can then provide further rapid groundwater transport via conduits. The latter of these seems a potentially plausible explanation, with 3 of the 5 minPD sites with detections being karstic springs which are known to be fed by conduit flow with rapid travel times. The influence of surface pathways in poorly drained soils is also reflected by the higher frequency of anthelmintic detections in the surface waters: almost half of the detections in surface streams occurred in catchments dominated by minPD soils.
As a consequence of small sample numbers, the observed relationships between anthelmintic detections and aquifer class (Pl and Ll), are likely being influenced more by the localised characteristics of these sites. For example, of the five sites with BZ detections that had ZOCs dominated by Ll aquifers, two sites (site 032 and 076) had an appreciable proportion of karstic aquifer within the ZOC (43% Rkd aquifer for site 076 and 42% Rkc aquifer for site 032), with these karstic aquifers supplying the borehole (076) and spring (032). In addition, three of the five sampling points were classified as monitoring BHs. It is therefore believed that the occurrence of anthelmintics at these sites is more likely related to a combination of such localised factors. Finally, although the relationship of anthelmintic detections with groundwater vulnerability was not statistically significant (Table S6), there is some evidence of a weak relationship with extreme (E) vulnerability. Whereas 34% of the 88 groundwater sites had ZOCs dominated by E vulnerability, 56% of the 16 groundwater sites with detections were dominated by E vulnerability. Further investigation of the relationship between detections and the percentage of the ZOC with E vulnerability showed only a weak relationship (p = 0.13), but it is possible that more targeted studies with fewer confounding factors could determine a clearer relationship with groundwater vulnerability.

Interpretation of temporal trends
Two of the most important factors likely to be influencing the temporal variation in occurrence of these drugs at the conduit dominated karstic sites are (a) the usage patterns and routes of entry of the drugs (source factor), and (b) meteorological events, primarily rainfall events controlling recharge (pathway factor). Considering the same usage patterns discussed for the spatial study (Section 4.1 and Supplementary  Fig. S1), and assessing the temporal trends observed as part of this study, it is evident that the spikes in detections during Feb/March and August/September of 2018 coincide with an increased source of the drugs entering the environment as a result of periods of heavier usage. The increased detections in Feb/March 2018 may be explained by similar usage patterns discussed for the spatial study; for 2018, 42% (national scale) of bovine slurry gathered during housing was landspread between January (mid-late) and April (Buckley et al., 2020). Both study catchments are also dominated by grassland, onto which 97% of bovine slurry is spread (based on national average (Hennessy et al., 2011)), so it is evident that the first spreading of slurry on recommencement after the closed period is also potentially a significant source of anthelmintics. This combined with the usage at turn-out provides reasonable explanation for the initial spike in detections in February/March. Further application (of the remaining 41% of stored slurry) through the summer months, in addition to the regular dosing of spring lambs and calves in pastures (through to the Autumn), is likely contributing to the spike in detections in August. The lower frequency of detections throughout the months preceding August, is interpreted as a result of meteorological factors.
As with the spatial study, ABZ and its TPs were commonly detected at several sites throughout the temporal study. When administered at housing, ABZ will be excreted primarily in urine (Wardhaugh, 2005), after which it will be land spread as slurry. This, along with direct excretion on pasture, could account for the initial early Feb/March detections of the ABZ parent drug, with the increasing concentrations of the ABZ TPs in the months thereafter likely due to a combination of both degradation and a potential increased source from animal excreta. The main sulphoxide metabolite that is excreted by animals (Junquera, 2015) can undergo further transformations to the amino-sulphone derivative, which is also the main environmental photodegradation product of ABZ (Horvat et al., 2012). The occurrence of OXYCLOZ is also probably attributable to land spreading of the drug administered during housing, given that it is only effective for treating mature liver fluke (Bloemhoff et al., 2014). LEV was typically only detected from the month of May onward, and this can be explained by the fact that LEV is used for treating Teladorsagia spp. and Trichostrongylus spp. roundworm, which pose the highest risk of infection from June to August. It is also notable that IVER was detected on a number of sampling occasions at one site (Clare D), with the highest concentrations detected during February and March 2018. IVER was not previously detected as part of the spatial study in March/April 2017.
Having established plausible sources of anthelmintics in the environment, assessing the meteorological conditions prior to the sampling events can provide further explanation for the trend in detections. Since the sites sampled as part of this temporal study have the potential for rapid transport of contaminants from sinking streams to springs in as little as 8 h, meteorological events during the hours and days just prior to the sampling events are likely to have the most influence on the transport of the contaminants. It is important to note however that this is not always the case, with the potential of longer lag times (from months to years) in other hydrogeological settings such as low permeability soils and less transmissive aquifers (Fenton et al., 2011), so recent meteorological events may be of less importance in such situations.
For the Roscommon sites, there is a clear trend of significant rainfall events occurring during the periods just prior to the February and March spike (Fig. 4), with most of the incident rainfall was recorded as effective rainfall, thus contributing to recharge. The mean daily rainfall and effective rainfall for the week prior to the sampling event was 2.8 and 2.3 mm per day, with as much as 12 mm of effective rainfall occurring on a given day during this period, which is likely sufficient to provide for the transport of the veterinary drugs to groundwater. It is also notable that January 2018 was particularly wet, with an average monthly rainfall of 177 mm (compared to a 30-year average of 92 mm). Consequently, although permitted by legislation, there is the potential that manure spreading was not carried out immediately after the end of the closed period, rather later into February, providing a potential increased source of the drug residues closer to the sampling event, with contaminants carried to groundwater with any subsequent rainfall. There were very few ER events occurring in Roscommon from April through July, with accumulated soil moisture deficits and higher evapotranspiration rates consuming any incident rainfall. Although there was a potential regular source of drugs throughout these summer months, these dry conditions may account for the decrease in frequency and concentration of the detected anthelmintics, due to the contaminants remaining on the surface or in the soil matrix. The spike in detections in August occurred after 11 mm ER within 2 days before sampling, which potentially resulted in flushing of a concentrate of the contaminants from the soil or sediment along the flow path.
Similar ER trends were observed for Clare sites (Fig. 5), to those described for Roscommon, with one main difference; the detections at Clare sites in August do not appear to coincide with any ER event, with large SMDs and no ER occurring from mid-April through to early September. These sites provide evidence of the importance of surface pathways in such karst hydrogeological settings. The meteorology data show up to 20 mm of rainfall occurring in Clare during the period prior to the sampling event. Although not contributing to ER due to the accumulated SMD, it is plausible that this high incidence of rainfall resulted in rapid run-off and/or direct infiltration to karstic aquifers.
Besides these usage and meteorological patterns, some interesting trends were observed for two of the Roscommon karst springs (C and E), with the same compounds detected on a number of sampling events, both at the spring and at a sinking stream known to be connected (based on previous dye tracing (GSI, 2014)). At Roscommon site C (sinking stream), the parent compound ABZ was detected initially during February and March, with an increasing concentration of TPs (particularly the amino sulphone TP) detected for several of the months thereafter (data provided in Table S7). Site observations indicated evidence of cattle grazing and sheep farming in the immediate vicinity of the sinking stream. At the associated spring (Roscommon D), the same ABZ TPs were detected during a number of sampling events (February, March, April, June and Sept), albeit at lower concentrations, which provides evidence of a potential poorly attenuated surface to groundwater pathway for anthelmintics. The lower concentrations combined with the lack of detection of parent drug and the prevalence of the amino-sulphone at the spring (compared to the sinking stream), does however indicates some degree of attenuation and/or degradation during transport. Another example is provided by Roscommon sites E and F: OXYCLOZ concentrations detected at the sinking stream (F) (10.8 and 11.7 ng L −1 ) and the spring (E) (11.8 and 13.3 ng L −1 ) in February and March, were quite comparable. While there is very little information reported on the fate of OXCLOZ in the environment, its relatively high LogK ow and LogK oc values (Table 1) indicate that it is likely to persist in the environment within the soil compartments. The occurrences at Roscommon spring E further indicate the potential for conduit fed systems to accommodate the unattenuated transport of some anthelmintics, which would otherwise persist and remain in the soils. While the suggested transport in both examples is not definitive, given that the groundwater arriving at a karstic spring is supplied by a network of conduits and fractures rather than a single sinking stream, it does clearly demonstrate the vulnerability of these systems and the need for further site-specific investigations.

Ecological significance
While some of the anthelmintic compounds are covered as "nematocides" under the pesticide definition specified in the EU Drinking Water Directive (European Commission, 1998) and Groundwater Directive (European Parliament, 2006), there is currently no regulatory monitoring of these contaminants in groundwater or drinking water. Furthermore, there are no environmental quality standards (EQS) set for these contaminants in environmental waters. Considering the concentrations of anthelmintics detected in the current study, and accounting for the TP concentrations, although there were no breaches of the pesticide parametric value of 100 ng L −1 (0.1 μg L −1 ), there were a number of occasions when the concentrations of anthelmintics (e.g. albendazoles) met or approached the pesticide threshold value of 75 ng L −1 (Government of Ireland, 2014), indicating that such contaminants may require additional investigation. Such investigations should also prioritise the assessment of the environmental impacts of anthelmintics. The majority of the benzimidazoles concentrations detected are not likely to pose an acute environmental risk given the reported ecotoxicological data have shown half maximum effective concentrations (EC 50 ) of the order of low μg L −1 for aquatic organisms such as Daphnia magna. However, there is a concern for the levels of the ML IVER, which was detected at concentrations up to 47.5 ng L −1 in the temporal study. Previous studies have shown that IVER is highly toxic to a number of aquatic organisms, including the Neomysis integer, Gammarus sp. and Daphnia magna, with reported lethal concentration (50%) (LC 50 ) values ranging from 25 to 70 ng L −1 (Hally et al., 1989;Horvat et al., 2012). Even higher acute toxicity to Daphnia magna has been reported by Garric et al. (2007) with even lower LC 50 values of 5.7 ng L −1 . The same authors report extremely high chronic toxicity to IVER at levels of 0.3-1 pg L −1 . Further research is needed to determine environmental quality standards for these emerging contaminants in water.

Conclusions
This study, the first of its kind in Ireland, reports on the most comprehensive investigation of an extensive suite of anthelmintic drugs (including 27 parent drugs and 13 transformation products) in groundwaters and a limited number of surface waters throughout Ireland. Up to sixteen different anthelmintic residues were detected across 18% of the groundwater sites, with up to eight compounds detected in 39% of the surface waters sampled. The overall detected concentrations of anthelmintics ranged from 1 to 41 ng L −1 . Of the 17 compounds detected, 13 belonged to the benzimidazole class of anthelmintics, which are commonly used in Ireland. Of these, albendazole and its three transformation products (albendazolesulphoxide, albendazole-sulphone and albendazole-amino-sulphone) were the most commonly detected, found at 52% (12 of 23) of the sites with detections, and at 11% of the total (106) sites, with occurrence in both groundwater and surface waters. Fenbendazole and its two metabolites were also frequently detected, in this case only occurring in groundwater.
Source factors, including sheep density and tillage land, and pathway factors within the zone of contribution (ZOC) of sampling sites, were statistically associated with anthelmintic occurrences. However, many of these associations were confounded by small sample numbers and small clusters of sites within the same catchment. Sampling point type and soil type were statistically related to anthelmintic occurrences, with detections more likely in monitoring boreholes compared to abstraction boreholes and in sites with ZOCs dominated by poorly drained soils. However, it was evident that the occurrence of anthelmintics in groundwater was not solely related to any one source or pathway factor at national scale, but rather to a combination of multiple factors on a more localised level, which varied on a site by site basis.
A temporal investigation indicated that the period with the highest anthelmintics occurrence was February/March following return of animals to pastures, pre-spring dosing and landspreading of manure, and again later during August/September. The outcome of this temporal investigation also highlighted the potential importance of unattenuated surface to groundwater pathways, such as sinking streams, for the transport of anthelmintics to groundwater.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.