Drivers of plant diversity in Bulgarian dry grasslands vary across spatial scales and functional-taxonomic groups

Questions: Studying dry grasslands in a previously unexplored region, we asked: (a) which environmental factors drive the diversity patterns in vegetation; (b) are taxonomic groups (vascular plants, bryophytes, lichens) and functional vascular plant groups differently affected; and (c) how is fine-grain beta diversity affected by environmental drivers? Location: Northwestern and Central Bulgaria. Methods: We sampled environmental data and vascular plant, terricolous bryophyte and lichen species in 97 10-m 2 plots and 15 nested-plot series with seven grain sizes (0.0001–100 m 2 ) of ten grassland sites within the two regions. We used species richness as measure of alpha-diversity and the z -value of the power-law species–area relationship as measure of beta-diversity. We analysed effects of landscape, topo-graphic, soil and land-use variables on the species richness of the different taxonomic and functional groups. We applied generalised linear models (GLMs) or, in the presence of spatial autocorrelation, generalised linear mixed-effect models (GLMMs) in a multi-model inference framework. Results: The main factors affecting total and vascular plant species richness in 10-m 2 plots were soil pH (unimodal) and inclination (negative). Species richness of bryophytes was positively affected by rock cover, sand proportion and negatively by inclination. Inclination and litter cover were also negative predictors of lichen species richness. Elevation negatively affected phanerophyte and therophyte richness, but positively that of cryptophytes. A major part of unexplained variance in species richness was associated with the grassland site. The z -values for total richness showed a positive relationship with elevation and inclination.


| INTRODUC TI ON
Global change, including climate and land-use change, is one of the major threats for biodiversity. However, the magnitude and direction of such changes varies among regions, ecosystems and taxonomic groups (Rounsevell et al., 2018). High plant species richness, i.e. the number of all plant species in a sample plot, is a biodiversity indicator linked to well-functioning ecosystems (Soliveres et al., 2016). Globally, semi-natural grasslands harbour the highest species richness at small spatial scales (Wilson et al., 2012). However, over the past 100 years in the Palaearctic region, species-rich semi-natural grasslands underwent a strong decline in area and diversity because of habitat conversion to arable and urban land and degradation due to land-use changes such as intensification and abandonment (Rounsevell et al., 2018;Dengler et al., 2020b). Land-use intensification is often associated with plant diversity loss due to the high loads of fertilisers and overgrazing.
Similarly, land-use abandonment negatively affects plant diversity and species composition (especially of typical grassland species) because successional grasses, shrubs and trees encroach the abandoned grasslands (Rounsevell et al., 2018;Dengler et al., 2020b).
The relevance of factors explaining diversity patterns varies greatly among studies. While anthropogenic factors were the most important in some studies (Maurer et al., 2006), in others abiotic factors had a greater effect (Merunková et al., 2014). Among the latter, high habitat heterogeneity (Stein et al., 2014;Boch et al., 2016), soil parameters (Moeslund et al., 2013), elevation (McCain & Grytnes, 2010) and climate (Chytrý et al., 2007;Palpurina et al., 2015) were found to shape diversity in grasslands. However, the direction of the effect and the importance of these factors can vary among studies and regions (Klaus et al., 2013;Palpurina et al., 2015). For instance, the relationship between soil pH and plant species richness in dry grasslands was reported as positive (Kuzemko et al., 2016), negative (Chytrý et al., 2007), unimodal (Mathar et al., 2016), not significant (Cachovanová et al., 2012), or varying, often because of the confounding effects of other factors such as climate (Palpurina et al., 2015(Palpurina et al., , 2017 or scale-dependency of species diversity (Turtureanu et al., 2014).
In general, drivers of species diversity and their relative importance vary across spatial scales (Shmida & Wilson, 1985;Siefert et al., 2012). For larger grain sizes, Field et al. (2009) demonstrated this phenomenon with a comprehensive meta-analysis. However, evidence for smaller spatial grains, where direct interaction between species can be assumed, is still sparse. Siefert et al. (2012) conducted a meta-analysis with vegetation data, for grain sizes from 0.01 m 2 to more than 1 km 2 , showing that the relative importance of edaphic vs. climatic variables decreases with increasing grain size. For plot-scale heterogeneity, Lundholm (2009) found no grain-size effect on plant diversity, while Tamme et al. (2010) showed that the positive effect of heterogeneity on species richness turns negative towards very small grain sizes. As all these studies were meta-analyses, there is a lot of variability in the data due to other factors and the interpretation therefore is difficult. This limitation could be overcome by studying scale effects using different grain sizes in the same system (Turtureanu et al., 2014;Kuzemko et al., 2016;Polyakova et al., 2016).
Moreover, different taxonomic groups like vascular plants, bryophytes and lichens can be differently affected by abiotic and anthropogenic factors (Löbel et al., 2006;Polyakova et al., 2016). However, this has only rarely been explored. In addition, dividing the overall vascular plant diversity into functional groups (e.g. grasses, legumes, herbs) or groups related to life-history traits (e.g. life forms) can give additional information, as these groups might respond differently (Chytrý et al., 2010;Palpurina et al., 2015).
Environmental factors also shape beta-diversity (Keil et al., 2012). For Palaearctic grasslands it has been shown with a comprehensive data set of nested-plot series that fine-grain species-area relationships (SARs) are adequately modelled with the power function (Dengler et al., 2020a). The slope parameter z of this function is an appropriate measure of multiplicative beta-diversity (Koleff et al., 2003;Jurasinski et al., 2009;Polyakova et al., 2016) and has been widely used to compare beta-diversity between different ecological situations (Drakare et al., 2006;de Bello et al., 2007).
Despite their high diversity with many unique vegetation types (Tzonev et al., 2009), Bulgarian grasslands have only recently received more attention (Velev & Vassilev, 2014;Palpurina et al., 2015). However, neither diversity patterns of different taxonomic and functional vascular plant groups (but see Palpurina et al., 2015), nor diversity patterns across different spatial scales (including beta-diversity) have previously been studied in Bulgarian grasslands.
We thus investigated fine-grain species richness patterns of vascular plants and non-vascular species in dry grasslands of two distinct mountain regions in Bulgaria. We asked: (a) how species rich are Bulgarian grasslands at different spatial grains and which environmental factors explain the patterns best; (b) are taxonomic groups (vascular plants, bryophytes, lichens) and functional groups of vascular plants differently affected; and (c) how is fine-grain beta-diversity of these grasslands affected by environmental drivers?

| Field sampling and lab measurements
The field sampling was conducted 14-24 August 2011 during the 3rd Research Expedition of the EDGG (www.edgg.org). We sampled four grassland sites in the Vratsa and six in the Koprivshtitsa region. Plots within sites were selected ensuring low within and high betweenplot variability in species composition (Ewald, 2003a). To address the scale dependence of biodiversity, we sampled 15 so-called "biodiversity plots", which are square plots of 100 m 2 , with two nested subplot series of 0.0001-, 0.001-, 0.01-, 0.1-, 1-and 10-m 2 plots in two opposite corners of the plot (Dengler et al., 2016b). To cover a wider range of grassland diversity, we sampled 67 additional 10-m 2 plots ( Figure 1).
In each plot and subplot, we recorded the shoot presence of all terricolous vascular plant, bryophyte, lichen and macroscopically visible "macroalgae" species (Dengler, 2008). Throughout the text, we use total species richness to refer to the combined richness of all taxonomic groups and non-vascular species to refer to bryophytes, lichens and macroalgae. Vascular plant species were further grouped by life forms and functional groups. Life forms follow Raunkiaer (1934): phanerophytes (trees and shrubs), chamaephytes (dwarf shrubs), hemicryptophytes (herbaceous perennial plants with buds near the soil surface), cryptophytes (perennial plants with subterranean buds) and therophytes (annuals). Functional groups were based on their taxonomic position: graminoids = Poaceae, Cyperaceae and Juncaceae; legumes = Fabaceae; other = all other families (Appendix S1).
Within each 10-m 2 plot, we measured coordinates and elevation with a handheld GPS. We also measured slope inclination, aspect and maximum microrelief (Dengler et al., 2016b) and estimated the percentage cover of stones and rocks (henceforth "rock cover"), fine soil and litter. Grazing intensity was estimated as four quasi-metric levels (0, not grazed, 1, low intensity, 2, medium intensity, 3, high intensity) based on dropping density and trampling disturbances.
For soil properties, we took a mixed soil sample from the uppermost 10 cm at five random points within each plot and measured pH and electrical conductivity with a standard glass electrode in a suspension of 10 g air-dried soil in 25 ml distilled water. Organic matter content was measured by loss on ignition at 550°C for 16 hr. Proportions of the soil particle fractions sand, silt and clay were measured using a Particle Size Distribution Analyzer (HORIBA LA-950V2, Kyoto, Japan). Due to the lack of clay fraction in most of the samples (76 plots out of 97) and very small clay contents (1-3%) in the rest of the samples, we did not use clay content for further analyses.
As a proxy of drought stress, we calculated the heat-load index using aspect and slope, assuming that slopes with 225° aspect receive the highest diurnal heat load in the northern hemisphere (Parker, 1988). Thus, the steepest SW-exposed slopes have maximum (positive) and the steepest NE-exposed slopes minimum (negative) values while flat areas have zero values. We also retrieved five climatic variables from WorldClim 2 at a spatial resolution of ~1 km (Fick & Hijmans, 2017; http://www.world clim.org/version2): annual F I G U R E 1 Map of the studied sites (n = 10), with localities of sampled "biodiversity plots" (n = 15) and "normal plots" (n = 67) [Colour figure can be viewed at wileyonlinelibrary.com] DEMBICZ Et al. mean temperature, annual mean precipitation, mean temperature of the warmest quarter, mean temperature of the coldest quarter and precipitation of the wettest month.
To account for potential variation in species pools between sites due to different management history, we calculated the distance of each plot to the centre of the nearest settlement (henceforth "distance to settlement"). To account for differences in landscape structure, we calculated the share of grasslands in a 4-km 2 circle surrounding each plot. We took the data on vegetation cover from CORINE Land Cover 2012 v. 20 (Copernicus Land Monitoring Service, 2012). The circle (buffer) had a 1,128-m radius (which had been used to explain grassland richness patterns by Janišová et al., 2014). For the purpose of this study, we considered the following non-forest CORINE land-cover types to represent grasslands: "pastures" (code 231), "land principally occupied by agriculture, with significant areas of natural vegetation" (243), "natural grasslands" (321), "beaches, dunes, sands" (331), "bare rocks" (332) and "sparsely vegetated areas" (333). The calculations were done in QGIS 3.12 software (QGIS Development Team, 2020).

| Statistical analyses of species richnessenvironment relationships
Prior to statistical analyses, we checked the distribution of all environmental variables for the 10-m 2 plots (Appendix S2) and transformed them if they strongly deviated from normal distribution (for details see Appendix S2). If predictors were highly correlated (Pearson's │r│ ≥0.7; Dormann et al., 2013), we selected the more ecologically relevant one to avoid multi-collinearity (see Appendix S3). Our final selection included eleven predictors: distance to settlement, share of grasslands in a 4-km 2 circle, elevation, inclination, heat-load index, pH, soil conductivity, sand proportion, rock cover, litter cover and grazing intensity (Appendix S2).
To avoid large differences in the variances among factors and to improve model convergence, we first rescaled all continuous variables to a range of 0-1 using a min-max scalar. For six variables (pH, heat-load index, grazing intensity, conductivity, sand proportion and elevation), the addition of their quadratic terms considerably improved model fit in the univariate generalised linear models (GLM) with Poisson error distribution (ΔAIC c > 2 between the model with and without quadratic term). However, to avoid overfitting in our relatively small data set, we only used the quadratic terms of the two variables that improved the model performance most (ΔAIC c > 10 between the model with and without the quadratic term), i.e. pH and grazing intensity (Appendix S2).
To evaluate the influence of the selected environmental factors on species richness, we calculated their relative importance by multimodel inference in multiple regression models (Burnham and Anderson, 2002), using the MuMIn package (Bartoń, 2018) in the R software (R Core Team, 2017). The relative importance of each predictor was quantified as the sum of Akaike weights over all possible models containing the predictor (Burnham and Anderson, 2002).
Altogether, we built 13 separate models with species richness in 10 m 2 as the dependent variable, i.e. for total species richness (number of all species across all taxa), richness of each of the four taxonomic groups (vascular plants, non-vascular species, bryophytes, lichens), the five Raunkiaer life forms and the three functional groups of vascular plants (Table 1). For these models, we used data from 97 10-m 2 plots (30 corners of "biodiversity plots" and 67 "normal plots").
To account for potential spatial autocorrelation in the model residuals due to the nested structure of our sampling, we applied GLMMs using lme4 (Bates et al., 2015). We first run a GLM for total richness and compared it with GLMMs including different random factors, with respect to spatial autocorrelation in the model residuals (Moran's test, using the ape package; Paradis et al., 2014) and AICc values. The GLM residuals were significantly autocorrelated (p = 0.001; Appendix S4). Thus, we selected the GLMM variant with site and plot as random intercepts because this model had the lowest AICc and effectively removed the spatial autocorrelation (p = 0.398) from the residuals (Appendix S4).To allow for comparisons, we used the same random effect structure also for the 12 GLMMs with subsets of species. Finally, we calculated the amount of variance in richness explained by the full random-intercept model (conditional: R 2 GLMMc ), by the fixed effects only (marginal: R 2 GLMMm ) and by the random effects (random: (Bartoń, 2018).
To compare results between different spatial grains, we built Poisson regression models for total species richness as a response variable separately for each grain size (n = 30 for grains 0.0001-10 m 2 ; n = 15 for 100 m 2 ). Due to the much smaller sample size and to avoid model overfitting, we excluded the following variables from TA B L E 1 Species richness in10-m 2 plots (n = 97) for all taxa, as well as for different taxonomic groups, Raunkiaer's life forms and functional groups of vascular plants

Mean ± SD
All taxa 10 62 38.5 ± 11.1 Taxonomic groups the former set of eleven predictors: distance to settlement, sand proportion, rock cover (these three variables had the lowest importance values within variable categories presented in Appendix S2), and grazing intensity (as this variable was very unequally distributed in the biodiversity plots, e.g. there were no plots in grazing category 2).
In this case, we implemented GLMs with Poisson distribution of errors, as there was no significant spatial autocorrelation of model residuals. For each GLM model, we calculated the proportion of explained variance as adjusted R 2 using the rsq package (Zhang, 2018).
We tested for overdispersion associated with counts in Poisson regression models using overdisp in the blmeco package (Korner-Nievergelt et al., 2015), but we found no significant overdispersion in any of our GLM models.

| Analyses of species-area relationships and beta-diversity
We used the exponent (z-value) of power-law species-area relationships (SARs) to quantify beta-diversity of total plant species richness.
The z-value can be considered as the log-transformed multiplicative beta-diversity standardised by the relative increase in area A between the α-and γ-level: z = log 10 (S γ /S α )/log 10 (A γ /A α ), where S γ and S α are the values of mean species richness at all small grain sizes (Jurasinski et al., 2009;Polyakova et al., 2016).
Here, "overall" z-values were calculated fitting a linear regression to each nested-plot series (n = 15) in a double-log space: log 10 S = log 10 c + z log 10 A, with S = species richness, A = area, c and z = fitted parameters (Drakare et al., 2006;Polyakova et al., 2016;Dengler et al., 2020a). The z-values were subjected to the same multimodel inference as the alpha-diversity measures from the biodiversity plots.
In addition, we calculated z-values between two subsequent grain sizes ("local" z-values) and tested their dependence on grain size in a one-way repeated-measures ANOVA. A significant scale dependence of local z-values would mean that the overall SAR deviates from a perfect power law (Turtureanu et al., 2014; see also de Bello et al., 2007). Overall and local z-values were calculated separately for each biodiversity plot, averaging richness values of the two subplots with sizes smaller than 100 m 2 .

| Species richness across scales
In total, we recorded 454 taxa, including 380 vascular plants (84%)  3.1 ± 3.6 Note: that the plots of areas ≤ 10m 2 in each biodiversity plot were considered as independent observations. n, number of plots (*, within biodiversity plots; **, all normal plots); SD, standard deviation.

DEMBICZ Et al.
Mean total species richness per plot ranged from 2.7 species in 0.0001 m 2 to 65.2 species in 100 m 2 ( Table 2). Regardless of plot size, mean vascular plant species richness was on average higher than that of non-vascular species, and mean bryophyte species richness was higher than that of lichens. Maximum total species richness was six in 0.0001 m 2 , 62 in 10 m 2 , and 89 in 100 m 2 (Table 2).

| Diversity-environment relationships of different taxonomic groups, life forms and functional groups
Soil pH (unimodal relationship, with a peak at pH 6.5; Figure 2, Appendix S5) and inclination (negative) were the two strongest predictors for total species richness (Table 3). For species richness of vascular plants, the strongest predictors were (in decreasing importance): soil pH (unimodal), conductivity (negative), litter cover (positive), inclination (negative) and grazing intensity (slightly unimodal). Richness of non-vascular species declined strongly with increasing inclination and litter cover, but increased with increasing rock cover, sand proportion in the soil and share of grasslands in the surrounding landscape. Species richness of bryophytes and lichens had similar sets of strong predictors, but conductivity and litter cover were less important for bryophytes and rock cover, sand content and share of grasslands for lichens (Table 3).
The relationship between pH and species richness of hemicryptophytes, therophytes, legumes and "others" was also unimodal (Table 4). A negative relationship between richness and inclination was found for therophytes, legumes and "others". Richness of hemicryptophytes, cryptophytes, legumes and "others" was positively related to conductivity. Litter cover was positively related to the richness of hemicryptophytes and graminoids. Hemicryptophyte and graminoid richness showed a unimodal relationship with grazing intensity. Elevation, which was not important for vascular plants, was negatively related to the richness of (juvenile) phanerophytes and therophytes and positively related to cryptophyte richness.
Among the three elements of spatial nesting, site was by far the most relevant (ΔAICc = 16.55 compared to GLM), followed by plot ID (ΔAICc = 13.27), while adding region even decreased model quality (ΔAICc = -2.80) (Appendix S4). The best-fitting model included both site and plot ID as random intercepts. The variation explained by these two random factors was particularly high for total species richness (R 2 GLMMr = 0.33) and vascular plant richness (R 2 GLMMr = 0.42) ( Table 3).

| Diversity-environment relationships across spatial scales
The explanatory power of the GLM models was highest for 0.1-m 2 plots (R 2 adj = 27%), followed by 100-m 2 and 0.001-m 2 plots, while it was below 15% for all other grain sizes ( Figure 3, Table 5). The importance of the predictors changed considerably across grain sizes: the importance of pH and elevation increased towards larger grain sizes (Figure 3). By contrast, heat-load index and inclination were most influential at intermediate grain sizes

| Factors influencing diversity of all plants and vascular plants at 10 m 2
Soil pH was the most influential factor for species richness of vascular plants and for the majority of functional groups, and mostly showed a unimodal relationship (maximum at pH ≈6.5). This finding is in agreement with a range of other studies that found soil pH to be highly influential for vascular plant diversity in grass-  Kuzemko et al., 2016 in Ukraine) was found. There are two major explanations for these regional differences: first, the chance of finding significant relationships depends on the length of the environmental gradient studied, which is particularly true for unimodal relationships. According to Dupré et al. (2002) unimodal relationships between species richness and soil pH are hardly ever found when the gradient length is less than three pH units. The studies in dry grasslands that found unimodal relationships also covered longer pH gradients (Bulgaria: 3.6 pH units; Sweden: 4.6; Czech Republic: 3.6) vs. those covering shorter gradients (Siberia: 1.5 pH units; Romania: 2.3; Ukraine: 2.9; see above). Second, the pH-richness relationship is modified by other environmental factors, as shown by Palpurina et al. (2017): while the unimodal relationship was pronounced in the more humid western parts of the Palaearctic, it became insignificant or even negatively linear in the most continental parts. Generally, the effect of soil pH on species richness can be attributed to its physiological effect: Few species are adapted to very high pH because of physiological stress and limited uptake of nutrients and few species can cope with very acidic soils because of poor nutrient availabity (leaching/fixation of nutrients). Most importantly, a very low soil pH increases the solubility of phytotoxic elements, such as aluminium (Tyler, 2003).
Even if generally in Europe the species pools of base-rich soils are Note: Importance values ≥ 0.5 are in bold.The variances explained by the fixed effects (R 2 GLMMm ),the random-factors plot ID and site (R 2 GLMMr ) and the whole model(R 2 GLMMc ) are also given. Positive (+) and negative (−) relationships are indicated.In case of quadratic terms, a negative sign means a unimodal relationship, a positive sign a U-shaped relationship. Transformations applied: sqrt, square root; log10, decimal logarithm; arcsin of sqrt, arcsine of square root.

TA B L E 3
Relative importance of the terms of fixed-effect factors obtained from a multi-model inference of generalised linear mixed-effect models (GLMMs) fitted to total, vascular, nonvascular, bryophyte and lichen species richness in 10-m 2 plots (n = 97) DEMBICZ Et al.
larger than those of acidic soils (Ewald, 2003b), the largest number of species' niches can be assumed to overlap at intermediate levels of the overall range of pH values allowing plant growth (Schuster & Diekmann, 2003).
The fact that conductivity was the second-most influential parameter for the richness of vascular plants was quite unexpected, as salinity generally negatively influences plant species richness (Garcia, Marañón, Moreno, and Clemente, 1993). However, in our study the conductivity gradient was relatively short, i.e. besides one sample being strongly saline, all others were rather not to moderately saline. Thus, we suppose that the "true" factor influencing vascular plants species richness was not salinity, The third-most relevant factor for vascular plant species richness was litter cover, having a positive impact. This was unexpected as accumulated litter can also reduce the generative reproduction of vascular plant species requiring light for germination (Facelli & Pickett, 1991). Accordingly, a negative effect was found by Kuzemko et al. (2016;Ukraine), Mathar et al. (2016;Siberia) and Turtureanu et al. (2014;Romania). One could assume that in extremely arid situations, a thin to moderately thick litter layer might protect seedlings from desiccation and thus facilitate their establishment (e.g. Polyakova et al., 2016; Siberia). However, this explanation is not plausible here as our sites were rather less arid. Perhaps not litter itself was causal, but productivity, reflected by higher litter accumulation: among dry grasslands, the slightly more productive ones, which are typically the semi-dry ones, have more litter, while they are also more species rich.
The fourth-most relevant factor for vascular plant species richness was inclination with a negative impact. Our sampling included very steep slopes (up to 75°) that are subject to strong erosion, making it hard for plants to establish and to survive if they are not particularly deep-rooted.
Last among the factors relevant for vascular plant species richness in the region was grazing intensity with a slight unimodal relationship. This was to be expected according to the intermediate disturbance hypothesis (IDH;Connell, 1978) and according to the dynamic equilibrium model (DEM; Huston, 2014), particularly at our sites with intermediate productivity. Note: The variances explained by the fixed effects (R 2 GLMMm ), the random-factors plot ID and site (R 2 GLMMr ) and the whole model (R 2 GLMMc ) are also given. Raunkiaer'slifeforms: Ph, phanerophytes; Ch, chamaephytes; H, hemicryptophytes; Th, therophytes; Cr, cryptophytes. Importance values ≥0.5 are set in bold. Positive (+) and negative (−) relationships are indicated.In case of quadratic terms, a negative sign means a unimodal relationship, a positive sign a U-shaped relationship. Transformations applied: sqrt, square root, log10, decimal logarithm; asin of sqrt, arcsine of square root.

| Differences between functional groups of vascular plants at 10 m 2
In many respects, the diversity-environment relationships of life forms and functional groups were similar to that of vascular plants as a whole, which logically is particularly true for those groups that accounted for the major fraction of vascular plants (i.e. hemicryptophytes and "others"). Thus, we will discuss here only major deviations from the results of vascular plants as a whole.
Elevation had a strong negative impact on phanerophyte and therophyte richness, but positively affected cryptophytes. Among all life forms, phanerophytes are most connected with lowlands and F I G U R E 3 Explained variance (expressed as adjusted R 2 ) of the generalised linear models (GLMs) for total species richness (top left plot) and relative importance of predictors for total species richness at seven spatial grains (plot sizes) (n = 30 for 0.0001-10 m 2 ; n = 15 for 100 m 2 ) [Colour figure can be viewed at wileyonlinelibrary.com]  Note: Importance values ≥ 0.5 are set in bold. Plus signs (+) indicate positive relationships, minus signs (−) negative relationships. If the linear term of a predictor has a positive sign, but its quadratic term a negative sign, the overall relationship is unimodal. Transformation applied: sqrt, square root; log 10 , decimal logarithm; asin of sqrt, arcsine of square root.
their species pool strongly decreases with elevation, so it is not unexpected that this is also reflected in plot-scale richness. A decrease of the fraction of therophytes with elevation (from 1,000 to 4,000 m a.s.l.) and an increase of the fraction of cryptophytes (from 1,000 to ca. 3,000 m a.s.l.) was also found by Mahdavi et al. (2013) for steppic to alpine grasslands in the Alborz Mts., Iran.
Rock cover had contrasting effects on four of the life forms: phanerophyte and therophyte richness was influenced positively, while hemicryptophyte and cryptophyte richness declined with increasing rock cover. For therophytes, which are often less competitive, small-statured species, this factor can similarly act as for non-vascular species (see below), supporting low-competition microhabitats.
In contrast, rocky outcrops provide shelter for emerging phanerophytes against grazing.

| Factors influencing diversity of bryophytes and lichens at 10 m 2
Drivers of diversity of non-vascular species were different from those of vascular plants, likely because of their different ecology.
The only strong predictor with the same direction was inclination.
Conductivity and litter cover had opposite effects for non-vascular species vs. vascular plants: the higher the conductivity (and, correlated with this, the humus content) and the higher the litter cover, the lower was the richness of non-vascular species and particularly of lichens. As discussed above for vascular plants, both parameters are related to more productive sites and later succes- Moreover, richness of non-vascular species and particularly bryophytes was positively related to rock cover and sand proportion. One explanation might be the worse growth conditions for vascular plants (lower water-holding capacity, less rooting space) in rocky and sandy soils that improve conditions of drought-adapted bryophytes and lichens. Moreover, some mainly saxicolous species (which were not counted here) might happen to grow also on soil beside the rocks. A similar pattern was found in Romania (Turtureanu et al., 2014) and Ukraine (Kuzemko et al., 2016), while Löbel et al. (2006) reported an optimum of bryophyte and lichen diversity at a rock cover of 50%.
Surprisingly, the share of grasslands in a 4-km 2 circle positively affected bryophyte richness, despite it being no influential factor for vascular plants. As was suggested by Löbel et al. (2006), who obtained similar results in Sweden, dispersal of bryophytes can be limited by distance (especially in pleurocarpous species, often reproducing vegetatively). A higher share of grasslands in the landscape increases grassland connectivity and promotes species with short-dispersal abilities.

| Random factors and their implications
In the three previous regional studies using the same sampling design (Turtureanu et al., 2014;Kuzemko et al., 2016;Polyakova et al., 2016), the authors also had tested for spatial autocorrelation, but never found significant patterns and thus did not account for it. By contrast, in our study we found very strong spatial autocorrelation for which we accounted by using GLMMs instead of GLMs with different random factors acknowledging our clustered (i.e. non-random) sampling. It turned out that adding "site" (i.e. the identity of the grassland patch of typically several hundred to a few thousand metres in diameter) improved the models (as measured with AICc) most among the three nesting levels and already removed the spatial autocorrelation in the residuals. Adding also "plot ID" (i.e. accounting for the particularly close proximity of the two corners of biodiversity plots) led to a slight further improvement of the model, while adding "region" (Vratsa vs. Koprivshtitsa region) did not improve the models. This means that the similarities of species richness within regions were already fully covered by the included environmental parameters. Also similar richness values in the two 10-m 2 corners of the biodiversity plots were largely explained by similar environmental conditions.
However, species richness within sites was more similar and among sites more dissimilar than expected by the measured environmental parameters. The differences were quite pronounced, with more than twofold higher richness under otherwise identi- correlated with precipitation and negatively with the temperature; Appendix S3) was highly important for the two largest grain sizes. This is in agreement with Siefert et al.'s (2012) results, which suggested that climatic variables become more influential towards larger grain sizes. At 0.1 and 0.01 m 2 different soil parameters became influential (rock cover, sand proportion and conductivity), which also agrees with findings by Siefert et al. (2012). The topographic factors, i.e. inclination and heat-load index, played a large role in grain sizes of 0.01 to 10 m 2 , not fitting well into Siefert et al.'s (2012) dichotomy of climatic vs. edaphic factors. Lastly, for grain sizes above 0.01 m 2 , the explanatory power of the models was relatively high for 0.1 and 100 m 2 , but showed a minimum at 1-10 m 2 . This contrasts with findings of three previous studies with the same sampling approach that found a peak in this grainsize range (Baumann et al., 2016;Kuzemko et al., 2016;Polyakova et al., 2016). Also, theoretically one would expect stochasticity to decrease with grain size and thus explanatory power to increase.
Since the current data set for scale dependence is quite small, we refrain from entering into speculation about possible causes, but suggest that this phenomenon should rather be studied with larger data sets from multiple sites.

| Species-area relationships
With a mean z-value of 0.240 for all taxa combined, the stud- We could not find a significant variation in z-values across spatial scales, indicating that indeed the power law was a suitable model to describe the SARs in these grasslands. This is in line with the main results of Dengler et al. (2020a) who found the (normal) power law to be the generally most appropriate model across more than 2,000 nested-plot data sets from Palaearctic grasslands.
In our study, we found a positive relationship of z-values with elevation. By contrast, most studies on beta-diversity along elevational gradients report negative relationships, explained by smaller species pools (Kraft et al., 2011) and different community assembly mechanisms (Lazarina et al., 2019) at higher elevations.
However, one should bear in mind that we studied beta-diversity at much smaller scales (centimetres to metres) than typical beta-diversity studies; thus different mechanisms and patterns are likely.
For example, small-scale microclimatic differentiation is known to be much higher at higher elevation than in the lowlands (Scherrer & Körner, 2009). In addition, we found a positive effect of inclination, which may increase erosion and thus small-scale disturbance. Such disturbances create open patches of soil or rock increasing heterogeneity and thus fine-grain beta-diversity.

| Conclusions and outlook
With regard to alpha-and beta-diversity at small scales, Bulgarian dry grasslands are intermediate compared to grasslands across the Palaearctic. We found pronounced differences in drivers of plotscale diversity among the three taxonomic groups, among the eight functional vascular plant groups and across the seven grain sizes.
Such systematic differences are theoretically expected (e.g. Shmida & Wilson, 1985;Siefert et al., 2012) and have been demonstrated before in other case studies (Löbel et al., 2006;Turtureanu et al., 2014), but they call for caution when generalising the results of biodiversity studies. One needs to account for such differences, e.g. when developing measures to maintain or increase biodiversity in grasslands, as species groups might be differently affected by the same measure.
Focusing on the main driving factors of grassland diversity, our study is in line with other studies (e.g. the general unimodal relationship of vascular plant diversity at plot scale with soil pH and the positive relationship of non-vascular species diversity with rock cover), but also deviated in some cases. For others, the results from regional studies are much less consistent or even contradictory, pointing to regional idiosyncrasies and/or interactions of different factors. A particularly important finding was that in the studied grasslands a large part of the variation in species richness was not explained by any of the included topographic, landscape, soil, climate and land-use variables, but was associated with unmeasured site-level differences. To understand the reasons of such differences, the next steps should be conducting a common analysis across the diverse grasslands of the Palaearctic while at the same time adding more potential predictors (e.g. historical factors).

ACK N OWLED G EM ENTS
We would like to thank Iva Apostolova for co-organising the 3rd EDGG Research Expedition on the data of which this paper is built, her and Aslan Ünal for participating in the field work and Anna Ganeva for revising critical bryophyte species.

AUTH O R CO NTR I B UTI O N S
The study was conceptually planned by JD and practically organised by HP, KV and NV. All authors, except ID and SB, were involved in the vegetation sampling. ID planned and conducted the statistical analyses, HP prepared the map. The manuscript was mainly written by ID and JD, with major contributions by NV. All authors critically revised the whole paper and approved its content.

DATA AVA I L A B I L I T Y S TAT E M E N T
The data used in this study are stored in and available from the GrassPlot database (GIVD ID EU-00-003; Dengler et al., 2018). The data of the 10-m 2 plots are additionally stored in the Balkan Dry