Present and historical landscape structure shapes current species richness in Central European grasslands

Current diversity and species composition of ecological communities can often not exclusively be explained by present land use and landscape structure. Historical land use may have considerably influenced ecosystems and their properties for decades and centuries. We analysed the effects of present and historical landscape structure on plant and arthropod species richness in temperate grasslands, using data from comprehensive plant and arthropod assessments across three regions in Germany and maps of current and historical land cover from three time periods between 1820 and 2016. We calculated local, grassland class and landscape scale metrics for 150 grassland plots. Class and landscape scale metrics were calculated in buffer zones of 100 to 2000 m around the plots. We considered effects on total species richness as well as on the richness of species subsets determined by taxonomy and functional traits related to habitat use, dispersal and feeding. Overall, models containing a combination of present and historical landscape metrics showed the best fit for several functional groups. Comparing three historical time periods, data from the 1820/50s was among the most frequent significant time periods in our models (29.7% of all significant variables). Our results suggest that the historical landscape structure is an important predictor of current species richness across different taxa and functional groups. This needs to be considered to better identify priority sites for conservation and to design biodiversity-friendly land use practices that will affect landscape structure in the future.


Introduction
Semi-natural grasslands are among the most species rich ecosystems in Central Europe (Habel et al. 2013). Since those grasslands are part of the cultivated landscape in Germany, abandoning them would lead to scrub encroachment and eventually to forests, culminating in the loss of culturally important and biologically diverse ecosystems. Over the last decades, European grasslands, their communities and related ecological processes have been facing various threats and challenges, such as management adaptations to changing socio-economic conditions, leading to intensification or abandonment. A consequence of intensifying fertilization practises on neighbouring fields is a decrease in soil moisture and an increase in soil nutrient availability. (e.g. Biró et al. 2013;Diekmann et al. 2019). Thus, species adapted to low-intensity land use over long time periods (Pärtel et al. 2005) struggle with these altered conditions, which now becomes evident through ongoing declines in diversity, abundance and biomass in grasslands (e.g. Seibold et al. 2019).
The general importance of present landscape structure (e.g. Sirami et al. 2019) and local land-use intensity (e.g. Allan et al. 2014) for the species richness of plants and arthropods are well studied. However, different spatial scales might be important for different taxonomic and functional groups (Honnay et al. 2003). In our study, landscape structure is comprised of landscape composition, i.e. the number and abundance of different land covers, and landscape configuration, i.e. the spatial arrangement of these land covers within a landscape (With 2019). Both can be quantified using landscape metrics (Walz 2011). Plant and arthropod species are suitable target groups for evaluating the relative importance of present and historical landscape structure for biodiversity, because of their essential contribution to many ecological processes and the general differences in their ability to react or be affected by disturbances at different spatiotemporal scales (e.g. active and passive dispersal), as well as their crucial role in biotic interactions such as plant-pollinator relationships (e.g. DiLeo et al. 2014). Plants with low dispersal ability, for instance, are mainly influenced by local habitat structures, whereas long-distance dispersed plant species are influenced more strongly by landscapescale structures (Damschen et al. 2014;DiLeo et al. 2014;Trakhtenbrot et al. 2014). Li et al. (2020), however, found the species richness of weak and strong dispersing species to be negatively affected by the area-weighted mean of shape metric and the nearest neighbour distribution at larger spatial scales. Complex local structures were found to be favourable for both dispersal types. For arthropod species, Seibold et al. (2019) found the biomass of weak dispersing species to decline with increasing cover of arable fields in the surrounding area, while the biomass of strong dispersing arthropod species was not affected by surrounding arable field cover. When Liu et al. (2014) investigated how landscape structure and local land-use intensity affected functional beetle diversity, they found herbivore, predator and decomposer abundances to be affected by a combination of landscapeand local-scale variables. At the same time, the negative effects of local land-use intensification on species abundance can be mitigated by increasing landscape structure diversity (e.g. Perović et al. 2015;Neff et al. 2019;Wintle et al. 2019).
Most studies focus only on the effects of presentday land use and landscape structure on species diversity, disregarding historical land management changes that may have altered both local land cover and the structure of the surrounding landscape. However, long lasting effects of historical land cover and landscape structure from as long as 150 years ago could still influence present-day species communities Lindborg and Eriksson 2004). For example, some studies showed that the current species richness of grasslands is influenced by past habitat amount (Finnish meadows, Raatikainen et al. 2018) as species show a time lagged response to habitat fragmentation caused by past management decisions (Fahrig 2003;Purschke et al. 2014) and the decrease of habitat amount (Fahrig 2017;Riibak et al. 2020). These examples show that changes in species richness are ongoing and may be strongly intertwined with historical changes of landscape structure. The extent and time frame of impacts of the historical land cover on current species richness are poorly understood. In an agricultural landscape, using data that extended 60 years into the past, Riibak et al. (2020) found that long-distance dispersing plant species occur more often in areas in which neighbouring patches were historically more persistently covered with grassland. Some calcareous grassland plant species may even be absent from suitable isolated patches due to low dispersal ability across unsuitable matrix habitats over extended time periods (Riibak et al. 2015). Thus, we also expect historical landscape structure to affect species richness in our study.
In this study we examined the effects of present land-use intensity, historical and present land cover and landscape structure on present-day plant and arthropod species richness in semi-natural grasslands of Central Europe. In particular, we analysed how strongly arthropod and plant communities are influenced by historical landscape metrics from three distinct points in time (1820/50, 1910/30, and 1960) relative to present local land-use intensity and landscape metrics measured at the time of plant and arthropod sampling (2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015)(2016). To better understand what shapes present-day species richness, we additionally analysed species richness of different subsets of plant and arthropod species. These were defined based on taxonomic categories as well as functional traits (i.e. dispersal ability, habitat and feeding specialization, guild, stratum use), which should help improve the understanding of the mechanisms underlying the impact of landscape structure on species' communities (Weiher and Keddy 1995).
We expect species richness to increase with increasing historical patch area the same way species richness is affected by present patch area as shown in other studies (e.g. Dembicz et al. 2020). Furthermore, we expect temporal habitat continuity to have a positive effect on species richness (Veldman et al. 2015;Le Provost et al. 2021). To investigate these expectations, we focused on the following research questions: (a) Do historical land cover and historical landscape structure affect present species richness across the studied grasslands and what is the relative importance of historical landscape structure compared to present landscape structure. (b) Which historical time periods are most important for the species richness of plants and arthropods? (c) Which functional groups of plants and arthropods are most strongly affected by historical land use and landscape structure?

Study area
The study was conducted within the framework of the Biodiversity Exploratories, comprising three production landscapes in Germany (Fischer et al. 2010). The sampling areas are located in the UNESCO Biosphere Reserve Swabian Alb (Schwäbische Alb, 48°20 0 28 00 -48°32 0 02 00 N, 9°10 0 49 00 -9°35 0 54 00 E, * 422 km 2 , 460-860 m asl) in south-western Germany, the National Park Hainich and its surrounding areas (Hainich-Dün, 50°47 0 25 00 -51°22 0 43 00 N, 10°10 0 24 00 -10°46 0 45 00 E, * 1300 km 2 , * 285-550 m asl) in central Germany, and the UNESCO Biosphere Reserve Schorfheide-Chorin (52°47 0 25 00 -53°13 0 26 00 N, 13°23 0 27 00 -14°08 0 53 00 E, * 1300 km, 3-140 m asl) in north-eastern Germany (for a map of the regions in Germany see Online Resource 1). In each region, 50 grasslands managed as meadows, pastures or mown pastures were chosen from 500 previously surveyed grid plots, following a stratified random selection (for details see Fischer et al. 2010). The plots were restricted to the two most common soil types in each region, thereby minimizing possible confounding effects in the data. The grasslands were designated as such for at least 10 years before the start of the project in 2006 and represent the regional-specific range of land use. On each grassland, a 50 m 9 50 m plot was established within a larger management unit for detailed species surveys. Details about the grassland management have been collected yearly since 2006 (Vogt et al. 2019).

Species richness and functional groups
Plant species richness was calculated from abundance data collected on 148 of the 150 Biodiversity Exploratories experimental grassland plots. All plant species were sampled in a 4 m 9 4 m subplot once in spring of 2009 (Socher et al. 2012). Of 347 recorded plant species, 304 were used for the analyses (Table 1). Highly uncharacteristic (e.g. Pinus sp.) and unidentified species (only identified to genus level) were excluded. Arthropod species richness was calculated from yearly sweep-net samples (20 double sweeps along three 50 m transects, i.e. the plot borders, once in spring and once in summer of each year) of all 150 grassland plots, cumulated over the years 2008 to 2015 to get more robust community data. In our analysis we focused only on the orders Hemiptera (restricted to Heteroptera and Auchenorrhyncha), Coleoptera and Araneae from these samples. Of 1214 recorded species, 1044 species were used in the analyses ( Table 1). Species that were highly atypical for grassland habitats were excluded for all arthropod related data . Additionally, to represent ground-active arthropods, Araneae and Coleoptera species caught with funnel pitfall traps of 15 cm diameter (filled with 3% copper sulphate solution) in the Hainich region were analysed. Three traps were installed on 43 plots of the 50 plots and were emptied monthly from May to October 2008. Two out of three traps were randomly selected per month and further analysed due to trap losses. Individuals were pre-sorted to taxonomic orders in the lab and later identified by taxonomic experts (see Acknowledgements). Of 236 recorded species, 228 species were used in the analyses (Table 1). Again, atypical species were excluded .
Besides testing the responses of total plant and arthropod species richness, we analysed species richness for subsets of plant and arthropod communities based on taxonomic (Araneae, Coleoptera, Hemiptera) and functional groups derived from a set of functional traits. For plants, the functional groups were based on species' dispersal abilities, i.e. short, medium, longdistance dispersers, derived from Hodgson et al. (1995), and whether they are mesic grassland species typical of meadows and pastures (phytosociological class Molinio-Arrhenatheretea), or dry grassland species, based on their phytosociological affiliation after Ellenberg et al. (1992). For arthropods, functional groups were based on species' dispersal ability using wing development and flying ability for Coleoptera and only wing development for Hemiptera, as well as on activity ranges and dispersal strategies for Araneae (weak, moderate, strong dispersers) as proxies. In addition, species were classified according to their feeding guild (carnivorous, herbivorous, omnivorous) and the main vegetation stratum usage (ground, herb layer), and herbivores were further classified according to their food specialization (polyphagous, oligophagous, monophagous). Trait classification was based on the data published by Gossner et al. (2015). For pitfall trap data, monophagous, oligophagous, herb layer dwelling, and herbivorous functional groups could not be analysed separately due to the low number of species.

Present local land-use-intensity
Present land-use intensity was quantified through a continuous land-use-intensity index (LUI, Table 2), which was calculated for each grassland plot based on main grassland management components according to Blüthgen et al. (2012): mowing (frequency per year), grazing (livestock units 9 days of grazing ha -1year -1 ) and fertilization (kg total nitrogen ha -1year -1 ). To calculate the LUI, mowing, grazing and fertilization intensity were standardized in relation to its mean depending on the region and then added up.
To calculate the compound LUI index over several years, the above components are first summarised, where for each management type the years are added up after standardising by the mean of management intensity of all focal years, divided by the number of years for each component. Finally the standardised values for the three management types averaged over the years were added and the result was square root transformed to achieve an even distribution (Blüthgen et al. 2012). The intensity of these different components was recorded using standardised questionnaires for the land owners (Vogt et al. 2019). LUI was calculated using the LUI calculation tool (Ostrowski et al. 2020) implemented in BExIS (http://doi.org/10. 17616/R32P9Q). In our modelling approach we calculated the LUI over the years 2006 to 2015 for the sweep net sampled arthropod species richness, 2006 to 2009 for the plant species richness, and 2006 to 2008 for the pitfall trap sampled arthropod species richness. We calculated the mean intensity beginning two or three years (depending on data availability) prior to the arthropod and vegetation sampling, because land-use intensity of the previous years likely affected the species communities. Historical and present landscape metrics For each region, historical maps were available for one period in the nineteenth century and two periods in the twentieth century. Additionally, we used present-day maps. For each of the four time periods, we calculated a set of landscape metrics. Historical maps from the Swabian Alb were digitized from cadastral maps (''Flurkarten'') of 1820, topographic maps (TK25; 1:25,000) from the German Empire of 1910, and topographic maps (TK25, 1:25,000) from the Federal Republic of Germany (BRD) of 1960. Historical maps from the Hainich were digitized from topographic maps (Urmesstischplatten, 1:25,000) of 1850, topographic maps (TK25, 1:25,000) from the German Empire of 1930, and topographic maps (TK10, 1:10,000) from the German Democratic Republic (DDR) of 1960. Historical maps from Schorfheide-Chorin were digitized from topographic maps (Urmesstischblatt, 1:25,000) of 1850, topographic maps (TK25, 1:25,000) from the German Empire of 1930, and topographic maps (TK25, 1:25,000) from the German Democratic Republic (DDR) of 1960. Present landcover features for each plot and the surrounding 2 km were recorded and mapped in 2008 in the field and polygon borders were defined from high resolution (40 cm) aerial photographs taken in 2009 to increase the accuracy of the maps. Polygon borders did not vary between 2008 and 2009 Perović et al. 2015).
Using a Geographical Information System (QGIS v. 2.18), every polygon of the historical and present land-cover layers was associated to one of the following land cover classes: grassland, forest, arable, and other according to the available information. For the analysis of the subset of dry grassland plant species, grasslands were further separated into dry and wet grasslands and only dry grasslands were used for modelling (for a comprehensive list of both land-cover keys and their original names and simplifications, see Online Resource 2).
The use of landscape metrics to quantify landscape structure has been well established (e.g. Turner and Gardner 2015) and is widely available as software packages (e.g. FRAGSTATS: McGarigal and Marks 1995, R Statistics: R Core Team 2019). Landscape metrics of the class and landscape spatial scale (Table 2) were calculated for each combination of plot patch and buffer zone (100 m, 250 m, 500 m, Index combining data on grazing, mowing and fertilization None Patch 1000 m and 2000 m around the plot centres) and time period (1820/50, 1910/30, 1960, 2008). Permanency and 'variation in field land use' (Table 2) were calculated for each buffer zone.

Statistical analysis
Using the R Statistics Software v. 3.6.1 (R Core Team 2019), we investigated the impact of present and historical land cover and landscape structure on plant and arthropod species richness. We fitted separate generalized linear models (GLMs) for the richness of each subset of species (response variables, i.e. for plants: species richness of all, short, medium, and long dispersing, dry and mesic grassland species; for arthropods: species richness of all, weak, moderate, strong dispersing, carnivorous, herbivorous, omnivorous, polyphagous, oligophagous, monophagous, ground, and herb layer specialist species) with the landscape variables as predictors (Table 2), as well as the regions (except for pitfall trap models) as covariates (base models; Eq. 1).
All continuous predictor variables were z-transformed to be able to compare effect sizes. To select the most important buffer radius and year of each landscape metric for use in the base models, GLMs were computed beforehand using only one variant of a certain landscape metric at a time, plus the plot patch area and the covariates LUI and region as predictors. We then selected the best fitting combinations of buffer radius and year of each landscape metric for each subset of species according to the models with the lowest Akaike Information Criterion value.
The final models were produced using model selection (step function in base R). We calculated effect sizes of the predictor variables on the original scale of species richness by raising e to the power of the estimate, which provides the ratio of change in species richness after a change of the predictor by one standard deviation (SD) unit, since the data were z-transformed.
We used Poisson distribution in the GLM, if possible, or negative binomial distribution (function glm.nb from the R package 'MASS' (Venables and Ripley 2002)) in cases where the response variables showed overdispersion. All models were fit with loglink function. Model validation was performed following Zuur et al. (2013). We checked linearity of the relationship between predictors and response variables using scatterplots and Pearson residuals. Further, collinearity was tested observing variance inflation factors (VIF). Variables with VIF-values higher than 5 were excluded (see Online Resource 3) Some significant correlations were observed among the predictor variables, particularly among the variants of the same landscape metrics (see Online Resource 4). However, there was no collinearity in the final models. Finally, we tested for spatial dependence of the model residuals by means of Moran's I and found no case of significant spatial auto-correlation (see Online Resource 5). Various other R packages were used for geographical data analysis, data wrangling and plotting (see Online Resource 6).
The significance of effects of predictor variables in the final models was tested with type-II Likelihood Ratio Tests using the 'Anova' function from the R-package 'car' (Fox and Weisberg 2019). Additionally, we tested whether the final models that included historical landscape metrics fit better to the data compared to models containing only present landscape metric variables using Vuong tests (Vuong 1989) in order to validate the significance of historical landscape effects.

Landscape change
Overall grassland area has increased since the nineteenth century by 7% in the Swabian Alb, by 59% in the Hainich, and by 24% in Schorfheide-Chorin, accompanied by an increase in forest area (25-36%) and a considerable decrease of arable land (38-59%) across all regions (Fig. 1). On average, all landscape metrics, i.e. patch density, Shannon evenness, proportional grassland area, edge density and shape index, slightly increased over the years, except for nearest grassland neighbour distance. The latter remained relatively constant in the Swabian Alb and the Schorfheide-Chorin but continuously decreased by 67.2% in the Hainich (Fig. 2).

Plant species richness
Of the 20 selected landscape metrics included in the final models of all the different plant species functional groups, 14 were historical variants. For shortand long-dispersing as well as mesic grassland plant species, models with mixed sets of historical and present landscape metrics fitted the data better compared to models consisting of only present variants (Vuong test, p \ 0.05). For the other functional groups (short distance dispersers and dry grassland species), there was no significant difference among 'mixed' nor 'present' models (Vuong test,p [ 0.05). The time periods most represented in the selected historical landscape variables were 1820/50 (29.7% of all significant variables) and 1960 (22.8%).
Among the functional groups, which were significantly more affected by a mixture of present and historical variables compared to models using only present variables, species richness of plants with short distance dispersal was negatively associated to the landscape patch density of 1820/50 in a radius of 250 m (p \ 0.01) and by the proportion of grassland area in the surrounding 2000 m of 1820/50 (p \ 0.01; Fig. 3). Furthermore, richness of long-distance dispersing plants was positively affected by the landscape patch density of 1820/50 in a radius of 2000 m (p \ 0.05). The proportional area of grasslands in the surrounding 2000 m of 1960 had a negative effect on long-distance dispersers (p \ 0.001). Mesic grassland specialists were positively correlated to the Shannon evenness of 1960 calculated for the surrounding 1000 m (p \ 0.001), but negatively to grassland permanency (p \ 0.001). As an example, the ratio of change in species richness for mesic grassland specialist species per change of one standard deviation unit (SD) for the Shannon evenness (e estimate ) was estimated as 1.097, which translates into a higher plant species richness of around 10% in grasslands with historically higher evenness in the surrounding landscape compared to grasslands with historically less even surrounding landscapes. Likewise, a change in grassland permanency by 1 SD would elicit a decrease in species richness of around 11%.

Arthropod species richness (sweep-net samples)
For the arthropod species collected via sweep-net sampling, 69% of all 58 landscape metric combinations selected in the final models were historical variants. For the herbivores (Vuong test, p \ 0.01), weakly and moderately dispersing species (Vuong test for both, p \ 0.05), as well as Coleoptera (Vuong test, p \ 0.05) and Araneae species (Vuong test, p \ 0.01) the models including a mix of historical and present landscape metrics showed a significantly better fit compared to models with exclusively present metrics (Fig. 4). The best performing models included historical landscape metric variants of different periods but the time period of 1960 was most often selected among the significant landscape metric variables of all sweep netted arthropod functional groups.
Species richness of herbivorous arthropods was positively affected by the landscape patch density of 1820/50 in a radius of 2000 m (p \ 0.001) and the nearest neighbour distance of 1910/30 within a radius of 2000 m (p \ 0.01). Species richness of arthropod

Arthropod species richness (pitfall traps)
For the pitfall trap samples of ground-active arthropods (Fig. 5), around 78% of the 63 selected landscape metrics were historical variants. Models of species richness of omnivorous species, and species with either low or high dispersal ability showed a significantly better fit (p \ 0.01) with a mixed selection of historical and present variables compared to models with only present variables. The time period most frequently present in the final models was 1820/50.
Omnivore richness was positively correlated with the proportional grassland area of 1820/50 within a radius of 500 m (p \ 0.001) and the shape index within radius of a 2000 m in 1910/30 (p \ 0.01). Species richness of weak dispersers positively correlated with the edge density within a radius of 500 m in 1910/30 (p \ 0.01) and negative correlated with the patch density within 500 m of 1820/50 (p \ 0.05). Species richness of strong dispersers was positively correlated with the Shannon evenness within 500 m of 1960 (p \ 0.05) and the edge density within 500 m of 1820/50 (p \ 0.01; see Online Resource 7 for a table with modelling parameters and statistical results for all models).  blue:1910/30, green: 1960, purple: 2008 Discussion Historical landscape structure strongly affected the present species richness of plants and arthropods in the studied grasslands in the three regions in Germany. Our study revealed that many of the tested plant and arthropod taxonomic and functional groups were not affected by present or historical metrics alone, but by a combination of variables at different time scales. This suggests that species richness cannot exclusively be explained by present landscape conditions, and historical landscape structure of the past 200 years is highly relevant for the species richness that we observe today. However, the impact of the time period and landscape metric differed for each functional group. Therefore, it is necessary to clearly evaluate the effects of the historical landscape from the point of view of the focal communities.

Effect of historical landscape structure on current species richness
The historical grassland area at the local patch level, i.e. the size of the grassland in which the plot is located, had no influence on the local richness of the different subsets of plant and arthropod species. This is in contrast to other studies on the influence of present grassland area which showed an increase in species richness with increasing historical or present patch area (e.g. Dembicz et al. 2020). However, many of the studied grasslands span large, connected, and unregularly shaped areas. This spatial situation is very different from the classical, more or less isolated patches, which are commonly studied (Arrhenius 1921; Ben-Hur and Kadmon 2020). A species-arearelationship, translating into higher plot-level richness on larger grasslands, might exist among the smaller and more compact grassland patches, but remains b Fig. 3 Forest plots of plant subsets (all three regions combined) of which the models with historical and present variables fit significantly better compared to models with only present variables (Vuong test undetected by our models. The habitat amount hypothesis (Fahrig 2013) suggests that species richness depends on habitat amount at landscape scale rather than on local patch size. We observed an impact of proportional grassland area at the landscape scale on some of the studied functional groups (see below: ''Differences between functional groups'' section). In summary, local grassland patch area had a neglectable influence on the studied groups, which is probably mainly due to the arrangement of the patches in the landscape. Historical landscape patch density played an important role for the studied functional groups. On the one hand, complex landscapes provide habitats for different species and opportunity for locally extinct species to be recruited from the surrounding species pool (Tilman 1997;Tscharntke et al. 2002). On the other hand, high patch density means a high degree of fragmentation of habitats, which may hinder colonisation of dispersal-limited species. These opposed effects of patch density are reflected in our results where, e.g., short dispersing plants were negatively affected, whereas long dispersers showed a positive response to high historical patch density. In general, the importance of historical patch density indicates a time lagged response to changes in landscape structure.
We observed landscape evenness, i.e. the distribution of different land cover classes in the landscape, to operate on a short time frame, as mainly present variables affected the species richness of the functional groups. Increasing evenness decreases disturbance spread and enhances recovery rates of species by providing refuge and facilitating a rescue effect, i.e. recolonization after disturbance (Turner and Gardner 2015). We found little variation in historical evenness corroborating this idea. Uniform landscape composition, for example with single, large arable patches, and a reduced grassland area may have short-term impacts as species dependent on heterogenous landscape not only lose their habitat, but also the ability to recover after heavy disturbances (Waldén et al. 2017).
A high nearest neighbour distance among habitat patches of the same class, i.e. larger distances between grassland patches, may affect species richness negatively by decreasing the likelihood of rescue effects and the chance of species to escape from disturbances and competition through affecting the connectivity of patches (Lindborg and Eriksson 2004;Cousins 2006). Yet nearest neighbour distance could also exert a positive effect by preventing the spread of disturbances (O'Neill et al. 1992). We observed almost exclusively positive effects of historical nearest neighbour distance, i.e. the farther away the grassland neighbours were, the higher the observed species richness. The results might be the consequence of disturbances spreading through better connected, closer patches, or due to a few species-rich remnant patches spread across the landscape with species-poor patches in-between, reducing the nearest neighbour distance and thus causing the observed negative effects.
We found no single most important spatial scale for any of the functional groups studied. Historical and present landscape metrics had significant effects at different radii depending on respective groups. Raatikainen et al. (2018) observed 1 km to be a meaningful connectivity threshold for meadow plant species, while investigating different connectivity thresholds in past and present landscape structures of Finnish meadows. Long distance dispersing plant species in our study were only impacted by the largest radius of 2 km, as expected for, e.g., zoochorous species (Helm et al. 2006). Short dispersing species are expected to operate on a shorter radius (Soons et al. 2005), which we could not observe consistently. Over longer time periods, different spatial scales may play a role and thus known effects of present landscape metrics on species richness cannot be applied to past situations. Short distance dispersing arthropod species in our study were mainly affected by the 500 to 1000 m radius, but long-distance dispersers by up to 2000 m. As opposed to plants, arthropod dispersal is usually an active process (informed dispersal, Clobert et al. 2009), where dispersal is dependent on habitat availability, detection ability of suitable habitats, as well as population density and resource availability in b Fig. 4 Forest plots of all sweep net sampled arthropod subsets (all three regions combined) of which the models with historical and present variables fit significantly better compared to models with only present variables (Vuong test (Raduła et al. 2020) and that saproxylic weevil frequency is higher in old woodland (Buse 2012). Colonization of habitat patches depends on both the connectivity of the patch and the conditions of the patch (i.e. land-use intensity in our study). As our grasslands fall along gradients of those confounding factors, this could explain the lack of clear patterns regarding habitat continuity in our study. The semi-natural grassland specialist plants of the phytosociological class Molinio-Arrhenatheretea were negatively affected by grassland permanency, i.e. grassland area changes over time. This is most b Fig. 5 Forest plots of all pitfall trap sampled arthropod subsets (only Hainich region) of which the models with historical and present variables fit significantly better compared to models with only present variables (Vuong test likely an indirect effect, facilitated by high historical land-use intensity on sites with high grassland permanency. However, due to the lack of data on historical land-use intensity, this assumption could not be tested in our study.

Importance of distinct historical time periods
The most important time period, i.e. the period most often included and significant as a coefficient in the GLMs, was 1820/50 (30%), followed by 1960 (23%) and finally 1910/30 (20%). Present variables were almost as important (28%) as our oldest time period.
Other studies have identified historical land use in the past 100 to 150 years to be an important driver shaping the composition of calcareous grassland communities (Heubes et al. 2011), plant community traits (Purschke et al. 2014), and arthropod seed predator communities (Stuhler and Orrock 2016). Our work underlines the importance of historical landscape metrics covering almost 200 years for studies on current species richness, and that a combination of past and present variables best explains species richness, but also that the influence of a distinct time period is communityand landscape metric-specific. However, it is not entirely clear why the landscape metrics from the 1820/50 period appeared to be more important than the other two periods included in our study.

Differences between functional groups
Plants of the phytosociological class of the Molinio-Arrhenatheretea, short and long dispersing plants, herbivorous, low to intermediately dispersing arthropods from the sweep net samples, beetles and spiders, as well as omnivorous, and weakly and strongly dispersing pitfall trapped arthropods were all more strongly affected by historical metrics than the other functional groups studied. In addition, the direction of effects of historical landscape metrics on group richness differed between differently sampled organisms. For example, the effect of proportional grassland area was different for the three species groups: plants, pitfall trapped arthropods and sweep net sampled arthropods. Increased habitat area should increase species richness, but if larger grassland areas were used more intensively in the past or present, as we suggested above, there may be a negative correlation. Pitfall trapped species were exclusively positively affected by historical and present proportional area of grasslands in the surrounding landscape. Potentially because these were positively affected by a higher land-use intensity in the surrounding landscape, which we cannot test with our data. In contrast, plant species richness correlated negatively with the historical grassland area, which we assume to be an indirect effect of higher intensity of grassland use at present and, possibly, already in the past. Sweep net sampled arthropods, except for monophagous species, which should be strongly related to the availability of plants, were not affected by the proportional grassland area. In our modelling approach we always included present land-use intensity as a candidate predictor variable to control for its strong effects on present species richness and, thus, to allow the other coefficients to show their effects independently of current intensity. Although we could not detect collinearity between predictor variables in the final models, we must assume some correlation between land-use intensity and proportional grassland area exists.

Conclusions
Our study showed that both historical and present landscape structure affect the richness of arthropods and plants, but effects differ between functional groups according to their habitat requirements and dispersal ability. The effects of historical landscapes appeared to be long lasting and to remain for up to almost 200 years. Historical landscape components remained important even when accounting for current land-use intensity. In particular, grassland specialist plant species (Molinio-Arrhenatheretea) and arthropod species with low to intermediate dispersal abilities showed strong responses to historical landscape structure. In contrast to the spatial variables, the chronological sequence of local land cover classes seems to have little impact on current species richness. The importance of historical landscapes for current species diversity indicates that historical data should, thus, be considered in conservation planning and management. Nature conservation projects should for example use information on past landscape development and consider its effects on species diversity when planning for suitable landscape structures and to improve management strategies. This might help to prevent an ongoing loss of diversity. Our results further suggest that responses to historical land use strongly depend on species specific traits, which need to be known to identify the optimal landscape structure for a conservation focus taxon. Future research should explore this aspect and the underlying mechanisms in more detail to provide support for halting the current biodiversity loss.
Data availability The datasets generated during and/or analysed during the current study are not publicly available due to an embargo period to protect young researchers but are available from the corresponding author on reasonable request. The datasets will eventually be available via https://www.bexis. uni-jena.de/sws/PublicDataLink.
Code availability The code created by the authors during the current study are available from the corresponding author on reasonable request. A comprehensive list of R packages used can be found in the Online Resource 6.

Declarations
Conflict of interest The authors have no conflicts of interest to declare that are relevant to the content of this article.
Ethical approval Not applicable.
Consent to participate Not applicable.

Consent for publication Not applicable.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.