Refining sampling protocols for cavefishes and cave crayfishes to account for environmental variation

Subterranean habitats represent focal habitats in many conservation strategies; however, these environments are some of the most difficult to sample. New sampling methods, such as environmental DNA (eDNA), show promise to improve stygobiont detection, but sources of sampling bias are poorly understood. Therefore, we determined the factors affecting detection probability using traditional visual surveys and eDNA surveys for both cavefishes and cave crayfishes and demonstrated how detection affects survey efforts for these taxa. We sampled 40 sites (179 visual and 183 eDNA surveys) across the Ozark Highlands ecoregion. We estimated the detection probability of cave crayfishes and cavefishes using both survey methods under varying environmental conditions. The effectiveness of eDNA or visual surveys varied by environmental conditions (i.e., water volume, prevailing substrate, and water velocity) and the target taxa. When sampling in areas with average water velocity, no flow, and coarse substrate, eDNA surveys had a higher detection probability (0.49) than visual surveys (0.35) for cavefishes and visual surveys (0.67) had a higher detection probability than eDNA surveys (0.40) for cave crayfishes. Under the same sampling conditions, 5 visual surveys compared to 10 eDNA surveys would be needed to confidently detect cave Subterranean Biology 39: 79–105 (2021) doi: 10.3897/subtbiol.39.64279 https://subtbiol.pensoft.net Copyright Joshua B. Mouser et al. This is an open access article distributed under the terms of the Creative Commons Attribution License (CC BY 4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited. RESEARCH ARTICLE Subterranean Biology Published by The International Society for Subterranean Biology A peer-reviewed open-access journal


Introduction
Variable species detection probability (i.e., the probability of detecting a species if present) is a fundamental sampling challenge when conducting ecological studies (Mac-Kenzie et al. 2018). Species detection can vary among habitats (Mollenhauer et al. 2018), sampling approaches (Pregler et al. 2015), species (McManamay et al. 2014;Mollenhauer et al. 2018), and over time (Hangsleben et al. 2013). The underlying species-environmental relationships of interest often do not emerge without consideration of variable sampling detection (Gwinn et al. 2016). Further, not accounting for variable detection can lead to incorrect estimates of extinction rates (Kéry et al. 2006;Pregler et al. 2015), species richness (Tingley and Beissinger 2013), and distributions (Chen et al. 2013;Lahoz-Monfort et al. 2014), largely due to false absences. For example, switching from seining to backpack electrofishing for sampling bridled shiner Notropis bifrenatus (Cope, 1867) in Connecticut led to underestimation of the species distribution due to differences in gear efficiencies (Pregler et al. 2015). As sampling design and statistics advances, ecologists had more options available to account for imperfect detection.
Variable species detection can be taken into account using appropriate study designs. Sampling standardization is useful for limiting some sampling variability (e.g., sampling at the same time of year; see also Bonar et al. 2009), but standardization alone does not account for environmental variability (i.e., flow or habitat) that is often of interest to ecologists (MacKenzie et al. 2004). For example, Mollenhauer et al. (2018) used standardized sampling to estimate the occupancy of Great Plains fishes and showed that sand shiner Notropis stramineus (Cope, 1865) occurrence was underestimated in one of two rivers due to differences in the environment (i.e., not controlled through standardization). Detection probability can be estimated by using a study design that uses repeated sampling while measuring environmental factors hypothesized to influence detection (MacKenzie et al. 2018). The concerns associated with sampling detection can be exacerbated in environments that are particularly difficult to sample or when species are rare.
Sampling difficulties in complex environments or where species are relatively rare create challenges for developing meaningful conservation actions. Large rivers, for example, are difficult to sample because of deep water, higher discharges (e.g., Detroit River, Lapointe et al. 2006), and vegetation cover (e.g., Niagara River, Crane and Kapuscinski 2018). Swampy streams, emblematic of complex habitat, make sampling using traditional approaches difficult (e.g., unconsolidated substrates, emergent vegetation, Jensen and Voukoun 2013). Even rivers with relatively homogenous substrates are difficult to sample when the target species is relatively rare (e.g., federally threatened Arkansas River Shiner Notropis girardi Hubbs & Ortenburger, 1929;Mollenhauer et al. 2018). In fact, many aquatic species are relatively rare or cryptic making adequate sampling problematic (e.g., bridle shiner, Jensen and Voukoun 2013; bull trout Salvelinus confluentus, Sepulveda et al. 2019). Despite advancements in sampling strategies, there remain notable examples of aquatic habitats that are difficult to sample but are considered especially important for both conservation and ecosystem services (e.g., subterranean environments, Mammola et al. 2019).
Sampling cavefishes and cave crayfishes can be difficult due to the challenges of traversing and sampling the subterranean environment. Cavefishes and cave crayfishes are typically surveyed by 1-3 people walking, crawling, or snorkeling slowly upstream in caves while recording the number of organisms observed (e.g., Graening et al. , 2010Bichuette and Trajano 2015;Behrmann-Godel et al. 2017). Stygobiotic organisms (groundwater obligates, Sket 2008) may go undetected during visual surveys due to similar environmental factors as surface aquatic environments (e.g., water depth, turbidity), but also because researchers can only access limited portions of the underground ecosystem and accessible areas are often difficult to traverse (Mammola et al. 2020). Additionally, the biology of cave organisms (e.g., low density due to k-selected life history and uneven distributions within caves) makes them difficult to detect (Mammola et al. 2020). Several cave surveys may be needed to detect stygobionts; thus, estimates of species occupancy and richness are skewed toward commonly sampled locations (Culver et al. 2004;Krejca and Weckerly 2008).
Sampling using environmental DNA (eDNA) is a relatively new technique in ecology and conservation biology that may improve detection of cavefishes and cave crayfishes; however, sources of variable detection probability are poorly understood. Environmental DNA surveys document species presence via the collection of DNA from the environment (Ficetola et al. 2008), which is derived from sources such as waste products, shed hair and skin, the slime coat of fishes and amphibians, shed exoskeletons of arthropods, and decomposing individuals (Tréguier et al. 2014;Thomsen and Willerslev 2015). Many taxa have been surveyed via eDNA in surface habitats, including fishes (e.g., Jerde et al. 2011), crayfishes (e.g., Tréguier et al. 2014, mollusks (e.g., Egan et al. 2013), and reptiles (e.g., Piaggio et al. 2014). Studies in subterranean habitats are few to date (reviewed in Gorički 2019) but include stygobiotic Proteus salamanders (Gorički et al. 2017;Vörös et al. 2017), Stygobromus amphipods (Niemiller et al. 2018), and two Cambarus species of cave crayfishes (Boyd 2019). Environmental DNA surveys can improve species detection when compared to traditional survey methods (Jerde et al. 2011;Smart et al. 2015;Schmelzle and Kinziger 2016). Further, eDNA surveys make it possible to survey karst environments without sampling entire caves.
Understanding how our sampling approaches relate to our ability to detect a species is important to developing meaningful conservation actions. In many cases, particularly with rare or cryptic organisms, sampling results in false absences (i.e., species was undetected when present). Therefore, our study objective was to determine some of the environmental factors associated with detection probability of cave crayfishes and cavefishes using both visual and eDNA surveys. Our overall goal was to assess how sampling bias related to the effort needed to adequately sample these taxa and obtain reliable presence or absence inferences. Results of this study will help managers choose the most efficient sampling approach for determining the presence of cavefishes and cave crayfishes and understand sources of detection error for both eDNA and visual surveys.

Study area
We conducted our study in the Ozark Highlands level-three ecoregion (hereafter referred to as the Ozark Highlands) of northeast Oklahoma, southwest Missouri, and northwest Arkansas (Figure 1). Average annual rainfall and air temperatures of the Ozark Highlands are 116 cm and 13.7 °C, respectively (30 yr climate normal for Springfield, Missouri; National Oceanic and Atmospheric Administration). The ecoregion was historically a mix of prairie, oak, hickory, and pine forests, but many lowland areas have been converted to agricultural uses (Woods et al. 2005). The lithology of the Ozark Highlands is primarily Mississippian limestone and Ordovician dolomite, which have been dissolved over time by groundwater, resulting in thousands of caves and springs (Unklesbay and Vineyard 1992).

Study species
We focused our study on two species of cavefishes, Ozark cavefish Troglichthys rosae (Eigenmann, 1898) and Eigenmann's cavefish Typhlichthys eigenmanni (Girard, 1859), and 5 species of cave crayfishes, Benton cave crayfish Cambarus aculabrum (Hobbs & Brown, 1987), bristly cave crayfish C. setosus (Faxon, 1889), Delaware county cave crayfish C. subterraneus (Hobbs, 1993), Oklahoma cave crayfish C. tartarus (Hobbs & Cooper, 1972), and Caney Mountain cave crayfish Orconectes stygocaneyi (Hobbs, 2001). The full distributions of many of our target species are unknown, though existing sampling data provide some insight. There is no known overlap in distributions among species within taxa (i.e., cave crayfishes or cavefishes; Figure 1). Troglichthys rosae is assumed to occur in the Springfield Plateau of northwest Arkansas, southwest Missouri, and northeast Oklahoma (Niemiller and Poulson 2010). Typhlichthys eigenmanni is considered endemic to the Ozark Highlands of central and southeast Missouri and northeast Arkansas (Niemiller et al. 2012); however, we only sampled along the western portion of the species estimated range. Cambarus aculabrum is known from only four locations in northwest Arkansas . Cambarus setosus is the widest-ranging cave crayfish of the Ozark Highlands and has been docu-Figure 1. We conducted eDNA and visual surveys for cavefishes and cave crayfishes at 40 caves, wells, and springs across the Ozarks Highlands ecoregion. The estimated ranges of the cave crayfishes (i.e., Cambarus aculabrum, C. setosus, C. subterraneus, C. tartarus, Orconectes stygocaneyi) are outlined using United States Geological Survey 12-digit watersheds that encompass locations where these species have been observed. Troglichthys rosae is thought to restricted to the Springfield Plateau (light-grey outline); however, we detected the species at some of the sites enclosed by the circle. Typhlichthys eigenmanni was only surveyed at our two northern-most sites. mented at 48 sites in southwest Missouri and two sites in Arkansas (Graening et al. 2006b). Cambarus subterraneus and C. tartarus have only been found in three and two caves in northeast Oklahoma, respectively (Graening and Fenolio 2005;Graening et al. 2006c). Orconectes stygocaneyi is assumed to be endemic to a single cave in southcentral Missouri (Hobbs III 2001). Little is known about the biology and ecology of these organisms; however, descriptions for each species are quite similar due to convergent evolution (e.g., cryptic behavior, habitat generalists, albinistic, and reduced eyes).

Study design
We conducted both eDNA and visual surveys for cavefishes and cave crayfishes at 21 caves, 12 springs, and 7 wells (Figure 1, Suppl. material 1: Table S1). We sampled caves, springs, and wells (hereafter referred to as sites) because they allow access to the groundwater habitat occupied by stygobionts. In fact, state agencies routinely sample hand-dug wells because they offer access to groundwater where caves may not be accessible and it is common to locate stygobionts at those locations (Doug Novinger, personal communication). We chose a combination of sites where some had previous documentation of cave crayfish and (or) cavefish occupancy (n = 24) and others had either never been sampled or no cavefishes or crayfishes had ever been identified (n = 16). Sites 15 and 16 occurred in the same cave but were considered different sites due to extreme differences in the hydrologic regime (Miller 2010). We selected 1-5 sampling units (n = 61) at each site based on presumed biological barriers (e.g., waterfalls or disconnected pools) and no sampling units were adjacent to one another. We chose to select multiple sampling units within caves because this allowed us to assess the spatial distribution of stygobionts. Sampling units were referenced by the site number, and then the sequential number of units within the site (e.g., 1.2 referred to the second sampling unit within site 1). For example, sampling unit 10.1 was a cave with a single pool of water and no discernible change in habitat, and sampling unit 16.1 was a pool within a cave bounded by a waterfall downstream and shallow riffle upstream. Sampling units were surveyed on 1-5 occasions (179 visual surveys and 183 eDNA surveys; Suppl. material 1: Table S2). Sampling was conducted during a relatively short time period (February-May 2017) to meet a closed-system assumption with respect to species occurrence (i.e., species neither colonize the sampling units nor go extinct during the survey period). Although there was some typical spring flooding at the end of our sampling period, we assumed there would be a lag between the initiation of high-water or low-water events before changes in species occupancy would occur (i.e., it would take time for species to recolonize when a sampling unit either became wet or dry again, Adams and Warren 2005). Further, defining our season to allow some changes in the physicochemical parameters at each sampling unit (Suppl. material 1: Table S1) was preferred to examine relationships between detection and a range of physicochemical parameters using both sampling methods.

eDNA surveys
We collected two water samples (≈ 1-L each) for eDNA analysis at each sampling unit during each visit. We collected two water samples to provide a replicate in case of error or contamination in subsequent steps. We immersed sampling equipment in 50% bleach for at least 30 s between sites and then rinsed it in deionized water to avoid contamination. If possible, we sterilized gear between sampling units, but some caves were too difficult to navigate with more than a single equipment set. We filtered distilled water in the field on four occasions to provide negative controls, which were treated the same as field samples in subsequent steps. Water was collected in two 1-L sample bottles (312187-0032, ThermoFisher Scientific, Waltham, Massachusetts) from approximately 10 cm above the substrate, where possible, without disturbing the substrate. Water was collected just above the substrate when water depth was < 10 cm. We did not sample the substrate to both avoid inhibitors (e.g., humic acid) and possibly sampling older DNA within the substrates that was not indicative of current occupancy. To collect water from wells, we lowered a Van Dorn sampler (3-1920-H62, Wildco, Yulee, Florida) to approximately 10 cm above the substrate, closed the sampler, returned it to the surface, and transferred the water to two 1-L sample bottles. We filtered the water immediately after collection, except for the samples collected from sampling units 4.1-4.4 on 21 March 2017, which were frozen and filtered later in the laboratory. While wearing nitrile gloves, we placed a 0.45-µm cellulose-nitrate filter (14-555-624, Fisher Scientific, Waltham, Massachusetts) inside a filter funnel (09745, Fisher Scientific, Waltham, Massachusetts) attached to a vacuum flask via a rubber stopper ( Figure 2). We used a hand pump (AC3310, Advance Auto Parts, Raleigh, North Carolina) to create a vacuum to pull water through the filter. Only one filter was typically needed to sample one L of water, but occasionally multiple filters (i.e., 2-6) were used due to clogging via sediment. Filters were stored at room temperature in vials of 900 µl of Longmire's buffer (Longmire et al. 1997), until extractions were completed (i.e., 1-18 mo after collection).

Visual surveys
Visual surveys for cavefishes and cave crayfishes occurred at most of the sampling units for later comparison to eDNA detection. We did not complete visual surveys at sampling units 10.1 and 18.1 on the last two survey dates due to local flooding. We did not visually survey the entirety of sampling units 5.1 and 6.2 due to sampling restrictions by the regulatory agency (i.e., safety concerns or concern for trampling crayfish). For springs and caves, two observers walked or crawled the entire sampling unit while carefully searching the whole wetted area for cave crayfishes or cavefishes by overturning rocks and examining crevices using headlamps to illuminate dark areas (e.g., Graening et al. , 2010. Hand-dug wells were surveyed in their entirety using a spotlight (QBeam Max Million III, The Brinkmann Corporation, Dallas, Texas) both before and after water samples were collected because disturbance from sampling sometimes caused stygobionts to emerge. We recorded the number of cavefishes and cave crayfishes observed and time spent observing (min).

Detection covariates
Our detection covariates were chosen based on a priori knowledge derived from the literature. We hypothesized that increased water turbidity (Thurow et al. 2006), greater water volume (Trajano 2001), flowing water (Thurow et al. 2006), and substrate (coarse or fine) (Albanese et al. 2011) would make it more difficult to detect stygobionts via visual surveys. The presence of light indicates surface connection and may affect detection via altered species abundance due to food availability (Simon et al. 2003) or predator abundance (Brown and Todd 1987). Increasing species abundance generally results in greater detection probability for various sampling methods (e.g., Pregler et al. 2015;Baldigo et al. 2017). For eDNA surveys, we hypothesized increased water turbidity would relate to more inhibitors in our samples (e.g., humic While wearing gloves, a 0.45-µm microbial filter was placed inside a filter funnel that was attached to a vacuum flask via a rubber stopper. A hand pump was used to create a vacuum and pull water through the filter. Filters were stored at room temperature in vials of 900 µl of Longmire's buffer (Longmire et al. 1997). acid, Jane et al. 2015), ultraviolet light could breakdown eDNA (Strickler et al. 2015), increased water volume would dilute eDNA (Rice et al. 2018), faster water would expel eDNA from the sampling unit (Jane et al. 2015), and fine substrates could easily be resuspended and lead to inhibition of eDNA PCR amplification (e.g., Buxton et al. 2017); all of which would decrease detection using eDNA surveys.
We estimated or measured (numbers in parentheses represent our measurement resolution): water turbidity (0.01 NTU), light (present or absent), water volume (1.0 m 3 ), water-column velocity (hereafter water velocity, 0.01 m/s), and substrate (coded as either coarse; or fine or bedrock) at each sampling unit to explain variable detection of cave biota. We collected 250-ml water samples before the start of each visual survey to measure water turbidity using a turbidity meter (AQUAfast AQ4500, Thermo Fisher Scientific, Waltham, Massachusetts). Light was recorded as ambient light visible (present) or not visible (absent) at the water-sample location. The water volume of each sampling unit was estimated by multiplying survey length (1.0 m), wetted width (0.1 m), and maximum water depth (0.1 m). Wetted width and maximum water depth were measured at 3-5 points along the sampling unit to represent average conditions. Water velocity was visually estimated at the same locations where we measured wetted width and maximum water depth. We visually estimated water velocity because it was unreasonable to bring a flow meter into many of the caves we sampled (e.g., narrow crawl spaces and deep water). Prior to the study, we compared our visual water velocity estimates to values measured with a Marsh-McBirney flow meter (Marsh-McBirney Inc., Frederick, Maryland) to ensure that our estimates were relatively accurate (i.e., ± 0.1 m/s). We also distinguished between the prevalence of clay, silt, or bedrock substrates (hereafter "fine"), or pebble substrate, cobble substrates, or woody debris of similar size or larger (hereafter "coarse") at each sampling unit (see Wentworth 1922 for sizes of each substrate). Substrate was combined into these two categories based on our ability to observe stygobionts in these habitats. Stygobionts are relatively easy to observe on clay and bedrock substrates because they cannot conceal themselves within either as they can in cobble or woody debris.

Primer and probe development
Primers and probes were designed to amplify DNA for each of our study species (i.e., a species-specific quantitative PCR [qPCR] Taqman assay). We acquired template DNA for each of our study species from various sources (Suppl. material 1: Table S3). Genomic DNA was extracted from tissue samples using the Qiagen DNeasy Blood and Tissue kit (69504, Qiagen, Hilden, Germany) according to the manufacturer's instructions. For cavefishes, a 500-bp fragment of the mitochondrial NADH dehydrogenase 2 (ND2) gene was PCR amplified using the forward primer MET: 5'-CAT-ACCCCAAACATGTTGGT-3' and reverse primer ND2B: 5'-TGGTTTAATCCGC-CTCAGCC-3' (Kocher et al. 1995). Each amplification reaction had a total volume of 30 µl, consisting of 1.0 µl of DNA, 2.4 µl of MgCl 2 (25mM), 4.8 µl of deoxynucleoside triphosphates (1 mM), 0.5 µl of forward primer (10 µM), 0.5 µl of reverse primer (10 µM), 2.4 µl of bovine serum albumin, 6.0 µl of GoTaq buffer, 0.12 µl of GoTaq DNA polymerase (M3001, Promega, Madison, Wisconsin), and 12.8 µl ddH 2 O. The thermal profile consisted of an initial denaturation step of 94 °C for 5 min followed by 35 cycles of denaturation at 94 °C for 30 s, annealing at 56 °C for 30 s, and extension at 72 °C for 30 s and elongation at 72 °C for 5 min. For cave crayfishes, a 710-bp frag-ment of the mitochondrial cytochrome c oxidase I (CO1) gene was amplified using the forward primer LCO1490: 5'-GGTCAACAAATCATAAAGATATTGG-3' and the reverse primer HCO2198: 5'-TAAACTTCAGGGTGACCAAAAAATCA-3' (Folmer et al. 1994). The amplification reaction consisted of the same reagents; however, the thermal profile was: an initial denaturation step of 94 °C for 5 min; 6 cycles of denaturation at 94 °C for 1 min, annealing at 50 °C for 1 min, and extension at 72 °C for 1.5 min; 35 cycles of denaturation at 94 °C for 1 min, annealing at 55 °C for 1 min, and extension at 72 °C for 1.5 min; and elongation at 72 °C for 5 min. We chose the CO1 and ND2 genes because they have a high copy number and are relatively easy to isolate and purify (Billington 2003), have rates of divergence that allow species to be distinguished (Billington 2003), and are commonly used to amplify DNA of cave crayfishes (e.g., Buhay et al. 2007) and cavefishes (e.g., Niemiller et al. 2012Niemiller et al. , 2013, respectively. PCR products were visualized on a 1.0% agarose gel then purified using a Wizard SV Gel and PCR Clean-Up System (A9281, Promega, Madison, Wisconsin). PCR products were Sanger sequenced and sequences were manually trimmed and aligned in Geneious (Version 11.1.5, Auckland, New Zealand) to generate a consensus sequence for the CO1 locus of each cave crayfish species and the ND2 locus of each cavefish species. The consensus sequences were entered in PrimerQuest (https://www.idtdna.com/primerquest/home/index) to generate species-specific qPCR Taqman assays (Table 1). Initial specificity of both the primers and probe was checked using Primer-Blast (https://www. ncbi.nlm.nih.gov/tools/primer-blast/) searching against the nr database on GenBank.
We performed in vitro validation and quantified the lower limit of detection for our assays. The lower limits of detection for C. setosus, C. subterraneus, C. tartarus, O. stygocaneyi, and T. rosae DNA were 1.5 × 10 -3 ng/µl, 3.9 × 10 -4 ng/µl, 1.5 × 10 -4 ng/ µl, 3.3 × 10 -4 ng/µl, and 2.5 × 10 -4 ng/µl, respectively. We were unable to test the assays in vitro for C. aculabrum and T. eigenmanni because we did not have genomic DNA for those species. We were unable to obtain samples of C. aculabrum DNA due to its rarity. We did not obtain samples of T. eigenmanni DNA because we only sampled a few sites, and many sequences were already available online. Not all assays developed were species-specific, but we confirmed species identity of field samples via Sanger sequencing of a subset of the positive samples.

eDNA extraction
We extracted eDNA from the filters using a DNeasy Blood and Tissue Kit by following the "purification of total DNA from crude lysates" protocol (Qiagen 2006) with the following modifications. We sterilized all laboratory surfaces and equipment with 10% bleach before extractions. DNA was initially extracted for only one filter collected at a sampling unit. Any additional filters were placed in fresh Longmire's buffer and set aside to use if the first filter was negative for the target species' DNA (see next section Quantitative PCR amplification). Using forceps, each filter was halved and torn into pieces. The pieces from each half were added to separate 2-ml microcentrifuge tubes. Forceps were sterilized between filters by immersion into 100% ethanol and Table 1. Taqman assays were designed to amplify DNA for each of our target species. The 5' end of the probe was labeled with the fluorescent dye (6-FAM), the 3' primer end with a quencher (Iowa Black ™ FQ), and there was an additional internal quencher (ZEN ™ ). Probes were doubled quenched to reduce background fluorescence and increase signal intensity. All primer and probe sequences are reported 5' to 3'.

Species
Forward flaming. The Longmire's buffer was split into two 1-ml tubes, and if the volume was < 360 µl, fresh buffer was added. The tubes of Longmire's buffer were then centrifuged at 8,000 g for 30 s. We then transferred 360 µl of the Longmire's buffer and the pellet to the respective tubes with the filter pieces. The above process resulted in a standard amount of filter pieces and buffer in each tube (i.e., 1 filter half and 360 µl). There were two tubes per sampling unit, and each tube was considered a subsample for that sampling unit. After samples were standardized, we followed the "purification of total DNA from crude lysates" protocol except we doubled the amounts of proteinase K, buffer AL, and 100% ethanol that were added to each tube. Further, we reduced the amount of buffer AE added in the final step to 125 µl. We stored our samples at 2 °C until amplification (i.e., up to 4 mo).

Quantitative PCR amplification
We amplified eDNA using quantitative Polymerase Chain Reaction (qPCR). Each amplification reaction had a total volume of 20 µl, consisting of 10 µl of TaqMan Environmental Master Mix 2.0 (4396838, ThermoFisher Scientific, Waltham, Massachusetts), 4.7 µl of ddH 2 O, 0.9 µl of forward primer (20 µM), 0.9 µl of reverse primer (20 µM), 0.5 µl of probe (10 µM), and 3.0 µl of template DNA. Samples were run in 96-well optical plates (BC3496, ThermoFisher Scientific, Waltham, Massachusetts) on a LightCycler 480 (Roche, Pleasanton, California). The thermal profile consisted of an initial denaturation step of 95 °C for 10 min followed by 40 cycles of denaturation at 95 °C for 15 s and annealing/extension at 60 °C for 1 min. Each subsample was run in triplicate, which resulted in an initial six pseudoreplicates for each sampling unit. If any pseudoreplicates amplified, then the sampling unit was considered positive for the species. If none of the pseodoreplicates amplified, we extracted eDNA from any re-maining filters from the sampling unit and ran another qPCR. We repeated the above process until all filters were processed, or until any pseudoreplicates amplified. If only one subsample amplified from a single survey date, then we processed that subsample again. If the subsample still amplified, then the survey was considered positive for the species and if it was negative, then the survey was considered negative. We also ran three negative controls during each qPCR in which the template DNA was replaced by ddH 2 O. If any of the negative controls amplified, then the qPCR run was discarded. A positive control was included that consisted of genomic DNA from the target taxa to ensure the reaction worked properly. We confirmed species identification of a subset of positive samples for each species using Sanger sequencing.

Statistical analysis
We modeled cave crayfish and cavefish detection probability when using both eDNA and visual surveys. Species of cave crayfishes and cavefishes tend to have narrow distributions (i.e., the species do not occur at all sites); thus, we modeled detection probability of all species of cave crayfishes as a single taxon and all species of cavefishes as a single taxon. Each taxon was either detected (1) or not detected (0) during each survey at a sampling unit, and the surveys were combined for sampling units to create capture histories for our response variable (i.e., a binomial response variable). For example, a capture history of 1010 would represent a taxon that was detected on the first and third surveys and undetected on the second and fourth surveys. Sampling units were included in the model twice if both taxa were detected, once if only one taxon was detected, and excluded if neither taxon was detected. Our approach required meeting the same assumptions for occupancy modeling with respect to the detection process: no false positives, sampling unit closure, and independent surveys. An alternative approach would be to make the individual surveys the outcome (i.e., 1 or 0, logistic regression), but we chose to use capture histories because it allowed us to evaluate model fit (see below). We assumed trait differences (e.g., morphology and behavior) among cavefish and cave crayfish species (i.e., within each taxon) would not influence detection probability. We excluded eDNA surveys for C. setosus because the assays did not amplify the subset of the field samples we tested with positive visual identification of the species. Our final model included 35 sampling units for cavefishes (105 visual surveys, 109 eDNA surveys) and 25 sampling units for cave crayfishes (77 visual surveys, 40 eDNA surveys). We modeled detection probability of cavefishes and cave crayfishes in relation to light, substrate, water volume, water velocity, and water turbidity, where each environmental variable varied by both taxa and sampling method. The continuous variables water volume and water turbidity were natural-log transformed due to rightskewed distributions. Water volume and water turbidity were standardized to a mean of zero and a variance of one to improve coefficient interpretation. The correlation level between water turbidity and water volume was low, indicating independence of these variables (Pearson's pairwise correlation coefficient = 0.11). We made velocity a category where 0 indicated no flow and 1 indicated flowing water. Light, substrate, and water velocity were treated as factors using a dummy variable approach (i.e., ambient light, no flow, and coarse substrate as the references). Independence between continuous and categorical variables was checked using point-biserial correlations and none were >0.22. Independence between categorical variables was checked using Cramer's V and none were >0.47. We also treated sampling method and taxa as factors using visual surveys and cave crayfish as the reference categories, respectively. The most complex model can be written as: for i = 1,2,…N, for j = 1,2,…J, where p ij is detection probability for survey j at sampling unit i, β 0 is the intercept, β 1 is the taxa main effect coefficient, β 2 is the method main effect coefficient, β 3 is the light main effect coefficient, β 4 is the turbidity main effect coefficient, β 5 is the velocity main effect coefficient, β 6 is the substrate main effect, β 7 is the volume main effect coefficient, β 8 is the taxa * method interaction term coefficient, β 9 is the taxa * light interaction term coefficient, β 10 is the taxa * turbidity interaction term coefficient, β 11 is the taxa * velocity interaction term coefficient, β 12 is the taxa * substrate interaction term coefficient, β 13 is the taxa * volume interaction term coefficient, β 14 is the method * light interaction term coefficient, β 15 is the method * turbidity interaction term coefficient, β 16 is the method * velocity interaction term coefficient, β 17 is the method * substrate interaction term coefficient, β 18 is the method * volume interaction term coefficient, β 19 is the taxa * method * light interaction term coefficient, β 20 is the taxa * method * turbidity interaction term coefficient, β 21 is the taxa * method * velocity interaction term coefficient, β 22 is the taxa * method * substrate interaction term coefficient, β 23 is the taxa * method * volume interaction term coefficient, X 1 is taxa, X 2 is method, X 3 is light, X 4 is turbidity, X 5 is velocity, X 6 is substrate, and X 7 is volume.
We fit our models using the program JAGS (Plummer 2003) called from the statistical software R (version 3.5.3; R Core Team 2019) using the package jagsUI (Kellner 2019). We used a broad uniform prior on the 0 to 1 scale for the detection probability intercept and broad uniform priors on the logit scale for other coefficients (Kéry and Royle 2016). Posterior distributions for coefficients were estimated using Markov chain Monte Carlo methods using 3 chains of 50,000 iterations each after a 10,000-iteration burn-in phase. We assessed convergence using the Brooks-Gelman-Rubin statistic (R , Gelman and Rubin 1992), where values < 1.1 for all model parameters indicates adequate mixing of chains (Kruschke 2015;Kellner 2019).
We used a three-step process to simplify our final model. We began by fitting the full model and simultaneously removed all three-way interaction terms with 95% highest density intervals (hereafter HDIs, Kruschke and Liddell 2018) that overlapped zero (i.e., were considered non-significant). The intervals are not interpreted in a traditional Frequentist sense (i.e., a 95% probability of containing the true value). Rather, the mean for the coefficient is the most plausible value, and the HDI contains credible values from the posterior distribution with a total probability of 95%. This use of a decision rule cut-off is analogous to hypothesis testing. However, an HDI that contains zero is not interpreted as failing to reject the null, but rather that an effect size of zero meets the minimum level of credibility. We then refit the model and used the aforementioned criteria to remove non-significant two-way interactions that were not retained in the three-way interactions. Finally, we repeated the above process to remove the main effects for environmental variables that were not significant. A modelselection process using HDIs has also been employed in similar studies (e.g., Kanno et al. 2015;Mihaljevic et al. 2015;White et al. 2020).
We examined model fit using posterior predictive distributions. The fit of the final model was assessed using a Bayesian p-value (Kéry and Royle 2016). A Bayesian pvalue closer to 0.5 suggests adequate fit, and extreme values (i.e., > 0.90 or < 0.10) indicate a lack of fit (Hobbs and Hooton 2015;Kéry and Royle 2016;Conn et al. 2018).
Because our goal was to assess how sampling bias related to the effort needed to adequately sample these taxa, we interpreted our results via cumulative detection plots. Cumulative detection probability (p c ) was calculated as: p c = (1 -(1 -p) k , where k is the number of surveys. We plotted the cumulative detection probability of each taxa for each significant relationship with method and environmental covariate.

eDNA and visual surveys
Environmental DNA surveys detected cavefishes at more sampling units than visual surveys, whereas visual surveys detected cave crayfishes at more sampling units compared to eDNA surveys. Environmental DNA surveys detected cavefishes at 33 of 61 sampling units, and visual surveys detected cavefishes at 14 of 61 sampling units. At 21 sampling units, we detected cavefish DNA but did not visually observe cavefishes. We detected cavefishes at six sites where they have never been detected using eDNA surveys, but did not detect any new populations using visual surveys. Environmental DNA surveys detected cave crayfishes at 10 of 61 sampling units, whereas visual surveys detected cave crayfishes at 17 of 61 sampling units. We detected cave crayfishes at one site where they have never been detected using eDNA surveys, but did not detect any new populations using visual surveys. Low eDNA detection could be the result of pseudogenes that we observed in the DNA of C. setosus and O. stygocaneyi. All of the negative controls collected in the field were negative, suggesting our decontamination protocol was adequate.

Detection covariates
The environmental factors we measured varied over the sample season (Suppl. material 1: Table S1). Water turbidity ranged from 0.20 to 41.50 NTU (mean ± SD = 2.97 ± 4.57 NTU). There was visible light at 26 sampling units, 34 sampling units were dark, and sampling unit 10.1 did not have visible light on the first 2 surveys but did on the last survey (i.e., we sampled at the cave entrance due to high water). We surveyed a range of water volumes across sampling units (0.06 m 3 -800.00 m 3 ; mean ± SD = 61.21 ± 132.00 m 3 ). Estimated water velocity ranged 0-0.53 m/s (mean ± SD = 0.06 ± 0.10 m/s), with 81 surveys classified as 0 (i.e., not flowing) and 102 surveys classified as 1 (flowing water). Substrate at 34 sampling units was classified as coarse substrate and 27 as fine substrate.

Statistical analysis
Detection probability of both cavefishes and cave crayfishes varied by survey method and was significantly related to water volume, substrate, and water velocity (Table 2). For cavefishes, detection probability at mean or reference levels of predictor variables was 0.35 (95% HDI: 0.19-0.55) using visual surveys and 0.49 (95% HDI: 0.32-0.67) using eDNA surveys. For cave crayfishes, detection probability at mean or reference levels of predictor variables was 0.67 (95% HDI: 0.47-0.84) using visual surveys and 0.40 (95% HDI: 0.19-0.65) using eDNA surveys. Cave crayfish and cavefish detection decreased sharply with increasing water volume using visual surveys. Cavefish detection decreased significantly when using visual surveys in sampling units classified by coarse rather than fine substrates. In contrast, cave crayfish detection decreased in when using eDNA surveys in fine compared to coarse substrates. Lastly, detection probability of cavefishes using visual surveys decreased significantly when water was flowing (i.e., water velocity > 0). R was < 1.1 for all coefficients. The calculated Bayesian p-value was 0.36.
The number of surveys needed to be confident the taxa were detected if present depended on sampling method and underlying environmental conditions. At mean levels of all predictor variables, approximately four eDNA surveys and nine visual surveys would be necessary to achieve a cumulative detection probability near one for cavefishes (i.e., confident the taxon was truly absent if undetected, Figure 3a). Alternatively, it would take approximately 5 visual surveys versus 10 eDNA surveys to achieve a cumulative detection probability near 1 for cave crayfishes when sampling under reference conditions (Figure 3b). When sampling in higher water volumes, greater than 10 visual surveys would be needed to confidently detect both cave crayfishes and cavefishes compared to less than 6 surveys in lower volume. Visually sampling for cavefishes at sites with fine substrates would require only four surveys to be confident of detection, whereas 10 surveys would be needed if the substrate was coarse. Seven surveys would be needed to confidently detect cavefish via eDNA surveys in both coarse and fine substrates. If we used eDNA sampling for cave crayfishes, then we would need 10 surveys in coarse substrate to be confident of detection versus more than 10 surveys in fine substrates. If we used visual surveys for cave crayfishes, then only five surveys  (1 -p) k , where p is detection probability at mean and reference levels (i.e., visual surveys, water not flowing, cave crayfishes, and coarse substrate) of predictor variables and k is the number of surveys. Table 2. Detection probability estimates from the final model for cavefishes and cave crayfishes using environmental DNA (eDNA) and visual surveys. Estimates for each parameter included in the detection model are reported on the logit scale as the mean ± standard deviation (SD) with a 95% high density interval (HDI). Mean values are reported as detection probabilities (Prob) by completing a logit transformation. The reference categories for categorical variables were visual surveys, water not flowing, cave crayfishes, and coarse substrate. would be needed in both coarse and fine substrates to achieve a detection probability near one. When water is flowing, it would take 5 visual surveys for cave crayfishes and > 10 surveys for cavefishes to achieve a cumulative detection probability near 1, whereas it would only take 4 surveys to detect both taxa using eDNA.

Discussion
We show detection probabilities for both cavefishes and cave crayfishes depend on both the sampling environment and method. Several studies have demonstrated that detection can be extremely low (< 0.01-0.18) for cave organisms when using visual surveys (Culver et al. 2004;Krejca and Weckerly 2008). We found that detection can be low for stygobionts and that it would take at least nine visual surveys to ensure cavefish detection with traditional visual survey methods. With relatively low detection probabilities for cavefishes, it would be possible to conclude that a cave is unoccupied when caves are often surveyed less than once per year in the Ozark Highlands (e.g., Graening et al. 2010). Because of the relatively low detection, managers would benefit from considering study designs that account for detection (i.e., accounting for the sampling bias). If repeat surveys are not possible, another option is to use the most efficient survey method for the target taxa under the prevailing environmental conditions realizing that underestimating occupancy would be likely. Detection probability via eDNA surveys can depend on the target species and its associated density. We observed that detection using eDNA surveys was typically higher for cavefishes than for cave crayfishes. Although some of the discrepancy in detection between cavefishes and cave crayfishes can be explained by the availability of genetic data, physiological differences may also play a role. For example, fish have a slime coat and release more DNA in the environment than crayfish that have a hard exoskeleton (Tréguier et al. 2014), thus making it easier to detect fishes. The abundance or biomass of the target organism also relates to how much DNA will be released into the environment (Takahara et al. 2012) and can influence detection (Dougherty et al. 2016;Baldigo et al. 2017). For example, we were unable to detect T. eigenmanni at sampling units where only one fish was observed across all surveys, but we detected them at sampling units where multiple individuals were observed. Other studies, however, have observed little relationship between target organism and detection (Rice et al. 2018). Water volume may interact with species abundance to further influence detection because eDNA may be diluted when there is more water. For example, we observed decreased detection with increased water volume.
The movement and persistence of eDNA in the environment can further complicate detection of aquatic organisms. In surface waters, eDNA flows downstream (up to 12.3 km; Deiner and Altermatt 2014) and can settle vertically (Turner et al. 2015). For example, Asian carp DNA was detected upstream of a fish barrier near the Great Lakes (Jerde et al. 2011), but flow reversals, not presence, were provided as the explanation (Song et al. 2017). In karst environments, water can flow in many directions due to gravity and topography (Aley and Kirkland 2012), which makes it difficult to understand the movement of eDNA in those environments. For example, we detected O. stygocaneyi DNA in sampling unit 20.1 which is approximately 100 m upslope from sampling unit 10.1 (i.e., the only location where it has been recorded), suggesting those locations may share water during flooding. We hypothesized O. stygocaneyi may not occur in that cave, but its DNA is present due to groundwater shared among systems during particularly wet periods. Environmental DNA can persist for up to 25 d in ex-perimental ponds (Dejean et al. 2011), in terrestrial soil for at least 6 y (Andersen et al. 2012), and in cave soils for thousands of years (Hofreiter et al. 2003). In relatively stable underground aquifers, eDNA may persist for months or even years resulting in detections that are not indicative of the current population status. Alternatively, large floods can quickly move sediment and organisms out of caves (Van Gundy and White 2009; Graening et al. 2010) resulting in quick expulsion of DNA. Our model results indicated flowing water increased detection for cavefishes and cave crayfishes via eDNA surveys, which would be expected because some flow would mix and transport eDNA that had been held in the soil or deeper groundwater, but the retention time of DNA is unknown.
Substrate and water velocity also influenced detection probability of cavefishes and cavefishes via visual surveys. Visual counts of stream fishes have been used for a variety of species in clear coldwater and warmwater streams (e.g., Lambert and Hansom 1989;Heggenes et al. 1991;Brewer and Ellersieck 2011). Sampling bias via visual surveys has been associated with water velocity (Heggenes et al. 1991), water depth (Brewer and Ellersieck 2011), surface glare, turbidity, and fish behavior (Bozek and Rahel 1991). Similarly, we found that coarse substrate and flowing water were negatively associated with our ability to detect cavefishes using visual observations. Both variables were represented as binary in our model, which does not provide a measure of the magnitude of the relationship (i.e., it is a shift in the intercept, rather than a slope). Nevertheless, our findings do suggest that substrate and water velocity are factors to consider when conducting traditional visual surveys.
We found false negative samples associated with cave crayfishes were often related to the presence of pseudogenes in some species' DNA. Pseudogenes are mitochondrial genes that have moved into the nucleus, become nonfunctional, and then acquire mutations (Buhay 2009). Pseudogenes can be identified by the presence of stop codons in the sequence and "messy" chromatograms (i.e., the presence of many PCR products; Buhay 2009). Therefore, the presence of pseudogenes can make it difficult or impossible to determine the species (Buhay 2009). We found pseudogenes in the DNA of O. stygocaneyi and C. setosus, which resulted in non-specific binding of the primers and probes and lower detection probability. Future efforts might attempt use of other genetic techniques to isolate the actual mitochondrial gene (e.g., cloning, RT-PCR, long PCR, mtDNA enrichment, sequencing mitochondrial rich tissues) or target different genes; however, all of these techniques have associated difficulties to overcome (e.g., expense and technicality; Song et al. 2008;Buhay 2009).
Our data suggest increasing the number and spatial distribution of cave crayfish DNA sequences would allow researchers to design better assays that might improve detection. Knowing the genetic variation of the population is critical when designing assays to successfully amplify the DNA of the target species while avoiding amplification of any non-target taxa (Furlan et al. 2015). We had access to 23 sequences for T. rosae, 21 sequences for T. eigenmanni, and 8 for C. tartarus to represent the genetic variation across the known distribution of these species. Consequently, the assays we developed for the aforementioned species worked well. Alternatively, we only had seven C. setosus DNA sequences to represent genetic variation for a species that is more broadly distributed than other cave taxa in this study (Suppl. material 1: Table S3). More samples of genetic material across the range of more broadly distributed species would be necessary to adequately capture the species' genetic variation (Niemiller et al. 2018). We also do not have a comprehensive understanding of the genetic variation and species designations among cave crayfish populations. For example, C. setosus individuals that were collected from opposite ends of their range were genetically different by almost 6% as reflected by their CO1 gene (Suppl. material 1: Table S3, accession numbers JX514464 and MN984899). We suggest future efforts focus on understanding the genetic variation among these species.

Conclusion
Environmental DNA is a useful tool; however, the limitations we identified indicate eDNA surveys for these taxa are currently not adequate to replace traditional surveys in subterranean environments. Environmental DNA is a viable option for sampling cavefishes from locations that provide access to groundwater but cannot be physically accessed easily (i.e., springs, wells, and flooded caves). In fact, we detected cavefishes' DNA in locations where they have not been previously identified (i.e., McDonald and Ozark counties, Missouri). Further, we show that fewer surveys using eDNA would be needed for cavefishes when compared to traditional visual surveys. Environmental DNA may serve as a useful initial surveillance method when followed up by focused, on-the-ground surveys or dye tracing to identify possible sources of DNA beyond the cave. Lastly, the life history and ecological data gained from traditional surveys provide important information necessary for developing conservation strategies though increasing survey effort to adequately capture species presence should be considered if that is the sampling goal. If eDNA surveys are to be used to supplement visual sampling in subterranean environments, it would be beneficial for future efforts to 1) examine DNA movement through karst environments, 2) evaluate the genetic diversity among the Ozark Highland cave crayfishes, and 3) attempt to isolate the actual CO1 (or other) gene of cave crayfishes to improve use of eDNA in these systems. laboratory assistance were provided by T. Dropps, M. Judkins, A. Miller, S. Schneider, D. Thomson, M. Wedgeworth, J. Wiggins, and C. Wood. We appreciate the constructive feedback from K. Kuklinski who provided comments on an earlier draft. An animal care and use protocol was not required for this research because the fish were not handled by the investigators. Any use of trade, firm, or product names is for descriptive purposes only and does not imply endorsement by the U.S. Government.
Woods AJ, Omernik JM, Butler DR, Ford JG, Henley JE, Hoagland BW, Arndt DS, Moran BC (2005) Ecoregions of Oklahoma (color poster with map, descriptive text, summary tables, and photographs). U.S. Geological Survey, Reston. are reported as ambient light visible (Yes) or not (No) and fine or coarse substrate, respectively. Sampling unit 10.1 did not have ambient light on the first survey, but did on later surveys due to cave flooding. The species of cavefish (Troglichthys rosae = ros, Typhlichthys eigenmanni = eig) or cave crayfish (Cambarus aculabrum = acu, C. setosus = set, C. subterraneus = sub, C. tartarus = tar, Orconectes stygocaneyi = sty, unknown = unk) known or thought to occur at each sampling unit are also reported. Table S2. Results of the environmental DNA (eDNA) and visual surveys (Vis) for each sampling unit (SU). Yes indicates that the species was detected, no indicates that the species was not detected, and NA indicates that a survey was not completed for that sampling unit. A dash indicates water samples were collected, but we did not complete the DNA analysis.