Statistical analyses
We used linear mixed models to analyse effects of ozone, nitrogen enrichment (i.e., including both the mean values of N fertilizer application on the agricultural fields and the NOx concentration in the air from satellite data), the risk of pesticide exposure and the proportion of (semi-) natural habitats and their two-way interactions on the abundance of pollinators and their contribution to crop production (see correlation matrix in Appendix S2).
The local abundance of honey bees is primarily determine by beekeeper behaviour rather than local effects of habitats (Büchler et al.2014; Wood et al. 2020). As managed species they are influenced differently by environmental pressures compared to wild pollinators, and we therefore analysed Apis mellifera separately from non-Apis pollinators (i.e., other bees and hoverflies).
To account for variation associated with the crop system on pollinators and pollination, crop identity was included as random effect in all models. Moreover, to remove potential confounding effects with study region or country, all explanatory variables included in each model were centered within study-year combinations (Van de Pol & Wright 2009).
As previous studies have also shown that densities of non-Apispollinators can in some circumstances be negatively affected by honey bee densities (e.g. Lindström et al., 2016; Geslin et al., 2017; Mallinger et al., 2017), honey bee abundance was included as explanatory variable in non-Apis pollinator models. For the analysis of the contribution of pollinators to crop production, in addition to sources of eutrophication, ozone pollution, pesticide risk and proportion of (semi-)natural habitats, we included abundance of honey bees (Apis mellifera ) and non-Apis pollinators as covariates.
First, to test for spatial autocorrelation, we compared models with different spatial correlation structure (exponential, Gaussian, Linear, rational quadratics and spherical spatial autocorrelation) and without spatial correlation structure, and defined the best random structure of the model based on their AICc scores (Akaike Information Criterion for small samples). Then, we applied model selection to the fixed terms of the model (∆AICc < 2 with the best model; Anderson et al., 2001). To not overfit the global model in relation to our sample size, the number of parameters in each tested model was restricted to 5 (including potential interaction effects). Selection of the best candidate models are presented in Supplementary material (see Appendix S3, S4 and S5 in Supporting Information). All analyses were computed using the ape (Paradis et al. 2019), nlme (Pinheiro et al.2020) and MuMIn (Bartoń 2011) packages in R software, version 3.4.2 (R Development Core Team 2018). All spatial extraction or landscape index calculation from shapefile and raster maps were made using QGIS software version 3.10 A Coruña (QGIS Development Team 2020).