Identification of factors associated with Fasciola hepatica infection risk areas on pastures via an environmental DNA survey of Galba truncatula distribution using droplet digital and quantitative real-time PCR assays using

Environmental DNA (eDNA) is a powerful tool for identifying the spatial and temporal presence and density of species in a range of aquatic habitats. The analysis of eDNA has a wide range of application, one of which may be to inform of Fasciola hepatica infection risk on pastures based on the detection of its eDNA as well as that of its intermediate snail host, Galba truncatula eDNA. Here, droplet digital PCR (ddPCR) and quantitative real- time PCR (qPCR) assays were developed to detect the eDNA of F. hepatica , and its intermediate snail host, G. truncatula in water samples collected from pastures grazed by cattle and/or sheep. Environmental factors associated with species presence, as detected via an eDNA survey, were identified using zero-


| INTRODUC TI ON
The detection and analysis of environmental DNA (eDNA) is recognized as a powerful methodology for evaluating species presence across space and time (Thomsen & Willerslev, 2015). Significant progress has been achieved over the past decade in optimizing analysis of water samples for eDNA of rare, invasive, and pathogenic species (Cristescu & Hebert, 2018). Vast surface areas can now be screened both quickly and efficiently for evidence of target species presence and to estimate population densities and biomass based on eDNA presence and concentration (Thomsen & Willerslev, 2015).
Surprisingly, despite the successful application of eDNA analysis of environmental water in human health (Larsen & Wigginton, 2020) and ecology (Beng & Corlett, 2020), it is yet to be widely applied in veterinary health despite multiple possibilities for aiding disease control.
One of the areas within the veterinary field that has shown promise for the successful application of eDNA analysis is parasite control, in particular the management of trematode parasites whose larval stages are found in water at key periods of their lifecycle (Jones et al., 2018). Trematodes are also dependent on an intermediate host to complete their lifecycle, many of which are either aquatic or amphibious snails. One of the main trematodes of concern for livestock producers is the liver fluke, Fasciola hepatica, whose infections are estimated to cost the European livestock industry € 635 million annually due to decreased milk yields and growth rates, effects which will also increase the carbon footprint of livestock production (Jonsson et al., 2022), and increase in mortality and treatment costs (Charlier et al., 2020). F. hepatica also infects humans, with its disease, fasciolosis regarded as a neglected tropical disease by the WHO (Mas-Coma et al., 2018). This trematode uses Galba truncatula as its primary intermediate host, an amphibious freshwater snail which also hosts the rumen fluke Calicophoron daubneyi (Jones et al., 2017), an emerging parasite in western Europe which has the potential to cause clinical disease and mortality in ruminants . Controlling these trematode parasites is becoming extremely challenging globally, due firstly to the impacts of climate change, which has allowed G. truncatula snails to thrive in temperate climates (Fox et al., 2011) and secondly due to the emergence of anthelmintic resistance which limits treatments options postinfection (Castro-Hermida et al., 2021). These challenges highlight the importance of an integrated control strategy for F. hepatica and C. daubneyi in livestock, whereby infection risk areas on farmland are managed proactively to limit infection opportunities alongside responsible treatment of infected animals (Crilly et al., 2015).
One of the key measures required to effectively assess current infection risk is the presence of G. truncatula snails on pastureland, with inclusion of spatial information on habitats known to significantly improve the strength of predictive fasciolosis risk models (Charlier et al., 2011). Challenges are associated with accurately mapping G. truncatula snail distribution, with the snails difficult to find and identify and may not always be present in seemingly suitable habitats (Charlier et al., 2011). Furthermore, our understanding of microenvironmental factors that influence G. truncatula presence and population densities is limited which further restricts the effectiveness of wide scale remote sensing spatial models of G. truncatula distribution (Charlier et al., 2014).
In recent years, in an attempt to overcome these technical restraints, eDNA assays have been developed to detect the presence of G. truncatula in a range of water saturated habitat types, including boggy areas, ponds, streams, and ditches, as well as the detection of eDNA of both F. hepatica and C. daubneyi in these areas (Jones et al., 2018). Conventional PCR and LAMP assays have previously been designed and optimized for G. truncatula eDNA detection (Davis et al., 2020;Jones et al., 2018). However, these tools, despite their simplicity, low cost, and practicality, can lack sensitivity and tolerance to inhibitors, and may not be capable of detecting the minute quantities of eDNA that is commonly observed in environmental water (Fediajevaite et al., 2021). These methods are also incapable of quantifying eDNA levels within samples which may be beneficial to estimate species biomass (Yates et al., 2019). New G. truncatula eDNA amplification methods, such as quantitative real-time PCR (qPCR) and droplet digital PCR (ddPCR), may therefore further enhance previously designed eDNA protocols to deliver superior target detection (Capo et al., 2021;Rathinasamy et al., 2018;Wood et al., 2019) that will subsequently support efficient mapping and modeling of G. truncatula distribution and trematode infection risk on farmland (Jones et al., 2018).
The aims of this study were to (1) develop qPCR and ddPCR assays to detect and quantify G. truncatula eDNA in environmental water samples collected across the pasture environment, and to (2) evaluate factors associated with the spatial distribution of G. truncatula snails on pastures in a temperate environment as identified using eDNA analysis.

| Primer and probe design
Primers and probes targeting a 90-bp amplicon of the G. truncatula ITS2 gene (GenBank accession no. AJ296271.1) were designed in silico using the PrimerQuest™ Tool (Integrated DNA Technologies, Coralville, USA). The forward 5′-CAAAGCATTCGTGTCCTTGC-3′ and reverse: 5′-CAGGGCGTATCCAACCATC-3′ primers and FAM/ BHQ1 labeled probe 5′-AAACCGAAGCCTTGCGGCGT-3′ were tested for their specificity in silico using Genious Prime 2019 (https://www.genei ous.com) and NCBI Blast (NCBI, Bethesda, USA). F. hepatica primers and probes designed and validated for qPCR by Rathinasamy et al. (2018), which target a 108-bp amplicon of the F. hepatica ITS2 gene, were also used in the study to analyze field samples. These primers and probe were incorporated into the G. truncatula ddPCR assay at concentrations of 900 and 250 nM, respectively. The probe was dual labeled with HEX and BHQ1, which allowed for multiplexing alongside the FAM labeled G. truncatula probe.

| qPCR protocol
The Q-qPCR machine (Quantabio, Beverly, USA) was used for qPCR analysis. The qPCR master mix composition and thermocycling profile of the qPCR assay were optimized within the manufacturer's recommended range. The PCR master mix was made in a final volume of 10 μl, and included 1× PerfeCTa® MultiPlex qPCR ToughMix®, 500 nM of each primer, 250 nM of Probe, and 2 μl of DNA extract.
The optimal primer and probe concentration was determined via a primer/probe concentration titration experiment (500/250, 250/125, 200/100, and 100/50 nM primer/probe concentrations were tested) where the lowest cycle threshold (CT) value for G. truncatula DNA observed would indicate optimal concentrations. qPCR was run with an initial denaturation at 95°C for 3 min, followed by 45 cycles of denaturation at 95°C for 15 s and an annealing/extensions step of 60°C for 45 s. The optimal annealing/extension temperature was determined via an annealing temperature increment experiment (56, 58, 60, 62, and 64°C) where the lowest CT value for G. truncatula DNA would indicate optimal conditions.

| ddPCR protocol
The QX200 Droplet Digital PCR System (BioRad, Hercules, USA) was used for ddPCR analysis. The ddPCR master mix composition and thermocycling profile of the ddPCR assay followed the manufacturer's recommended protocols, with the exception that the number of cycles was increased to 45. The PCR master mix was made in a final volume of 20 μl, and included 1× supermix for probes (BioRad), 900 nM of each primer, 250 nM of Probe, and 4 μl of DNA extract.
The master mix was placed in sample wells of the droplet generation cartridge, with 70 μl of droplet generation oil (Bio-Rad) added to each of the oil wells. These cartridges were placed in a QX200 Droplet Generator (Bio-Rad) which generated microdroplets, which were transferred into a 96-well PCR plate (Bio-Rad), heat-sealed with foil using the PX1 Plate Sealer (Bio-Rad). PCR was run with an initial denaturation at 95°C for 10 min, followed by 45 cycles of denaturation at 94°C for 30 s, annealing/extensions step of 60°C for 60s, and a final extension step at 72°C for 5 min. The optimal annealing/extension temperature was determined via an annealing temperature increment experiment (57.4, 58.5, 60, 61.8, and 63.3°C) where the greatest separation between positive and negative droplets and the highest copy number of G. truncatula DNA detected would indicate optimal conditions. The droplets were analyzed individually using the QX200 Droplet Reader (Bio-Rad), with the data analyzed using QuantaSoft software (Bio-Rad).

| Assay optimization, sensitivity, and specificity
To test the sensitivity of both assays, a series of 10-fold dilutions of DNA extracted from G. truncatula snails was analyzed. This DNA was extracted following Jones et al. (2017), with the addition of a final DNA purification step using AMPure XP magnetic beads (Beckman Coulter, Brea, USA) following standard protocols. The DNA within each assay ranged from 2 log 10 ng to −5 log 10 ng (as measured by NanoDrop™ 2000 Spectrophotometer) across x8 ten-fold dilutions.
This analysis was conducted in technical triplicates, with the limit of detection (LOD) of each assay determined by the lowest detectable DNA concentration where all three replicates were positive. To determine the limit of quantification (LOQ) of the qPCR and ddPCR assays, the coefficient of variation among the three replicates at each dilution was calculated, with the LOQ determined as the last dilution where the LOQ was below 35% (Forootan et al., 2017). To further refine the LOD/LOQ estimates, a 2-fold dilution of initially calculated LOD DNA concentrations was analyzed in triplicate. To test the specificity of both assays, DNA extracted from Omphiscola glabra, Stagnicola palustris, Stagnicola fuscus, Lymnaea stagnalis, and Radix balthica snails (Davis et al., 2020) were run in biological duplicates for both assays, as was five environmental DNA extracts from farmland habitats known not to contain G. truncatula snails.

| Field eDNA sample collection
During April and May 2021 and May 2022, ten 50-ml samples of water were collected from 64 water saturated areas (500 ml in total) on six beef and/or sheep farms in Ceredigion, Wales, UK, with known cases of F. hepatica and C. daubneyi in their flocks/herds. The number of habitats sampled on each farm ranged from 8 to 16 (mean = 11).
These habitats were either watercourses, poached areas on pasture (wet areas where soil has been mechanically disturbed by animals and/or vehicles), or ponds and were present in fields grazed by livestock. The environmental and physical nature of each habitat was recorded ( Table 1). These data were later used for statistical analysis of factors associated with G. truncatula eDNA detection. The pH of water in each habitat was also measured using an electronic pH meter (Hanna Instruments, Rhode Island, USA). During water collection, a concurrent search for G. truncatula snails was completed where snails within 50 cm radius of each 10 collection spots were counted. All visits, sample, and data collection were conducted by the same individual on each occasion.
Each 500-ml water sample was transported in chilled conditions to the laboratory, where the samples were filtered through 2.7μm micro-glass fiber filters (Whatman, Maidstone, UK) using a faucet vacuum suction pump, Büchner flask, and a funnel. Any soil sediments present at the bottom of the water sample were not filtered to limit the presence of inhibitors in the samples. DNA from half of each filter was extracted using an adapted DNeasy® PowerSoil® kit (Qiagen, Hilden, Germany) protocol, where the initial bead-beating step was undertaken in a 7-ml bijou tube. Extracts were eluted in 70 μl of elution buffer. Negative controls included double distilled water samples that were transported alongside field samples and processed following the same protocol to identify potential DNA contamination during the eDNA collection and filtering process, and blank filters to identify potential DNA contamination during the DNA extraction process. Water filtering, DNA extraction, and assay setup/DNA amplification were completed in separate laboratories (Goldberg et al., 2016).

| eDNA analysis
Field samples were analyzed in technical duplicates using both the qPCR and ddPCR assays, respectively. The qPCR assay incorporated a standard curve run in duplicate in each run, where a 10-fold dilution series (1 log 10 ng -−3 log 10 ng) was used as a reference to quantify DNA concentration in each sample. For the second ddPCR replicate, the multiplex assay was used to enable detection of both G. truncatula and F. hepatica eDNA. This multiplex assay followed the same protocol that was optimized in this study for the singleplex G. truncatula assay. To test if significant differences in G. truncatula eDNA concentrations were identified whilst using the singleplex and multiplex assays, a paired t-test was used to compare the log 10 transformed eDNA concentrations identified by both assays.
Each ddPCR and qPCR run incorporated a positive control (G. truncatula DNA extract) and non-template control samples. An automatic threshold was set in the Q-qPCR software to identify positive samples in the qPCR assay, whilst the ddPCR threshold for detection was set at 2000 for both G. truncatula and F. hepatica detection as these thresholds attained optimal separation of negative and positive droplets. For data reporting and visualization purposes, the concentration of eDNA in 1 L of field water per habitat was calculated following the equation described by Agersnap et al. (2017).
As the ddPCR assay identified more positive G. truncatula samples compared to the qPCR assay, all negative qPCR samples were reanalyzed, with an adaptation to the protocol, where 4 μl of DNA extract was added to each reaction. Furthermore, to test if inhibitors were impacting eDNA detection, all negative field samples were reanalyzed using the qPCR assay whilst spiked with −3 log 10 ng of G. truncatula DNA extract.

| Statistical analysis
All statistical analyses were conducted in R (R Core team, 2019), primarily using the glmmTMB package (Brooks et al., 2017). Preliminary analysis was conducted to establish which generalized linear mixed model (GLMM) family would be most appropriate to analyze field eDNA data. Here, the deviance of simulated residuals of both normal and zero-inflated GLMMs of the Gaussian family, using log 10 TA B L E 1 Environmental and physical factors recorded for each habitat during eDNA survey, along with the number of habitats that exhibited each factor transformed copy concentration as the dependent variable, and the Poisson and negative binomial (nbinom1 and nbinom2) families, using raw copy count data as the dependent variable, were compared using the DHARMa package in R (Hartig, 2020). Zero-inflated GLMM of the Gaussian family was deemed to be superior to other models tested due to an observed lower deviance of simulated residuals, and thus models belonging to this family were created to identify factors asso- The study also aimed to identify factors associated with F. hepatica eDNA detection using a binomial mixed model where farm and field were included as random factors. However, no significant associations between F. hepatica eDNA presence and any factor were detected (p > 0.05).

| Assay optimization and performance
The optimal annealing/extension temperature for both qPCR and ddPCR assays was determined to be 60°C. For ddPCR, an annealing/extension step of 60°C yielded highest copy numbers and the greatest separation between negative and positive droplets habitats, a mean eDNA concentration of 4.76 log 10 copies/L was detected by ddPCR analysis, a significantly higher eDNA concentration compared to positive habitats where G. truncatula snails were not physically detected (3.24 log 10 copies/L) (p = 0.013). ddPCR also identified 5 of 64 samples to be positive for F. hepatica eDNA, with a mean eDNA concentration of 2.24 log 10 copies/L. There was no significant difference between log 10 transformed G. truncatula DNA concentrations observed in samples when analyzing using the singleplex ddPCR assay and the multiplex ddPCR assay (p = 0.959).
Examples of habitats sampled and the eDNA concentrations detected can be seen in Figure 6.
Factors identified as being significantly associated with G. truncatula eDNA detection and the log 10 number of G. truncatula eDNA copies per liter of field water in each habitat as identified by ddPCR can be seen in Table 2. Habitats with black or dark brown soil, those that contained deep water pools and those that were fully shaded were significantly less likely to contain G. truncatula eDNA (p < 0.05), whilst higher eDNA concentrations were found in temporary G. Mean log 10 target DNA copies detected three technical replicates via ddPCR are plotted against G. truncatula genomic DNA concentrations (log 10 ng/assay) used in the dilution series. 1 and 2 log 10 ng led to oversaturation of positive droplets in ddPCR. X markers indicate mean log 10 copies/CT value of DNA concentration falling below LOD threshold as one replicate did not amplify in each respective assay.

F I G U R E 3
Association between eDNA concentrations of G. truncatula in eDNA extracts as measured using qPCR and ddPCR where a significant positive association was present (p < 0.001). Red markers indicate samples that were positive via ddPCR but negative via qPCR analysis.
for unicellular algae, the main food source of G. truncatula, which grows on bare mud surfaces (Dreyfuss et al., 2015). Shaded areas may, therefore, lack the capability of supporting populations of G.
truncatula snails due to limited food (Petzold, 1989). Shaded areas are also likely to be cooler which may impact snail population size as G. truncatula snails are increasingly active, grow quicker, tend to live longer, and produce more offspring as temperatures increase in temperate conditions (Hodasi, 1976) and traits that may also lead to increased eDNA excretion (Jones et al., 2021). However, contradictory positive associations between trees and G. truncatula presence were identified by Roessler et al. (2022) in Switzerland. It was concluded there that the presence of single trees or bushes as well as a certain proximity to woodland is beneficial to the snails as trees offer protection from direct sunlight and drying out. In our work, no significant difference was observed between eDNA detection rates in habitats partially shaded and those not shaded. This suggests that some shading is not completely detrimental to G. truncatula populations, and when considering the findings of Roessler et al. (2022), some shading may be beneficial in certain circumstances to protect habitats from desiccation. Furthermore, a near significant difference in eDNA detection was also observed when comparing habitats partially shaded and those fully shaded which further emphasizes that the degree of shading is vital to determining habitat suitability to G truncatula. Currently, major drives are underway to plant trees on agricultural land with the aim of increasing carbon storage in efforts to reach net zero carbon emissions within the agriculture sector.
Within this drive to plant more trees, it has been suggested that strategically planted trees could be employed to limit habitat suitability to G. truncatula, especially considering their ability to improve soil drainage (Crilly et al., 2015). However, questions remain regarding the impact of strategic planting of trees on G. truncatula populations and the subsequent dynamics of trematode transmission in those areas. To answer these questions, eDNA analysis could be employed in future studies to monitor changes in G. truncatula population dynamics over time as trees, or any other mitigations are employed.
Soil color was also associated with eDNA detection, with eDNA detection less likely in soils categorized as dark brown or black.
Darker soils often contain high levels of organic matter and peat and

F I G U R E 5
Boxplot of mean (×) and median (−) G. truncatula (n = 42) and F. hepatica (n = 5) eDNA concentrations across habitats sampled in study as measured by ddPCR.
are often associated with lower pH values. Studies have shown that G. truncatula snails can tolerate a wide range of soil pH, although they tend to avoid highly acidic soils (pH <5.1) (Dreyfuss et al., 2015).
Increasing G. truncatula population densities are often positively correlated with soil pH (Caron et al., 2014;Charlier et al., 2014) and this could explain the negative associations between these soils and the presence of detectable eDNA; however, there was no significant association between eDNA detection and sampled water pH. It has been suggested that peaty soils are less suited to G. truncatula snails due to an absence of suitable food sources (Heppleston, 1972) and its texture (Walton & Wright, 1926), rather than its acidity, and the darker soils observed in this study could be less suited to G. truncatula due to these alternative reasons. Furthermore, darker soils may have high concentrations of specific inorganic compounds such as manganese oxide (Dixon et al., 1990), yet the impact of soil mineral content on G. truncatula has yet to be studied in detail despite the Significantly higher eDNA concentrations were also found in samples collected from temporary G. truncatula habitats (i.e., those that tend to dry up for periods in the summer). This finding contradicts the widespread notion that G. truncatula snails are found in higher densities in permanent habitats as snails in temporary habitats have increased mortality rates due to temporal variations in moisture levels which may vary from standing water to extremely dry ground F I G U R E 6 A selection of habitats surveyed in this study. (a) The habitat with highest G. truncatula eDNA concentration in study (424,000 copies/L); (b) a poached area on pasture (83,000 copies/L); (c) a Tire rut habitat (124 copies/L), which also had the highest concentration of F. hepatica eDNA in the study (373 copies/L); (d) a wet flush on pasture which had only recently appeared due to the failure of underground drains (6767copies/L); (e) drainage ditch habitat where G. truncatula snails were physically detected (98,000 copies/L); (f) a negative habitat, note the that this habitat was shaded by trees which were present within the habitat. (Heppleston, 1972). Nevertheless, it has been found that G. truncatula snails inhabiting temporary habitats tend to be larger (Chapuis et al., 2007;Heppleston, 1972) which may cause an increase in eDNA levels expelled per individual snail which may explain the increased eDNA concentrations observed in this study in temporary habitats. Furthermore, fluctuations in habitat water volume may dilute or concentrate eDNA (Katz et al., 2021;Seeber et al., 2019), and the increased eDNA concentrations observed in temporary G. truncatula habitats may be, in part, caused by the drier than usual conditions observed during the study period. This hypothesis is supported by the fact that G. truncatula eDNA was less likely to be detected in habitats with deeper water pools, although in this instance it is widely accepted that G. truncatula snails avoid deeper pools due to their requirement to access air unlike other lymnaeid species (Skuce et al., 2014). This highlights a major challenge in quantifying and comparing eDNA concentrations in habitats of varying types, sizes, and water volume as is observed for G. truncatula habitats which may take the form of permanent wetlands, spring heads, ponds, streams, ditches, and poached soils on pastures in both temperate and arid climates (Dreyfuss et al., 2015).  (Parr & Gray, 2000). Further longterm solutions were also identified for future application which will aim to limit G. truncatula habitat suitability. For example, on two of the farms, G. truncatula habitats were identified in wet flushes where field drainage systems had failed ( Figure 6) and fixing these drains would likely make these habitats uninhabitable to G. truncatula snails (Dreyfuss et al., 2015). Furthermore, multiple watercourses were identified across the farms ( Figure 6) and fencing to limit livestock access could be a viable option in some of these areas.  (Cuthill, 2020). However, major practical challenges are associated with screening non-liquid matrices such as plants and soil for eDNA, which include limited sample volume capacity per extraction and the presence of PCR inhibitors (Cuthill, 2020;Thomsen & Willerslev, 2015). Considering that water is the optimal screening matrix for eDNA, other molecular markers also present in water and that may be specific for distinct F. hepatica life stages, such as protein or RNA (Pawluk et al., 2018;Veilleux et al., 2021), could be targeted in future environmental water-based surveys, although substantial research would need to be undertaken initially to characterize proteins and/or RNA markers that are specific to each F. hepatica life cycle stage (Cwiklinski & Dalton, 2018).
In conclusion, this study developed novel ddPCR and qPCR assays to detect G. truncatula eDNA in water from pastureland. These assays can be utilized to map F. hepatica infection risk areas on farmland which can guide disease control strategies in flock and herds.
Surveillance of G. truncatula eDNA and potentially F. hepatica eDNA can also be used in future to further our understanding about G.
truncatula ecology and the F. hepatica lifecycle, and to monitor the effectiveness of management strategies aimed at limiting habitat suitability for G. truncatula or to limit contact between these snail and grazing livestock. The tools developed may also be applied in a human health context, specifically in regions of the world where F. hepatica commonly infects humans.

CO N FLI C T O F I NTE R E S T
The authors declared that they have no conflict of interest.

DATA AVA I L A B I L I T Y S TAT E M E N T
The data that support the findings of this study are available from the corresponding author upon reasonable request.