Occupancy in dynamic systems: accounting for multiple scales and false positives using environmental DNA to inform monitoring

1 –––––––––––––––––––––––––––––––––––––––– © 2019 The Authors. Ecography published by John Wiley & Sons Ltd on behalf of Nordic Society Oikos This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited. Subject Editor: Simon Creer Editor-in-Chief: Miguel Araújo Accepted 13 October 2019 42: 1–11, 2019 doi: 10.1111/ecog.04743 doi: 10.1111/ecog.04743 42 1–11


Introduction
Landscape-scale monitoring is important to understand species occupancy through time and monitor for the impacts of habitat and climate change (Noon et al. 2012). Many survey and monitoring programs focus on repeated surveys of historically-occupied sites or sentinel sites, although habitat is often spatially and temporally dynamic (Yoccoz et al. 2001, Buckland et al. 2005. Occupancy modeling has improved estimation of species occurrence by accounting for imperfect detection (MacKenzie et al. 2002, Miller et al. 2011), but choice of sampling unit remains difficult and is further complicated in landscapes with dynamic habitat (Noon et al. 2012, Gould et al. 2019. Studies can implement dynamic occupancy models that account for changes in occupancy through time (MacKenzie et al. 2003), but these require extensive survey efforts yearly and focus on changes in occupancy of existing sites; they do not account for changes in the availability of sites as is common in highly dynamic systems. Multi-scale occupancy parameterization can allow for inferences of occupancy and habitat relationships at multiple spatial scales (Nichols et al. 2008, Mordecai et al. 2011, Pavlacky et al. 2012, Lipsey et al. 2017); however, the increased costs of sampling at a landscape level can be prohibitive, especially for rare and difficult-to-detect species.
Due to the variability in climate now and in the future, aquatic species in arid systems could benefit from efficient landscape level monitoring that considers occupancy at multiple spatial scales (George andZack 2001, Hagen et al. 2016). The need for landscape level monitoring of amphibians has been an important point of emphasis in the literature (Hamer and Mahoney 2010) and has been implemented at broad scales (Gould et al. 2012, Adams et al. 2013, Hossack et al. 2015; however, in some systems application has been challenged by low detection probabilities, which increase the replication and effort required for accurate estimation (Bailey et al. 2004, Fairman et al. 2013. Environmental DNA (eDNA) has emerged as a promising tool to survey and monitor aquatic species that are rare or elusive (Goldberg et al. 2011, Jerde et al. 2011, Thomsen et al. 2012. As eDNA is detected imperfectly, an occupancy modeling framework has been increasingly applied to interpret eDNA datasets (Schmidt et al. 2013, Hunter et al. 2015, De Souza et al. 2016, Schmelzle and Kinziger 2016 and is beginning to be applied to estimate occupancy at multiple ecological scales and at broad geographical extents (Sutter and Kinziger 2019). The use of highly-sensitive eDNA detection techniques could allow for an effective approach to model occupancy at site-level and landscape-scales.
One challenge in applying occupancy modeling to eDNA detections is that most models assume that no false positives occur, but at the PCR replicate level, false positives are widely recognized (Wilcox et al. 2013, Ficetola et al. 2015. These issues have been handled either by assuming false positives did not occur (Rees et al. 2014) or by setting ad hoc thresholds (e.g. retesting ambiguous samples and considering only samples that test positive on the second run to be positive; Goldberg et al. 2013), but these methods can bias estimates (Ficetola et al. 2015, Lahoz-Monfort et al. 2015. To account for false positives in eDNA qPCR replicates, a site confirmation approach that relies on an additional detection method (that does not have false positives) can provide confirmation of site occupancy and be used to estimate the rate of false positives (Miller et al. 2011, Chambert et al. 2015. False positive detection models have been further extended to incorporate a probability of false positive detections at two hierarchical levels (Guillera-Arroita et al. 2017). False positives have been estimated with eDNA data, but to date they have not been used in modeling habitat occupancy relationships or incorporated to assess influence on model selection and ecological inference (Guillera-Arroita et al. 2017).
Riparian areas make up 1-3% of arid landscapes; aquatic species in these landscapes rely on these isolated areas of habitat (Patten 1998). These limited patches are further influenced by seasonal and between year variation that produce dynamic conditions for survival and connectivity for aquatic species (Naiman and Decamps 1997). Consequently, changes in aquatic species occupancy in arid systems can be difficult to detect and interpret at the scale of individual site, presenting challenges for species monitoring into the future. Hydrologic variability in arid environments makes watershed occupancy a useful scale to monitor populations that rely on permanent aquatic environments. A loss of a single site could be negligible if connectivity to other occupied sites within the watershed remain, however, if watershed occupancy decreases substantially due to warmer and drier climate, connectivity and site-level occupancy could be severely reduced. However, factors influencing occupancy at sites can also have a large impact on overall population viability. Therefore, monitoring at multiple scales may be the most effective approach for understanding the status of populations, especially in highly dynamic systems.
The Great Basin clade of the Columbia spotted frog Rana luteiventris lives in isolated populations in semi-arid portions of southeastern Oregon, southwestern Idaho and Nevada, USA. Columbia spotted frogs are highly aquatic at all life stages (Reaser and Pilliod 2005). Under future climate scenarios, the Great Basin is predicted to become hotter and have increased variability in precipitation (Hurd et al. 1999, Cubashi et al. 2001) and the Columbia spotted frog is predicted to have a reduced area of suitable climate . Previous surveys in the northern Great Basin did not detect Columbia spotted frogs at a large portion of historical sites, but a relatively high proportion of proximal sites were occupied, leading to difficulties in the interpretation of changes in occupancy and population declines (Wente et al. 2005).
We propose to alleviate common barriers to the implementation of long-term monitoring programs in dynamic systems by demonstrating an approach to modeling occupancy at multiple spatial scales simultaneously and employing a rapid and simple survey method. We conducted paired eDNA and traditional field surveys (visual observation and dipnet survey) across the Great Basin range of the Columbia spotted frog and applied a multi-scale occupancy modeling framework accounting for false positives to investigate factors influencing the occupancy and detection of this amphibian in a highly dynamic system. In addition, we further investigated how eDNA methods can improve detection probabilities compared to traditional field methods and, evaluated the occurrence of false positives in eDNA sampling, influence on model selection and interpretation of our multi-scale occupancy model. These results can be used to inform future monitoring design for difficult to detect species in dynamic aquatic systems.

Site selection
We selected 220 survey sites throughout the range of the Great Basin clade of the Columbia spotted frog in southeastern Oregon, southwestern Idaho and Nevada. To gather information on potential locations of this rare amphibian, we used the historical dataset and predictive maps of Pilliod et al. (2015), the National Wetland Inventory (NWI; US Fish and Wildlife Service), and two meetings with biologists from the Bureau and Land Management, United States Fish and Wildlife Service, United States Forest Service and Nevada Dept of Wildlife. For site selection, we first randomly selected 21 known locations from historical records compiled by Pilliod et al. (2015), stratified by geographic area (Fig. 1). To identify new sampling sites near these known locations, we then selected two 12th level hydrologic units (Watershed Boundary Dataset; US Geological Survey) adjacent or one hydrologic unit away from the units containing historical sites because all adjacent hydrologic units were already being sampled (i.e. hydrologic unit was already selected for sampling because another historical site was randomly selected for sampling in the same or adjacent watershed). When multiple watersheds were available for selection, we selected watersheds in the following order of priority: those with sites suggested by resource managers, those with the highest proportion of predicted suitable climate , and those with the greatest wetland area (from NWI) (n = 42 watersheds). Three individual sites within each watershed were then selected (in order of priority) from: on the ground knowledge of collaborators, random selection from semipermanent wetlands (NWI), and those identified through Google Earth (Google) when there were no NWI wetlands, but because of water availability and limitations to public access, the number of sites varied (n = 126; Fig. 2). When the pre-selected sites where unavailable for sampling during site visit these were replaced when suitable habitat was observed in the area. The choice of watershed for neighboring population sampling was not influenced by the presence of additional historical locations. We implemented this categorized sampling approach to provide a site level understanding of Columbia spotted frog populations for resource managers and account for the dynamism of the landscape that is a challenge for landscape-scale monitoring of aquatic systems in arid landscapes.
For additional survey sites, we randomly selected semipermanent wetlands that had a probability of climatic suitability above 0.20 (n = 49) and below 0.20 (n = 24) from a climatic suitability model developed for Columbia spotted frogs in the Great Basin ; for flowchart of site selection see Fig. 2). The additional sites were added to potentially identify new populations of Columbia spotted frogs in a highly fragmented landscape. We selected 0.20 climatic suitability as a cut-off to include sites unlikely to have Columbia spotted frog to better characterize factors influencing occupancy. The additional sites provided the opportunity for an increased chance of identifying new populations while accounting for the variability of aquatic resources in the study area. The approach outlined was not a completely randomized design; however, taking a more targeted approach was necessary due to highly disjunct known populations and to account for changes in a dynamic landscape with infrequent aquatic resources (i.e. if we had used random sampling, we would have had very few detections). We did not explicitly test for spatial autocorrelation between sites, focused on including sites across the range of climatic conditions in the Great Basin, including sites unlikely to have Columbia spotted frogs. Through this design, our goal was to collect a representative sample with enough power to gain understanding of the covariates important to Columbia spotted frog occupancy at multiple scales while balancing the risk of overestimating precision through correlated samples. All sites were selected from public land. Sites were surveyed in May-August of 2015 and 2016 to target tadpole season, when breeding sites for the species can be most clearly identified. Surveys were conducted starting at lower elevations earlier in the summer to track the gradient in breeding season timing. We surveyed geographic areas A, B, C and G in 2015 and B, D, E and F in 2016 ( Fig. 1).

Environmental DNA sample collection and visual/dipnet surveys
We performed visual/dipnet surveys and collected eDNA samples at each site, defined as the specific wetland, pond or pooled water that was selected for sampling based on our sampling approach. Upon initial approach to a site, we conducted a visual survey where two biologists surveyed different portions of the site until the entire site was observed. After the visual survey, we collected eDNA samples before any other equipment or clothing touched the water to prevent cross-contamination between sites. Following the collection of water and eDNA filtering, we performed dipnetting around the observed site moving along the shoreline approximately every 1 m when accessible.
We identified areas where frogs or tadpoles would most likely occur to collect water samples for filtering, focusing on areas without canopy cover, with emergent vegetation, and 4 along south facing shorelines. We collected three temporally replicated filter samples by filtering 250 ml of water for each filter at a site; per sample detection rate in a pilot study using these methods with Columbia spotted frogs was found to be 0.66, indicating that > 0.95 detection can be achieved with three samples (Goldberg and Pilliod unpubl.). Most sites had 1 sampling location (i.e. physical location within the site where the filter sample was taken); however, we collected samples from 2 to 6 additional locations within the site and combined equal portions from each location into a single sample for 14 sites, when sites were approximated to be > 1200 m 2 .
We used single-use whirl-paks (Nasco) to collect the individual water samples. We attempted to filter 250 ml for each sample; when the filter clogged before reaching this amount, we recorded the amount successfully filtered. We filtered eDNA samples through a 0.45 µm cellulose nitrate filter within single use sterile analytical test filter funnel (Nalgene). After filtration was complete, we folded the filter in half inwards and placed into an individual coin envelope marked with unique identification, location and date. The coin envelope was then either placed into a larger bag containing silica desiccant until completely dried (2015) or placed into an individual bag with silica desiccant (2016); this approach was adapted from Carim et al. (2014). There was no evidence of contamination from either method, but to further safeguard against possible contamination observed in other studies we took the more careful approach in 2016. Any contact with the filter was either with a new latex glove or with forceps sterilized in 50% bleach and rinsed with distilled water. We also collected one negative control at each site by filtering distilled water on site using the same materials and methods as the field samples. We decontaminated boots with 10% commercial bleach solution between each site and filtration materials downstream of filter once per week or sooner if moving to new geographical area.

eDNA lab protocol and analysis
We extracted and prepared all eDNA samples for quantitative PCR (qPCR) in a lab dedicated to low quality DNA extraction where high quality DNA and PCR product are restricted. Access is limited and researchers must shower and change clothes before entering after being exposed to high quality DNA or PCR product. We extracted DNA from samples using the Qiagen DNeasy Blood and Tissue Kit with the addition of the Qiagen QIAshredder as described in Goldberg et al. (2011). Each extraction batch included a negative control to test for contamination.
Great Basin Columbia spotted frog presence was assessed using targeted sequence detection through quantitative PCR (qPCR). The assay was designed from sequences from Funk et al. (2008) obtained through GenBank. Sequences for the Great Basin clade were aligned using Clustal X2 (Larkin et al. 2007) and assay designed using Primer Express 3.0 (Life Technologies). Primers were analyzed for speciesspecificity through Primer-BLAST (NCBI) and a set was chosen that would not cross-amplify with other species found throughout the range of the Columbia spotted frog. Of western North American species, this database indicated that the primers may create products for Rana cascadae, R. pretiosa, R. aurora and R. sierrae under some conditions, however these species are not found within the range of this clade. We validated this assay (Supplementary material Appendix 1 Table  A1) by testing it against 10 Great Basin Columbia spotted frog tissues and five of each of the following: northern leopard frog Lithobates pipiens, American bullfrog L. catesbeianus, Woodhouse's toad Anaxyrus woodhousii, Pacific treefrog Pseudacris sierrae, western toad Bufo boreas and Great Basin spadefoot Spea intermontana. We analyzed each sample in triplicate. Each qPCR included 3 µl of DNA extract in a total volume of 15 µl. Reactions were run using 1× QuantiTect Multiplex PCR Mix (Qiagen), 0.2 µM of each primer and 0.2 µM of probe on a Bio-rad CFX96 Touch Real-Time PCR Detection System. To test for inhibition of the qPCR, each well included an exogenous internal positive control (IPC; Applied Biosystems). Reactions activated for 15 min at 95°C then ran for 50 cycles of 94°C for 60 s followed by 60°C for 60 s. All qPCR plates included negative controls to test for contamination. Positive standards consisted of a serial dilution of 10 −3 through 10 −6 of extracted DNA from a toe clip of a Great Basin Columbia spotted frog quantified on a Qubit 3.0 Fluorometer (Life Technologies; 67.4 ng µl −1 ) and were included in duplicate on each qPCR plate. All of the samples from the target species tested positive and all of the samples from non-target species tested negative.

False positive estimation
Even though the Columbia spotted frog assay showed high species specificity, the sensitivity required to analyze rare and/or degraded templates such as eDNA (i.e. extended PCR cycle counts) can lead to false positives (Wilcox et al. 2013). replicates (1 = exponential amplification of target DNA or 0 = no amplification) and detection (1) or non-detection (0) using visual/dipnet survey. We modeled a constant parameterization (ψ(.) fp(.) p(.)) of the false positive occupancy modeled in the Unmarked R package (Fiske and Chandler 2011). Identifiability issues can occur in a site confirmation design (Chambert et al. 2015, Guillera-Arroita et al. 2017. To address this concern, we used the two-stage false positive model developed by Guillera-Arroita et al. (2017) that produces completely identifiable parameter estimates with the addition of two calibration experiments. We ran 50 qPCR blanks that consisted of the Columbia spotted frog eDNA assay with 3 µl of deionized water and included our field negatives to account for false positives at the sample level. We examined profile likelihoods for multiple local maxima and assessed parameter estimates for agreement with the Miller et al. (2011) false positive occupancy model. We subsequently applied the estimated probability of a false positive at the qPCR replicate level from the Miller et al. (2011) model to our qPCR detections that informs the filter sample detection in the following multi-scale occupancy analysis.

Multi-scale occupancy
The multi-scale occupancy detection histories were made up of three temporal replicate eDNA filter samples that were determined to be a detection (1) or non-detection (0) based on the amplification of Columbia spotted frog mtDNA through qPCR. If any of the three qPCR replicates per filter sample was positive, then the filter sample was considered a detection. We combined survey seasons (2015 and 2016) into a single season. Occupancy models require independent repeat sampling to account for detection, however, there is no requirement that sites are physically visited on separate occasions (Mackenzie and Royle 2005). In fact, depending on the system or species, time between survey events can potentially lead to a failure in assuming closure of sites to changes in true occupancy and in general to provide a snapshot of the system surveys should be conducted as quickly as possible (Mackenzie and Royle 2005). Environmental DNA sampling at the site (i.e. wetland, pond or pooled water) level as outlined here represents three independent filter samples from a single sampling location (i.e. not spatially replicated) taken in short succession (seconds) that provide temporal replication. We estimated watershed occupancy using spatial replication of watersheds with two or more sites (n = 62). Sites were either isolated within the watershed due to intermittent flow or at a distance and flow rate unlikely to violate the assumption of independence, however, if independence was violated it would likely lead to an overestimation our detection probability that would then underestimate watershed occupancy. We used the Nichols et al. (2008) multi-scale model and similar to Schmidt et al. (2013), but redefined our scales of interest as: where, z i , describes the presence or absence of the species at watershed i, as a function of the occupancy probably ψ i . The next hierarchal level, a ij , describes the presence or absence of the species at site (i.e. wetland, pond or pooled water) j in watershed i, as a function of the occurrence state at watershed i and the occupancy probability of a site, θ ij . The observed data, y ijk , describes the presence of eDNA in filter sample k from site j in watershed i, as a function a ij and the detection probability given eDNA is in the filter sample and is amplified by qPCR, p ijk . Therefore, the model parameters consisted of a broad-scale occupancy probability (ψ), which represented the probability of occupancy within a watershed, a fine-scale occupancy probability (θ) of a site (i.e. wetland, pond or pooled water) being occupied conditional on the watershed being occupied, and the probability of detection (p) of eDNA given it is in the filter sample and amplified by qPCR that is conditional on occupancy of the site and watershed. An example detection history of 101 000 would indicate that Columbia spotted frogs were detected in filter sample one and filter sample three at site one and were not detected at site two with both sites located within the same watershed. We note that this parameterization differs from previous hierarchal modeling with eDNA (Schmidt et al. 2013) by combining availability within the filter sample and qPCR detection into a single detection parameter (p). We did this in order to model occupancy at multiple ecological scales within the contemporary software limitation of three hierarchical levels. We made this compromise to identify important ecological differences in occupancy across multiple scales and account for the dynamism of landscapes when being monitored for a species of interest. Future developments in modeling multi-scale occupancy with eDNA would benefit from the addition of a 4th hierarchical level that would disentangle availability within the filter sample and qPCR detection while still estimating occupancy at multiple ecological scales.
We hypothesized that watershed occupancy (ψ) was influenced by up to four different watershed covariates: average precipitation, average temperature, total length of perennial streams and average slope of the stream network (Supplementary material Appendix 1 Table A2a). We hypothesized that average precipitation and total length of perennial streams would have a positive relationship and average temperature and average slope would have a negative relationship on watershed occupancy. We considered average slope as a potential influence on occupancy because Columbia spotted frogs occurred mostly in stream habitats in the study area, but in slow or still side channels or margins of streams (Reaser and Pilliod 2005). We hypothesized that more gradual slopes within the stream network would be associated with these aquatic site characteristics. Average precipitation included 7 rainfall and snow melt. We used a combination of watershed covariates to test five potential hypotheses on the effect on watershed occupancy (Supplementary material Appendix 1  Table A3). We hypothesized site occupancy (θ) was influenced by up to six covariates: emergent vegetation, water type, livestock impact, average precipitation, average temperature and elevation (Supplementary material Appendix 1 Table A2b). We used a combination of site covariates to test five potential hypotheses of the drivers of site level occupancy (Supplementary material Appendix 1 Table A3). Finally, we hypothesized that eDNA detection at a site was influenced by up to two covariates: temperature (°C) in a negative relationship and pH in a positive relationship (Supplementary material Appendix 1 Table A2c). At each site, we measured these parameters just under the surface of the water after eDNA samples were collected with a multiparameter meter (OAKTON Instruments, Vernon Hills, IL). We modeled detection probability (p) for eDNA filter samples as a constant, covariate of pH or covariate of water temperature. We calculated multicollinearity of covariates within each level of the multi-scale occupancy model (ψ, θ and p) by calculating variance inflation factors.
We ran all combinations of hypotheses at each level for a total of 75 models with each multi-scale occupancy level having five or fewer potential models (Supplementary material Appendix 1 Table A3). Multi-scale occupancy models were fitted using the Nichols et al. (2008) multiscale occupancy parametrization in Program Mark 8.1 (White and Burnham 1999). Models were ranked according to Akaike's information criterion corrected for small sample size (AICc) and assessed for support based on AICc weights and delta AICc.
We estimated the potential impact of false positives on inference in our dataset by applying the estimated probability of a false positive to each positive qPCR replicate in our multi-scale occupancy dataset to create 30 full-dataset replicates. In this approach positive detections were classified as a false detection and changed from a 1 to a 0 in the detection histories used for the multi-scale dataset if a random number between 0 and 1 was lower than the probability of a false positive for that technical replicate. The reclassified qPCR replicates for each sample were then converted to a sample level detection (1) if any of the three replicates were positive or no detection (0) if none. We applied the adjusted encounter histories to our multi-scale occupancy model set to assess differences in model selection and parameter inference.

Results
We analyzed a total of 660 eDNA samples from 220 sites with a mean filtered volume of 248 ml (SD = 15 ml). Two sites had evidence of on-site contamination of target DNA (field negatives testing positive) and were excluded from occupancy analyses. The contaminated field negatives were from 2016 and we believe that the contamination was due to an error in handling in the field and was isolated to these two sites. Target DNA was not detected in any extraction negatives or qPCR negatives. We detected Columbia spotted frogs at a greater proportion of sites using eDNA than traditional surveys, with a naïve occupancy estimate of 0.37 with eDNA and 0.20 with visual/dipnet surveys (Table 1). We detected Columbia spotted frog at 86% of historical sites and 36% of neighboring sites. Sites with < 0.20 climatic suitability (n = 24) had a single Columbia spotted frog detection while 41% of sites > 0.20 climatic suitability had detections. Of the sites where eDNA detections occurred, qPCR replicates were highly consistent, with only 17% of positive sites having only 1 positive qPCR replicate and 76% having all three positive replicates.

False positive estimation
False positive detection was modeled using each individual filter sample (n = 654) as the sampling unit with an encounter history. The full dataset had 504 qPCRs that amplified for our target sequence. None of the 50 qPCR blanks amplified for Columbia spotted frog.

Multi-scale occupancy
Out of the full dataset, sixty-one watersheds had > 1 (range 2-8) sites per watershed for a total of 168 sites available for use in the multi-scale analysis. The top two candidate models were nested and together had 0.82 of the evidence weight ( Table 2). The most supported model from the candidate set included climatic covariates of precipitation and temperature at the watershed occupancy scale, emergent vegetation and water type as covariates of site occupancy, and water temperature as a covariate of detection ( Table 2). The next model included all covariates of the top model with the addition of stream length as a covariate of watershed occupancy but was considered uninformative due to its 95% confidence interval overlapping zero and the extra parameter did not improve AIC over the top model (Table 2). Water type was present in the top two candidate models, but was also considered uninformative due to 95% confidence intervals including zero (Arnold 2010). Variance inflation factors did not exceed 4.39.
Precipitation had a positive association with watershed occupancy (Fig. 3A) while temperature had a negative impact on watershed occupancy (Fig. 3B). Aquatic covariates of emergent vegetation had a positive influence on site occupancy (Supplementary material Appendix 1 Table A4). Water temperature had a positive effect on eDNA detection (Fig. 4), opposite of our hypothesis. At the mean water temperature 8 (18.4°C), eDNA sampling had a high probability of detection per sample, p = 0.78 (SE = 0.03; 95% CI: 0.71-0.84).
The resampled datasets that accounted for false positives had an average of 379.4 positive qPCRs (n = 30) out of the original 389 positive qPCRs. Twenty-seven resampled datasets that accounted for false positives did not lead to any change in model selection with a combined model weight of the top two models ranging from 0.66 to 0.88 (mean = 0.80). Two resampled datasets showed slight changes in model selection compared to when false positives were not accounted for with additional support for a third model that no longer included water temperature as a covariate of detection. An additional resampled dataset also had the top two models switched in most supported, with the top model including streams as a covariate of ψ with an AICc weight of 0.37. The parameter estimates of the resampled datasets for the top model remained similar to the original dataset and accounting for false positives did not change parameter interpretation (Supplementary material Appendix 1 Table A4).

Discussion
Landscape-scale monitoring of species occupancy in variable systems can be highly challenging, but improvements in accuracy across scales for aquatic species can be made through new methodologies. We found that multi-scale modeling of occupancy of an aquatic species in a dynamic arid system informed by eDNA sampling could lead to improved monitoring through increased detections, including when the probability of false positives was incorporated. Using eDNA increased our total detections in this system and proved to be more sensitive than our visual/dipnet surveys, allowing for an effective sampling strategy across the range of the Great Basin clade of Columbia spotted frogs. Similarly, greater detection sensitivity of eDNA compared with field surveys has been observed for other species and systems (Dejean et al. 2012, Pilliod et al. 2013, Schmelzle and Kinziger 2016.
Contrary to expectations, we found that eDNA detection probabilities increased with increasing temperatures, indicating that degradation was not the primary influence on detection in this study. Increased eDNA detection of amphibians in warmer water could be caused by several factors, including increased activity in warmer water, increased metabolism, faster tadpole growth and greater frequency of time spent in the water during hot conditions. Warmer water temperatures can increase amphibian activity (Brattstrom 1963), leading to the possibility of greater eDNA shedding and higher detection. However, previous studies observed no effect of water temperature on eDNA production for fish in microcosm experiments (Takahara et al. 2012, Klymus et al. 2015. It is likely that the increase in detection probability at warmer temperatures is being influenced by behavioral choices or a physiological response of Columbia spotted frogs in this system. Accounting for this variation in detection improved estimates of watershed and site occupancy, indicating that water temperature is an important measure to collect during landscape-scale monitoring efforts based on eDNA data. One challenge with using eDNA for occupancy estimation is that the probability of a false positive at the technical replicate scale is non-negligible (Ficetola et al. 2015, Lahoz-Monfort et al. 2015. Quantitative PCR replication can help improve parameter estimates when a degree of uncertainty is assigned to the number of positive replications and to a greater degree using an additional detection method (Ficetola et al. 2015, Lahoz-Monfort et al. 2015. This study used empirical data to estimate the rate of false positives for eDNA qPCRs using multiple detection methods and found that it had a minimal influence on model selection results in this dataset, which was likely influenced by the consistency among replicates of the eDNA data and high detection probability. When available, an estimation of the probability of false positives from eDNA samples is relatively straightforward and can improve or support the inferences being made of occupancy and detection (Lahoz-Monfort et al. 2015). It is more likely, but not always the case, that a site with fewer positive detections (i.e. a single sample with a single positive qPCR) represents a false detection and when removed will influence model outcomes. It is important to acknowledge that this method does not identify true false positive sites, but rather applies a probability that any single detection may be a false positive. It is still necessary to consider variability within specific systems because the impact of false positives is likely to be greater when detection probability is low, and the target species is rare. We found little impact of the removal of potential false positives on modeling results because positive sites tended to have multiple positive wells across multiple samples. Future monitoring using this sampling and analysis design could use this false positive rate (0.023 per technical replicate) to categorize sites with one out of nine technical replicates positive as uncertain, as these would have a 0.17 probability of being a false positive (dropping below 0.02 if two technical replicates tested positive). If the status of a site characterized as uncertain is needed for management or conservation, additional sampling would be necessary. Arid and semi-arid landscapes like the Great Basin of North America have a high degree of variation in water availability and an increase in variability of precipitation as well as predicted hotter temperatures under future climate scenarios (Hurd et al. 1999, Cubashi et al. 2001, likely decreasing available habitat for aquatic species. The results from our eDNA surveys were consistent with previous studies that described emergent vegetation and climatic variables as important influences on Columbia spotted frog presence in the Great Basin (Arkle andPilliod 2015, Pilliod et al. 2015). The multi-scale analysis in this study highlights the importance of climatic factors influencing Columbia spotted frog occupancy at a watershed scale while site presence was influenced by finer scale aquatic characteristics. By using this multi-scale approach, habitat relationships at multiple scales can be incorporated into landscape level monitoring while accounting for future variability across watersheds and sites. The specific modeling results here can help resource managers identify 12th level hydrologic units that have a higher probability of occupancy to identify newly occupied watersheds or target watersheds for restoration or conservation. This research provides the necessary framework to design and implement a long-term monitoring program for the Great Basin clade of the Columbia spotted frog that takes into account the changing aquatic resources in this system.
A primary goal of long-term monitoring is to identify declines and use that information to allocate scarce resources. Assessments of declines from historically occupied sites through long-term monitoring can be complicated when spatial distribution of occupied sites shifts due to changing resource availability and/or population dynamics (Hecnar andM'Closkey 1996, Cruickshank et al. 2016). This leaves the potential that identified declines in occupancy are just variation in the spatial locations that are occupied   (Yoccoz et al. 2001, Skelly et al. 2003. The application of a multi-scale occupancy model could help account for this variation and the spatial distribution when monitoring occupancy of patchily distributed species. For long-term monitoring of occupancy, the approach to choosing sampling sites within watersheds as well as watersheds to be sampled would have to be standardized to allow for comparison of occupancy estimates across years. Alternatively, all sites within the watershed or study area could be surveyed, but this is not feasible for landscape-scale or broader species assessments and monitoring. This study demonstrates the utility of eDNA sampling across a large geographical area to assess ecological questions of species occurrence at multiple scales while accounting for false positives. This approach brings together multiple advances in eDNA and occupancy modeling and lays the groundwork for improved applications of eDNA detection to address ecological and conservation questions. We implemented a robust eDNA sampling scheme that had greater detection sensitivity over traditional detection methods. The multi-scale approach we applied allowed us to examine factors at the individual habitat patch but also the watershed scale that is more relevant to monitoring the status of this species in this dynamic landscape, while accounting for the false positives that have the potential to derail eDNA-based decision-making. We demonstrated that eDNA can be a highly sensitive and reliable method to assess amphibian occupancy across a broad geographic area. The approach we took for estimating false positives can be undertaken whenever an unambiguous survey method (e.g. dipnetting or trapping) can be conducted alongside with replicate eDNA samples. The resulting probability can then be used to create a threshold for the actionability of an eDNA result. The greater sensitivity and efficiency afforded by eDNA detection can improve landscape scale monitoring of aquatic species and create greater confidence in methods needed to continue to monitoring species trends into the future.