Focusing only on native species assemblages, here we apply robust richness estimation methods to evaluate richness changes during the past 80 years in three European countries (Great Britain, the Netherlands and Belgium; ca . 29.2 million records in total) at multiple spatial scales (from 10 × 10 km cells to the whole of each country), and also evaluate changes in patterns of biotic homogenisation (due to the spread of native species or loss of spatially restricted species). We find that extensive species richness loss and biotic homogenisation occurred before 1990. However, such negative trends have slowed down (or recovery has set in) for several taxa during the most recent decades.

Biesmeijer et al . ( 2006 ) use rarefaction methods to estimate changes in diversity in plant and pollinator taxa based on historical collections before and after 1980. They found evidence of declines in bees and pollinator‐dependent plants in Great Britain and the Netherlands. However, whereas it is well known that during the past decades the studied countries underwent substantial climatic (EEA, European Environment Agency 2004 ), land‐use (Haines‐Young et al . 2003 ; EEA, European Environment Agency 2010 ; FAO, Food & Agriculture Organization of the United Nations 2012 ) and environmental policy changes (Kleijn & Sutherland 2003 ), it remains unclear how the rate of biodiversity loss has changed over time.

Ongoing species richness loss and biotic homogenisation are a growing concern for society (Pereira et al . 2010 ; Cardinale et al . 2012 ). The increasing awareness of the role of human activities in these declines and their potential impacts on our health, food supply, ecosystem services and well being (e.g. Isbell et al . 2011 ) have led to policy changes towards the protection of biodiversity, many implemented after 1990 (Kleijn & Sutherland 2003 ). The Convention on Biological Diversity's main target is to substantially reduce the rate of biodiversity loss at global, regional and national levels (CBD, Secretariat of the Convention on Biological Diversity 2010 ), and the EU committed itself to an even more ambitious target of halting biodiversity loss. Most agree that these targets have not been met (e.g. Butchart et al . 2010 ); but in the absence of standardised long‐term monitoring programmes (which exist for only a few taxa and countries) it is difficult to assess if biodiversity loss trends are changing.

Materials and methods

For each of the studied countries (the United Kingdom, the Netherlands and Belgium), four 20‐year time periods were compared: 1950–1969 (hereafter P1), 1970–1989 (P2), 1990–2009 (P3) and whenever possible 1930–1949 (P0), covering a period during which the studied countries were subjected to substantial changes in land‐use, climate and environmental policies. The analysis of data collected prior to 1950 (1930–1949) was limited to taxonomic groups and countries for which data quality allowed this.

Three groups of flower visiting insects which contribute to plant pollination were studied here: bees (Hymenoptera: Apoidea), hoverflies (Diptera: Syrphidae) and butterflies (Lepidoptera: Papilionoidea and Hesperioidea). Given the importance of bees for pollination, and the recognised high vulnerability of bumble bees (Apidae, Bombini) (Carvell et al. 2007; Williams & Osborne 2009), we analyse bumble bees separately from other bees (the honeybee, Apis mellifera L. was not considered, as its presence in the focal region is largely dependent on apiculture practices). Plants were divided into three groups according to their dependence on flower visitors for pollination: fully dependent, intermediate dependency and independent, as described below. Exotic species (archaeophytes and neophytes) represent a relatively small minority of species, and analyses including exotics generated similar result patterns for most spatial scales; therefore, we focus only on the results obtained using data on native plant species. For details on the databases sources see Table S1.

Information on the dependency of plants on insects for pollination was obtained for ca. 50% of the plant species (2481 of 4943 species) from the Ecological Flora database (see Table S1). When information was not consistent among the three databases a given species was classified as intermediately dependent. On the basis of this information we detected that only 6% of the genera (23 of 377 genera, single species genera excluded) and only 35% of the families (36 of 103 multi‐species families) included both insect‐dependent and non–insect‐dependent species. Therefore, for the ca. 50% of species for which information was not available in any of the databases, we used information from genus and family level to classify the species (e.g. all Poaceae were classified as pollinator independent). If species within a genus had different dependencies, species would be classified as intermediately dependent. As many tree species are planted for ornamental or reforestation purposes, we repeated the analyses with and without all tree species to distinguish between effects of natural vs. man‐induced range changes. As no major differences were found between analyses with and without trees, we presented the results including all native plant species.

Analyses of species richness change In the absence of formal monitoring programmes, long‐term databases containing validated species records collected at different times and by many different recorders provide a valuable source of information on past and present species occurrences. As these records were usually not collected following standardised sampling schemes, their analysis faces several methodological challenges. First, improvement of taxonomic and collection tools/skills through time can lead to more species being registered in more recent time periods simply due to detectability changes. 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 past or present time periods were lumped into aggregate species, based on information provided by the specialists on each of the taxa (Bees: DMichez, PR, SR; Hoverflies: MR, FVM, Butterflies: RF, DMaes, MW, Plants: QG, SH, WVL, BO, JS). Plant species with multiple subspecies or varieties were always aggregated under a unique name. Secondly, estimates of richness change are generally dependent on the spatial scale at which they are evaluated (Keil et al. 2010). However, such scale‐related differences in estimates provide valuable information to help understand the patterns of change, and for conservation management (Whittaker et al. 2001). For example, range expansions affect richness values of multiple fine‐scale cells, and hence have a substantial effect on the mean change value at finer spatial scales (Gotelli & Colwell 2001), but no effect at country level (if the species was present in both periods, and simply changed its distribution pattern). Conversely, country‐level extinctions of spatially restricted species and species introductions will affect coarse (e.g. national)‐scale richness, while only influencing a few fine‐scale cells (see Cassey et al. 2006). Therefore, whenever data quality allowed, we repeated richness change analyses at multiple spatial scales (e.g. 10 × 10, 20 × 20, 40 × 40, 80 × 80 or 160 × 160 km grid cells as well as for the whole country). As only one 160 × 160 km grid was available for the Netherlands and Belgium, this scale was not considered in these countries. Finally, unequal sampling intensity between grid cells or time periods and oversampling of rare species can bias richness estimates (Hellmann & Fowler 1999; Garcillan & Ezcurra 2011; ter Steege et al. 2011). Here, we combine robust richness estimation methods with a meta‐analysis approach to deal with the unstandardised nature of historical collections. When comparing two periods with unequal sampling, previous studies have used rarefaction to estimate (or rather, interpolate) the value of richness in the more sampled period that would be expected if sampling effort had been equal to that of the less sampled period (e.g. Biesmeijer et al. 2006; Keil et al. 2010). As an alternative, Colwell et al. (2012) have shown that species accumulation curves can be extrapolated (up to threefold of the real sampling effort), producing reasonable estimates of richness further in the accumulation curve. This approach has been used here to estimate richness of the least‐sampled period, extrapolating to a point where sampling effort is equal between the two periods (pre‐period X 1 [n] and post‐period X 2 [n]). However, where the more sampled period had more than threefold the number of records of the less sampled period, we combined extrapolation and interpolation in our comparisons: extrapolating the smaller sample's accumulation curve up to three times its sample size, and rarifying the larger sample down to this same ‘comparison point’. We repeat this process for every cell, calculating relative richness change between the pre‐ and the post‐period as . For analyses of richness change, we applied a log transformation ( , hereafter termed ‘logratio’) to normalise residuals. Its sampling variance (VAR QX ) is approximated as As in the absence of singletons, the assemblage is considered to be well sampled, and SD of X will be zero (Colwell et al. 2012), to account for any under‐ or overestimation of singletons and doubletons (collectors may put more effort on registering the rarest species, Garcillan & Ezcurra 2011, or try to obtain two specimens to capture morphological diversity of the species, e.g. males and females; flowering and fruiting plant specimens), we a priori excluded cells with very poor quality of sampling and we used bootstrapping to calculate SD of X (by re‐sampling the data we create a corrected estimator of the number of unseen species, where the effect of under/oversampling of singletons is reduced). For details on selection criteria applied in richness change analyses, see legend of Table S2. Meta‐analysis techniques were then applied to obtain an overall weighted value of richness change (Qw) and assess if the mean value of change at a given scale was significantly different from zero. Using the rma.uni function of the R package metaphor (Viechtbauer 2010), each grid cell was weighted based on the inverse of the variance, so that cells with more reliable estimations have a higher weight in the analyses (Hartung et al. 2008). To check if the method completely corrected for bias due to differences in sampling efforts, we included the log of the relative difference in the number of records between the two time periods as a covariate. If accumulation curve estimations did not completely remove the bias due to sampling effort (i.e. whenever ΔR had a significant effect on estimated richness change across sites), we calculated the partial residuals after removing the effect of sampling effort for each cell to obtain unbiased estimates of richness change for each grid cell. We then assessed in which cells richness had significantly increased or decreased using a z‐test, for each species group, country and time period. Finally, we evaluated if values of richness change were spatially autocorrelated by comparing a model with spatial autocorrelation structure (linear or exponential) with a null model using a log‐likelihood ratio test and compared the Akaike Information Criterion values of each model. For verification of the sensitivity of our results, we repeated all richness analyses using only interpolation and only extrapolation and we compare the results using re‐sampled and non–re‐sampled SD and variance when analysing data.