Development of highly sensitive environmental DNA methods for the detection of Bull Sharks, Carcharhinus leucas (Müller and Henle, 1839), using Droplet Digital™ PCR

Background: As apex and mesopredators, elasmobranchs play a crucial role in main‐ taining ecosystem function and balance in marine systems. Elasmobranch populations worldwide are in decline as a result of exploitation via direct and indirect fisheries mortalities and habitat degradation; however, a lack of information on distribution, abundance, and population biology for most species hinders their effective manage‐ ment. Environmental DNA analysis has emerged as a cost‐effective and non‐invasive technique to fill some of these data gaps, but often requires the development of spe‐ cies‐specific methodologies. Aims: Here, we established eDNA methodology appropriate for targeted species de‐ tections of Bull Sharks, Carcharhinus leucas , in estuarine waters in the northern Gulf of Mexico.


| INTRODUC TI ON
Elasmobranchs (sharks, skates, and rays) play a crucial role in marine ecosystems as apex and mesopredators, influencing prey abundance, behavior, and trophic interactions across multiple trophic levels in marine food webs (Ferretti, Worm, Britten, Heithaus, and Lotze 2010;Ritchie et al. 2012). Healthy elasmobranch populations help to maintain ecosystem function, increase biodiversity, and buffer against invasive species and transmission of diseases (Heithaus, Frid, Wirsing, and Worm 2008;Ritchie et al. 2012). However, many elasmobranch populations are in decline as a result of exploitation via direct and indirect fisheries mortalities and habitat degradation (Dulvy et al. 2014). The life history strategies of many elasmobranchs are characterized by late maturity, longevity, and low fecundity, making the recovery of exploited populations a biologically slow process (Garcia et al., 2008;Hoenig and Gruber 1990). According to the International Union for Conservation of Nature (IUCN) Red List of Threatened Species, one-quarter of elasmobranch species are estimated to be threatened with extinction and almost one-half are categorized as Data Deficient, meaning there are insufficient data to properly assess their conservation status (Dulvy et al. 2014). Robust data on species distribution, abundance, biology, and population biology are necessary to enact appropriate conservation strategies for the maintenance of healthy elasmobranch populations; unfortunately, such data are often incomplete or lacking for many species (Dulvy et al. 2014).
Analysis of environmental DNA (eDNA) has recently emerged as an alternative, powerful approach to fill data gaps on the distribution, habitat use, abundance, and population biology of aquatic species (Ficetola, Miaud, Pompanon, and Taberlet 2008), including elasmobranchs (Sigsgaard et al. 2016). All organisms leave traces of DNA in the environment through shedding of cellular debris, skin cells, blood, and biological waste, all of which can be collected in water samples (Rees, Maddison, Middleditch, Patmore, and Gough 2014); however, differences in how organisms shed DNA (i.e., mucus, scales, feces) suggest that eDNA accumulation may differ across species (Le Port, Bakker, Cooper, Huerlimann, and Mariani 2018), requiring taxon-specific research. In targeted species detections, water samples are typically filtered, DNA extractions are performed on the resulting particulate material, and extracted DNA samples are analyzed using a quantitative real-time polymerase chain reaction (qRT-PCR) platform with species-specific primers, developed to amplify a small DNA fragment in the target species (Foote et al. 2012;Taberlet, Coissac, Hajibabaei, and Rieseberg 2012). The collection of water samples is a cost-effective and efficient method of surveying elasmobranch populations when compared to traditional survey methods involving setting nets or lines, which can have high incidence of bycatch and inflict varying degrees of stress to both target and nontarget species (Larson et al. 2017;Lewison, Crowder, Read, and Freeman 2004). Post-release recovery and survival tends to vary widely across species, with some species being particularly sensitive to net capture and handling (Stobutzki, Millter, Heales, and Brewer 2002). With a well-designed sampling scheme, eDNA methodologies offer increased sensitivity for detecting the presence of rare species while negating the need to capture, handle, or even observe the target species (Port et al. 2016;Rees et al. 2014). In elasmobranchs, eDNA methods have been used in targeted species detections for the Critically Endangered Largetooth Sawfish, Pristis pristis (Simpfendorfer et al. 2016), the Endangered Maugean Skate, Zearaja maugeana (Weltz et al. 2017), the Vulnerable Chilean Devil Ray, Mobula tarapacana (Gargan et al. 2017), and the Vulnerable Great White Shark, Carcharodon carcharias (Lafferty, Benesh, Mahon, Jerde, and Lowe 2018). Furthermore, eDNA has been used to assess population characteristics in the Endangered Whale shark, Rhincodon typus (Sigsgaard et al. 2016) and to estimate shark diversity in tropical habitats using metabarcoding (Bakker et al. 2017;Boussarie et al. 2018).
Bull Sharks, Carcharhinus leucas (Müller and Henle, 1839), are found in temperate, subtropical, and tropical latitudes globally and are distinctive as one of only a few sharks that can use freshwater for extended periods of time (Thorson 1962;Thorson 1971;Thorson, Cowan, and Watson 1973). As upper trophic level predators, they play a crucial role in maintaining ecosystem health across both marine and freshwater habitats (Every, Pethybridge, Fulton, Kyne, and Crook 2017;Polovina, Abecassis, Howell, and Woodworth 2009;Ritchie et al. 2012). Using acoustic telemetry data to examine the habitat use of C. leucas in northern Gulf of Mexico waters, Drymon et al. (2014) found C. leucas may preferentially select higher-quality, less-urbanized rivers, although a spatially limited acoustic array hindered a full evaluation of this pattern. Targeted eDNA surveys of C. leucas could provide a cost-effective, sensitive method to examine this pattern more widely, as there could be substantial ecological implications of such habitat preference. Here, we establish an eDNA methodology appropriate for targeted species detections of C. leucas in estuarine waters in the northern Gulf of Mexico. Specifically, we compare total eDNA yields for different QIAGEN ® DNeasy ® DNA extraction kit protocols and develop a species-specific C. leucas eDNA assay using a relatively novel, Bio-Rad® Droplet Digital™ PCR (ddPCR), platform to detect low quantities of target DNA.
Finally, we apply these methods to investigate the detectability of C. leucas eDNA in known habitat in the northern Gulf of Mexico and in ex situ closed and flow-through environments containing a single C. leucas individual.

| Laboratory controls
Strict laboratory controls were implemented throughout this study to reduce the risk of cross-contamination and contamination by exogenous DNA (see Deiner, Walser, Mächler, and Altermatt 2015;Goldberg et al. 2016). Water processing, DNA extractions, and PCR amplifications were conducted in physically separated laboratory spaces to prevent cross-contamination between stages.
Negative controls were incorporated into every stage of sample processing, and PCR was performed on them to check for potential contamination. Filter negatives contained target-free, autoclaved deionized water, DNA extraction negatives contained no filtered particulate material, and PCR amplification negatives contained no DNA; all negative controls produced negative results, indicating no contamination had occurred. The ddPCR assay conditions used to carry out these negative control tests are described below.

| Water sample collection and filtration
Water samples throughout this study were collected just below the surface of the water in 1 L high-density polyethylene Nalgene ® bottles precleaned in a 10% bleach solution and sanitized under ultraviolet (UV) light for 20 min. New gloves were used to collect each water sample and samples were stored on ice in a cooler until filtration using a vacuum pump could take place, which occurred within 24 hr of collection (see Pilliod et al. 2013), except where otherwise noted. Water samples were filtered in a dedicated, precleaned laboratory space that had never had C. leucas tissue or total genomic DNA (gDNA) present.
Each 1 L water sample was inverted at least three times to ensure homogenization of particulate matter and was then vacuum-filtered using 47-mm-diameter, 0.8-μm nylon filters, which were replaced when clogging occurred every ~350 ml (e.g., three filters per 1 L) and preserved in 95% ethanol at room temperature, unless noted otherwise (see Appendix S1). During all water filtration, filters were handled with designated sterile forceps for each sample and gloves were changed in between samples to avoid cross-contamination.

| DNA extraction methods
Due to the wide variety of DNA extraction methods used in eDNA literature (Renshaw, Olds, Jerde, McVeigh, and Lodge 2015), we compared eDNA extraction kits to establish an appropriate method for the nylon filters used to filter water samples in this study. The QIAGEN ® DNeasy ® Blood & Tissue Kit is a frequent choice for DNA extractions from filters in eDNA studies, but with numerous variations (see Rees et al. 2014). The performance of this kit using the Goldberg et al. (2011) variation incorporating QIAshredder™ spin columns was compared to that of an extraction kit designed specifically for water samples, the QIAGEN ® DNeasy ® PowerWater ® Kit. The Goldberg et al. (2011) protocol incorporating QIAshredder™ spin columns was selected because in preliminary trials, it yielded higher relative quantities of DNA compared to some other variations (Appendix S2). Additionally, four variations of physical disruption methods to dislodge the particulate matter from the filters prior to digestion were tested with each extraction method: (a) no physical disruption, (b) bead beating, (c) filter scraping, and (d) freezing filters with liquid nitrogen and crushing them using an autoclaved mortar and pestle. The QIAGEN ® DNeasy ® PowerWater ® Kit contained bead beating as part of the standard manufacturer's protocol, so this step was eliminated for the no physical disruption variation to determine if this step was a critical factor in DNA yields. Three × 1 L water sample replicates were used in each extraction/physical disruption treatment, collected from Mobile Bay, Alabama using the water collection and filtration protocols described. To eliminate the filter preservation step, the filters for each 1 L sample were immediately placed into the appropriate lysis buffers (see Hinlo et al. 2017). The DNA extracts for each 1 L water sample were combined and the DNA qualities were assessed using 2% agarose gel and the relative quantities were measured using Thermo Fisher Scientific NanoDrop™ spectrophotometer technology, with each extract measured four times.

| Development of a species-specific assay
To develop a species-specific assay, primers and an internal probe were manually designed in conserved regions of the mitochondrial   (Table 1)  (Instrument no. 773BR1456) and thermal cycling conditions were as follows, using a ramp rate of 1°C/s: initial denaturation at 95°C for 10 min, followed by 35 cycles of 94°C for 30 s and 56°C for 2 min, followed by enzyme deactivation at 98°C for 10 min, and a final hold at 4°C. To ensure the optimized assay was species-specific for C. leucas using the ddPCR platform, the primers and probe were tested using these ddPCR reaction and cycling conditions, in replicates of three, with 0.2 ng/μl of gDNA extracted from five C. leucas individuals and one individual of each of 18 other genetically similar, local exclusion species, collected from the Gulf of Mexico (Table 1).
All ddPCR data were analyzed with the Bio-Rad ® QX200™ Droplet Reader and QuantaSoft™ software using the Rare Event Detection (RED) analysis, a manual detection threshold of 3,000 amplitude (Figure 1), and a limit of detection (LoD) of the developed assay. The LoD is considered the lowest concentration of C. leucas DNA that can reliably be detected using the optimized assay conditions. The lower LoD was determined by conducting ddPCRs with gDNA from two C. leucas individuals using a sixfold series of 10X dilutions (e.g., 1:10 to 1:1,000,000), from a starting concentration of 25.0 ng/μl. Means and standard errors of detected DNA concentration (copies/μl) were calculated for each individual, across the three ddPCR replicates for each dilution. The flow rate of the mesocosm was designed to mimic flow in a coastal system at ~30 cm 3 /hr, with complete system turnover at approximately 2 hr. One wild-caught juvenile male C. leucas, ~930 mm total length, was introduced to this system and five × F I G U R E 1 Raw output of the optimized Droplet Digital™ PCR (ddPCR) for the designed Carcharhinus leucas specific assay showing one ddPCR replicate for one individual (0.2 ng/μl of genomic DNA) and one replicate for the ddPCR negative from the Bio-Rad ® QX200™ Droplet Reader. Each droplet in each well was classified as either positive (blue droplets) or negative (gray droplets) for target DNA, based on a manual detection threshold set to 3,000 amplitude (the horizontal pink line) using the QuantaSoft™ Rare Event Detection analysis. Event number refers to the number of droplet events generated for a given well or sample; Ch 1 amplitude measurement refers to the level of fluorescence emitted by a droplet event; and each column is a single well 1 L water samples were collected immediately (time 0.0), spanning the diameter of the mesocosm; this sampling regime was repeated every 0.5 hr for 3 hr, allowing for complete turnover of the system.

| Collection of positive water samples
Water samples were stored in a −20°C freezer for 1 month, due to laboratory equipment constraints, similar to Bakker et al. (2017) and Gargan et al. (2017), and were thawed at room temperature prior to filtration.
Water samples from these experiments were vacuum-filtered using 47-mm-diameter nylon 0.8-μm filters (three per 1 L), which were preserved in 95% ethanol at room temperature (Appendix S1) and DNA extractions followed the Goldberg et al. (2011) protocol incorporating the QIAshredder™ spin columns (Appendix S2).
DdPCR amplifications were carried out in replicates of five, using the optimized C. leucas assay previously described in this study. All ddPCR reactions were set up using aerosol barrier filter pipette tips and designated pipettes, separate from those used in setting up PCR reactions, were used to add eDNA extracts to the reactions. DdPCR results were analyzed using the Bio-Rad ® QX200™ Droplet Reader and QuantaSoft™ RED analysis, a manual detection threshold of 3,000 amplitude, and the LoD.

| Optimal eDNA methods
The

| Analysis of water samples
Using the developed ddPCR assay and the QuantaSoft™ RED analysis with a manual detection threshold of 3,000 amplitude, an average of 1.62 copies/μl (SE = 0.12) of C. leucas DNA was detectable in the ddPCR reactions from water samples collected from known habitat, Mobile Bay, without visually confirming the presence of C. leucas ( Figure 4). In the ex situ positive eDNA experiment, 30 min after a C. leucas was added to the closed tank containing this water, large amounts of target eDNA were present, with an average concentration of 166.6 copies/μl (SE = 3.01) in the ddPCR reactions ( Figure 4).
In the flow-through mesocosm experiment, when applying a lower LoD of 0.6 copies/μl to the data analysis, target C. leucas DNA was not detectable in any of the ddPCR replicates at time 0.0 but was detectable in all ddPCR replicates 0.5 hr after the shark was added

| D ISCUSS I ON
The use of eDNA as a tool to study the distribution and ecology of marine species has increased substantially in recent years (Bakker et al. 2017;Foote et al. 2012;Lafferty et al. 2018;Port et al. 2016).
However, careful consideration and optimization of the methods employed in such studies are necessary, ultimately allowing for an appropriate interpretation of the results. Here, we found filtering water with nylon 0.8-μm filters, preserving the filters in 95% ethanol (Appendix S1), and then performing DNA extractions using the Goldberg et al. (2011)  LoD determined in this study shows the sensitivity and detection capability of the developed assay and was demonstrated to be sufficient for C. leucas eDNA detection in Mobile Bay and in ex situ positive samples. However, the LoD may require further refinement through additional dilution series between the 1:10,000 and 1:100,000 dilutions before being used in data analysis for large numbers of field samples. Furthermore, due to potential differences across ddPCR machines, we recommend the LoD to be refined independently for each machine, using the LoD here as a starting reference point for this assay.
The ability of ddPCR to detect low concentrations of target DNA, for example, 2.5 pg of C. leucas DNA in this study, means this platform may be less likely to produce false negatives when used alongside an appropriate sampling regime and water processing methods (e.g., spatial and depth coverage, volume collected, filter pore size). False negatives can occur when target DNA is captured in water samples but is not detected due to limitations of the genetic assays employed (  Limit of detection (LoD) tests using a 6-fold 10X dilution series (1:10-1:1,000,000) of total genomic DNA (gDNA) from two Carcharhinus leucas individuals from the northern Gulf of Mexico. (a) The mean DNA concentrations (copy number/ μl) and standard error bars were calculated from three Droplet Digital™ PCR (ddPCR) replicates for each of two individuals, using a manual detection threshold of 3,000 amplitude and the Rare Event Detection analysis setting on the Bio-Rad ® QX200™ Droplet Reader and QuantaSoft™ software. The 1:10 and 1:1,000,000 were not graphed due to oversaturation of the PCR product, and the lack of DNA copies present showing no positive droplet detections, respectively. The LoD (0.6 copies/μl) is represented by a dotted line. (b) Raw droplet output of ddPCR serial dilution products from one ddPCR replicate of one C. leucas individual detected by the Bio-Rad ® QX200™ Droplet Reader and QuantaSoft™ software. Each droplet in each well was classified as either positive (blue droplets) or negative (gray droplets) for target DNA. Each well is separated by yellow bars and corresponds to the same dilution concentrations graphed in Figure 3a, labeled with each dilution series it represents difference in detection abilities between the two PCR platforms is likely due to fundamental differences in how they quantify target DNA. DdPCR quantifies the starting DNA copy number present in a sample using end-point PCR without reference to a standard (absolute quantification) (Whale et al. 2012), making it a more sensitive and precise assay, ideal for eDNA applications targeting a single target species. Additionally, the RED analysis setting using the Bio-Rad ® QuantaSoft™ software is designed to identify low copy numbers of target DNA in a background largely composed of nontarget DNA copies (Bio-Rad ® Droplet Digital™ PCR Applications Guide).
Given the ability of ddPCR to detect such low quantities of DNA, it may replace qRT-PCR in eDNA research Nathan, Simmons, Wegleitner, Jerde, and Mahon 2014) assessing the distribution, habitat use, and abundance of species found in low abundance and/or are of conservation concern (Baker et al. 2018; Hunter et al. 2018;Tréguier et al. 2014), including elasmobranchs (Bohmann et al. 2014;Lafferty et al. 2018). However, we caution that the ability to detect such low quantities of DNA also increases the potential for false positives (Goldberg et al. 2016;Huggett, Cowen, and Foy 2015). All eDNA studies, but especially those using ddPCR, require strict field and laboratory controls and procedures be in place to reduce the potential for false positives, typically the result of contamination by exogenous DNA or cross-contamination of samples (see Ficetola, Taberlet, and Coissac 2016). In addition to the contamination controls described by Goldberg et al. (2016), Deiner et al. (2015), and Port et al. (2016), when using ddPCR, we also suggest: (a) using two cleaning methods for decontamination of all field and water filtration equipment (e.g., a bleach wash, plus autoclaving, and/or UV light exposure), (b) that water filtration is conducted in a laboratory space that has never had tissue or gDNA from the target species present, (c) that gloves and any tools are changed in between samples during water filtration (see Pilliod, Goldberg, Arkle, and Waits 2013), (d) that negatives be incorporated into field collection, water filtration, DNA extraction, and PCR, with each negative run through to PCR (see Bakker et al. 2017;Jerde, Mahon, Chadderton, and Lodge 2011), (e) that a designated pipette, separate from that used to set up reactions, be used to add DNA extracts to ddPCR reactions, and (f) that multiple replicates for each sample are run during ddPCR (see Rees et al. 2014).
Strict field and laboratory controls will ensure the authenticity and reliability of eDNA results, which is increasingly critical in eDNA research using highly sensitive technologies, such as ddPCR, especially when the results of such studies will be used to inform conservation and management initiatives (Hunter et al. 2017;Hunter et al. 2018).
F I G U R E 4 Raw Droplet Digital™ PCR (ddPCR) output from the ambient water sample in Mobile Bay, the Carcharhinus leucas eDNA positive water sample taken from a closed system 30 min after adding the shark, and each negative control from the Bio-Rad ® QX200™ Droplet Reader. Each droplet in each well was classified as either positive (blue droplets) or negative (gray droplets) for target DNA based on a manual detection threshold set to 3,000 amplitude (the horizontal pink line) using the QuantaSoft™ Rare Event Detection analysis. Event number refers to the number of droplet events generated for a given well or sample; Ch 1 amplitude measurement refers to the level of fluorescence emitted by a droplet event; and each column is a single well. Columns, or wells, are separated by yellow bars; Column D01 corresponds to one ddPCR replicate from the ambient Mobile Bay water sample and F01 corresponds to one ddPCR replicate from the C. leucas eDNA positive water sample. Columns B11, A12, and B12 correspond to one ddPCR replicate from each negative control incorporated and shows no contamination occurred during any stage of this experiment F I G U R E 5 Carcharhinus leucas mean eDNA concentrations (unit of measure) in a flow-through mesocosm detected using the Bio-Rad ® QX200™ Droplet Reader and QuantaSoft™ using a manual detection threshold of 3,000 amplitude with the Rare Event Detection analysis setting. Each time point sample was run in Droplet Digital™ PCR (ddPCR) replicates of five, and standard error bars were used to show the variation in concentration estimates across the five ddPCR replicates for each sample. The lower limit of detection, found to be at least 0.6 copies/μl in this study, is indicated by a dotted line Fundamental research on the accumulation, persistence, and degradation of elasmobranch eDNA is necessary to improve the interpretation of results in eDNA field research. Here, we have shown that after adding a shark into closed and flow-through systems, target eDNA was detectable within 30 min. In the flowthrough system, the initial spike in target eDNA that occurred between 0.5 and 1.0 hr could be due to initial stress experienced by the shark after being added to the mesocosm, causing it to expel more DNA (e.g., Barnes et al. 2014 with the target species present (e.g., Simpfendorfer et al. 2016) or collecting water samples from known habitats, but without visually confirming the presence of the target species (e.g., Weltz et al. 2017). Future studies should assess DNA accumulation over different timescales than presented here, as well as how altered flow rates, water conditions (pH, temperature), weather conditions (photoperiod, cloud cover), and number and size of target species impact the accumulation and persistence of elasmobranch eDNA in marine systems.

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

AUTH O R CO NTR I B UTI O N S
All authors contributed to the conception and design of the study, the acquisition and the interpretation of the data, and writing the manuscript.

DATA AVA I L A B I L I T Y S TAT E M E N T
The datasets generated during and/or analyzed during the current study are available from N.M. Phillips on reasonable request.