Movement, demographics, and occupancy dynamics of a federally-threatened salamander: evaluating the adequacy of critical habitat

Study area

Our study was conducted in the Bull Creek watershed, Travis County, Austin, Texas. The spring population case study was conducted near Lanier Spring, on the Balcones Cayonlands Preserve, a local preserve system dedicated for the protection of endangered avian and karst invertebrate fauna. Lanier Spring emerges from alluvial deposits and discharges into Bull Creek several meters from the spring opening. The presence of potential habitat upstream and downstream of the Lanier Spring confluence in the form of shallow, clear water with abundant rock cover (Bowles, Sanders & Hansen, 2006) made it an ideal site to study aspects of population biology and individual movement of E. tonkawae in the context of a spring-stream habitat gradient.

Headwater streams in Bull Creek are spring-fed, gaining streams (i.e., discharge increases with stream length) with flows that undulate between the surface and subsurface, particularly during dry periods. We conducted occupancy surveys within five tributaries of Bull Creek, including those whose catchments were predominantly within the Balcones Canyonlands Preserve (Mainstem/Trib. 8 & Trib. 7) and streams surrounded mostly by urban development (Trib. 4, Barrow Hollow, & Trib. 2).

Study 1: movement and distribution around a spring

We surveyed Lanier Spring every two weeks from January through April 2013 and delimited nine 5 m survey sections around the spring outlet (four upstream, four downstream, and one at the spring outlet) every 15 m. The wetted channel boundaries determined the survey width, as each 5 m section was surveyed from bank to bank. Based on previous capture-recapture work at this site (Bendik et al., 2013b), we suspected the number of salamanders available for capture (and recapture) would be enough to detect upstream or downstream movements of individuals between sections. Detection of movement was restricted to the stream channel because this species is aquatic and overland movements are not possible. Sections were exhaustively searched for salamanders by flipping cover objects such as rocks and leaf litter. Once detected, salamanders were caught using dip nets and photographed on a standardized 0.5 cm grid. Measurements of body length were obtained from photographs (following Bendik & Gluesenkamp, 2013) using Image J (Rasband, 1997). We cropped photos to include only the head and used Wild-ID (Bolger et al., 2012) for individual identification based on pigmentation patterns. A previous study demonstrated that identification errors using this technique were lower (false rejection rate = 0.76%) compared to visible implant elastomers (VIE; false rejection rate = 1.90%) for E. tonkawae (Bendik et al., 2013b). Therefore, we were confident in our ability to correctly identify individuals. In addition to documenting intra-seasonal movement, we were interested in the presence and location of VIE individuals (N = 1, 115 marked from 2007 to 2011) from a prior study at this site. Their presence in areas beyond the spring outlet may indicate dispersal from the spring.

To compare abundances in and around the spring, we estimated superpopulation sizes for each 5 m section during the four-month sampling period using the POPAN model in program MARK 8.0 (White, 2015). The superpopulation is the abundance of all animals associated with the study system (including those within and outside of the sample site). We used model-averaged values to estimate superpopulation sizes at each section. Our model set included combinations of either time-varying or constant detection probability (p) and apparent survival (ϕ); we fixed probability of entry to vary by time (but not section) and superpopulation size by section. We also investigated the effect of body length on p and ϕ.

To calculate the probability of movement while accounting for detection error, we used a multi-state model with live-recaptures in program MARK 8.0 (White, 2015). This is an extension of the Cormack–Jolly–Seber model for calculating apparent survival and detection probability in open populations that allows for movement between strata with transition probabilities (Hestbeck, Nichols & Malecki, 1991; Brownie et al., 1993). We considered models whereby salamanders could either remain in their current position, or move to a new section. Transition probabilities estimate the rate movement between sections, an average distance of 20 m. Movement probability was calculated as twice the transition probability, since transitions in both directions indicate movement. Our model set included all combinations of time variance for apparent survival (ϕ), detection (p) and transition probabilities (ψ), but we did not incorporate differences among states. We also tested for the influence of body length on apparent survival, detection, and transition probability. Program U-CARE (Choquet et al., 2009) was used to estimate the $\hat{c}$ ratio, a measure of goodness-of-fit, for both the multi-state and POPAN models.

Study 2: habitat-occupancy dynamics

Starting at the downstream confluence below the nearest known E. tonkawae locality, we surveyed 97 sites in 2013, and surveyed an additional 30 sites in 2014 and 2015, for a total of 127 sites (encompassing over 12 km of linear stream channel). Of those, 118 contained wetted habitat at some point during our sampling. The distribution of sites was nearly evenly split between those occurring within preserves (n = 57) and those in urbanized tributaries (n = 61). With the exception of avoidance of deep pools (>0.6 m) such as man-made impoundments, we selected sites in a systematic, random fashion. Site selection started at a random point 0–100 m from the confluence and proceeded upstream, while maintaining equal distances among sites. Because we were targeting an even distribution of sites among tributaries and urban/non-urban catchments, sites in Barrow Hollow, Trib. 2, Trib. 4 and Trib. 7 were spaced 70 m apart, while those in Mainstem/Trib. 8 were spaced 140 m apart due to the larger size of that tributary. Each site was 10 m in length (from downstream to upstream) and the wetted channel boundaries determined the width of each site. Sites were surveyed each year from 2013 through 2015; our survey seasons generally started in late March and ended in early May. Because sampling occurred once per year, we use the term “year” to refer to occupancy seasons. Two experienced observers simultaneously performed a time-limited search removing cover objects for 5 min at each site, three times per year, at a minimum interval of one week. During our surveys, we documented new spring localities where we observed water issuing from conduits within or near the stream channel and included previously documented springs (City of Austin, 2015) to compare to the distribution of occupied sites.

We also collected a suite of environmental data at each site, which was reduced to a smaller set of covariates based on our specific hypotheses of interest (see below) and an assessment of multicollinearity using eigenvector analysis. Scatter plots among candidate predictive variables and an assessment of multicollinearity are included in Supplemental Information 1. We conducted one-way ANOVAs among continuous predictor variables and site groupings (tributary or urbanization) and Chi-squared tests among categorical variables and groupings. The program R (R Core Team, 2014) was used for data preparation and preliminary statistical analyses. Significance was evaluated at α = 0.05.

We used a model that integrates both habitat and occupancy dynamics to fit detection/non-detection data and changes in habitat suitability (Mackenzie et al., 2011). We used a simplified version of this model, with only three possible states: suitable, but species undetected; suitable with species detected; or unsuitable and unoccupied. Because neotenic Eurycea are stream-dwelling, aquatic organisms, we classified sites that had flowing water during some portion of our sampling season as suitable (S) and those that were either dry during the entire period, or initially stagnant and later became dry, as unsuitable (U). In this context, “suitable” habitat refers to a minimum condition for the presence of an aquatic species, and beyond that does not imply that such habitat is necessarily of good quality. Following the notation of Mackenzie et al. (2011), the parameter π^[S] represents the probability of a site being suitable habitat during the first year. The probability of habitat changing state between years is represented by two parameters: ${\eta }_t^{\left[S,S\right]}$ represents the probability of habitat remaining suitable between year t and t + 1; ${\eta }_t^{\left[U,S\right]}$ is the probability of habitat transitioning from unsuitable to suitable between year t and t + 1. The probabilities of habitat to transition from suitable to unsuitable, or to remain unsuitable are 1 $-{\eta }_t^{\left[S,S\right]}$ and 1 $-{\eta }_t^{\left[U,S\right]}$, respectively. The probability the species was initially present at a site is represented by ψ^[S] (initial occupancy probability for unsuitable sites is fixed at zero). The probability a species colonizes a site between years t and t + 1 is described by two parameters, ${\gamma }_t^{\left[S,S\right]}$ and ${\gamma }_t^{\left[U,S\right]}$. The probability a species goes extinct at a site between years t and t + 1 is ${\epsilon }_t^{\left[S,S\right]}$. By definition, when habitat becomes unsuitable, colonization does not occur (${\gamma }_t^{\left[S,U\right]}={\gamma }_t^{\left[U,U\right]}=0$) and extinction does (${\epsilon }_t^{\left[S,U\right]}=1$). Detection probability of habitat state was assumed to be perfect and we also assumed that species occupancy had no bearing on habitat suitability. A primary assumption of the occupancy model we used required that both habitat suitability and occupancy remained constant within each sampling season.

With regard to habitat dynamics, we predicted that habitat suitability would be dependent upon the prior state because spring locations and other geomorphological factors that influence hydrologic periodicity are assumed to be relatively constant. We therefore compared models where habitat suitability exhibited a first-order Markovian pattern of state dependence (represented by two parameters, ${\eta }_t^{\left[S,S\right]}$ and ${\eta }_t^{\left[U,S\right]}$) to those where habitat suitability was independent of the previous state (${\eta }_t^{\left[S,S\right]}={\eta }_t^{\left[U,S\right]}$). Our first year of sampling was also the driest, so we expected that urban tributaries would have more wet sites than non-urban tributaries because of artificial groundwater recharge from irrigation and leaking treated water and wastewater infrastructure (Christian, Banner & Mack, 2011). To test for this effect, we considered models where the initial habitat state π^[S] was a function of the group variables tributary (TRIB) or urbanization (URB). Sites are nested within tributaries, and tributaries are nested within the urbanization classification. URB is a tributary-wide categorical variable based upon predominant land use within the catchment. Compared to those within preserves, urbanized tributaries are characterized by having markedly higher water conductivity (Supplemental Information 1 and Figure S4) and are surrounded by residential development.

We also tested the importance of TRIB and URB groupings with regard to site occupancy and colonization. The initial descriptive analysis we performed indicated both tributary and urbanization groupings could be explained by certain environmental covariates. Therefore, we modeled the effects of group and environmental explanatory variables on the parameters of interest separately.

To determine what habitat conditions were associated with site occupancy, we compared models where initial occupancy and colonization were dependent upon a combination of habitat covariates. We initially collected a suite of environmental variables and narrowed down this set to a candidate list of possible predictors based on biological hypotheses and elimination of redundant or collinear variables (see Supplemental Information 1). Because previous studies of central Texas Eurycea have demonstrated the importance of groundwater habitats, we were interested in variables that could be reliable predictors of groundwater influence at a given site. Seasonal variation in groundwater temperature is much less than that of surface water, so sites with lower variation in temperature within sampling seasons could indicate groundwater-influenced habitats such as groundwater springs and stream-fed “springs” (created when upstream flow sinks underground, reemerging downstream). We measured water temperature during each survey and computed temperature variation as the standard deviation of temperature (TS) within each year. Additionally, we also used the bankside presence of maidenhair ferns (MH), Adiantum sp., as potential indicators of groundwater or water permanence. These plants are commonly found near springs and seeps in Texas (Brune, 1981). We predicted low water temperature variation and fern presence would be positively correlated with salamander occupancy and colonization. Site occupancy within these stream reaches may also be predicted by other aspects of their habitat. Rock cover (RC) is an important habitat characteristic that is often correlated with presence and abundance of Eurycea (Bowles, Sanders & Hansen, 2006; Pierce et al., 2010; Diaz et al., 2015). RC was measured as the mean of three visual estimates of the percentage of rocks available (i.e., unembedded) for cover (>8 mm) at the substrate surface. Water depth (WD) may be important for exclusion of fish predators (although we purposely avoided sampling deep pools) so we expected it to negatively affect occupancy. However, WD may also influence dry risk and temperature variation as well. WD was taken as the mean of three evenly spaced measurements at the thalweg. Sites with high calcium carbonate (CA) cementation may not be preferred habitat (Sweet, 1982), lacking enough interstitial spaces for salamanders or their prey items. CA was visually estimated as being in one of four categories: (1) no carbonate build-up; (2) low to moderate build up (but <50% area); (3) moderate to high build up (>50% area); and (4) complete covering of area with carbonate and tufa. For analyses, these categories were collapsed into low (1 or 2) and high (3 or 4); 43% of observations had no carbonate cementation and only two observations were ranked 4. We incorporated yearly differences in TS, WD and RC in our analysis, but assumed that CA and MH did not change (we used the single, maximum value recorded among all years). Because missing covariate data are not allowed in occupancy models, we used the within-tributary average for 5 missing values of WD and 1 of TS (1% of measurements for those covariates) in lieu of excluding those sites from our analysis. Twelve survey events with available occupancy data were missing most covariate values due to having a single visit (e.g., temperature variation measurements required at least two visits), and were dropped from the analysis.

To understand how habitat instability influences occupancy dynamics, we compared models where colonization was dependent on the prior habitat state (${\gamma }_t^{\left[U,S\right]},{\gamma }_t^{\left[S,S\right]}$) or was random with respect to prior habitat (${\gamma }_t^{\left[U,S\right]}={\gamma }_t^{\left[S,S\right]}$). Estimates of ${\gamma }_t^{\left[U,S\right]}$ are a measure of the ability to colonize a previously dry site, which may indicate migration between stable spring habitats and/or an ability to take advantage of habitats throughout the stream during wetter years. Our sampling occurred over a period of progressively wetter conditions, so we were able to determine whether dry sites were eventually colonized or not. Covariates were mirrored among ψ and γ_t to simplify the model set, and because we were interested in the set of covariates that most consistently predicted whether a site could be occupied or not. For example, if we included WD (from year 1) as a covariate on ψ^[S], γ₁ and γ₂ were also modeled with WD (corresponding to year 2 and year 3, respectively). The effect of environmental variables on colonization was assumed to be consistent regardless of the prior habitat state. Continuous covariates were standardized for analysis by subtracting the arithmetic mean from each value.

The program PRESENCE 9.7 (Hines, 2006) was used to fit models using maximum likelihood. We employed a multi-stage approach to model selection similar to those used in Dugger, Anthony & Andrews (2011), Lee & Bond (2015), and MacKenzie et al. (2012), as follows: for each stage we held fixed a general model structure while varying the structure for the focal parameter(s); then, from each model set we determined the best model using small-sample Akaike’s Information Criterion (AICc) values, keeping the best structure for the focal parameter(s) constant for each subsequent step. We estimated the effective sample size as the average value of the total number of sites and the total number of individual surveys. First, we held habitat and occupancy parameters fixed with full time and state dependencies while we modeled detection as either constant, time variant among (but not within) years (YEAR), constant among tributaries (TRIB), or as additive, YEAR + TRIB. Proceeding from the best detection structure, we then tested habitat structures (as above), keeping occupancy time and state dependent. For the third step we compared models with and without state dependency (but with time variance), including a random occupancy model, imposing the constraint γ_t = 1 − ε_t, whereby changes in occupancy were independent of the prior state (MacKenzie et al., 2006). We then tested whether group effects TRIB or URB better explained variation in ψ^[S] and γ_t compared to the more general unstructured model. We included site-level covariates (TS, WD, RC, CA, MH) as explanatory variables in place of the group effects on ψ^[S] and γ_t to distinguish among hypotheses of habitat selection. We first determined the best structure of continuous, time-varying covariates (TS, WD, RC), to account for environmental variation in occupancy dynamics among years. Using the best combination of time-varying covariates, we then included variables MH and CA to determine if they added explanatory value to the model. We used the sum of AICc variable weights to quantify the relative importance of covariates. We model-averaged parameter estimates on the real scale and model averaged covariate coefficients (following Lukacs, Burnham & Anderson, 2010). Finally, because we used a multi-stage approach to model building, we re-tested detection and habitat model structures based on the top model to see if they were consistent with the initial decisions made in the model testing process.

Study 1: movement and distribution around a spring

Unique individuals were recaptured in all nine survey sections. Individuals moved frequently between sections along the spring-stream gradient. Throughout the four-month period we captured 215 unique individuals (>25 mm total length) and recaptured 81 at least once. Of the recaptured animals, 21 were identified in two different sections, implying a minimum movement distance of 15 m at least once. Fourteen moved one section (an average distance of 20 m), four moved two sections (40 m), two moved three sections (60 m) and one moved four sections (80 m). We did not observe strong directionality in movement: 10 moved downstream, 11 moved upstream. We observed VIE-tagged individuals in all but two sections, including both sections 80 m from the spring (Fig. 1). Incidentally, during the 2014 occupancy surveys, we observed a single VIE-tagged individual over 500 m upstream of Lanier Spring. This individual had been last observed in 2010.

The single-state goodness-of-fit analysis indicated a significant effect of transience, which is expected given the movement between sections we observed. There was also some evidence of overdispersion in the POPAN model ($\hat{c}=1.58$), therefore adjusted QAICc values were used for model selection and model-averaging. Model-averaged estimates of superpopulation size were greatest near the spring and 20 m downstream, although the section farthest downstream had the third highest estimate (Fig. 1). Small juveniles (<25 mm total length) were most abundant downstream, with the largest proportion (32%) occurring in the most distant section (Fig. 1), possibly accounting for the high superpopulation size estimate in that section. Body length was an important covariate in our model set, contributing to 98% of the QAICc weights (including any model where body length was a covariate), although the most optimal models only included body length on either ϕ or p, not both (Table 1). Model-averaged estimates of bi-weekly survival and detection probability based on mean body length (24.8 mm) ranged from 0.44–0.66 and 0.35–0.53, respectively.

For the multi-state capture-recapture model, goodness-of-fit was adequate and we did not find evidence of overdispersion ($\hat{c}=0.91$). Body length was also an important predictor of both survival and movement probability (Table 2), and was positively correlated with both (Fig. 2). The sum of AICc weights for models including body length as a covariate for φ and ψ was 0.92 and 0.74, respectively. Model-averaged estimates of bi-weekly survival and detection probability based on mean body length (24.8 mm) ranged from 0.61–0.66 and 0.11–0.55, respectively. Movement probability, corrected for imperfect detection, was 0.15 (SE = 0.07) for an average length salamander. A full list of model selection results and model-averaged parameter estimates for both capture-recapture analyses are provided in Supplemental Information 2.

Study 2: habitat-occupancy dynamics

We documented E. tonkawae nearby and well beyond previously recognized spring locations among 118 sites and five tributaries of Bull Creek (Fig. 3). The first year was the driest, and the proportion of wet sites progressively increased each year. These resulted in an overall positive trend in WD and a negative trend in TS each year (Table 3). We found significant effects of both tributary- and urbanization-based groupings for TS, RC (Table 3), MH, and CA (for results of statistical analyses, see Supplemental Information 1).

Based on the multi-stage modeling approach, the best model prior to including environmental covariates had a model weight of 0.94 and the following structure: tributary-dependent detection, p(TRIB); year- and state-dependent changes in habitat suitability (${\eta }_t^{\left[S,S\right]}$, ${\eta }_t^{\left[U,S\right]}$); initial habitat suitability, π^[S](URB), where sites grouped by urbanization category outperformed a tributary-based structure; Markovian site occupancy, dependent upon both the prior habitat and prior occupancy states for each year (${\gamma }_t^{\left[U,S\right]}$, ${\gamma }_t^{\left[S,S\right]}$, ${\epsilon }_t^{\left[S,S\right]}$); initial occupancy, ψ(TRIB), where sites grouped by tributary outperformed those grouped by urbanization category (evidence ratio for ψ (TRIB) vs. ψ(URB) models = 23). Based on this model structure, we compared models with environmental covariates on initial occupancy (in place of the group structure due to redundancies; Table 3) and colonization, resulting in two top models with a combined weight of 0.99 (Table 4). Re-testing all detection- and habitat-based model structures on the top model did not alter the model hierarchy. The following results are based on model-averaged estimates from the top two models.

Initial habitat suitability was higher in urbanized streams (${\hat{\pi }}^{\left[S\right]}=0.89$, SE = 0.05) compared to those in preserves (${\hat{\pi }}^{\left[S\right]}=0.63$, SE = 0.07). The probability of a site remaining wet from 2013 to 2014 was high (${\hat{\eta }}_1^{\left[S,S\right]}=0.96$, SE = 0.02), while the probability of a site transitioning from dry to wet during the same time period was lower (${\hat{\eta }}_1^{\left[U,S\right]}=0.26$, SE = 0.10). All sites were wet during the final year of sampling, resulting in both estimates of η₂ at the boundary of 1.00.

The best AICc (effective sample size = 620) time-varying covariate models of initial occupancy and colonization contained covariates TS and WD (Supplemental Information 2 and Table S13), each with 100% of the relative variable weights compared to 3% relative weight for RC (calculated from the model set including all combinations of TS, WD, and RC). Consistent with our hypotheses, estimates of TS and WD coefficients had a significant negative effect on initial occupancy with 95% confidence intervals excluding zero (Table 5). Inclusion of CA and MH categorical site-level covariates improved upon the continuous covariates model, decreasing the AIC value by over 8 units. MH was positively correlated with salamander presence, while CA was negatively correlated, though not significantly (Table 5).

In contrast to initial occupancy probabilities, the magnitude and direction of these covariate effects were inconsistent among colonization probabilities (Table 5). Previously unoccupied sites were approximately three times more likely to be colonized if the prior habitat state was dry, compared to sites that remained wet (e.g., ${\hat{\gamma }}_1^{\left[S,S\right]}=0.11$, SE = 0.06; ${\hat{\gamma }}_1^{\left[U,S\right]}=0.34$, SE = 0.06; Table 6). Extinction probabilities among sites that remained wet declined by approximately two-thirds for 2015 compared to the prior year (Table 6).

Odds and odds ratios are a convenient way to interpret the effects of covariates in occupancy models that utilize logistic regression (MacKenzie et al., 2006). Odds ratios based on the top two models did not differ in terms of the ranks of their magnitude (high values and low values were consistent among models), indicating that covariate effects were consistent among models. We therefore present results based on the model-averaged estimates. Sites with average WD and TS, but with MH were less likely to be occupied than unoccupied during 2013, with odds of 0.76:1 (calculated as e^−4.35 × e^4.07×1). Similarly, sites 1 standard deviation (SD) below average WD (SD = 7.5 cm) without MH were also less likely to be occupied (odds = 0.15:1). Sites 1 SD below average TS (1.62 °C), but with mean WD and no MH were slightly more likely to be occupied than unoccupied (odds = 1.37:1). However, when sites had both MH and low TS (1 SD below average), the overall odds of occupancy vastly improved at 80:1. Even sites with moderate, negative deviations from average covariate values and with MH dramatically improved the odds of occupancy (e.g.,−0.5 SD = 26:1 odds in favor of occupancy).

A full list of model selection results for occupancy analyses is provided in Supplemental Information 3.