The contribution of land cover change to the decline of honey yields in the Northern Great Plains

Decreased availability of forage, as well as increased pesticide exposure, are important factors in the decline of honey bee health. Here, we isolate land cover transitions and their effect on honey production at 160 commercial apiaries in the Northern Great Plains. We found that land cover changes from 2008 to 2012 caused an annual decline in honey yields of 0.9% in the study area. Transitions from grassland to soybean (but not corn) were particularly detrimental to honey yields, potentially due to bee contact with pesticides within and around agricultural fields. When our results are applied to known apiary locations across all of North Dakota (U.S.A.), we estimate a 2.5% (1.6 million USD) decline in 2012 honey yields due to land cover changes occurring between 2008 and 2012. Even when controlling for changes in land cover, we found that on average colonies in the study area experienced a 14% annual decline in honey yields. We discuss possible explanations for these non-land-cover-related honey yield declines, including changing economic conditions (e.g. pollination services), changes in land management (e.g. pesticides), and increases in pests or diseases.


Introduction
Crop pollination is an economically critical ecosystem service. Bees are the predominant pollinators for a variety of specialty crops, such as fruits, vegetables, nuts, and oilseeds [1,2]. In 2019, annual revenues from managed honey bees (Apis mellifera) for pollination services in the United States were over $300 million [3]. In that same year, honey bees produced over $385 million worth of honey and other products [3].
Land covers with high-quality forage contain flowering plants that provide abundant and diverse nectar and nutritious pollen [11]. Land meeting those criteria could be wildland, restored land, or land cultivated with pollinator-friendly crops (e.g. alfalfa). For example, the U.S. Department of Agriculture's Conservation Reserve Program (CRP) incentivizes farmers to remove land from agricultural production and plant grasses or pollinator-friendly covers that have been shown to improve colony health [12].
A review of 30 years of land cover changes in the U.S. found that forage resources in the Northern Great Plains, a region of great importance to . Each apiary has a 2 km radius buffer around it to mimic the foraging range of a honey bee (Apis mellifera). In general, colonies that are producing considerably more than the mean honey production that year are mostly located in areas with pasture and grassland cover, yet that area of our study region is undergoing rapid land cover change with the westward expansion of the U.S. corn and soybean belt. Data for 2012 are not shown as there is no discernable pattern that year. Source for land cover data: USDA Cropscape. honey bee rearing [13], increased from 1982 to 2002 but declined from 2002 to 2012 [13]. Analyses using more recent data support the finding that broad-scale changes in land cover may have detrimental impacts on pollinators [14][15][16][17][18] and that habitat loss resulting from land conversion drives wild-bee declines [5,19,20].
Previous research has focused on the influence of land cover on bee health [16] and the suitability of sites for apiaries [13,15]. A recent study by Smart et al [18] estimated the impact of land cover on the survival and productivity of honey bee colonies at six commercial apiary sites over 3 years. The small sample size (necessary due to field data collection constraints) enabled the authors to assess the types of pollen collected and the presence of pesticide residues, in addition to honey yields. While our study does not include this kind of detailed field data collection, using data from a commercial beekeeping operation allows us to assess the role of land cover change on honey yields for a much larger sample than encountered in previous studies. We empirically estimate the relationship between land cover and honey yields for over 5000 colonies at 160 commercial apiary sites over 4 years (figure 1). With this unique dataset, we can control for unobserved apiary-level confounding factors. This approach allows us to validate forage and insecticide toxicity indices (ITI) introduced in recent studies [19,[21][22][23][24] and to determine the effect of land cover changes on honey yields. Because we have multiple years of data from the same locations over a relatively broad geography, we can ascertain the effect of specific land cover transitions on honey production and scale the results to this major honey production region of the U.S.

Study area
The Northern Great Plains is of great importance to managed honey bees as beekeepers rely on the region for honey production and refuge for their bees. During the summer, approximately one-third of managed honey bees in the U.S. reside and produce honey in the Northern Great Plains [13,25]. North Dakota in particular is the state with the most honey-producing colonies and the highest density of apiaries in recent years [3]. We focus our work on a five-county area in southcentral North Dakota (Kidder, Logan, LaMoure, McIntosh, and Stutsman counties; figure S1 (available online at stacks.iop.org/ ERL/16/064050/mmedia)), for which we have apiary location and colony count data, as well as associated honey-production records for 2008, 2009, 2010, and 2012.

Apiary data
We obtained honey-production data from a commercial honey-production and pollination business 6 . Data include locations of between 106 to 158 honey bee apiaries each year with spring, summer, and fall colony counts at each apiary, and the number of honey supers collected. The beekeeper placed between 27 and 56 colonies (40.9 ± 4.6 mean ± 1 SD) at each apiary over the 4 year period. Honey supers consist of frames holding combs used to store honey in a hive, the managed boxes in which a single colony lives [26]. Precise honey weights tied to a specific location are not available as supers are mixed during the extraction process. The number of supers collected is used by commercial beekeepers to estimate profits and while there may be some variation in the honey stored on each frame, supers are generally consistent in their weight. Our commercial honey producer exclusively uses shallow supers, which hold 28 lbs of honey on average, accounting for the fact that the last super is not always full when the honey is harvested.
Our analysis uses the number of supers collected divided by the number of colonies at each apiary in each year (table S1), which we refer to as honey yield. Records for 2011 were unavailable for this analysis because they had been misplaced by the beekeeper. We exclude from our analyses colonies that were moved to a new apiary location during a particular year. If the total number of colonies placed at a yard in the spring decreased by more than 35% when the season ended in the fall, the yard was removed from the subsequent analysis as it was deemed a failing apiary (most likely due to mite infestation or other diseases). This resulted in the removal of six apiaries in 2009, three in 2010, and ten in 2012. Public land survey system (PLSS) data were used to geographically locate the apiaries. The township, range, and section from the honey production data were matched to the PLSS dataset and the PLSS sections were converted to points for each apiary. Thus, the apiary locations are only approximate. Given that within a PLSS section the land cover in our study area is rather homogeneous and that we consider land cover well beyond the PLSS section in our analyses, we expect this issue to be minor.

Land cover data
Land cover data at 30 m spatial resolution were down- (CDL) are raster datasets classified into 40-43 land cover types, depending on the year.

Forage resources index
We build our model based on InVEST (Integrated Valuation of Ecosystem Services and Tradeoffs), which is a suite of models for valuing ecosystems services [21]. Here, we modified the pollination module of InVEST, developed for wild bees [21], so that it is applicable to managed bees in that only forage resources, not nesting resources, are considered. We call this model ABEILLES, which is French for bees, and stands for 'Apis, other Bees, and Environmental Indicators in a Land-use Land cover Suitability' model. The model creates a spatially explicit forage resources index (FRI), using a land cover map, information about floral resources in each land cover class, and an estimate of honey bee foraging range.
Each land cover class is assigned a value from zero to one that represents the estimated availability of nectar and pollen within that land cover by season. For example, water provides no floral resources for bees and therefore is assigned a score of zero in all seasons, while sunflower is assigned a summer score of 0.434. These scores are based on work by Koh et al [19] modeling wild-bee abundance using national land cover data. For the land covers present in our study area but absent from Koh et al [19], e.g. CRP categories, we conducted our own expert assessment of floral resources (table S2). We set the foraging range to 2.5 km based on Steffan-Dewenter and Kuhn [27]. We calculate the mean FRI surrounding each pixel within the foraging range, to which we apply a distance decay function.
The model is run for each season, which is given a relative weight depending on its importance to the provision of floral resources that are attractive to bees. In North Dakota, floral resources for honey bees peak during the summer and decline into the fall. As honey bee colonies in the study sample were not present on the landscape in the spring, we consider only the summer and fall seasons, weighted as contributing 60% and 40% of the final model output, respectively. These weights match those used in similar modeling studies [19,21,23]. The output of the model is representative of the suitability of forage within the landscape for honey bees.

Insecticide toxicity indices
The insecticide toxicty indices (ITI) are created similarly to the forage resource indices. For these indices, each land cover class is assigned a value from zero to one that represents the honey bee toxicity weighted insecticide use. Our scores are based on work by Douglas et al [22] on toxic load. The insecticide use data is from the USGS Pesticide National Synthesis Project and includes both foliar and seed treatment applications 7 . The toxicity indices are created by dividing the insecticide usage by the honey bee acute toxicity values. These toxicity values are LD50s (dose producing 50% mortality) compiled from the US-EPA's ECOTOX database and the Pesticide Properties Database. The compiled toxicity indices (i.e. toxic loads) are downloaded from the Insecticide Explorer website. Data reclass tables for the USDA Cropland Data Layer are downloaded by year for North Dakota (see table S3). These scores are then assigned to the CDL pixels and the distance weighted index is calculated for each apiary in our study. There are two versions of the index: contact-ITI and oral-ITI. Oral exposure occurs when the insect eats pollen or nectar laced with insecticide. Contact exposure occurs when the insect is thought to visit the field but may not eat the nectar or pollen. An index of zero would be an apiary located in a landscape without any insecticides applied to any crops within the 2.5 km foraging radius. An index of one would be the highest possible toxic load. This would occur in our study area and time-period if an apiary was surrounded exclusively by soybeans, and insecticides were applied at the rate they were in 2008, which was the year with maximal insecticide application during our study period.

Vegetation growth index
Aside from changes to land cover, annual changes in precipitation and temperature have the capacity to alter floral environments and the timing of floral emergence, which may in turn impact where beekeepers place their hives [28]. In this study, we use the Normalized Difference Vegetation Index (NDVI) to account for seasonal and annual fluctuations in floral production.
The NDVI data relate to the growth condition of vegetation on a landscape and is akin to an index of 'greenness.' It varies from negative one to one with values less than zero for areas of open water, values greater than 0.2 generally representing vegetation, and the higher NDVI values indicating areas of actively growing plants [30,31]. Land-cover factors that influence plant growth, such as temperature and precipitation, affect NDVI over time.

Empirical methods
To assess the impact of changes in land cover, FRI, and ITI on honey yield, we use a linear mixed modeling approach. This approach takes advantage of the repeated measurements of honey production at each apiary to control for any time-invariant factors affecting honey yields (e.g. soil). In addition to these apiary-level controls we also include a year effect to control for possible confounding factors that are correlated across our study area during each honey production season (e.g. temperature, honey prices). Several random-effects structures were tested and only the best random-effect structure is reported.
We calculate the changes in honey yield due only to land cover changes over time by fixing the years and NDVI and estimating honey yield. This isolates the land cover effect from the year and NDVI effects. In addition, we project changes in honey production to the state of North Dakota using apiary location data obtained for the entire state.

Results
Land cover varied greatly over our study area, with some of the apiaries in landscapes dominated by cultivated crops while others are placed in landscapes dominated by pasture and grassland (figure 1). We first examine the correlation between honey yields and the distance-weighted proportion of different land covers surrounding an apiary. Apiaries near more corn and soybeans had lower honey yields (figure 2(a)), while apiaries near more pasture and grassland had higher honey yields ( figure 2(b)). The correlation between flowering crops and other flowering land covers (i.e. alfalfa and sunflower, clover, and wildflowers) and honey yields was even more pronounced ( figure 2(b)). However, these flowering land covers were uncommon near the apiaries (on average only about 1% of the landscape surrounding the apiaries).
Honey yields were higher in landscapes with more forage resources. Simple correlations between honey yields and FRI, ITI, as well as NDVI are shown in figure 3. Using our linear mixed model of forage resource availability, which controls for apiary-level confounding factors and year effects, we found that land covers with a higher FRI had significantly higher honey yields (model 1a in column 2 of table 1). If FRI increased by 0.01, honey yields increased by 0.05 honey supers (1.4%). These results validate the use in the literature of FRI as an index of forage resources available to honey bees [19,21,23,24].
We also found that apiaries near landscapes with a higher contact-ITI had lower honey yields (model 2b in column 4 of table 1). If contact-ITI increased by 0.01, honey yields decreased by 0.04 supers (1.1%; table 1). We found no significant effect of the oral-ITI on honey yields (model 2a in column 3 of table 1). Since we do not have field data on pesticide uptake by the bees in our sample, we cannot rule out other mechanisms that could explain our finding that higher contact-ITI reduces honey yields. For example, pest presence could increase pesticide use but also interfere with crop blooming and therefore nectar production.
Conversion of grasslands to other land covers over time affected honey production. In our third set of models, we regressed honey yields on distanceweighted land cover proportions within the foraging range of honey bees, as well as NDVI and year indicator variables (models 3a and 3b columns 5 and 6 of table 1). To avoid collinearity, the most common land cover (i.e. grasslands) was not included in the model. Therefore, our estimates of the effect of land cover on honey yields are all relative to grasslands, our base land cover type. We estimated that land cover changes between 2008 and 2012 reduced per-hive honey yields for our study area on average 0.9% per year. This impact was large considering that broad land cover changes tend to occur over longer time frames than the 4 years in our study [24]. Soybean production had a significantly negative impact on honey yields with a 1% change from grassland to soybean associated with a decline in honey yields of 0.03 supers (0.9%; model 3b in column 6 of table 1).
Other factors such as annual variations in floral resources due to weather and unexplained declines over time also impacted honey yields. Across all our models, the NDVI had a positive, significant effect on honey yields. In the first model, a 0.01 increase in summer and fall NDVI was associated with an increase in honey yields of 0.02 supers (1.4%; model 1a in column 2 of table 1).
Lastly, for the hives in this study, we find a 14% annual decline in honey yields independent of land cover changes and weather variability captured by NDVI. This could be due to several factors, which we explore in the discussion. For comparison, honey yields in all of North Dakota fell from 90 to 69 pounds per colony from 2008 to 2012, an annual decline of 6% [3].
Declines in commercial honey yields in the Northern Great Plains are due, at least in part, to land cover changes. The Northern Great Plains serve as a summer refuge for one-third of commercial honey bee colonies in the U.S. and are an important region for honey production [15]. North Dakota, in particular, is the leading honey production state in the U.S., accounting for just under a quarter of total U.S. honey production in any given year [3]. In 2012, honey production in North Dakota was valued at $64 million. Using our estimates of the impacts of land cover change on honey yields and the location of apiaries in North Dakota, we estimated the effect of land cover changes between 2008 and 2012 on honey production

Discussion
We find that grassland conversion to soybeans had a significantly negative effect on honey yields, yet that is not the case when grassland transitioned to corn. This finding is somewhat unexpected, as soybean fields provide marginally more floral resources for bees than corn [32,33]. A potential explanation might lie in differences in pest management [34], especially since our best model includes the contact-ITI.
As of 2009, 65% of the corn crop in the U.S. was planted with genetically modified, plantincorporated protectants [35] and to date, these are not commercially available in soybeans. This is important because crops with plant-incorporated protectants are associated with a decrease in foliar and soil insecticide applications (see figure 4 in [26]). In addition, Coupe and Capel [35] reported a dramatic increase in insecticide applications on soybeans from 2004 to 2009 to treat a new pest, the soybean aphid (Aphis glycines Matsumura), that had been found in the U.S., including in our study area. When economic thresholds of pest damage are reached, which can be during the flowering stage, soybean aphid is treated with insecticides [36]. In addition to the exposure of pollinators to these insecticides within fields, foliar application of insecticides can drift from the field to non-target flowering plants [37]; both of which would result in 'contact' insecticide exposure. This difference in the management of pests in corn and soybeans may explain some of the differences in the effect of these land cover changes on honey yield.
Additionally, during our study period, the use of seed-treated insecticides became more common. While this was the case for both corn and soybeans, the active ingredients used in the study area across these two crops differed. For example, clothianidin was commonly used in corn as a seed treatment, while its use in soybean is largely absent [38]. Instead, imidacloprid was commonly used in soybeans [38]. While all of the neonicotinoids have the potential to impact honey bees, the effects of each active ingredient can differ [39]. Furthermore, because the primary Table 1. Impact of forage resources, insecticide toxicity, and land cover on honey yields. Honey yields are the mean number of supers (honey collection box with 8-10 frames) per colony. All indices are unitless and range from 0 to 1. Details of how the indices are calculated can be found in the methods section. The p-values were obtained by bootstrapping the estimates (1,000 replications) at the apiary-level (i.e. clustered) to obtain a non-parametric estimate of the uncertainty distribution of the coefficients. Asterisks * , * * , and * * * denote 10%, 5%, and 1% levels of significance, respectively when testing the null hypothesis that the coefficients equal zero. The dots (.) indicate the base year or land cover in the model. AIC is the Akaike information criterion. route of exposure to bees from neonicotinoid seed treatment is through off-target airborne dusting during planting, differences in the planting methods (e.g. vacuum planter, seed flow lubricant), planting date, and the ability for the coating to remain on the seed has the potential to affect non-target flowering plants near corn and soybean fields differently [40,41]. While this route can cause direct exposure for nontarget insects near these fields during planting, the commercial honey bees in our study are moved into the area after corn and soybean planting. Insecticides can also act synergistically with temporal and spatial floral-resource scarcity. While soybean fields provide nectar for honey bees during the bloom period [33], pesticide exposure during this period can have lasting effects on honey bees' ability to collect nectar and pollen in later periods. Tosi et al [6] report that when sub-lethal neonicotinoid exposure was associated with food deprivation, honey bee mortality increased by 50% compared to mortality from exposure alone and others show that 'Chronic exposure [to thiamethoxam, a neonicotinoid insecticide] significantly decreased flight duration (−54%), distance (−56%), and average velocity (−7%) after either one or two days of continuous exposure that resulted in [honey]bees ingesting fieldrelevant thiamethoxam doses' [42]. These sub-lethal foraging impacts could reduce honey yields without causing colony mortality. Aside from acting as a potential source of pesticide exposure to honey bees, areas with substantial soybean or corn cultivation reduce the availability of plants that flower during later periods, which affects the ability of honey bees to prepare adequate honey stores for the winter [33].
The effect of herbicides is more insidious. Herbicide applications are meant to suppress the growth of weeds within fields, which are known to provide nutritional resources for bees in agricultural systems [43]. Given that from 2008 to 2012, glyphosate applications steadily increased on soybean fields in the U.S. [44] and the average use of glyphosate per acre is higher in soybeans than in corn [45], this may also serve to explain why we see a large decrease in honey yield when land cover transitions from grassland to soybean, but not to corn.
On the whole, the Northern Great Plains is a rapidly changing landscape [14,15,46,53]. The area has undergone a conversion from grassland to corn and soybean crops (1%-5.4% annually from 2006 to 2011 in the Northern Great Plains [31]) and continues to do so (figure S2). The loss of grasslands is associated with the loss of foraging resources [15], which affects honey bee health [17,18]. Interestingly, others have found that honey bees 'that did not undergo nutritional stress were not significantly impaired' by sublethal doses of the two neonicotinoids (clothianidin and thiamethoxam) tested in their study [6]. Similarly, in laboratory experiments involving honey bees, 'high-quality diets (polyfloral pollen and highquality single-source pollen) have the potential to reduce mortality in the face of infection with Israeli acute paralysis virus' [47]. It is likely that the combination of changes in the forage area, with the increased insecticide and herbicide applications, acts together to decrease honey yields.
Many studies that consider the effect of landscape composition on honey bee health [16,18,48,49] did not differentiate between corn and soybean fields. Our findings might explain seemingly incongruent results: e.g. where Sponsler and Johnson [49] reported that the proportion of cropland within a 1 km radius of honey bee colonies in the midwestern USA was positively related to food accumulation in the colony, yet Smart et al [18] found a negative impact of combined cultivated acreage on honey bee-colony health in the Northern Great Plains region of the U.S.A.
We also found a substantial (14%) annual percolony decline in honey production unrelated to land cover changes. While weather is known to affect honey production [50,51], and flowering plants that honey bees depend on fare poorly in excessively dry or rainy seasons [29,52], we aimed to control most of that variability using the NDVI. Factors other than weather must be at play, such as a worsening of mite infestations or nosema infections, or a decrease in the effectiveness of pesticides and fungicides used to treat them. Our lack of field data prevents us from probing these questions further, though a recent study detected resistance to the miticide amitraz and reduced efficacy of Varroa mite (Varroa destructor) treatment in commercial beekeeping operations [53]. In addition to environmental factors, honey yields can also be impacted by changing economic conditions, as honey bees are used for both crop pollination and honey production with competing demands (detailed further in SI). Of course, combinations of all the factors described above acting in synergy are the likely cause of honey yield decline [4], but, to the best of our knowledge, our study constitutes the first attempt to isolate specific land cover transitions' effect on commercial honey yields.
While our study did estimate robust and consistent results, observational studies of commercial hives have limitations. For example, because we did not weigh each super, and thus we had to assume that the supers were of uniform weight. We did confirm with the beekeeper that he and his crew only remove full supers except during the last harvest when the uppermost super may be only partially full of honey. This may lead to some unexplained variation that could bias our results. However, in the absence of a state-or regional-level dataset on colony location, honey production, and beekeeper management practices, focusing on data provided by an individual beekeeping operation is the next best option. The commercial beekeeper who provided the data for our study implemented a similar management strategy across all apiaries. This included control of mite infestations when they occurred. This would reduce variation caused by other non-land cover change factors. Other factors, such as pesticide exposure, should be validated with residue analysis.

Conclusion
Our research provides empirical evidence that recent land cover changes in the Northern Great Plains have caused a substantial and sustained decline in honey yields over a short time period. Of particular significance are the honey yield declines associated with a land cover shift from grasslands to soybeans (Glycine max), which is not observed with transitions from grassland to corn (Zea mays). This research also helps to validate both floral resource and insecticide toxicity indices currently used in the literature. Our results provide evidence of the value to honey bees of grasslands, such as those established under the U.S. Department of Agriculture's CRP. However, because our results are based on observed commercial honey production records, future validation of our results in the field is necessary.
The main strength of our study is our ability to isolate the effect of specific land cover transitions, such as grasslands to corn or grasslands to soybeans. We are also, to the best of our knowledge, the first to explicitly account for annual fluctuations in floral production associated with weather using NDVI, which we find to be a significant predictor of honey yields.
Lastly, we validated the output of ABEILLES, a model of forage resources for bees that is related to the pollination module of InVEST. Our results indicate ABEILLES can be helpful to generate forage resource metrics for other areas of the U.S.

Data availability statement
The data that support the findings of this study are available upon reasonable request from the authors.