One step ahead of the plow: Using cropland conversion risk to guide Sprague's Pipit conservation in the northern Great Plains

a Wildlife Biology Program, University of Montana, 32 Campus Drive, Missoula, MT 59812 USA b United States Fish and Wildlife Service, 134 Union Blvd., Lakewood, CO 80228, USA c United States Fish and Wildlife Service, Region 6 HAPET Office, 922 Bootlegger Trail, Great Falls, MT 59404, USA d The Nature Conservancy, 117 E Mountain Ave Ste 201, Fort Collins, CO 80524, USA e Department of Zoology and Physiology, University of Wyoming, 1000 E. University Ave, Laramie, WY 82071, USA f Environment Canada — Canadian Wildlife Service, 300, 2365 Albert Street, Regina, SK S4P 4 K1, Canada g Natural Resources Institute, University of Manitoba, 70 Dysart Rd. Winnipeg, MB R3T 2 M6, Canada


Introduction
Grasslands are among the most imperiled ecosystems worldwide (Hoekstra et al., 2004) because their soils provide some of the most productive farmland on earth. As rising global food demand surpasses improvements in yields on existing cropland, additional grassland conversion will be required to feed a projected 11 billion people by 2050 (Foley et al., 2011;Ray et al., 2013). Rising commodity prices exacerbated by demand for biofuels threatens to further expand cropland agriculture (Fargione et al., 2009;Wright and Wimberly, 2013). In temperate North America, historic grassland losses total approximately 70%, including complete conversion of the most productive areas where nothing but remnant tracts persist (Samson et al., 2004). In the northern Great Plains where most grassland remains, agricultural conversion is happening five times faster than grasslands can be protected Walker et al., 2013).
A steep and consistent decline in songbird populations reflects eroding ecosystem integrity in North American grasslands (Brennan and Kuvlesky, 2005;Sauer et al., 2014). Of high concern is Sprague's Pipit (Anthus spragueii; herein "pipit"), a grassland obligate species that breeds in the native mixed prairie of Saskatchewan, Alberta, Montana, and the Dakotas . The pipit that has been declining N3% annually across North America since 1966 (Sauer et al., 2014), is listed as globally vulnerable by the International Union for Conservation of Nature (IUCN, 2014), is federally threatened in Canada (Environment Canada, 2012) and is being considered for federal protection under the U.S. Endangered Species Act (ESA, 1973;USFWS, 2010). The ESA status assessment focuses attention on pipits and underscores the urgency for conservation of northern grasslands.
Broad-scale planning enables systematic targeting of scarce conservation resources (Bottrill et al., 2008), and sensitive species provide a useful lens for delineating landscapes of high conservation value and Biological Conservation 191 (2015) [739][740][741][742][743][744][745][746][747][748][749] identifying impacts of human activity (Sanderson et al., 2002). Spatially explicit tools enable practitioners to target implementation where populations will benefit most (Margules and Pressey, 2000). We present a three-part analytical framework for strategic conservation of pipits in North America. First, we depict a range-wide distribution model by integrating survey efforts across a million-km 2 area of Canada and the United States. Using our model, we describe factors shaping pipits' continental distribution and delineate core areas of high bird abundance. Second, we assess vulnerability to future habitat loss using soil capability for agriculture as an index of conversion risk. For the U.S. portion of the range, we employ a quantitative risk model to develop future scenarios of cropland expansion and assess potential impacts to populations. Finally, we explore the relationship between land tenure and population distribution to guide appropriate conservation approaches.

Study area
Our study area includes the intersection of the Breeding Bird Survey range for pipits (Sauer et al., 2014) and the Plains and Prairie Pothole Landscape Conservation Cooperative (PPPLCC; Millard et al., 2012), a consortium of public and private conservation partners (Fig. A1). The region covers portions of Alberta, Saskatchewan, Montana and the Dakotas. This area is made up of diverse mixed grass prairie with level to rolling terrain. It encompasses interspersed badlands and sagebrush steppe in the west and pothole wetlands and prairie parklands in the east and north. It includes portions of Great Plains-Palouse Dry Steppe (331), Great Plains Steppe (332) and Prairie Parkland (251) provinces as described by Bailey (1995).

Bird survey data
Range-wide perspectives are required for conservation of migratory songbirds. Data limitations and inconsistent collection methods have hindered efforts to model bird distributions across broad scales. To describe the distribution of pipits in their breeding range, we combined data from 76,623 point counts (2007-2012; Table A1) into one integrated analysis. Integration allowed us to achieve spatial coverage that made a continental perspective possible. To assess the influence of point count distance and duration on pipit detectability, we conducted a sensitivity analysis in the heart of the range (northeast Montana). The test dataset, collected in 2012-2013, contained known distance and time intervals for evaluation (author ML, unpublished data). We truncated data by 1-min time interval (0-1, 0-2, 0-3, etc.) and distance interval estimated to the nearest 10-m. We then used linear models to estimate the effect of time and distance on observed detection probability and abundance. Detection probability remained relatively insensitive to count duration (1% increase per minute) and distance (4% per 100-m). Bird abundance was more strongly affected (3 and 8% respectively) so we limited subsequent modeling to presence/absence data.
We removed repeated and overlapping records, keeping the most recent surveys within 200-m of one another based on average pointcount radius. Surveys were not targeted for pipits and data were heavily skewed towards absence, so we stratified records to ensure appropriate class balance with 40% occurrence. Because survey locations were highly clumped in some regions, we thinned the dataset to 10,000 records (approximately 30%) using a random sampling algorithm weighted by the inverse proportional kernel density estimate of sampling intensity using the R package spatialEco (Evans, 2015). Thinning directly addresses sample independence by modeling the spatial process of the observed data and effectively smoothing the sample distribution to a more uniform spatial process. However, some portions of the study region, including the Dakotas and Saskatchewan, had lower data availability and contributed proportionally fewer records even after thinning (Fig. A1).

Environmental predictors
Climate has a strong relationship with bird habitat in North America and long-term averages reflect envelopes that shape geographic ranges (Jiménez-Valverde et al., 2011;Thomas, 2010). We obtained long-term averages for North America from USDA (2012). Developed from monthly climatic surfaces spanning , this dataset uses thin-plate splines to depict continuous variables over a stable climatic period that likely reflects establishment patterns for plant communities (McKenney et al., 2001;Rehfeldt, 2006). We included five long-term averages in model selection: mean annual precipitation (mm), mean annual temperature (C°/10), total growing season precipitation (mm), summer precipitation balance and average frost free period in days. Climate variables were highly correlated, so we chose those most relevant to herbaceous vegetation growth and that had correlations ≤0.8.
Because grassland birds are sensitive to vegetation structure (Fisher and Davis, 2010), we included three shorter-term measures of vegetation growth and moisture, including Gross Primary Productivity (GPP), maximum annual snowfall, and the Palmer Drought Severity Index (PDSI). We averaged short-term measures across 2006-2010, including the year preceding bird surveys because residual vegetation is an important component of grassland bird habitat (Ahlering et al., 2009). Comparable data for 2011-2012 were not available at the time of analysis. GPP provides an index of the amount of vegetation growth and is derived from Moderate Resolution Imaging Spectrometer satellite imagery at 8-day intervals (Reeves et al., 2006). We represented GPP as the maximum measurement during April-July using values obtained from NASA (2012). Maximum snow depth for winter between October and April was obtained from Snow Data Assimilation System (National Operational Hydrologic Remote Sensing Center, 2004). For PDSI we used global 2.5°gridded monthly data for May self-calibrated with the Penman-Monteith potential evapotranspiration formulation, 1900-2010 (Dai and NCAR, 2014).
We included attributes for land cover because pipits require relatively large and intact grassland on their breeding grounds (Davis et al., 2006;Sliwinski and Koper, 2012). We used four variables that describe patterns in land-cover: proportion of cropland, forest, grassland and grassland aggregation index. We derived 400-m resolution binary layers of crop, forest and grassland from 30-m products created by Agriculture Agri-Food Canada in 2010(2015 and level II of the United States Geological Survey 2011 National Land Cover Dataset (Homer et al., 2015). Land use classes had overall accuracies of 78-79% in the U.S. (Wickham et al., 2010(Wickham et al., , 2013 and about 93% in Canada (Agriculture and Agri-Food Canada, 2015). For each 400-m cell, we estimated the proportion land cover and aggregation in the surrounding neighborhood using a moving window average. We used the program FRAGSTATS (McGarigal et al., 2012) to estimate grassland aggregation index, calculated as the proportion of within-class adjacencies among neighboring pixels and indicating degree of fragmentation. We tested a range of window sizes (0.026 to 100-km 2 ) with binomial regression and used Akaike's Information Criterion (AIC) to determine the most predictive scale (Fig. A2). Grassland cover was most predictive 1-km 2 around surveys (ΔAIC = 3-400) and aggregation index was most predictive at 10.4-km 2 (ΔAIC = 3-224). Grassland cover was closely related to aggregation index (r = 0.88) but correlation coefficients among other land cover variables were moderate (−0.78 to −0.13).

Statistical methods
We specified a binominal model with a probabilistic outcome using the nonparametric model Random Forest (Breiman, 2001) in program R (Liaw and Wiener, 2002;R Development Core Team, 2013). Random Forest is a bootstrapped Classification and Regression Tree (CART) approach that is based on the principle of weak learning (Hastie et al., 2008), where a set of weak subsample models converge on a stable global model. This method has been shown to provide stable estimates while being robust to many of the issues associated with spatial data (e.g., autocorrelation, nonstationarity). It fits complex, nonlinear relationships and accounts for high dimensional interactions (Cutler et al., 2007;Evans et al., 2011). We assessed competing models by comparing model importance, calculated as smallest out-of-bag (OOB) error, smallest maximum within-class error and fewest parameters (Murphy et al., 2010) using package rfUtilities (Evans and Murphy, 2014). Parsimony in Random Forests reduces noise, produces a more interpretable model and results in better model fit (Evans et al., 2011;Murphy et al., 2010). To assess model fit, we calculated OOB error, root mean square error (RMSE; Willmott, 1981) and the calculated area under the receiver operating curve (AUC; Metz, 1978) in two ways: first, with data used in model fitting, and second, with data withheld from fitting after thinning (see Section 2.2.1). To assess model performance, we conducted independent cross-validation using 1000 replicated predictions with 10% withheld data from each replicate. We quantified error as the cumulative error rate across bootstrap replicates.
Imperfect detections (false absences) in bird survey data can be a significant source of error in models of distribution. Available estimates of detection probability for pipits range between 0.70-0.82 (S. Davis [0.70] and M. Lipsey [0.82], unpublished data). We assessed the effect of false absences on model fit by conducting a sensitivity analysis in the package rfUtilities by randomly changing a proportion (p = 0.05) of presences to absences and running a series of perturbed models. This test also reflects model sensitivity to lack of independence in predictors. We observed a small standard deviation (δ = 0.0005), across n = 999 simulations, indicating model stability. This can be partially attributed to the ability of Random Forests to predict through noise and is an advantage of weak learners (Breiman, 2001).

Population core areas
To estimate the regional distribution of populations, we first resampled the model prediction raster from the arbitrary resolution of environmental layers (400 × 400 m or 16 ha) to units that approximated territory size for male pipits (160 × 160 m or about 2.6 ha; Fisher and Davis, 2011) using bilinear interpolation in ArcGIS (ESRI, 2010). We summed the probability of occurrence across all pixels in the study region to generate an index of total population. We placed each grid cell prediction in the context of the study area by dividing the individual pixel probability by the total index. Starting with the highest-value pixels, we cumulatively summed the probabilities until a given threshold was met. We set 25, 50 and 75% thresholds to delineate cores as the smallest possible areas containing the highest concentrations of predicted pipits. We estimated proportion of the population within multiple political and ownership boundaries by dividing the sum of occurrence probabilities in each class by our total population index.

Continental cropland risk
To estimate future conversion risk of grassland to cropland, we used existing soil databases to overlay soil capability for agriculture on the pipit distribution. Soil capability classes are ranked 1 to 8, with 1 being the most suitable for crops and 8 the least. We accounted for slight differences in soil classifications between the two countries by combining categories 1-2 (most arable), 3-4 (some limitations) and 5-8 (least arable). Conversion rates tracked soil capability, with the most arable land (classes 1-2) largely already converted (70%). By comparison, only 47 and 5% of the moderate (classes 3-4) and least arable soils (5-8), respectively, have already been converted. Using the species distribution model probability surface, we calculated the simple proportion of the predicted population on untilled land in each soil class. We also calculated the proportion of land and population that are legally protected from agricultural conversion within each soil class. We obtained soil data from the Natural Resources Conservation Service web soil survey database (NRCS, 2014) in the U.S. and from the Canada Land Inventory in Canada (1998).

U.S. cropland risk scenarios
To identify regions and populations at risk from conversion in the U.S., for which we had more detailed data than we did for Canada, we used a cropland suitability model described in Smith et al. (in preparation). The model provides a probability surface with values from 0-1, representing the relative suitability of each grid cell for conversion to cropland. We used this surface to develop three potential build-out scenarios, a-c. In each scenario, land above a given probability cut-point was assumed to be converted. Pattern and rate of future grassland loss is difficult to predict; therefore, we use scenarios only as reference points for planning. Scenarios do not reflect variation in rates of conversion, nor do they refer to a given time horizon. They represent the spectrum of plausible absolute losses in grassland area due to cropland expansion based upon observed rates of loss ( Fig. A3; Doherty et al., 2013;GAO, 2007).
Scenario (a) represents minimal, or background conversion, scenario (b) represents a constrained growth scenario for cropland and scenario (c) represents unconstrained cropland growth. We reclassified the probability surface raster to produce predicted conversion layers. Probability cut-points that defined scenarios were selected as (a) 0.98, (b) 0.7 and (c) 0.3 (Table 2) after visual inspection of the area accumulation curve derived from the tillage model (Fig. A3). For each scenario, we removed pixels of predicted new cropland from the original land cover layer of grassland. All federal land and state lands were considered protected from cropland conversion except state school trust lands. Tribal lands included in the Bureau of Indian Affairs (BIA) databases were treated as private, thus not protected from conversion. Existing cropland was also excluded. We converted altered grassland layers to proportion and aggregation variables and substituted them into the original model to re-predict pipit distribution under each scenario.

Ownership
To estimate the composition of land tenure and conservation status of the pipit population in the U.S., we used ownership and protection data compiled for Doherty et al. (2013). In Canada, we built an ownership layer by combining boundaries from provincial, federal, and private conservation areas. To quantify areas protected from cropland in Canada, we obtained parcel boundaries of provincial land that were legally protected from cultivation from the Alberta, Saskatchewan, and Manitoba provincial governments and those of federal lands from Environment Canada. We also obtained information on lands that were privately owned and legally protected from cultivation from private conservation agencies (e.g., Nature Conservancy of Canada, Ducks Unlimited Canada). We considered lands with perpetual conservation easements to be protected from conversion but lands under volunteer or management agreement to be available for conversion. For each ownership class, we summed the value of the species distribution model probabilities and divided by the total to produce a proportional estimate of population density by ownership.

Species distribution model
The most supported and parsimonious model of pipit occurrence included nine predictors: proportion of grassland, grassland aggregation index, PDSI, average maximum snowfall, growing season precipitation, summer precipitation balance, average frost-free period, mean annual precipitation and mean annual temperature (Fig. 1). RMSE was 0.06 with a 14.9% OOB error rate. The AUC was 0.99 when predicting data used in model fitting and 0.88 when predicting data that were not  used (see Section 2.2.1). The model performed well based on error estimated by cross-validation (15.6%). All assessments indicated good model performance with high predictive accuracy (Fawcett, 2006).
Landscapes with a high proportion of aggregated grassland and with relatively cool, moist climates were most likely to contain pipits. Effects of environmental and climatic predictors on pipit distribution were nonlinear (Fig. 1). Strongest predictors were moisture variables (maximum snowfall, PDSI, growing season precipitation and summer precipitation balance) combined with proportion of grassland 1-km 2 around survey points.

Population core areas
Breeding pipits were unevenly distributed across their range and were concentrated in core areas characterized by grassland (Fig. 2a). The relationship between population density and area was steep, with 25% of the population within 5% of the study area and 75% of birds within 30% of the study area (Fig. 2b). Regions of highest pipit density were predicted in southeast Alberta, southwest and south-central Saskatchewan, and northeast Montana. Our model also indicated several smaller core areas in southwest Manitoba and in western portions of the Dakotas. About 60% of the population occurred in Canada, with 40% in the U.S. Alberta and Saskatchewan together contained nearly all (97%) of the Canadian population and 57% of the global population. Montana contained 63% of the U.S. population, with most of the remainder in the Dakotas (Table 1).

Continental cropland risk
Observed frequency of pipits was three times lower in cropland (13%) than across all other land uses combined (40%). Continentally, we estimate that 21% of breeding pipits occupied grasslands that are legally protected from conversion to cropland. Conversely, a quarter of the population occupied unprotected grasslands at risk of future conversion (soil capability classes 1-4; Fig. 3b). Pipits occupied protected grasslands underlain by arable soils (classes 3-4) more than expected ( Fig. 3a; χ2 test, p = 0.06). In contrast, they avoided the most arable, unprotected landscapes (classes 1-2) where widespread conversion has already impacted grasslands (χ2 test, p = 0.1). Protection from conversion was inversely related to soil capability, with grasslands on more arable soils less protected. Protection status was low with 2% of the most arable soils protected (classes 1-2), 8% of classes 3-4 and 23% of classes 5-8 protected.

U.S. cropland risk scenarios
Within the U.S., potential population declines varied from 2-27% across three conversion scenarios (Fig. 4). Our model indicated a 2% population decline with background growth (Fig. 4; scenario a), a 13% decline with constrained growth (b) and 27% decline with unconstrained growth (c; Table 2). Background rate of conversion in scenario (a) predicted few grassland losses with only the easternmost fringe of core areas affected. Under the constrained conversion scenario (b), additional core populations were at risk, particularly in smaller habitat blocks along margins of the pipit range. Scenario (b) also predicted habitat loss in the largest core area in the U.S. in northern Montana (Fig. 4). Unconstrained cropland expansion in scenario (c) resulted in habitat losses across most of the eastern portion of the range and the western margins along the Rocky Mountain Front. Intact grasslands were predicted to remain in the south and central portions of the U.S. distribution (Fig. 4).

Ownership analysis
Land tenure was heavily skewed to private ownership amidst a mosaic of federal, tribal, and state/provincial lands. Our model suggests that 70% of the global breeding population was located on lands under private ownership. We also document that state/provincial, federal and Tribal/First Nation lands each contained considerable portions of the population (Table 3).

Discussion
A broad-scale perspective can inform systematic approaches to achieving conservation with limited resources. Anchored within core areas of high abundance, our approach links populations to landscape conservation at a continental scale. In western North America, core areas are being used to guide investments for high-profile and at-risk species like woodland caribou (Schneider et al., 2010) and greater sage-grouse (Copeland et al., 2013;Doherty et al., 2010;USFWS, 2013). Often, core areas for focal species coincide with important habitat for other species of interest, as recently demonstrated for mule deer and sage-grouse in Wyoming (Copeland et al., 2014). Indeed, predicted core populations of pipits in northeast Montana and southern Saskatchewan overlapped qualitatively with important migratory corridors for pronghorn (Poor et al., 2012) and sage-grouse (Tack et al., 2011). At state and provincial scales, efficiency of conservation for pipits would be maximized by focusing initial investments in southeast Alberta, southwest and south-central Saskatchewan, and northeast Montana (Fig. 2).
Our distribution model suggests that broad-scale climate patterns strongly influence pipit habitat selection (see also George et al., 1992;Wiens et al., 2008). Climate variables, especially those related to precipitation, were highly predictive and pipits selected an envelope of moderate moisture at a continental scale and in a non-linear fashion. Whether this moisture envelope produces vegetation structure that is relatively sparse or dense depends on geographic context as well as management (Bakker et al., 2002;Madden et al., 2000). This may explain why a regionally based model developed by Niemuth et al. (2008) did not report a relationship between climate and pipit abundance. Variability captured in our range-wide approach boosted power to detect climate relationships, and the random forest approach is well suited to characterize non-linear relationships. Our continentalscale analysis explains variation in local-scale studies by capturing the range of environmental conditions that shape populations.
Pipits are known to avoid exotic and planted grasses (Davis et al., 2013), but to the best of our knowledge no continental data exists for grassland composition and we were unable to include it in this analysis. Moisture and climate variables relevant to plant communities were included to capture some variability in herbaceous plant structure and composition. Although coarse land cover data limited our model's ability to discriminate fine scale ecological processes such as vegetation composition (second order; Johnson, 1980) or inform site level decision making, it helps us to understand broad scale (first order) habitat selection.
Analysis of soil capability for agriculture demonstrates that the distribution of pipits has contracted in response to cumulative impacts of tillage in arable grasslands. This pattern is supported by spatial variability in trend estimates in the Breeding Bird Survey (Fig. A4; Sauer et al., 2014). Current evidence suggests that pipits avoid cropland (Davis et al., 1999;Owens and Myres, 1973) and future research is needed to evaluate the contribution of individuals occurring in or around cropland to population growth (Donovan and Thompson, 2001;Pulliam, 1988). Despite avoidance, low predicted densities across 56.4 million hectares of cultivated land could represent up to 29% of the breeding population.
Unconstrained cropland growth predictions suggest that risk to the population is moderate and that potential losses are comparable at U.S. (27%) and continental (25%) scales. However, tillage of an additional 12.5 million hectares in the U.S. is unlikely and losses can be mediated through proactive and targeted action. Conservation of this species depends on a shared vision for sustainable ranching as 70% of pipits rely on privately owned grasslands, often maintained as rangelands for livestock production. Moreover, another 12% of the population breeds on provincial lands in Canada that are privately managed. Public and tribal lands support less than a third of populations, though tillage risk scenarios suggest that these are continentally important for insulating against increased cropland expansion.
Results presented here are specific to effects of cropland conversion on the breeding grounds and similar analyses for winter range would inform a more holistic strategy throughout the pipit life cycle. Nextgeneration analyses could also incorporate other potential risks such as climate change (Skagen and Yackel Adams, 2012) and energy development. Effects of energy infrastructure on pipit abundance are variable (Hamilton et al., 2011;Kalyn Bogard and Davis, 2014) but linear features and establishment of exotic grasses associated with oil and gas development may reduce reproductive success (Ludlow et al., 2015). An additional 50,000 new oil and gas wells are added annually in central North America (Allred et al., 2015) and if drilling continues as anticipated, a regional understanding of potential impacts to pipits is warranted.

Conclusion: a roadmap for pipit conservation
We show that conservation of pipits depends upon a systematic approach that invests heavily in private land partnerships. Identification of priority landscapes can enhance decision-making and promote immediate and effective conservation actions. Private landowners are willing to implement beneficial practices for wildlife (Henderson et al., 2014) and the capacity to do so is growing as coordinated approaches become available (Neudecker et al., 2011). Voluntary incentives can help offset high economic returns from cropland by compensating producers for the conservation value of native grasslands. Partnerships should work to develop a portfolio of incentives relevant to the diverse needs of landowners. Some examples include conservation easements (Fishburn et al., 2009), compensation for ecological goods and services, drought mitigation and marketing of livestock products raised on native grasslands. Acceleration of pipit conservation will likely require additional coordinated funding.
Improvements in agricultural policy that incentivize ranching would also curb tillage expansion. For example, the new 'Sodsaver' provision in the 2014 Farm Bill (U.S. Agricultural Act of 2014; H.R. 2642) denies full federal insurance subsidy to recently converted cropland (Miao et al., 2014). Additional modification of subsidies could further reduce conversion of marginal land as incentives still favor farming over ranching (GAO, 2007). Higher returns from cropland also entice policy makers to lease public lands for farming where permitted. For example, prohibiting tillage on state school-trust lands would remove a primary threat to grassland cores in Montana. Other policy incentives that generate and maintain interest in grassland conservation should also be considered. One approach would be to modify the U.S. Conservation Reserve Program allowing more frequent grazing, mirroring the Permanent Cover program in Canada (McMaster and Davis, 2001).

Acknowledgments
This work was supported by the U.S. Bureau of Land Management and the U.S. Fish and Wildlife Service Plains and Prairie Potholes Landscape Conservation Cooperative, the Prairie Pothole Joint Venture and the Habitat and Population Evaluation Team. We thank the large      A2. Comparative strength of relationship between Sprague's Pipit occurrence and grassland proportion (a) or grassland aggregation (b) measured across spatial scales. Variables were calculated for rectangular neighborhoods surrounding survey locations using area indicated on the x-axis (log scale). Relationships estimated using binomial regression and compared using Akaike's Information Criterion (AIC; y-axis).