Conservation of bird populations is increasingly focused on landscapes. We combined data collected in 2005–2011 from 16,250 North American Breeding Bird Survey (BBS) survey points with local and remotely sensed environmental data to model the distribution of 7 grassland bird species in the Northern Great Plains of the United States. We analyzed data at the survey point level, which is consistent with the scale of conservation treatments that we apply, and avoided information loss caused by pooling data at the BBS route level. By accounting for observer effects, nesting of survey points within routes, and sequence of survey points, we accommodated BBS survey design, refined estimates of important habitat predictors, improved model fit, and reduced or eliminated positive spatial autocorrelation in model residuals. The predictive power of models was greatly increased by including variables that characterized annual and long-term precipitation, as well as local land cover attributes not available from satellite-derived land cover data. Occurrence models from survey-point-level BBS data and environmental data with high thematic resolution were able to describe habitat relationships that are often associated with fine-grained, local studies, but across broad spatial extents and at scales relevant to local conservation actions. Predicted occurrence was strongly correlated with observed numbers, suggesting that occurrence models may be useful indicators of density. Relationships derived from models allowed us to develop spatially explicit decision support tools, which can be used to target areas for conservation treatments and to assess the conservation actions of multiple conservation programs and joint ventures (e.g., Prairie Pothole, Rainwater Basin, and Northern Great Plains joint ventures) in the U.S. Northern Great Plains.
INTRODUCTION
Concern over decreasing bird populations has stimulated a variety of bird conservation plans, many of which (e.g., North American Waterfowl Management Plan, Partners In Flight, The Nature Conservancy's Migratory Bird Program) explicitly promote a landscape approach to bird conservation. Increasing awareness of the importance of landscape composition to avian ecology and conservation, in conjunction with a recent upsurge in the availability of spatial analysis software and data, has led to increased development and application of spatially explicit models to direct conservation actions (Carroll et al. 1996, Askins 2000, Wiens 2002). These models, often referred to as spatial planning or conservation assessment tools, are used for a variety of purposes, including identification of habitat and lands for protection, prioritization of funding, and identification of opportunities for restoration.
Spatial tools for guiding bird conservation may be particularly important in the Northern Great Plains, which have the highest diversity of grassland bird species in North America (Peterjohn and Sauer 1999) as well as 6 species of endemic passerine (Samson et al. 1998). Rich soils and limited topographic relief also make the Northern Great Plains an important area for crop production, and native grasslands in the region are among the most threatened ecosystems in the world, especially in the eastern portion of the region where precipitation supports more crop varieties (Licht 1997, Hoekstra et al. 2005). Conversion of grassland, particularly native prairie, to cropland in the region is extensive and ongoing as agricultural subsidies, new crop varieties, and altered climate enable the planting of lands that were previously considered unsuitable for crop production (Stephens et al. 2008, Rashford et al. 2011, Lark et al. 2015). Habitat loss for grassland birds is exacerbated by roads, shelterbelts, wind turbines, and oil and gas infrastructure that fragment the landscape and reduce habitat suitability for grassland birds (Grant et al. 2004, Shaffer and Buhl 2016, Thompson et al. 2015). As a consequence of habitat loss and degradation, grassland birds have a larger proportion of species that are decreasing than any other bird group in North America (Askins 1993, Peterjohn and Sauer 1999, Sauer et al. 2017).
The need for spatial tools that can be used to evaluate, allocate resources to, and increase efficiency of conservation actions in the Northern Great Plains is magnified by the sheer size of the region, extensive private land ownership, and the variety of available conservation treatments. The Great Plains region of North America covers ∼162 million ha (Samson and Knopf 1994) and exhibits considerable variation in climate, topography, soil quality, and land use (Licht 1997, Samson et al. 2004). Most land in the Great Plains region is privately owned, and many conservation programs address the differing needs and interests that motivate people who own the land (Heard 2000, Ryan et al. 2003, Ernst and Wallace 2008). Because of varying interests of landowners and the diversity of land types and uses in the region, an array of conservation treatments is available to benefit grassland birds, including acquisition of perpetual conservation easements to preserve existing grasslands, as well as delayed haying, planting of cropland to grassland, tree and brush removal, prescribed burns, and grazing management to enhance or restore habitat (Gray et al. 2005, Johnson 2005, USFWS 2012).
We used data from the North American Breeding Bird Survey (BBS) in conjunction with environmental predictors to develop comprehensive, species-specific spatial planning tools for guiding grassland bird conservation in the U.S. Northern Great Plains. The BBS is an annual, continent-wide survey that is the primary source of information regarding populations of many North American bird species, thanks to the efforts of thousands of volunteer observers combined with the scientific rigor of the survey and analysis of resulting data (Bystrak 1981, Sauer et al. 2013). Despite not being intended for the development of spatial models, the consistent sampling framework, long timeframe, widespread distribution of survey routes, and variety of habitat types and land uses that the BBS encounters make BBS data valuable for developing spatial models as well as for monitoring avian population trends (Niemuth et al. 2005, Thogmartin et al. 2006a, Sauer et al. 2013, Gorzo et al. 2016, Sauer et al. 2017). Our study had 3 main objectives: (1) to identify factors, especially landscape characteristics, associated with the detection of grassland birds at BBS ‘stops' (individual survey points along a BBS route); (2) to create maps showing predicted occurrences of grassland birds across our study region; and (3) to use relationships identified in models to create additional decision support tools to guide conservation actions for grassland birds in the U.S. Northern Great Plains.
METHODS
Study Area
Our study area included the states of North Dakota, South Dakota, Nebraska, Montana, and Kansas, as well as those portions of Colorado and Wyoming east of the Rocky Mountains (Figure 1). This region, which we refer to as the U.S. Northern Great Plains, is characterized by relatively flat topography and limited rainfall that follows an east–west gradient, with higher precipitation in the east (Wiens 1974). Because water is generally limiting in this semiarid landscape, the precipitation gradient greatly influences land use, vegetation composition and structure, and bird communities (Wiens 1974, Samson et al. 1998, Niemuth et al. 2008). Much native grassland has been converted to crop production, with losses of native prairie exceeding 99% in the eastern portion of the region (Samson and Knopf 1994, Licht 1997). In addition to cropland, the study area has more trees and woody vegetation than it did historically as a result of fire suppression, altered grazing regimes, tree planting, and alteration of hydrologic regimes following settlement by Euro-American immigrants in the 1800s (Licht 1997, Courtwright 2007). The U.S. Northern Great Plains region also encompasses millions of hectares of grasslands that have been enrolled in the U.S. Department of Agriculture (USDA) Conservation Reserve Program (CRP), which substantially benefits the populations of many grassland bird species (Johnson and Igl 1995, O'Connor et al. 1999, Johnson 2005).
BBS Data
We downloaded stop-level BBS data from 2005 to 2011 for routes within our study area from the U.S. Geological Survey Patuxent Wildlife Research Center, Laurel, Maryland, USA (Pardieck et al. 2014). Each 40-km route contained 50 stops, or survey points, ∼0.8 km apart; details of route placement and sampling are described by Bystrak (1981). We assigned the resultant 16,250 stops to geographic coordinates using a variety of techniques, with 55% of locations coming from observer-provided information, including GPS locations, field descriptions, and digitization of stops marked on aerial photographs by observers, and 45% of locations coming from automated generation of points at 0.8-km intervals from the starting point along individual survey routes. The accuracy of locations of stops assigned at 0.8-km intervals was likely aided in our study area by the fact that many of the survey routes followed roads laid out on a 1.6-km grid based on the public land survey. We selected the 2005–2011 timeframe for bird survey data as it overlapped with the time period of land cover data collection and provided a broad range of precipitation conditions. We analyzed data from 83,500 counts collected at the 16,250 stops along 325 routes by 264 observers, only using data that passed BBS quality criteria for weather conditions, daily timing, and seasonal timing (see www.pwrc.usgs.gov/BBS/help/BBS_Run_Type_Codes.txt for more information).
We analyzed data for the Upland Sandpiper (Bartramia longicauda), Sprague's Pipit (Anthus spragueii), Lark Bunting (Calamospiza melanocorys), Savannah Sparrow (Passerculus sandwichensis), Grasshopper Sparrow (Ammodramus savannarum), Bobolink (Dolichonyx oryzivorus), and Eastern Meadowlark (Sturnella magna), as these species have been identified as conservation priorities (Rosenberg et al. 2016), exhibit a variety of grassland habitat preferences and geographic distributions, and had sufficient observations with which to develop models.
Predictor Variables
Because many factors affect observations of birds, we developed models from a suite of candidate predictor variables that characterized landscape composition and configuration, weather and climate, daily and seasonal changes in bird activity and detectability, topographic variation, and survey structure, all of which have been well supported by previous research (Table 1). Land cover data were derived in part from the National Land Cover Database 2006 (NLCD 2006; Fry et al. 2011). NLCD 2006 has overall agreement of 78% between classified satellite data and a primary or alternate land cover class visually interpreted from aerial photography, although accuracy has been consistently lower among grass-dominated classes (Wickham et al. 2013). To improve thematic resolution and classification accuracy of grass-associated land cover data, we incorporated spatial data from the USDA National Agricultural Statistics Service identifying alfalfa (Medicago sativa) fields (Boryan et al. 2011), as well as data delineating 3.8 million ha of land enrolled in CRP grasslands, which were mapped rather than interpreted from remotely sensed imagery. All predictor data were processed at a spatial resolution of 30 m.
TABLE 1.
Predictor variables considered in the development of models predicting the occurrence of grassland birds in the U.S. Northern Great Plains were selected based on documented associations with bird presence, density, or detection. All predictors were treated as continuous variables unless otherwise noted. NLCD class = land cover class derived from the National Land Cover Database 2006 (Fry et al. 2011).
We obtained precipitation and temperature data from the PRISM (Parameter-elevation Regressions on Independent Slopes Model) climate mapping database, which uses weather station data to model precipitation and temperature across space (Daly et al. 2008). Previous-year and current-year precipitation were strongly correlated with long-term precipitation because they generally followed the same east–west gradient as long-term precipitation. Therefore, we subtracted previous-year and current-year precipitation from the long-term mean to create a variable reflecting the precipitation anomaly for each time period.
Because changes in topography may influence the detection (Dawson 1981) or densities of birds (Renfrew and Ribic 2002), we included the standard deviation of elevation around each survey stop as an index of topographic roughness; we also included mean elevation to capture gradients that might be associated with topography or soil characteristics. We included the number of each stop along an individual BBS route as an index to the time of day, thereby mitigating daily time-related changes in bird detection, which varies during the interval required to run a BBS route (Robbins 1981, Rosenberg and Blancher 2005). Similarly, we included ordinal date as a covariate to explain seasonal changes in bird detectability (Anderson et al. 1981, Skirvin 1981). We also included year as an indicator variable to account for interannual variation in population size and distribution not attributable to changes in observers or patterns of annual precipitation. Even though our objective was to develop spatial models to predict occurrence across a regional landscape, we included nonlandscape factors such as annual precipitation, daily timing, and seasonal timing to explain additional variation in the data, thus improving our ability to make inferences about landscape-level habitat selection. Because of repeated observations along routes across multiple years and differences in the skills of observers, some of whom ran multiple routes over multiple years, we treated route, observer, and year as random effects to address changes in variance associated with these variables (Crawley 2007). All other variables were treated as fixed effects. We did not include a first-year observer effect (Kendall et al. 1996) because we assumed that this effect would be less problematic for detecting presence than estimating population trends, and not doing so resulted in a simpler model.
Because many bird species are influenced by the landscape beyond the area included in the point-count circle (Bakker et al. 2002, Cunningham and Johnson 2006, Greer et al. 2016), we sampled the habitat around each BBS stop at 7 scales using circular moving window analysis, which summarizes data within a ‘window' of a selected size around each cell in a GIS data layer. Landscape data were in raster format, and the area within each moving window was ∼50, 200, 450, 800, 1,250, 1,800, and 3,200 ha, respectively, for circles with radii of ∼400, 800, 1,200, 1,600, 2,000, 2,400, and 3,200 m. We chose these scales as they were multiples of the maximum survey distance used in the BBS and also coincided with distances commonly used in land acquisition and management in the region. We chose not to use finer scales for 2 reasons. First, even though many species have detection distances much less than the nominal 400-m sampling window of the BBS (Thogmartin et al. 2006b), the locations of recorded individuals within the window were unknown and may have been outside a circle with a smaller radius. Second, the locations of some BBS stops were imperfectly known, and maintaining a broad sampling window helped to ensure that stop locations at which bird data were collected coincided with sampled environmental predictors. We standardized all continuous variables to a mean of 0 and standard deviation of 1 to improve convergence of the model-fitting algorithm. We analyzed spatial data using ArcMap (Environmental Systems Research Institute, Redlands, California, USA).
Model Development
Because most of the species that we evaluated were not detected or were detected in low numbers at BBS stops, we used logistic regression to model apparent occurrence. Even though our models were biologically justified and well supported by past research, we used model selection to develop a parsimonious model that suitably balanced bias and variance (Burnham and Anderson 2002), as well as to evaluate models developed with different combinations of correlated variables. Prior to developing models, we assessed collinearity among predictor variables to ensure that highly correlated (Pearson's r > 0.7) variables were not considered simultaneously. We began by developing a null model that included only the intercept and random effects (Crawley 2007), which served as a baseline for assessing improvements in model fit based on changes in Akaike's Information Criterion (AIC), followed by a full model containing all predictor variables. The full model was run at each of the candidate spatial scales, and we selected the scale with the lowest AIC value. We discriminated among reduced versions of the full model, holding out one parameter or set of parameters at a time and assessing improvements in AIC values to select a best approximating model (Burnham and Anderson 2002, Crawley 2007). When the full model would not converge, we used different subsets of the full model to evaluate predictors and identify the model best supported by the data. We calculated Akaike weights (wi) for each model within 4 AIC units of the model with the lowest AIC value, which is a useful rule of thumb for identifying the set of models plausibly supported by the data (Burnham and Anderson 1998). Akaike weights provide an indication of the relative likelihood of competing models best fitting the data, and thus enable evaluation of the relative strength of evidence for models relating bird observations to predictor variables.
In an attempt to develop a parsimonious model and avoid spurious correlations, we only evaluated main effects of linear relationships, except for quadratic relationships that characterized climatic envelopes and nonlinear relationships with the amount of cropland in the landscape or seasonal changes in detection (Table 1). We conducted statistical analyses in the R environment (R Core Team 2013), specifically the generalized linear mixed models capacity of the lme4 package (Bates et al. 2015), using a binomial distribution. We used the bound optimization by quadratic approximation option to improve convergence of the maximum likelihood estimator.
Because geographic distributions varied among species, we did not use the same analysis extent for all species. For species that were not distributed across our entire study area, we selected analysis areas by states, as that is the level at which many conservation programs in the region are implemented. By excluding large areas where species did not occur, we were able to reduce the preponderance of zeroes and resulting overestimation of model performance metrics (Lobo et al. 2008, Barve et al. 2011, Zuur et al. 2012) while retaining sufficient observations where birds were not detected to model biologically important climatic factors influencing species distributions (Guisan and Thuiller 2005).
Analyzing BBS data at the stop level allows inferences to be made at a much finer spatial resolution than using BBS data at the route level, but increases the potential for positive spatial autocorrelation, which, if ignored, can lead to overestimation of the precision of parameter estimates, obscure ecological patterns, and reduce model performance (Legendre 1993, Carroll and Pearson 2000, Lennon 2000, Lichstein et al. 2002). We included climatic and land cover variables to account for broad-scale gradients that may influence bird distribution, as well as observer and time-of-day variables to account for local spatial autocorrelation. When spline correlograms (Bjornstad 2015) indicated that positive spatial autocorrelation remained in model residuals, we reran the best-supported model with an autologistic term that indicated the presence of the target species at adjacent stops to improve model fit and reduce local autocorrelation (Augustin et al. 1996, Klute et al. 2002). When creating correlograms, some of which used correlation matrices resulting from 83,500 observations, computing limitations required that we thinned residuals from 4 of the models by year. We evaluated models by calculating the area under the curve (AUC) of receiver operating characteristics (ROC; Hosmer and Lemeshow 2000) using 10-fold cross validation (Stone 1974).
We created maps showing the relative predicted occurrence of each species throughout the study region by incorporating corresponding GIS data into the logistic regression equation, using coefficients estimated from all folds of the data used to develop the model. We used the mean value of nonlandscape variables (i.e. those related to detection or annual weather conditions) when applying models to landscape data. Because of the difficulty of applying the autologistic term across the landscape, particularly in an environment as variable as our study region, we used the autologistic term to improve statistical inference but did not apply it to create predictive surfaces (Boyce 2006, Dormann et al. 2007). We also created plots, by species, of bird response to the amount of perennial cover (i.e., pasture and hay, grassland and herbaceous, CRP, and alfalfa cover classes) and forest in the sampling window, again holding other variables at their mean value. These plots were used to compare species' responses to these factors and to identify thresholds in responses to landscape characteristics that could be addressed at local scales through conservation treatments such as tree removal and grassland restoration. Finally, we assessed simple correlations between predicted occurrence and number of birds observed at each stop, by species, to determine whether occurrence models were useful predictors of density.
RESULTS
Landscapes surrounding BBS stops throughout our study region varied considerably in type and distribution of land cover (Table 2). High correlations between forest cover and topographic roughness (r = 0.72) and long-term January (minimum) and August (maximum) temperatures (r = 0.70) precluded simultaneous consideration of these variables in models using the complete dataset. In the subset of data from Kansas and Nebraska, cropland and the grassland and herbaceous cover class were strongly correlated (r = −0.70), precluding their simultaneous consideration in the Eastern Meadowlark model, which was constrained to those 2 states. Data were dominated by zeroes for all species, although prevalence varied among species (Table 3). Given the complexity of our models and the number of variables considered, some models that we considered did not successfully converge, even when the number of maximum likelihood iterations was increased to 500,000.
TABLE 2.
Means and standard deviations (SD) for continuous predictor variables at 83,500 Breeding Bird Survey (BBS) counts conducted at 16,250 stops (individual survey points). Values for land cover and digital elevation model data were derived from a sampling window with 800-m radius. Land cover data were static, but climatic and temporal data varied among years. See Table 1 for variable definitions.
TABLE 3.
Species, scale of model, model performance (area under curve [AUC] of receiver operator characteristics), difference in Akaike's Information Criterion (AIC) from null model (Δn AIC), correlation between predicted occurrence and individuals actually observed (Cor), U.S. states included in analysis, number of Breeding Bird Survey (BBS) stop (survey point) counts included in analysis (n), and number of counts during which each species was detected (Detections) for best-supported models predicting the occurrence of 7 grassland bird species in the U.S. Northern Great Plains, 2005–2011. Variables are defined in Table 1.
Habitat and observed bird numbers showed strong positive spatial autocorrelation, but spatial autocorrelation was eliminated in model residuals (Figure 2) for 4 of the 7 species that we assessed. Climatic and land cover variables accounted for much spatial autocorrelation (Figure 2C), but observer and time variables were necessary to remove remnant spatial autocorrelation (Figure 2D). Models for the Upland Sandpiper, Lark Bunting, and Grasshopper Sparrow also required the addition of an autologistic term to remove remnant positive spatial autocorrelation from model residuals.
FIGURE 2.
Positive spatial autocorrelation was evident in (A) the amount of grassland in the landscape surrounding Breeding Bird Survey (BBS) stops (individual survey points) and (B) the number of Eastern Meadowlarks recorded at BBS stops in Kansas and Nebraska, USA. Positive spatial autocorrelation was (C) substantially reduced in residuals of a model predicting the occurrence of Eastern Meadowlarks in Kansas and Nebraska that included only habitat, climatic, and topographic variables, and (D) eliminated from residuals of a model predicting Eastern Meadowlark occurrence that also included observer, BBS stop, and date. The heavy solid line represents estimated autocorrelation, and the thin dashed line indicates the 95% confidence envelope. Data and models for other species and geographic extents showed similar patterns.
The best-supported models characterizing bird–environment relationships indicated that the occurrence of all species was influenced by climate, weather, or topography, as well as landscape composition and configuration (Table 4). Improvements in AIC values over the null model indicated substantial support for the best-supported model for all species, and AUC values ranged from 0.80 to 0.95 (Table 3), indicating excellent to outstanding discrimination (Hosmer and Lemeshow 2000). Model uncertainty varied among species, but competing models were nested and often differed from the best-supported model due to the inclusion or exclusion of only one variable (Appendix Table 5). The focal species showed similar responses to landscape characteristics, with consistent negative associations with trees, positive and varying associations with grassland cover classes, and negative, weak positive, or curvilinear responses to cropland (Table 4, Figures 3 and 4).
TABLE 4.
Variables and estimated coefficients (and standard errors) for landscape models predicting the occurrence of 7 grassland bird species in the U.S. Northern Great Plains, 2005–2011. Variables are defined in Table 1, except Autologistic = a binary term indicating the presence or absence of the target species at one or both adjacent survey points, added to improve model fit and reduce local spatial autocorrelation.
FIGURE 3.
Response to the area of forest in the sampling window varied among species, with Lark Bunting and Upland Sandpiper showing the strongest avoidance of trees, and Grasshopper Sparrow showing the lowest avoidance of trees. Response curves were scaled to a common unit for all species. The order of species in the legend follows the order of species in the figure: GRSP = Grasshopper Sparrow, BOBO = Bobolink, EAME = Eastern Meadowlark, SAVS = Savannah Sparrow, LARB = Lark Bunting, and UPSA = Upland Sandpiper.
FIGURE 4.
Response to the amount of grassland in the sampling window varied among species, with Grasshopper Sparrow and Savannah Sparrow showing the greatest and Lark Bunting and Sprague's Pipit showing the smallest increases in occurrence as the amount of grassland in the landscape increased. The amount of each grassland type (grassland and herbaceous, pasture and hay, Conservation Reserve Program [CRP] grassland, and alfalfa cover classes) was equally divided among the cover classes included in the best-supported model for each species. The order of species in the legend follows the order of species in the figure: GRSP = Grasshopper Sparrow, SAVS = Savannah Sparrow, UPSA = Upland Sandpiper, BOBO = Bobolink, EAME = Eastern Meadowlark, LARB = Lark Bunting, and SPPI = Sprague's Pipit.
Precipitation strongly influenced the occurrence of all 7 species, with 5 of the 7 species influenced by short-term (either current-year or previous-year) precipitation and 6 influenced by long-term (30-yr mean) precipitation. The occurrence of Upland Sandpipers, Lark Buntings, and Eastern Meadowlarks was more strongly associated with mean long-term January (minimum) temperature than mean long-term August (maximum) temperature (Table 4). The detection of all species but Sprague's Pipit was influenced by the daily and/or seasonal timing of surveys, as well as survey structure, including observer, year, and route effects (Table 4).
Spatial patterns in the predicted occurrence and numbers of grassland birds reflected regional climatic patterns, land forms, and cover classes, with Sprague's Pipits and Savannah Sparrows selecting dry and moist portions, respectively, of northern states; Upland Sandpipers, Lark Buntings, and Grasshopper Sparrows found throughout much of the study area; and Eastern Meadowlarks occurring most frequently in the moister, eastern portion of Nebraska and Kansas (Figure 5). Consistent with these patterns, the best-supported models showed a negative relationship between the occurrence of Lark Buntings and the area of emergent herbaceous wetlands, and a positive relationship between the area of emergent herbaceous wetlands and occurrence of Upland Sandpipers, Savannah Sparrows, Bobolinks, and Eastern Meadowlarks (Table 4). Predicted occurrence at each BBS stop was strongly correlated with observed numbers for all species (Table 3).
DISCUSSION
Our results demonstrate that analyses using stop-level BBS data and environmental data with high thematic resolution are able to describe habitat relationships often associated with fine-grained, local studies but across broad spatial extents and at scales relevant to local conservation actions. For example, our models indicated that the Savannah Sparrow was positively associated with pasture and hay, which was found primarily in the northeastern, or tallgrass prairie, portion of our study region, CRP grasslands, and emergent wetlands, all of which are consistent with previous findings of selection for mesic sites, tall, dense cover, and exotic grasses (Davis and Duncan 1999, Madden et al. 2000, Davis et al. 2016). Bobolinks showed a similar response, but were also associated with alfalfa, again consistent with previous findings showing selection for exotic grasses and legumes (Renken and Dinsmore 1987, Delisle and Savidge 1997). Conversely, the strong association of Sprague's Pipits, Grasshopper Sparrows, Lark Buntings, and Upland Sandpipers with the grassland and herbaceous cover class, which was found primarily in the central and western portion of our study region, is consistent with previous findings that these species generally select drier sites with short or sparse vegetation (Davis et al. 1999, Madden et al. 2000, Lueders et al. 2006).
The association between the area of land enrolled in CRP grasslands and the occurrence of the Lark Bunting, Savannah Sparrow, Grasshopper Sparrow, Bobolink, and Eastern Meadowlark reinforces previous findings as well as the importance of the CRP program to grassland bird populations in the Great Plains (Johnson and Igl 1995, Delisle and Savidge 1997, Johnson 2005). A lack of association between CRP grassland and the occurrence of Sprague's Pipits and Upland Sandpipers reflects the Sprague's Pipit's selection for native grasslands of short to intermediate stature (Davis and Duncan 1999, Davis et al. 1999, Madden et al. 2000, Davis et al. 2016) and the Upland Sandpiper's frequent selection of sites with short, sparse vegetation (Renken and Dinsmore 1987, Sandercock et al. 2015). As expected, some of the species that we assessed showed quadratic or weak positive associations with cropland, which is consistent with previous findings of lower density or likelihood of occurrence in cropland than in grasslands (Johnson and Igl 1995).
Responses to climate varied among species but, similarly to other studies (i.e. Thogmartin et al. 2006a, Ahlering et al. 2009, Albright et al. 2010, Lipsey et al. 2015), precipitation and/or temperature were strong predictors of occurrence for all species. The biological significance of climatic variables is unclear, as they may be correlates of other factors (e.g., plant community composition, primary and secondary productivity) that more directly influence species occurrence, likely in concert with other factors such as soils and landform (Guisan and Zimmerman 2000, Niemuth et al. 2008). The occurrence of 4 of the species that we assessed was more strongly associated with long-term mean August temperatures, while the occurrence of the remaining 3 species was more strongly associated with long-term mean January temperatures, but the mechanism responsible for that difference, and whether the difference was real and not an artifact of correlation between the 2 variables, is unknown. Regardless of mechanism, weather and climate in our study region are highly variable and strongly affect bird occurrence, whether directly or indirectly.
We did not find support for an association between the occurrence of Sprague's Pipits and the number of patches in the landscape, even though previous analyses have found Sprague's Pipits to be sensitive to landscape fragmentation (Davis 2004, Lipsey et al. 2015), nor did we find associations between Sprague's Pipit occurrence and stop number or ordinal date, which were present for all other species that we considered. Lack of support for these relationships may be a function of the small number of observations of Sprague's Pipits, which had <10% of the detections of the other species that we considered. The Sprague's Pipit is simply an uncommon species throughout much of its range, but the problem of the small number of detections was addressed in part by the 2015 addition of 42 BBS routes in Montana, which had the lowest BBS route density (1 route per degree block) and highest Sprague's Pipit density in the United States.
The BBS only provides an index to bird presence and numbers, as existing protocols provide no mechanism for assessing and correcting for detectability, and some unknown fraction of the birds present at each stop is not recorded (Sauer et al. 2013). Nevertheless, uncorrected data can still provide useful estimates of relative density or probability of occurrence (Johnson 2008, Etterson et al. 2009, Leston et al. 2015), and spatial models developed from BBS data have been useful for providing ecological insights, guiding conservation, and providing spatially explicit minimum estimates of population size and distribution (e.g., Newbold and Eadie 2004, Thogmartin et al. 2006a, Hudson et al. 2017, Rosenberg et al. 2017). Predicted occurrence was positively and significantly correlated with observed counts for all of the species that we considered, suggesting that the occurrence models that we present are also useful for identifying areas of high density.
Our models included several variables (i.e. stop number, ordinal date, current-year precipitation, previous-year precipitation, and autologistic) that weren't applied to spatial data to create maps showing relative probability of occurrence. These variables explained spatiotemporal or fine-grained spatial variation in bird occurrence that improved estimates of variables that were in line with our goal of developing landscape-scale predictive models over broad spatial and temporal extents. Models that include variables to accommodate observer and route effects as well as daily and seasonal timing can have AIC values >100 points lower than models without such variables (data not shown), indicating that models that do not accommodate sampling and design issues have essentially zero support for adequately describing the data relative to models that contain these variables (Burnham and Anderson 1998). In addition, the elimination of spatial autocorrelation of residuals when timing and observer variables were included suggests that our modeling process accounted for spatiotemporal patterns in detection caused by observer and timing effects.
Interestingly, of the 3 species that required an autologistic term to reduce spatial autocorrelation in model residuals, 2 species, the Lark Bunting and Upland Sandpiper, are thought to be colonial or semicolonial nesters (Shane 2000, Casey et al. 2011). This suggests that some of the spatial autocorrelation that we observed may have been rooted in bird behavior rather than habitat or sampling, which reinforces the appropriateness of an autologistic term to capture such dynamics. However, autologistic regression contains a degree of circularity and reduces the size of coefficient estimates for habitat variables (Dormann 2007), which complicates the application of the models to conservation. In our analyses, confidence intervals for some environmental variables included zero due to a reduced size of coefficient estimates and/or increased standard errors after the autologistic term was added. We chose to retain these variables, given their biological importance and selection in the nonautologistic models; alternatively, one could simply use models without the autologistic term, treating remaining autocorrelation as a behavioral artifact beyond the scope of management actions.
The radius of the sampling window at which landscape data best described bird occurrence was ≤800 m for 5 of the 7 species that we evaluated, but extended to 1,200 m for the Sprague's Pipit and 3,200 m for the Lark Bunting. Our findings are consistent with other studies which have shown that landscape characteristics influence the occurrence or density of grassland birds and that the scale of the landscape influence varies among species (Ribic and Sample 2001, Cunningham and Johnson 2006, Thogmartin et al. 2006a). Birds likely respond to different landscape features (e.g., trees vs. wetlands) at different scales, but we did not assess landscape characteristics at multiple scales within individual species' models due to the absence of a priori information about selection preferences of each species.
Management Implications
Spatially explicit models provide a biological foundation for identifying landscapes suitable for protection or management, as well as for assessing the effects of conservation programs and investigating the potential effects of changes in climate and land use. The relative scarcity and limited distribution of some species reinforce the importance of using spatial models to direct conservation efforts, as conservation treatments in areas without the appropriate climatic envelope or landscape characteristics will provide little benefit for target species. The models presented in this paper are of sufficiently fine spatial and thematic resolution to assess individual land parcels, unlike models developed using coarse-grained response data (i.e. entire BBS routes) or predictor variables. However, even with relatively fine resolution, management may be necessary to ensure that appropriate fine-grained habitat features (i.e, absence of trees, appropriate vegetation structure and composition) are present (Grant et al. 2004, Derner et al. 2009, Greer et al. 2016).
In our study region, spatial models and decision support tools derived from those models are widely used to guide conservation efforts (Niemuth et al. 2009, RWBJV 2013). Paper and digital copies of occurrence or density models are distributed to conservation practitioners, who use them to evaluate landscapes for conservation treatments. For example, value to grassland birds is one of the criteria for assessing candidate land parcels for acquisition of perpetual grassland easements in the Prairie Pothole Region (USFWS 2010), where tens of millions of dollars are spent annually to conserve habitat for upland-nesting waterfowl. As a result of these efforts, ∼2 million ha of wetlands and grasslands have been protected in the Prairie Pothole Region through the acquisition of perpetual easements for waterfowl conservation (Niemuth et al. 2014). In the Flint Hills of Kansas, the responses of grassland birds to tree and grassland cover depicted in Figures 3 and 4 were used to develop spatially explicit decision support tools showing areas where tree removal and grassland restoration would provide the greatest benefits to grassland birds (M. Estey personal observation).
In the absence of direct access to applied models, the varying responses that we have documented for response to the amount of grassland or forest cover in the surrounding landscape may provide a framework for providing benefits for multiple species. Whereas many species richness models focus on areas of distributional overlap without considering species requirements and conservation treatments, the relationships that we have identified allow practitioners to identify portions of the landscape needing treatment (e.g., tree removal or grassland restoration) and, by meeting the needs of the most restrictive species, to provide habitat for multiple species. The occurrence of 4 species was positively associated with wetlands in the landscape, which provides justification for the restoration of grassland and wetland complexes for migratory bird conservation in the U.S. Northern Great Plains. Finally, negative responses by grassland birds to urban areas and grassland fragmentation provide justification for conservation easements and grassland restoration that prevent development and reduce fragmentation, respectively.
Climatic conditions in the Northern Great Plains are highly variable, with the result that the distributions and numbers of birds can change greatly from one year to the next (Cody 1985, George et al. 1992, Niemuth et al. 2008). Variability in distributions reinforces the importance of broad spatial extents and long timeframes in conservation planning and action; the BBS is well suited for providing data to help guide conservation actions for many species across much of North America.
ACKNOWLEDGMENTS
We thank the many BBS observers who collected data and helped to identify stop locations; K. L. Pardieck and D. Ziolkowski of USGS Patuxent Wildlife Research Center for providing BBS stop data and route information; J. R. Sauer of USGS Patuxent Wildlife Research Center for discussion about observer biases; T. L. Shaffer for discussion about modeling; K. L. Pardieck, R. D. Pritchert, J. A. Shaffer, H. M. Specht, and 2 anonymous reviewers for comments on an earlier draft of this paper; K. W. Barnes and C. F. Jorgensen for providing R code and advice; and the USDA for providing CRP data. The findings and conclusions in this article are those of the authors and do not necessarily represent the views of the U.S. Fish and Wildlife Service.
Funding statement: This study was undertaken by all authors as part of their normal duties.
Ethics statement: We followed the Code of Ethics of the Ecological Society of America.
Author contributions: N.D.N., M.E.E., A.A.B., and S.P.F. conceived the idea; S.P.F., B.W., M.E.E., and N.D.N. developed methods; N.D.N., M.E.E., B.W., P.J.M., R.C.G., and A.J.R. analyzed the data; and N.D.N. wrote the paper.
LITERATURE CITED
Appendices
APPENDIX TABLE 5.
Constituent variables (preceded by a sign indicating the direction of the relationship), differences in Akaike's Information Criterion (ΔAIC), and Akaike weights (wi) for candidate models with AIC differences <4 relating the apparent occurrence of 7 grassland bird species to environmental and survey predictors in the U.S. Northern Great Plains. All models contain random effects for observer, Breeding Bird Survey (BBS) route, and year. Variables are defined in Table 1.
