Land-use change reduces habitat suitability for supporting managed honey bee colonies in the Northern Great Plains

In 2006, biofuel crops surrounding commercial apiary locations were generally confined to far eastern portions of ND and SD (Fig. 2A). In 2014, biofuel crop area surrounding apiaries generally expanded west and northward across the study region, with continued intensification in eastern ND and SD and southern SD (Fig. 2B). Our trend analysis revealed significant annual gains in biofuel crop area around registered apiary locations from 2006 to 2014 [ ${\widehat{\beta }}_{YEAR}$ = 9.1 ha annually, 95% credible interval (CI) 8.9–9.3]. Across ND and SD, between 2006 and 2014, there were an additional 1.2 Mha of biofuel crops surrounding registered apiary locations. Much of the increase in biofuel crop area around apiaries was focused in the Prairie Pothole Region (PPR) of the Dakotas, a region extending east and north of the Missouri River in ND and SD (Fig. 3A). Average annual gains in biofuel cropping area were four times greater among registered apiaries in the PPR [ $\overline{x}$ = 10.3 ha ± 11.3 (1 SD)] than in apiaries west or south of the Missouri River, a region also known as the Badlands and Plains Region (BPR) ( $\overline{x}$ = 2.5 ha ± 5.7). There were 13,038 and 5,325 registered apiary locations in the PPR and BPR, respectively. Of the 432 apiary locations exhibiting an annual increase in biofuel crops of >30 ha, 98% were located east or north of the Missouri River, in the PPR. In general, counties with greater gains in biofuel crop area tended to have higher densities of registered apiary locations, suggesting that recent expansion of corn and soybean plantings may be encroaching into the core area of Dakota beekeepers (Fig. 3A).

The grassland trend analysis revealed a systematic decrease in grassland land cover surrounding registered apiary locations from 2006 to 2014 ( ${\widehat{\beta }}_{YEAR}$ = −0.8 ha annually, 95% CI −0.59 to −0.97). Our interpolation model of grassland change showed that apiaries with larger gains in biofuel cropping area also lost more grassland (Fig. 3B). Of the 3,105 apiary locations exhibiting a >10-ha annual decrease in grassland, 81% were located east or north of the Missouri River, in the PPR. Areas that exhibited high levels of grassland loss and high apiary density were generally confined to central and southern ND and the eastern half of SD (Fig. 3B).

Relationships among our land-cover and land-use covariates were highly varied, with Grassland and Biofuels exhibiting the strongest negative correlation (Fig. S2). All covariates included in the same model had correlation coefficients <0.3. Grassland was the most common land cover surrounding apiaries in this region, followed by biofuel crops, small grains, and open water (Fig. S3). Our COMMODITY crop model revealed that the probability of a site being used as a commercial apiary was negatively related to our commodity crop covariates (Fig. 4A). In general, Biofuels (−0.64; 95% CI −0.77 to −0.50) exhibited a stronger negative correlation with site use than Sm_Grains (−0.43; CI −0.58 to −0.28), suggesting a slightly stronger avoidance of biofuel crops than small grain fields by commercial beekeepers. Our HABITAT model estimated a strong positive relationship between apiary site use probability and grassland area (Grassland, 0.70; CI 0.56 to 0.83), alfalfa (Alfalfa, 0.25; CI 0.13 to 0.28), and open water (Water, 0.29; CI 0.17 and 0.42) (Fig. 4B). The model revealed equivocal results for associations between apiary site use and woodlands (Forest, −0.016; CI −0.45 to 0.13) and sunflower fields (Sunflower, −0.04; CI −0.18 to 0.11), with both parameters having credible intervals that overlapped zero. Results from our CONSERVATION model show that commercial beekeepers were more likely to use sites with larger areas of CRP land (CRP, 0.19; CI 0.08 to 0.31) (Fig. 4C). This model also demonstrated a weak positive relationship between other state and federal lands and apiary site selection probability (Fed_State, 0.08; CI −0.03 to 0.20); however, the credible intervals overlapped zero.

Model validation showed that all models performed better than random in predicting use of 196 sites (Fig. S4). Our HABITAT and COMMODITY models yielded similar discriminatory results, with both models having comparable area under the curve (AUC) values and correctly discerned a higher number of validation sites than our CONSERVATION model.

As global demand for resources and sustainable energy increases, there is a pressing need for a holistic examination of the impact of land-use change on a suite of ecosystem services, environmental tradeoffs, and biodiversity impacts (25, 37, 38). Here, we examined the impact of biofuel crop production on honey bee habitat; however, other impacts could also be evaluated to better understand how socioeconomic factors and global markets drive land-use change and affect multiple ecosystem service outputs. Whereas considerable investments have been directed toward developing commodity crops on private lands, few studies have evaluated how these investments have affected ecosystem services that benefit the public (39, 40). Pollinators serve as effective model organisms for evaluating ecosystem service tradeoffs because their service to humans is directly quantifiable (2, 41) and their health and provided pollination services can be linked with land management activities. Conservation efforts designed to promote habitat for pollinators in the NGP will likely benefit other ecosystem services, including conservation of biodiversity; however, these added benefits need to be quantified so that informed policy decisions can be made that maximize ecosystem service delivery while reducing ecosystem disservices from specific types of agricultural practices.

Future research is needed to understand how land-use change affects honey bee colony health, productivity, and pollination services. Similar to life cycle analyses conducted for naturally migrating species (42), models are needed to guide conservation investment throughout the migratory range of commercial honey bees. To maximize conservation investments, land management activities designed to benefit pollinators should be developed within an adaptive management framework so that management uncertainties can be addressed during the early stages of program development. In addition to quantifying large-scale habitat features that pollinators require, finer scale studies are also needed to investigate floral resources that maximize benefit to pollinators and will grow readily in agricultural landscapes (43 –45). This information can be useful for designing and evaluating conservation seed mixes that are cost-effective for implementing across large spatial extents and regional programs. Integrated ecological and economic models are also needed to evaluate how land-use change in one part of the country affects ecosystem service delivery elsewhere in the United States. Development of such models would be useful for identifying stakeholders who may directly benefit from pollinator habitat enhancement in the NGP because healthy honey bee colonies are required for agricultural crop pollination elsewhere in the United States.

We created maps highlighting (i) the spatial distribution of biofuel crops (i.e., corn and soybeans) and (ii) changes in biofuel crops and grassland area surrounding commercial apiaries in North and South Dakota from 2006 to 2014. We focused on these years because they represent a period of significant land-use change in this region, including loss of CRP and expansion of biofuel crop production (22). We obtained spatial locations of 18,363 registered apiaries from the North Dakota (number of apiaries, 11,629) and South Dakota (n = 6,734) Departments of Agriculture (Fig. 1) (data accessed January 12, 2015). In a Geographic Information System (R Core Team 2015, packages rgdal, rgeos, raster, sp) (SI Appendix), we georeferenced and placed a 1.6-km (∼1.0 mile) buffer around each apiary location and quantified annual land covers as classified in the Cropland Data Layer (CDL) (46) within each buffer. We used 1.6 km as buffer distance because commercial beekeepers generally maintain a distance of >3.2 km between apiary locations to minimize colony competition for floral resources. We extracted pixel counts of each CDL land-cover category in Geospatial Modeling Environment, Version 0.7.4.0 (47) and converted these counts to area using annual CDL resolution. We created two new land-cover classes, biofuel and grassland (Table S1), and summed the area values of contributing land-cover categories to calculate area (ha) in biofuel crop and grassland land covers for each registered apiary location and year. We then calculated the annual gains or losses of biofuels and grassland area for each apiary and calculated mean annual change for both land covers for each apiary across the entire study time span. We used inverse distance weighting interpolation in ArcGIS Desktop (48) to create spatial maps of biofuel crop production in 2006 and 2014, the annual rate of change in biofuel cropping area from 2006 to 2014, and the annual rate of change in grassland area from 2006 to 2014. All models used land-cover data from the 18,363 registered apiary locations to create the interpolation surface across the Dakotas. To estimate annual rates of biofuel change from 2006 to 2014, we constructed a linear trend model within a Bayesian framework, with either biofuel cropping area or grassland area at apiary i in year x as the response and YEAR as the predictor variable. Trend models were fitted using WinBUGS (49) and the R2WinBUGS package (SI Appendix) in R (50). For both models, we used normally distributed priors with zero means and large variances (i.e., diffuse priors). We report parameter estimates for the YEAR regression coefficient and associated 95% Bayesian credible intervals. YEAR regression coefficients that do not overlap zero would suggest a systematic trend in biofuel crop or grassland area change surrounding apiaries from 2006 to 2014. To highlight areas of potential landscape stress for managed pollinators, we overlaid a county-level apiary density map with each one of our interpolated land-use change maps. These maps revealed areas within the Dakotas that have a high density of registered apiaries and significant rates of increase or decrease in biofuels and grassland area change from 2006 to 2014. For each county, we report the number of registered apiaries per 10,000 ha (38.6 mi²). Because of differences in the apiary registration process for each state, apiary density should be interpreted as relative density within each state, rather than as comparisons across states.

We developed apiary selection models by identifying commercial beekeepers who operate across a broad geographic area, including portions of the Dakotas that have experienced significant gains in biofuel cropping area and grassland loss, as determined from our land-cover trend analysis. We focused our analysis on three large-scale commercial beekeepers who operate in eastern and central portions of the Dakotas. We used the North and South Dakota Departments of Agriculture apiary registration databases to delineate the operating domain of individual beekeepers. Within these domains, we conducted aerial photograph interpretation of all registered apiary locations (∼1,500 apiaries) to verify that the apiaries were used from 2005 to 2007. We used 2005 to 2007 as our study period to correspond with the 2006 CRP enrollment data we obtained through a memorandum of understanding with the USDA Farm Service Agency (FSA). This time period represented the height of CRP participation in the Dakotas, when ∼1.3 Mha (3.2 million acres) were enrolled. We used a combination of aerial images from Google Earth (2016 Google) and the USDA National Agricultural Imagery Program (USDA Geospatial Data Gateway, https://gdg.sc.egov.usda.gov) to determine whether an apiary was used during a given year. Our aerial interpretation revealed 644 apiaries that were verifiably used from 2005 to 2007. We removed all sites that were within 3.2 km of another known apiary to avoid pseudoreplication. Of these remaining sites (n = 583), we set aside one-sixth of the occupied sites (n = 98) for model validation and used the remaining 485 sites for model training. Within the operating domain of the commercial beekeepers, we randomly generated 800 points to represent unused apiary locations. Unused sites had the following selection criteria: (i) could not be located in a body of water or an urban center or on restricted federal lands, (ii) had to be within 50 m of an access road, and (iii) could not be within 0.8 km of each other or any known apiary locations. We applied this separation distance to minimize overlap with other used or unused sites. We set aside 98 randomly generated unused sites for model validation and used the remaining points for model training. This selection process yielded 1,183 and 196 sites for model training and validation, respectively. For all used and unused apiary sites, we used the National Agricultural Statistics Service (NASS) 2006 CDL to quantify land-cover and land-use features within 1.6 km of the point location. Similar to the land-cover trend analysis, we combined various land-cover and land-use categories into broader classes to reflect their hypothesized relationship for supporting commercial apiaries (Table S1). We reclassified the CDL by reassigning the original pixel value of each land-cover category to a new value representing one of the broader classes. Land-cover categories that occupied <0.5% of the landscape and were not easily assigned to our broader land-use categories, such as double crop classes, were excluded from quantification (Table S1). We determined the area of each land-cover class by extracting pixel counts within each apiary buffer in Geospatial Modeling Environment, Version 0.7.4.0 (47) and converting counts to area using the 2006 CDL spatial resolution. We also calculated the area of CRP and other private lands enrolled in federal conservation programs and all federal- and state-owned lands (Table S2). Shapefiles of federal and state lands were merged into a single layer, and both this layer and CRP were rasterized and reclassified to binary rasters to reduce processing time. Rasterization can cause loss of polygon edge definition; however, our 900-m² pixel resolution minimized potential edge effects. Area of CRP and federal and state land within each apiary buffer were determined by the same method used for the CDL.

We developed three logistic models to quantify how apiary site selection was influenced by land cover and land use. The first model (COMMODITY) was designed to assess how apiary site selection was affected by major commodity crops, which we classified into two broad categories: biofuel crops (covariate Biofuels, corn and soybeans) and small grains (Sm_Grains) (Table S1). The second model (HABITAT) included a Grassland covariate (Table S1). We also included cultivated crops and other land covers with suspected benefit to honey bees: Alfalfa, Forest, Water, and Sunflower (12). We did not include canola fields in the model because of a general lack of canola in our study region. The third model we created, CONSERVATION, quantified the role federal and state conservation lands play in influencing apiary site selection. We included CRP as a separate covariate because of the sizable amount of private land enrolled in CRP in the Dakotas. All other federal and state lands were combined into a single Fed_State covariate. We constructed a Pearson’s correlation matrix of all raw covariates before analysis; covariates with a correlation coefficient >0.3 were not included in the same model. All covariates were then scaled to have a mean of zero to allow for comparison of slope parameters generated from the regression models.

We developed all models within a Bayesian framework to allow for posterior prediction of used and unused sites during model validation. We fitted logistic regression models using WinBUGS (49) and R2WinBUGS (SI Appendix) in R (50). Logistic regression was used because our response variable was binary (i.e., 1 = Used apiary, 0 = Unused, randomly generated point) and land-cover predictor variables were continuous (see SI Appendix for model code and covariates). For all models, we used normally distributed priors with zero means and large variances (i.e., diffuse priors). We evaluated the 95% credible intervals of the slope coefficients to determine association between site use and habitat covariates.

We used the inverse of the logit-link function to predict apiary use probability for all 196 validation sites based on the slope parameter estimates generated from each model. We used the package pROC (SI Appendix) in R to calculate receiver operating characteristic (ROC) curves and integrated the area under the curve (AUC) to assess model performance and predictive capabilities (51). A model with perfect predictive power would yield an AUC of 1.0, and a model with no predictive power would yield an AUC of 0.5.