Long-term data reveal unimodal responses of ground beetle abundance to precipitation and land use but no changes in taxonomic and functional diversity

While much of global biodiversity is undoubtedly under threat, the responses of ecological communities to changing climate, land use intensification, and long-term changes in both taxonomic and functional diversity over time, has still not been fully explored for many taxonomic groups, especially invertebrates. We compiled time series of ground beetles covering the past two decades from 40 sites located in five regions across Germany. We calculated site-based trends for 21 community metrics representing taxonomic and functional diversity of ground beetles, activity density (a proxy for abundance), and activity densities of functional groups. We assessed both overall and regional temporal trends and the influence of the global change drivers of temperature, precipitation, and land use on ground beetle communities. While we did not detect overall temporal changes in ground beetle taxonomic and functional diversity, taxonomic turnover changed within two regions, illustrating that community change at the local scale does not always correspond to patterns at broader spatial scales. Additionally, ground beetle activity density had a unimodal response to both annual precipitation and land use. Limited temporal change in ground beetle communities may indicate a shifting baseline, where community degradation was reached prior to the start of our observation in 1999. In addition, nonlinear responses of animal communities to environmental change present a challenge when quantifying temporal trends.

The world is run by the little things 1 , with insects being incommensurately under-described, yet comprising over 1 million of the 1.7 million named living species on Earth 2,3 . However, disproportionately few long-term ecological monitoring programs have targeted insects, and large knowledge gaps remain for many insect groups regarding their temporal trends and potential drivers. For example, BioTIME, a database of published biodiversity time series from 362 studies, contains only 22 studies (6%) primarily targeting insects 4 .
Ground beetle sampling. Ground beetle sampling followed standard procedures using either unbaited pitfall traps (either six or eight traps per site) 59,60 or both unbaited pitfall traps and hand collections 61,62 . In the Ruhr region, pitfall trapping was supplemented with a total of 204 hand captures following standardized procedures 61,62 . Beetles collected by hand captures were only used in analyses based on within-site comparisons and were excluded from analyses of activity density based on across-site comparisons. Sampling methodology was identical across years within each site but varied between regions. All sites (1) had an observation period (= study length) of at least ten years except one site with nine years (averaging 13.9 yrs ± 0.59 SE), (2) had at least four years of sampling within the observation period (averaging 7 yrs ± 0.33 SE), and (3)  www.nature.com/scientificreports/ the same season within the time series 63 . To maintain a minimum sample size of at least four sites per region, we included one site (in the Ruhr region) with an observation period covering nine years (and eight sampling years). All regions were sampled between May and July except the RMO region, which was sampled between August and September (accounted for in statistics). Captured individuals were stored in ethanol 64,65 and identified to species 66 .
Community metrics. For each sampling event (per site and year), we calculated 21 taxonomic and functional community metrics (Table 1) representing taxonomic diversity, taxonomic turnover, activity density, activity densities of functional groups and functional diversity.
Taxonomic diversity and activity density. To assess taxonomic diversity, we calculated five community metrics: species richness, Shannon's diversity 67,68 , Simpson's diversity 69 , evenness (Pielou's J) 70 , and 'Evar' evenness 71 . We opted to assess Evar in addition to Pielou's J as Evar is equally sensitive to rare and abundant species and therefore considered more appropriate than Pielou's J 71 , which on the other hand is more commonly used in other studies and therefore offers better comparability. To assess taxonomic turnover, we calculated four community metrics: turnover (number of species 'appearances' and 'disappearances' divided by the total number of species) and its components 'species appearances' (immigrations) and 'species disappearances' (extinctions) and its complementary species exchange rate (SERa) according to the formula provided in Hillebrand et al. 72 . SERa takes into account species proportional abundances and is therefore considered a more robust measure to detect compositional change through time compared to turnover 72 . The number of ground beetles captured with pitfall traps are considered measures of 'activity density' rather than true abundance as they reflect both the density of individuals and their level of locomotion facilitating their capture 51 .
Functional traits. We considered four trait groups representing seven functional traits 73 , hereafter functional groups: habitat preference (specialists and generalists; species inhabiting only one habitat class were considered specialists, all others generalists), flying ability (winged and dimorphic), hibernation stage (imago and larva), and feeding habits (predators and herbivores, only predators common enough to examine individual trends in activity densities). The selected traits address physiological and behavioral aspects of a species survival, particularly regarding land-use intensity (habitat preference) 42 , resource acquisition (feeding habits) 55 , dispersal and dissemination (flying ability) 56,57 , and resistance to harsh conditions (hibernation stage) 54 . Functional traits were retrieved from www. carab ids. org-an online database of ground beetle species traits 37 (accessed 29.07.2020).
Functional diversity. We calculated four distance-based functional diversity metrics based on Gower dissimilarity according to Villéger et al. 74    www.nature.com/scientificreports/ space, weighted by the relative abundance 76,78 ). Prior to analyses, species-specific traits (functional groups) were fuzzy-coded between 0 and 3 following Chevenet et al. 79 . Dimensionality reduction was required during the calculation of functional richness; the final quality of the multidimensional trait-space was 0.61.
Climate data. We retrieved daily temperature and precipitation data from the gridded observational European dataset (E-OBS Temperature and Precipitation data set) 80,81 with a spatial resolution of 0.1 degrees. We then calculated the mean annual daily temperature and annual cumulative precipitation of the 12 months preceding each sampling event to account for climatic conditions across the full range of ontogeny 82 .
Land use data. We extracted land use data around each sampling site from the CORINE land cover dataset 83 , with the highest resolution (5 ha) available from the year 2012. CORINE land cover data provide spatial coverage of land use types and were available for the years 2012 and 2018. We selected the year 2012 as it represents the mean of the timespan covered by our time series. Land use types were assumed to be consistent during the observation period as all sites are situated in biosphere reserves or along river floodplains. As an indicator for land use intensity, we calculated the land use index (LUI) 84 based on the coverage of land use categories within a radius of 500 m (LUI_500), 1000 m (LUI_1000) and 2000 m (LUI_2000) around ground beetle sampling sites.
The three different buffers of land use intensities were all highly correlated (R 2 = 0.88, p < 0.001). Hence, we only included LUI_1000 as site-specific estimates of LUI.

Data analysis. Standardization of community metrics.
To avoid bias in richness metrics (taxonomic and functional diversity) due to variable exposure times of pitfall traps within sites over time, we used rarefied indices. All richness-related community metrics were rarified within a given site to the year with the shortest exposure time 85,86 . Rarefied community metrics are only used in analyses based on within-site comparisons. We refrained from applying rarefaction across sites to prevent rarefaction-induced bias of richness-related community metrics at sites with shorter exposure times (i.e., substantially reducing exposure time and thus removal of rarely captured species) and because the shape of rarefaction curves may vary across sites. LUI = %pasture + 2 x %arable land + 4 x %urban area Table 1. Community metrics, regions and their abbreviations evaluated in this manuscript. 1 = taxonomic diversity, 2 = taxonomic turnover, 3 = activity density 4 = functional traits, 5 = functional diversity. www.nature.com/scientificreports/ Activity density was standardized by dividing by the total exposure time of pitfall traps per site and year (units are the number of captured beetles per day) and used in analyses of across-site comparisons. Analyses based on within-sites comparisons also include hand captures (Ruhr region only) as those had a standardized (equal) sampling effort in each sampled year within the sites.
Statistical approach. We used two approaches to examine climate and land-use effects on ground beetle communities. The first approach focuses on temporal trends and their drivers and is based on within-site comparisons. Approach one is independent of local scale environmental heterogeneities between sites. The second approach takes advantage of temporal variation in community metrics and climate variables (of each sampled year at each site) in addition to spatial variation. Approach two examines the effects of climate and land use on ground beetle communities, but it is limited to changes in activity densities, which were comparable across sites. We summarize both statistical approaches briefly below and provide more details in the supplementary materials.
Approach one uses a meta-analytical approach based on Mann Kendall trend tests of time series within each site. It accounts for both temporal and spatial autocorrelation and tests for (1) overall and regional temporal trends in all ground beetle community metrics, (2) overall trends in temperature and precipitation, and (3) the effects of trends in climate variables and site-specific LUI's on the overall trends in community metrics. Following Maire et al. 87 , we report the model coefficients from meta-analytical mixed effects models as trend mean effect sizes (TMES) and standard errors (SE). Observation period was additionally included in models to account for time series length. The major advantage of this approach is the ability to assess overall and regional trends based on heterogeneous site-specific time series (covering heterogeneous habitat types) for all community metrics.
With the second approach, we investigate whether site/year-specific variation in climate and site-specific variation in land use intensity explain variation in ground beetle activity densities. This approach tests for relationships between the drivers of site/year-specific temperature, precipitation, site-specific LUI's and the responses of activity density, functional groups, and the dominance of functional groups (using axes of a Principal Component Analysis [PCA] to eliminate multicollinearity, appendix Fig. A1 & Table A1) using generalized least squares regression (gls) models. The major advantage of this approach is the ability to incorporate both temporal and spatial variation in predictors and responses. However, this approach could only be conducted for community metrics based on activity density as only activity density could be realistically standardized across sites.
Statistical software. All analyses were conducted in R 3.6.1 88 . Species richness, Shannon and Simpson diversity indices, evenness (Pielou's J), turnover, species appearances and disappearances, species rarefaction, and PCA were calculated using package vegan, version 2.5-5 85 . The Evar metric was calculated using package microbiome, version 1.8.0 89 . Functional diversity metrics were calculated using the FD package, version 1.0-12 90 . Package metafor, version 2.4-0 91 was used for meta-analyses and package glmulti, version 1.0.8 for multi model selection and inference 80 . Generalized least squares models were calculated using the package nlme, version 3.1-140 92 .

Results
A total of 65,874 ground beetles from 194 species were captured in 46,566 days of pitfall trapping and 204 hand captures. Rarefaction was based on a total exposure time of 40,916 pitfall trapping days and hand captures, reducing the initial dataset to estimates of 59,548 ground beetles from 190 species. Analyses of spatio-temporal variation in community metrics based on activity density included 61,826 ground beetles from 184 species from pitfall traps only.
Climate trends and land use gradient. Trends in mean temperature and cumulative precipitation over the 12 months prior to sampling indicated an increase in temperature (TMES ± SE: 11.07 ± 4.61; P = 0.016) and a marginal increase in precipitation (4.46 ± 2.40; P = 0.063) across all sites during the observation period (Appendix Fig. A2). Across all 40 sites and sampling years, annual mean temperatures varied from 7.8 to 12.6 °C and annual cumulative precipitation ranged from 424.5 to 878.4 mm (Appendix Fig. A3). The index of land use intensities (LUI) at buffers of 1000 m around sites averaged 144 ± 88 (mean ± SD, range: 10-375, Appendix Fig. A4).

Temporal trends in community metrics and their drivers.
On average across the 40 sites, none of the 21 assessed community metrics changed significantly over time (Fig. 2). Within regions, APP decreased (p = 0.007) whereas DIS increased (p < 0.001) in the CHO region and SERa increased (p = 0.002) in the SPW region (Appendix Fig. A5). No significant trends were evident in any of the other community metrics and the five regions (Appendix Figs. A5, A6, A7 & A8).

Drivers of spatio-temporal variation in activity density and functional groups.
Activity densities from pitfalls were higher at the three northeastern regions (CHO, Elbe, SPW) than at the two southwestern regions (RMO, Ruhr; Fig. 3A). Functional group activity densities tended to parallel total activity densities, with the notable exception of much higher captures of predators in the Ruhr region (Appendix Fig. A9). Total activity density had a unimodal response to both precipitation (peaking at around 600 mm; www.nature.com/scientificreports/ land use intensity (peaking at around LUI ~ 160 Table 3, Fig. 3C). All seven functional groups showed similar unimodal responses to LUIs (Appendix Table A2 & Fig. A10). The activity densities of three functional groups were driven by annual precipitation (Appendix Table A2, Fig. A11). Activity densities of generalists and winged beetles were highest at ~ 600 mm of annual precipitation and then declined at higher precipitation rates. Predators increased linearly with precipitation. Three functional groups: dimorphic beetles, larval-, and imago-hibernators had significant second order polynomial responses to temperature (Appendix Table A2E,F,G), though the fitted relationship was more indicative of linear declines in activity densities with temperature (Appendix Fig. A11D,E,F).
Ground beetle functional group dominance varied with land use intensity and annual precipitation. The first principal component (PC1) of the functional dominance PCA represents communities with high proportions of dimorphic beetles and low proportions of specialists, winged beetles, and imago-hibernators (Appendix ,  Table A1). PC1 had a u-shaped response to land use intensity, indicating specialists, winged beetles and imagohibernators drove the pattern of increasing activity densities at intermediate land use intensities (Appendix Fig. A12). The second component (PC2) represents communities with high proportions of larval-hibernators and low proportions of predators. PC2 was negatively correlated with precipitation, indicating predators had increased dominance and larval-hibernators were proportionally less common with increasing annual precipitation (Appendix Fig. A12).

Discussion
We did not detect any net directional temporal changes in ground beetle communities across the 40 sites over the last two decades in Germany. However, across all sites and sampling years, ground beetle activity densities peaked at intermediate annual precipitation and land use intensity, a result supported by other studies reporting compositional shifts across land-use gradients 45,46,[93][94][95] . The unimodal response of activity densities to precipitation was most pronounced for habitat generalists and winged beetles, which may be indicative of better dispersion abilities of these groups to more suitable habitats 49 following exceeded tolerance limits to precipitation 31 . Increases in activity density and functional dominance of predators with precipitation may be a consequence of rainfall-driven increases in habitat volumes resulting in extended food chain length 96 .
In contrast to our expectations of increased taxonomic and functional diversity (Hypothesis 1), we detected no temporal trends in the 21 analyzed community metrics representing taxonomic and functional diversity,  Table 1. Table 2. Drivers with a significant influence on the overall trends in community metrics. Arrows indicate the direction of effect sizes (z values) of influential drivers identified in the model selection and multi-model inference procedures. Bold font indicates significant drivers and their effect sizes. Abbreviations of community metrics are explained in Table 1. www.nature.com/scientificreports/ taxonomic turnover, activity density, and activity densities of functional groups. However, we cannot reject our prediction that trends in functional diversity paralleled trends in taxonomic diversity as both taxonomic and functional diversity were fairly invariant over time. Our results are in line with other large-scale biodiversity assessments of various taxonomic groups reporting few long-term net changes in taxonomic diversity metrics 12,16,72,97 . These studies have generally reported overall increases in temporal turnover, whereas in our study, turnover changes were restricted to two regions (CHO and SPW, both in northeast Germany). Our results Figure 3. Spatio-temporal variation of activity densities across the five regions (A) and activity density responses to significant drivers from the overall model (Table 3): Precipitation (B) and land use intensity (LUI) around 1000 m of sampling sites (C). Error bars in Panel A represent one standard error. Points in Panels B and C represent activity densities within each site and year. Table 3. Spatio-temporal drivers of overall ground beetle activity density. The generalized least squares model included an autoregression term to account for temporal autocorrelation, a grouping factor of site nested within region to account for repeated sampling within sites. Estimates are shown in bold when significant. Second order polynomial terms are denoted as "2nd poly". www.nature.com/scientificreports/ also contrast several local scale studies from the UK, Northern Germany and the Netherlands reporting declines in ground beetle taxonomic diversity and abundance 47,48,50 , and to recently reported declines in activity density and species richness in ground beetles at a single German site 98 . However, these studies examine earlier time periods, starting between 15 to 30 years before the time series assessed here. We speculate that this may indicate a shifting baseline effect, whereby the ground beetle communities investigated here may have been altered prior to the start of our observation period, with only taxa robust against anthropogenic disturbances remaining. In partial support of Hypothesis 2, taxonomic turnover (SERa) increased with temperature and tended to increase with precipitation across all sites. Regionally, we found an increasing species exchange rate in the SPW region, while in the CHO region, species disappearance decreased and species appearance increased over time. These results illustrate that long-term ecological trends can vary across spatial scales 99 . Although research focusing on single 50,98 or few sites 40,47 are indispensable to understand the influence of local scale drivers, our results illustrate that localized changes in diversity may not hold at broader spatial scales 97 .

Community metric Predictors
Precipitation was the most important environmental driver of site-based trends in ground beetle communities. Activity density of habitat specialists increased with increasing precipitation while Simpson's diversity, evenness (EVE), and functional dispersion, all declined. This suggests that habitat specialists became the dominant ground beetle group in moist conditions, at the expense of habitat generalists. High rainfall increases ground beetle habitat specialization, a consequence of specific adaptations to moist conditions 56,57 . Site-based trends in the species exchange rate (SERa) increased with site-based trends in temperature and declined with land use intensities. Consequently, variable environmental conditions at the investigated sites may increase selection for species that are better adapted to more extreme living conditions [40][41][42]52 . Accordingly, high land use intensities might lower the exchange rate of species. However, our meta-analytic approach, while necessary for a standardized assessment of all examined responses besides activity densities, collapses climate and ground beetle response values from each site into one overall trend, and thus has reduced power to detect the influence of environmental drivers, especially if sites have high variability across years.
Considering temporal variation in addition to spatial differences revealed climate-driven influences on the overall activity density of all ground beetles and on activity densities of individual functional groups that were not evident from our site-based meta-analytic approach. High activity densities at intermediate precipitation suggests moderate rainfall either provided more resources for ground beetles or increased ground beetle movement relative to extreme rainfall conditions, in turn increasing trap catch. While temperature did not have any influence on the overall activity density, three functional groups (dimorphic beetles, imago hibernators, larval hibernators) tended to decline with increases in annual temperature. This suggests decreasing overwintering survival as a consequence of desiccation 50 and a decline of beetles with reduced abilities to relocate to habitats that are within their temperature tolerance range. Microhabitat conditions and local habitat structure can also drive ground beetle activities at smaller spatial scales 36,100 . The sampling sites assessed in our study covered many heterogeneous habitat types and corresponding microclimatic conditions. Disentangling how interactions between these conditions contributed to the observed and unobserved influences of climate on activity density of functional groups is a fruitful avenue for future work.
We did not find support for our third hypothesis that LUI would have a linear negative effect on temporal trends in taxonomic and functional diversity. However, across an increasing gradient of LUI, ground beetles varied non-linearly both in activity density and by functional groups. Activity densities of all functional groups had unimodal responses to LUIs. Ground beetle taxonomic reorganization (SERa) was highest at low LUIs, potentially due to a larger species pool. At intermediate LUIs, specialists, winged beetles, and imago-hibernators proportionally dominated the ground beetle communities, while very high land use intensities reduced activity densities of ground beetles independently of their functional traits.
Our study is subject to several common challenges in long-term research that may limit our ability to capture trends in the sampled ground beetle communities. While all but one of our 40 sites cover periods of at least 10 years (average 13.9 years), a caveat of our temporal analysis is the comparably lower number of sampling years ranging from 4 to 10 (averaging 7 years). This particularly applies to the Elbe region, which had the fewest sampling years. However, regional trends in the Elbe region were similar to those from the other four regions, indicating that overall, our analyses reflect well the temporal patterns in the ground beetle communities across the last two decades. Additionally, the peak of human pressures may have been reached before the onset of our observation period 24 starting in 1999. The sampled community was potentially already at the "bottom of the barrel"-that is major long-term changes in ground beetle communities prior to 1999 may have filtered out particular taxa resulting in no evidence of recent trends 52,101 .
In comparison to previous long-term studies on ground beetle communities 40,41,98 , a major strength of our study is its high number of sites from multiple regions. While heterogeneous habitat types covering different local environmental conditions might result in averaged-out trends at larger scales, similar patterns across all five regions suggest rather low temporal variation in the assessed ground beetle communities. However, and importantly, responses of functional group activity densities had several unimodal or u-shaped responses to climate and land use. Such non-linear responses challenge trend detection, especially when sample number is limiting. Additionally, incorporating both local scale microclimate and habitat structure poses a challenge in long-term and large scale-studies of invertebrates, but remains a key consideration for future work. Finally, we echo recent calls for distributed and standardized long-term monitoring schemes to unravel temporal changes in biotic communities and their driving forces 102  www.nature.com/scientificreports/