Occupancy modeling of Black-backed Woodpeckers on burned Sierra Nevada forests

The Black-backed Woodpecker (Picoides arcticus) has been designated by the USDA Forest Service as a management indicator species for snags in burned conifer forests of the Sierra Nevada of California, USA. However, little is known about the characteristics that affect between-fire and within-fire habitat selection by the species in the region. Here we report on the first broad-scale survey of Black-backed Woodpeckers on wildfire-affected forests of the Sierra Nevada. We implemented a Bayesian hierarchical model to: 1) estimate Black-backed Woodpecker occupancy probability in fire areas burned within a time window of 1–10 years; 2) identify relationships between occupancy probability and habitat covariates (fire age, change in canopy cover pre-to-post fire, and snag basal area), elevation, and latitude; and 3) estimate detection probability and relate it to survey interval length and survey type (passive v. broadcast). We included random fire-area effects in our model of occupancy probability to accommodate clusters of nonindependent points surveyed within the larger set of fire areas. Mean occupancy probability was estimated to be 0.097. Elevation (after controlling for latitude) had the strongest effect on occupancy probability (higher occupancy at higher elevation) followed by latitude (higher occupancy at northerly sites). Fire age was also important; occupancy probability was about 43 higher on the youngest compared to oldest fires. Although the direction of regression coefficients were in the expected direction (positive), snag basal area and canopy cover change were of minor importance in affecting occupancy probability. There was some indication, however, that the importance of snag basal area increased with fire age. Weak links between occupancy and canopy cover change suggested the species uses a range of burn severities, and such heterogeneity may promote habitat longevity. Our estimate of overall detection probability (across all survey intervals) was 0.772. We found strong effects of survey interval length (higher for longer interval) and, especially survey type (higher for broadcast survey) on detection probability. Our modeling framework and implementation illustrates the flexibility of the Bayesian hierarchical approach for handling complexities such as estimating derived parameters (and variances) and modeling random effects, and should prove generally useful for occupancy studies.


INTRODUCTION
The use of indicator species to assess environmental conditions has a long history in ecology, conservation, and land management (Simberloff 1998, Lindenmayer et al. 2000, Carignan and Villard 2002, Niemi and McDonald 2004).Although the indicator species concept has been criticized, it remains popular for a variety of historical (e.g., the National Forest Management Act of 1976 in the United States) and practical (e.g., cost and ease of designing and implementing effective monitoring) reasons.In addition, scientific debate over the use of indicator species sharpened the focus of many modern applications; greater attention is now typically paid to stating clear objectives, selecting appropriate species, and testing assumptions related to how well the species actually reflects the processes or conditions it was chosen to indicate.
The indicator species concept is often employed to monitor effects of land management activities (Lindenmayer et al. 2000).Wildfire and salvage logging are priority management issues in coniferous forests of western North America (Hutto 2006), and the Black-backed Woodpecker (Picoides arcticus) was recently selected by the USDA Forest Service as a management indicator species (MIS) for snags in burned forests of the Sierra Nevada (USDA Forest Service 2007).Highseverity stand-replacing wildfires have been suggested to be particularly important for this species (Hutto 1995, 2008, Murphy and Lehnhausen 1998, Kotliar et al. 2002, Smucker et al. 2005, Hanson and North 2008), with habitat quality peaking within the first few years of a high-severity fire (Murphy and Lehnhausen 1998, Saab et al. 2007, Nappi and Drapeau 2009).However, the species may also respond positively to lower intensity fires such as controlled burns (Russell et al. 2009), and it also occurs to some extent in unburned forests (Dixon andSaab 2000, Tremblay et al. 2009).Habitat heterogeneity, with low-severity and unburned areas interspersed with more severely burned patches, may be important in maintaining the quality of woodpecker habitat beyond the first few years post fire (Nappi et al. 2010).Clearly, a better understanding of region-specific responses of Black-backed Woodpeckers to characteristics of burned forests is needed to refine their use as an indicator of conditions within recently burned forest stands.
Occupancy probability, or the proportion of sites occupied, is a common metric for wildlife monitoring studies (including MIS monitoring) and the assessment of wildlife-habitat relationships.It is now widely appreciated that estimates of occupancy will be biased low if detection probability is less than one, and that inferences about spatial and temporal variation in occupan-cy could be misleading unless detection probability is properly accounted for within the statistical modeling framework (MacKenzie et al. 2002, 2005, MacKenzie 2005, Bailey et al. 2007, Guillera-Arroita et al. 2010).Nevertheless, many occupancy studies still routinely ignore issues of detectability despite employing field methodologies that would permit its estimation.In addition, although often regarded as a nuisance parameter, detection probability itself may be of keen interest to land managers.This is especially true in the case of MIS or species of conservation concern, because detection probability estimation allows assessment of how much survey effort is needed to state with a given level of confidence the presence or absence of that species in the area of interest.Modeling detection probability can also allow evaluation of field methodologies and assist survey planning (Guillera-Arroita et al. 2010).
Here we develop and implement a hierarchical model to estimate detection and occupancy probabilities and identify habitat relationships of Black-backed Woodpeckers in burned forest stands of the Sierra Nevada of California, USA.Our model and implementation illustrate the flexibility of the approach for handling complexities such as estimating derived parameters (and variances) and modeling random effects.Random effects are particularly important to consider for surveys that follow a typical survey design whereby a series of (likely non-independent) sites (often survey points) is surveyed within a larger set of discrete study areas (i.e., cluster sampling; Thompson 1992).Our specific objectives were to: 1) estimate Black-backed Woodpecker point occupancy probability in Sierra Nevada forest stands that burned within a time window of 1-10 years; 2) identify relationships between occupancy probability and habitat covariates (time since fire, change in canopy cover pre-to-post fire, and snag basal area), elevation, and latitude; and 3) estimate detection probability and relate it to survey interval length and survey technique (passive v. call-broadcast survey) in an effort to evaluate the efficacy of our survey methodology and inform future surveys.With respect to the habitat modeling of occupancy (objective 2), we expected the species to occur relatively more frequently in more recent fire areas, at sites with greater loss of canopy cover, and at sites with higher snag basal area.

Sampling design
We randomly selected 51 forest stands from a pool of 72 fire areas containing !50 ha of forest that burned at mid-or high-severity (severity classes indicating moderate to complete mortality of the tree layer) in conifer forest habitats on 10 national forests in the greater Sierra Nevada region, California, USA, that burned between 1999 and 2008 (based on fire boundaries delineated in the geodatabase VegBurnSeveri-ty08_1.mdb; downloaded from hhttp://www.fs.fed.us/r5/rsl/clearinghouse/gis-downloadi).Here we define the greater Sierra Nevada region broadly to include not just the physiographic Sierra Nevada, but also spatially contiguous areas of the southern Cascades and Modoc Plateau within northeastern California, an area that is often thought of as an ecologically cohesive unit for management and conservation purposes (e.g., USDA Forest Service 2007).Selected fire areas ranged from 96 ha (2001 White Fire on Stanislaus National Forest) to 56,683 ha (2002 McNally Fire on Sequoia National Forest) in size.The total burned area (any severity) within the 72-fire sampling frame was 323,358 ha.
We conducted Black-backed Woodpecker surveys at 5-21 survey points in each fire area during a single visit between 16 May and 7 July 2009, resulting in a total of 899 points sampled.The number of points surveyed within a fire area was determined based on fire size, accessibility of the fire area, and the number of surveys that could be completed by 2 observers during the morning hours 0530-0930 PST.Survey points were established along transects that followed trails or secondary roads whenever possible, with the first survey point established as near as possible to a randomly selected point within the fire area.Transects along roads were offset from the road 50 m in a random direction [unless only one side of the road was accessible (e.g., due to a cliff ) or only one side of a road was burned].In cases where there were no roads or trails within the fire area, observers initiated the first survey point as near as possible to a randomly selected point and then continued in a random direction whenever possible.Survey points within a fire area were separated by approximately 250 m.

Black-backed Woodpecker surveys
A common survey methodology to allow simultaneous modeling of occupancy and detection probabilities is based on a repeated measures design.Under this design, replication can be achieved through repeated surveys (Mac-Kenzie and Royle 2005) or by dividing single surveys into sequential time intervals (Alldredge et al. 2007).The latter approach is a cost-effective, efficient alternative whenever all surveys can be completed within a time frame over which occupancy probabilities can be reasonably considered to be stationary (MacKenzie and Royle 2005).A disadvantage of the sequential interval sampling, however, is that the assumption of the most general occupancy model of independent and equal detection probabilities among sampling intervals is almost certainly violated.In particular, detection in one interval is likely to yield higher probability of detection in subsequent intervals (i.e., an observer-induced ''traphappy'' response; Riddle et al. 2010).This is especially true when calls are broadcast to enhance detection probabilities (Sliwa and Sherry 1992, Conway and Gibbs 2005, Wintle et al. 2005, Nadeau et al. 2008).In such circumstances, most information for estimating detection probability is contained within the first detection, such that a survey at a given site can be discontinued once the animal is detected without sacrificing precision of detection or occupancy probability estimates.Such surveys have been termed 'removal surveys', because a site is removed from further sampling once a detection is recorded (Farnsworth et al. 2002, MacKenzie andRoyle 2005).
Here we employed what might be termed a double-removal methodology based on singlevisit interval-survey data.At each survey point, we completed 1-2 passive (i.e., no broadcast used) survey intervals followed by three callbroadcast survey intervals.Passive surveys were initiated 1-2 minutes after the observer arrived at the survey point to allow local birds to settle.For both the passive and broadcast surveys we recorded a zero for intervals without detections and, for points with detections, a one for the interval of first detection.The number of passive surveys (one or two) was alternated between successive survey points.For survey points with one passive survey interval, interval duration was 3-min.For survey points with two passive intervals, a 3-min interval was followed by a 2min passive interval.The two-interval (5-min total) passive surveys were part of a broader multi-species study for which we recorded all bird species detected during those intervals; the methodology we employed is recommended and commonly used in such surveys (Ralph et al. 1995, Siegel et al. 2007).Passive surveys were followed by three 2-min broadcast survey intervals intended to enhance detection probability by eliciting responses from Black-backed Woodpeckers.Each broadcast survey interval consisted of playing Black-backed Woodpecker vocalizations (scream-rattle-snarl and pik calls) and territorial drumming sounds (recording by G. A. Keller, Macaulay Library of Natural Sounds, Cornell Laboratory of Ornithology) for approximately 25 s followed by 95 s of listening and watching for Black-backed Woodpeckers.Recordings were broadcast at a standardized volume using FOXPROt ZR2 digital game callers (FOXPRO Inc., Lewistown, PA, USA).In addition to recording intervals of detection for passive and broadcast surveys, we recorded numbers of individuals, ages and sexes of individuals, initial and minimum detection distances, and whether or not birds drummed or vocalized during the survey.Although we did not limit the point survey area, .90% of individuals were detected ,120 m from survey points.

Habitat data and exploratory analysis
Prior to formal modeling (see below), we examined the distribution of survey points with detections and those without detections (nondetection stations) in relation to covariates that we expected to influence Black-backed Woodpecker occupancy rates, using bean plots generated with the beanplot package (Kampstra 2008) in R (R Development Core Team 2009).Bean plots facilitate comparison of distributions of data points by simultaneously displaying the data along with normal density traces of the data.Our preconceptions of important predictors of woodpecker occupancy were based on our own field observations (including pilot field data collected in 2008) as well as previous research conducted on the species (Hanson andNorth 2008, Hutto 2008).We expected woodpecker occupancy to be influenced by fire age (higher occupancy in more recently burned areas), fire severity (higher occupancy in more severely burned areas), and pre-and post-fire forest composition and structure.Because forest composition is closely allied with elevation, we also examined the relationship between detections and elevation.Elevational distributions of birds in the Sierra Nevada vary from the northern to southern extents of the study region with individual species occurring at higher elevations in the south (Siegel et al. 2011); thus we also examined the distributions of detections and non-detections in relation to latitude and to the residuals of a linear regression of elevation on latitude.The latter essentially represents the response of birds to elevation after controlling for latitude, and it is what we included as a covariate in the logit-linear model for occupancy probability (see The model below).
Specific variables we examined included: (1) California Wildlife Habitat Relationships (CWHR) classification (California Department of Fish and Game 2005), as determined on the ground (2) elevation (m) determined from hand-held Geographic Positioning System (GPS) units and United States Geological Survey topographic maps; (3) snag (standing dead tree) basal area (m 2 /ha) estimated at each survey point by applying the Bitterlich variable-radius plot method (Mueller-Dombois and Ellenberg 1974) to stem counts made using a 10 basal area factor (BAF) timber-cruising crutch; (4) pre-fire % tree cover calculated using the Spatial Analyst Toolbox in ArcGIS (Ver.9.2, Environmental Systems Research Institute, Redlands, CA) by averaging midpoints of the mean % tree cover variable (WHRDEN-SITY) from 10,000-m 2 (100 3 100 m pixel) California Multi-source Land Cover Data hhttp:// frap.cdf.ca.gov/data/frapgisdata/download.asp?spatialdist¼1&rec¼fveg02_2i at 90,000-m 2 (300 3 300 m pixel) resolution; (5) number of years since fire (range ¼ 1 to 10 years based on USDA Forest Service data; see Sampling design above); (6) change in percent canopy cover, a measure of burn severity, based on satellite derived relativized difference normalized burn ratio score RdNBR (Miller et al. 2009) summarized at 8,100-m 2 (90 m 3 90 m pixel) resolution by averaging 900-m 2 (30 m 3 30 m pixel) values from Geographic Information System (GIS) coverages provided by the USDA Forest Service (J.D. Miller) using the raster package hhttp://cran.r-project.org/web/packages/raster/vignettes/Raster.pdfi in R (R Development Core Team 2009); and ( 7) latitude (in decimal degrees) recorded from USGS topographic maps.

The model
We developed our occupancy model based on i ¼ 1, . . ., N ¼ 899 survey points, j ¼ 1, . . ., M ¼ 51 fire areas, and k ¼ 1, . . ., K ¼ 5 survey intervals.The designation of two location indices (i and j ) reflects a typical study design where a series of sites, in our case survey points, is surveyed within a set of discrete study units, here fire areas (Thompson 1992).Often, covariate data of interest are available at these two different scales.
We followed the basic modeling approach outlined in Royle and Dorazio (2008; Chapter 3) and references therein.We modeled detections, y(i, j, k), conditional on occupancy, z(i, j ), such that: yði; j; kÞjzði; jÞ ; Bern zði; jÞp k where y(i, j, k) ¼ 1 if at least one Black-backed Woodpecker was detected at survey point i in fire area j during sampling interval k and y(i, j, k) ¼ 0 otherwise; and z(i, j ) represents the true unobserved occupancy state of the survey point, such that z(i, j ) ¼ 1 if one or more woodpeckers were at the survey point and z(i, j ) ¼ 0 if no woodpeckers were present.The probability of detecting at least one individual at an occupied survey point i in fire area j and interval k (i.e., Pr(y(i, j, k) ¼ 1 j z(i, j ) ¼ 1)) is a Bernoulli trial with success (i.e., detection) probability p k .The model assumes no false detections (i.e., Pr(y(i, j, k This assumption is reasonable because the species is easily distinguishable from similar species by both appearance and vocalizations. As noted above (Methods: Black-backed Woodpecker surveys), the y(i, j, k) observations consisted of recording the passive interval and broadcast interval of first detection (the double removal design).Note that although the K intervals varied among survey points in our study, depending on the number of passive survey intervals and whether or not a woodpecker was detected prior to the last broadcast interval, we input data as encounter histories with a fixed K ¼ 5 intervals (see Supplement), and treated missed intervals (due to conducting only one passive interval and/ or detecting a woodpecker prior to the last broadcast interval) as missing data (NAs).This data structure resulted in 16 possible detection histories: 00000, 00001, 0001NA, 001NANA, 01000, 01001, 0101NA, 011NANA, 0NA000, 0NA001, 0NA01NA, 0NA1NANA, 1NA000, 1NA001, 1NA01NA, 1NA1NANA.
We modeled the latent occupancy state, z(i, j ), as a Bernoulli random variable: such that is the probability of survey point i at fire area j being occupied by at least one Blackbacked Woodpecker.
We defined logit-linear models for both p k and w ij .We modeled detection probability, p k , as: where the regression parameters are represented by a's, such that a 0 is the intercept, or logittransformed detection probability for a 2-min passive survey interval; a 1 is the regression coefficient for the variable effort k , which is an indicator for interval length where effort k ¼ 1 for 3-min intervals ( just the first passive interval; i.e., k ¼ 1), and zero otherwise; and a 2 is the coefficient for the indicator variable type k , which denotes whether the count interval was passive (type k ¼ 0) or broadcast (type k ¼ 1).In a separate analysis (not reported here) we considered a model with a covariate representing basal area of live trees; however, there was little indication that this variable affected detectability.Here we consider detection probability to be homogeneous across space.
We defined the model for occupancy probability as: where b 0 is the intercept, or the overall mean logit-transformed occupancy probability for the average survey point at mean covariate values; fire j represents a random fire-area effect; and the remaining b.s represent the logit-linear coefficients for the model covariates.By including the random fire-area effect, we admit the likely similarity in occupancy probability (i.e., nonindependence) among points within a fire area.
Covariates in the model for w ij were defined as follows (see Methods: Habitat data and exploratory analysis for detail): fire.agej was the number of years since fire for fire area j (range ¼ 1 to 10 years); latitude i was the latitude in decimal degrees for survey point i; snag.bai was the basal area of snags at survey point i ; cc i was the change in percent canopy cover at survey point i; and elev.resi was the residuals from a regression of elevation on latitude for survey point i.By including the elevation residuals rather than elevation itself, we avoided problematic collinearity between these two variables (r ¼À0.681;N ¼ 899; P , 0.0001) in the logit-linear model (Graham 2003), and acknowledge the known dependence of bird responses to elevation on latitude (Siegel et al. 2011).
Although the intercept from our logit-linear model (Eq.2) allows inference about occupancy probability at a point within an average fire at mean values of covariates, we were also interested in estimating the mean proportion of survey points in the 2009 sample that were occupied by Black-backed Woodpeckers.This latter quantity is a derived parameter, involving both a known component, the number of sites with woodpeckers detected [i.e., y(i, j, k) ¼ 1] and an unknown component based on Eq. 1.It can be calculated as: Because we were interested in the relative efficacy of our two survey types, we also derived quantities representing detection probabilities for the combined passive intervals, the combined broadcast-survey intervals, and the overall detection probability (i.e., all intervals combined).We estimated detection probabilities for the combined passive intervals, p P , as: where 1 À p k represents the probability of not detecting an occurrence at a point during time k.
Thus, the probability detecting an occurrence at a point during the two passive intervals is one minus the probability of not detecting an occurrence during the two intervals.Similarly, we estimated detection probability for the broadcast intervals, p B as: and an overall (passive þ broadcast) detection probability, p., as: We estimated parameters using Bayesian Markov chain Monte Carlo (MCMC) methods (Gilks et al. 1996).We provided non-informative prior distributions for all model parameters.For the inverse-logit transformed intercepts from the logit-linear models for occupancy and detection probabilities [i.e., inverse-logit(a 0 ) and inverselogit(b 0 )] we used U(0, 1) priors; for the detection probability coefficients a 1 and a 2 we used U(À10, 10) priors; for the occupancy probability covariates b 1 -b 5 we used Norm(0, 0.001) priors; for the fire j random effect in the occupancy linear model we used a Norm(0, s fire ) prior; and for s fire , the precision parameter of fire j , we used a Gamma(0.001,0.001) prior.We standardized covariates of the logit-linear model for occupancy probability to facilitate parameter estimation and interpretation.Models were run using Win-BUGS (Spiegelhalter et al. 2003) from R (R Development Core Team 2009) via the R2Win-BUGS package (Sturtz et al. 2005).Posterior distributions were based on 30,000 iterations of 2 chains after discarding the first 20,000 simulations and thinning by 2. We assessed convergence by examining trace and density plots of posterior distributions of the two chains and from values of the potential scale reduction factor, R, which were ,1.1 for all model parameters (Gelman et al. 2003).We provide data and R and WinBUGS code in the Supplement.

Black-backed Woodpecker detections and exploratory analysis
We detected Black-backed Woodpeckers at 169 of 899 survey points distributed across 28 of the 51 fire areas.We detected Black-backed Wood-peckers on both the west and east sides of the Sierra crest, and across nearly the full latitudinal range of our study area, including the most northerly fire area we surveyed and the third most southerly fire area surveyed (Fig. 1).
Black-backed Woodpecker detections varied somewhat by major habitat type (Fig. 2A).Five points in mixed deciduous-coniferous habitats (montane riparian and mixed coniferous forest) had no detections and are not shown in Fig. 2. Pinyon-juniper forest (PJN) was notable among coniferous habitats for having 49 points surveyed, but no Black-backed Woodpecker detections.White fir (WFR) and Sierra mixed conifer (SMC) habitats had relatively few detections with 6% and 12% of the survey points in those habitats with detections, respectively.The remaining coniferous forest habitats sampled had 21-33% of survey points with detections.Much of the habitat-related variation in detections could be explained by the residuals of a regression of elevation on latitude (Fig. 2B).In nearly all forest types, detections tended to occur at survey points that were at relatively high elevations, and detections tended to occur at about the same relative elevation across forest types (elevation residual values of ;þ100-200 m above the mean of 0).
We also found a latitudinal gradient in Blackbacked Woodpecker detections, with detections more common at higher latitude fire areas (Fig. 3).There seemed to be some evidence of a relationship between detections and each of the remaining predictor variables considered except pre-fire canopy cover (Fig. 3).Because of the lack of relationship between pre-fire canopy cover and detections, we did not consider this variable in our linear model for occupancy probability (below).

Occupancy modeling
Occupancy probability.-Weestimated an overall mean occupancy probability (i.e., mean occupancy probability for the average fire area at mean covariate values) of 0.097, with a central 95% credible interval (hereafter CRI) of 0.049-0.162.Our estimate of the proportion of points surveyed during the 2009 that were occupied by woodpeckers was much higher: mean ¼ 0.252; 95% CRI: 0.219-0.299,suggesting that most occurrences were clustered within a few sites and/or extreme covariate values.
Three of five covariates included in the logitlinear model for w ij had 95% CRIs for regression coefficients that did not include zero, and thus appeared to be important predictors of occupancy probability: fire.agej , latitude i , elev.res i (Table 1).Standardized regression coefficients and plots of predicted occupancy probability across observed covariate values (Table 1, Fig. 4) suggest that elev.resi (i.e., elevation adjusted for latitude) had the strongest effect on occupancy probability, followed by latitude i and fire.agej .Mean predicted occupancy probability was higher for survey points at higher elevations (for a given latitude), at more northerly latitudes, in more recent fire areas.The smaller standardized regression coefficient and 95% CRI that included zero suggested that snag basal area, snag.bai , had a relatively minor effect on occupancy probability (Table 1); however, the predicted change in occupancy probability over the range of observed snag.bai values suggested a potential role for this variable (Fig. 4).High snag basal area may be especially important in maintaining viability of woodpecker habitat in older fire areas.For example, snag basal area was relatively high for survey points with detections compared to non-detection survey points in older fire areas (6-10 years old); but was similar for detection and non-detection survey points in younger (1-5 year-old) fire areas (Fig. 5).The standardized regression coefficient for cc i was similar in magnitude to the coefficient for snag.bai and the 95% CRI for this parameter also included zero; the direction of this effect was, however, as expected (positive), suggesting that Black-backed Woodpecker occupancy rates were slightly higher at survey points that burned at higher severity (Table 1).
Detection probability.-Botheffort k (interval duration) and, especially, type k (passive v. broadcast interval) were important in affecting detection probability, p k (Table 1).Our estimated mean detection probability for the 2-minute passive count interval, p 2 , was just 0.065 (95% CRI: 0.024-0.126),and for the 3-minute passive interval, p 1 , 0.178 (95% CRI: 0.125-0.239).Our estimate of the mean overall probability of detection for the two passive intervals combined, p P , was 0.231 (95% CRI: 0.164-0.309).We estimated the mean detection probability for a call broadcast interval, p 3-5 to be: 0.337 (95% CRI: 0.253-0.424),and the overall detection probability for the three call-broadcast intervals combined, p B , of 0.705 (95% CRI: 0.584-0.809).Most of the birds detected during passive count intervals (33 of 45 total passive interval detections) were also detected during the broadcast intervals, and there was little difference between the overall broadcast detectability of 0.705 and the overall detectability, p., based on the combined passive and broadcast surveys: mean ¼ 0.772 (95% CRI: 0.663-0.860).Extending these estimates to include multiple site visits, at least 9 visits to an occupied site would be required to achieve a 90% chance of detecting a Black-backed Woodpecker if we were only using the twointerval 5-min passive count technique, compared with only two visits using either the combined passive-broadcast or broadcast-only methods (Appendix).

DISCUSSION
Management indicator species (MIS) are commonly monitored by land management agencies to assess effects of management activities (Lindenmayer et al. 2000).The Black-backed Wood-pecker (Picoides arcticus) has been designated by the USDA Forest Service as a MIS for snags in burned forests of the Sierra Nevada (USDA Forest Service 2007).However, little information exists on this species' response to characteristics of burned areas within the Sierra Nevada region v www.esajournals.org(Hanson and North 2008).Our study represents the first broad-scale survey of Black-backed Woodpeckers across recently burned forests of the Sierra Nevada and the first to assess responses of the species to habitat characteristics relevant to managing habitat for the species in the region and refining its use as an MIS.
An important aspect of our study was the monitoring application of providing an estimate of mean occupancy probability for Black-backed Woodpeckers on 1-10 year-old fire areas on national forests during the 2009 breeding season.
We estimated the mean occupancy probability for an average fire area at mean values of covariates to be 0.097.Our estimate of the mean proportion of survey points occupied by woodpeckers, however, was 0.252, suggesting considerable overdispersion in the data with most occupied points clustered within relatively few fire areas.This discrepancy highlights the importance of including cluster-level effects in models (here our random effect fire j ) to account for such overdispersion.If we assume our sample to be representative of 1-10 year-old fire areas in the   v www.esajournals.orgto suggest a rough estimate of about 81,486 ha (or a range based on the 95% CRI of 70,815-96,684 ha) of area occupied by Black-backed Woodpeckers.Continued monitoring will enable assessment of changes in habitat availability, occupancy probability, area occupied, and the distribution of woodpeckers.
An additional goal of our study was to assess habitat effects on occupancy probability.Our exploratory analysis of the distribution of detections among pre-fire habitat types and pre-fire canopy cover suggested that Black-backed Woodpeckers did not distinguish among fire areas based on pre-fire forest composition.This finding is in agreement with other broad-scale studies of the Black-backed Woodpeckers (e.g., Hutto 2008), and with findings of smaller scale studies from various parts of the species' range that report foraging and nesting on fire-killed snags of variety of coniferous tree species (Bull et al. 1986, Villard and Beninger 1993, Villard 1994, Murphy and Lehnhausen 1998).One exception to the apparent lack of pre-fire habitat-type selectivity in our study was pinyon-juniper forest The mean difference in snag basal area between stations with detections and those without detections was greater for older fire areas (mean of 16.17 m 2 /ha for sites with detections vs. 7.75 m 2 /ha for sites without detections) than for newer fire areas (12.41 vs. 9.98 m 2 /ha for detections and non-detections, respectively).v www.esajournals.org(PJN in Fig. 2) for which we had no detections; this habitat may not have enough large-boled trees to support this species (Nappi et al. 2003).
Of the three habitat covariates considered in our occupancy model related to post-fire habitat, the strongest response of Black-backed Woodpecker occupancy probability was to fire age.Our model predicted woodpecker occupancy probability dropping off markedly in the first few years following a fire.This finding is consistent with other studies that report the strongest positive responses of Black-backed Woodpecker populations to fire disturbance 1-3 years post-fire (Murphy and Lehnhausen 1998, Saab et al. 2007, Nappi and Drapeau 2009), during which time birds are likely responding to increased abundance of wood-boring beetle larvae (Villard andBeninger 1993, Murphy andLehnhausen 1998).Yet, the relationship between fire age and occupancy may be more complicated than suggested by our model.For example, our sample showed detections peaking 2-3 years post-fire, a smaller peak 6 years post-fire, and at least a few occupied survey points (on two of seven fire areas) out to 10 years post-fire, the temporal extent of our sampling frame.Because fire ages in our sample were from a single sampling season, it is impossible to say whether this pattern is pervasive or simply a reflection of characteristics of fires that burned in those particular years.Others have, however, also shown the species to remain fairly common 6-8 years post-fire (Hoyt andHannon 2002, Nappi et al. 2010), and within-fire heterogeneity in fire severity may create conditions that promote a steady supply of dying trees, diverse prey resources, and habitat longevity extending beyond the first few years post-fire (Nappi et al. 2010).
Contrary to our expectations, we did not find a strong positive response of occupancy probability to the covariate that described fire severity (cc i ).The direction of the response, was however, in the expected direction (higher occupancy where canopy cover loss was greatest).The lack of strong response of woodpeckers to canopy cover change could reflect the fact that our sampling frame was limited to points within fire areas that contained at least 50 ha of forest that burned at mid-or high-severity, and woodpecker survey points were typically not very far from a severely burned patch.Thus, woodpecker responses to burn severity at a larger landscape scale may have obscured patterns summarized at individual sampling points (in our case at the scale of 90 3 90 m pixels).Studies that have shown strong positive responses of Black-backed Woodpeckers to high-severity stand-replacing wildfires have sampled a wider range of conditions, comparing sites between burned and unburned forest stands or individual sites prior to and following a fire (Hutto 1995, 2008, Murphy and Lehnhausen 1998, Kotliar et al. 2002, Smucker et al. 2005, Hanson and North 2008).Nevertheless, our finding of a relatively weak fire-severity effect could indicate that Black-backed Woodpeckers in the Sierra Nevada may require, or at least be more tolerant of, more heterogeneous landscapes than has been suggested for other regions (e.g., Hutto 2008).
The size of the effect of snag basal area (snag.bai ) on occupancy probability, as measured by the standardized regression coefficient, was similar to the effect of canopy cover.However, because of the skewed distribution of this variable, the potential effect of snag basal area over the range observed values (although extreme values were rare), suggested its potentially greater importance in affecting occupancy patterns.It should also be noted, however, that the canopy cover change and snag basal area were positively correlated (r ¼ 0.254) and so including both as covariates in the occupancy model may have diluted the apparent effects of each to some extent.Comparison of the distributions of survey points with and without detections for young fire areas (1-5 years post-fire) to distributions of detection and non-detection points for older (6-10 year-old) fire areas, suggested that the importance of snag basal area for maintaining local woodpecker populations increases with fire age.We plan to explore this possibility further in future modeling efforts (e.g., by including a fire.agei 3 snag.bai interaction term in the logitlinear model for occupancy probability).If this inference holds, it suggests that salvage logging following wildfire may not just negatively affect short-term responses of woodpeckers and other wildlife (Koivula andSchmiegelow 2007, Hanson andNorth 2008), but it may have particularly strong negative effects on the longevity of fire areas as habitat for these species.
A final objective of our study was to estimate and model detection probability.Our logit-linear model for detection probability allowed us to assess effects of survey interval type (type k ) and duration (effort k ), both of which (especially survey interval type) positively affected detection probability.Furthermore, by summarizing MCMC output for model parameters we were able to easily provide point estimates and credible intervals for detection probability at the level of the survey interval (either 2-or 3-minute interval), survey interval type (passive or broadcast), and the complete survey (all five intervals combined).Our survey objectives did not include meeting any particular detection probability threshold.However, in some instances land managers might need to determine with a known level of certainty whether Black-backed Woodpeckers are present in a project area.If we extend our estimates of detection probability for survey interval types to include multiple site visits, we would need to visit an occupied point at least 9 times to achieve a 90% chance of detecting a Black-backed Woodpecker if we were only using the two-interval 5-min passive count technique.On the other hand, if we used either the combined passive-call-broadcast or broadcastonly approaches, we would only need to visit an occupied point twice to achieve this detectability threshold.This suggests that we could eliminate the passive survey altogether with little effect on our ability to estimate detection probability or occupancy of Black-backed Woodpeckers.
Our analysis of detectability suggests that broadcast surveys may be important for reliably confirming the presence of Black-backed Woodpeckers at a survey point with any reliability.However, broadcast surveys have certain drawbacks.In some instances woodpeckers may respond to broadcasts (even if only by silently approaching closer) from perhaps hundreds of meters from a survey point.Furthermore, the same individuals may respond to broadcasts at successive points.This complicates interpretation of occupancy somewhat, which should in this case be interpreted as 'use'.Although we account for lack of independence of survey points with the random fire area effect in our occupancy model, we are currently initiating a radiotelemetry to provide better understanding of habitat use and scales of responses of woodpeckers to playback.Such data will provide a more thorough evaluation of our survey methodology and more meaningful biological interpretation of our occupancy parameter.
One consequence of using call broadcasts to enhance detectability is that once a bird responds to the broadcast and is detected by the observer, detection probability in subsequent intervals is enhanced, suggesting that for complete-detection-history data a 'behavioral response' model is needed (analogous to the M b closed-population capture-recapture model of Otis et al. 1978).Under this model, detections beyond the first detection interval yield no new information about detectability, and suggest that a removal design (Farnsworth et al. 2002, MacKenzie andRoyle 2005), such as we implemented here, may be preferred.From a sampling standpoint, the removal design is more efficient (at least in the case of a single species survey), because the survey is discontinued at a site once an individual is detected.This can allow time to survey more sites or add more visits/intervals at sites without prior detections, which should improve the ability to model and estimate occupancy probabilities.Although data from removal surveys yield less overall information about detectability and afford fewer options for modeling this parameter than complete detection histories (e.g., time-varying models; Alldredge et al. 2007), removal designs may in some circumstances require fewer overall surveys (where survey ¼ number of sites 3 number of visits or intervals) than complete detection history designs to achieve a given level of precision in occupancy estimates (MacKenzie and Royle 2005).
Our study, the first spatially extensive survey of Black-Backed Woodpecker across the greater Sierra Nevada region, provides managementrelevant information on the spatial distribution and habitat relationships of this designated MIS for snags in burned forest.As has been reported elsewhere, the species is closely allied with areas recently affected by wildfire; however, heterogeneity of burn severities within fire areas may promote habitat longevity.Future modeling efforts will use additional years of monitoring data to understand the implications of dependence on an ephemeral habitat-recent fire areas-on occupancy dynamics.Our analysis of detection probability also has important implications for designing surveys for this species; in particular, use of call broadcasts was clearly an important tool for allowing effective modeling of detection and occupancy probabilities with realistic sample sizes.Beyond the scope of our study, our occupancy modeling framework and implementation illustrate the flexibility of the Bayesian hierarchical approach for handling complexities such as estimating derived parameters (and variances) and modeling random effects.Given the surge of interest in occupancy modeling with imperfect detection (MacKenzie et al. 2005, Bailey et al. 2007), and the often complicated survey designs and objectives of these studies, the Bayesian hierarchical modeling approach should prove broadly useful.

Fig. 1 .
Fig. 1.Distribution of 51 fires surveyed for Black-backed Woodpeckers on National Forests of the Sierra Nevada during 2009 and proportions of points with detections.The range map of this species, based on California Department of Fish and Game (2005), is shown to highlight the spatial extent of sampling relative to the species' distribution in California.Note that sampling and detections occurred outside of the described yearround range.

Fig. 2 .
Fig. 2. (A) Distribution of detections and non-detections by habitat type.Habitat types follow the California Wildlife Habitat Relationships (CWHR) classification system (California Department of Fish and Game 2005) and include: SMC ¼ Sierra mixed conifer, JPN ¼ Jeffrey pine, EPN ¼ eastside pine, PPN ¼ ponderosa pine, RFR ¼ red fir, PJN ¼ pinyon-juniper, WFR ¼ white fir, JUN ¼ juniper, LPN ¼ lodgepole pine.Data from five points in mixed hardwood-conifer and montane riparian habitat for which no detections were recorded are not shown.Asymmetric bean plots in (B) show distributions of survey points by habitat type for which detections were recorded (left of vertical center lines) and for which no detections were recorded (right of vertical center lines) in relation to residuals of a regression of elevation on latitude.Data points are shown as tick marks, with longer tick marks representing multiple points at the same elevation.Shaded regions delineate density traces of the data.Black horizontal lines show mean elevations for detection (left of center lines) and non-detection (right of center lines) point survey stations.The dashed line at zero shows the mean residual value for all habitat types.Habitat related variation in detections appeared to be due largely to differences in the distributions of elevations among habitats.

Fig. 3 .
Fig. 3. Bean plots showing distributions of detections and non-detections in relation to geographic and habitat variables.Shaded (gray) regions show density traces of the data; the individual data points are represented by white (inside trace densities) or gray (outside trace densities) line segments.Means are indicated by bold black lines.Survey points with Black-backed Woodpecker detections were characterized by greater fire-induced change in canopy cover, greater snag basal area, and fewer years since fire than survey stations where Black-backed Woodpecker was not detected.There was no apparent difference between detection and non-detection points in pre-fire percent canopy cover.

Fig. 4 .
Fig. 4. Predicted mean occupancy probability and 95% credible intervals in relation to individual covariates included in the hierarchical occupancy model.See Methods: The Model for variable definitions.

Fig. 5 .
Fig. 5. Bean plots showing the distributions of detections (left sides of beans in dark gray) versus nondetections (right sides in light gray) for recent (1-5 year-old) fire areas and older (6-10 year-old) fire areas.Shaded (gray) regions show density traces of the data; the individual data points are represented by white (inside trace densities) or gray (outside trace densities) line segments.Means are indicated by bold black lines.The mean difference in snag basal area between stations with detections and those without detections was greater for older fire areas (mean of 16.17 m 2 /ha for sites with detections vs. 7.75 m 2 /ha for sites without detections) than for newer fire areas (12.41 vs. 9.98 m 2 /ha for detections and non-detections, respectively).

Table 1 .
Posterior summaries (means, standard deviations [SD], and 2.5% and 97.5% credible interval boundaries [lower and upper 95%]) for intercepts and regression coefficients (continuous variables for coefficients b 1-5 standardized to facilitate comparison) from the models for occupancy probability (w ij ) and for detection probability ( p k ) applied to the Black-backed Woodpecker survey data from 899 points on 51 burned forest stands.