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).

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.

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 (Graening et al. 2006a). Cambarus setosus is the widest-ranging cave crayfish of the Ozark Highlands and has been documented 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 south-central 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).

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.

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).

Figure 2. 

Filtration setup for eDNA collection. 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).

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. 2006a, 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).

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 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³), 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.

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’-CATACCCCAAACATGTTGGT-3’ and reverse primer ND2B: 5’-TGGTTTAATCCGCCTCAGCC-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₂ (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₂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 fragment 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. 2012, 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.

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 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).

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₂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 remaining 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₂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.

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 right-skewed 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:

$\begin{array}{l}logit\left(p_{ij}\right)={\beta }_0+{\beta }_1{X}_{1[i,j]}+{\beta }_2{X}_{2[i,j]}+{\beta }_3{X}_{3[i,j]}+{\beta }_4{X}_{4[i,j]}+{\beta }_5{X}_{5[i,j]}+{\beta }_6{X}_{6[i,j]}+{\beta }_7{X}_{7[i,j]}+\\ {\beta }_8{X}_{1[i,j]}{X}_{2[i,j]}+{\beta }_9{X}_{1[i,j]}{X}_{3[i,j]}+{\beta }_{10}{X}_{1[i,j]}{X}_{4[i,j]}+{\beta }_{11}{X}_{1[i,j]}{X}_{5[i,j]}+{\beta }_{12}{X}_{1[i,j]}{X}_{6[i,j]}+\\ {\beta }_{13}{X}_{1[i,j]}{X}_{7[i,j]}+{\beta }_{14}{X}_{2[i,j]}{X}_{3[i,j]}+{\beta }_{15}{X}_{2[i,j]}{X}_{4[i,j]}+{\beta }_{16}{X}_{2[i,j]}{X}_{5[i,j]}+{\beta }_{17}{X}_{2[i,j]}{X}_{6[i,j]}+\\ {\beta }_{18}{X}_{2[i,j]}{X}_{7[i,j]}+{\beta }_{19}{X}_{1[i,j]}{X}_{2[i,j]}{X}_{3[i,j]}+{\beta }_{20}{X}_{1[i,j]}{X}_{2[i,j]}{X}_{4[i,j]}+{\beta }_{21}{X}_{1[i,j]}{X}_{2[i,j]}{X}_{5[i,j]}+\\ {\beta }_{22}{X}_{1[i,j]}{X}_{2[i,j]}{X}_{6[i,j]}+{\beta }_{23}{X}_{1[i,j]}{X}_{2[i,j]}{X}_{7[i,j]}\end{array}$

for i = 1,2,…N, for j = 1,2,…J,

where p_ij is detection probability for survey j at sampling unit i, β₀ is the intercept, β₁ is the taxa main effect coefficient, β₂ is the method main effect coefficient, β₃ is the light main effect coefficient, β₄ is the turbidity main effect coefficient, β₅ is the velocity main effect coefficient, β₆ is the substrate main effect, β₇ is the volume main effect coefficient, β₈ is the taxa * method interaction term coefficient, β₉ is the taxa * light interaction term coefficient, β₁₀ is the taxa * turbidity interaction term coefficient, β₁₁ is the taxa * velocity interaction term coefficient, β₁₂ is the taxa * substrate interaction term coefficient, β₁₃ is the taxa * volume interaction term coefficient, β₁₄ is the method * light interaction term coefficient, β₁₅ is the method * turbidity interaction term coefficient, β₁₆ is the method * velocity interaction term coefficient, β₁₇ is the method * substrate interaction term coefficient, β₁₈ is the method * volume interaction term coefficient, β₁₉ is the taxa * method * light interaction term coefficient, β₂₀ is the taxa * method * turbidity interaction term coefficient, β₂₁ is the taxa * method * velocity interaction term coefficient, β₂₂ is the taxa * method * substrate interaction term coefficient, β₂₃ is the taxa * method * volume interaction term coefficient, X₁ is taxa, X₂ is method, X₃ is light, X₄ is turbidity, X₅ is velocity, X₆ is substrate, and X₇ 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 model-selection 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 p-value 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.

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.

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³–800.00 m³; mean ± SD = 61.21 ± 132.00 m³). 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.

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 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.

Figure 3. 

Cumulative detection probability as a function of the number of surveys for cave crayfishes (panel A) and cavefishes (panel B). Solid lines eDNA surveys and dashed lines are visual surveys. Cumulative detection probability (p_c) was calculated as: p_c = (1 - (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.