Combining ddPCR and environmental DNA to improve detection capabilities of a critically endangered freshwater invertebrate

Isogenus nubecula is a critically endangered Plecoptera species. Considered extinct in the UK, I. nubecula was recently rediscovered (in one location of the River Dee, Wales), after 22 years of absence. In a similar way to many other species of Perlodidae, I. nubecula could be utilised as a bio-indicator, for assessing water quality and health status of a given freshwater system. However, conventional monitoring of invertebrates via kick-sampling, is invasive and expensive (time consuming). Further, such methods require a high level of taxonomic expertise. Here, we compared the traditional kick-sampling method with the use of eDNA detection using qPCR and ddPCR-analyses. In spring 2018, we sampled eDNA from twelve locations on the River Dee. I. nubecula was detected using kick-sampling in five of these locations, three locations using both eDNA detection and kick-sampling and one location using eDNA detection alone – resulting in a total of six known and distinct populations of this critically endangered species. Interestingly, despite the eDNA assay being validated in vitro and in silico, and results indicating high sensitivity, qPCR analysis of the eDNA samples proved to be ineffective. In contrast, ddPCR analyses resulted in a clear detection of I. nubecula at four locations suggesting that inhibition most likely explains the large discrepancy between the obtained qPCR and ddPCR results. It is therefore important to explore inhibition effects on any new eDNA assay. We also highlight that ddPCR may well be the best option for the detection of aquatic organisms which are either rare or likely to shed low levels of eDNA into their environment.


Results
Specificity and validation of eDNA assay using PCR, qPCR and ddPCR. The eDNA assay (primers and probe) designed in this study were species-specific in-silico and in-vitro with both conventional PCR and qPCR. The negative controls or samples with DNA from non-target species did not amplify with either method. For qPCR, we analysed the standard curve and compiled the limit of detection (LOD) and limit of quantification (LOQ) as per the MIQE guidelines 33,34 . The LOD was 6.82 × 10 −6 ng DNA µL −1 at 39.29 ± 2.00 Ct (i.e. Cycle threshold) and the LOQ was 6.82 × 10 −4 ng DNA µL −1 at 34.48 ± 0.95 Ct (Slope = −3.86, Y inter = 19.52, R 2 = 0.97, Eff% = 81.63) (Fig. 1). With ddPCR, five of the replicates from the dilution, equating to 0.08 pg of DNA yielded a positive detection for I. nubecula (mean 0.05 copy per µL −1 ). Interestingly, only one replicate from the next dilution series (0.016 pg) yielded a positive detection (0.08 copy per µL −1 ). All further dilutions and the negative controls were negative. However, as shown in other studies 35 , the lower LOD readings (with a 95% confidence level) can sometimes overlap with apparent artefacts seen in the negative controls. For this reason, we considered 0.08 pg of DNA to be the lowest amount able to be detected using ddPCR and only considered samples >0.5 copy per µL −1 (as in 35 ) to meet the threshold for a positive detection.
Kick sampling assessment. Populations of I. nubecula were identified at five different locations along the River Dee, and apparently absent at a further five (Fig. 1, Table 1). Abundance ranged from just one individual at two of the sites (W7 and W8), up to a highest density of 30 individuals at W3. Two of the sites (surveyed using the eDNA assay) were unable to be assessed via kick sampling due to dangerous access and weather conditions at the time of sampling ( Table 1).

comparison of qpcR versus ddpcR analyses. Despite the success of the assay (in-silico and in-vitro),
we were unable to amplify DNA via qPCR in any of the eDNA samples (Table 2). This was even true of the 'positive control eDNA sample' which consisted of 11 I. nubecula individuals housed in a 1 litre mesocosm for a period of one hour before filtering (see methods). During each run, the positive dilution range indicated the assay ran without any issue (Slope = −3.65/−4.05, Y inter = 19.22/26.46, R 2 = 0.98/0.99, Eff% = 76.46/88.03). In contrast, the ddPCR analysis revealed a positive detection of I. nubecula at four sampling locations (Fig. 2, Table 2). Interestingly, only two of the sites (W4 and W5) were positive with both undiluted DNA and diluted template (1:2). W2 was positive only when we utilised undiluted eDNA, and W7 only when we diluted 1 in 2 ( Table 2). Concentration of eDNA was relatively low and ranged from 0.6 to 0.14 copies per µL −1 across all samples. The 'positive eDNA' sample generated a much higher DNA concentration of 5.4 copies per µL −1 (undiluted) and 8.2 copies per µL −1 (diluted).
The site occupancy modelling approach did not reveal any significant effect of the environmental variables on the presence of eDNA or on the probability of detection (Tables 3 and 4). Probabilities of I. nubecula occurrence were relatively low and ranged from 0.45 to 0.53 (Table 4). Probabilities of eDNA detection at each sampling site ranged from 0.59 at site W5 (where all 'natural replicates' where found to be positive using ddPCR) to 0.27 at site W10, a site with high turbidity where no stonefly were found.

Discussion
In our study, we compared the use of kick-sampling and eDNA detection for monitoring a critically endangered bioindicator macroinvertebrate. While our eDNA detection approach using qPCR showed high sensitivity ( Fig. 1), with no false positive results during the validation process and assessment of the MIQE guidelines (Appendix 2),we were, however, not able to amplify DNA traces of I. nubecula in any of the eDNA samples. This  www.nature.com/scientificreports www.nature.com/scientificreports/ is surprising as one should expect positive detection at least in the five locations where we found the species via kick-sampling, and especially in the'positive eDNA' sample. These observations thus clearly pose doubts on the concept of eDNA using the qPCR methodology. Potential explanations for these false negative observations might be (i) an incorrect sampling protocol, (ii) the presence of PCR inhibitors in the DNA extracts, or (iii) a very limited shedding rate of the targeted species 36 . As previously shown, the sampling design of any eDNA based study can affect the reliability of detection 34 . In this case, we accounted for this by taking three natural replicates at each site and incorporating six technical PCR replicates per sample. Therefore, we believed our sampling protocol to be sufficient.   www.nature.com/scientificreports www.nature.com/scientificreports/ This left inhibition of the qPCR assay as the most likely reason for the false negative detections, as has been shown in a number of other studies 37,38 . One can assess for inhibition via the use of internal positive controls, such as spiked synthetic DNA or DNA from organisms different than the targeted species 36 . Limited detection or complete failure of such internal controls may then clearly show the occurrence of inhibition factors. If there is inhibition, two methods can be utilised to overcome this issue. The first method is to dilute the DNA extracted from the field sample 36 , whilst the second is the use of an inhibitor removal kit 36,39 . However, both methods have been shown to reduce the yield of target DNA in the extracted sample 36 . In our study, qPCR showed no results from the eDNA samples and so we hypothesised that inhibition may be an important driver for the false negative observations in this assay. We did not use an inhibitor removal kit, as we wanted to avoid reducing the amount of DNA extracted from the field samples. Instead, we opted to run the samples on a ddPCR with two different dilutions. The use of ddPCR worked and four sites revealed a positive signal, three of which mapped with the results from our kick-sampling survey. However, the influence of diluting the eDNA extracts also became apparent when analysing the results. eDNA was shown to be positive for only two of the sites (W4 and W5) regardless of dilution, one (W2) was positive only with undiluted eDNA template and another (W7) only when the extract was diluted 1:2 ( Table 2). This result indicates care should be taken with regard to dilution of extract in future eDNA-based studies (both for qPCR and ddPCR) and where possible multiple dilutions (starting from zero) should be run to give greater confidence in the results.

NR A NR B NR C NR A NR B NR C NR A NR B NR C
It is not surprising that ddPCR outperformed qPCR in this study and similar results have been shown before [40][41][42][43] . This is simply because ddPCR partitions any given sample into thousands of droplets, performing independent end-point PCR amplification on each droplet, thereby enabling the detection and quantification of very low amounts of DNA 41 . After amplification, the fluorescence of each droplet is measured allowing quantification of the targeted DNA (Bio-Rad's QuantaSoft software version 1.7.4.0917). This is in contrast to qPCR, which relies on the detection of PCR amplification, rather than amplification efficiency. Interestingly, the analysis of the 'positive eDNA sample' (11 I. nubecula in a 1 litre mesocosm for a period of one hour before filtering) showed an increase from 5.4 copies per µL −1 (undiluted) to 8.2 copies per µL −1 (diluted). This indicates that inhibition is still affecting the ddPCR (although not strong enough to block amplification in this instance). Future research should therefore explore the role of inhibition in eDNA based methods as a matter of urgency to ensure confidence in these tools remains high.
The very low I. nubecula eDNA concentrations in the samples (0.6 to 0.14 copies per µL −1 ) and the still relatively low concentration of eDNA in the 'positive eDNA' sample (8.2 copies per µL −1 ) indicates that this species may have very low shedding rates. However, as this is the first study to utilise ddPCR for detecting low populations of an endangered invertebrate in a fast-flowing river, we are unable to compare our results with other studies to date. Besides the fact that invertebrates are generally found to shed only limited amounts of eDNA in the water, potential other explanatory variables of these low levels of eDNA for this species could be the high flow rate of the river and low temperature during sampling. Sample were collected at the end of winter/beginning of spring, when environmental conditions such as high flow rates or flood events could have decreased and diluted the quantity of DNA traces. However, this was unavoidable for this species as I. nubecula emerges from March onwards 31,44 and so sampling time could not be altered.
Finally, when sampling for any eDNA study, it is useful to have an understanding of the ecology of the species under study, such as the species habits and preferred habitat in which it occurs. However, again, as I. nubecula was only recently rediscovered in the UK, there is very little information on this species 32 . Our site occupancy modelling approach was also unable to identify any specific variable that would have a significant effect on the probability of detection of this species using our eDNA assay (Fig. 3, Tables 3 and 4). This was not surprising however, as all the individual sites were on the same river system. Indeed, occupancy modelling is known to have certain limitations, mainly driven by the number of locations sampled and restricted range of environmental values collected 43 . In addition, the rarity of our target organism in this study and its likely high stochasticity (with regard to population distribution) would influence the models' outputs and ultimate usefulness. Thereby exploring the www.nature.com/scientificreports www.nature.com/scientificreports/ effects of the underlying environmental drivers on the distribution of I. nubecula remains difficult at the current time. Further work will therefore be necessary in order to increase our understanding of the ecology of I. nubecula if we want to optimize the sampling protocol and conservation plans for this species. However, the application of site occupancy modelling can become beneficial when prior survey data is combined with a more intensive survey effort. In such cases, a more informed post experimental understanding will be obtained 45,46 .
In conclusion, even if the highest standards of validation are undertaken in the design and implementation of an eDNA based PCR or qPCR assay [28][29][30] , false negative results can and do appear, driven by inhibition factors 36 , low shedding rates from the target species 18,47 and/or low population sizes 20 . In this case we are dealing with an extreme scenario, in which none of our eDNA samples showed any amplification via qPCR despite the fact that populations of I. nubecula were known to be present. However, we were able to get positive detection (using ddPCR) at five independent sites, three of which mapped against a physical detection of the species using kick sampling. Less than ten studies have (at the time of writing) utilised this technology for eDNA assays 35,[40][41][42][43][48][49][50] , but this is likely to increase significantly due to the apparent benefits observed in this study for example. We end by highlighting that negative results, derived from assays reliant solely on qPCR should be viewed with caution, for the reasons given above.

Methods primers and probe design.
A species-specific set of primers and probe, targeting the COI gene (Cytochrome C Oxidase subunit 1 mitochondrial gene) of I. nubecula was designed using the Geneious Pro R10 Software (https://www.geneious.com 51 . This assay amplifies a 124 bp fragment using the forward primer (5′-CCAGAAGCCTTGTAGAAAAC-3′), the reverse primer (5′-ACCCCGGCTAGATGAAGAGA-3′) and a probe (6-FAM-CCCCACTCTCTGCTGGAATT-BHQ-1). Specificity of the assay was assessed in-silico by comparing against sequences from 21 genetically similar invertebrate species, previously submitted to the NCBI (National Centre for Biotechnology Information; https://www.ncbi.nlm.nih.gov/) see Appendix 1 for full list. The specificity of the assay was tested in-vitro using PCR and qPCR, with DNA extracted from the nine invertebrate species (closely related or likely to be present in the same ecosystem). These included; I. nubecula, Leuctra hippopus (Kempny, 1899), Perlodes mortoni (Klapálek, 1906), Nemoura lacustris (Pictet, 1865), Leuctra geniculata (Stephens, 1836), Nemoura erratica (Claassen, 1936), Taeniopteryx nebulosa (Linnaeus, 1758), Diura bicaudata (Linnaeus, 1758) and L. fusca (Linnaeus, 1758). eDnA samples. 12 locations from the River Dee, were sampled for eDNA between 9 th March 2018 and 1 st of April 2018 ( Fig. 1 and Table 1). These locations were chosen following previous knowledge of historical observations in 1981 and 1982 32 . At each location, three independent (i.e. A, B and C) 1 L water samples (referred to here after as natural replicates) were collected using a 40 mL sterile polypropylene ladle and placed into a sterile plastic bag (Whirl-Pak ® 1242 ml Stand-Up Bag Merck ® , Darmstadt, Germany) 34 . Sub-samples were regularly collected from surface water downstream to upstream (to avoid disturbing sediments), across the width or the bank of the river, depending on the access and weather conditions following the method outlined in 52 . Each independent 1 L water sample was then filtered with a sterile 50 mL syringe (sterile Luer-Lock ™ BD Plastipak ™ , Ireland) through a sterile 0.45 μm Sterivex ™ HV filter (Sterivex ™ filter unit, HV with luer-lock outlet, Merck ® , Millipore ® , Germany). Sterivex filters were immediately placed in a freezer bag and stored at −80 °C until further analysis. At each location, new sterile equipment and disposable nitrile gloves were used during the sampling process to avoid contamination. A 'positive' eDNA sample was collected by creating an isolated mesocosm onsite, which consisted of river water from site W4 and 11 specimens of I. nubecula stored for 1 hour. Two negative control samples were additionally filtered in the field with sterile ddH 2 O in parallel with the natural samples, to control for potential cross-contamination during the workflow.
DnA extraction. DNA extraction from both the eDNA samples and the tissue samples (utilised for validating the assay) was done using the Qiagen DNeasy ® Blood and Tissue Kit. We followed the manufacturer's instructions for performing DNA extraction from tissue samples. Sterivex filters were extracted following the methods outlined in 53 ). All laboratory equipment was disinfected and decontaminated using UV-treatment prior to conducting any laboratory work. Laboratory equipment and surfaces were regularly disinfected using 10% bleach and absolute ethanol before conducting analyses. ddpcR. Digital droplet PCR was conducted using the Bio-Rad QX200 ddPCR system in a 20-μl total volume. Each reaction contained 10 μL Bio-Rad ddPCR supermix for probes (no dUTP), 750 nM of each primer, 375 nM probe, 3 µL DEPC water, and 4 µl template DNA. Twenty microlitres of the PCR mix was pipetted into the sample chambers of a Droplet Generator DG8 Cartridge (Bio-Rad, cat no. 1864008), and 70 μL of the Droplet Generation Oil for Probes (Bio-Rad, Cat No. 186-4005) was added to the appropriate wells. The cartridges were covered with DG8 Gaskets (Bio-Rad, cat no. 1863009) and placed in a QX200 Droplet Generator (Bio-Rad) to generate the droplets. After droplet generation, the droplets (40 μL) were carefully transferred to a ddPCR 96-well plate (Bio-Rad, Cat No. 12001925). The PCR plates were sealed with pierceable foil (Bio-Rad, Cat No. 181-4040). PCRs were performed using a C1000 Touch TM Thermal Cycler with a 96-well Deep Reaction Module (Bio-Rad). PCR conditions were 10 min at 95 °C, followed by 40 cycles of denaturation for 30 s at 94 °C and extension at 60 °C for 1 min, with ramp rate of 2 °C s-1, followed by 10 min at 98 °C and a hold at 12 °C. Droplets were then read on a QX200 droplet reader (Bio-Rad). All droplets were checked for fluorescence and the Bio-Rad's QuantaSoft software version 1.7.4.0917 was used to quantify the number of I. nubecula copies per µL. Thresholds for positive signals were determined according to QuantaSoft software instructions. All droplets beyond the fluorescence threshold (3500) were counted as positive events, and those below it as negative events. All eDNA samples were analysed in duplicate (one replicate undiluted and one replicate diluted 1:2). One positive control (i.e. DNA extracted from I. nubecula at a concentration of 1 ng/µL diluted 1:100 (Nanodrop 2000 Spectrophotometer, ThermoFisher Scientific)), one No Template Control (i.e., IDTE pH 5.0) and the two negative field controls were additionally included. The LOD using the ddPCR was assessed following the method outlined in 35 . We conducted a serial dilution of a DNA extracted from I. nubecula. The starting point was an initial 1: 100 dilution of extracted genomic DNA from I. nubecula at 1 ng/µL, followed by a serial 1:5 dilution. The serial dilution included ten replicate of each dilution. estimation of the LoD and LoQ. To attain estimates of the LOD and LOQ for the primers and probe for both qPCR and ddPCR, we set-up a dilution range from 10 −1 to 10 −9 with 10 technical replicates used for each of the dilution steps. Following 34 , the LOD was defined as the last standard dilution when the targeted DNA was detected and quantified in at least one replicate with a threshold cycle under 45. The LOQ was defined as the last standard dilution in which the targeted DNA was detected and quantified in at least 90% of positive samples 34,52 . All eDNA samples were then analysed with six technical replicates 34,52 on a qPCR plate, with six negative controls and a positive control dilution series from 10 −1 to 10 −6 in duplicate.
Kick-sampling. Kick-sampling was performed using the standard of the Freshwater Biological Association (UK), i.e. using a kick-sampling net with a 1 mm mesh (see detailed protocol: https://www.fba.org.uk/ practical-guidance-sampling-and-collecting). Sampling duration was recorded at each site but varied dependant on access, depth, flow rate, and/or weather conditions (Table 1). Perlodidae specimens found during kick-sampling were either preserved in 99% ethanol or kept alive as a part or a separate rearing experiment. Specimens were identified in the laboratory by two independent taxonomy experts (John Davy-Bowker & Michael Hammett) using a low-power binocular microscope with cold light source and using an identification key 44,54 . Statistical analysis. A site occupancy modelling approach 45,55,56 was utilised to assess the effect of environmental covariates on the presence of eDNA of I. nubecula and to estimate the detection probability of the new assay. This hierarchical modelling framework has the advantage of accounting for the risk of false negative results when estimating the probability of detection. This analysis was run with the ddPCR data (Appendix 3). Covariates tested included: (i) turbidity (likely to inhibit the PCR reaction, with the volume of filtered water being used as a proxy), (ii) pH, (iii) dissolved oxygen concentration, (iv) amount of time spent at each location (for both eDNA sampling and kick-sampling)used as a proxy for field conditions and (v) site accessibility as a binary indicator (possible to perform kick-sampling/absence of kick-sampling survey) (Appendix 4). Analyses were performed using the 'eDNAoccupancy' package 43,57 in the R statistical programming environment (R Core Team, 2018). Model selection and interpretation followed procedures given in 43,57 . We fitted our model using the 'occModel' function from the described package. MCMC chains ran for 11,000 iterations, with 10,000 retained for obtaining parameter estimates and credible intervals.