Forest birds respond to the spatial pattern of exurban development in the Mid-Atlantic region, USA

Study area

Our study area encompassed approximately 4,300 km² and included nine counties in north-central Virginia (Clarke, Culpeper, Fauquier, Frederick, Madison, Page, Rappahannock, Shenandoah, and Warren) and two in western Maryland (Washington and most of Frederick; Fig. 1), USA. The region has experienced high population growth rates, ranging from 4% (Page County) to 36% (Culpeper County) in the past decade (US Census Bureau, 2013). The region has also experienced an increase in exurban settlements over the same time period (Suarez-Rubio, Lookingbill & Elmore, 2012), stimulated at least in part by the close proximity of natural amenities (Suarez-Rubio, Lookingbill & Wainger, 2012).

Breeding bird survey

Using the North America Breeding Bird Survey (BBS) (Peterjohn & Sauer, 1994; Sauer, Fallon & Johnson, 2003), a large-scale annual roadside survey to monitor the status and trend of breeding bird populations in the USA and southern Canada, we selected two groups of species that represent contrasting habitat preferences (forest vs. edge). Forest species—Ovenbird (Seiurus aurocapilla), Red-eyed Vireo (Vireo olivaceus), American Redstart (Setophaga ruticilla), Wood Thrush (Hylocichla mustelina), Scarlet Tanager (Piranga olivacea), and Eastern Wood-Pewee (Contopus virens) (Poole, 2005)—were defined as birds that use a wide variety of deciduous and mixed deciduous-coniferous forests and that might favor interior forested habitats (Mikusiñski, Gromadzki & Chylarecki, 2001). Forest-edge species—Eastern Towhee (Pipilo erythrophthalmus), Eastern Phoebe (Sayornis phoebe), Gray Catbird (Dumetella carolinensis), Northern Cardinal (Cardinalis cardinalis), and Indigo Bunting (Passerina cyanea) (Poole, 2005)—are those species strongly associated with forest edges and open habitats (Mikusiñski, Gromadzki & Chylarecki, 2001). These 11 species were also selected because they were detected on at least 5% of surveys during the 1986–2009 interval. In addition, many of the species are reported to have experienced population declines or reduced fecundity due to habitat loss or fragmentation (Donovan & Flather, 2002; HaganIII, 1993; Sherry & Holmes, 1997; US NABCI Committee, 2009).

BBS routes involve 39.4 km-long road transects, with 3-minute point count surveys conducted at stops every 0.8 km. From each BBS route located in the study area, we selected every fifth stop along the route to reduce overlap between adjacent areas around survey stops and decrease the likelihood of spatial autocorrelation (Moran’s I = 0.108, p = 0.182). We only considered survey stops that had detailed direction descriptions (i.e., geocoding information and characterization of site-specific features) and fell within the study region (125 survey points in total) (Fig. 1). We focused our analysis on survey stops instead of the entire route because of our interest in local variability of breeding habitats.

To characterize local characteristics of breeding habitats, we established potential zones of influence (Glennon & Kretser, 2013) of 400-m and 1-km radius around the selected BBS stops. These areas represented both breeding bird territories (Bowman, 2003; Mazerolle & Hobson, 2004), which were assumed to be in the immediate surroundings of survey stops, and areas feasibly visited during bird daily movements (Krementz & Powell, 2000; Lang et al., 2002). Within these areas, we quantified the proportion of forest and exurban development and the spatial pattern of exurban development from 1986 to 2009.

We used a hierarchical Bayesian model to adjust BBS counts (Suarez-Rubio et al., 2013) and account for BBS sources of variability such as observer differences (Sauer, Peterjohn & Link, 1994), first-year observers’ skills (Erskine, 1978; Kendall, Peterjohn & Sauer, 1996), environmental conditions (Robbins, Bystrak & Geissler, 1986), and habitat features (Sauer, Pendleton & Orsillo, 1995). We modeled count data as hierarchical over-dispersed Poisson and fit models using Markov Chain Monte Carlo (MCMC) methods in WinBUGS 1.4.3 (Lunn et al., 2000). We specified C_it as the count for each species on stop i and time t where i = 1, …, N; t = 1, …, T; and N and T were the number of stops and the number of years species were observed, respectively. C_it was assumed to be Poisson distributed with mean μ_it ${\displaystyle C_{it}\sim \text{Pois}\left({\mu }_{it}\right)}$ and the full model was: ${\displaystyle log\left({\mu }_{it}\right)={\beta }_{0\mathrm{stop}}+{\beta }_{1\mathrm{stop}}\times {\mathit{Year}}_t+{\beta }_2\times {\mathit{FirstYear}}_{it}+{\mathit{Route}}_{it}+{\mathit{Observer}}_{it}+{\mathit{Error}}_{it}}$ where each stop was assumed to have a separate intercept (β₀) and time trend (β₁). The model included several sources of variability including unknown route environmental conditions and habitat features (Route_it), observer effects (Observer_it), first-year observer effects (FirstYear_it) and over-dispersion effects (Error_it). Given that route conditions could also change among years, we also included year into the model. We used two Markov chains for each model and examined model convergence and performance through Gelman–Rubin diagnostics (Gelman, Carlin & Rubin, 2004; Link & Barker, 2010). Once the model reached convergence, we derived estimates of the count at each stop and in each year which were then used for the threshold analysis.

Defining exurban development

To characterize the land cover in the areas around survey stops, we classified Landsat 5 TM images (pixel size: 30 m) for 1986, 1993, 2000, and 2009. We performed standard pre-processing procedures (atmospheric and topographic correction) prior to image classification and conducted a supervised classification of areas of exurban development using a training dataset generated from aerial photos. Exurban development was defined as areas with housing densities between 1 unit per 0.4 ha and 1 unit per 16.3 ha (e.g., 6–250 houses per km²) (Brown et al., 2005). We identified exurban development using both spectral and structural characteristics following the methods outlined in Suarez-Rubio, Lookingbill & Elmore (2012). We derived spectral characteristics from spectral mixture analysis (Adams, Smith & Johnson, 1986) of the corrected Landsat images to estimate the fractional cover of vegetation, substrate, non-photosynthetic vegetation, and shade within each image. Based on spectral mixture analysis outputs, we built decision trees to classify exurban development for each of the four image dates.

To further analyze pixels belonging to branches of the decision trees that could not discriminate between exurban and urban areas based on spectral characteristics alone, we used morphological spatial pattern analysis (MSPA) (Soille, 2003; Vogt et al., 2007). The analysis evaluates map geometry by applying mathematical morphological operators to allocate each pixel to one of a mutually exclusive set of classes. We used an 8-neighbor rule as our structural element (i.e., both cardinal directions and diagonal neighbors are considered) and edge width of one. Pixels that fell into the MSPA-Islet (representative of isolated housing units), Bridge, Branch, and Loop classes (representative of associated roads) were considered exurban development. All other MSPA classes were considered urban development. Lastly, all cells originally designated as exurban development in the decision tree were then added back to attain the final exurban development maps. Overall classification accuracy for the final exurban development maps ranged from 93 to 98% (kappa: 0.87–0.96) (Suarez-Rubio, Lookingbill & Elmore, 2012).

Analyzing the spatial pattern of exurban development

To examine the spatial pattern of exurban development, we used the final exurban development maps as foreground and analyzed them using MSPA. Here, we focused specifically on the Islet class which represented scattered, isolated housing units. Using the MSPA classification output, we developed a compactness index to describe how clustered exurban development was in the area surrounding each survey stop at each time period considered. The compactness index was a measure of the proportion of exurban development within any MSPA classes other than the Islet class (i.e., 1−(Exurban Development islets/Exurban Development all classes)) and ranged from 0% (all Islets) to 100% (no Islets). Survey stops lacking exurban development within the potential zone of influence were excluded from the analysis (28 and 20 survey stops for the 400-m and 1-km radius buffers, respectively were excluded). Hence, dispersed exurban development was represented by 0% and maximally clumped exurban development by 100% compactness (see example in Fig. 2).

Identifying species response to compactness of exurban development

To examine the relationship between compactness of exurban development and bird species at the survey stops, we fitted a non-parametric locally weighted polynomial regression (loess) (Cleveland & Devlin, 1988). When the loess regression highlighted nonlinearity in the relationship, then a change-point analysis was used to test for a nonlinear threshold response.

We estimated potential species threshold responses to compactness of exurban development using Threshold Indicator Taxa ANalysis (TITAN) (Baker & King, 2010). TITAN allows the identification of change points in both occurrence frequency and relative abundance of individual species along an environmental gradient. It distinguishes responses of individual species with low occurrence frequencies or highly variable abundances and does not assume a linear response along all or part of an environmental gradient. TITAN uses normalized species scores (z) to establish a change-point location that separates the data into two groups and maximizes association of each species with one side of the partition. Z scores measure the association of a species’ abundance weighted by their occurrence and are normalized to facilitate cross-species comparison. Thus, TITAN distinguishes if a species responds to an environmental stressor (in this case compactness of exurban development) and whether the response is negative (z−) or positive (z+).

To measure quality of the response and assess uncertainty around change-point locations, TITAN bootstraps the original dataset and recalculates change points with each simulation. Uncertainty is expressed as quantiles of the change-point distribution. Narrow intervals between upper and lower change-point quantiles (i.e., 5 and 95%) indicate a nonlinear response in species abundance whereas broad quantile intervals are characteristic of species with a linear or more gradual response. Diagnostic indices of the quality of the response are purity and reliability. Purity is the proportion of bootstrap replicates that agree with the direction of the change-point for the observed response. Pure indicators (purity ≥ 0.95) are those that consistently assign the same response direction during the resampling procedure. Reliability is the proportion of change-point individual value scores (IndVal) among the bootstrap replicates that consistently have p-values below defined probability levels (0.05). Reliable indicators (reliability ≥ 0.95) are those with consistently large IndVal.

We ran TITAN (R package: TITAN2) (Baker & King, 2010) for the 11 selected bird species and compactness index in R 3.1.1 (R Development Core Team, 2013). We used the minimum number of observations on each side of the threshold split that is required by TITAN (n = 5) and specified 250 permutations to compute z scores and diagnostic indices as suggested by Baker & King (2010).

Evaluating species responses to forest loss and forest fragmentation in relation to compactness of exurban development

To evaluate the effects of compactness of exurban development in relation to other factors known to affect birds (i.e., forest loss and forest fragmentation), we used generalized additive models (GAMs) (Hastie & Tibshirani, 1990). GAMs were used to better account for potential non-linear trends between the response and predictor variables (e.g., Guisan, Edwards Jr & Hastie, 2002; Zuur et al., 2009). GAMs require fewer assumptions of data distributions and error structures, assuming only that functions are additive and components can be smoothed by local fitting to subsets of the data.

The models used adjusted counts for each bird species as dependent variables and compactness of exurban development, proportion of exurban development, proportion of forest, number of forest patches greater than 0.45 ha, and forest edge as predictor variables. The latter variables were estimated following Suarez-Rubio et al. (2013). Gaussian errors and an identity link were used, and smoothing parameters were automatically selected based on the effective degrees of freedom and a generalized cross validation criterion in R package mgcv (Wood, 2001; Wood, 2006). We did a multi-model comparison using a stepwise backwards selection process and calculated the Akaike information criterion (AIC_i) and the $\mathrm{\Delta }{\mathrm{AIC}}_i$ to rank models and select a best-fitted model (Zuur et al., 2009). We used the results to strengthen the inference regarding factors affecting birds in forested environments. Models were evaluated based on graphical diagnostic plots and the explanatory power of a model was assessed by examining the amount of the explained deviance. Predictors of the best-fitted model with high significance levels (p < 0.01) were identified as key factors that have strong effects on bird species.

Landscape composition and compactness of exurban development around survey stops

Landscape composition around survey stops changed through time during the time period studied, except for the 21% of stops that were inside protected areas (Table 1). The inclusion here of MSPA classes that represented associated roads (i.e., Bridge, Branch, and Loop) in addition to scattered isolated pixels (i.e., Islets) in the definition of exurban development differed from other operational definitions of exurban development used in previous work; as a result, the total amount of development that was classified as exurban was higher for our study than was reported for more restrictive definitions (e.g., Suarez-Rubio, Lookingbill & Elmore, 2012). For both the 400-m and 1-km radius buffers, there was a 6% increase in exurban development from 1986 to 2009 (Table 1).

Compactness of exurban development also increased over time (Table 1). For the 400-m radius buffer, compactness increased from 18% in 1986 to 39% in 2009. For the 1-km radius buffer, compactness increased even more, from 11% in 1986 to 44% in 2009. For both extents, the increase was higher between 2000 and 2009 than for any other time period. Compactness was slightly correlated with the amount of exurban development (Pearson’s correlation coefficient for 400-m buffer: 0.38, and 1-km buffer: 0.46) and not correlated with forest at either extent (Pearson’s correlation coefficient for 400-m buffer: −0.15, 1-km buffer: 0.04).

Response of bird species to compactness of exurban development

Non-parametric locally weighted polynomial regression (loess) models indicated a non-linear relationship between the compactness index and abundance of selected bird species (Fig. 3). Forest species differed in their threshold response to compactness of exurban development (Fig. 4). For the 400-m radius buffer, only one of the six forest species (i.e., Scarlet Tanager) showed a significant and reliable threshold response to compactness. Although Wood Thrush also responded negatively, the quality of the indicator was less reliable (0.80) (Table 2). In contrast, for the 1-km radius buffer, almost all forest species responded positively and reliably to the compactness of exurban development (Table 2).

Forest-edge species also had significant though less consistent threshold responses to compactness of exurban development at both extents (Fig. 4). For the 400-m radius buffer, Eastern Phoebe and Gray Catbird had a significant positive response to the compactness metric, while Eastern Towhee responded negatively to compactness. For the 1-km radius buffer, Eastern Phoebe, Gray Catbird, and Indigo Bunting responded positively to compactness, with reliability values and change points spanning a wide range of compactness values, similar to the finding for forest birds (e.g., Red-eyed Vireo, Eastern Wood-Pewee; Fig. 4).

In general, reliability information was redundant with purity (i.e., species with ≥0.95 purity were usually also reliable) (Table 2). In some instances, the direction of the response changed with extent of analysis. Wood Thrush responded positively to compactness of exurban development for the 1-km radius buffer. Although the direction of the response changed for the 400-m radius buffer, the indicator was not reliable at this extent (reliability = 0.80). For other species (e.g., Scarlet Tanager and Eastern Towhee), wide confidence bands and low z scores at the 400-m extent, highlighted uncertainty when the abundance distributions did not show a clear response. Therefore, where there were differences in the reliability and direction of response at different extents, the 1-km relationships were more reliable.

Most species (both forest and forest-edge) had relatively broad bootstrapped change-point distributions indicating that there were not sharp threshold responses to compactness of exurban development (Fig. 4). In addition, the width of the bootstrapped change-point distributions varied between the two buffer distances for only a few species. For example, Eastern Phoebe was one of the few species with a sharp response to compactness, but this occurred only at the 400-m radius buffer.

Response of bird species to forest loss and forest fragmentation in relation to compactness of exurban development

When forest loss and forest fragmentation were included as predictor variables in addition to the exurban development measures (i.e., proportion and compactness), forest had a highly significant effect on all forest species modeled and most forest-edge species at the 1-km radius buffer (Table 3). Number of forest patches had a significant influence on Red-eyed Vireo and Scarlet Tanager, and forest edge did not affect any of the forest species. The effect of exurban development varied among forest species. Only Red-eyed Vireo was significantly influenced by both proportion of exurban development and compactness of exurban development. Eastern Wood-Pewee and Wood Thrush were influenced by compactness of exurban development, whereas Scarlet Tanager was only influenced by proportion of exurban development.

None of the forest-edge species were influenced by compactness of exurban development at the 1-km radius buffer, although Eastern Phoebe, Eastern Towhee, Indigo Bunting, and Northern Cardinal were affected by its proportion. Regarding forest fragmentation, Indigo Bunting and Northern Cardinal were influenced by number of forest patches, whereas Eastern Phoebe, Eastern Towhee, and Gray Catbird were affected by forest edge. Models at the 400-m buffer and for American Redstart and Ovenbird at the 1-km buffer did not converge.