Soil eutrophication shaped the composition of pollinator assemblages during the past century

www.ecography.org ECOGRAPHY Ecography 209 –––––––––––––––––––––––––––––––––––––––– © 2019 The Authors. Ecography published by John Wiley & Sons Ltd on behalf of Nordic Society Oikos This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited. Subject Editor: Jian Zhang Editor-in-Chief: Miguel Araújo Accepted 27 September 2019 43: 209–221, 2020 doi: 10.1111/ecog.04656 doi: 10.1111/ecog.04656 43 209–221 Atmospheric nitrogen deposition and other sources of environmental eutrophication have increased substantially over the past century worldwide, notwithstanding the recent declining trends in Europe. Despite the recognized susceptibility of plants to eutrophication, few studies evaluated how impacts propagate to consumers, such as pollinators. Here we aim to test if soil eutrophication contributes to the temporal dynamics of pollinators and their larval resources. We used a temporally and spatially explicit historical dataset with information on species occurrences to test if soil eutrophication, and more specifically nitrogen deposition, contributes to the patterns of change of plant and pollinator richness in the Netherlands over an 80 yr period. We focus on bees and butterflies, two groups for which we have good knowledge of larval resources that allowed us to define groups of species with different nitrogen related diet preferences. For each group we estimated richness changes between different 20-yr periods at local, regional and national scale, using analytical methods developed for analyzing richness changes based on collection data. Our findings suggest that the impacts of soil eutrophication on plant communities propagate to higher trophic levels, but with a time-lag. Pollinators with nitrogenrelated diet preferences were particularly affected, in turn potentially impairing the performance of pollinator-dependent plants. Pollinator declines continued even after their focal plants started to recover. In addition, our results suggest that current levels Soil eutrophication shaped the composition of pollinator assemblages during the past century Luísa G.Carvalheiro, Jacobus C.Biesmeijer, MarkusFranzén, JesúsAguirre-Gutiérrez, Lucas A.Garibaldi, AveliinaHelm, DenisMichez, JuhaPöyry, MennoReemer, OliverSchweiger, Leon van den Berg, Michiel F.WallisDeVries and William E.Kunin L. G. Carvalheiro (https://orcid.org/0000-0001-7655-979X) ✉ ([email protected]), Depto de Ecologia, Univ. Federal de Goias, Goiania-GO, Brazil, and Centre for Ecology, Evolution and Environmental Changes (CE3C), Faculdade de Ciências da Univ. de Lisboa, Lisboa, Portugal. – J. C. Biesmeijer (https://orcid.org/0000-0003-0328-0573), Naturalis Biodiversity Center, Leiden, the Netherlands, and Inst. for Environmental Sciences (CML), Leiden Univ., Leiden, the Netherlands. – M. Franzén (https://orcid.org/0000-0001-8022-5004) and O. Schweiger, UFZ Centre for Environmental Research, Dept of Community Ecology, Halle, Germany. MF also at: Ecology and Evolution in Microbial Model Systems, EEMIS, Dept of Biology and Environmental Science, Linnaeus Univ., Kalmar, Sweden. – J. Aguirre-Gutiérrez (https://orcid.org/0000-0001-9190-3229), Environmental Change Inst., School of Geography and the Environment, Univ. of Oxford, Oxford, UK, and Naturalis Biodiversity Center, Leiden, the Netherlands. – L. A. Garibaldi, Inst. de Investigaciones en Recursos Naturales, Agroecología y Desarrollo Rural (IRNAD), Sede Andina, Univ. Nacional de Río Negro (UNRN) y Consejo Nacional de Investigaciones Científicas y Técnicas (CONICET), San Carlos de Bariloche, Río Negro, Argentina. – A. Helm, Inst. of Ecology and Earth Sciences, Univ. of Tartu, Tartu, Estonia. – D. Michez (https://orcid.org/0000-0001-8880-1838), Research Inst. of Biosciences, Laboratory of Zoology, Univ. of Mons, Mons, Belgium. – J. Pöyry (https://orcid.org/0000-0002-5548-1046), Finnish Environment Inst. (SYKE), Biodiversity Centre, Helsinki, Finland. – M. Reemer, European Invertebrate Survey – the Netherlands/Naturalis Biodiversity Center, Leiden, the Netherlands. – L. van den Berg, Bosgroep Zuid Nederland, Geldrop, the Netherlands, and Radboud Univ. Nijmegen, Nijmegen, the Netherlands. – M. F. WallisDeVries (https://orcid.org/0000-00033808-2999), De Vlinderstichting/Dutch Butterfly Conservation, Wageningen, the Netherlands, and Wageningen Univ., Plant Ecology and Nature Conservation Group, Wageningen, the Netherlands. – W. E. Kunin, School of Biology, Univ. of Leeds, Leeds, UK. Research 210 of nitrogen deposition still have a negative impact on most groups here analyzed, constraining richness recoveries and accentuating declines. Our results indicate that the global increase in nitrogen availability plays an important role in the ongoing pollinator decline. Consequently, species tolerances to soil nitrogen levels should be considered across all trophic levels in management plans that aim to halt biodiversity loss and enhance ecosystems services worldwide. Keywords: extinction debt, herbivory, historical biodiversity changes, nitrogen deposition, nitrophily, pollinator communities Introduction Ecosystems are undergoing rapid changes worldwide (IPBES 2019), resulting in biodiversity losses and biotic homogenization across many taxa (Dornelasetal. 2014), including plants and their pollinators (Biesmeijeretal. 2006, Carvalheiroetal. 2013). While there is strong evidence of overall declines in insect abundance and diversity worldwide (Hallmanet al. 2017, Lister and Garcia 2018), a previous study has shown that, in some regions of NW Europe, the negative trends in species richness become substantially lessened after 1990 (Carvalheiroetal. 2013). Land use and climate change partly explain these trends (Aguirre-Gutiérrezetal. 2015, 2016), but the contributions of other drivers, e.g. soil eutrophication, are still unclear. Nitrogen is essential for many biological processes (Elseretal. 2007), but due to combustion of fossil fuels and intensification of fertilizer use in agriculture and forestry (Gallowayetal. 2008, Suttonetal. 2011), excessive nitrogen availability is a current threat for biodiversity (Öckingeretal. 2006, WallisDeVries and Bobbink 2017). Impacts on plant community are well-recognized, including losses in species diversity (Stevensetal. 2004, Bobbinketal. 2010), and these may propagate to primary consumers (Stevensetal. 2018). While there has been recent research into effects on herbivores (Nijssenetal. 2017, Pöyryetal. 2017, WallisDeVries and van Swaay 2017), the potential effects of soil eutrophication on those feeding on pollen and nectar (potential pollinators) are largely unexplored (Stevensetal. 2018, but see Betzholtzetal. 2013, Tamburinietal. 2017, Ramosetal. 2018). Such potential bottom–up effects could be mediated by reduction in nectar or pollen availability or quality (Petanidouetal. 1999, Vanderplancketal. 2017). Also, insects may be deterred from visiting, or ovipositing, on plants whose chemical composition has been changed by nitrogen enrichment (Abbas et al. 2014, Audusseau et al. 2015, Kurzeetal. 2017). Even if visitation is not deterred, population fitness can be negatively impacted by the above mechanisms, potentially leading to local losses of diversity (as reviewed in Elliottetal. 2008). Recent implementation of environmental policies aiming to reduce deleterious impacts of agricultural practices in Europe has led to significant reductions of nitrous oxide emissions, which, combined with changes in industry strategies (e.g. offshoring), has contributed to an overall decrease of local atmospheric nitrogen deposition (EEA 2007). While leaching and denitrification can lead to rapid nitrogen loss from soil systems, eutrophication effects may still be detectable in soil and plant tissues’ chemical composition up to 10–20 yr (Dörretal. 2009, Kandeleretal. 2009). Past and ongoing changes in diversity may hence be mediated by species’ sensitivity to nitrogen. In this study, we take advantage of the detailed spatially explicit information on nitrogen deposition (< www.RIVM.nl >), and plant and pollinator richness changes (Carvalheiroetal. 2013) in the Netherlands to explore if nitrogen deposition and tolerance of species to nitrogen is influencing the patterns of richness change of plants and of insects that feed on them (bees and butterflies). Historical data frequently only allows to run distribution and species richness analyses (Biesmeijer et al. 2006, Carvalheiroetal. 2013), not being adequate for analyzing species abundance dynamics. When doing short-term evaluations focused solely on species richness patterns, negative effects on biodiversity may be missed because community change is slow to detect on the basis of presence–absence data. However, in the long term, if more species show a contracting spatial range (i.e. less presence datapoints) rather than range expansion, then negative effects in richness will be detected. In addition, changes in species spatial distribution regulate richness change patterns, with species’ ranges contractions or expansions leading to more accentuated richness changes at finer spatial scales (e.g. 10 × 10 km) than at broader scales (with no change at country level) (Thomas and Abery 1995). As high nitrogen availability favors plant species that are commonly associated with nitrogen-rich habitats (McClean et al. 2011), we predict that the range expansion of nitrophilous plants drove the previously detected (Carvalheiroetal. 2013) increase in local plant richness during periods of increased nitrogen deposition. Moreover, most plant species have multiple pollinator species to assure the levels of pollen deposition essential to their reproduction (Waseretal. 1996, Johnson and Steiner 2000). Therefore, we predict that the level of nitrophily of plant species will play a more important role in defining overall patterns of plant richness change than species’ pollinator dependence. The response of consumers is likely to depend on their ability to adjust their diet (i.e. on the level of specialization, Carvalheiroetal. 2010), potentially lagging behind plants. Consequently, we predict that during periods of increased dominance of nitrophilous plants, declines will be more accentuated for more specialized consumers particularly if their major plant resources are nitrophobous species. Finally, given that many landscapes in our study region are still subjected to high levels of nitrogen deposition (Fig. 1), we predict that recent nitrogen deposition will be negatively spatially related with richness changes 211 of nitrophobous plants and their associated insects, and positively related with richness changes of nitrophilous plants and their associated insects. Material and methods We evaluate if the patterns of richness changes (i.e. local and regional declines or increases) of Dutch plants, bees and butterflies over an 80-yr period depended on species’ resource requirements (plant preference for nitrogen rich or poor habitats, and insect diet specialization). Bee (Hymenoptera: Apoidea: Anthophila), butterfly (Lepidoptera: Papilionoidea) and vascular plant data sources are described in Supplementary material Appendix 1 Table A1. Species traits Plant species were classified into two pollinator dependence groups: pollinator-dependent and non-pollinator-dependent (abioticor self-pollinated). For each plant species we then obtained the Ellenberg N-value (Ellenbergetal. 1991), which refers to the species’ soil nitrogen preference, from 1 (low) to 9 (high). Following the classification used by Öckingeretal. (2006), plant species were then divided into two classes of nitrogen responses: nitrophobous (Ellenberg values from 1 to 5) and nitrophilous (Ellenberg values from 6 to 9). Information on the dependency of plants on insects for pollination was obtained from several databases (Supplementary material Appendix 1 Table A1). For additional analyses focusing on plants that are considered to be important bee flower resources, we considered all plants that are classified as ‘bee flower’, ‘bumblebee flower’, or as a transition type from bee or bumblebee flowers to any other groups of flower in the BIOFLOR database (Müller classification, Kühnetal. 2004). Extensive information on Dutch species’ diet allowed us to classify bee (Supplementary material Appendix 1 Table A2) and butterfly (Supplementary material Appendix 1 Table A3) species into four diet classes (Fig. 2). We first considered three classes related to the level of specialization of their larval diet (i.e. host plants for butterflies; and pollen sources for bees) on plants with a specific N requirement: 1) species whose major larval diet sources (i.e. those making more than 75% of the diet) consist of nitrophilous plants; 2) species whose major larval diet sources consist of nitrophobous plants; and 3) nitrogen generalists (i.e. species without indication of strong preference for nitrophilous or non-nitrophilous plants). Figure 1. Nitrogen deposition changes through time in the Netherlands. Mean (± standard deviation) nitrogen deposition (kg N ha−1 yr−1) per 10 × 10 km grid cell of the Netherlands 20 yr–1 period since 1930 until 2009 are presented. Grey horizontal bar represents the nitrogen critical load range (quantitative estimate of an exposure to one or more pollutants below which significant harmful effects on specified sensitive elements of the environment do not occur according to present knowledge) determined in van Dobbenetal. (2013) for each of the Dutch habitats classified as suitable for pollinators (i.e. all grasslands, moors/peat, forests and sandy soils, Vogiatzakisetal. 2015, Aguirre-Gutierrezetal. 2016). Overall there were considerable changes in nitrogen deposition across the four periods compared in this study. Nitrogen deposition increased between the first (P0: 1930–1949), second (P1: 1950– 1969) and the third (P2: 1970–1989) period, reaching its peak in P3. Between the third and fourth period (P4: 1990–2009) nitrogen deposition slightly decreased. Figure2. Plants and pollinators can be grouped according to their nitrogen related resource preferences. Plant species with Ellenberg values from 1 to 5 were considered nitrophobous; plant species with Ellenberg values from 6 to 9 were considered nitrophilous. Pollinators were divided in four diet related classes 1) species whose major larval diet (i.e. those making more than 75% of the diet) consists of nitrophilous plants; 2) species whose major larval diet consists of nitrophobous plants; 3) polyphagous nitrogen generalists and 4) non polyphagous nitrogen generalists. Information on insect larval diet was obtained from several data sources for butterflies (Bink 1992, Stoltze 1996, Eliasson et al. 2005, van Swaay 2006, Dennis 2012), and bees (Westrich 1990, Peetersetal. 1999, 2004, Peeters and Reemer 2003, Petterssonetal. 2004, Müller and Kuhlmann 2008, Davisetal. 2011, Scheperetal. 2014). Additional information on diet specialization was obtained from the BWARS trait database (the Bees, Wasps and Ants Recording Society, < www. bwars.com >). For parasitic bees we used information on the diet of their main host(s) in the Netherlands. 212 Nitrogen generalists were further divided into polyphagous (i.e. having less than 90% of records on the same plant family, Müller and Kuhlmann 2008); and non-polyphagous species (i.e. feeding on only one family of plants, including mono and oligophagous species). For the groups of nitrogen specialists (i.e. species that show a preference towards plants with a given Ellenberg N value) the number of species was not sufficient to analyze separately polyphagous and nonpolyphagous species: polyphagous species make up 13% of the bees preferring nitrophilous plants, 19% of the bees preferring nitrophobous plants, 40% of the butterflies preferring nitrophilous plants and 5% of the butterflies preferring nitrophobous plants. For details on number of species and records per group and time period see Supplementary material Appendix 1 Table A1, and for species list per each class see Supplementary material Appendix 1 Table A2, A3. Only native species were considered for richness change calculations (i.e. for plants all neophytes were excluded, for bees and butterflies all species present in a region as an intentional or accidental result of human activity were excluded, see Richardsonetal. 2000). For plants, we focused on herbs to minimize the influence of range changes caused by recent restoration programs in the Netherlands, which involved the plantation of many shrub and tree species. However, it should be noted that several bees and butterflies can feed on shrubs and trees or on non-native plants (Kleijn and Raemakers 2008), but those insect species were still included. To ensure that taxonomic changes and collection tools/skills would not affect estimates of richness change, species that could not be easily distinguished in the analyzed time periods were lumped into species aggregates, based on information provided by specialists on each of the taxa (Supplementary material Appendix 1 Table A2, A3). Plant species with multiple subspecies or varieties were always aggregated under a unique name at species level. Richness change calculations We compare four 20-yr time periods (P0: 1930–1949, P1: 1950–1969, P2: 1970–1989, P3: 1990–2009), for which we have extensive geographically and temporally explicit databases of species records, and across which substantial changes in diversity have occurred (Carvalheiroetal. 2013). We characterize these periods in terms of nitrogen deposition, using information from the Netherlands Environmental Assessment Agency (<www.rivm.nl>; Veldersetal. 2010; see also Monteny and Hartung 2007). Using an approach developed and tested in a previous work (Carvalheiroetal. 2013), and applied for multiple taxa in different regions of the World (Hendriksetal. 2013, Eskildsenetal. 2015), we estimated the mean change in species richness for each group between P0 and P1, P1 and P2, and P2 and P3 at three spatial scales: local (10 × 10 km grid cells); regional (40 × 40 km); and for the country as a whole. We selected these three scales because of data availability and their relevance for spatial dynamics of species populations (the smallest scale covers the maximum foraging of the studied bees and butterflies, e.g. Hagenetal. 2011, but see Rao and Strange 2012), and for conservation management and policy making, which typically focus at local (e.g. when managing ecosystem services), regional/county and country level (e.g. re-introduction programs). In addition, comparisons of results at these different spatial scales give us a more complete picture of the overall patterns of change. For example, if large changes are detected for the whole country but not at regional and local scales, this suggests that such whole country changes are driven by species that have very limited spatial distributions, and hence only affect a small number of regional and local scale cells. Following Carvalheiro et al. (2013), for each cell that matched the selection criteria, richness change was estimated using a combination of both interpolation and extrapolation (for cell selection criteria and a detailed explanation of the method see Supplementary material Appendix 1). This approach allowed us to deal with unequal sampling intensity between grid cells or time periods, and with oversampling of rare species which may bias richness estimates. The approach leads to richness change estimates that are not correlated with sampling effort change (Carvalheiroet al. 2013), and has shown to be robust to extrapolations of sampling effort up to three-fold in range (Colwelletal. 2012). Subsequent analyses used a meta-analytical approach where the contribution of each datapoint is weighted based on the error associated with the richness change estimate. More specifically, the richness change estimates were log-transformed and analyzed using weighted general linear models (GLMw), with the inverse of variance (bootstrapped to correct for under/over-representation of singletons) applied as weight. Scripts were developed in R (R Core team) are available at <https://github.com/ lgcarvalheiro/richness.change>). Effect of nitrogen deposition on recent richness change patterns Geographically explicit information on nitrogen deposition (sum of reduced and oxidized nitrogen) for each 10 × 10 km grid cell in the Netherlands for the year 2009 were obtained based on a model developed by PBL (Planbureau voor de Leefomgeving) and the Netherlands Environmental Assessment Agency (< www.rivm.nl >; Veldersetal. 2010). While the intensity of nitrogen deposition has changed over time (Fig. 1), the spatial distribution of its sources within the Netherlands over recent decades has not changed greatly (Monteny and Hartung 2007, Veldersetal. 2010). Therefore, the values of nitrogen deposition obtained in 2009 can be considered as a proxy of the increase in nitrogen availability during the past several decades (i.e. regions with greatest nitrogen deposition values in 2009 are assumed to be the ones where nitrogen levels most changed from 1970 to 2009). While the focus of this study is on the effect of nitrogen availability, we also checked for possible confounding effects of landuse and climate. For estimating land-use change, we considered percentage of habitat considered suitable for a large number of pollinator species (as defined by Aguirre-Gutierrezetal. 2015 following Vogiatzakis et al.’s 2015 classification). Grassland, 213 moors/peat, forest and sandy soils were classified as ‘suitable for most species’, and agriculture (in Netherlands, it is mostly highly intensified monocultures), urban, water and swamps as ‘non-suitable for most species’. Land-use data were obtained from the geo-information department of Wageningen Univ. with an original resolution of 25 by 25 m pixels for 1980 (P2) and 2008 (P3). We then calculated the difference in percentage of habitat considered suitable to pollinators for each cell between P3 (1990–2009) and P2 (1970–1989). For climate we considered nine ‘bioclim’ variables (Ficketal. 2017) known to have an important effect on Dutch pollinators (Aguirre-Gutiérrezetal. 2016 and Supplementary material Appendix 1 Table A5). For each climatic variable we estimated change as the difference between the average of the annual values within P3 and P2. We then selected variables that had a Pearson correlation value lower than 0.5 between them, so that the number of excluded variables was minimized, and maintaining a similar number of temperature and humidity related variables. To test for interactive effects between nitrogen deposition and resource preferences on the patterns of richness change of plants, bees and butterflies we used general linear mixedeffects models (GLMM), with resource preference classes nested within cell identity as random structure. We tested if richness change values were spatially autocorrelated by comparing a model with linear and exponential correlation structures with a null model (R package ‘nlme’, Pinheiroetal. 2019). As no spatial effect was detected, we proceeded with the analyses with the R package ‘metafor’ weighing each data point based on the inverse of its variance. Nitrogen deposition (log-transformed to linearize relationship) and resource preference group were included as explanatory variables, also considering their interaction. To inspect if effects of climate or land use change were confounding the effects of nitrogen deposition, we included the two-way interactions between nitrogen deposition and several climatic and land use change variables in our models. Interactions between trait group and all environmental variables were also considered. We then selected the most parsimonious model for each group, based on the lowest Bayesian information criterion (BIC), using a minimum model that included nitrogen deposition, trait group and the interaction between both, and allowing a maximum of two climatic variables in each model. The R package ‘multcomp’ (Hothorn et al. 2008) was used to run post-hoc contrast test on the most parsimonious model detected for each group, and extract an estimate (and 95% confidence intervals) of the effect of nitrogen deposition for each group, after accounting for the effect of climate and land use change). Results Do resource preferences modulate historical patterns of biodiversity richness change? Our results show that, in addition to the expected impacts of other drivers (climate, land use), species’ nitrogen-related resource preferences do influence the patterns of richness changes of plants and pollinators. Local and regional (10 × 10 km and 40 × 40 km grid cells, respectively) plant richness declines between P0 (1930 and 1949) and P1 (1950 and 1969), a period during which nitrogen deposition levels increased (Fig. 1), were most accentuated for nitrophobous plants, particularly if they depended on insects for pollination (Fig. 3a, d). At country level no significant declines were detected, with even slight but significant increases in richness being detected for plants dependent on pollinators. This is an indication that although many species reduced their spatial range during this period, they did not go nationally extinct, and concurrently some species which were not detected in P0 (possibly due to very restricted spatial ranges) were detected or had colonized new locations in P1 (Fig. 3a). Between P1 and P2 (when nitrogen deposition had its most accentuated increase) there was, on average, an increase in local (10 × 10 km grid cells) plant richness. For pollinatordependent plants, these increases seem to have been independent of the plants’ nitrophily classes (Fig. 3b, e). However, for plants used by bees as food sources (including those not fully dependent on insects for pollination, Fig. 3g–i), increases were only detected for nitrophilous plants, while large-scale declines were detected for nitrophobous plants (possibly as a consequence of local scale declines detected during the previous time-period for this group of plants). As described in Carvalheiroetal. (2013), in recent time periods (i.e. between P2 versus P3, periods between which nitrogen deposition levels slightly decreased), average changes in plant local richness were less accentuated. Here we show that increases at large spatial scales were more pronounced for nitrophobous plants (Fig. 3c, f), suggesting that some species from this trait group, whose distributions were previously so restricted that they were not registered in the surveys (of several 40 × 40 km cells or even at country level), expanded their ranges. For flower visitors, declines detected between P1 and P2 (Fig. 4b, e), periods between which nitrogen deposition greatly increased and nitrophobous bee–plants lost space for nitrophilous plants, were most accentuated for nitrogen specialists (i.e. for either species feeding solely on nitrophobous or solely on nitrophilous plants). Country-level declines were particularly strong for bee species whose major larval resources are nitrophilous plants (rather than non-nitrophilous, as predicted). During the same period, butterfly richness also declined substantially, but the most susceptible butterfly species were the non-polyphagous species which were also nitrogen generalists. Among butterfly nitrogen specialists (i.e. specialized on either nitrophilous or nitrophobous plants), substantial declines were also detected, these being most accentuated for nitrophobous species and at finer scales (10 × 10 km). This pattern suggests both that some previously widespread species contracted their ranges, and that some species became extinct (or were no longer detectable) throughout the country. For species specialized on nitrophilous plants, declines were only detected at larger scales, suggesting that declines occurred mostly for species that were already spatially restricted. 214 Between P2 and P3, significant recoveries were detected in most bee groups at different spatial scales (e.g. at local level for polyphagous species, and at country level for nitrophilous bees), with no clear relation to diet choice detected. For butterflies, increases were restricted to species specialized on nitrophilous species, and declines were still detected for other groups between P2 and P3, these being significant at all scales for species specialized on nitrophobous plants. When repeating the analyses excluding polyphagous species from the nitrogen specialist groups, results were maintained (Supplementary material Appendix 1 Fig. A1). For polyphagous nitrogen generalists and for species whose major larval resources are nitrophilous plants, mild recovery signs were found at fine spatial scales. Are recent biodiversity change patterns still modulated by nitrogen availability? Although average changes in richness became less accentuated between P2 and P3, spatial distribution of nitrogen deposition played an important role in explaining the variability of recent patterns of change for both plants and pollinators (Fig. 5, see also Supplementary material Appendix 1 Fig. A2, A3). For all groups analyzed here, current levels of nitrogen deposition have negative effects on richness change patterns (i.e. either constraining recoveries or accentuating declines). For plants, such effects were significant only for those that are non-pollinator-dependent. For bees, the negative trends (associated with the very low number of 10 × 10 km cells with sufficient data for analysis for each group) were non-significant. For butterflies, nitrogen deposition effects depended on species resource preference. Only species whose major larval resources were nitrophobous plants did not show significant negative effects of N deposition, i.e. the negative effects detected in Fig. 4f were similar in areas with low and high current N deposition. While climate and land use changes do have a significant effect on many of the groups analyzed here (Supplementary material Appendix 1 Table A5 and Aguirre-Gutiérrezetal. 2016, 2017), the effects of nitrogen deposition were not Figure3. Nitrogen preferences and pollination mode affected plant richness change patterns. Change of species richness (estimated weighted mean ± 95% confidence intervals) of plants at different spatial scales (local scale: 10 × 10 km; regional scale 40 × 40 km; whole country) are presented for plants dependent on pollinators (a–c), plants not dependent on pollinators (d–f) and plants that are used by bees as food sources (g–i). Grey symbols represent nitrophobous species (i.e. plants that do not thrive in N rich environments); black symbols represent nitrophilous plants. Red dashed line represents no change (0%). Filled symbols indicate that change was significantly different from zero, otherwise symbols are open. All statistical details are presented in Supplementary material Appendix 1 Table A4. 215 significantly influenced by other drivers of change (i.e. no interactive effect between N deposition and other drivers were detected, Supplementary material Appendix 1 Table A5). However, effects of climate warming on pollinators depended on their nitrogen related diet preferences, with positive effects only being detected for bees that prefer nitrophobous plants (i.e. increases were more accentuated in regions subjected to greater climate warming) and negative effects being most accentuated for butterflies preferring nitrophilous plants (i.e. declines were most accentuated in regions where warming was greater) (Supplementary material Appendix 1 Fig. A4). Discussion Increases in nitrogen deposition have important impacts on plant communities and can propagate to higher trophic levels (Mansonetal. 2013, Pöyryetal. 2017). However, until now, little has been known about how such changes affect the spatial and temporal dynamics of diversity patterns of pollinators, and what role they played in past declines and more recent partial recovery of plant and pollinator local richness in NW-Europe (Carvalheiroetal. 2013). Here we show that historical and ongoing declines in pollinator richness were, at least partially, associated with nitrogen levels and larval diet preferences. Overall, our findings suggest that while other drivers (e.g. climate or land use changes) did play a strong role, soil eutrophication (i.e. increases in nitrogen availability) mediated how biodiversity has changed. Below we discuss the possible causes and implications of these findings. Limitations associated to richness change analyses In the absence of formal monitoring programmes, long-term databases containing validated species presence records collected at different times and by many different recorders, provide a valuable source (frequently the only source) of information on past and present species occurrences. However, such databases can have considerable bias, e.g. unstandardized sampling effort and overrepresentation of rare species. While the analytical approach here applied allows to deal with such exceptional databases (Carvalheiroetal. 2013), it is Figure4. Nitrogen related larval diet preferences affected pollinator richness change patterns. Change of species richness (estimated weighted mean ± 95% confidence intervals through time at different spatial scales (local scale: 10 × 10 km; regional scale 40 × 40 km; whole country) are presented. Red dashed line represents no change (0%). Triangles represent species with a more generalized diet in terms of plant nitrophily (i.e. major larval resources include both nitrophilous and nitrophobous plants), these being either polyphagous (black triangles) or nonpolyphagous (grey triangles). Circles represent nitrogen specialist species, i.e. whose major larval resources are either nitrophilous plants (black circles) or nitrophobous plants (grey circles). Filled symbols indicate that change was significantly different from zero, otherwise symbols are open. For bees whose major larval resources are nitrophilous plants, we had insufficient data (less than three cells) to run the analyses at 10 × 10 km and 40 × 40 km. All statistical details are presented in Supplementary material Appendix 1 Table A4. 216 important to consider the sensitivity of the estimates obtained to small changes in the analytical method. While the overall comparative conclusions (i.e. between groups or time periods) of this manuscript were not affected by such sensitivity tests (Supplementary material Appendix 1 Table A7, e.g. nitrophilous butterflies outperform nitrophobous butterflies independently of the type of method and weight), and we do have a recommended method and weight (applying a combination of interpolation and extrapolation and use of bootstrapped inverse of variance as weight, Carvalheiroetal. 2013), estimates of mean values of change in isolation need to be interpreted with caution. We also emphasize that changes in species abundance are not detectable with our approach (changes are only detected in case of local extinction or colonization) and information on species abundance is essential to better understand community dynamics over time. Effect of resource preferences on historical patterns of biodiversity richness change As discussed in a previous study (Carvalheiroetal. 2013) the accentuated variation in richness change patterns across time periods, detected for most groups, is likely related to accentuated changes in land-use patterns during the studied period (EEA 2010, FAO 2012), as well as recent increased investment in environmentally–friendly practices (Kleijn and Sutherland 2003). Here we go a step further and show that, as expected, nitrogen availability did play an important role in defining the overall patterns of plant richness, particularly between 1930 and 1970. Vascular plants are susceptible to increases in nitrogen levels (and associated acidification, Goulding 2016), with known impacts on plant diversity (Bobbinketal. 2010, Isbelletal. 2013), and physiology (Santiagoetal. 2011, Barbosaetal. 2014). Species adapted to low nutrient levels (nitrophobous) can be very efficient in extracting nutrients, and an increase in nutrient levels does not always lead to a higher uptake (Fichtner and Schulze 1992, Barbosaetal. 2014). On the other hand, species adapted to nutrient-rich conditions (nitrophilous plants) tend to react much more strongly to increases in the soil N availability (Fichtner and Schulze 1992). Consequently, the initial response to such increased N availability can lead to a greater biomass of nitrophilous plants (Aberetal. 1989, Bobbink 1991), outcompeting co-occurring nitrophobous plants (Tilman 1993, Dirnböcketal. 2014) and contributing to their local extinction (Fig. 3a, d). The accentuated parallel changes that were detected for most plant groups (e.g. increases at country level before the 1970s, Fig. 3a, followed by increases of local richness, Fig. 3b) during the period of increased nitrogen deposition, are likely due to other environmental factors that are known to affect plants independently of their dependence on pollinators (e.g. land use, climatic changes). The fact that the advantage of nitrophilous over nitrophobous plants became less accentuated over time for pollinator-dependent plants suggests that other environmental changes played a more important role in this particular group. For example, nitrophobous plants might be favoured in restoration programs (e.g. plantation of native flower margin plantations in rural and suburban areas has been done throughout the whole country since early Figure5. Recent recoveries in local (10 × 10 km) richness change rates were limited by current N deposition levels. Effect of nitrogen deposition on different functional groups of plant, bee and butterfly richness change between P2 (1970–1989) and P3 (1990–2009) is presented (see statistical details in Supplementary material Appendix 1 Table A5, A6). All explanatory variables had a Pearson correlation value lower than 0.5. Red dashed line represents no effect of nitrogen deposition. For plants, black and grey circles represent nitrophilous and nitrophobous plants, respectively. For pollinators, circles represent nitrogen specialist species, i.e. whose major larval resources are either nitrophilous plants (black circles) or nitrophobous plants (grey circles). Triangles represent species with a more generalized diet in terms of plant nitrophily (i.e. major larval resources include both nitrophilous and nitrophobous plants), these being either polyphagous (black triangles) or nonpolyphagous (grey triangles). Filled symbols indicate that change was significantly different from zero, otherwise symbols are open. For bees, whose major larval resources are nitrophilous plants, we had insufficient data to run the analyses. 217 1990s, Kleijn and Sutherland 2003), whereas nitrophilous plants that tend to grow in agricultural habitats might have suffered more from increasing intensity of farm management. Before 1950, differences in local richness declines between nitrophilous and nitrophobous plants were more accentuated for plants that depend on pollinators (Fig. 3a), suggesting that pollinators that feed on nitrophobous plants were more affected than those with different feeding preferences. Although no substantial pollinator richness declines were detected during these early time periods, local density declines or changes in their foraging behavior might have already occurred, preceding the richness declines detected later on (Fig. 4b, e). This suggests a time-lag for pollinator richness changes in relation to plants. Insects have faster life cycles than most plants, and changes in density (WallisDeVries and van Swaay 2017) or visitation patterns (Ramosetal. 2018) in response to changes in nitrogen can be quickly detected. However, that does not necessarily lead to changes in richness at landscape level. Given their mobility, insects can rapidly change their local distribution and foraging patterns to avoid the negative effects of excess nitrogen in localized areas. Consequently, richness changes at landscape or regional scales take longer to occur. On the other hand, plants can only change their distribution over generations (which are typically longer than those of insects), so they are less able to counteract local and regional extinctions, i.e. impacts can be detected earlier. In addition, the time lag detected for insects could be a response to changes in resource quality (Adler and Irwin 2011, Hooveretal. 2012, Nijssenetal. 2017), as discussed below. Increased availability of nitrogen can change the quality of floral resources (e.g. changes in sugar, amino acid or secondary compound content, Petanidouetal. 1999, Gardener and Gillman 2001, Gonthieretal. 2011), which can affect visitation or oviposition patterns (Chenetal. 2008). Further studies testing if such changes in resource quality are more accentuated for plant species capable of using the excess nitrogen (i.e. nitrophilous plants, Chapin 1980), or if some parts of the plants (e.g. leaves versus pollen) are more affected than others, would be needed to fully understand the patterns of change here detected. For instance, why are declines at coarse scales more accentuated for bees that were specialized on nitrophilous plants (Fig. 4b)? Or why were effects more accentuated for butterflies than for bees? Also, further attention should be paid to potential interactive effects between nitrogen availability and changes in other biogeochemical cycles (e.g. water, carbon, phosphorous, Sardansetal. 2011, Myersetal. 2014). Such interactive effects could lead to further changes in plant stoichiometry and nutrient content. While many mechanisms request further investigation, overall, our results suggest that eutrophication of the environment, and consequent changes of the flora, has contributed substantially to the decline of pollinators (i.e. global domino effect; Schleuningetal. 2016, Pöyryetal. 2017). Among the three species groups that include polyphagous species, we found consistently for bees and butterflies that nitrogen generalists were much more resilient to environmental change through time than either group of nitrogen specialists. The fact that pollinators with more diversified diets were less affected could be explained by their flexibility and ability to forage in a diverse set of habitats (Rogeretal. 2017), thus reducing or even excluding the contribution of plants with toxic compounds to their diets. Alternatively, these species might be better able to detoxify these compounds, a process that requires the production of specific enzymes from dietary proteins (Manson and Thomson 2009, Slansky 1992). However, larval habits are much more specialized in butterflies than bees, with non-polyphagous species frequently having strong preferences for a single plant genus; whereas a large number of non-polyphagous bees use a diverse set of plants from the same family as preferred resources (Vanderplancketal. 2017). While the high mobility of both insect groups may have delayed the detection of community level impacts, this may explain why among species with a more specialized diet (a trait known to be linked with high susceptibility to environmental changes, Carvalheiro et al. 2010), butterflies were more susceptible to declines than were bees (Fig. 4b, e), even if they were nitrogen generalists. Are recent biodiversity change patterns still modulated by nitrogen availability? During recent decades, when nitrogen deposition has started to decrease, nitrophobous plants outperformed nitrophilous plants (Fig. 3c, f), suggesting a change in vegetation composition. While previous studies suggest that impacts of nitrogen on plants can last for decades (Isbelletal. 2013), the slightly positive changes detected here in recent decades could partly be a consequence of the reduction in nitrogen deposition, despite levels being still far above natural levels (Fig. 1). However, the fact that, despite increases in richness, plants were affected negatively by N deposition (Fig. 5, i.e. recoveries were less likely to occur or less accentuated in areas with high N deposition), suggests that other drivers (e.g. changes in land use, habitat restoration, including practices that aim to reduce N soil levels, WalliesDeVries and Bobbink 2017) are responsible for recoveries detected in plant richness, and that N deposition does constrain such positive effects. As for pollinators, even though nitrophobous plants are recovering, pollinators that prefer nitrophilous plants continue to outperform those that prefer nitrophobous ones (Fig. 4c, f). For bees, increases of nitrophilous species were even detected at country level. Such a pattern is partially driven by species that are expanding their range. For example, the wide distribution of the nitrophilous host plant of the wild bee Colletes hederae (Hedera helix) allowed the quick expansion of C. hederae throughout Europe (Dellicouretal. 2014), being first detected in the Netherlands in 1997. Several butterfly species with a preference for nitrophilous diet are also known to be expanding in NW Europe (Betzholtzetal. 2013). However, even for nitrophilous groups that show signs of recovery, and for species without nitrogen-related