Wild pollinators and honeybees respond differently to landscape-scale organic farming and increase sunflower yields

Wild pollinators play a critical role in crop production, yet they are increasingly threatened by agricultural intensification and habitat loss. Hence, identifying effective measures to support pollinators at landscape and field scale is crucial for maintaining pollination services and ensuring sustainable food production. We assessed how landscape composition (area of organic farming, semi-natural habitats and mass-flowering crops) and field management (farming system, weed cover and weed richness) influence wild pollinators and honeybees in sunflower fields. Additionally, we used a pollinator exclusion experiment to assess the effects of landscape composition, field management and pollinators on seed weight, seed number, pollination services and overall yield. Bumblebee abundance increased with organic farming area in the landscape, while solitary bee richness increased with semi-natural habitat area. Both bumblebees and hoverflies declined in abundance with increasing mass-flowering crop area in the landscape. At field level, the abundance and richness of solitary bees and hoverflies increased with weed richness. Insect pollination in open compared to pollinator-excluded treatments increased yields on average by 25%. Pollination services and overall yields were not affected by weeds. Overall yields did not differ between conventional and organic fields, while pollination services were marginally higher in organic fields. Our findings underscore the need for multi-scale conservation strategies to sustain pollinators and pollination services. Increasing organic farming at the landscape scale can support pollinators across both organic and conventional systems but cannot replace semi-natural habitats, which remain essential to enhance solitary bees in crop fields. Landscape management should therefore promote both organic farming and semi-natural habitats. Tolerating moderate weed levels within fields can further enhance wild pollinators without reducing yields. Farmers should also consider the amount of simultaneously mass-flowering crops in the landscape, to avoid dilution effects. Our findings provide practical strategies to support different groups of wild pollinators through integrative landscape and field management and strengthen pollination services in agroecosystems.

Methods
Site selection

The study was carried out in 2022 on 14 conventional and 15 organic sunflower fields in agricultural landscapes around Würzburg (Northern Bavaria, Germany, Fig. 1A). The fields were selected in such way that the surrounding landscape composition covered the broadest possible gradients in the area of 1) annual organic farming, 2) semi-natural habitats and 3) sunflower fields, whereby maximizing gradient 1 had highest and gradient 3 had lowest priority (Fig. 1B). Conventional farmers applied mineral fertilizer and herbicides, but no insecticides in the sunflower fields. Organic farmers applied no pesticides or fertilizers, with two exceptions where manure or compost were used. The average number of different crops in the typical rotation was 2.75 ± 0.48 (mean ± SE) for conventional and 4.3 ± 0.33 (mean ± SE) for organic sunflower fields. Minimum distance between study fields was 5.8 ± 0.75 km (mean ± SE). In two cases, distance between fields was less than 2 km (958 and 1671 m). Field size ranged between 0.8 and 39.9 ha and was on average 5.6 ± 1.3 ha (mean ± SE). Neither field size nor sunflower coverage differed between conventional and organic sunflower fields (t-test, all P > 0.1).

Landscape composition

The study region is dominated by annual crop fields, interspersed with semi-natural habitat fragments, villages and forests. The most dominant crop in the study region is winter wheat, with other typical crops being barley, spelt and oilseed rape (Supporting Information, Fig. S1). We calculated landscape composition around each field in a 500 m and 1 km radius. These radiuses were chosen to cover the typical foraging distance of the studied groups of pollinators (Greenleaf et al., 2007; Kendall et al., 2022; Meyer et al., 2009; Zurbuchen et al., 2010). Using ArcGIS Pro version 3.2.1 (ESRI, 2023), we calculated the proportion of annual organic crop fields and sunflower fields in the study year (hereafter “organic farming area” and “sunflower area”), using data from the Integrated Administration and Control System (Federal Ministry of Food and Agriculture, 2022). To calculate proportions of semi-natural habitats (“semi-natural habitat area”), we created a map combining data from different sources (Supporting information) which was validated and updated via landscape surveys in the study year. The following habitats were characterized as semi-natural habitats: orchard meadows, extensively managed grasslands, fallows, reed, hedges and forest edges. Landscape variables were not correlated, with the exception of organic farming area and sunflower area, which were positively correlated (R = 0.54, Fig. S2), because sunflower is more often part of organic than of conventional crop rotations. We included both, organic farming area and sunflower area, as variables in the same models since they were expected to have opposite effects on pollinators.

Pollinator, flower and yield surveys

We sampled bees and hoverflies in two rounds in July 2022 during sunflower bloom. We differentiated between honeybees, bumblebees, all other wild bees (hereafter referred to as “solitary bees”) and hoverflies since these groups differ considerably in their ecology. Pollinators were surveyed using transect walks, standardized by time and area. Each transect consisted of 30 minutes searching time and covered 300 m2 (300 m length, 1 m width). One half of the transect (15 minutes, 150 m2) was walked along the edge of the field, the other half was walked in the center of the field, along the sowing row lines. Individuals were identified in the field or caught using a net for identification in the laboratory. For each flower-visiting pollinator recorded during the transect walks, we recorded whether it visited a sunflower or a weed. Additionally, we conducted flower surveys of weeds and sunflowers on the same transects in each round. Weed richness was assessed by counting the number of distinct non-sunflower morphospecies per transect. Weed and sunflower coverage were measured by the size (i.e. area) of one representative flower or floral unit per species and multiplying it by the number of flowers or floral units in the transect. If flower size varied within a transect, we calculated flower coverage for each section based on the respective representative flower size and summed them up for total coverage. We summed flower coverage across all weed species per transect. We also measured sun and shade temperature during each round. Pollinator surveys were always conducted between 9:00 am and 5:30 pm, when shade temperatures were at least 18° C and wind speeds did not exceed 3 Bft. The order (i.e. time of the day) in which sunflower fields were visited was randomized between rounds to avoid effects of sampling time. Pollinator abundances were not affected by sampling temperature (Table S1, Fig. S3). To assess the contribution of insect pollination to sunflower yield components, we conducted pollinator exclusion experiments on 28 of our study fields (14 conventional, 14 organic). In mid-June (i.e. before flowering), we marked eight sunflower heads of different plants along the edge and in the center of each field. All plants had a minimum distance of 2 m to each other and plants in the center of the field had a minimum distance of 20 m to all edges. Half the heads per location were bagged using a fine mesh net before bud opening, to allow wind and self-pollination only. After sunflower bloom, we removed the nets. Shortly before harvest, we hand-harvested the marked flowers. Seeds were extracted from the heads; empty seeds were sorted out and the head diameter was measured. Then, the seeds were dried at 55° C for 48 hours. Seeds per head were weighed using a fine scale and counted using a counting machine. Seed number and seed weight per plant were averaged for each field and location (center vs. edge). We calculated the yield of open and bagged flowers as average seed weight per flower per location, multiplied by flower density per location. Pollination services were calculated as the difference in yield between open and bagged flowers per location and field. To estimate overall yield, we used the yield calculated from open pollinated flowers in the field center.

Statistical analyses

We summed pollinator abundances per group across rounds and calculated cumulative species richness for pollinators and weeds per field. Weed and sunflower coverage were averaged per field across rounds. We calculated Chao1 estimators for species diversity for each wild pollinator group using the ‘vegan’ package (Oksanen et al., 2022). Since the Chao1 estimator was highly correlated with measured species richness for all pollinator groups (r > 0.7, P < 0.001, Fig. S4), we used measured species richness in all our models.

We used generalized linear models (GLMs) to identify effects of landscape composition and field-scale management on pollinators, weeds, pollination services and overall yield. For pollinator abundance models, we used negative binomial distributions, to account for overdispersion (Venables & Ripley, 2002). We used Gamma distributions for pollinator richness and Poisson distribution for solitary bee richness models. For weed, pollination services and overall yield models, we used Gaussian distributions. We fitted separate models for abundance and richness of each pollinator group, weed coverage (log-transformed), weed richness, pollination services, and overall yield as responses.

We implemented a multi-step modelling approach, because field-scale predictors were partly moderated by landscape-scale predictors, which could obscure the actual relationships between predictor and response variables (Arif & MacNeil, 2023). In the first model set, we included landscape-composition variables (organic farming area, semi-natural habitat area and sunflower area within a 1 km radius), farming system (conventional vs. organic), as well as the interactions of each landscape variable with farming system. For solitary bees, we ran an additional model set using landscape variables at a 500 m radius, because foraging distances in this pollinator group can be smaller than 1 km (Zurbuchen et al., 2010).

To assess the effects of field-scale management on pollinators, we calculated a second set of models with the predictors sunflower coverage, weed coverage, weed richness, and honeybee abundance as predictors. To identify effects of weeds on pollination services and overall yield, we calculated models with weed coverage and weed richness as predictors. In all models from the second set, we included all predictors from the first set of landscape models as control variables to obtain unbiased causal estimates for the field-scale variables (Arif & MacNeil, 2023). In the models for pollination services we included field ID as random effect.

To assess pollinator resource use (sunflower vs. weeds), we summed the number of pollinators recorded on sunflowers and the number of pollinators recorded on weeds per pollinator group per field. We then calculated negative binomial models using flower resource (i.e. sunflower vs. weeds), farming system and their interaction as predictors.

To test whether seed number, seed weight, and yield differ between bagged and open flowers,  we calculated GLMs (Brooks et al., 2017) with treatment (i.e. bagged vs. open), farming system and their interaction as fixed effects, and field ID as random effect. For models with seed weight or seed number, we used normal distribution and for models with yield, we used Gamma distribution. To evaluate the effects of pollinator abundance on pollination services, we calculated linear mixed effects models (Bates et al., 2015) with pollinator abundance as fixed effect, and field ID as random effect.

We assessed the significance of predictors for all models using a Type II ANOVA from the ‘car’ package, since our focus was to test the main effects of predictors, while still accounting for potential interactions (Fox & Weisberg, 2019). Residuals of all models were checked using ‘DHARMa’ (Hartig, 2022). Statistical analyses were done in R version 4.3.1 (R Core Team, 2023).