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

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.


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 andGarcia 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, 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 andBobbink 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 andSteiner 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 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 (abiotic-or 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). 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, 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.
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.
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 andStrange 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 andHartung 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, 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) 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.
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, the effects of nitrogen deposition were not  Table A4. 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. 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 andSchulze 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. 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 andThomson 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 preferences, recent nitrogen deposition continues to have a negative effect (Fig. 5), constraining pollinator richness recovery. This suggests that any range expansion that may be driving the increases in local richness is likely occurring towards areas with lower N deposition. Furthermore, the recent declines in nitrophobous butterflies (Fig. 4f ) are not spatially related to N deposition patterns (Fig. 5), which suggests that this specific group of butterflies is not as able to take advantage of the reductions in N deposition as other groups. A possible explanation for this unexpected result is that increases in host-plant distribution of this butterfly group (reflected by the plant richness increases shown in Fig. 3) did not translate into increases in local abundance to the degree necessary to support consumer populations. Unfortunately, the sort of standardized monitoring data required to assess abundance changes is not available for most insect groups, particularly before 1990. Also, if N deposition levels still exceed the critical levels for nitrophobous butterflies at landscape scale (WallisDeVries and van Swaay 2017), the recovery detected for nitrophobous plants at local scales might not be sufficient to support the recovery of its consumers.
As the effect of temperature on plant richness was similar across plant trait groups (i.e. no interaction detected, Supplementary material Appendix 1 Table A5), the fact that for pollinators climate warming effects depended on nitrogen related diet preferences (Supplementary material Appendix 1 Fig. A4) is unlikely explained by the effects of climate change on the distribution of their host plants. Further studies exploring interactive effects between climate and nitrogen availability in more detail are essential to understand if the success of conservation practices aiming to restore soil nutrient conditions may be affected by global warming.

Concluding remarks
Our results suggest that anthropogenic nitrogen deposition has long-term consequences for the composition (functional group representativeness) of plant and flower-visitor assemblages. Species tolerance for raised nitrogen levels and dietary preferences should be considered in future policy and management plans that aim to halt biodiversity loss and increase the provision of ecosystems services (e.g. guaranteeing that levels of fertilizers applied to crop fields do not exceed recommended dosages for a specific crop in a specific region). While insects often respond faster than plants to other environmental drivers (Krauss et al. 2010), our results suggest that community-level effects mediated by plants will take much longer to be detected (lagging behind changes in plant assemblages, see also Sang et al. 2010). This lends support to the idea that long-term effects of eutrophication (and the spatial scale at which changes occur) should be considered when evaluating and predicting impacts of global environmental changes.
The effects mentioned above may also depend on plant (e.g. symbiotic relationships with nitrogen fixing bacteria) and pollinator (e.g. sociality, nesting or larval microhabitat) species traits may also modulate effects of nutrient enrichment on biodiversity (WallisDeVries and van Swaay 2006, Elliott et al. 2008, Manson andThomson 2009). Future studies looking in more detail into species functional traits and into changes in pollinator foraging habits will allow to better disentangle the mechanisms regulating the propagation of nutrient enrichment from plants to pollinators.
Finally, interactive effects between eutrophication and other environmental changes (e.g. climate or changes in habitat quality due to fragmentation or pollution, see Gonzalez-Varo et al. 2013 for a review) need to be better explored to fully understand the dynamics of plant-pollinator relationships under global change.
Funding -This project was partly funded by The European Commission Framework Program (FP) 7 via the Integrated Project STEP (grant 244090). LGC was supported by Fundação para Ciência e Tecnologia (FCT) and European Union via the programa operacional regional de Lisboa 2014/2020 (project EUCLIPO-028360) and by the Brazilian National Council for Scientific and Technological Development (CNPq) via the project BJT 300005/2015-6. JB was supported by COST action FA-1307 Super-B (Sustainable Pollination in Europe). AH was supported by the Estonian Ministry of Education and Research (IUT20-29), and by the European Union through the European Regional Development Fund (Centre of Excellence EcolChange). JAG was supported by the Netherlands Organisation for Scientific Research (NWO) project number 019.162LW.010.