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).