Distribution of priority grassland bird habitats in the Prairie Pothole Region of Canada

Grassland ecosystems and the species that rely on them are one of the most urgent habitat conservation concerns in North America. Fundamental to any landscape conservation efforts is the identification of priority habitats to help target management efforts. Many avian species associated with prairie ecosystems have experienced population declines along with continued loss of prairie habitats. Additionally, given the long history of research in avian systems and the close grassland associations of some species, birds are excellent candidate taxa for the identification of priority habitats and can provide an informed starting point for multispecies assessments. We used data from the North American Breeding Bird Survey (1997-2014) to develop species distribution models for 15 grassland bird species across the Prairie Pothole Region of Canada. Model performance varied widely across species. Ten species demonstrated good model performance (average Boyce Index > 0.64 across 5-fold cross validation). We used these 10 species to assess the influence of habitat covariates on the relative probability of occurrence, to compare the spatial scales of selection, and to generate multispecies habitat priority maps. Of the nine habitat covariates considered, most species predictably demonstrated positive associations with grassland habitats and avoidance of areas of high tree and shrub cover. Two covariates representing wetland abundance were also frequently included in the top models. The area covered by wetlands (w.area) was present in the top model for 5 of 10 species with a consistently estimated negative coefficient. However, a covariate, which represented the number of wetland basins (w.basins), was present in the top model for 8 of 10 species with an estimated positive coefficient for all but 1 species, representing a preference for more heterogeneous wetland landscapes. The larger spatial scales we considered tended to have greater explanatory power than smaller spatial scales and were thus more prevalent in the top models. The multispecies priority habitat maps that we produced can be used for future assessments of potential habitat management actions. Our work provides a critical foundation for the incorporation of grassland bird conservation goals into on-going landscape-planning initiatives in the Prairie Pothole Region of Canada. Répartition de milieux prioritaires pour la conservation d'oiseaux de prairie dans la région des cuvettes des Prairies au Canada RÉSUMÉ. Les écosystèmes de prairie et les espèces qui en dépendent représentent une des préoccupations de conservation d'habitat les plus importantes en Amérique du Nord. Au moment d'établir les activités d'aménagement, la détermination des milieux prioritaires est fondamentale à tout effort de conservation de paysages. De nombreuses espèces aviaires associées aux écosystèmes de prairie ont subi des baisses de population et des pertes continuelles d'habitat. De plus, étant donné la longue tradition de recherche en ornithologie et l'association étroite de certaines espèces avec les milieux de prairie, les oiseaux représentent un excellent taxon pour la détermination de milieux prioritaires et peuvent également offrir un point de départ éclairé pour les évaluations multispécifiques. Nous avons utilisé les données du Relevé des oiseaux nicheurs (ou BBS) en Amérique du Nord (1997-2014) afin d'élaborer des modèles de répartition pour 15 espèces d'oiseaux de prairie habitant la région des cuvettes des Prairies au Canada. La performance des modèles a beaucoup varié d'une espèce à l'autre. La performance du modèle était bonne pour dix espèces (moyenne de l'indice de Boyce > 0,64 lors d'une validation croisée répétée 5 fois). Nous nous sommes servis de ces dix espèces pour évaluer l'influence des covariables d'habitat sur la probabilité de présence relative, comparer les échelles spatiales de sélection et produire des cartes de milieux prioritaires pour plusieurs espèces. Neuf covariables d'habitat ont été considérées. Sans surprise, la plupart des espèces ont montré une association positive avec les milieux de prairie et un évitement des secteurs où le couvert en arbres et en arbustes était élevé. Deux covariables liées aux milieux humides se retrouvaient souvent dans les meilleurs modèles. La superficie de milieux humides (w.area) figurait dans le meilleur modèle pour cinq des dix espèces, avec un coefficient constamment négatif. Toutefois, la covariable représentant le nombre de milieux humides (w. basins) était présente dans le meilleur modèle pour huit des dix espèces, avec un coefficient positif pour toutes les espèces sauf une, ce qui indique une préférence pour des paysages de milieux humides plus hétérogènes. Les plus grandes échelles spatiales que nous avons considérées avaient en général une capacité d'explication plus élevée que les plus petites échelles spatiales, et se retrouvaient, par conséquent, plus souvent dans les meilleurs modèles. Les cartes de milieux prioritaires que nous avons produites pour plusieurs espèces à la fois peuvent être utilisées dans le cadre d'une éventuelle évaluation d'activités potentielles d'aménagement d'habitat. Nos travaux établissent une assise fondamentale pour que soient inclus des objectifs de conservation d'oiseaux de prairie dans les initiatives de planification à l'échelle du paysage ayant cours dans la région des cuvettes des Prairies du Canada.


INTRODUCTION
Grassland ecosystems are some of the most imperilled in the world (White et al. 2000, Hoekstra et al. 2005 with less than 1% of original prairie habitat remaining in many regions of the North American plains (Samson and Knopf 1994). In Canada, estimates published in 2003 indicated that 25-30% of native grasslands remained, primarily concentrated in southern Alberta and Saskatchewan (Gauthier and Wiken 2003). Loss of remaining native grassland in the prairies of Canada continues at an annual rate of approximately 0.7% (Watmough and Schmoll 2007). Although primary threats are conversion and loss due to crop production (Coupland 1979, Vitousek et al. 1997, Fargione et al. 2008, additional threats include fire suppression, poor grazing management, and industrial and urban development (Askins et al. 2007). Wildlife relying on these systems for their life history are also experiencing widespread population declines (Brennan andKuvlesky 2005, Askins et al. 2007). This decline has led to multiple efforts in North America to conserve and restore grasslands and the associated wildlife species. For example, in prairie Canada, conservation partners in the Prairie Habitat Joint Venture (PHJV; http://www.phjv.ca/) and Saskatchewan Prairie Conservation Action Plan (SK PCAP; http://www.pcap-sk.org/) have been delivering grassland conservation and restoration, and conducting conservation research, since the late 1980s.
Negative population trends are of concern for many avian species within grassland systems (Peterjohn andSauer 1999, Brennan andKuvlesky 2005). Declines in both abundance and diversity (i.e., species richness) of grassland birds are greater than for other avian subgroups based on breeding habitat affinity (e.g., woodland, wetland, shrubland; Schipper et al. 2016). The general pattern of declines in grassland birds seems to continue despite much of the conversion of native grasslands to agriculture having been completed by the 1940s, particularly in the United States (Waisanen and Bliss 2002). Continued loss of native grasslands, however, is well documented on the Canadian prairies (Watmough and Schmoll 2007) and the compounding effects of habitat loss, fragmentation, declining patch size, industrial disturbance, and the potential impacts of climate change still threaten grassland species through direct and indirect pathways (Davis 2004, Ludlow et al. 2015, Jarzyna et al. 2016. Thus, there remains a need to identify important habitats for grassland birds in Canada.  Identifying important habitats to species of interest is a critical step in guiding wildlife conservation efforts (Thogmartin et al. 2006, Forcey et al. 2007). Over broad geographic ranges, this goal is commonly achieved by developing species distribution models (Guisan et al. 2013). These models can then be used to guide management within the extent of the models to target conservation action for optimal outcomes (Thogmartin et al. 2014).
Publicly available data have been used successfully to develop species distribution models (SDMs) and elucidate relationships between habitat and the likelihood of species presence (Seoane et al. 2004, Brotons et al. 2007). In North America, the Breeding Bird Survey (BBS) has documented species presence and abundance along set routes during the breeding season since 1966. These data have been used to develop species distribution models for many species (Peterson and Robins 2003, Matthews et al. 2011, McCarthy et al. 2012, Barbet-Massin and Jetz 2014, Steen et al. 2014). These and other similar studies have made important contributions to avian ecology, and the SDMs developed have the potential to make substantial contributions to spatial conservation decision making (Guisan et al. 2013).
We developed SDMs for a suite of grassland bird species using BBS data from the Prairie Pothole Region of Canada (PPR; Fig.  1). We had three main objectives. First, we developed SDMs using resource selection functions (RSFs) for 15 grassland bird species. Second, we assessed model performance and examined which factors influenced model performance. Third, we quantified niche overlap for a suite of species by combining the SDMs for those species with good model validation to identify multispecies priority areas across the region. Our work provides a critical foundation for landscape-planning initiatives that incorporate grassland bird conservation goals into ongoing conservation efforts in the PPR.

Study area
Our study area comprises the Canadian portion of the Prairie Potholes Bird Conservation Region 11 (BCR 11; Environment Canada 2013). We restricted our analysis to within Canada because of a lack of comparable wetland data spanning the border and to avoid duplication with previous and ongoing U.S. efforts (Drum et al. 2015). The Canadian region covers approximately 467,000 km² of the northernmost portion of North America's Great Plains and encompasses the southern portions of three Canadian provinces (southwestern Manitoba, southern Saskatchewan, and southern Alberta; Fig. 1). The region derives its name from millions of shallow depressional "pothole" wetlands that formed as subterranean masses of ice, which melted following the retreat of glaciers at the end of the last ice-age (Doherty et al. 2017).
The region is characterized by a generally dry climate supporting temperate grasslands. Remaining native grasslands are primarily fescue prairie in the western extents and mixed-grass prairie centrally. Historic tall grass prairie in the eastern extents has been all but extirpated. The northern extents of the region represent a transition from grassland to boreal forest and are interspersed with stands of aspen (Populous tremuloides) and other woody species (Doherty et al. 2017).
The region is primarily agricultural; land uses included cropland, predominantly for cereal grain and oil-seed production, and introduced and native grass forage lands (pasture and haylands) for cattle production (Ecological Stratification Working Group 1995). Loss of native grassland to cultivation remains a primary challenge for bird conservation in this region. Less than 43, 20, and 1% of native grasslands remain in Alberta, Saskatchewan, and Manitoba, respectively (Samson and Knopf 1994).

Study species
We developed species distribution models for 15 bird species in the PPR that use grassland habitats (Table 1). We focused on avian species having an affinity for grassland habitats during the breeding season. We initially considered a broader array of potential species; however, the list was rarefied based on data availability for each species of interest as described below. Six species of interest have limited range edges within the PPR. For these six species (Table 1) we limited all analyses to within their breeding range extent as defined by BirdLife International and NatureServe (2014). The BBS is an annual North American roadside survey that has been conducted for over 40 years. Observers follow set routes that are 39.4 km and conduct 3-minute point counts at each of 50 stops (i.e., one stop every 0.8 km) following standardized protocols. These data are used extensively to assess population trends (e.g., Amano et al. 2012, Bled et al. 2013, Schipper et al. 2016) and for the development of species distribution models (e.g., Barbet-Massin and Jetz 2014, Goetz et al. 2014, Steen et al. 2014. We obtained Breeding Bird Survey route and stop locations by request from the Canadian Wildlife Service (CWS). These spatial data included route beginning and end locations and the locations of all 50 stops along most routes ( Fig. 1). However, some routes contained locations for only the start and end locations ( Fig. 1). We selected all stops within the PPR boundary (or within the species range), regardless of whether the entire route was contained within the boundary. Provincial count data (50-stop) were obtained from the United States Geological Survey (https://www.pwrc.usgs. gov/bbs/ downloaded on 15.06.12). Both datasets were filtered to include only those routes classified as "Active" and "run type 1." We retained counts from 1997-2014 and summed the number of birds counted for all focal species at each stop across those years. For example, if a single bird were detected in each of five years, then the count at the stop would sum to five. Count data were georeferenced with the spatial locations by a unique combination of province number (i.e., "statenum"), route identifier, and stop identifier (e.g., 1-50).
We developed datasets for each species by identifying all stops in which a species was detected on ≥ 1 occasion between 1997 and 2014 in the count data set and assigned these locations a value of 1 (Table 1). We assigned all stops within the PPR in which the target species was not identified a value of zero for counts. We then calculated the proportion of stops at which a target species was detected over the time series and across the PPR. A species needed to be detected at ≥ 5% of stops to be considered for the development of species distribution models.

Covariate development
We developed spatial habitat covariates reflecting land-cover and wetland abundance with demonstrated or hypothesized habitat associations for each species of interest. To facilitate comparison across species, we considered the same suite of spatial covariates (n = 9; Table 2) for each species. Each covariate was summarized in neighborhoods with radii of 400 m, 800 m, 1600 m, and 3200 m around each stop for a total of 4 scales. These scales were chosen to represent multiple levels of habitat selection and for consistency with previous literature examining habitat selection in grassland birds (Bakker et al. 2002, Cunningham and. We used two geospatial datasets to develop our covariates. The first was used to generate seven variables to capture land cover at the BBS survey stop level. These variables were estimated from the 30-m resolution Annual Crop Inventory digital data layer produced by Agriculture and Agri-Food Canada (AAFC; AAFC 2014). Variables included all grassland classes combined (grass), all annual crops (crop), putative native grass (n.grass), seeded grassland for hay and pasture (t.grass), deciduous, mixed and coniferous trees (tree), and shrubs (shrub) and trees combined (tree.shrub; Table 2). These variables were static to 2014 and assumed that summarization in that year was generally representative of land-cover composition during the period of study.
The final two variables described the wetland area (w.area) and number of wetland basins (w.basins) in surrounding 41 km² landscapes (Table 2). These variables were extracted from a Ducks Unlimited Canada (DUC) data layer derived from publicly available CanVec hydrography data created from best available sources ranging in scale from 1:10,000-1:50,000 (Natural Resources Canada 2011). Because the CanVec wetland layer is known to vary in wetland capture (e.g., missing small wetland basins), the DUC database adjusts wetland area and basin count using a spatial statistical model built from overlapping DUC wetland inventory data (digitized at imagery resolution 0.5 m-2.5 m), CanVec data, and Soil Landscapes of Canada landscape covariates (Soil Landscapes of Canada Working Group 2011). The DUC data layer thus estimates wetland area and wetland basin count by applying adjustment factors to CanVec data across the Canadian PPR.

Model development and evaluation
We developed resource selection function (RSF; Manly et al. 2002) models using logistic regression (Hosmer and Lemeshow 2000) to characterize habitat selection and develop species distribution predictions for each species. We selected bestapproximating models using Akaike's Information Criterion (AIC; Burnham and Anderson 2002). Each variable was summarized at each of the four scales (i.e., neighborhood sizes). We identified the best scale for each variable by comparing single variable models at each scale and an intercept-only model. Thus, there were five univariate models compared for each variable by species comparison. The scale with the lowest AIC was selected and carried forward into the complete candidate set of models.
If the variable was not more predictive than the intercept-only model, the variable was excluded from all subsequent analyses.
Our goal was to integrate variables at the scale at which they make the highest contribution to the explained variance in the univariate model, which generally results in better model performance (Graf et al. 2005).
After excluding covariates with low explanatory power, we assessed correlation (Pearson's r ≥ |0.65|) among all potential predictor variables to avoid issues of multicollinearity. When variables were highly correlated, we chose the most predictive single variable based on AIC to carry forward into the multivariable models. We subsequently ran all possible model combinations of the remaining variables and selected the top model (i.e., the model with the lowest AIC). This model was mapped spatially using the model coefficients and covariates to create a raster of predicted probability of occurrence across the species range within the PPR. All variables were standardized (zscore transformation) prior to analysis and models were thus applied to standardized raster surfaces.
We evaluated model performance using a k-fold cross validation technique developed by Boyce et al. (2002). We used five folds of the data and iteratively fit the top model to each set of training folds (i.e., 80% of the observations) and calculated the areaadjusted observed number of observations (i.e., remaining test 20% of observations) falling into five binned RSF classes. We calculated the Spearman rank correlation between the RSF score and the area-adjusted frequency of validation points for each of the five folds and the mean area-adjusted frequency across folds.
High correlation values between RSF scores and area-adjusted frequencies suggest a model that is good at predicting occurrence of the focal species ). We assessed model performance based on the average and range of the Spearman correlation values. To examine the possibility of spatial autocorrelation, we examined residuals of the top model for spatial patterns.
To identify the most important habitat (or regions) and aid in conservation planning, we amalgamated predicted surfaces across species using several approaches. First, we summarized the average estimated coefficient value across species for each covariate and all scales. Next, we summed the number of species in which the top model contained each covariate and scale. This allowed us to assess patterns in the relative frequency of covariates and scales across all species. Finally, we produced two maps, one of which displayed the average RSF value (0-1) for each pixel across the PPR. The second map was more restrictive because it summarized the best habitat across all species. To generate this figure, we used the predicted RSF values (i.e., rather than the equal-interval bins) and selected the top 25% of pixels on the landscape for each species thus creating a binary surface representing whether a pixel was contained in the top 25%. We then summed the number of species represented in each pixel using the top 25% surface for each species. All analyses were conducted using R statistical software ver. 3.3.0 (R Core Team 2016).

RESULTS
We considered 9118 point-count locations in our analyses for all species that had a range commensurate with the PPR boundaries. Sample sizes for the six species with limited ranges were smaller and are presented in Table 1. Observations at point-count locations varied widely by species from our minimum of 5% to > 70% for the most common species, i.e., Clay-colored Sparrow (Spizella pallida), Savannah Sparrow (Passerculus sandwichensis), and Vesper Sparrow (Pooecetes gramineus; Table 1). As expected, the structure of the top models varied across species (Table 3).
Model performance based on Spearman rank correlation values between the area-adjusted frequency of validation points and RSF bins across the five folds varied widely among species (Table  3). The best models demonstrated excellent predictive capabilities of the relative probability of occurrence and the worst demonstrated negative Spearman correlations. We did not summarize and interpret models with poor validation. Therefore, we excluded models for species with an average Spearman correlation of < 0.5 or a range of correlation values > 0.5. This resulted in the retention of 10 species (Fig. 2). All subsequent result summaries and synthesis only include those species in which the model validation met these criteria.
The global model for each species included the top covariates at the scales carried forward from the univariate comparisons. The total of all possible model combinations for each species was 512. The top models for species included 3-6 covariates ( Table 3).
All variables in the top models were summarized at the selected scale and used to generate the final predicted surfaces (Appendix 1). There were no obvious spatial patterns in the residuals and we interpreted this as a qualitative indication of an absence of spatial autocorrelation in the model-predicted values (Appendix 1). Patterns in the standardized coefficients for the 10 species with good model validation suggested generally positive or neutral (i.e., β = 0) coefficients for most variables included across scales (Fig. 3). Across species, the most common variables included in the top models were wetland area (w.area; n = 5 species), wetland basin (w. basins; n = 8 species), tree shrub (n = 7), and the grass variables (t. grass n = 6, grass n = 5, n.grass n = 4). Generally, model results suggested that birds strongly avoided areas of high tree and treeshrub cover (Fig. 3). Positive selection was estimated for grass, native grass, and tame grass covariates (Fig. 3). A slight negative association was found for the wetland area count when it was included in the top model for a species (n = 5 species; Table 3). Most covariates included in the top models occurred at the larger spatial scales (Fig. 3). Unsurprisingly, one of the three grass variables (grass, n.grass, t.grass) was included in the top model for all species. The other most common covariates included in the top models were wetland basin count (n = 8 species), tree-shrub (n = 7 species), and the wetland area count (n = 5 species; Table 3; Fig. 3). The meanbinned value (Fig. 4) and top-habitat predicted surfaces (Fig. 5) both identified large areas of priority habitat for the suite of grassland birds in the Prairie Pothole Region. Figure 5 depicts a more restrictive approach to identifying priority habitats and implicates less of the region as priority habitat for the 10 species of grassland birds considered here.

DISCUSSION
We developed SDMs for 15 avian species having an affinity for grassland habitats during the breeding season using data from the BBS (1997)(1998)(1999)(2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014). As expected, models for these species generally depicted positive associations with grass habitat, negative associations with tree habitat, and mixed associations with wetland habitat, depending on how wetland was depicted. Model performance, assessed using k-fold cross validation, varied widely among the 15 species. Models with poor validation were typically characterized by data sets with a low proportion of detections across sites. However, a low proportion of detections (i.e., data availability) did not explain all of the variation in model performance. Four species with high average Spearman rank correlations had a proportion of detections < 0.20 (Baird's Sparrow, Ammodramus bairdii; Chestnut-collared Longspur, Calcarius ornatus; Grasshopper Sparrow, Ammodramus savannarum; and Sprague's Pipit, Anthus spragueii). These species all shared the common characteristic of limited distributions concentrated in large areas of grassland habitats in Alberta (Appendix 1). Because we used both detection and nondetection data, it was possible for a species to be detected (i.e., 1) and not detected (i.e., 0) on the same BBS route. Compared to other species, there were limited occurrences of the routes in which the species were both present and undetected over the 1997-2014 time series for these four species. It is interesting to note that for these species, data are absent (i.e., no BBS routes) throughout much of the best-predicted habitat in Alberta. Significant patterns emerged from the predicted surfaces associated with the six species having both low detections and poor model performance (Bobolink, Dolichonyx oryzivorus; Common Yellowthroat, Geothlypis trichas; Grasshopper Sparrow; Lark Bunting, Calamospiza melanocorys; LeConte's Sparrow, Ammodramus leconteii; and Sedge Wren, Cistothorus platensis). The three species with the most limited distributions (Bobolink, Sedge Wren, and LeConte's Sparrow) did not produce well-validated models likely because of their limited extent and data availability. For example, when producing the Spearman's rank correlations based on the Boyce index, very few of the use locations for these species fell in the top two bins. This occurred because the distribution of RSF values across the study extent was strongly skewed toward low values and resulted in the negative Spearman rank correlations estimated for each species. More focused studies that limit the spatial extent and do not predict across such large areas of mostly unsuitable habitat for these species could result in better performing models. For example, of the six species in which we limited model development to within their restricted ranges in the PPR, two species (Chestnut-collared Longspur, Grasshopper Sparrow) had better model performance when the models were developed across the restricted range (B. Fedy, unpublished data).
The coefficient estimates for w.area and w.basin were modest, but consistent in their directionality. When included in the top model, w.area was consistently negative whereas w.basins was largely positive. This difference is likely due to the different landscape patterns captured by the two covariates. Higher values and the positive association with the w.basins covariate typically represent landscapes with many smaller wetlands and a more heterogeneous landscape. This heterogeneity in the context of wetland areas, could be more beneficial to grassland birds than landscapes containing large expansive wetlands as represented in the w.area covariate (Skinner and Clark 2008). Additionally, smaller wetlands tend to be more ephemeral and may provide ideal conditions for certain species over the breeding season (e.g., Baird's Sparrow). The strongest coefficients (when included in the top models) were positive associations with native grass (n.grass) and grass and avoidance of tree-shrub and trees, particularly at the larger scales. These relationships are in the expected direction and suggest that our models are capturing ecologically relevant variation in habitat preferences across the study extent.
Larger scales (e.g., 1600 m and 3200 m radii) were consistently selected for all variables with the exception of crop and tame grass variables. Thus, using these data on species presence and absence, landscape context (i.e., configuration at larger spatial scales) is more important to the relative probability of use than smaller, patch-scale, metrics of habitat composition as estimated from our spatial layers. The scales used in our study ranged from a 400 m to 3200 m radius and thus encompassed areas surrounding locations from 0.5 km² to 32.2 km². These findings suggest that larger scales are consistently more predictive than smaller scales for the 10 grassland bird species modeled here. Shahan et al. (2017) also documented that landscape variables were more predictive than local variables for eight grassland songbird species and suggested the importance of considering the landscape context of prairie fragments (Shahan et al. 2017). Our findings regarding the importance of landscape scales are contradictory to the conclusions of some previous studies of grassland birds. For example, Koper and Schmiegelow (2006) and Davis et al. (2013) highlighted the importance of and argued that local-scale (e.g., 400 m radius) summaries of environmental covariates had greater predictive capabilities than larger scale summaries of landscape composition. There are several important differences between our study and previous efforts, however. First, multiple previous studies did not explicitly address the predictive capabilities of landscape composition summarized across more than two window sizes (Brotons et al. 2005, Koper and Schmiegelow 2006, Davis et al. 2013). Additionally, two earlier studies highlighting smaller scales included a covariate for "habitat type" in the models. Though it is not entirely clear from the manuscripts, the habitat type classification likely summarized habitat composition at relatively large spatial scales and this covariate was common in the top-ranked models for many species (Ribic et al. 2009, Davis et al. 2013. For Sprague's Pipit, the selection of grassland variables at 1600 m radius (8.01 km) was similar to the scale of a grassland aggregate covariate (10.4 km²) that was the most predictive of Sprague's Pipit occurrence in a previous study developing SDMs using different data and modeling approaches than ours (Lipsey et al. 2015). Habitat selection is clearly a hierarchical process and there is not one scale that will capture all relevant information for a species (Johnson 1980). For grassland species in particular, Cunningham and Johnson (2006) found that combining proximate and landscapelevel predictors resulted in the best-performing models for multiple grassland bird species. Johnson and Igl (2001) (Johnson and Igl 2001). These results are similar to Davis (2004) who found area sensitivity in Baird's Sparrows, Chestnut-collared Longspurs, Grasshopper Sparrows, and Sprague's Pipits and the inverse relationship (i.e., area insensitivity) for Clay-colored Sparrows and Western Meadowlarks. Our study did not explicitly address area sensitivity, but we documented positive association with grass variables at the two largest scales (1600 m and 3200 m radii) for Baird's Sparrow, Brown-headed Cowbird, Chestnut-collared Longspur, Grasshopper Sparrow, and Western Meadowlark, whereas Savannah Sparrow and Vesper Sparrow were the only species for which the smallest scale (400 m) for grass or cropland was selected. For Baird's Sparrow the strongest covariate was for the avoidance of trees at the largest spatial scale. Generally, the consistently better performance of the larger scales may suggest some form of area sensitivity. As highlighted by Johnson and Igl (2001), replication in space is important and results from one area may not apply to others because of a suite of factors including study design, analytical methods, location relative to range of the species, and surrounding landscapes. Our results are consistent, though, with other studies on Baird's Sparrow that have found higher occurrence in grassland habitat (McMaster and Davis 2001) and avoidance of cultivation (Owens and Myres 1973), shrubs, and high visual obstruction vegetation (Madden et al. 2000). The questions surrounding area sensitivity and requirements in grassland birds are complex and likely somewhat location and species specific. For example, With et al. (2008) concluded that even relatively large grassland areas may be insufficient for the persistence of some grassland species in the Flint Hills region of Kansas and Oklahoma. Fisher and Davis (2010) reviewed 57 studies of grassland bird habitat relationships and identified 9 of the most important variables based on previous research. As noted by the authors and others, identification of remaining grassland habitat is a critical first step to population conservation and recovery Herkert 2001, Fisher andDavis 2010). Our development of SDMs for multiple grassland bird species represents a significant step toward that goal. The challenge with the incorporation of the nine variables identified by Fisher and Davis (2010) into spatially explicit models is that most of the covariates are not currently available in GIS (e.g., bare ground, forbs, litter) and several variables will likely prove too difficult to ever model efficiently in a GIS environment (e.g., vegetation height, vegetation volume, litter depth, dead vegetation). We suspect that other variables that were measured in previous studies, but not at the landscape scale addressed here (e.g., shrub and tree cover), will prove important in future efforts (Thompson et al. 2014).
Aside from our crop variable, we did not directly include other anthropogenic variables that could affect habitat use such as energy development (Ludlow et al. 2015) or roads (Sutter et al. 2000). Multiple studies have documented the potential negative effects of both conventional and renewable energy development on grassland birds. For example, Shaffer and Buhl (2016) documented displacement of several species (Grasshopper Sparrow, Western Meadowlark, Bobolink, Savannah Sparrow, Vesper Sparrow, Chestnut-collared Longspur, and Clay-colored Sparrow) in relation to wind turbines. Thompson et al. (2015) documented avoidance of habitats well beyond the area occupied by the infrastructure for many grassland species in the Bakken. However, as in other studies (Ludlow et al. 2015), the presence, strength, and mechanisms of negative effects were not consistent across all species. For example, Ludlow et al. (2015) found negative influences on density and reproductive success for Sprague's Pipit and Baird's Sparrow, but not for Vesper Sparrow, Western Meadowlark, or Brown-headed Cowbird, with more Brown-headed Cowbirds occurring in close proximity to energy infrastructures. In a report prepared for the Petroleum Technology Alliance Canada, Linnen (2008) examined the impacts of energy development on most of our 10 species (with the exception of Clay-colored Sparrow) and concluded that Chestnut-collared Sparrows, Sprague's Pipits, and Baird's Sparrows were the most sensitive to gas development and edge avoidance whereas Vesper Sparrows were more abundant near traditional oil and gas development. Species distribution models can be combined with information on energy development potential to help guide the siting of infrastructure (Tack and Fedy 2015), identify potential areas for mitigation (Kreitler et al. 2015), and explore the impact of various development scenarios (McGowan et al. 2013, Lipsey et al. 2015. We suggest that our models could be used in similar applications, particularly for those species expected to respond negatively to energy development (e.g., Baird's Sparrow).
As with most species distribution models, it is important to recognize that predictions are static and assume the relationships between species and environment are at equilibrium (Araújo and Pearson 2005). In light of recent global changes in climate and habitats, species-environment relationships may become more transient (Yackulic et al. 2015). Our models represent a baseline for the development of more sophisticated dynamic occupancy models that could include time varying conditions for both habitat and climate. Future updates of our models with new data on both bird occurrences and temporally varying covariates could prove illuminative regarding temporal variation in the relative probability of habitat use (Gorzo et al. 2016) and in particular, the potential impacts of climate change (Steen et al. 2014, Conrey et al. 2016).
One of our goals was to identify important habitat as a first step toward informing spatial conservation decisions regarding grassland birds within the PPR in Canada. Our models can be incorporated on a species by species basis or through the use of our two summary surfaces that were intended to represent the best habitat across all 10 species. We recommend future research combine our results with other sources of information, particularly those related to bird abundance and stressors affecting abundance, to help define critical habitat. The integration of SDMs and spatially explicit estimates of abundance can be a valuable tool for the identification of priority habitats for species of conservation concern (Thogmartin et al. 2014, Doherty et al. 2016. Additionally, these surfaces could help inform priority areas for reserve design when combined with other relevant data (Moilanen et al. 2005). Ultimately, SDMs can be a valuable tool within structured and transparent decision-making processes (Guisan et al. 2013).
Responses to this article can be read online at: http://www.ace-eco.org/issues/responses.php/1143