Evaluating artificial shelter arrays as a minimally invasive monitoring tool for the hellbender Cryptobranchus alleganiensis

Hellbenders Cryptobranchus alleganiensis are critically imperiled amphibians throughout the eastern USA. Rock-lifting is widely used to monitor hellbenders but can severely disturb habitat. We asked whether artificial shelter occupancy (the proportion of occupied shelters in an array) would function as a proxy for hellbender abundance and thereby serve as a viable alternative to rock-lifting. We hypothesized that shelter occupancy would vary spatially in response to hellbender density, natural shelter density, or both, and would vary temporally with hellbender seasonal activity patterns and time since shelter deployment. We established shelter arrays (n = 30 shelters each) in 6 stream reaches and monitored them monthly for up to 2 yr. We used Bayesian mixed logistic regression and model ranking criteria to assess support for hypotheses concerning drivers of shelter occupancy. In all reaches, shelter occupancy was highest from June–August each year and was higher in Year 2 relative to Year 1. Our best-supported model indicated that the extent of boulder and bedrock (hereafter, natural shelter) in a reach mediated the relationship between hellbender abundance and shelter occupancy. More explicitly, shelter occupancy was positively correlated with abundance when natural shelter covered <20% of a reach, but uncorrelated with abundance when natural shelter was more abundant. While shelter occupancy should not be used to infer variation in hellbender relative abundance when substrate composition varies among reaches, we showed that artificial shelters can function as valuable monitoring tools when reaches meet certain criteria, though regular shelter maintenance is critical.


INTRODUCTION
A major challenge associated with monitoring rare herpetofauna is balancing the tradeoff between survey rigor and disturbance to organisms and their habitat. Reptiles and amphibians often have cryptic habits and many imperiled species occur in low density. Because of these traits, multiple, rigorous surveys are often required to obtain accurate and precise estimates of population status (MacKenzie et al. 2002, Dorazio & Royle 2005, Mazerolle et al. 2007. A valid concern for many biologists is the increased disturbance to organisms and their habitat that is often inherent in increased survey effort. Capture and handling can have short-term effects on physio logy of vertebrates (Romero et al. 2008), and re peated cap-ture can yield long-term sub-lethal effects (e.g. myopathy) in some species (Cattet et al. 2003). Aside from direct effects of capture, physical search efforts can lead to disturbance of critical microhabitat that may adversely affect populations (Pike et al. 2010). Thus, there is a need for minimally invasive survey and monitoring tools for many imperiled herpeto fauna.
The hellbender Cryptobranchus alleganiensis is an amphibian of growing conservation concern that is widely suspected to be negatively impacted by traditional survey methods. Hellbenders are large (up to 74 cm total length), fully aquatic, stream-dwelling salamanders that are endemic to the eastern United States (Nickerson & Mays 1973). Hellbenders have undergone precipitous declines since the 1970s (Wheeler et al. 2003, Burgmeier et al. 2011, Graham et al. 2011. The eastern subspecies C. a. alleganiensis is listed as endangered and/or a species of special conservation concern in most states where it occurs (USFWS 2011a), and a distinct population segment in Missouri was recently recommended for federal listing under the US Endangered Species Act (ESA; USFWS 2019). The Ozark subspecies C. a. bishopi is federally listed as an endangered species under the ESA (USFWS 2011b). Obtaining baseline estimates of, and monitoring changes in, occurrence and abundance of hellbenders is a primary goal of biologists throughout the species' range (Quinn et al. 2013, Pugh et al. 2016, Pitt et al. 2017. Hellbenders rely heavily on rock crevices for shelter and nesting and are notoriously difficult to detect. Snorkeling while rock-lifting is recognized as the most effective method to detect all age classes and is, by far, the most common sampling method (Nickerson & Krysko 2003). However, disturbing large cover rocks can destroy critical microhabitats that hellbenders and their prey depend on and place hellbenders and surveyors at risk of injury (Browne et al. 2011). Furthermore, because detectability while rock-lifting can vary with site characteristics, survey conditions, and hellbender size class, multiple surveys per site are necessary to properly account for detectability . Increasingly, biologists have become concerned about the potential damage that rock lifting surveys inflict on hellbenders and their habitat and must decide between increasing survey intensity to maximize data quality or minimizing risks to hellbenders.
Artificial cover objects have been used to minimize survey disturbance for a wide range of reptiles and amphibians (Willson & Gibbons 2010) and might offer a minimally invasive alternative to rock-lifting for hellbenders. Examples of artificial cover for her-petofauna include metal and wooden cover boards for fossorial salamanders, frogs, snakes, and lizards (Grant et al. 1992), arboreal cover boards for lizards (Nord berg & Schwarzkopf 2015), PVC pipes for treefrogs (Boughton et al. 2000), leaf-litter bags for stream-dwelling salamanders (Waldron et al. 2003), and subterranean hibernacula for newts (Dervo et al. 2018). Hellbender artificial shelters, also referred to as 'nest boxes', were recently designed as a conservation tool (Briggler & Ackerson 2012). Artificial shelters were originally designed to supplement nesting cavities and have been primarily deployed in captive-breeding facilities and free-flowing streams to attract spawning hellbenders for the purpose of collecting eggs for captive propagation (Ettling et al. 2013). However, hellbenders use artificial shelters outside of the breeding season (Briggler & Ackerson 2012), suggesting that artificial shelter arrays may be a more broadly applicable monitoring tool. Specifically, searching established artificial shelter arrays to detect hellbenders may function as a less invasive surrogate for rock lifting while snorkeling, similar to using cover board surveys as a surrogate for natural cover searches that re quire lifting rocks, fallen logs, and tree bark to detect terrestrial herpetofauna (Grant et al. 1992).
The overarching goal of our study was to investigate the efficacy of artificial shelter arrays as a monitoring tool for hellbenders. Our specific objectives were to investigate how temporal and site-specific factors might influence hellbender use of artificial shelter and determine whether the proportion of occupied shelters on a given occasion (hereafter, shelter occupancy) might function as a reliable index of local hellbender abundance. We hypothesized that shelter occupancy in all arrays would vary temporally due to seasonal shifts in hellbender activity patterns and the length of time allowed for hellbenders to discover recently deployed shelters. We hypo thesized that shelter occupancy would vary among arrays due to site-specific hellbender population density, natural shelter density, or both, and we explicitly evaluated relative support for 3 non-mutually exclusive hypotheses to explain among-site variation in shelter occupancy. We refer to these as the (1) 'population density driven occupancy', (2) 'natural shelter driven occupancy', and (3) 'natural shelter per capita driven occupancy' hypotheses. A key assumption for all 3 hypo theses is that hellbenders recognized artificial shelter, when encountered, as a viable alternative to natural shelter.
Under 'population density driven occupancy' we would predict that artificial shelter occupancy in -creases as population density of the target species increases (Fig. 1A). The log-linear (i.e. threshold) trend in our prediction is used to represent the idea that shelter occupancy will reach a saturation point (i.e. 100% occupancy) at and above some threshold of population density. A key assumption specific to this hypothesis is that artificial shelter occupancy is independent of natural shelter availability. It is worth noting that 'population density driven occupancy' is an implicit assumption in many herpetofauna studies, exemplified when raw counts obtained from artificial cover arrays are used to infer relative abundance among sites or between treatments. However, to our knowledge, this hypothesis has not been explicitly tested.
Under 'natural shelter driven occupancy', we would predict that artificial shelter occupancy de creases as natural shelter density increases (Fig. 1B) because the probability of animals encountering, and thus using, artificial shelter should decrease as the total number of shelters in an area increases (Lindenmayer et al. 2003). A key assumption of 'natural shelter driven occupancy' is that population density has a negligible effect on shelter occupancy, thereby allowing for the scenario where population density is high but shelter occupancy remains low (i.e. shelter occupancy does not function as a proxy for population density).
The 2 hypotheses described above can be merged and extended to yield the hypothesis of 'natural shelter per capita driven occupancy'. Under this hypothesis, we would predict that animals are most likely to use artificial shelter when natural shelter is limited, where limitation might be due to low abundance of natural shelter, high conspecific competition for natural shelter, or both. As a result, artificial shelter use should increase as the number of natural shelters per capita decreases (Fig. 1C). A key assump tion of this third alternative is that natural shelter density mediates the relationship between population density and artificial shelter occupancy, also allowing for the scenario where population density is relatively high but artificial shelter occupancy remains low.

Study area
Our research took place in 6 stream reaches (hereafter R1-R6) located in 3 streams (n = 1-4 reaches stream -1 ) in the upper Tennessee River basin in southwest Virginia, USA. We selected our stream reaches because pilot data suggested they represented a relatively wide gradient of hellbender density and natural shelter density. Stream reaches had an average water depth of 44 cm (range: 11-125 cm) and an average wetted width of 13 m (range: 8-17 m). In order to standardize the area sampled by each artificial shelter array, we defined reach length (range: 272-384 m) such that each reach consisted of a wetted spatial extent of approximately 5040 m 2 . We considered all wetted portions of the stream channel in our calculation of spatial extent, regardless of whether natural shelter was present.

Target species
The eastern hellbender Cryptobranchus alleganiensis alleganiensis is a large (up to 74 cm), longlived (25+ yr), fully aquatic amphibian (Taber et al. 1975). Hellbenders are typically associated with cool, fast-flowing streams and are benthic habitat specialists that heavily depend on rocky crevices for daily shelter and nesting. Sub-adults and adults exhibit high site fidelity to specific stream reaches and rock cavities (Nickerson & Mays 1973). Sexual maturity is reached at around 6-8 yr or approximately 300 mm total length (Peterson et al. 1988). Spawning occurs annually (Smith 1907) in late summer through early fall, when females may deposit 500+ eggs (Topping & Ingersol 1981) in a nest cavity that is guarded by a single male for several months until larvae emerge. Larvae and juveniles are rarely encountered and thus much of their ecology is poorly understood.

Establishing artificial shelter arrays
We deployed artificial shelter arrays (n = 30 shelters array -1 ) across 6 stream reaches, yielding 180 shelters in total. We constructed artificial shelters by hand according to the 'modified boot design' of Briggler & Ackerson (2012). We built internal frames out of galvanized hex-mesh and hardware cloth and covered each frame with an equal-part mixture of sand, Portland cement, and quick-setting concrete. Each shelter consisted of a rectangular cavity with a single entrance via a tunnel and a removable lid to allow researchers access to the main chamber (Fig. 2). Chambers averaged 40 cm long (range: 29-49 cm), 38 cm wide (range: 23-48 cm), and 14 cm high (range: 7-23 cm). Tunnels averaged 24 cm long (range: 20-34 cm), 11 cm wide (range: 8-18 cm), and 10 cm high (range: 6-14 cm). Unlike Briggler & Ackerson (2012), we outfitted most shelters with a PVC cleanout plug (5 cm diameter) in one corner of the 170 Fig. 2. Artificial shelters installed into natural stream habitat and occupied by hellbenders Cryptobranchus alleganiensis main chamber to function as a camera port and retained a solid floor in each shelter. Once cement dried, we drilled 3 holes (4-5 mm diameter) through the upper rear wall of each shelter to facilitate water circulation through the main chamber and soaked shelters outdoors in cattle tanks for ≥30 d prior to deployment. Our initial shelter design relied on the use of boulders and cobble to anchor lids once shelters were deployed. However, due to relatively high rates of lid displacement during seasonal peaks in stream discharge early in our study (see Section 3.1) we retrofitted all shelter lids with a locking mechanism (see  in May 2015, about two-thirds of the way into our study. We stratified shelters evenly throughout each stream reach such that shelter density in each array was approximately 0.6 shelters per 100 m 2 of wetted in-stream habitat. We assumed that each array effectively sampled the entire extent of wetted habitat in the stream reach (ca. 5040 m 2 ). Due to the time required to construct shelters and our desire to fieldtest artificial shelters before deploying 6 full arrays, we staggered shelter deployments. We deployed an initial 30 shelters between June and September 2013 in 2 stream reaches (see Fig. 3) and deployed the final 150 shelters about 1 yr later (June-August 2014). To install shelters, we excavated a hole in existing substrate, positioned the shelter such that the tunnel opening was oriented 75-180° from oncoming current, surrounded the periphery of the shelter with gravel, sand and pebbles, and lined the tunnel and chamber with sand and gravel to a depth of approximately 2 cm. Finally, we placed one or more large rocks on top of each shelter for concealment and to anchor lids (Fig. 2). We used numbered aluminum tags to uniquely identify shelters. When a shelter was relocated for any reason, we assigned it a new number (hereafter, shelter ID). Thus, shelter ID reflected a unique combination of shelter and location.

Estimating hellbender density
We estimated hellbender density in each stream reach using robust design mark-recapture surveys (Pollock 1982). We conducted these surveys as part of a separate study to understand effects of land cover on hellbender populations. Detailed field and analytical methods are described by Bodinof Jachowski & Hopkins (2018). Briefly, we used rock-lifting while snorkeling to exhaustively search a randomly selected sub-reach within each reach; where each subreach represented one-third of the area covered by a given artificial shelter array. We did not survey the entire 5040 m 2 covered by each array because we wanted to limit microhabitat destruction caused by rock-lifting. Habitat appeared similar throughout each array, leaving us confident that hellbender density in randomly selected sub-reaches was representative of density throughout the extent of wetted stream (5040 m 2 ) covered by the corresponding array. We conducted robust design surveys such that 2 consecutive surveys of each sub-reach occurred in both 2014 and 2015. During surveys, we marked hellbenders using passive integrated transponder (PIT) tags injected subcutaneously along the lateral tail musculature and used capture histories for individuals to estimate abundance (while accounting for imperfect detection) in each reach for 2014 and 2015, separately. Larvae and young juveniles (<131 mm total length) were rarely encountered during these surveys and thus we ignored them when estimating hellbender abundance. As a result, our density estimates are specific to the pooled group of sub-adults and adults (hereafter individuals or hellbenders).
We used the extent of wetted stream in a sub-reach (ca. 1680 m 2 ) to convert our abundance estimates to density estimates. In 2014, density among reaches ranged from 0.02 ± 0.002 to 0.29 ± 0.015 ind. per 10 m 2 . In 2015, density among reaches ranged from 0.04 ± 0.003 to 0.31 ± 0.015 ind. per 10 m 2 . To place these values in context, our lowest density reach was estimated to contain 4-6 (2014: N = 4.12 ± 0.38; 2015: N = 6.20 ± 0.50) ind. per 100 linear meters of stream, while our highest density reach was estimated to contain 47-51 (2014: N = 47.66 ± 2.43; 2015: N = 51.98 ± 2.56) ind. per 100 linear meters. For any given reach, density estimates were highly similar be tween years (see Bodinof Jachowski & Hopkins 2018). For simplicity, we averaged mean density estimates for 2014 and 2015 to obtain a single mean estimate of hellbender density for each reach to use in our analysis addressing shelter occupancy (see Table 1).

Estimating natural shelter density
To estimate the proportion of each reach characterized by natural shelter, we characterized substrate composition in each reach using a modified Wolman Pebble Count with 100 random observations (Wolman 1954). We conducted our pebble count in the same randomly selected sub-reach used to estimate hellbender density. Again, habitat quality appeared consistent throughout each reach, leaving us confident that substrate composition in randomly selected sub-reaches was representative of the broader reach over which an array was deployed.
We used the proportion of each reach characterized by boulder (rocks with a secondary axis > 256 mm) and bedrock as a proxy for natural shelter density in each reach. We used boulder and bedrock cover to define natural shelter because the majority of hellbender captures (1006 of 1067 or 94%) recorded from our study system between 2007 and 2015 were associated with bedrock or boulders. Among our 6 stream reaches, percent natural shelter spanned about an order of magnitude (range: 5-51% of observations from each sub-reach). Notably, hellbender density and the percent of each reach characterized by boulder and bedrock were not noticeably correlated (see Table 1).

Artificial shelter surveys
We used snorkeling to inspect shelters approximately monthly from deployment through early September 2015. On each visit, we classified the status of each shelter as either (1) 'available', (2) 'unavailable', or (3) 'censored'. We classified shelters as available if the shelter was in place and the tunnel was at least partially open. Each time a shelter was classified as available, we recorded whether the lid was on or off and whether the shelter was occupied by a hellbender. We classified shelters as unavailable if the shelter had been dislodged since the last time it was surveyed or if the tunnel was entirely blocked by substrate. When possible, we restored unavailable shelters to an available condition by clearing any tunnel blockage and re-installing dislodged shelters that were still intact in a more sheltered location within ~15 m of the original location. When shelters were destroyed or permanently lost following dislodgment, we removed them from the study but retained all observations up to that point for our analysis. Finally, we classified shelters as censored if we were unable to inspect the shelter due to poor visibility or high flows.
When we encountered occupied shelters, we removed hellbenders by hand, weighed, measured, and sexed them, permanently marked them with a PIT tag (model HPT 12; Biomark), and released them in the shelter where they were captured. The only exception to our processing protocol occurred during our final monthly survey (September 2015), when we recorded occupancy status of each shelter but did not capture or process hellbenders to avoid disrupting breeding.

Statistical analyses
We used a multi-model approach and Bayesian mixed effect logistic regression models to examine relative support for our hypotheses regarding factors influencing shelter occupancy. We used observations of individual shelters as the sampling unit in our analysis, and only considered observations when a shelter was deemed available. Our analysis assumed that the state of each shelter (occupied/unoccupied) on each occasion was determined without error since shelters were searched exhaustively by tactile and visual means on each occasion. Importantly, our analysis allowed for the state of each shelter to change between successive surveys. We considered hellbender density and proportion of the site characterized by natural shelter as drivers of among-site variation in shelter occupancy and considered day of year, days since deployment, and lid status as covariates explaining temporal (i.e. within-site) variation in artificial shelter occupancy (see Table S1 in Supplement 1 at www. int-res. com/ articles/ suppl/ n041 p167 _ supp. pdf). Prior to model fitting, we scaled continuous covariates to have a mean value of zero and screened them for collinearity.
We used a candidate set of 7 models representing our hypotheses about factors driving occupancy of artificial shelters by hellbenders (see Table 2). To represent 'population density driven occupancy', we used a model that included only population density. To represent 'natural shelter driven occupancy', we used a model that included only natural shelter. To represent 'natural shelter per capita driven occupancy', we modeled shelter occupancy as an interaction between population density and natural shelter. We considered models that assumed all of the variation in shelter occupancy was driven by spatial factors, a model that assumed all variation was driven by temporal factors, and several models that included both spatial and temporal factors (see Table 2).
For temporal factors, we used a quadratic (x + x 2 ) form of day of year to represent our prediction that a peak in artificial shelter occupancy would coincide with the peak in hellbender movement just prior to spawning (late August -mid September). We used a pseudothreshold (ln [x]) form of days since shelter deployment to represent our prediction that artificial shelter occupancy rates would gradually increase post-deployment up to some threshold and thereafter remain stable. We included lid status to represent our hypothesis that hellbenders would be less likely to use shelters when the lid was ajar or off. Finally, to account for temporal autocorrelation in observations from the same shelter, we included an additive random effect term (ε[ID]) that assumed a first-order autoregressive (AR1) correlation structure to our data. This term as sumed that the status of a shelter at time t was dependent on the status of the same shelter at time t -1. Correlation among sequential observations from the same shelter was indicated in a preliminary analysis by the fact that all models in cluding the AR1 term outranked models that lacked the term.
We coded the outcome of each observation as binomial response, where a '1' indicated a shelter was occupied and '0' indicated a shelter was not occupied. We used Markov chain Monte Carlo (MCMC) and a Gibbs sampler in JAGS (Plummer 2013) with the package 'runjags' (Denwood 2016) in R v.3.2.1 (R Core Team 2013) to obtain posterior distributions for model parameters (see Supplement 2). We ran 3 chains for 10 000 iterations with a burn-in of 50 000. Trace plots of each parameter and Gelman-Rubin convergence statistics (Gelman et al. 2013) were used to determine convergence of chains (R < 1.1). We ranked candidate models according to Watanabe-Akaike's information criterion (hereafter, WAIC; Watanabe 2013). Similar to other model ranking criteria, WAIC represents a measure of model fit corrected by a penalty for model complexity, and models with lower scores outrank those with higher scores. However, unlike other criteria, the WAIC assesses fit based on full posterior predictive distributions rather than point estimates and thus is recognized as the only fully Bayesian model ranking criterion (Hooten & Hobbs 2015). We considered the hypothesis represented by our top-ranked model to be the best supported in our candidate set. We report mean estimates of parameter coefficients and mean estimated effects of covariates in our top-ranked model and their associated 95% credible intervals unless otherwise noted. We interpreted there to be strong evidence of effects from a covariate if the 95% credible intervals for the associated parameter coefficients did not overlap zero.
While WAIC can be used to rank relative support among candidate models, it provides little information regarding model performance. To assess performance of our top-ranked model, we conducted posterior predictive checks using Bayesian p-values (Conn et al. 2018). Briefly, posterior predictive checks entail a comparison between data simulated under the model to data used to fit the model. For each MCMC iteration, a goodness of fit statistic was used to assess divergence of the true data from their expected value according to the model. We compared a χ 2 test statistic computed from the actual data to a χ 2 test statistic computed from the replicated data to obtain a Bayesian p-value, where pvalues close to 0.5 indicate a well-fitting model and values close to 0 or 1 indicate poor fit (Kéry & Schaub 2011).

Rates of shelter availability, unavailability, and censoring
We recorded a total of 2060 observations from artificial shelters over the course of our study. Because we assigned new identities to shelters that were relocated during our study, our observations were distributed among 249 uniquely identified shelters, where each shelter was surveyed on an average (± SE) of 8.2 ± 0.28 occasions (range: 1-20 occasions).
Shelters were available on 87% of all survey occasions (range: 77-95% of occasions array -1 ; Table 1). Among occasions where a shelter was deemed available, lids were ajar/off 8% of the time (range: 5-17% of occasions array -1 ; Table 1). Disturbance to lids was about 5 times more common prior to installation of lid-locking mechanisms (lids missing on 144 of 1372 or 10% of occasions) than after (lids missing on 10 of 429 or 2% of occasions). Shelters were deemed un available on 12% of all occasions (range: 5-21% of occasions array -1 ; Table 1). Tunnel blockage by sand and gravel explained 81% of unavailable cases (range: 66-100% of unavailable occasions array -1 ). Dislodgement explained 19% of unavailable cases (range: 0-34% of unavailable cases array -1 ) but was rare overall (2% of all occasions; range: 1-6% of occasions array -1 ; Table 1). The majority of dislodgements occurred during heavy rain events in November 2014 (n = 18 dislodgments) and May 2015 (n = 21 dislodgements). In 19 cases (range: 0-10 cases array -1 ; Table 1), dislodged shelters were permanently de stroyed. Anecdotally, dislodgment of shelters tended to be most common in reaches with ex tensive amounts of exposed bedrock and in reaches where stream discharge exhibited acute spikes following rain events (i.e. flashier flows). Shelters were censored on approximately 1% of occasions (range: 0-2% of occasions array -1 ), typically due to high turbidity, flow, and water depth.

Characteristics of hellbenders using artificial shelter
Collectively, we detected 161 unique hellbenders (62 F, 83 M, 16 unknown sex) from our arrays (range: 6-69 unique ind. array -1 ; Table 1) over the duration of our study. There was no evidence that one sex was consistently more likely to use artificial shelter than the other. In total, 77 individuals were captured only once, 35 were captured twice, and 44 were captured ≥3 times from shelters. One female was captured on 18 occasions from only 2 shelters located about 10 m apart. Hellbenders occupying shelters were all mature adults (≥291 mm total length), except for larvae that originated from nests established in shelters. We detected 33 nests in our arrays over the 3 breeding seasons that our study encompassed.

Spatial patterns in shelter occupancy
Hellbenders used artificial shelters in all 6 arrays, but shelter occupancy varied considerably among arrays (range: 7-48% of available occasions array -1 where characterized by occupancy; Table 1) and over time (Fig. 3). Among the hypotheses we considered to explain among-site variation in artificial shelter occupancy, 'natural shelter per capita driven occupancy' was the best supported (Table 2). Our topranked model included an interaction between population density and natural shelter and included all 3 temporal covariates. Our posterior predictive check (Bayesian p = 0.57) indicated that the top-ranked model fit our data well.
Model coefficients and their credible intervals provided strong evidence that shelter occupancy varied among reaches in response to population density, natural shelter, and their interaction (Figs. 4 & 5). Our top-ranked model indicated that the direction and intensity of the association between population density and shelter occupancy was heavily dependent on the density of natural shelter in the site. For example, shelter occupancy was estimated to double with every 0.10 unit increase in individuals per 10 m 2 when natural shelter covered just 5-10% of a reach (Fig. 5). In contrast, occupancy was estimated to increase by only 10-15% for every 0.10 unit increase in population density when natural shelter covered 16-20% of a reach, and was estimated to remain consistently low (0-10%), regardless of hellbender density, when natural shelter characterized > 20% of the reach.

Temporal patterns in shelter occupancy
Credible intervals for model coefficients indicated strong evidence that shelter occupancy rates varied in response to day of the year, days since shelter deployment, and whether or not a shelter lid was intact (Fig. 4). Across all reaches, the probability of a shelter being occupied was highest each year from early June through August (Fig. 6A). The probability of a shelter being occupied also increased gradually throughout the first 365 d post-deployment (Fig. 6B); such that season-specific rates were higher in year 2 relative to year 1. There was considerable uncertainty whether shelter occupancy continued to increase after 1 yr post-deployment (Fig. 6B). Odds ratios (e β ) indicated that absence of a lid from an otherwise available shelter decreased the odds of occupancy by around 85% but did not entirely eliminate the probability that hellbenders would use a shelter (Fig. 4).

DISCUSSION
Across space, we found that artificial shelter occupancy was highest in stream reaches where the ratio of hellbender density to natural shelter cover was highest. Thus, our findings were most consistent with the hypothesis of 'natural shelter per capita driven occupancy', and suggest that hellbenders are most likely to use artificial shelter when natural shelter is most limited, where limitations may be driven by either a low abundance of natural shelter items or relatively high conspecific competition for natural shelter. As a result, shelter occupancy rates did not function as a reliable index of relative abundance for the 6 stream reaches that we considered.
Given the ecology of hellbenders, it is not surprising that artificial shelter was most readily occupied in areas where population densities were high and natural shelter was relatively rare. Hellbenders are benthic habitat specialists and their dependence on large rocks and bedrock crevices for daily shelter, nesting, and to support a prey base is well documented (Nickerson & Mays 1973, Bodinof et al. 2012a. Adults are typically solitary (but see Jachowski & Hopkins 2013) and aggressively compete for shelters and nest sites, especially near the time of spawning (Smith 1907). For all these reasons, the availability of cover is suspected to function as a common limiting factor for hellbender populations and may have been a limiting factor for hellbender populations that used our artificial shelters most frequently. However, our findings suggest that 176 Table 2. Candidate set of models representing alternative hypotheses regarding the factors driving hellbender occupancy of artificial shelter, arranged from most to least supported according to the Watanabe-Akaike information criterion (WAIC). β: model coefficient; ε: error term to account for temporally autocorrelated repeated measures on unique shelters; ID: factor corresponding to unique shelter identity Rank Hypothesis Model WAIC ΔWAIC 1 Natural shelter per capita driven occupancy with temporal effects β 0 + β 1 × density + β 2 × natural shelter + β 3 × (density × natural shelter) + β 4 × ln(days deployed) + β 5 × day of year + β 6 × day of year 2 + β 7 × lid + ε(ID) 1398 0 2 Population density driven occupancy with temporal effects β 0 + β 1 × density + β 2 × ln(days deployed) + β 3 × day of year + β 4 × day of year 2 + β 5 × lid + ε(ID) 1407 9 3 Natural shelter driven occupancy with temporal effects β 0 + β 1 × natural shelter + β 2 × ln(days deployed) + β 3 × day of year + β 4 × day of year 2 + β 5 × lid + ε(ID) natural shelter was not a limiting factor for hellbender density in several of our study reaches. For example, 2 of our low-density (< 0.1 ind. per 10 m 2 ) populations inhabited stream reaches with 3-10 times higher levels of natural shelter than what was observed in the reach hosting our highest population density (Table 1). This is a particularly salient point, because it emphasizes the potential for ex tremely low rates of artificial shelter use in stream reaches where use of such structures as a man agement tool is most desired (i.e. those harboring low-density hellbender populations). More work is needed to assess the degree to which natural shelter density might mediate the relationship between hellbender density and artificial shelter occupancy in other portions of the species' range. While our study indicated that the relationship between hellbender density and artificial shelter occupancy depended on the extent of bedrock and boulder within a stream reach, a more recent investigation into patterns of shelter occupancy using 10 shelter arrays (including 5 arrays from the current study) failed to detect an effect of natural shelter and concluded that artificial shelter occupancy is driven primarily by hellbender density (Button 2019). We suspect that the contradiction between our findings and those of Button (2019) can be largely explained by differences in the way natural cover was defined in each study and differences in study reaches considered. It is worth noting that the nature of the interaction between hellbender density and natural cover in our top-ranked model indicated that shelter occupancy was positively associated with hellbender density when boulder and bedrock covered ≤20% of a stream reach, but uncorrelated with hellbender density when boulder and bedrock was more common. Therefore, our results suggest that the power to detect substrate effects on artificial shelter occupancy can depend heavily on substrate composition in the suite of reaches considered.
The ecological processes driving the spatial variation in shelter occupancy that we report here remain unclear. Low rates of artificial shelter occupancy in areas where natural shelter was relatively abundant could have been due to hellbenders perceiving artificial shelter as a less-viable alternative to natural shelter or, alternatively, could have simply been due to a lower frequency of encounters between hellbenders and artificial shelter. A previous study found that hellbender home range size and the average dis- tance of daily movements was negatively correlated with boulder density (Bodinof et al. 2012b). If similar relationships between hellbender movement and natural cover occurred in our system, the high rates of artificial shelter use in areas where natural cover was relatively rare might be at least partially attributable to relatively high rates of hellbender encounters with artificial shelter during typical movements. Temporal variation in artificial shelter occupancy within a site, like what we observed for hellbenders, is a common observation among herpetofauna (Grant et al. 1992, Monti et al. 2000. Seasonal variation in artificial shelter occupancy is common among birds in particular and has been associated with variation in resource demand and availability, population fluctuations, or some combination of the 2 (Mainwaring 2011). We suspect that hellbender behavior related to spawning is the most likely explanation for the increase in artificial shelter occupancy that we observed each summer. Competition for shelter sites and aggressive interactions between hellbenders increase during the weeks leading up to spawning (Smith 1907), suggesting that many individuals may become displaced during summer months and may be searching for cover.
We attribute the increase in shelter occupancy over time in our study to additional opportunities hellbenders had to encounter new shelters. However, it is also possible that shelters became more attractive to hellbenders as they weathered. The gradual increase in shelter occupancy rates during the first year postdeployment might also be attributable to the hellbenders' activity patterns. We observed several hellbenders moving along the stream bottom in April and May and therefore suspect that spring foraging movements may have been particularly important for facilitating discovery of shelters that were installed 6-12 mo prior. The effect of time since deployment that we describe may at least partially explain why Messerman (2014) failed to detect any hellbenders using 54 artificial shelters that were installed into North Carolina streams, given that shelters were deployed in August and only monitored for 5 mo.
Our study is the first of which we are aware to explicitly evaluate support for 'population density driven occupancy' in a herpetofauna system. Under this hypothesis, shelter occupancy or counts of individuals using artificial shelter should increase proportionally with abundance, which is implicitly assumed in many studies that utilize artificial cover searches to monitor herpetofauna (Hesed 2012). For example, uncorrected count data from artificial cover arrays have commonly been used to investigate eco-logical relationships for herpetofauna and evaluate effects of alternative habitat management scenarios on local abundance (Degraaf & Yamasaki 1992, Morneault et al. 2004, Semlitsch et al. 2007, Mathewson 2009, Tilghman et al. 2012). However, our work showed that artificial shelter occupancy failed to reflect variation in relative hellbender abundance among the 6 reaches considered in our study, which has important implications for future studies. Additional work is needed to investigate how often the assumption of 'population density driven occupancy' might be violated for other herpetofauna.

MANAGEMENT AND CONSERVATION IMPLICATIONS
Our study highlights both the potential value and important limitations of using artificial shelters to monitor hellbenders. Our findings tentatively suggest that shelter occupancy should function as a rough index of relative abundance when the following criteria are met: (1) abundance of natural cover is similar among stream reaches involved, and fewer than 20% of observations in a Wolman Pebble Count fall within the categories of boulder or bedrock in any reach; (2) shelter surveys are conducted at the same time of year at each site, preferably mid-summer; and (3) estimates of shelter occupancy are adjusted to account for shelter availability (i.e. dislodged shelters or those blocked by sediment are not included in the estimate). Our findings also suggest that shelter occupancy rates may be an informative proxy for temporal variation in relative abundance within a single site once arrays have been established, which recent work suggests make take 2-3 yr (Button 2019). Finally, given that artificial shelters facilitated a large number of hellbender detections and evidence of breeding success in some stream reaches, they may be particularly valuable for facilitating long-term capture-mark-recapture work or reproductive ecology research. When applied effectively, artificial shelters can facilitate minimally invasive monitoring and research, reduce personnel needs during field surveys, and offer a method for standardizing survey effort over time. In contrast, one important limitation of artificial shelters is their inability to facilitate detection of immature hellbenders, except for eggs and nestling larvae that originated from nests established in the shelters themselves. We suspect that adjustments to shelter design and/or the microhabitat where shelters are deployed might improve attractiveness to sub-adults and juve-nile animals and encourage work to explore alternative shelter designs for immature age classes. Secondly, we show that raw rates of artificial shelter occupancy can be driven by a more complex suite of factors than hellbender abundance alone. Because we found that shelter occupancy rates were not a reliable proxy for relative abundance among the 6 reaches considered in our study, we urge caution against using raw shelter occupancy rates to infer patterns of relative abundance among reaches when substrate composition is highly variable.
While artificial shelters can be useful for some applications, deploying arrays requires overhead costs associated with time and labor for shelter construction, installation, and maintenance. Aside from materials (~$30 US per shelter), the largest cost associated with construction as per Briggler & Ackerson (2012) is labor. We did not explicitly quantify person-hours spent constructing shelters but estimate that an experienced individual can construct a shelter frame in about 1 h and complete concrete work for the shelter body and lid in about 3 h. However, efficiency of the construction process generally improves as shelters are produced in bulk. Shelters weigh ≥50 lb (~23 kg), which can make deployment physically challenging, especially in remote locations. Perhaps most importantly, regular maintenance is critical to ensure that shelters remain in place, intact, and unblocked by sediment. Encouragingly, a recent study demonstrated that a slight modification to shelter and lid de sign paired with more targeted placement of shelters in certain microhabitats virtually eliminated lid displacement and shelter dislodgement events in streams throughout southwest Virginia (Button 2019). Sediment buildup in tunnels, however, may be more difficult to overcome in many streams. Evidence reported here is consistent with more recent work (Button 2019) and suggests that maintenance every 30-40 d and targeted visits following high flow events are necessary to ensure that most shelters re main at least partially open. While Button (2019) found that slight modifications in shelter orientation relative to water flow can help reduce frequency of tunnel blockage, the study also suggests that sediment loads in some streams may be high enough to preclude the practicality of shelter deployment entirely. For all these reasons, we encourage others to use pilot deploy ments to validate suitability of stream reaches before deploying extensive shelter arrays.