Analysis of Pollination Services Provided by Wild and Managed Bees (Apoidea) in Wild Blueberry (Vaccinium angustifolium Aiton) Production in Maine, USA, with a Literature Review

Abstract

Maine is the largest producer of wild blueberry (Vaccinium angustifolium Aiton) in the United States. Pollination comes from combinations of honey bees (Apis mellifera (L.)), commercial bumble bees (Bombus impatiens Cresson), and wild bees. This study addresses (1) previous research addressing wild-blueberry pollination, (2) effects of wild-bee and honey-bee activity densities on fruit set, yield, and crop value, (3) the economic value of wild-bee communities, and (4) economic consequences of pollinator loss. Bee communities were sampled in 40 fields over three years (2010–2012) and bee activity densities were estimated for bumble bees, honey bees, and other wild bees. These data were applied to an economic model to estimate the value of bee taxa. Bumble bees and honey bees predicted fruit set and reduced its spatial heterogeneity. Other wild bees were not significant predictors of fruit set. Yield was predicted by fruit set and field size, but not pest management tactics. Our analysis showed that disruption in supply of honey bees would result in nearly a 30% decrease in crop yield, buffered in part by wild bees that provide “background” levels of pollination. Honey-bee stocking density and, thus, the activity density of honey bees was greater in larger fields, but not for wild bees. Therefore, a decrease in crop yield would be greater than 30% for large fields due to the proportionally greater investment in honey bees in large fields and a relatively lower contribution by wild bees.

1. Introduction

Lowbush blueberry (Ericaceae: Vaccinium angustifolium Aiton), the dominant plant in the wild blueberry species complex (see Appendix A), contributes prominently to Maine’s economy. Maine is the largest producer of wild blueberry in the United States [1], a crop that, from 2012–2018, averaged 39.3 ± 9.2 (SD) million kg berries/year and was worth an average of USD 3.85/kg [2]. With over 8000 harvested hectares each year [2], blueberry crop production values are only exceeded by northern Maine’s potato industry [3].

Wild blueberry is primarily an outcrossing native North American fruit crop that relies on insects to vector pollen from one clone (genetically unique individual plant) to another [4]. Bees are the primary pollinators of wild blueberry (see Appendix A). Bushmann and Drummond [5] documented more than 120 bee species associated with wild-blueberry fields. This is approximately 43% of the total 278 bee species recorded in Maine [6]. On the basis of a single bee flower visit, bumble-bee queens (Bombus) are the most efficient pollinator of wild blueberry, followed in efficiency by bees of the genera Andrena and Osmia [7,8]. Before 1960, most wild-blueberry farmers relied upon resident wild bees for pollination [9]. Honey bees (Apis mellifera (L.), which are not native to North America, were evaluated as pollinators in the 1950s–1980s [9] and are currently used by the majority of wild-blueberry farmers in the US and Canada despite many growers relying on wild bees for at least part of their pollination [10,11,12] (also see Appendix A). The commercially available common eastern bumble bee (Bombus impatiens Cresson) and the alfalfa leafcutting bee (Megachile rotundata Fabr) also provide pollination that supplements that provided by wild bees [13,14,15,16,17]. Currently, in Maine, a substantial proportion of farmers employ these alternative commercial bees or simply rely on resident wild pollinators [10,11,12].

A few studies compared the field-level pollination efficacy of wild bees and the managed honey bee. Asare et al. [18] found that, on a per bee basis, honey bees and wild bees increase fruit set by 0.8% and 11.0%, respectively, for each increase in a unit bee. On a field basis, they estimated that honey bees contribute more to fruit set than wild bees owing to their greater numbers. An economic analysis of pollination by Hoshide et al. [19] showed that the net income for Maine wild blueberry farmers attributable to wild bees was USD 613/ha and that for honey bees was USD 923/ha. They also estimated attributable net income per bee, calculated using observed wild-bee densities and honey-bee hive bee counts. For rented honey bees, pollination value was USD 0.012/bee. Wild-bee pollination value per bee was estimated at USD 0.311/bee. Attributable net income values per wild bee were greater than those for rented honey bees due to wild bees’ greater pollination efficiency.

In Section 1 and Appendix A, we review published studies on wild-blueberry reproductive biology, pollination, and bee communities. Despite all of these studies (>80), there is only limited knowledge of the role that the wild-bee community plays in wild-blueberry pollination. Our study of wild-bee communities in blueberry fields conducted between 2010 and 2012 will expand understanding of wild-blueberry pollination. Specifically, we ask the following questions: (1) How does wild-bee abundance (bumble-bee queens and other wild bees) and species richness affect mean fruit set and its variance and fruit yield, relative to the honey bee? (2) Can we estimate the economic value of the wild native bee community as a whole, and which taxonomic group(s) most contributes to this value?

We ask these questions about a production system that depends on commercial pollinators while wild bees forage freely in fields. Given adequate populations, wild pollinators can provide sufficient pollination services for wild-blueberry production [18,20] or serve as additional pollinators in fields that use commercial bees [14,17,18]. However, given the past, current, and predicted future decline of wild-bee populations, densities sufficient for primary or even background pollination services cannot be assumed in the future [21,22]. What is the economic value of wild bees for an insect-dependent crop? The answer to that question will help us understand the value or cost of considering pollination strategies [23,24,25].

2. Materials and Methods

2.1. Fields

With this study, we chose 40 Maine wild-blueberry fields located in Hancock and Washington Counties in 2010–2012. Nearly 88% of the wild blueberries produced in Maine originate from these two counties [2]. These are the northeastern counties that comprise Maine’s “Downeast” region, created by retreating glaciers and characterized by sand and silt glacial deposits [26,27,28]. Fields were visited over a three-year period: 2010 (12 fields), 2011 (12 fields), and 2012 (16 fields). Fields ranged from 0.5 to 16.3 ha with a median field size of 3.4 ha. Several field characteristics were related to field size (see Section 3); thus, in addition to using the continuous measure of field size, fields were categorized as small and large in relation to the estimated median field size (n = 40) in order to minimize correlation among predictor variables. Small fields (n = 22) included fields less than or equal to the median field size, while large fields (n = 18) were those that exceeded the median. Fields of both size categories were included in all years. The 2010 and 2011 fields were located in Hancock County. In 2010, fields ranged from 0.2 to 19.0 km (mean = 8.6 km) between fields. In 2011, fields ranged from 0.1 to 22.3 km (mean = 6.2 km) between fields. In 2012, the study area included northern regions of Hancock County and Washington County, which increased the distance between the fields. The 2012 fields ranged from 0.9 to 77.0 km (mean = 30.8 km) between fields. The distances between fields were measured on Google Earth (©2013 Google, https://www.google.com/earth/versions/#download-pro) from points originating in the middle of the field, without compensation for curvature of the earth.

Most of the fields were distinct, isolated blueberry fields surrounded by boreal spruce/fir forest [29] or non-crop open land, with two exceptions. In one case, a field opened on one edge to a non-study blueberry field that was in a different production cycle (i.e., bordering field was pruned and not in bloom) and, in the second case, a road separated two study fields that were under different management regimes.

Pollination management for most fields relied on honey bees rented by farmers (24 fields). Six fields were stocked with commercial bumble bees (Bombus impatiens) purchased by farmers, whereas seven relied on wild bees including rental honey bees stocked in neighboring fields (see [18] for prevalence). One field used a combination of rental honey bees and commercial bumble bees and two fields were serviced by two colonies of stationary, resident honey bees.

Each study year, we obtained yield measures (kg berries/field) from the cooperating farmers. For all unharvested fields (n = 2), we harvested berries by hand, raking in a measured area of the field and estimated kg berries from these subsamples of harvested yield. We determined field size from digitized records using ArcGIS version 10 [30]. With these data, we estimated kg berries/ha for 34 of the 40 fields. For all analyses pertaining to yield or crop value, the six fields without yield data were omitted.

2.2. Fruit Set

In early to mid-May of each year, nine representative clones (genetically distinct plants, see Appendix A) were selected. Clones typically range in size from <1 m² to >500 m², mean = 43.3 m² [31]. Three clones were selected near the center of the field, three about halfway between the edge and center, and three near the field edge. In order to compensate for microregional differences in growth, clones were selected throughout the field, although no clone was selected directly adjacent to honey-bee hives. Four stems were randomly chosen and tagged in each clone. Stems were visited three times. The first visit resulted in a count of the flowers per stem. The second visit took place in early June during which the ripening fruit was counted as a measure of fruit set, which we refer to as early fruit set. The last visit occurred in late June for a final count of developing fruit. The proportion of fruit at the final count that developed from the initial flowers is a measure of late fruit set, after some of the premature fruit drop has occurred [32]. Late fruit set is not necessarily a measure of harvestable fruit as harvest typically occurs 3–6 weeks later, a time period during which environmental conditions can further reduce fruit number. We also estimated the coefficient of variation (CV) of late fruit set for each field. This metric was used to measure the evenness of pollination between the nine selected clones [33]. Using least-squares linear multiple regression, we examined the influence of bee activity densities for all three bee groups (see Section 2.8) on measures of fruit set [34].

2.3. Economic Modeling

Losey and Vaughan [35] developed a mathematical expression to estimate the economic value of wild bees as pollinators of vegetable and fruit crops produced nationwide. Their formula is based upon United States Department of Agriculture reported national monetary values of vegetable and fruit crops, the reported dependency of that crop on insect pollination for fruit set, and the proportion of pollinating insects that are wild bees in order to calculate the level of pollination services attributable to wild bees.

Their equation states

V_np = Σ [V × D × (1 − P)],

(1)

where V_np is the annual value of the crop attributable to native pollinators, V is the annual value of the crop as determined by the US Department of Agriculture, D is the dependency of each crop on insect pollination, and 1 − P is the proportion of the effective insect pollinators (in abundance) that are wild native bee species.

We used the Losey and Vaughan [35] model by considering that, for wild blueberry, D = 1, since it is assumed no wild-blueberry fruit is produced without insect pollinators [4,9] (see also Appendix A). Therefore, D was eliminated from the equation. We also changed V_np and V to be equal to the value in US dollars found per blueberry field, not as a total of the entire blueberry industry, and we used the representation V_bee instead of V_np to stand for any bee taxon or group.

The modified equation states, for each field,

V_bee = V × (1 − P),

(2)

where V_bee is the value in US dollars of the blueberry crop for a given field for a given year that is attributable to a target bee group, V is the annual gross value in US dollars of the blueberry crop for a given field for a given year, 1 − P is the proportion of bee abundance represented by the target bee group found in a given field and year, and P is the proportion of bee abundance that belongs to the non-target bee groups. The gross value of crop (V) was farmer-reported and the proportions of bee groups were determined as described in Section 2.4.

We chose to consider three target bee groups: honey bees (HB), bumble bees (BB), and other wild bees (OB). The OB group includes all non-Apis, non-Bombus bee genera which are not commercially managed, including the cleptoparasitic bees (e.g., Sphecodes) and species that may not be native to North America [36], but these non-native species are not common in Maine wild-blueberry fields [6]. We did not consider pollinating insects other than bees in our assessments of total pollinator community, which excludes suspected but undocumented pollination by some Lepidoptera, Diptera, Coleoptera, and Formicidae [37,38].

We initially calculated the value of the three bee groups by giving equal weight to the pollination ability of honey bees vs. wild bees. This likely underestimates the value of the wild-bee community, as those are shown to be more efficient pollinators of blueberry than honey bees [7,8,18]. We then recalculated the value of each bee group by adjusting the proportions of the groups by pollination equivalencies determined by Javorek et al. [7] (see Section 2.6).

2.4. Determining Proportions of Bee Groups

The proportion of each bee group was based on timed activity-density counts of foraging bees observed in the field. For all three study years, over the course of blueberry bloom (early May–early June) and up to 10 days beyond bloom, each of the 40 fields was visited 3–4 times for measuring bee activity density. Observations were made only when bees were most likely to be active and foraging: mostly sunny conditions, with wind speeds below 1.6 kph and temperature above 12° C [8]. The activity-density counts of all three bee groups were obtained on the same visit.

In all years of the study, a bumble-bee (Bombus spp.) activity-density measure was obtained by one observer, working in a randomly chosen section of the field, counting foraging bumble bees while walking through the field along a straight line at a slow pace for 15 min. Care was taken to make sure bees were not recounted as they moved from one foraging spot to another or were disturbed by the observer. As wild-blueberry plants are less than one-half meter tall and queens (usually >1.5 cm in length) mostly forage near the top of the plant, this method of observation is likely to reveal all bees in the vicinity. The bumble bees were not identified to species, but workers from fields that stocked with Bombus impatiens were not counted. Workers of other species (determined by visual inspection without capturing the bees) and B. impatiens from fields not stocked with commercial colonies were counted when they appeared toward the end of bloom.

In 2011, foraging honey-bee activity density was assessed by trapping single bees in individual lidded cups during a 5–10 min period. The bees were released after counting. In 2012, honey-bee activity density was measured by placing a quadrat, 0.9 m per side, around a blooming patch and, after waiting one minute to allow the bees to settle, counting the number of bees foraging on flowers in one minute. This was done for 15 randomly chosen locations in the field. Field assessments of honey-bee activity densities were not recorded in 2010, but estimated from our observations in 2011 and 2012 (see Section 2.5).

In all years of the study, wild-bee collection involved groups of 2–3 people hand-catching bees found foraging on flowers. Each collection bout lasted 15 min. During the bouts, any non-Bombus bee found foraging on any flower in the field (including non-blueberry) was collected in a clean, lidded cup and placed on ice in a cooler until transfer to a freezer (−20 °C) where they were stored until pinning. Bumble bees were specifically avoided as only queens forage during blueberry bloom. All bees were identified to genus using the online key at Discover Life (https://www.discoverlife.org/bees/) and later identified to species [5]. From this collection, we calculated the numbers of bees caught per minute for each non-Apis and each non-Bombus genus. These measures include bees caught during bloom and those caught on non-crop flowers as bloom ended.

The activity-density measures for all three bee groups, expressed as bees/minute, account for all bee taxa found in the fields. From these activity-density measures, we calculated the proportion of honey bees (HB), bumble bees (BB), and other wild bees (OB) foraging in the field.

2.5. Determining Honey-Bee Activity and the Proportion of Honey Bees Foraging for 2010

We used linear regression to estimate values for honey-bee activity and the proportion of honey bees foraging in 2010 (see Section 3.1). The number of honey-bee hives placed in the fields was known for 2010. This model, described by Equation (4) (see Section 3), assumes a background proportion of 0.519 honey bees/field (y-intercept). This proportion was then used to calculate the value of the bee groups for 2010 fields.

2.6. Adjusting the Economic Value of the Bee Groups on the Basis of a Pollination Equivalency

Javorek et al. [7] calculated a pollination equivalency, i.e., the number of wild-blueberry flowers pollinated in the time it takes a honey bee to pollinate one flower, for five genera including Bombus (bumble bees) and Andrena. They determined that, for each flower pollinated by a honey bee, bumble bees pollinate 5.9 flowers and that an Andrena pollinates 3.6 flowers. This pollination equivalency for bumble bees reflects an average of values presented for queens and workers [7]. We used this average since a few workers can appear at the end of bloom, even though most foraging bees during blueberry bloom in Maine are queens. Using the pollination equivalencies, field measures of bees per minute for both Bombus and Andrena were recalculated, and then the proportions of honey bees, bumble bees, and “other bees” foraging in each field were recalculated. We considered pollination equivalency for Andrena because species of this genus were common in our collection and several species carried high levels of Vaccinium pollen [5]. We only considered the 34 fields from the full dataset for which we had yield values. For this analysis, we used the estimated honey-bee activity densities for 2010.

2.7. Determining Species Richness

The hand-caught solitary bee collection was supplemented with solitary bees caught in bowl traps during blueberry bloom and with a summer bumble bee collection that occurred after bloom. Details of this collection can be found in Bushmann and Drummond [5]. With this larger collection of bees, a measure of species richness defined as the total number of species was determined for all 40 fields.

2.8. Statistical Analysis

Data from all 40 fields were used for most analyses. All statistical tests were performed with JMP 14.0 [34]. We were not able to obtain yield for six fields (four fields from 2010, two from 2012). These fields are excluded from all analyses of yield and crop value. Linear Pearson correlation analysis was used to assess the association between the activity density of bee groups and measures of fruit set (early, late, and CV) over all years. To explicitly model the effect of year on bee group activity density and our three measures of fruit set, we used general linear models with year as a categorical independent variable. General linear models were constructed with year as a categorical independent variable and field size as both a categorical and a continuous independent variable to assess the effect of field size on bee group activity density, fruit set, and yield.

Using generalized linear mixed models, the influence of bee activity density was examined for all three bee groups (BB, HB, OB) and bee species richness on the three fruit set measures on yield. In these models, year was considered a random effect and all bee group activity density and bee richness predictors were considered fixed effects. The marginal r² (proportion of total variance explained by fixed effects only) was calculated to determine the total variance explained by the predictors. Model fit and appropriateness was evaluated with the chi-square convergence test. Fixed effects were tested for significance using Satterthwaite’s method for estimating denominator degrees of freedom in fixed effects F-tests in a mixed model [39].

Predicting Fruit Set—Two general linear models were constructed to examine wild-bee activity density in conjunction with honey-bee activity-density measures on the influence of percentage late fruit set and fruit set CV. For the first set of models, all wild bees were combined into a total wild-bee (WB) activity-density measure (“other bees” + bumble bees) and, in the second set of models, bumble bees and “other bees” as separate bee groups were used. These models were constructed both with and without the estimated honey-bee activity density of 2010.

Predicting Crop Yield—Using a general linear model, we constructed a predictive relationship of yield (kg berries/ha) by assessing the potential predictors of late fruit set, pest management, and field size (ha) as fixed effects and year and its interaction with fixed effects as random effects. We did not use any of the bee activity-density measures to model yield as we previously showed that bee activity densities were significant predictors of late fruit set. Fields were categorized into two types of pest management. Management category 1 consisted of those fields not treated with pesticides during the spring (May–June) of the study year (n = 16) and management category 2 included those that were treated with pesticides during the five months (May–September) of the study year (n = 24). Both pest management category fields were included in all study years.

Predicting Value of Bees—To determine the effect of field size on honey-bee and bumble-bee monetary value, we used a Wilcoxon nonparametric test because of heterogeneity and non-normality of the data.

3. Results

3.1. Honey-Bee Activity Density and Hive Stocking Rate

The different honey-bee counting protocols of 2011 and 2012 did not reveal significantly different mean values of bee activity density among years (t_(20.1) = 1.150, p = 0.264). When considering 2011 and 2012 fields for these two years, we found a positive relationship between numbers of honey bee hives in fields and honey-bee activity density with 27% of the variation in honey bee foraging explained by the model (F_(1,26) = 9.532, p = 0.005, r² = 0.27). The relationship shown in Figure 1 is

HB/min = 0.750 + 0.214 × HB hives/ha.

(3)

Similarly, we found that, in 2011 and 2012, the proportion of foraging honey bees and the number of honey-bee hives in the field exhibited a positive relationship with 26% of the variation in proportion of honey bees explained by the model (F_(1,26) = 8.942, p = 0.006, r² = 0.26). This linear relationship is

Proportion HB = 0.519 + 0.036 × HB hives/ha.

(4)

3.2. Study-Year Variation in Bee Activity Density

Figure 2 illustrates the distribution of the activity density for all three bee groups for all 40 fields for the three years of the study. Mean bumble bee (BB) activity density differed significantly between years (F_(2,37) = 25.248, p < 0.001). Activity density of BB declined in 2012, with a mean of 0.05 BB sightings per minute, compared to a mean of 0.27 bees/minute in both 2010 and 2011. Likewise, OB (other bees) activity density differed between years, but with an opposite trend; OB activity density significantly increased in 2012 compared to 2010 and 2011 (F_(2,37) = 11.66, p < 0.001). Honey-bee (HB) activity density did not vary between the years 2011 and 2012, the only years with HB activity-density measures recorded. Across all three years, BB activity density was weakly negatively correlated to OB activity density (Pearson’s r = −0.295, p = 0.066), but independent of HB activity. Activity density of OB was weakly positively correlated with HB activity density at the α = 0.10 significance level (Pearson’s r = 0.322, p = 0.095). There were no indications of a bee activity density by year interaction (p > 0.05).

3.3. Species Richness and Bee Activity Density

Total bee species richness averaged 29.3 ± 1.2/field (range = 17–50 species). Total bee species richness was positively correlated with OB activity density (r = 0.572, p < 0.001) and HB activity density (r = 0.398, p < 0.036). Total bee species richness was negatively correlated with BB activity density (r = −0.474, p < 0.002).

3.4. Field Size

Seventeen of the 18 large fields were stocked with honey bees, and only one field relied upon background (“other bees” + bumble bees) wild pollinators. The pollination of smaller fields was evenly distributed: six fields with bumble bees, six fields relying on background bees, and seven fields with honey bees. When considering field size as a continuous measure and not as a categorical variable, then bumble-bee stocked fields and background-bee fields were of similar sizes (mean size = 2.08 and 2.34 ha, respectively), which differed significantly from the honey-bee stocked fields (mean size = 6.68 ha) (F_(2,34) = 8.041, p = 0.001). Yield (kg berries/ha) was different for each size category. Smaller fields (n = 18) had lower yields than larger fields (F_(1,32) = 4.669, p = 0.038, r² = 0.13). Mean yield for small fields was 1914.3 kg/ha, while, for larger fields, it was 3201.0 kg/ha. Therefore, larger fields have a greater likelihood of honey-bee use for pollination and higher yield.

Total wild-bee (WB) activity-density measures (BB + OB) did not vary with field size class, and neither did BB or OB considered separately (WB: F_(1,38) = 0.155, p = 0.696). Activity density of HB, however, was significantly lower in the smaller fields (F_(1,26) = 10.871, p = 0.003, r² = 0.29), when only considering the 28 fields of 2011 and 2012. This relationship did not change when the estimated honey-bee activity measures of 2010 were included (F_(1,37) = 15.499, p < 0.001, r² = 0.30).

3.5. Fruit Set

Fruit set (all measures) did not vary according to field size (early fruit set: F_(1,38) = 0.612, p = 0.439; late fruit set: F_(1,38) = 0.416, p = 0.523; CV late fruit set: F_(1,38) = 0.551, p = 0.463). Early fruit set and the CV of late fruit set were determined by year (early fruit set: F_(2,37) = 1.32, p = 0.002, r² = 0.28; CV fruit set: F_(2,37) = 3.28, p = 0.049, r² = 0.15). In 2012, early fruit set was lower than the two previous years, while the CV of late fruit set was higher in this same year. Late fruit set followed the same annual pattern as did early fruit set (lowest value in 2012), but the year effect was not statistically supported at the α = 0.05 level (F_(2,37) = 2.356, p = 0.109, r² = 0.11). Late fruit set for fields relying solely on wild bees for pollination was 43.8%. Late fruit set for all fields (n = 40) was 48.3%. Early fruit set was highly correlated to late fruit set (r = 0.669, p < 0.001). Therefore, we only looked at late fruit set and fruit set CV as measures of pollination to interpret the effectiveness of the three bee groups that we observed (HB, BB, and OB).

The summary of analyses in

Table 1. Predictive models for percentage late fruit set and coefficient of variation (CV) of fruit set describing relationship between late fruit set and CV of fruit set, and bee activity density for honey bees, two wild-bee groups, total wild native bee activity density, and bee species richness. Activity density was measured as bees/minute; WB = total wild bees (bumblebees + “other bees”), BB = bumble bees), HB = honey bees.

Fruit Set Model	Predictors	Degrees of Freedom	F Ratio	Probability 1	Coefficient
Late fruit set HB, WB, Richness
Fixed effects Marginal r2: 0.336	HB WB Richness	1, 24 1, 24 1, 24	6.604 6.416 4.046	0.017 0.018 0.056	5.18 ± 2.02 28.56 ± 11.28 −0.75 ± 0.37
Random effects	Year Year × HB Year × WB Year × Richness			<0.001 0.164 0.133 0.376
Late fruit set HB, BB, OB, Richness
Fixed effects Marginal r2: 0.398	HB BB OB Richness	1, 23 1, 23 1, 23 1, 23	5.982 18.523 2.098 3.029	0.023 <0.001 0.161 0.095	5.14 ± 2.10 28.45 ± 6.61 29.37 ± 20.28 −0.77 ± 0.44
Random effects	Year Year × HB Year × BB Year × OB Year × Richness			<0.001 0.210 0.311 0.402 0.719
CV of fruit set HB, WB, Richness
Fixed effects Marginal r2: 0.203	HB WB Richness	1, 24 1, 24 1, 24	3.328 4.493 1.115	0.086 0.045 0.301	−6.25 ± 3.43 − 40.62 ± 19.16 0.67 ± 0.64
Random effects	Year Year × HB Year × WB Year × Richness			<0.001 0.095 0.115 0.299
CV of fruit set HB, BB, OB, Richness
Fixed effects Marginal r2: 0.278	HB BB OB Richness	1, 22 1, 22 1, 22 1, 22	159.366 2731.374 0.119 0.102	<0.001 <0.001 0.733 0.752	−8.15 ± 0.65 −52.26 ± 5.98 −31.22 ± 90.50 0.32 ± 0.35
Random effects	Year Year × HB Year × BB Year × OB Year × Richness			<0.001 0.177 0.226 0.326 0.564

¹ Probabilities for fixed effects are based upon Satterthwaite’s F statistic, and probabilities for random effects are based upon Wald’s statistic. Bold entries are significant at p < 0.05.

indicates that the random effect of year in the four predictive models for the two fruit set measures was significant, but none of the year × bee activity-density measures or bee richness were significant. The predictive model of late fruit set with HB, WB, and bee richness suggests that HB along with WB explained 33.6% of the variance in late fruit set. In this model, bee richness was marginally significant (α = 0.10), but with a negative coefficient. However, the model with all three bee groups (HB, BB, OB) and bee species richness was the best predictor of late fruit set explaining 39.8% of the total variance in late fruit set. In this model, only HB and BB were significant, while bee species richness was again only marginally significant with a negative coefficient (α = 0.10).

The bee group OB was not a good predictor of late fruit set except when combined with BB as WB. The poor performance of OB as a predictor may be due to its correlation with bee species richness and BB; thus, much of the variance explained by OB is already accounted for by BB. A similar pattern was observed with CV of fruit set. The mode with HB, BB, OB, and bee species richness was the best predictive model, explaining 27.8% of the variance in the CV of fruit set. In this model, only HB and BB were significant and both had negative coefficients, suggesting that the presence of these two bee taxa reduced the spatial variation in fruit set across a field, providing a more uniform potential yield.

3.6. Yield and Crop Value

The two measures of crop success, yield (kg berries/ha) and the crop value per hectare, were strongly correlated (r = 0.826, p < 0.001) as expected. Additionally, both measures were affected by year (crop value: F_(2,31) = 4.32, p = 0.022, r² = 0.22; kg/ha: F_(2,31) = 4.725, p = 0.016, r² = 0.23). The year 2012 was characterized by the greatest yield (3522 kg/ha), followed by 2010 (2189 kg/ha) and 2011 (1571 kg/ha). Similarly, 2012 experienced the greatest value per hectare even though, in that year, the price of berries was not the highest of the three years [40].

Table 2. Mixed model results for prediction of yield (kg berries/ha). Yield is a farmer reported value, based on 34 fields, 2010–2012.

Fruit Set Model	Predictors	Degrees of Freedom	F Ratio	Probability 1	Coefficient
Yield (kg/ha) Late fruit set, field size (ha), field management
Fixed effects Marginal r2: 0.238	Late fruit set Field size Management	1, 22040 1, 100771 1, 385689	6.555 4.023 0.887	0.011 0.049 0.346	48.72 ± 19.03 −549.46 ± 198.08 −301.89 ± 320.56
Random effects	Year Year × Late fruit set Year × Field size (ha) Year × Management			<0.001 0.764 0.589 0.677

¹ Probabilities for fixed effects are based upon Satterthwaite’s F statistic and probabilities for random effects are based upon Wald’s statistic. Bold entries are significant at P ≤ 0.05.

summarizes the predictive model of yield (kg berries/ha). These results indicate that late fruit set and field size were significant predictors of yield. Field pest management was not a significant predictor. The variance in yield explained by pest management may have been masked by the annual variation in yield and would have been accounted for by the random year effect. No indication of an interaction between late fruit set, field size, or management and year was detected (p > 0.05). The CV of fruit set was not included as a predictor because it was negatively correlated with late fruit set (r = −0.911, p < 0.001). Although, if the CV of fruit set was used as a predictor variable, we found it to be a significant predictor of yield (F_(1,12410) = 8.577, p = 0.003) and, thus, crop value. The negative coefficient (−35.49 ± 12.11) indicates that, as the clonal variation in fruit set decreases, the total field yield increases. Because crop value is an annual price multiplier of yield, we did not attempt to develop a predictive model for crop value.

3.7. Economic Value of Wild Bees and Honey Bees

Following the modified Losey and Vaughan [35] model defined in Equation (2), the value of a bee group is related to the proportion in which that group is found in the field. As shown above, the activity density of all bees (commercial and wild) predicts fruit set and fruit set predicts yield and, thus, crop value. We determined that the mean proportion of total wild bees (WB) found in all 40 study fields was 0.36. With application of the model, this suggests that WB, as a large group encompassing multiple taxa, are on average related to up to 36% of the value of the wild-blueberry crop.

Table 3. Mean value of the pollination services provided by three bee groups when (A) each group is considered to be an equally efficient pollinator of wild blueberry and (B) the activity density of wild bees is adjusted for efficiency in pollination using the pollination rate equivalency developed by Javorek et al. [7]. HB = honeybees, BB = bumblebees, OB = other wild bees (i.e., all non-Bombus) bees; n = 34 fields.

Bee Group	Mean Value of Pollination (USD)	Mean Value of Pollination (Percentage of Gross in USD)
A: bee groups considered as equivalent pollinators
HB	$19,370.55	64%
BB	$2033.85	14%
OB	$4239.15	22%
B: bee group activity density adjusted by pollination equivalencies
HB	$15,024.61	43%
BB	$5614.49	33%
OB	$5016.96	24%

reports the mean value of pollination and the percentage of the value of the crop attributable to pollination for the three target bee groups as determined by application of the Losey and Vaughan model [35]. The values given reflect two conditions: when the bee groups are considered equally effective pollinators and when the pollinating power is adjusted by pollination equivalencies as determined by Javorek et al. [7]. The OB group includes the genus Andrena whose abundance was adjusted by a factor of 3.6 (see Section 2.6).

3.8. Pollination without Rental Honeybees

Using the data from 34 of our study fields, we built a model that describes the relationship between yield (kg/berries/ha) and HB hives/ha (F_(1,32) = 4.310, p = 0.046, r² = 0.12).

kg/ha = 1754.85 + 216.39 × HB hives/ha.

(5)

Another model illustrates the relationship between the value of the crop and the density of HB stocking (F_(1,32) = 6.102, P = 0.019, r² = 0.16).

USD crop value/ha = 2759.12 + 417.80 × HB hives/ha.

(6)

Both of these models (see Figure 3A,B) can be used to estimate the effect of the absence of rental HB in wild-blueberry production, given that all other aspects of production and market forces remain the same. However, the intercepts can be seen as the strength of background pollination which is mostly due to wild bees (WB) or the fruit set in the absence of stocked HB. As the number of hives/acre drops to zero, the expected yield is 1754.9 kg berries/ha (95% confidence interval (CI) = 966.1) and the expected value attains USD 2759/ha (95% CI = USD 1567.6). This represents nearly a 30% decrease in comparison to the mean yield.

For all three years, one isolated field in this study had no honey-bee activity. Since this field was not harvested by the manager, we estimated yield by hand-raking berries. The estimated yield was 2222 kg/ha with an associated value of USD 3742/ha.

4. Discussion

4.1. Bee Abundance

In this study, we provide an analysis of the wild and managed bee community that influences fruit set, yield, and, ultimately, economic value of the wild-blueberry crop in Maine. We measured bumble-bee (BB) and other wild-bee (OB) activity in all 40 study fields and found that, on average, the WB abundance (OB + BB) made up 36% of the total bee abundance in fields, although, as seen in Figure 2, the abundance is not consistent across fields.

In our study, BB activity density dropped dramatically in 2012 to a mean of 0.05 bees/min from a mean of 0.27 bees/min recorded for the two previous years. We suspect this drop is due to unusually warm weather in March 2012 that triggered abnormally early queen foraging activity. By the time wild blueberry came into bloom, the temperatures were more normal, but bumble-bee queens were not commonly found in fields. Our study revealed six fields that registered 0.0 BB/min, all of which occurred in 2012. While the sample size is too small for significance tests of a reliable power, late fruit set for these six fields was among the lowest recorded in this study.

The OB activity density significantly increased during 2012, although the overall activity of OB was only weakly negatively associated with BB activity. We do not know why the OB abundance increased in 2012, but this may have been an artefact of sampling methodology. The activity-density measure derives from active bee-hunting with visual scans on flowers. It is possible that the lower BB abundance changed the visual environment such that the OB were easier to detect. However, taken together, the activity-density measure indicates that wild-blueberry fields have a persistent wild-bee (WB) community that varies over time and space.

4.2. Bee Abundance and Crop Fruit Set, Yield, and Crop Value

This persistent abundance of wild bees positively relates to late fruit set, which is a key predictor of the harvestable berry measure of yield (Table 2). In our strongest model, 39.8% of the variance in late fruit set was due to HB, BB, and OB activity, as well as species richness (Table 1). Honey-bee and BB abundance were strongly significantly associated with increased fruit set and negatively associated with the CV of fruit set, which suggests that these bee groups reduce spatial variation in fruit set across the field. Other wild-bee abundance (OB) and species richness were not significantly related to fruit set. One might expect OB abundance to have a positive relationship with fruit set, as many of the bees we collected were known pollinators of wild blueberry and were often carrying pure loads of ericaceous pollen [5]. The lack of a relationship could be due to the fact that this group includes all wild bees, including many (e.g., Sphecodes) suspected as being weak pollinators at best [5]. Furthermore, while most of the bees we caught were actively foraging on wild-blueberry flowers, representatives of other flower taxa located within the fields were examined for bees and those bees were a component of the OB activity-density measure.

Honey-bee and BB relationships with fruit set were shown in other studies. In a two-year study, Eaton and Murray [41] documented a significant relationship between both HB and BB activity densities and fruit set, with a strong relationship for HB during both years of the study and a relationship for BB during one year. Unlike our study, they also found a strong correlation between OB activity densities and yield in one year of the study, indicating that non-Bombus wild bees can significantly contribute to the value of a crop, but not consistently. Assessing fruit set during a wild-blueberry study in fields stocked either with commercial BB or with HB, Drummond [17] found a positive relationship between BB foraging density and fruit set. Eaton and Nams [42] showed yield and fruit set to increase with greater activity densities of honey bees and fruit-set measures to improve with increasing densities of non-Apis wild bees. However, unlike our study, their research indicated an interaction between HB and OB, such that evidence of competition for floral resources emerged at high honey-bee activity densities. Considering wild bees in general, in a long-term study, Asare et al. [18] assessed HB and WB activity densities on 162 wild-blueberry fields from 1993–2015 and found that both HB and WB activity densities were related to fruit set. They found that fruit set was impacted on a per bee basis 1.6 times more by wild bees than honey bees. Similar to Eaton and Nams [42], a negative interaction was observed between HB and WB activity densities on fruit set (but at p = 0.089), also indicating a potential competition effect for floral resources between the two bee groups.

One field from our study illustrates the power of wild pollinators as this field was not stocked with any commercial bees and was isolated from human activity such that we did not record any HB at any time in the field. Final fruit set for this field was 42%, which we can assume was due primarily to WB pollination, although other pollinating insect taxa cannot be discounted [38]. Pollination research at the University of Maine in the 1950s suggests a similar level of fruit set due to WB visits [43]. Researchers at this time were interested in determining the efficacy of HB as pollinators and compared fruit set in cages with trapped HB to the fruit set in adjacent open areas without HB but with WB visitors. At two sites in 1959, the research documented an average of 30% fruit set outside the HB cages. Desjardins and De Oliveira [44] documented a fruit set of 42% for wild blueberry due to native, “background” pollinators. Drummond [17] found a fruit set of 28% for wild-blueberry fields in Maine without commercial pollinators, and Stubbs and Drummond [14] found the fruit set due to wild bees in 10 wild-blueberry fields in the absence of commercial pollinators to be 55.5% ± 6.5% (standard error, SE). Asare et al. [18] found that fields with only wild bees had an average fruit set of 33.7% and, for every unit (bees/m²/min) increase of wild bee density, fruit set increased by 11%. Similarly, we recorded the mean final fruit set for fields relying on background pollination (some of which may be provided by honey bees placed on nearby farms) to be 43.8%. Taken together, our data along with previous research indicate that wild-bee pollination achieves 30–55% fruit set in wild blueberry without the addition of commercial pollinators. It must be acknowledged, however, that other factors including environmental conditions [40,45], genotype [4], bloom phenology [46], and possibly nectar production [47] also affect fruit set. Asare et al. [18] presented a conceptual model for the multitude of factors that affect fruit set and yield including stem density, flowers/stem, bee species community diversity, bee species efficiency, bee density, weather conditions, clone genotype, frost, fertility, pruning, weeds, disease, and insect pests.

Yield and crop value are also related to factors beyond pollination. In wild-blueberry production, inputs from farmers include, among others, resources needed for pruning, fertilizing, pest monitoring, pest control, and harvesting, in addition to pollination purchases. For a farmer seeking an average yield, complete costs (fixed and variable) are estimated at USD 3370/ha [48], and those inputs are not trivial compared to total crop value, which in our study ranged from about USD 1000 to USD 10,000/ ha. However, given that no berries will result without pollinators [9], the value of these non-pollination inputs can be seen as a “risk avoidance” investment to enhance the likelihood of developing fruit becoming an economically viable fruit.

Therefore, our study suggests that farmers might expect increased fruit set and yield and, thus, crop value, with increasing abundance of wild bees generally and bumble bees specifically. Applying the Losey and Vaughan [35] model to our data suggests that Maine farmers can expect up to 36% of crop value to be due to WB activity density if all pollinators are considered equivalent in terms of power of pollination. However, research indicates that many species of WB and especially BB are far more efficient pollinators (on a per bee basis) than HB [7,8,17,31,44]. Therefore, the value of wild pollinators, according to our modified Losey and Vaughan [35] model, suggests WB value may approach 57% of the crop even when rental HBs are used in the field (Table 3).

4.3. Bee Abundance and Field Size

Neither BB nor WB activity density varied by field size. Given that none of our study fields exceeded 17 hectares, we cannot speak directly on the role of WB in fields larger than those we examined. Our study fields were small compared to the “blueberry barrens” of Downeast Maine where many fields are larger than 200 ha. However, more than 75% of farmers harvest less than 20 ha per year, suggesting that our data represent field sizes typical of the blueberry agroecosystem in Maine [10]. Given the size of Maine BB queens and workers measured from our collections [5], we expect Bombus ternarius Say queens (the most common species) to forage up to 5.8 km from the nest site and B. ternarius workers to travel up to 0.9 km [49,50]. We would, therefore, expect the proportions to diminish only with field sizes approaching double that of our study’s largest field size. Considering that foraging bumble-bee queens (which are predominant during blueberry bloom) are also engaged in nest site selection, we cannot say that queen foraging is constrained by distance to the nest.

Honey-bee activity density decreased significantly with decreasing field size, due to fewer hives/hectare for the smaller fields. Why do farmers of fields the size of our large fields increasingly rely upon HB when we recorded no differences in WB or BB activity density within larger and smaller fields, and given the superior pollinating abilities of WB or BB? We first suggest that different survey methods could reveal a greater WB abundance along the field edges as reported in other studies [51,52] even in fields limited to those of the size we sampled [53,54]. If this were the case, increased reliance upon imported HB in larger fields approaching 20 ha would be justified. Additionally, HB may ensure pollination at field interior regions, especially if the hives are not placed along the field edge. Aras et al. [55] showed that HBs initially forage close to where the hives are placed in wild-blueberry fields, but quickly expand their foraging distance until they move across even the largest of fields. In wild-blueberry fields in Canada, Cutler et al. [52] found that most wild bees decreased in abundance as the distance from the field edge approached 200 m; however, larger-bodied bees were found further from the field edges than the smaller bees. Therefore, BB [17,44] and WB such as Andrena spp. that nest within blueberry fields [56] may be responsible for much of the pollination in a field interior. Validation of this hypothesis that HBs play a role in increasing pollination toward field interiors would suggest that farmers wishing to decrease reliance on imported bees could use wild-bee habitat plantings [57] and adjust crop rotation cycles so that larger fields function effectively as smaller fields [1].

4.4. Conclusions

Our estimates of the percentage of crop value attributable to a given bee group is analogous to a measure of “vulnerability” [58] and indicates vulnerability for both the decline of WB and the lack of HB. It is not unrealistic to consider what a pollination regime would look like given a lack of HB or a decline in the abundance of WB [23,59,60]. Our results indicate that, all else being equal, a farmer could expect a little over 1755 kg/ha without HB (Figure 3A), although this value can also be considered as the strength of WB pollination. If WB and HB pollination services are complementary and not interchangeable [33], then declines in wild pollinators may result in farmers losing 1755 kg/ha of berries that cannot be recovered with increasing HB colonies. If the pollination services of wild bees can be compensated for by honey bees, then the pollination costs to farmers will increase with a loss of WB populations. For each additional HB hive/ha, a farmer can expect an additional 1971.39 kg of berries/ha (Equation (5)), which, at a current value of USD 100–150/hive (F.A. Drummond, pers. com.), will increase costs for larger fields substantially.

We applied our analysis to wild blueberry in Maine. Our goal was to illustrate how WB pollinators, especially BB, are related to increased wild-blueberry yield through increasing fruit set. Our hope is that a farmer who understands the economic value of WB may be motivated to alter management regimes in order to conserve this natural resource, especially if the bee’s presence in the field provides a sufficient economic value. Following this reasoning, it could be argued that low economic value of a wild resource or changing economic values associated with land use may motivate a farmer to not conserve a resource [61,62]. However, given the dependency of most fruits and vegetables on pollination [61,63], we feel that sufficient evidence exists to conclude that significant levels of pollination are performed by WB [33], which would, therefore, provide a financial boon for farmers in the form of a highly valuable ecosystem service. This finding also provides a strong justification for bee and bee-habitat conservation by all wild-blueberry farmers.