Introduction

Temperate grasslands are the most diverse plant communities globally with regard to vascular plants for plot sizes below 100 m2 (Wilson et al. 2012). Except for two examples in Argentina, all known world-records of plant species richness of grasslands are located in the Palaearctic realm. While the Palaearctic steppes with an original extent of 13 million km2 form the largest grassland biome on Earth (Wesche and Treiber 2012), the documented extremely diverse grasslands are all located within the nemoral biome of Europe, meaning that they are semi-natural communities that have emerged and are maintained by low-intensity land use (Dengler et al. 2014; Roleček et al. 2014). It is challenging to understand why semi-natural analogues to the natural steppe vegetation should host higher species densities than the natural steppes. However, this finding could be an artefact caused by the disproportionally higher number of studies on plant diversity in Palaearctic semi-natural grasslands (e.g. Löbel et al. 2006; Merunková et al. 2012; Janišová et al. 2014; Turtureanu et al. 2014; Palpurina et al. 2015; reviews by Dengler 2005; Dengler et al. 2014) than in natural steppes (e.g. Chytrý et al. 2007; Mathar et al. 2016; Kuzemko et al. 2016, in this issue). Indeed, there have also been sporadic reports of species richness close to the global maxima from natural meadow steppes (Walter and Breckle 1986; Lysenko 2007), which might indicate that more comprehensive studies could find hotspots of small-scale diversity also in natural steppes.

Palaearctic steppes form one of the largest biomes on earth. They are very diverse concerning the physical habitat and floristic composition, ranging from the relatively mesic steppes of Eastern Europe to dry-continental steppes (including semi-desert and desert steppes) widely distributed across the Asian continent (Lavrenko et al. 1991; Korolyuk 2002). Many steppes were used for nomadic pastoralism over centuries or even millennia. Over the last few decades, however, vast areas formerly covered by natural grasslands were converted into arable fields (Hoekstra et al. 2005; Smelansky and Tishkov 2012; Kuzemko et al. 2016, in this issue). Today, Asian steppes are among the most threatened ecosystems of the world (Henwood 1998; Wesche and Treiber 2012). Studying patterns and processes of Eurasian steppe biodiversity is of great scientific relevance, because they may represent a direct continuation (or at least a modern analogue) of the Pleistocene glacial steppe that was widespread across the Palaearctic, and from which much of the floristic composition of present-day secondary grasslands is derived (Bredenkamp et al. 2002; Chytrý et al. 2007; Kuneš et al. 2015; Pokorný et al. 2015).

We still lack a general theory of the drivers of plant richness at fine spatial scales in steppe and grassland ecosystems. The environmental variables most frequently analysed to explain richness patterns at plot scale in grasslands and steppes include soil pH (e.g. Löbel et al. 2006; Chytrý et al. 2007; Palpurina et al. 2015), soil fertility (e.g. Becker and Brändel 2007), grazing intensity and land use (e.g. de Bello et al. 2007; Turtureanu et al. 2014; Mathar et al. 2016, in this issue), elevation (Austrheim 2002), microrelief (Löbel et al. 2006; Deák et al. 2015), topographic wetness and soil humidity (Moeslund et al. 2013; Mathar et al. 2016 in this issue) and climate (Chytrý et al. 2007; Palpurina et al. 2015). The effect and importance of these factors, however, significantly varied between studies. For instance, soil pH can have a positive (Kuzemko et al. 2016 in this issue and literature cited therein; see also Chytrý et al. 2010), negative (Chytrý et al. 2007), hump-shaped (Mathar et al. 2016 in this issue) or no effect (Cachovanová et al. 2012) on plant species richness of dry grasslands. Another example is the role of management and disturbance factors, which were found to be the most explanatory variables in some studies (e.g. Maurer et al. 2006; Pedashenko et al. 2013; Turtureanu et al. 2014) while in others it was physical habitat that had the greatest effect (Löbel et al. 2006; Chytrý et al. 2007; Moeslund et al. 2013; Merunková et al. 2014). One reason for heterogeneous findings may be the existence of interactions among environmental variables (Bakker et al. 2006; Chytrý et al. 2007; de Bello et al. 2007; Huston 2014). A second reason could be that many studies did not cover a sufficiently long gradient of the independent variable, thus missing the peak of hump-shaped relationships (Fraser et al. 2015). Further, the patterns of biodiversity–environment relationships may strongly depend on spatial scale (Auestad et al. 2008; de Bello et al. 2007; Turtureanu et al. 2014) or differ between α- and β-diversity (de Bello et al. 2007; Hanke et al. 2014). Finally, it should be stressed that while non-vascular plants (bryophytes and lichens) often contribute significantly to biodiversity levels of Eurasian dry grasslands (Dengler 2005; Löbel et al. 2006), their diversity–environment relationships are much less well-known and they seem to show contrasting responses to environmental variables compared to vascular plants (see Kuzemko et al. 2016 in this issue, for a brief review).

The aim of our study thus was to contribute to the general knowledge of plant diversity patterns in natural steppes of the Palaearctic. For this purpose, we selected the Republic of Khakassia in Central Siberia, Russia, which is particularly interesting as a study area because here the two main ecological-geographical types of Palaearctic steppes—the West Palaearctic steppes and the Central Asian steppes—co-occur in the same landscape(Lavrenko et al. 1991). The steppes of Khakassia still cover huge areas and comprise, among others, specific cryo-petrophytic communities (rocky steppes strongly influenced by frost in winter) with a high proportion of glacial relicts (Polyakova and Larionov 2012; Ermakov et al. 2014). The conservation relevance of these communities is also enhanced by the occurrence of a large number of red-list taxa at regional and national scale (Ankipovich et al. 2012). Previous studies on the vegetation of Khakassian steppes were mainly devoted to plant community classification (Kuminova et al. 1976; Makunina 2006; Ermakov 2012), conservation assessment of communities with Red Data Book species (Polyakova 2013) and large-scale mapping (Ermakov et al. 2014). One study addressed diversity–environment relationships of steppe, forest and tundra communities across along bioclimatic and edaphic gradient in a region adjacent to our study area (Chytrý et al. 2007). However this study was restricted to vascular plants and had a fixed plot size of 100 m2.

We sampled vascular plant, bryophyte and lichen composition in nested plots of seven different sizes (0.0001–100 m2) along topographic, edaphic and bioclimatic gradients, using the well-established multi-scale sampling approach of the EDGG Research Expeditions (see Dengler 2009b for a description of the methodology), to answer the following questions:(i) How rich in plant species are the natural steppes of Southern Siberia compared to natural and semi-natural grasslands in other regions of the Palaearctic? (ii) What are the main environmental drivers of the diversity patterns in these steppes? (iii) What are the diversity–environment relationships and do they vary between spatial scales and among different taxonomic groups?

Methods

Study area

Our study was conducted in the northern part of the Minusinskaya Basin (54°–55°N, 88°–90°E; Fig. 1), located in the Republic of Khakassia, Russia. This intermountain basin extends 180 km in northeast-southwest direction with a width of 60–100 km and is situated between the large mountain systems of the Kuznetskij Alatau in the West, Western Sayan in the South and Eastern Sayan in the East. The area has gently undulating landforms and elevation is between 150 and 700 m.a.s.l.

Fig. 1
figure 1

Location of Khakassia within Russia (a) and of the four sampling sites b within Khakassia: c1 surroundings of Itkul Lake, c2 around the village of Gol’dzha, c3 around the village of Katyushkino, c4 Podzaploty area. Squares indicate all sampled nested-plot series (n = 39) and circles additional 10-m2 plots (n = 54). Images c2 and c3 have the same scale as c1. Background images for c1c4 obtained from open-source Google Earth

Due to the “rain shadow” effect by the Kuznetskij Alatau, little of the westerly humid air masses reaches the basin, resulting in an ultra continental climate. Winters are extremely cold, while summers are warm (see Table 1 for mean temperatures).Annual precipitation is 300–550 mm (modelled data based on Gidrometeoizdat 1966–1970), with up to 80 % falling between April and October (Gavlina 1954). In winter, snow cover does not exceed 10–20 cm, and is unevenly distributed: wind-exposed slopes are often snow-free, leading to deeply frozen soils. The growing season is short, spanning from late May to mid-September, when there is the onset of frost events (Kuminova et al. 1976).

Table 1 Summary data of all environmental variables used in analyses prior to modelling and those included in the further modelling of species richness of the 10-m2 plots (n = 126) with an indication of addition of squared terms to modelling and the of data transformation applied to improve normality, if needed

The soils are characterized by low humus content, but high carbonate content (Nikolskaya 1968; Kuminova et al. 1976). The most common soil types in Khakhassia include haplic kastanozems chromic and haplic chernozemspachic, below and above 300 m.a.s.l., respectively (Gradoboev 1954; Kolyago 1967).

The natural steppe vegetation occupy 21.7 % of total area of Khakassia and is distributed in the steppe and the forest-steppe belts. The rest of the area is covered by forest (46 %), floodplain and semi-natural meadow (11.5 %), hight-mountain (10.8 %) and other types (10 %) of vegetation (Kuminova et al. 1976).The various plant communities of Khakassian steppes are traditionally assigned to two different phytosociological classes, the Euro-Siberian steppes of the Festuco-Brometea and the Central Asian steppes of the Cleistogenetea squarrosae (Korolyuk 2002). Cleistogenetea communities are characterized by a predominance of xerophytic Central Asian and Mongolian species, while in Khakassia, Festuco-Brometea communities represent meadow-steppes with prevailing mesophytic Eurasian species. Both classes co-occur in the landscapes of the Minusinskaya Basin, but occupy sites with different ecological and topographical characteristics. Zonal steppes of the Cleistogenetea predominate on flat areas and on the south-facing slopes. The Festuco-Brometea steppes prevail within locally more humid sites on northern slopes.

Khakassian steppes have been used as pastoral land by nomadic tribes since ancient times. The conversion of steppe into arable land began in the early twentieth century (Smelansky and Tishkov 2012), reaching a peak during 1960’s. For the last two decades, the proportion of the arable land in the region has decreased as many fields were abandoned during 1990’s. Nowadays, steppes, including “petrophytic communities” (i.e. rocky grasslands), occupy 47.5 % of the steppe belt and nearly 50 % of forest-steppe belt of Khakassia (Kuminova et al. 1976). Today grazing and mowing are the main land-use types within the remaining steppe areas (Kandalova and Lysanova 2010). There is no agricultural land use within the “Khakassky” nature reserve (2680 km2), in which seven of its nine subunits are located in the steppe and forest steppe zones (273 km2).

Plot selection and arrangement

Field work took place during the 6th Research Expedition of the European Dry Grassland Group (EDGG; see Vrahnakis et al. 2013) in the second half of July 2013 at the peak development of vegetation. Four main areas were sampled to cover as much as possible of the different habitats of steppe communities (Fig. 1). Within each of these areas we subjectively selected homogenous stands of the different physiognomic and ecological types of steppe vegetation that could be recognized in the field; replicates of the same type were sampled only if clearly separated, to avoid pseudo-replication. We only sampled in the steppe and forest-steppe zones. Our plots were concentrated in protected, unmanaged areas within the nature reserve, but some of them were located in moderately grazed land outside the protection zone.

We sampled 39 nested-plot series (“biodiversity plots”) following the basic concept proposed by Dengler (2009b) and implemented in previous EDGG Research Expeditions, such as Romania (Turtureanu et al. 2014), Ukraine (Kuzemko et al. 2016, in this issue) and Bulgaria (Pedashenko et al. 2013). Each biodiversity plot was a square of 100 m2, with square-shaped subplots of 0.0001, 0.001, 0.01, 0.1, 1 and 10 m2 arranged in nested series in two opposite corners. In addition, we sampled 54 quadrats (“normal plots”) of 10 m2, resulting in 132 plots at the 10-m2 scale for richness data. Due to missing environmental data, one 10-m2 plot in one “biodiversity plot” and four “normal plots” were excluded from the analyses of diversity–environment relations. For each nested plot and normal plot, we recorded geographic coordinates and elevation with a GPS device.

Field sampling

For plots of all sizes we recorded all species of vascular plants, terricolous bryophytes, lichens and macroscopic taxa of ‘algae’ using the any-part present system (see Dengler 2008), while critical specimens (mainly small cryptogams) were collected for further determination. Nomenclature follows Czerepanov (1995) for vascular plants, Ignatov et al. (2006) for mosses and Andreev et al. (1996) for lichens.

For all 10-m2 plots (normal plots + those placed within nested plots) also environmental parameters were recorded: slope aspect, slope inclination, relief position, microrelief, soil depth. We visually estimated percentage cover of litter. Microrelief was measured as the vertical deviation of the surface from an idealized plane across the plot. Soil depth was measured as the maximum depth reached with an iron pin (diameter of 10 mm, length of 80 cm) at five random points within the plot. The median and the range of soil depth (the difference of maximum and minimum measurement within plot) were used in the analyses. Soil sample was taken to measure pH (in water), content of humus, CaCO3 and soil skeleton content (%) of a mixed sample from five randomly placed sub-samples from the upper 10 cm layer. We are aware that the microhabitat conditions in plots of sizes other than 10 m2 could differ from the measured values, but we nevertheless use those values to model richness patterns within all sampled spatial scales. We assume that within-plot variation of measured factors is smaller than between-plot variation. Moreover, measuring all microhabitat variables within all spatial grains would not have been feasible.

The plot data are stored in and available from the Database SpeciesArea Relationships in Palaearctic Grasslands (EU-00-003; Dengler et al. 2012), registered in the Global Index of Vegetation-Plot Databases (GIVD; Dengler et al. 2011).

Extraction of climatic variables

Heat load index was calculated for each plot as a proxy of microclimate (especially of dryness), using slope aspect and inclination (Olsson et al. 2009). Mesoclimate did not vary much between our sites despite the geographical distance, thus we used only five basic parameters mean annual temperature, mean annual precipitation, January mean temperature, July mean temperature and evaporation during the warm period, all extracted in ArcGIS 10.2 from a climatic model. The model was obtained by integrating weather data from Russian climate stations (Gidrometeoizdat 1966–1970) with elevation information from a digitized 1:200,000 topographic map. Temperature values were interpolated across the elevation gradient according to the adiabatic lapse rate of 0.65 K per 100 m of elevation. Precipitation was computed from precipitation-altitude charts compiled for each of the bioclimatic sectors of the Altai-Sayan region, defined on the basis of their precipitation regime (Polikarpov et al. 1986). Evapotranspiration was calculated by the climatic model on the basis of temperature and air humidity during the growing season. We decided against using the WorldClim precipitation data (Hijmans et al. 2005; http://www.worldclim.org/), because they appeared unrealistic for our study area when compared to local sources.

Analyses of richness–environment relationships

Before assessing richness–environment relationships, were selected 13 variables out of 18 that did not exhibiting strong collinearity among each other (Pearson linear correlation |r| < 0.7, Table 1).To test whether data were normally distributed, data on each environmental variable was checked visually using histogram, QQ plot and box plot graphs. If the distribution deviated strongly from normality, an appropriate data transformation was applied (for details see Table 1). One more 10-m2 plot was removed from further analyses at the stage of model preparation as it disproportionally influenced model results, as indicated by partial leverage plots (the heat load index value of this plot was an extreme outlier and had both high leverage and high influence in the regression). For two explanatory variables, humus content and mean annual precipitation, squared terms were additionally included in the analyses as those improved model performance in the univariate generalized linear models with Poisson error distribution that we conducted prior to the main analyses (ΔAICc > 2 when comparing model with and without quadratic term, Table 1).

In the main analyses, environmental variables explaining species richness patterns were analysed using a multi model interference based on generalized linear models (GLMs) with Poisson error distribution. We used Akaike information criterion for small sample size (AICc; Burnham and Anderson 2002) as measure of model quality. The estimation of variable importance was based on multimodel inference (Burnham and Anderson 2002) as implemented in the R-package MuMIn version 1.14.; Bartoń 2015). This approach calculates AICc values for the full model and all partial models derived by any possible combination of 1 to all variables of the full model. The relative importance of each variable is then calculated by summing the Akaike weights over all models that contain the specific variable (Johnson and Omland 2004). Importance values range from 0 to 1, with values close to 1 indicating that a variable is included in all statistically plausible models.

Analyses of richness–environment relationships were finally carried out for data collected within plots of the area 10 m2 (n = 126) for overall species richness and separately for different taxonomic groups as vascular plants, non-vascular plants (bryophytes and lichens together) and bryophytes and lichens separately. Moreover, to identify differences in important predictors between spatial scales, nested-plot series (data collected from spatial scales from 0.0001 m2 to 100 m2, n = 77 per scale but 39 for 100 m2) were analysed for total species richness, vascular and non-vascular plants.

The assumption of equidispersion was evaluated using function dispersion test R-package AER version 1.2-4 tests (Kleiber and Zeileis 2008). As significant overdispersion was found particularly for overall species richness, all analyses were repeated using negative binomial generalized linear models. Results did not differ qualitatively from those of the GLMs (Online Resource 1). No significant spatial autocorrelation was detected in the residuals of any model when tested using function moran.test in R-package spdep version 0.5-88 (Bivand and Piras 2015). McFadden’s pseudo-R 2 (1 − (log likelihood of the full model/log likelihood of the null model)) was calculated for the “best” (in terms of AICc) GLMs in order to compare the amount of variability explained by the response variables between different scales. Note that McFadden’s pseudo-R 2 tends to report much smaller values when compared to common R 2 measures (see Online Resource 2 for demonstration of similarity in patterns of different R 2 measures).

Analyses of species–area relationships and β-diversity

We used the “slopes” of the species–area relationships (SARs), to quantify β-diversity (Riccota et al. 2002; Jurasinski et al. 2009). Specifically, we used the so-called z-values, one of the parameters of power-law SARs. The power law (S = c A z, with S = species richness, A = area, c and z being fitted parameters) has been shown to be to a good general approximation for the increase of species richness with area at any spatial scale (Connor and McCoy 1979; Dengler 2009a), and the z-value is widely used to compare spatial species turnover between different ecological situations (Drakare et al. 2006; de Bello et al. 2007). The z-values can be seen as the logarithmic transformation of multiplicative β-diversity (β mult = S γ /S α , where S γ and S α are the mean species richness values at the coarse and at the fine grain sizes, respectively; Jurasinski et al. 2009), standardised by the increase in grain size between the α- and γ-level, i.e.:

$$z = {{\log_{10} \left( {\frac{{S_{\gamma } }}{{S_{\alpha } }}} \right)} \mathord{\left/ {\vphantom {{\log_{10} \left( {\frac{{S_{\gamma } }}{{S_{\alpha } }}} \right)} {\log_{10} \left( {\frac{{A_{\gamma } }}{{A_{\alpha } }}} \right)}}} \right. \kern-0pt} {\log_{10} \left( {\frac{{A_{\gamma } }}{{A_{\alpha } }}} \right)}}$$

We calculated the z-values in double-log space, i.e. log10 S = log10 c + z log10 A with linear regression (Drakare et al. 2006; Dengler 2009a), resulting in an overall z-value. Second, we calculated “local”z-values between two subsequent plot sizes to test whether β-diversity depends on spatial scale, i.e. the actual relationship deviates from a perfect power law (Turtureanu et al. 2014; see also de Bello et al. 2007). Both calculations were performed separately for each of the 39 biodiversity plots, using for each grain size the mean richness values of the two corners.

Results

Species richness at different spatial scales

In the whole dataset, we identified 498 taxa: 352 vascular plants (71 %), 75 bryophytes (15 %), 69 lichens (14 %) and two macroscopic taxa of terricolous “algae”. Among vascular plants, the most frequent graminoid species were Stipa krylovii, Helictotrichon schellianum and Carex humilis, and the most frequent forbs were Thalictrum foetidum, Aster alpinus, Bupleurum scorzonerifolium, Galium verum, Hedysarum gmelinii and Schizonepeta multifida (Online Resource 3). Regarding non-vascular plants, the most frequent species were Cladonia pyxidata, Peltulaeuploca, Barbula unguiculata, Bryum caespiticium and Weissia condensa.

Mean total species richness varied between 2.3 species at the smallest spatial scale (0.0001 m2) and 73.4 species at 100 m2 scale (Table 2). Species richness of vascular plants was evenly distributed among the plots, in contrast to non-vascular plants, for which the standard deviation of species richness exceeded the mean up to a plot size of 0.1 m2. The maximum recorded species richness for all taxonomic groups was 5 species in 0.001 m2, 80 species in 10 m2, and 106 species on 100 m2 plots (for vascular plants the numbers are 5, 72 and 94, and for non-vascular plants 3, 22 and 26, respectively) (Table 2). At the plot scale, the average proportion of vascular plants was higher and varied between 84 % (in the 0.01-m2 plots) and 92 % (in the 0.0001-m2 plots). The species richness of non-vascular plants was lower than that of vascular plants at all spatial scales. The species richness of bryophytes exceeded that of lichens only at smaller spatial scales (Table 2).

Table 2 Overview of the species richness encountered in the all nested-plot series and normal plots of the study

Taxon-dependence of species richness–environment relationships

Species richness in the 10-m2 plots (n = 126) was controlled predominantly by heat load index (negative linear relationship), mean annual precipitation (hump-shaped relationship) humus content (hump-shaped relationship) and soil skeleton content (positive linear relationship) (Fig. 2a; Table 3).

Fig. 2
figure 2

Effect of heat load index, humus, pH and mean annual precipitation on a total richness, b vascular plant, c bryophyte and d lichen richness within 10-m2 plots (n = 126). Response surfaces were calculated based on full models (GLM with Poisson family error) using function visreg2d in R package visreg version 2.2-0 (Breheny and Burchett 2015)

Table 3 Estimates of the relative importance of predictor variables in the multiple regression models of total, vascular, non-vascular, bryophyte and lichen species richness within 10-m2 plots (n = 126) assessed with multimodel inference

Vascular plant richness followed the same patterns as total richness with respect to humus content, mean annual precipitation, heat load index and soil skeleton content. However, vascular plant species richness was also affected by relief (higher species richness at the ‘bottom’ than on ‘slopes’ and ‘ridges’) and was positively correlated with cover of litter, pH and range of soil depth (Fig. 2b; Table 3).

Species richness of bryophytes, similar to total richness and vascular plants, was also affected by heat load index (negatively) and relief (but with higher species richness on ‘slopes’ and ‘ridges’ than at the ‘bottoms’ of valleys). On the other hand it was positively associated with mean annual precipitation and mean annual temperature. The relationship of bryophyte species richness with humus content was hump-shaped (Fig. 2c; Table 3).

Species richness of lichens was strongly negatively associated with soil depth, cover of litter, microrelief, and humus content. Similarly to bryophytes, lichens species richness was higher on ridges than in valley bottoms. Their relationship with mean annual precipitation was hump-shaped, with the maximum richness at mean precipitation values (around 450 mm) (Fig. 2d; Table 3).

Scale-dependence of richness–environment relationships

The explanatory power of the regression models for the 77 series of nested plots was highest at a plot size of 100 m2, with a secondary peak for 10-m2 plots, gradually decreasing towards the smallest scales (Fig. 3).

Fig. 3
figure 3

Change of the explanatory power of the regression model for total species richness across the studied spatial scales from 0.0001 to 100 m2 (n = 77 for 0.0001–10 m2; n = 39 for 100 m2 and z-values). Note that McFadden’s pseudo-R 2 is not comparable to common R 2 measures as it tends to be small

Heat load index was the most important factor (with a negative correlation), for total species richness across most of the spatial scales with a particularly high significance for 0.01, 1, and 10 m2 (Table 4, Online Resource 4). Microrelief was important at scales from 0.01 to 0.1 m2 (with a negative correlation), while soil skeleton content had an increasingly positive effect with increasing plot size, and became very important at scales from 10 to 100 m2. Relief was most important at the scale of 0.1 m2, and was of negligible importance for the total richness at the largest and smallest spatial scales. As for climatic variables, annual precipitation was important at most spatial scales. It exhibited a hump-shaped relationship with total species richness at small scales and a strongly positive relationship at scales from 10 to 100 m2. Mean annual temperature was (moderately) important only at the scale of 1 m2. Presence of grazing had a positive influence on total species richness only at the scale of 0.1 m2 (Table 4, Online Resource 4).

Table 4 Estimates of relative parameter importance for the seven spatial scales of α-diversity and of the overall z-value as measure of β-diversity based on multimodel inference (n = 77 for 0.0001–10 m2; n = 39 for 100 m2 and overall z-values)

Species–area relationships and β-diversity

The overall z-values for all taxonomic groups and across all plot sizes ranged from 0.197 to 0.345, with a mean of 0.258 in the 39 nested-plot series. Among the studied factors, the overall z-values were positively related to microrelief had a u-shaped relationship with differences in mean annual precipitation (Table 4).The local z-values were significantly scale-dependent (repeated measures ANOVA: p < 0.001). The highest mean z-values were found at the smallest spatial scales (with a maximum of 0.322 at the 0.01–0.1 m2 transition) and gradually decreased above these plot sizes (reaching only 0.166 for the transition 10–100 m2) (Fig. 4). Maximum variance in the local z-values was found at the smallest plot sizes, where also the highest absolute z-values were recorded (Fig. 4).

Fig. 4
figure 4

Changes of species turnover or β-diversity (expressed asz-values) with spatial scale in the 39 biodiversity plots. Box plots of local z-values for each pair of subsequent plot sizes are shown. Letters indicate homogenous groups according to Tukey’s post hoc test. Medians and interquartile ranges (IQR) are shown by boxes; whiskers indicate ranges from Q1 − 1.5(IQR) to Q3 + 1.5(IQR)

Discussion

Diversity of Khakassian steppes in comparison with other regions

In general, our mean and maximum species richness values were in the range found in other Palaearctic dry grasslands (reviewed by Kuzemko et al. 2016 in this issue). At all spatial scales between 0.01 and 100 m2, our maximum as well as mean vascular plant species richness were higher than those of semi-natural dry grasslands of Europe, with the exception of those in Transylvania (Turtureanu et al. 2014). Compared to other natural Palaearctic steppes, those in Khakassia were much richer at 100 m2 than those of the Russian province of Tyuman (27 species reported by Mathar et al. 2016, this issue).

Unlike dry grasslands in Central Podolia (Kuzemko et al. 2016 in this issue) non-vascular plants were present at all scales. Non-vascular plants were generally most frequent and abundant in sites where abundance of vascular plants was limited by environmental stress, especially in the petrophytic steppe, which is in accordance with other studies demonstrating the competitive superiority of vascular over the non-vascular plants (Wal et al. 2005).

Richness–environment relationships

In our dataset, the environmental variable most markedly controlling species richness of vascular plants (and, less strongly, of non-vascular plants, especially lichens) was heat load index, with a negative correlation. This is very similar to the findings by Kuzemko et al. (2016, in this issue) for the Podolian steppes in Ukraine and by Turtureanu et al. (2014) for Transylvanian semi-natural grasslands. For Podolia, this finding was explained in that south-facing, steep slopes in a summer-warm climate give rise to a habitat with a marked drought stress that prevents colonization by many species of the regional species pool (Kuzemko et al. 2016, in this issue). In the extreme continental climate of Khakassia, the causal mechanism for the negative correlation between richness and heat load could also take place in spring, when solar radiation on south-facing slopes is already high but temperatures are still low and the soil is still frozen, leading to strong physiological drought (see Jarvis and Linder 2000).

Vascular plant richness had also a very strong negative correlation with convex topographic positions: diversity was much higher in bottoms than on ridges. Contrary to a positive correlation with convex topographic position found in Danish grassland diversity (Moeslund et al. 2013), our results can be explained again with the ultra continental climate of the study area, characterized by very low winter temperatures and a very little amount of snowfall. During the winter, the wind can blow the snow cover from ridges, leaving the soil and vegetation exposed to a long period of severe frost (about five months), with the formation of a layer of frozen ground (up to 2.5 m deep) with distinct cryogenic processes. This leads to patches of cryo-petrophytic steppe poorer in vascular plants and richer in bryophytes and lichens. In contrast, the presence of even a shallow snow layer in valley bottoms and concave micro-topographic position is favourable for the formation of moderately thermophilous steppes rich in forbs (Ermakov et al. 2014). Moreover, ridges are strongly exposed to erosion, resulting in shallow and nutrient-poor soil leading to unfavourable habitats for many vascular plants, while soil will accumulate in depressions.

The relationship between humus content and vascular plant richness was hump-shaped with the highest richness shifted towards high humus content. This finding might be connected with the productivity–richness relationship, as the humus content is a proxy of soil fertility, and of soil water retention, and thus aboveground productivity. The hump-shaped relationship between productivity and vascular plant species richness has recently been confirmed by Fraser et al. (2015). Given the relatively low productivity of the Russian steppes on a global scale (Titlyanova 1991), we can assume that most of our samples were still below the tipping-point of the hump-shaped response that can be found over a large productivity gradient (Fraser et al. 2015). They were thus below levels where competition for light starts to limit vascular plant richness (see Grime 2002; see also Mathar et al. 2016, in this issue), except on the most humus-rich soils (see Fig. 2). The above-mentioned explanations linking diversity patterns of vascular plants within Khakassian dry grasslands with productivity–diversity relationship are additionally supported by the observed strong positive relationship of litter cover with vascular plants species richness, as production of litter is a direct effect of primary productivity (Fraser et al. 2015).

Interestingly, soil pH had a strong positive effect on vascular plants species richness. The positive effect of the pH on vascular plants species richness is also remarkable (and currently unexplained) because pH was, in general, very high (mean 8.57, max. 9.25). The richness-pH relationship has to be inherently unimodal (Chytrý et al. 2007), and for such high pH values, other studies found that vascular plants richness decreased with increasing pH, in European secondary grasslands (e.g. Becker and Brändel 2007), in South Siberian steppes (Chytrý et al. 2007) and in West Siberian meadows and steppes (Mathar et al. 2016, in this issue).

Among the considered mesoclimatic factors, the only one that was important for most of the studied groups was mean annual precipitation, contrary to other studies that underlined the role of annual temperature (e.g. Palpurina et al. 2015), this can probably be explained because the temperature gradient within our study area was very short. The strong hump-shaped relationship of total richness (and that of vascular plants) with precipitation can stem from the same principle as the relationship with humus content. The primary productivity of temperate dry grasslands can be controlled by water to the same extent as nutrient deficiency (Knapp et al. 2001). At high precipitation values, the most competitive and large-sized species (suppressed elsewhere by drought stress), can take dominance, limiting plot richness. As within Khakassian grasslands, the flattening of the precipitation-vascular plant species richness curve happened only at the highest values of precipitation: over most of the gradient, the present results can be compared with those of Adler and Levine (2007) or Chytrý et al. (2007) (for steppe vegetation), who described positive relationship of mean annual precipitation with plant species richness.

Scale-dependence of richness–environment relationships

Heat load index was more important at medium to large plot sizes, which is in accordance with the results of Turtureanu et al. (2014) from Transylvania and Kuzemko et al. (2016, in this issue) from Ukraine. Humus content was the most important factor between 0.0001 to 0.1 m2 plot sizes in Transylvania, and between 0.1 and 1 m2 in Ukraine; in Khakassia, this predictor variable was most important between 0.1 and 100 m2. These results may contrast with the explanation given by Turtureanu et al. (2014) that resource availability is most important at small spatial scales (see Siefert et al. 2012).

An important positive role of microrelief in the richness model was noted in Transylvania between 0.1–100 m2, which turned negative only for the smallest plot size, consistently with heterogeneity-diversity theory (Tamme et al. 2010). By contrast, in the Khakassian steppes, the influence of microrelief was weakly positive at the smallest scales and strongly negative at plot scales 0.01–1 m2, becoming positive only above 10 m2 (Table 4, Online Resource 4).

Another interesting feature of our results is the strong importance of mean annual precipitation for all plots larger than 0.0001 m2 (Table 4, Online Resource 4). By contrast, both Kuzemko et al. (2016, in this issue) and Turtureanu et al. (2014) found mesoclimatic factors important only at larger spatial scales and gaining importance with scale, which is in accordance with results of meta-analytical studies by Siefert et al. (2012) and Field et al. (2009). The contrasting results of our study can be attributed to the harsh climatic conditions of the study area, where drought and associated lack of snow cover during winter are important stress factors (see Chytrý et al. 2007), which can act at different spatial scales (for example through cryogenic processes).

It is interesting to note that disturbance (grazing, mowing) was found in similar studies as an important predictor variable (Löbel et al. 2006; Turtureanu et al. 2014; Kuzemko et al. 2016, in this issue) positively influencing species richness at different spatial scales. For instance, in Transylvania it had a strong positive effect across all plot sizes between 0.1 and 100 m2, but in the Khakassian dataset it was important for total species richness only at 0.1 m2. At the 10-m2 scale its influence was rather weak, but still positive, for all taxonomic groups except lichens. The influence of grazing on grassland species richness strongly depends on its intensity (Dupré and Diekmann 2001). Within the area sampled in the present study, pastoral management was very low intensity (with free-ranging cows). Such type of management can resemble grazing by the wild animals that shaped the steppes of Khakassia from their very beginnings (Akhmetyev et al. 2005). On the other hand, as steppe is the zonal vegetation type within Khakassian lowlands the lack of human management does not result in rapid successional changes and loss of biodiversity, as is the case in semi-natural grasslands of Western and Central Europe (Pykälä et al. 2005). Additionally, as management influence was not the focal question of our study, preferential plot selection avoiding intensively grazed areas shortened the gradient of land use intensity, which could reduce its importance as a factor influencing species richness.

Species–area relationships and β-diversity

The z-values as measure of β-diversity were higher in Khakassia than those found in a review of European dry grassland associations (Dengler 2005) or in Podolian steppes (Kuzemko et al. 2016, in this issue), and nearly as high as in Transylvanian dry grasslands (Turtureanu et al. 2014). The scattered knowledge about z-values in different grassland ecosystem currently hinders an explanation of this outstanding position. In contrast to Estonian dry grasslands(Dengler and Boch 2008) and Ukrainian steppes (Kuzemko et al. 2016 in this issue), we found a significant scale-dependency of the local z-values, with a shape of the relationship (Fig. 4) similar to that reported by Turtureanu et al. (2014) for Transylvania. We observed an already high value at the smallest transition (0.0001–0.001 m2), gradually increasing up to the maximum at the 0.01–0.1 m2 transition (this is the same spatial scale showing the maximum z-value in Transylvania).At larger spatial scales, the mean z-values showed a steep decrease, reaching 0.166 for the transition 10–100 m 2(strikingly similar to the 0.165 value found in Transylvania for the same transition).

As expected, the environmental variable that had a significant positive influence on overall z-values was microrelief. This can be viewed as a proxy of micro-ecological heterogeneity within the larger-scale plots (although it was measured only within the 10-m2 plot),i.e. of the probability that new, different microhabitats will be included in sampling with increasing plot size, leading to the inclusion of plant species linked to those microsites.

Differences in diversity patterns between vascular plants and non-vascular plants

Consistent with previous studies (see Turtureanu et al. 2014, and references therein), we found that the richness patterns and drivers of bryophytes and lichens were different than those of vascular plants. For instance, lichen species richness was negatively affected by litter cover and humus content, while vascular plants species richness was generally positively correlated with these factors. Both bryophytes and lichens were more diverse on ridges, while vascular plants were more diverse in valley bottoms. These patterns can be explained by the lack of snow-cover on ridges, exposing the soil to extremely cold temperatures in winter and to drought stress in spring, thus favouring non-vascular plants. The latter instead outcompete bryophytes and lichens on the snow-covered, deep and humus-rich soils of the valley bottoms.

We found that soil pH had a strong positive correlation with vascular plant richness but, contrary to some previous studies (e.g. Löbel et al. 2006), only very weak relationship with non-vascular plant diversity. However, this latter result is similar to the findings of Turtureanu et al. (2014) in Transylvanian secondary grasslands. That phenomenon can be explained by a relatively short pH gradient, where all sampled sites could be regarded as calcium rich, as is the case in the Transylvanian grasslands (Turtureanu et al. 2014).

Bryophytes were the only group for which species richness was strongly correlated with mean annual temperature (positively) and at the same time strongly positively influenced by mean annual precipitation. As bryophytes are not rooted and thus susceptible to desiccation, their positive relationship with annual precipitation is understandable (see Pharo et al. 2005; Turtureanu et al. 2014). Low winter temperatures have also been found to indirectly negatively influence bryophyte richness through low vapour pressure (Pharo et al. 2005).

Conclusions

In comparison to other Palaearctic grasslands studied with a comparable methodology, Khakassian steppes ranked quite highly for both α-diversity (measured by the mean and maximum richness at each plot size) and β-diversity (measured as the overall and local z-values). We also found a clear peak of local z-values between 0.001 and 0.1 m2, which is currently unexplained.

Consistent with theory and previous studies, the importance of different factors in the regression models explaining richness varied systematically with spatial grain. Moreover, the analysed drivers varied in their importance and sometimes even had contrasting effects on vascular versus non-vascular plants, and partly also between bryophytes and lichens. This indicates the importance of considering groups other than just vascular plants, especially in fundamental ecological studies. Our research confirmed that the steppes of Southern Siberia are, at small spatial scales, one of the plant biodiversity hotspots of the Palearctic, which strengthens the need for the conservation of natural ecosystems in this region.