Inferring presence of the western toad (Anaxyrus boreas) species complex using environmental DNA

Western toads (species complex comprised of Anaxyrus boreas, A. canorus, A. exsul, and A. nelsoni) are widely distributed in the western United States but are declining, particularly in the southeastern extent of their range. The subspecies A. b. boreas is listed as a Species of Greatest Conservation Need in New Mexico, Colorado, Utah, and Wyoming. Reliable and sensitive methods for delineating distributions of western toads are critical for monitoring the status of the species and prioritizing conservation efforts. We developed two qPCR assays for detecting western toad DNA in environmental DNA samples. Both markers efficiently and reliably detect low concentrations of western toad DNA across their range in the conterminous U. S. without detecting non-target, sympatric species. To determine the optimal annual sampling period, we then tested these markers using repeated sampling in ponds where western toads were known to be present. Quantities of collected eDNA varied widely across samples, but sample-level detections across sites exceeded 80% for June sampling. In the later summer, detection dropped off sharply with only a single detection in the ten samples collected throughout August. © 2018 Published by Elsevier B.V. This is an open access article under the CC BY-NC-ND license (http://creativecommons.org/licenses/by-nc-nd/4.0/).


Introduction
Toads in the genus Anaxyrus represent a widely distributed and sometimes cryptic complex of species, subspecies, and lineages in North America (Frost et al., 2017). In the western United States, a number of forms are regarded as taxa of conservation concern because they are highly local endemics (Forrest et al., 2017;Gordon et al., 2017) or have undergone widespread declines leading to consideration for listing under the U.S. Endangered Species Act (U.S. Fish and Wildlife Service 2017). To prioritize conservation efforts for this species complex, reliable methods for assessing their presence and distributions are needed.
Environmental DNA (eDNA) has emerged as an efficient and useful tool for detecting rare or invasive aquatic species (Dejean et al., 2012;Wilcox et al., 2013;Sigsgaard et al., 2015) and delimiting distributions of rare species McKelvey et al., 2016). By integrating minor-groove-binding probes into quantitative PCR (qPCR), eDNA sample analysis has proven effective in detecting low concentrations of targeted DNA and can be more sensitive than traditional PCR methods (Kutyavin et al., 2000;Turner et al., 2014). Development of eDNA assays, however, can be problematic for species complexes that exhibit incomplete lineage sorting or minor levels of divergence. In such cases, it may only be feasible to develop assays that are specific and effective for groups of species within a genus.
Here, we describe the development of two qPCR assays to detect DNA of the western toad species complex. We follow Goebel et al. (2009;also see Pyron and Wiens, 2011) and use western toad to collectively refer to A. boreas, A. canorus, A. exsul, and A. nelsoni. Multiple independent assays can be useful in validating eDNA detections as they increase the likelihood of positive detections while reducing the probabilities of false positive detections (Carim et al., 2016a).

Assay development
To develop environmental DNA markers for western toad, we chose the cytochrome b (cytb) and cytochrome oxidase c subunit I (COI) regions of the mitochondrial genome because these regions provided the best combination of publically available reference samples for target and non-target species across a comprehensive geographic range. These regions were also used for in silico candidate primer and probe design because they provided sufficient nucleotide differences to distinguish western toads from most other species. After initially screening all publically available sequences for western toads and other closely related and sympatric toads, we found that sequences for all species within the western toad species complex were identical or nearly identical in these regions (Table A.1). Consequently, we pursued eDNA markers that would detect all species within the western toad species complex. We obtained candidate primers in silico using the DECIPHER package (Wright et al., 2014) in R v. 3.2.3 (R Core Development Team 2015, and aligned them with the genetic sequence data using MEGA 7.0 (Kumar et al., 2016). We then manually adjusted primer lengths and positions to optimize annealing temperatures and maximize base-pair mismatches with non-target species. Using the MEGA sequence alignments, we visually identified regions unique to western toad and designed TaqMan™ MGB probes (Applied Biosystems) with 6-carboxyfluorescein (FAM)labeled 5' ends and minor-groove-binding, non-fluorescent quenchers (MGB-NFQ) within the target amplicon sequence of each gene (Table A.2). We assessed annealing temperatures for each primer-probe set in Primer Express 3.0.1 (Life Technologies; Table A.2) and screened each marker for secondary structures using IDT OligoAnalyzer (https://www.idtdna.com/ calc/analyzer).
To test the specificity of each marker, we screened DNA extracted from 101 western toad tissues from 46 locations in the western U.S. (100 A. boreas, 1 A. canorus; Fig. 1), 9 non-target amphibian species, and 31 non-target fish species (Table A.3). All samples used in this study were from existing collections acquired under appropriate sampling permits. Toad tissues included toe clips, muscle tissues, and liver tissues from adult specimens and tail clips from tadpoles; fish tissues were fin clips. Tissues were extracted with the DNeasy Tissue and Blood Kit (Qiagen, Inc.) following the manufacturer's protocol. We screened each marker in vitro against tissue-derived DNA in a single qPCR reaction. Screening was performed on a StepOne Plus Real-time PCR Instrument (Life Technologies) in 15 ml reactions containing 7.5 ml Environmental Master Mix 2.0 (Life Technologies), 900 nM forward primer, 900 nM reverse primer, 250 nM probe, 4 ml DNA template (~0.12e0.88 ng), and 2.75 ml deionized water. Thermocycler conditions included 95 C for 10 min followed by 45 cycles of denaturation at 95 C for 15 s and annealing and extension at 60 C for 1 min. To minimize the risk of sample cross contamination, all qPCR tests were set up inside of a UV hood where consumables and pipettes were irradiated with UV light for 1 h prior to each test. Each test included a notemplate control with distilled water used in place of DNA template.
We optimized primer concentrations following methods outlined in Wilcox et al. (2015;Table A.2). Using the optimized concentrations and cycling conditions above, we tested sensitivity of the markers by performing standard curve experiments created from western toad (specifically, A. boreas) qPCR product. For each marker, qPCR product was purified using PureLink™ PCR Micro Kit (Invitrogen) and quantified on a Qubit 2.0 fluorometer (ThermoFisher Scientific). From this stock, we prepared a six-level standard curve dilution series (6 250, 1 250, 250, 50, 10, and 2 copies per 4 ml) in sterile TE. We ran six replicates of each dilution.
Finally, we applied both assays to eDNA samples collected from two sites in Utah with known presence of western toads based on visual encounter surveys (Lower Pond at Left Fork Kanab Creek and Snow Lake; Table 1) and a site in Montana believed to be occupied by amphibians other than western toads (Hidden Lake, 47.14311 N, À113.57065 W). For the assay validation samples, eDNA was collected from a 300-ml water sample at Snow Lake on 11 September 2016, a 1000-ml water sample at Left Fork Kanab Creek on 31 August 2016, and a 5000-ml water sample at Hidden Lake on 28 October 2015 using a peristaltic pump and a 1.5 mm glass mircrofiber filter as described in Carim et al. (2016b). The sampling method for these differed from Carim et al. (2016b) in that sampling concluded when a filter was clogged, preventing water passage. The eDNA samples were then extracted with the DNeasy Blood and Tissue Kit (Qiagen, Inc.) using a modified protocol (Carim et al., 2016c) in a room dedicated solely to this practice where the work area, scissors, and forceps were cleaned with 50% bleach prior to extractions. Negative extraction controls were paired with eDNA extraction sets and filtered pipette tips were used for all laboratory based methods. All extracts were stored at À20 C until analyzed. Environmental DNA samples and extraction controls were analyzed in 15 ml volumes using optimized concentrations (Table A.2) and the same PCR cycling conditions and recipe above except a TaqMan Exogenous Internal Positive Control (IPC; Life Technologies) consisting of 1.5 ml of 10X IPC assay and 0.30 ml 50X IPC DNA, was included in place of 1.8 ml distilled water to test for inhibition. All reactions were run in triplicate, along with triplicate no-template controls and positive controls for each analysis, on a QuantStudio 3 Realtime PCR System (Life Technologies). For all qPCR experiments, a reaction was considered positive if the amplification curve crossed the threshold during the exponential phase; a sample was considered positive if at least one reaction in the triplicate was positive. A sample was considered inhibited if there was a > 1 cycle-threshold (Ct) shift in the IPC relative to the notemplate control.

Field trials
Like most amphibians, western toads are exclusively aquatic from the egg stage to metamorphosis, after which juveniles and adults are often found in terrestrial habitats. To assess temporal variation in detection rates in potential breeding habitats, we sampled at roughly 2-week intervals from May to August 2017 at five sites in southern Utah collecting a single sample per visit. These included the two sites where samples were collected for in vivo assay tests and three additional sites: a stock pond on Dry Creek, the upper pond at Left Fork Kanab Creek, and a wetland 2.25 km north of Deep Creek Lake. We used a modified Carim et al. (2016b) sampling protocol to account for the wetland system in which we stopped filtering after 5000-mL or after three filters were clogged per a sample, whichever came first. If a sample required multiple filters, they were combined during extraction and analyzed as a single sample.
Because both assays detected western toad DNA with robust specificity and sensitivity during in vitro and in vivo validation (see results), we opted to analyze these samples using the COI assay. If any samples appeared inhibited, DNA was extracted from the second half of the filter and DNA from both filter halves was combined and cleaned by running through an inhibitor removal column (Zymo Research; www.zymoresearch.com) and then re-analyzed. To estimate the number of target DNA copies in each positive sample, samples were analyzed alongside the same six-level standard curve used for assay sensitivity analysis (above). To estimate copy number, we generated a regression line of the amplification Ct values on the known starting quantities of DNA in the standard curve. All reactions were set up inside of a hood which was irradiated with UV for at least 1 h prior to PCR set-up. Locations <5 km apart were combined (see Table A.3 for a list of specific sample locations). The hatching represents the range of the western toad species complex in the United States (derived from IUCN Red List of Threatened Species 2018). (For interpretation of the references to colour in this figure legend, the reader is referred to the Web version of this article.)

Assay development
The in silico design phase identified primers amplifying a 92-and 111-base-pair fragment in the cytb and COI genes, respectively. There was a minimum of nine mismatches between the cytb assay, and twelve mismatches between the COI assay and the sequences of non-target species (Table A.1); both probes are identical to the reference sequences from the western toad species complex that are of sufficient length to provide coverage of the probe hybridization region. Both markers detected DNA in vitro from all target samples and did not detect DNA from any non-target species (Table A.3). The standard curve analysis resulted in an amplification efficiency of 92.7% (standard curve y-intercept ¼ 39.383, r 2 ¼ 0.986) and 96.1% (standard curve y-intercept ¼ 39.674, r 2 ¼ 0.995) for the cytb and COI markers, respectively. The limit of detection (lowest concentration with >95% amplification success; Bustin et al., 2009) for each marker was 10 mitochondrial DNA copies/reaction (mtDNA copies/rxn) with successful detection of target DNA in all six replicates; however, each marker also detected target DNA in five out of six replicates at 2 mtDNA copies/rxn (83.3% detection success). Using both the cytb and COI markers, we detected western toad DNA in the eDNA samples from both locations known to be occupied in Utah, and did not detect this DNA in the presumed unoccupied location in Montana.

Field trials
A total of 34 samples were taken across five sites from May to August 2017 with filter volumes ranging from 350-to 5000-mL. Filter clogging was common and occurred in 30 of 34 samples. Inhibition was initially observed in 11 of 34 samples, but treatment successfully removed inhibition in all 11 samples. Detection rates varied by site. The Dry Creek stock pond was lowest (1/7) and the upper pond on Left Fork Kanab Creek was the highest (6/8; Table 1). Detection rates were generally greater in early summer and declined by late August (Fig. 2). Lastly, there were no detections of western toad DNA in the extraction blanks or the qPCR no-template controls using either assay.

Discussion
The eDNA markers described here are highly sensitive and reliably detect low concentrations of western toad DNA. Because the in vitro tests successfully detected western toads collected across the western conterminous U.S., we expect this assay to perform reliably within this region. Although the markers were only field validated at A. b. boreas locations, both markers are likely to detect all species of the western toad complex because their reference sequences are identical in our primer and probe regions. This is important to note because if sampling at locations where multiple species/subspecies of the complex co-occur, follow-up efforts involving physical surveys may be needed to resolve distribution boundaries of individual species/subspecies if that level of resolution is necessary. Nonetheless, given the often local distribution and uncertain divergence of particular lineages, we suggest screening assays and/or sequencing the cytb and COI genes of specimens of these other taxa to ensure that the assays will be effective for all individuals of the western toad complex. Although we could not obtain tissue from all congeneric non-target species, specifically Arizona toad (Anaxyrus microscaphus) and Great Plains toad (Anaxyrus cognatus), the COI assay boast at least 13 basepair mismatches with these species (Table A.1). Given the number and location of mismatches with the COI assay to both the Arizona toad and Great Plains toad, we are confident there would be no non-target amplification. This is supported by the lack of amplification for all of the non-target Anaxyrus tissues tested in which they all had comparable numbers of mismatches (see Table A.1 below). Unfortunately, there are no cytb sequences available for Arizona toad and Great Plains toad that overlap with the assay region. If analyzing eDNA samples at locations where Arizona toads and Great Plains toads could occur, we recommend applying the COI assay since the cytB assay could not be validated in vitro.
These assays have the potential to form the basis for eDNA-based surveys of western toads across their range, but must account for the issues encountered in sampling highly productive, often turbid, and sometimes ephemeral habitats for the presence of this species (Keinath and McGee, 2005). Filter clogging, which lengthens the sampling process and complicates eDNA extraction, as well as sample inhibition, which can potentially compromise results, are primary concerns when sampling in these habitats. Loss of DNA during inhibitor removal is possible (on average less than 10% loss in 100e200 ml elutions; see http://www.zymoresearch.com for more details), and could potentially result in false negative detections when target DNA is present at low quantities. However, this treatment process still provides the best results for inhibitor removal and successful analysis of inhibited samples (McKee et al., 2015). By adopting a modified protocol in which we filter 5000-mL or clog three filters, we are maximizing the amount of water, and therefore DNA, filtered while not overwhelming the extraction process with too many filters. Other studies have shown that increasing the pore size of the filter in turbid systems can reduce clogging, increase the volume of water filtered, and, therefore, increase the amount of DNA captured on the filter even though smaller particles of eDNA may be lost through the larger pore sizes (Turner et al., 2014;Goldberg et al., 2018). Combining these techniques with our fixed filtration protocols may enhance the between-sample consistency required to make reliable inferences about species presence.
Environmental DNA sampling has been used to infer seasonal activities of amphibians based on eDNA concentrations in both lotic and lentic systems Buxton et al., 2017). Similarly, our field sampling demonstrated the expected temporal variation of western toad presence in potential breeding habitats. Following breeding season in late spring, larval forms are expected to be present until metamorphosis in early to mid-summer, the timing of which varies with elevation of air and water temperature (Keinath and McGee, 2005). At that point, juvenile and adult western toads may only occasionally use standing water habitats, or may be more likely to frequent streams that promote movements across an individual's home range (Young and Schmetterling, 2009). Because western toads are extremely difficult to detect outside of the breeding season and during larval development, the ease of geographically distributed, repeat sampling for eDNA (McKelvey et al., 2016) suggest that this tool could prove valuable for evaluating seasonal occupancy of breeding, dispersal, and overwintering habitats. In this study, the only positive results associated with eDNA surveys in ponds lacking visual evidence of toads occurred in the lower site on Left Fork Kanab Creek, results which could be attributable to downstream drift of DNA. We did, however, obtain positive eDNA detections in the Dry Creek Stock Pond, where at most a single adult toad was observed, indicating that the method is capable of detecting toads at low densities if present. The likely drift of DNA into the lower site in Left Fork Kanab Creek indicates both a need for caution and an opportunity: caution because detections could indicate upstream toad presence and opportunity because sampling at downstream locations may allow efficient detection in complex wetlands such as are commonly produced by beaver dams. Given the high per-sample detection rates, abilities to detect toads a low densities, and to detect toads downstream from wetlands where they live, a robust detection protocol for this species appears feasible.

Declaration of interest
The authors declare no financial or other conflict of interests.

Funding
This study was funded by the USDA Forest Service Intermountain Region and the USDA Forest Service Intermountain Region BeSMART Microgrant.

Role of funding
The funding sources had no role in the study design; in the collection, analysis and interpretation of data; in the writing of the report; or in the decision to submit the article for publication. * Represents species included in the western toad complex (as described by Goebel et al., 2009) used as "targets" for assay development. ** The Dixie Valley toad is designated as Anaxyrus boreas in GenBank and is considered part of the A. boreas species complex as reported by Forrest et al. (2017). The species name has since been designated as Bufo (Anaxyrus) williamsoni as reported by Gordon et al. (2017).