Landscape- and local-scale habitat influences on occurrence and detection probability of Clark’s nutcrackers: Implications for conservation

Whitebark pine (Pinus albicaulis), a keystone species and an obligate mutualist of the Clark’s nutcracker (Nucifraga columbiana), is rapidly declining throughout its range. Evidence suggests this decline is leading to a downward trend in local nutcracker populations, which would in-turn decrease whitebark pine regeneration. Our objectives were to (1) evaluate temporal variation in nutcracker habitat use as a function of whitebark pine and Douglas-fir (Pseudotsuga menziesii) habitat, at local and landscape scales, (2) develop metrics for predicting when whitebark pine communities require intervention to sustain nutcracker visitation, and (3) test McKinney et al. (2009) and Barringer et al.’s (2012) models predicting nutcracker occurrence. Between 2009 and 2013, we carried out 3,135 audio-visual Clark’s nutcracker surveys at 238 random points in the southern Greater Yellowstone Ecosystem. Using Bayesian occupancy models and cross-product model selection, we evaluated the association between nutcracker occurrence and habitat variables during five stages of the nutcracker annual cycle, while accounting for imperfect detection. Nutcracker occurrence was most strongly associated with the presence of cone-bearing whitebark pine trees (rather than cone density) and the area of whitebark pine on the landscape. To promote a high, >75%, probability of occurrence at a site within the study area, we recommend a management plan that achieves a landscape composed of a minimum of 12,500–25,000 ha of cone-bearing whitebark pine habitat within a 32.6 km radius. Additionally, an optimal habitat mosaic includes moderate levels of Douglas-fir habitat. Models currently used to guide whitebark pine management strategies underpredicted nutcracker occurrence in our study area, suggesting these strategies may not be appropriate in the region. We cannot predict how this mutualistic relationship will change as the population density of each species shifts. We therefore suggest conducting periodic surveys to re-evaluate the relationship as the environment changes and management strategies are implemented.


Overview
In this appendix, we provide some additional details about the final occupancy model used in this study, i.e., the model used to produce the results presented in the main text. The inference objective here is to model Clark's nutcracker occurrence in each of s = 1, . . . , 5 biological seasons (BR: breeding season, ES: early and LS: late summer, FH: fall seed harvest, and PH: post-harvest).

Observation process -detectability
Rarely can a species be detected with absolute certainty and conducting multiple visits to sites in close succession in a period that is assumed to be closed provides the 'repeated visit' data required to estimate, and thus account for, imperfect detection of the true occupancy state.
Preliminary analysis showed no support for any covariate effects on detectability based on the 0.5 posterior model support threshold cut-off recommended by Barbieri and Berger (2004). In the absence of covariates, but to allow for variation in detectability, we modeled year-(t) and season-(s) specific detection probability as a normally distributed random effect with a single (global) mean (µ θ ) and standard deviation (σ θ ): Due to a lack of fit of the detection model (see below), we allowed for additional variation at the level of the unique site-season-year combination also as a normally distributed random effect with a mean of 0 and a standard deviation (σ ) to be estimated: Binary observations (y = 1 when detected and y = 0 otherwise) at site i, during survey j, in season s in year t, are then modeled as Bernoulli random variables with success probability p which is multiplied by z i,s,t , the latent occurrence state of site i, in season s, in year t, setting detection probability to p = 0 when sites are unoccupied and p when occupied: The detection parameters to be estimated, and the (uninformative) priors used in the Bayesian analysis, are provided in Table 1.

Ecological process -occurrence
The interest is in estimating ψ ist , the probability that site i is occupied in season s in year t, and specifically, to relate these probabilities to site specific covariates. Again, we conducted a preliminary analysis and, using the 0.5 posterior model support threshold cut-off recommended by Barbieri and Berger (2004), selected the final set of covariates which were: whitebark pine cones per hectare (cone crop: cc), the proportion of whitebark pine within 32.6 km radius (pw), and the proportion of Douglas-fir within 3.2 km radius (pd). In addition, and to account for extreme zero-inflation in cone counts, cone crop was included both as both a binary covariate (cc I = 1 if cone count is 0 and cc I = 0 otherwise), and as continuous covariate of observed cone crop counts (cc). This allows the intercept of the occurrence model to vary depending on whether or not cones are present.
We modeled year and season covariate effects as covariate-specific (k) random effects, i.e., with a covariate-specific mean (µ β,k ) and standard deviation (σ β,k ): for the k = 5 parameters of interest, i.e., the intercept, and the four covariate effects (cc, cc I , pw,and pd).
The resulting global model is: , ω = ω 1 , . . . , ω 5 , which includes or removes the k th term depending on whether ω k = 1 or ω k = 0, respectively. Thus, the general model can be reformulated to allow for cross-product model selection as follows: Note that the intercept term, β 1 , is always included in the model.
There are two sets of priors to be considered: the priors for the parameters to be estimated (Table 2), and the season-specific prior model probabilities (Table 3).
Although prior model probabilities were set to be equal, i.e., each model had an equal probability of being selected, there were certain covariate combinations that, based on the biology of Clark's nutcrackers, were considered not plausible a priori. Thus, prior model probabilities are equal across all plausible model

Model Fitting
To fit these models, we use JAGS (Plummer 2003), called from within R (R Core Development Team) using the library jagsUI (Kellner 2016). The data are available as an R workspace available as the supplement (nutcracker.RData).
Contained in this workspace are two data objects: • 100m radius data nutData.100 • infinite radius data nutData.all Within each object, all the data required to each model is contained: • the detection/non-detection data, y • the covariates (see Table 1) • the dimensions variables for the number of rows in the data, R

Goodness of fit
Below is the code for computing the posterior distribution of the Pearson's chi-square statistic comparing the observed (y) and expected encounter history frequencies (see JAGS code above).