Meta-analysis of multidecadal biodiversity trends in Europe

Local biodiversity trends over time are likely to be decoupled from global trends, as local processes may compensate or counteract global change. We analyze 161 long-term biological time series (15–91 years) collected across Europe, using a comprehensive dataset comprising ~6,200 marine, freshwater and terrestrial taxa. We test whether (i) local long-term biodiversity trends are consistent among biogeoregions, realms and taxonomic groups, and (ii) changes in biodiversity correlate with regional climate and local conditions. Our results reveal that local trends of abundance, richness and diversity differ among biogeoregions, realms and taxonomic groups, demonstrating that biodiversity changes at local scale are often complex and cannot be easily generalized. However, we find increases in richness and abundance with increasing temperature and naturalness as well as a clear spatial pattern in changes in community composition (i.e. temporal taxonomic turnover) in most biogeoregions of Northern and Eastern Europe.

T he current biodiversity crisis, manifested in a global decline of species, affects many taxonomic groups and biotic realms [1][2][3][4] . These changes may be less evident at specific locations, since local factors, such as small-scale colonization and species turnover may compensate or even counteract trends occurring at larger spatial scales. Heterogeneity in patterns of change in biodiversity observed at a local scale 5,6 has been described in several studies. For example, strong declines in local biomass and distribution have been reported for terrestrial insects 7-10 and birds [11][12][13] , but reports of local increases in biodiversity are also widespread and span multiple taxonomic groups 5 including freshwater invertebrates 14,15 , fishes 16 , birds 12,13 and plants [17][18][19] .
Ecosystem functions and their benefits to people at local to global scales ultimately depend on the taxonomic and functional diversity of local communities 5,20 . The same relationship applies to conservation measures: while they need to be harmonized at larger scales, most measures need to be tailored to local conditions. Therefore, it is crucial to understand the variation in biodiversity trends between localities and how and why these may vary across biota, regions and local conditions 21 .
Analyses of the trends in local biodiversity over large spatial scales and multiple taxonomic groups are needed to fully understand the patterns of local biodiversity change and the discrepancies between local and global biodiversity trends. Unfortunately such syntheses are rare. Current evidence from multi-taxa biodiversity trend comparisons is limited and equivocal, and direct comparison across studies is hampered by substantial differences in the temporal and spatial coverage and resolution of the data. Dornelas et al. 6 studied 100 time series of terrestrial, freshwater and marine taxonomic groups around the world and found no systematic temporal changes in α-diversity of local communities. However, the authors detected significant increases in β-diversity and increasing trends in species richness for terrestrial plants in the temperate region 6 . The most comprehensive study to date to our knowledge 22 showed that species turnover, i.e., a measure of temporal community variability, is stronger in marine than freshwater and terrestrial assemblages and that this is often decoupled from changes in species richness. However, these results were based on a relatively coarse spatial and temporal resolution, and mostly short time series. Other studies focused on a smaller number of sites or more restricted biotic scope. A study of 22 sites in different realms from Central Europe detected stronger effects of temperature on population trends in terrestrial than aquatic habitats (i.e. populations of terrestrial "species with warmer temperature preferences increased more than [terrestrial] species with colder temperature preferences" (page 3), but found no such relationship for aquatic taxa) 23 . Finally, Gibson-Reinemer et al. 24 found that the increasing species turnover in mountain communities was stronger in ectotherm communities and in tropical compared to temperate regions. These results point towards region, biotic, and realm specific patterns in local biodiversity trends, but a comprehensive overview directly comparing trends among complementary time series has been lacking so far.
The present study analyzes trends in various taxonomic groups measured at specific locations in nine different biogeoregions across the European continent. We ask: (1) are long-term trends in biodiversity detectable at individual localities across Europe? (2) If so, to what extent are such trends attributable to changes in climate at a regional scale and/or changes in local conditions? (3) Do observed trends in biodiversity vary predictably among biogeoregions, realms and taxonomic groups and, thus, inform our understanding and prediction of larger scale patterns? Due to the ongoing global change we expect to observe: (1) long-term reductions in biodiversity indicated by declining species richness and abundances (2) increasing variability in community composition, indicated by higher temporal species turnover and (3) differential responses across taxonomic groups and biogeoregions, reflecting differences in the extent of and vulnerability to climate change (i.e. Southern Europe should be more negatively affected than Northern Europe 25 ).
We compile 161 long-term (minimum 15 years) biomonitoring time series from 115 sites, mostly belonging to the International Long-Term Ecological Research network (ILTER 26 ), in 21 European countries ( Fig. 1), covering nine biogeoregions, three realms and eight taxonomic groups (Fig. 1 We then apply a meta-analytic approach to identify the patterns among biogeoregions, realms and taxonomic groups. Our time series have an unbalanced distribution across biogeoregions, realms and taxonomic groups, which is a common issue in macroecological studies 6,22,23 . We account for the unbalanced design by testing the robustness of our results using a sensitivity analysis (Supplementary notes) and when interpreting the results. Our results show that trends in common biodiversity metrics differ among biogeoregions, realms, and taxonomic groups. Particularly, we find stronger changes in community composition in Northern and Eastern Europe and increases in richness and abundance with increasing temperature and naturalness.

Results
Overall trends. We observed overall increasing trends in taxonomic richness, diversity and turnover, across all time series, while there was no significant trend in abundance (Table 1).
Biodiversity trends. The trends of all biodiversity metrics differed among biogeoregions (Table 1). Abundances declined in the Atlantic region and increased in the North Sea over time. Species richness and diversity increased in the Black Sea, Boreal region and North Sea. Species turnover increased over time in time series from the Alpine, Boreal and Continental regions, and decreased in time series from the Black Sea and North Sea ( Table 1, Fig. 2). We did not detect any clear trends in abundance in any realm, while trends in taxon richness, diversity and turnover varied across realms. We recorded increasing taxon richness in the freshwater realm, increasing taxon richness and diversity in the marine realm, and increasing turnover in the terrestrial realm (Table 1; Fig. 3).
We found a significant decline in the abundance of terrestrial invertebrates. Species richness, diversity and turnover trends differed among taxonomic groups (Table 1). Species richness and diversity increased in birds and aquatic invertebrates. This contrasted with the decreasing diversity in benthic algae (note that algae data were only available for sites within a single river catchment). Turnover trends significantly increased for plants. Other taxa did not show a trend in the respective biodiversity metrics ( Table 1, Fig. 4).
Influence of climatic trends and site characteristics. The most important predictors (i.e., highest absolute values of z-scores) of abundance and richness trends were site naturalness (i.e., a measure of local anthropogenic pressure), temperature trend (i.e., S-statistics of air or water temperature trend, according to the realm) and their interaction ( Table 2). Increases in temperature and site naturalness were correlated with increased abundance and richness. The negative interaction between site naturalness and temperature trend indicates stronger increases in abundance and richness with increasing temperature at sites with lower naturalness than at sites with high naturalness ( Supplementary  Fig. 1). The most important predictor of diversity trends was longitude (positive correlation). Species turnover trends were mostly affected by elevation (positive correlation) and by the positive interaction between site naturalness and temperature trend, indicating stronger increases in turnover with increasing temperature at sites with higher naturalness than at sites with low naturalness ( Table 2, Supplementary Fig. 1). Study length had much less influence on the trends, indicating that the differences in study intervals among time series did not bias our results ( Table 2).

Discussion
Our results demonstrate considerable heterogeneity in the extent and direction of change in biodiversity metrics in recent decades between biogeoregions, realms and taxonomic groups in Europe. We identified increases in taxon richness in Northern and Eastern Europe. Records in these regions are primarily represented by terrestrial and aquatic invertebrate datasets. These observed patterns are in line with recent modelled predictions of future responses to global changes 25 and are likely due to climatechange induced poleward range shifts across taxa [27][28][29][30] . We also detected declines in species abundances for terrestrial invertebrates and in the Atlantic biogeoregion (from which most of the Biodiversity trends for the whole dataset (overall trends) and within the different biogeoregions, realms and taxonomic groups, as resulting from meta-analysis mixed models. Note that only significant results (p ≤ 0.05) are reported for the biogeoregion, realm and taxonomic groupspecific analysis.
terrestrial invertebrate time series in our dataset are derived). Our results, based on data with larger spatial, temporal and taxonomic coverage, and finer temporal resolution, corroborate recent reports of worldwide declines of local terrestrial insect communities [31][32][33] . Several studied biogeoregions, realms and biotic groups showed no significant trends in biodiversity metrics. Other studies on local changes in biodiversity also detected no overall changes 5,22,34 in apparent contradiction to the documented global-scale biodiversity loss (e.g. IPBES 4 ). However, the extinction of rare species, by definition, is often restricted to the very few local places where these species occur and may thus go undetected in large extent quantitative studies 35 . Furthermore, the loss of specialist taxa could be compensated locally by the colonization of more tolerant and generalist species 36 and/or invasive taxa, resulting in biotic homogenization despite no change in overall species richness 36,37 . Such patterns were, for example, observed in vascular plants at coarse scales in Europe due to the extinction of rare native species and spread of alien species 38 . Thus, temporal changes in taxonomic composition, i.e., turnover, are likely to be more sensitive than simply taxonomic richness (i.e., alpha-diversity) and abundance in responses to global change 39,40 . Accordingly, we show that temporal changes in taxon turnover are more pertinent across biogeoregions than the other three studied biodiversity metrics. The observed increases in temporal taxon turnover are spatially structured across Europe, involving mostly biogeoregions of Northern and Eastern Europe. Such changes in community composition can reflect non-equilibrium dynamics (including time-lags and transient phenomena 41 ) with, for example, climate change 42 , pollution or the introduction and spread of alien species 43 . Moreover, over the past 30-40 years (i.e. the period covered by most time series in our dataset) positive impacts of environmental regulation, e.g., reductions in atmospheric emissions and hence acid rain, could also be important drivers of biotic change (see e.g., Monteith et al. 44 ). In a recent analysis of UK vegetation data, an overall increase in vegetation species richness was linked to recovery from acidification 45 . Accordingly, increases in species turnover could also reflect a process of biological recovery from past disturbances. Another possible interacting factor is the change in land use, which follows different temporal trajectories in different European regions, and thus could concur in explaining regional differences in biodiversity trends 46 . Future research should clarify to what extent the increasing taxon turnover is led by the spread of generalist and invasive species compared to declines in rare species, and whether the observed consistent with other studies, e.g. 14,25 . However, we could also identify a combined effect of naturalness and temperature on abundance, richness and turnover trends. More specifically, and counterintuitive to common expectation, we found that sites considered to be in a less natural state are those experiencing the  strongest changes in biodiversity metrics with increasing temperature (i.e. steeper increase in abundance and richness and steep decline in turnover). We speculate that degraded sites are more prone to invasion of generalist and invasive species, while natural sites may be more resilient 49 . A major obstacle in assessing biodiversity trends is that longterm data are unevenly distributed among taxonomic groups and are biased towards charismatic taxa, such as birds and mammals, and towards groups with long study traditions, such as vascular plants and marine fishes, while invertebrates are relatively neglected, except butterflies and bees in some countries. Our study encompasses an unusually large variety of taxonomic groups, including those largely overlooked in previous large-scale biodiversity studies 6,20 . Some of these overlooked groups that we included in our analysis, such as aquatic invertebrates, showed unexpected increases in richness and diversity, likely due to the aforementioned processes (e.g. recovery from stressors, spread of generalist or invasive taxa and taxa adapted to warmer temperatures). On the other hand, we recorded declines in species abundances for terrestrial invertebrates in the Atlantic biogeoregion, consistent with previous findings 31,32,50 . Our results therefore emphasize that patterns of change in the biodiversity of the most studied 'iconic' groups cannot be extrapolated across other taxa.
These findings reiterate the need to not only maintain but to increase the numbers of long-term monitoring schemes of local ecosystems. Such schemes can provide unique insights for ecologists and conservationists, yet they are often threatened by a lack of support as they do not fit into the temporal extent of most funding schemes. In contrast to space-for-time or aggregated snapshot (e.g. opportunistic data or atlases of species distributions) approaches, the long-term tracking of communities in specific locations minimizes the risk of biases related to shifts in sampling locations, sampling areas and protocols 51 .
Our dataset is, to the best of our knowledge, the most comprehensive in terms of the spatial and temporal extent and taxonomic and realm representation of high temporal resolution long-term biodiversity monitoring records in Europe. Most sites included in this study are part of the global ILTER network and the vast majority are characterized by low anthropogenic pressure, which may have led to an underestimation of the true scale of biodiversity changes at continental scale. Perhaps most importantly, most studied sites are shielded from direct effects of changes in land-use and loss of habitat, e.g., conversion to intensive agriculture and urbanization. To overcome such sources of bias, long-term monitoring programs should include a larger representation of more intensively used (e.g. agricultural) areas and incorporate sites vulnerable to significant anthropogenic perturbations 51 . On the other hand, our approach reduces the potential risk of tracking immediate biodiversity responses to localized disturbance or successional recovery processes, which can bias the estimates of biodiversity change 52 .
Although most of our time series do not predate the 1980s (only one study goes back to the 1920s), we were still able to detect evidence of substantial reorganization of communities within relatively short time frames. However, our data, as well as most of the data used in other studies 6,10,22,23 , is being considered in the absence of the longer historical perspective required to capture the overall changes in biodiversity in the Anthropocene, as we lack baseline data from times when human impacts were lower (e.g. pre-industrial era). This limitation is a common issue in studies of long-term biodiversity change (but see ref. 53 ) that can have a detrimental effect in restoration ecology 54,55 . The lack of baseline data could be overcome by the integration of ecological and paleobiological approaches (e.g. Battarbee et al. 56 ), which could extend the temporal dimension back to, e.g., the end of the last ice age. Such approaches could have important conservation implications 57 and allow for comparisons between impacted and unimpacted sites within the same geographical area, which could reveal differences between changes driven by climate and by direct anthropogenic pressures (e.g. changes in land use or pollution) 58 . Looking into the future, there is an urgent need to harmonize biodiversity monitoring schemes 59,60 that would also allow improved up-and downscaling of trends as well as integrating cross-domain feedback loops 61 .
The inherent complexity of ecological systems manifested by diverging long-term responses of local biodiversity still hamper any upscaling of these trends to a continental or even global scale. This might explain the partly contradictory results not only within but also among large-scale studies that are based on local biodiversity data 6,22 . For example, our study revealed higher taxon turnover in the terrestrial realm at European scale while, at a global scale, turnover seems higher in marine assemblages 22 . However, these studies do agree on the lack of an overall decline in species richness and on the increase in taxon turnover over time. We argue that these contradictory results can be driven by the insufficient number and quality of systemic and harmonized biodiversity monitoring activities at a local scale and by the insufficient length of the underlying time series. Regarding the latter, we should bear in mind that prior to the start of most of our biodiversity time series (mainly in the 1980s), many species had already declined in abundance or gone extinct. Moreover, the vast majority of larger scale studies describing biodiversity changes have not been able to clearly identify the environmental drivers of these changes. This demonstrates the urgent need to complement biodiversity time series with environmental data that can be used to explain the observed patterns but are often not collected on a routine basis (with exceptions in limnology and oceanography). The complexity of biodiversity dynamics that our results and those of other studies highlight is not to be interpreted as an obstacle to the development of conservation measures. Rather we argue that a better understanding of the patterns, trends and changes in biodiversity and environmental conditions will allow conservation measures to be better tailored to specific locations and taxonomic groups in order to significantly retard or even reverse further biodiversity loss 62 .

Methods
Data compilation. We circulated a call for biodiversity data within the ILTER network and additional partners to fill in geographical gaps. The criteria for data selection were: (1) each time series covers at least 15 years, (2) with preferably at least ten survey events during that time, (3)  Biodiversity data. Biodiversity data were expressed as abundance or biomass of surveyed taxa at each survey occasion, in some instances as percent coverage (e.g. for some time series of benthic algae and plants). Survey methods and season varied among time series, but were kept constant within each time series throughout the entire study period. In most cases, surveys were carried out once a year, but some time series had gaps or more frequent survey intervals (e.g. weekly resolution for phytoplankton, zooplankton and moths from the Finnish and Czech sites). In the latter case, we filtered the data in the time series to select only the months/seasons that were consistently surveyed throughout the whole study period, and we pooled the data (as sums or averages) within each year. Taxonomic resolution was kept constant for each time series throughout the entire study period, generally at the species or genus level, with a few exceptions (i.e. a few macroinvertebrates groups were identified at family resolution).
For each year within each time series, we computed four biodiversity metrics: total number of organisms or biomass (hereby referred to as abundance), taxonomic richness (i.e. number of taxa), Simpson´s diversity and temporal taxon turnover. The latter was computed as the proportion of taxa gained and lost between two subsequent years relative to the total number of taxa observed, using Eq. (1), as implemented in the R package codyn 63 .
Number of taxa gained þ number of taxa lost total number of taxa : We classified the time series into nine biogeoregions, three realms and eight taxonomic groups (Fig. 1).
Abiotic variables. For each study site we extracted the daily mean temperature and daily total precipitation data from the gridded observational dataset for precipitation, temperature and sea level pressure in Europe (spatial resolution: 0.25 degrees 64 ), and computed the mean annual temperature and total annual precipitation. For aquatic ecosystems we used in situ water temperature measured at the surface and calculated the mean temperature across the yearly monitoring periods. We gathered information on local anthropogenic pressures in a standardized questionnaire (Supplementary Methods). The questionnaire asked the data providers to assess the impact (from 1 = no to 4 = strong impact) of a series of pressures at the site (e.g. urbanization, sources of pollution, agriculture, etc.; for a full list see Supplementary Methods) and indicate whether impacts were constant or changed throughout the study period. Data providers were also asked to estimate the overall environmental quality of the site (i.e. naturalness; from 1 = low to 5 = high) and to state whether it was constant or changed throughout the study period. We used this information to define the quality-class of the sites, based on the overall assessment. The majority of sites (n = 67) scored 5 (i.e. high naturalness), 42 sites scored 4, 37 sites scored 3 and 15 sites scored 2.
Data analysis. We used a two-step procedure that allowed us to combine very heterogeneous original datasets (see paragraph Biodiversity data). First, we analyzed each time series separately to quantify time-series-specific biodiversity trends. Second, we used the effect sizes of the individual time-series-specific biodiversity trends to synthesize the overall trends and identify common patterns and drivers. The results of this second (meta-)analysis are reported in the paper.
For the first step of the analysis, we used the Mann-Kendall trend test to identify monotonic trends in each biodiversity and climate time series over the study period 65,66 . We detected serially correlated time series using auto-and crosscovariance and correlation functions 67 , and we applied the modified Mann-Kendall with Hamed and Rao 68 variance correction approach. We used Sstatistic and its variance as effect size of the trend 65 for the next step of the analysis. A similar meta-analytical approach has already been applied in a previous study on ecological time series 69 .
At 23 study sites, multiple time series were available for a given taxonomic group, e.g., when surveys were conducted at multiple transects or plots or with multiple traps. To avoid pseudoreplication, we combined those time series using meta-analysis mixed models (using the R package metafor 70 ) and extracted the cumulative effect sizes and their variances prior to the second step of the analysis. The total of 161 time series reported above refers to the aggregated final set of time series.
The second step of the analysis aimed at synthesizing the trends across the different time series. For that, we fitted meta-analysis mixed models to account for random effects. To compute the overall biodiversity trends and to explore how the trends varied among biogeoregions, realms and taxonomic groups, we included biogeoregion (nine levels), realm (three levels: freshwater, marine and transitional zone, and terrestrial) and taxonomic group (eight levels) as explanatory variables in the models with no intercept. We did so separately because biogeoregion, realm and taxonomic group were not independent (see Supplementary Table 1). The results of these models hence show whether or not the overall S-statistics for individual trends in the groups differ significantly from zero.
Additional meta-analysis mixed models were used to test the influence of selected abiotic variables describing site characteristics and climate on the biodiversity trends. We included the following variables: latitude, longitude, elevation, site naturalness, S-statistics 65 of temperature and precipitation trends, and the length of each time series. These explanatory variables showed little collinearity (|r| < 0.6), and were thus all retained as potential predictors. Similar to Everaert et al. 71 , we applied an information-theoretic approach to model selection and multi-model inference 72 , to determine the relative importance of those explanatory variables on the trends in biodiversity metrics. For this we created a candidate set of models with all possible linear combinations of explanatory variables and we extracted the corrected Akaike's information criterion (AICc) of each model using the R package glmulti 73 . We retained only plausible models with δAICc ≤ 2 of the best model (i.e. the one(s) with the lowest AICc 72 ) and computed the relative importance of each predictor variable as the sum of the Akaike's weights of all the selected models in which that variable was included. We computed the model-averaged coefficients (±95% C.I.) for each predictor variable in each selected model, weighted by the Akaike's weights (Supplementary Table 2). To evaluate the effect of the interaction among climatic and local stressors, we added to the selected models the interactions between site naturalness and temperature trends, and the interaction between site naturalness and precipitation trends. We then compared the resulting models (without interaction, with single interaction, with both interactions) and chose the one with the lowest AIC value. We have implemented this procedure because, to the best of our knowledge, it is not possible to include selected interactions into the information-theoretic approach to model selection and including all interactions would have resulted in an overly complex model.
To account for biases in biogeoregions, realms and taxonomic groups, we also performed sensitivity analyses (Supplementary notes and Supplementary Table 3), which confirmed that the results were robust.
Reporting summary. Further information on research design is available in the Nature Research Reporting Summary linked to this article.

Data availability
All datasets analyzed during the current study have been deposited in public repositories, at the links reported in Supplementary Data. Please note that some of the datasets will be publicly accessible after a period of embargo, from 1st January 2021. The source data underlying all figures are provided as a Source Data file. Source data are provided with this paper.
Code availability