Contemporary spatial and environmental factors determine vascular plant species richness on highly fragmented meadows in Central Finland

Habitat loss is a major threat to biodiversity. It can create temporal lags in decline of species in relation to destruction of habitat coverage. Plant species specialized in semi-natural grasslands, especially meadows, often express such extinction debt. We studied habitat loss and fragmentation of meadows and examined whether the changes in meadow coverage had caused an extinction debt on vascular plants. We also studied whether historical or present landscape patterns or contemporary environmental factors were more important determinants of species occurrence. We surveyed the plant species assemblages of 12 grazed and 12 mown meadows in Central Finland and detected the meadow coverages from their surroundings on two spatial scales and on three time steps. We modelled the effects of functional connectivity, habitat amount, and isolation on species richness and community composition. We observed drastic and dynamic meadow loss in landscapes surrounding our study sites during the last 150 years. However, we did not find explicit evidence for an extinction debt in meadow plants. The observed species richness correlated with contemporary factors, whereas both contemporary factors and habitat availability during the 1960s affected community composition. Effective conservation management of meadow biodiversity builds on accurate understanding of the relative importance of past and present factors on species assemblages. Both mown and grazed meadows with high species richness need to be managed in the future. The management effort should preferably be targeted to sites located near to each other.


Introduction
Habitat loss and degradation caused by anthropogenic land-use changes are major global threats to biodiversity (Pimm et al. 1995;Sala et al. 2000;MEA 2005;Fischer and Lindenmayer 2007). Concurrently with habitat loss, habitat fragmentation is often implicated in the loss of biodiversity (e.g. Haddad et al. 2015;Wilson et al. 2016). A recent synthesis (Haddad et al. 2015) revealed that decreased fragment area, increased isolation, and the creation of habitat edges consistently degraded ecosystems, reducing species persistence, species richness, nutrient retention, trophic dynamics, and, in more isolated fragments, movement. Such ecological consequences are mediated through three broadly defined mechanisms: loss of habitat amount, changes in the spatial configuration of the landscape, and the different interaction effects of area and configuration and the matrix surrounding the habitat patches (Wilson et al. 2016).
There is vast empirical evidence demonstrating that habitat loss has large and consistently negative effects on biodiversity (Fahrig 2003(Fahrig , 2017Villard and Metzger 2014;Haddad et al. 2015). Several studies have also shown that an increase in the amount of habitat in the local landscape increases species richness (Fahrig 2013). The extent to which the spatial configuration of the habitat patches within the landscape contributes to these outcomes has been a subject to intensive debate (Villard and Metzger 2014). This especially concerns the question whether fragmentation, i.e. the degree to which a given amount of habitat is broken apart, has significant effects on biodiversity atop habitat amount (Fahrig 2003;Villard and Metzger 2014).
On general level, the ecological effects of habitat fragmentation per se are not significant when compared to habitat loss (Fahrig 2003(Fahrig , 2017. Unfortunately, a majority of research thus far does not separate the effects of habitat loss from the configurational effects of fragmentation, i.e. changes in number of patches, patch sizes, and between-patch distances (Fahrig 2003). Yet separation of habitat loss vs. fragmentation effects is important for conservation practices, as they differ from each other. However, those studies that have made this distinction have not been consistent in their conclusions (Fahrig 2003(Fahrig , 2017. Firstly, it is predicted that habitat configuration does not matter and the number of species is simply a function of the total amount of habitat in the local landscape (Fahrig 2013). Secondly, adverse consequences of fragmentation have been demonstrated in cases where there is little habitat left across large areas (Hanski 2011(Hanski , 2015Rybicki and Hanski 2013). Thirdly, habitat fragmentation per se is suggested to have positive ecological effects as it increases the functional connectivity of the landscape by creating larger number of smaller habitat patches with smaller distances between them (Fahrig 2017).
Functional connectivity considers organisms' behavioral responses to individual landscape elements and the spatial configuration of the entire landscape (Kindlmann and Burel 2008). Connectivity is important in determining species distributions as it predicts colonization of habitat patches (Moilanen and Nieminen 2002). For example, when landscape structure is more favorable to movements, a species may require a smaller total area of habitat to persist (Villard and Metzger 2014). For highly fragmented habitats, the best and most consistent performance is found for connectivity measures that take into account the size of the focal patch and the sizes of and distances to all potential source populations (Moilanen and Nieminen 2002). Therefore connectivity as a concept ties together the habitat amount, its spatial configuration, and possible matrix effects, and gives these an ecological meaning.
There is a consensus that higher functional connectivity results in higher probability of species persistence (Villard and Metzger 2014). However, as fragmentation proceeds with habitat loss, functional connectivity deteriorates and eventually the total amount of habitat declines to the point where the species cannot sustain itself (extinction threshold, see Fahrig 2003;Villard and Metzger 2014). Thus, the ecological response of the species to habitat loss and fragmentation is nonlinear (Hanski 2011), and strongest ecological effects of habitat fragmentation should be observable at the intermediate levels of habitat amount (Pardini et al. 2010;Villard and Metzger 2014). When the habitat amount is high, only habitat loss will matter (above so-called fragmentation threshold; see Andrén 1994;Fahrig 2017). If the amount falls below the extinction threshold, the species will not persist (Villard and Metzger 2014). According to metapopulation theory, the extinction threshold refers to the point along a gradient of habitat loss and fragmentation where the metapopulation loses its viability as local colonizations do not suffice to compensate for local extinctions (Hanski 2011).
There is evidence for an interaction between the total amount of habitat and the degree of fragmentation: for more habitat, more fragmentation is possible without additional effect on the number of surviving species (Rybicki and Hanski 2013). Several studies have demonstrated fragmentation effects when the total amount of habitat falls under a threshold of 10-30% of the entire landscape (Andrén 1994;Pardini et al. 2010;Hanski 2011;Rybicki and Hanski 2013). Yet, the fragmentation threshold hypothesis remains to be under continuous testing, as a recent review by Fahrig (2017) found no consistent empirical evidence for it.
In general, the fraction of available habitat that is occupied by a species is an important indicator of population viability (Hanski 2011). In practice, understanding the interactive effects of habitat amount and configuration is challenging because nonlinear and synergistic or antagonistic effects are generally present (Villard and Metzger 2014). Species also respond to landscape changes in different ways depending on their degree of specialization with the focal habitat type (Andrén et al. 1997). Habitat generalists may seem to survive in very small patches, when they actually utilize resources also from the surrounding matrix (Andrén 1994). The situation is further complicated when multiple species are considered simultaneously (Mönkkönen and Reunanen 1999). Species that occur within the same landscape differ in their habitat requirements and functional responses to land cover composition, making analyses of fragmentation effects on multispecies communities a true challenge (Hanski 2015).
The effects of habitat loss and fragmentation on communities and ecosystems can take years to decades before becoming apparent, suggesting that patches will continue to lose species for considerable time periods (Haddad et al. 2015;Wilson et al. 2016). Several studies have demonstrated such delayed loss of species, a phenomenon called extinction debt (Kuussaari et al. 2009;Haddad et al. 2015). If extinction debt exists, the remnant populations are deemed for local extinction after species-specific temporal lags without habitat restoration measures (Kuussaari et al. 2009;Krauss et al. 2010;Haddad et al. 2015). Extinction debt can arise through several mechanisms that are related to species' traits and population dynamics (Hylander and Ehrlén 2013;Haddad et al. 2015). Its existence is a profound challenge to the conservation of species and their habitats (Kuussaari et al. 2009;Hanski 2011).
A common procedure to detect extinction debt involves mapping current species richness and examining to what extent species richness is explained by past and present habitat extent (Kuussaari et al. 2009). Extinction debt exists if excessive number of species, especially habitat specialists, is observed after the loss of habitat: i.e. species richness is higher than what would be expected in relation to the present habitat amount or connectivity (Cousins 2009;Kuussaari et al. 2009;Krauss et al. 2010). For example, species richness of semi-natural grassland plants is often shown to reflect historical rather than contemporary habitat configuration (Lindborg and Eriksson 2004a, b;Helm et al. 2005;Cousins 2009;Krauss et al. 2010;but see Ö ster et al. 2007;Cousins and Vanhoenacker 2011). Although species richness often decreases significantly with increasing habitat loss and fragmentation, other measures of community structure, such as species abundance and community composition, may more accurately inform how these processes affect biotic communities (Lampila et al. 2005;Wilson et al. 2016;Haddad et al. 2015). For example, population size is a more sensitive measure of extinction debt than species richness, as populations decline in species-specific manner before going extinct (Lampila et al. 2005;Hylander and Ehrlén 2013).
In Finland, meadows have faced severe habitat loss and fragmentation due to agricultural modernization during the last 100 years (Salminen and Kekäläinen 2000;Luoto et al. 2003b;Raunio et al. 2008). Meadow habitats depend on constant management actions such as mowing and extensive grazing that create intermediate disturbances, which maintain high biodiversity (Pykälä 2000). During the twentieth century, traditional management practices were abandoned or replaced by intensive agriculture, and meadow cover decreased to less than one percent (Salminen and Kekäläinen 2000;Luoto et al. 2003b). This and suboptimal grazing or mowing of the remaining meadows have led to the endangerment of hundreds of species (Rassi et al. 2010).
The termination of management of meadows has caused significant changes in landscape structure and a decline in the overall number and population sizes of vascular plant species (Luoto et al. 2003b;Rassi et al. 2010). Therefore meadow habitats are a central target for conservation action on national level (Heikkinen 2007). Most of the Finnish studies exploring the ecological effects of past land uses have been conducted in southern parts of the country (Luoto et al. 2003a;Pykälä et al. 2005;Pykälä 2007;Pitkänen et al. 2014Pitkänen et al. , 2016Lampinen et al. 2017). However, due to regional differences in biogeography and agricultural history (Ahti et al. 1968;Soininen 1974), the biological legacies may express alternative patterns in other parts of the country.
In this paper, we study whether meadow loss and fragmentation from the mid-nineteenth century onwards has caused extinction debt on contemporary meadow species in Central Finland. We compare the effects of past and present habitat amount, connectivity, and isolation together with contemporary environmental factors on vascular plant richness and community composition. We hypothesize that the meadow coverage has become strongly fragmented from the mid-nineteenth century to the present day, and thus remnant species communities may witness extinction debt, expressed as a causal relationship between species and historical amount and/or connectivity of meadow habitat. Our results give insights into what criteria (e.g. past or present spatial factors, or current environmental site attributes) are most influential in optimizing future conservation actions concerning these threatened habitats. The results also indicate if the total amount of meadow habitat, historical or contemporary, is more important than habitat configuration in explaining the current meadow species assemblages.

Study sites
We studied 24 meadow sites located in the province of Central Finland (Fig. 1). All study sites were located in the southern half of the province in order to reduce biogeographical turn-over in their species assemblages. We included moist, mesic, or dry meadows under either mowing or grazing management (n = 12 and 12, respectively). The latter were grazed by horses, cattle, and/or sheep. All meadows had been managed yearly through mowing or grazing for at least 10 years, and their minimum size was 0.1 ha.

Plant and soil data
We conducted vegetation surveys on the study sites between late June and early July 2014. On each site, we established a 44 m edge-to-center study transect with a randomized compass course determining the Fig. 1 The locations of the study sites (IDs 1-24) within biogeographical vegetation zones in southern Central Finland. Mown sites (1-12, square symbol) and grazed sites (13-24, circle symbol) are differentiated on the map. Note that on three occasions (24/9; 10/15; 11/23) symbols of grazed and mown sites are overlapping on the map since the sites are located on the same farms. Index map (below) shows the overall location of the study area in relation to the cities of Helsinki and Jyväskylä in Finland starting point. On small or narrow meadows with limited space for transects, we split transects into two parts so that the first part covered the original transect and the remaining length was positioned perpendicularly to the first, crossing it at its mid-point. In every ten meters, we located a 2 m 9 2 m plot to the right side of the transect, starting at 2 m point and ending at 42 m. We sampled five plots on each meadow, comprising together 20 m 2 . Within each plot all vascular plant species and their abundances (as percentage cover) were recorded. The plot-wise species occurrences were pooled in order to determine site-level species occurrences. Since we were interested in the responses of species that typically grow in meadows, we extracted meadow species from the list of all species based on the list provided in Appendix 1 in Pykälä (2001). In this list, Pykälä has included all vascular plant species that have been characteristic for Finnish meadows in the late nineteenth century. In the subsequent analyses we included only the number of these meadow species as a response variable and excluded all other species.
Average vegetation height (cm) was assessed in the study plots. This was repeated three times during the growing season in 2014: In June-July and in September-October. For each site, the average of recorded values was used in the analyses.
We collected soil samples from the plots in May 2014. We used a soil corer of 3 cm diameter to take two samples from the opposite corners of each plot, totaling 10 samples per site. Samples were taken from soil surface to 5 cm depth, excluding undecayed litter. The samples taken from each site were mixed, homogenized with a 4 mm mesh sieve, and preserved in the freezer. Soil pH was measured from a suspension of 6 ml soil and 30 ml 0.01 M CaCl 2 . We used a median of three pH measurements in our analyses.

Landscape variables
We created 1 and 2 km radius buffer zones (hereafter study landscapes) around each study site for historical and contemporary meadow coverage detection. We used two distance thresholds (the radii), since the exact scale of ecological effect was uncertain (Fahrig 2013). These selected distances correspond to reasonable dispersal distances for plant species living on meadow habitats (Lindborg and Eriksson 2004a;Helm et al. 2005;Cousins and Lindborg 2008). The average dispersal ranges from less than 1 to 3 km (Helm et al. 2005;Soons et al. 2005).
To map meadow coverages, we acquired historical maps and contemporary geospatial data covering the study landscapes from three time periods: 1841-1855 (time step 1), 1961-1972 (time step 2) and year 2010 (time step 3). For time step 1 we acquired cadastral maps as photographed map sheets from the National Archives of Finland. Time step 2 consisted of scanned topographical map sheets from Central Finland Centre for Economic Development, Transport and the Environment. Geospatial data for time step 3 was obtained from the vector-format topographical database of the National Land Survey of Finland. Topographical and cadastral maps were georeferenced and rectified in ArcGIS to match spatially with the 2010 topographic database (using ESRI Ò ArcMAP TM version 10.3).
All meadow patches within the study landscapes were interpreted and digitized from the maps (Fig. 2). The digital conversion process into geospatial meadow layers was done retrospectively beginning with the latest time step 3 (Bender et al. 2005;Mökkönen 2006). As the cadastral maps had gaps and missing land cover information in some places, we classified both meadows and missing data areas for the oldest time step 1. In topographical maps (time step 2) such problems were not encountered. The interpretation and digitizing was done in scale of 1:6000 with a minimum mapping unit of 10 9 10 m, determined from the scale of the original maps.
Past and present habitat connectivity in each study landscape were calculated from the meadow coverage data according to the time layers. We used a graphtheoretic approach to functional connectivity that is based on a patch-corridor-matrix model and measures habitat availability (Pascual-Hortal and Saura 2006; Laita et al. 2011). The chosen measure, integral index of connectivity (IIC), is shown to provide reliable information on habitat amount and the degree of connectivity between patches (Pascual-Hortal and Saura 2006). It recognizes increasing topological distances between patches as lower connectivity, and favors habitat located in a single large patch (Laita et al. 2011). Because it treats connectedness of patches as binary, it is suitable for exploring scale-dependency in connectivity effects (i.e. which is the threshold distance for observed ecological response).
In ArcGIS, rasterized meadow coverages for each of the three time steps were extracted and re-vectorized. Patch sizes and all edge-to-edge distances were calculated for every meadow patch within the study landscapes. These data were analyzed using Conefor 2.6 (an updated version of Conefor Sensinode 2.2; Saura and Torné 2009), which is a software designed for graph-theoretic connectivity analyses. We calculated IIC within each study landscape in every time step separately, given by: where a i is the area of each patch and nl ij is the number of links in the shortest path (measured as Euclidean distance) between patches i and j. A L is the total analyzed landscape area, comprising both habitat and non-habitat patches. In forming links, respective to the radius of the studied landscape, we used a distance threshold of either 1 or 2 km. Patches that were less apart from each other than the threshold distance were considered as connected to each other. In calculating time step 1 connectivity, we adjusted the A L value to correspond to the area where data was available. In more recent maps A L was equaled to the total area of the study landscape. This opportunity to adjust the value of the total landscape area was an important reason why we chose to use IIC.
For each study site, each time step, and both landscape radii, we calculated two functional connectivity index values: a regular IIC with information on patch sizes and a distance-based IIC with a fixed patch  Fig. 2 Example of the geospatial data from one of the study sites (ID 22; located in the center of the circle). Top maps show the original map data sets from the three time steps, and below are their corresponding meadow patches. Buffers with 1 km (dashed line) and 2 km (solid line) radii exemplify how the data was delineated into study landscapes for the analyses. For a full color version of the figure, the reader is referred to the online version of the article size attribute (a i = a j = 1 ha). This allowed us to control for the effect of inclusion of patch size information into the index, giving an index value that was based on the variation of inter-patch distances and number of patches within the threshold distance.
To compare the functional connectivity indices to the total amount of habitat, the areas of all meadow patches were pooled within the study landscapes. Such buffer measures are commonly used in connectivity studies and they measure structural rather than functional connectivity (Moilanen and Nieminen 2002). For time step 1, the amount of meadow area was corrected according to the missing data (small destroyed patches or white areas in the maps). We assumed that the amount of meadows that we had no data on was proportional to the recorded meadow coverage within each study landscape. For each landscape, we divided the combined amount of habitat with the A L value and added a respective proportional area of missing meadows to the sum of habitat amount. In further analyses we square-root transformed the total amount of meadows to linearize the cumulative species-area relationship (MacArthur and Wilson 1967).
In addition, we analyzed buffer-independent patch isolation for each study site and each time step. This was measured as the nearest-neighbor distance to the closest meadow patch edge from the starting point of the study transect. Distance to nearest-neighboring population is the simplest form of connectivity measures and it is widely applied both in spatial ecology and in reserve site selection literature in conservation biology (Moilanen and Nieminen 2002). However, it is a poor predictor of species viability, and therefore it is not advisable to rely solely on it when measuring connectivity (Moilanen and Nieminen 2002).
Finally, to derive coarse information on the landuse history of the study sites, we calculated a simple age depth variable. We categorized the binary occurrence of meadow cover on each site for every time step and added up the values. Thus the variable included information on the number of time steps when the site had been meadow.

Statistical analyses
The richness of meadow species was analyzed using generalized linear models (GLMs) with Poisson distribution. We had seven landscape variables, with data from the three time steps for each of the variables, and given the small sample size (n = 24) it was not reasonable to fit all of the variables in one global model. Furthermore, some of the variables correlated strongly with each other. By separating the landscape variables into different models we were able to compare their relative support over each other. Thus, in each model, the effect of past vs. present state of one landscape variable on meadow species richness was analyzed. As a result, we built seven GLMs: 1. regular connectivity (IIC) of meadows within 2 km radius; 2. distance-based connectivity (IIC distance) within 2 km radius; 3. the amount of meadow habitat within 2 km radius, square-root transformed; 4. regular connectivity (IIC) of meadows within 1 km radius; 5. distance-based connectivity (IIC distance) within 1 km radius; 6. the amount of meadow habitat within 1 km radius, square-root transformed; and 7. the distance to the nearest neighboring meadow.
All of the landscape variables were measured on all three time layers. In every GLM, we also included other environmental variables: Management type (grazed or mown), soil pH, the quadratic effect of soil pH, the interaction between management type and soil pH, the age depth of the study meadow, and the East and North coordinates of the site. Vegetation height was not included in the models because it correlated strongly with soil pH. Before analyses, all continuous explanatory variables were standardized into zero mean and unit variance. The lack of spatial autocorrelation in the full model residuals was confirmed with Moran's I test (using 2-and 4-nearest-neighbor structures).
We subsequently compared the likelihoods of the separate models given the data on meadow species richness. The models were ranked based on their second-order Akaike Information Criterion (AICc) values, and those having empirical support in terms of AICc differences were simplified and parameterized through model averaging (Burnham and Anderson 2002;Grueber et al. 2011). For every full model, all possible submodels were derived, including also an intercept-only model. From this set of models, those fulfilling a cut-off criterion of top 2 AICc were extracted and averaged in order to estimate the final model parameters, using the so-called zero method (Burnham and Anderson 2002).
The site-wise abundances of meadow species were analyzed with canonical correspondence analysis (CCA). The constrained ordination method enabled the use of environmental factors in explaining composition of species communities, thus complementing the GLMs that were used to model meadow species richness. In CCA, we built seven models that were similar to GLMs, except that the quadratic effect of pH and the interaction between management and pH were excluded. Model-building in CCA was based on permutation tests (with significant terms added sequentially using 999 permutations). As the information-theoretic approach is not yet compatible with ordination analyses (Burnham and Anderson 2002), we utilized statistical approach based on null-hypothesis testing in the community analysis.

Meadow loss and fragmentation
We observed declines in the functional connectivity (Fig. 3a, b and d, e) and the total amount of meadow habitat (Fig. 3c, f; Table 1) before year 2010 (time step 3). In the mid-nineteenth century (time step 1), the proportional meadow coverages surrounding our study sites ranged from zero to 17.2% (within study landscapes defined by the 1 km radius). The current meadow coverage was 2.3% in maximum, and in twenty cases less than one percent (for mean values, see Table 1). A median of 7.7% of the meadow coverage on time step 1 remained on time step 3. Even though the general trend during the time interval was a decrease in coverage, the magnitude of the change ranged widely and we noticed an increase in meadow coverage on three specific occasions.
The decline in the amount and connectivity of meadows were similar on both studied spatial levels (study landscapes with 1 and 2 km radii). In unison, the distances from study sites to their nearest-neighboring meadow patch increased (Fig. 3g). Thus the study sites have become increasingly isolated as the meadow coverage decreased. This general observation, however, was not seen in the mean edge-to-edge distances that were measured between all meadow patches located within the 1 km radius. On those occasions where there was more than one meadow patch within the study landscape, mean between-patch distances halved during the studied time period (Table 1). This indicates that the few remnant meadow patches are on average near to each other. Yet in six study landscapes the contemporary map data included only one meadow patch, and one even zero, because one of the study sites was mislabeled as forest.
Comparison of the three time steps indicates that the loss of meadow habitat has been pronounced before the 1960s (time step 2; Fig. 3c, f, Table 1). This initial decline in habitat amount did not increase nearest-neighbor distances; the isolation effect became observable only after the middle time step (Fig. 3g). Similarly, the number of meadow patches did not decline linearly, but a drastic drop was observed before year 2010 (time step 3; Table 1).
On some occasions habitat loss was tied to fragmentation, which is reflected in peaks of distance-based IIC values on time step 2. Within these study landscapes, the increase in distance-based connectivity between time steps 1 and 2 indicates more patches with shorter inter-patch distances (Fig. 3b, e). For example, distance-based IIC calculated for the landscape surrounding study site 22 on time step 2 is the highest in Fig. 3b and second highest in Fig. 3e (the data are depicted in Fig. 2). In general, the meadow patches have been small throughout the study period, but they were at their smallest on time step 2 (Table 1). This, also, is in accordance with the fragmentation effect (more and smaller patches nearer to each other). However, fragmentation per se was not always observed: within several landscapes habitat loss occurred together with declining number of patches and increasing between-patch distances (Fig. 3b, e).

Vascular plant species richness
A total of 156 vascular plant species occurred on the study sites, with a mean species richness of 43.6 species per site. Seventy-four of the species were classified as meadow species and their mean richness was 30.3 species per site. No threatened species were observed, but two near-threatened species were recorded: Botrychium lunaria and Dianthus deltoides. All observed species are listed in Online Resource 1.
Comparison of the full models indicated that the amount of meadow habitat within 1 km radius best explained the number of meadow species when compared to the other landscape variables (Table 2). Three models that were very close to each other in  terms of model selection statistics followed it: regular connectivity (IIC) within 2 km radius, nearest-neighbor distance, and distance-based connectivity within 1 km radius. The differences in the AICc values of these models fell within four units from the best model, giving a total sum of Akaike weights of 0.87. Other models had considerably less empirical support (with AICc differences [ 4; Burnham and Anderson 2002) and were omitted from model averaging.
In the simplified models, variables that explained meadow species richness were mostly contemporary (time step 3) factors (Table 3, Fig. 4). Of these, the most important variables were the current habitat amount and distance-based connectivity, both being present in all selected models (Table 3, Fig. 4d, e). Both were measured within 1 km radius and associated with increasing species richness. When the variables were compared with each other, habitat amount at 1 km radius had the highest predictive power with the data (Table 2). Two historical variables were left in the models: regular connectivity (2 km radius IIC) at time step 2 (Fig. 4c) and habitat amount at time step 1 (within 1 km radius; models 1 and 6 in Tables 2 and 3). Their effects, however, were weak. The distance to the nearest neighboring meadow on time step 3 had a slight effect on species richness: When the distance increased, species richness decreased (Table 3, Fig. 4f). Also other contemporary factors, excluding management and management*pH interaction, were left in the simplified model 7. Among these, the most important ones were the observed negative effect of the East coordinate (Fig. 4a) and the quadratic effect of soil pH (Fig. 4b). These variables were also retained in most of the other models (Table 3).

Community composition of meadow species
In Canonical Correspondence Analysis, the models that had significant impact on the community composition of meadow species included regular connectivity within 2 km radius and the combined amount of meadow habitat within both 1 and 2 km radii (Table 4). During simplifying the models, contemporary landscape variables were removed and only time step 2 values were retained including regular connectivity and habitat amount within 2 km radius (Table 5). In addition to these remnant effects from 1960s, the simplified models showed that the community composition of meadow plants was associated with variation in soil pH and management type (Table 5, Online Resource 2). The best fit model, with 28.8% of the total inertia explained, included time step 2 connectivity (2 km regular IIC), management type, and soil pH (Table 5). AICc gives second-order AIC that is corrected for small sample size D i are AICc differences relative to smallest AICc among the models Log-likelihood of each model is given in log(L) K is the number of estimable parameters in the model The rank order of the models from best fit to least fit is based on Akaike weights (w i ) and AICc differences D 2 gives the ratio between deviance explained by the model and null deviance

Discussion
Strong habitat loss but weak evidence for extinction debt We observed drastic changes in meadow coverage in landscapes surrounding our study sites during the last 150 years. As we hypothesized, meadows have faced severe habitat loss, but this overall process seems to have been locally highly dynamic. We observed that on different locations, the rate and magnitude of change varied; in few study landscapes even the direction was opposite, resulting in increased instead of decreased habitat coverage. This indicates that the temporal continuity of meadow patches is likely to be low in Central Finland. Given the magnitude of overall meadow loss since the mid-nineteenth century, it was surprising that we did not observe explicit historical effects on presentday vascular plant species richness. Observed numbers of meadow species resulted from contemporary factors rather than from historical habitat amount and/ or configuration, suggesting that extinction likely plays a minor role. Based on historical records, populations of several meadow species have declined in Central Finland during the late twentieth and the early twenty-first century, and this is mainly due to the cessation of meadow management (Uusitalo 2007). The overgrowing of abandoned meadows is the second most important reason behind both regional extinctions and endangerment of species in Finland (Rassi et al. 2010).
In our data, the greatest decline in habitat coverage occurred after the mid-nineteenth century, but we did not observe that this change had effects on species richness or communities. This conclusion may be partially due to restrictions in our analytical approach. Unfortunately, there were no historical species occurrence data available from our study sites, and therefore we could not analyze temporal trends in species richness. Instead, we used a common but perhaps simplistic method that seeks for correlations between current species richness and past landscape structure (Lindborg and Eriksson 2004b;Kuussaari et al. 2009). However, the use of species abundances in the community analysis seemed to complement this approach as it revealed that current community composition was affected by the amount and connectivity of meadow habitat during the 1960s. Certain species, whose occurrence was connected to the For soil pH, both linear and quadratic (x 2 ) effects were included in the models

Management*pH refers to an interaction between the two variables
In model 6, habitat amount was square-root transformed on all time steps (1: 1841-1855, 2: 1961-1972, and 3: 2010) historical landscape structure, are considered indicators of high meadow quality (Bistorta vivipara, Pilosella lactucella, and Trifolium medium; see Online Resource 2 and Pykälä 2001). This suggests that there are certain meadow species that express temporal lags in their decline, but this is not reflected to the overall species richness. Thus, the extinction debt as described by Kuussaari et al. (2009) was not evident in our data. In a recent study, Lampinen et al. (2017) surveyed 74 grassland generalist and 41 specialist species from power line clearings in Southern Finland and noted that both groups responded positively to late nineteenth century connectivity of nearby meadows. Using the same species categorization, we observed 74 generalist but only four specialist species (of which one was B. vivipara), and a strong relation to contemporary, but not to historical habitat pattern. Lack of specialists in our data may be partially caused by survey sampling; an inventory of total species richness of the sites would have resulted in more exhaustive data (Ö ster et al. 2007). On the other hand, we were interested in capturing the average number and abundance of species on the study sites. As the study meadows differed in size, sampling was needed to produce comparable data sets. As we observed strict meadow specialists only rarely, we suggest that the populations of meadow species may either already have become extinct or faced significant decline on our study sites, or that the extinction debt never developed on the sites.
The first alternative conclusion is supported by a comparison to other empirical studies. It is common to argue that extinction debt is already paid, if the contemporary habitat distribution is sparse and species richness follows its pattern (Cousins 2009;Kuussaari et al. 2009;Kolk and Naaf 2015). A Swedish metastudy demonstrated that extinction debt was evident in landscapes where over 10% of the original grassland habitat was still left; in more fragmented landscapes species richness was related to contemporary habitat pattern (Cousins 2009). Furthermore, in landscapes with small and isolated patches extinction debt will be paid off faster (Kuussaari et al. 2009). Of our study landscapes, 58.3% were below the ten percent threshold, when the loss of habitat between the oldest and the most recent time steps was calculated within 1 km radius. Generally low connectivity values reflected sparse and isolated habitat distribution. Low connectivity and habitat amount was observed throughout the temporal extent of our study, leading to a suspicion on whether meadow coverage has ever been extensive in Central Finland. If this is the case, it is well possible that the meadows we studied have never been colonized by some of the specialist species and thus there has not been extinction debt to begin with. This second alternative conclusion has grounds on the theory underlying the concept of extinction debt. Extinction debt is inherently based on an assumption of an equilibrium state where the number of species is not changing because the rate of local extinctions equals the rate of local colonizations (Kuussaari et al. 2009). In case of long-lived vascular plants, development of such equilibrium requires temporal continuity from the habitat patches, and local habitat history can modulate the existence of extinction debt (Koyanagi et al. 2017). Our observation on the landscape-specific temporal changes in the meadow coverage patterns indicates that the process of meadow loss is not straightforward in Central Finland. It rather follows local land-use dynamism; this spatiotemporal instability of habitat patches can prevent the species richness from reaching the equilibrium state.
This latter alternative is, on our opinion, the more likely explanation of the lack of evidence for extinction debt on our study sites. Based on the observation that habitat availability during the 1960s was reflected in meadow species' community structure, but not on site-level species richness, we suggest that determining extinction debt in terms of species' population declines would be a more sophisticated method. Ultimately, extinction debts arise from delayed extinctions of many different species, each of which will respond individualistically to the observed destructive change in habitat (Hylander and Ehrlén 2013). Furthermore, the dynamic local context calls for new methodologies in measuring such changes within landscapes. Treating present and historical patterns of habitat cover as distinct and separate categories ignores the level of shared history of habitat patches (Ewers et al. 2013). Even within the same landscape, and considering single species, extinction debt can vary in response to local habitat history (Koyanagi et al. 2017).
Local contemporary factors impact meadow species richness and community composition Based on our results, contemporary habitat amount is an important but scale-dependent driver for richness of meadow plants in Central Finland. This supports earlier findings that present-day grassland proportion positively influences species richness at the local scale (e.g. Cousins and Vanhoenacker 2011). We found that increasing total habitat amount increased species richness in landscapes with 1 km, but not with 2 km radius. At the same time, nearest-neighbor distances between meadow patches decreased. The connection between nearest-neighbor distance and habitat amount is based on the fact that both measure the availability of habitat in the immediate landscape surrounding the focal patch (Fahrig 2003). They are, however, nonlinearly connected with each other (Andrén 1994). This may explain why the response of species richness was observed only on the smaller spatial level: the 1 km distance may be a meaningful connectivity threshold for meadow species in Central Finland. One kilometer distance is often used in studies concerning forest plants (Kolk and Naaf 2015). In wind-dispersing grassland plants seed dispersal distances are unlikely to be more than 0.1 km (Soons et al. 2005), but for zoochorous species that are dispersed by grazers average distances of 3 km are realistic (Helm et al. 2005). In addition to habitat amount, current distancebased connectivity explained meadow species richness. This underlines the importance of the number and reachability of available habitat patches for meadow plants. We observed that from the midnineteenth century onwards, the overall number of meadow habitat patches declined, the patches decreased in size, and their edge-to-edge distances shortened. Finally, due to the pronounced habitat loss and fragmentation, the remnant meadows have become isolated and this may be reflected in the realized dispersal distances. Whereas the contemporary habitat amount and configuration effects were observed on the 1 km radius, the historical effects were more prevalent on the 2 km radius. This may suggest that the dispersal distance of meadow species has become shorter. This could be possible if the number of dispersal propagules has decreased along with decreasing population sizes.
In addition to habitat availability, we found soil pH and East coordinate to be important in determining vascular plant species richness on our study sites. More westerly located sites had more species than easterly located sites. We suspect that this results from a well-known climatic west-to-east gradient across Finland that affects vegetation and species richness (Ahti et al. 1968).
The importance of soil pH for plant diversity in Europe is well known, but usually the relationship between the two is positive (Dupré and Ehrlén 2002;Pärtel et al. 2004). Interestingly, our finding partially contrasts this, as we observed higher species numbers at intermediate soil pH values. An earlier study conducted on wood-pastures located in the same geographical region found that vascular plant richness increased with increasing pH, but also in this case the effect was humped, with a peak around pH of 4.5 (Oldén et al. 2016). In this wood-pasture data the pH values ranged from 3.1 to 4.9 and in our meadow data the range was from 3.7 to 5.2. This indicates that in boreal semi-natural habitats the relationship between vascular plant species richness and soil pH may be humped, and peaks somewhere between values 4.0 and 4.5, but testing this assumption needs further effort.
Another mechanism that may explain both the nonlinear effect of soil pH and the previously mentioned observed lack of meadow specialists is habitat deterioration due to low management quality. We found that pH was positively correlated with average vegetation height, and it affected species communities. For example, two tall-growing species associated with high soil pH were Alopecurus pratensis and Anthriscus sylvestris, both of which are indicators of poor management (Online resource 2, Pykälä 2001). Tall vegetation is generally considered as a sign of inadequate management on dry and mesic meadows (Raatikainen 2009). Ineffective removal of biomass may lead to proliferation of rapidly growing strong competitors, which reduces the overall species richness (Grime 1973). According to our field observations, some of the study sites suffered from low-quality management. However, as we did not directly analyze the effect of management intensity on species richness, the question on the possible effects of habitat deterioration remains open and needs further exploration.
Although management regime did not affect species richness, mown and grazed meadows differed from each other in terms of their communities. A recent meta-analysis on the biodiversity benefits of grassland management revealed only small differences in the effect of grazing and annual mowing, but indicated that grazing had an overall more positive effect (Tälle et al. 2016). Our results complement this conclusion by highlighting that mown and grazed meadows host different species. We agree with Tälle et al. (2016) that the mechanisms involved in the effects of grazing and mowing are not clear and need further clarification. The studies would benefit of including also other environmental and historical factors, as we found out that community structure was also significantly affected by soil pH and past habitat availability.

Conclusion
Based on our results, we suggest that in order to effectively conserve meadow biodiversity in landscapes where meadows are already highly fragmented, a focus on contemporary factors such as location, soil acidity, and management method is of essence. The existence of possible biological legacies depends on local to regional land-use history, and because of landscape-specific differences it is advisable to confirm whether extinction debt is present or not in each case. We found that remnant meadows' species communities respond to habitat amount and connectivity 50 years ago, but this effect did not express extinction debt, which is usually examined in terms of species richness.
We recommend usage of multiple levels of analysis as species' responses to habitat amount and configuration are scale-dependent and individual. In Central Finland, both mown and grazed meadows with high species richness need to be managed in the future. The management effort should preferably be targeted to sites located near to each other. This would maximize the richness of species benefiting from the management actions and reduce the isolation effect that most likely delimits the dispersal of meadow species within the current fragmented landscape.
We suggest that future studies concerning highly fragmented meadows could explore which other key environmental factors, in addition to soil pH and management regime, affect their biodiversity. Considering spatiotemporal factors, biodiversity legacies in terms of possible connections between patchspecific land-use history (so-called terragenies) and species' populations' persistence would be an interesting research subject.