Species–area relationships in continuous vegetation: Evidence from Palaearctic grasslands

Species–area relationships (SARs) are fundamental scaling laws in ecology although their shape is still disputed. At larger areas, power laws best represent SARs. Yet, it remains unclear whether SARs follow other shapes at finer spatial grains in continuous vegetation. We asked which function describes SARs best at small grains and explored how sampling methodology or the environment influence SAR shape.

Historically, studies of SARs have largely been restricted to two functions, (a) the power function (often called the power law; Arrhenius, 1921;Preston, 1962) and (b) the logarithmic function (sometimes erroneously termed the 'exponential' function; Gleason, 1922). This was mainly because the fit of these two functions was easily explored using least squares linear regression techniques.
A comparison of a broader set of functions became possible with the advent of nonlinear regression techniques (e.g. Dengler, 2009;Flather, 1996;Guilhaumon, Gimenez, Gaston, & Mouillot, 2008;Stiles & Scheiner, 2007). In recent years, a wide array of different function types has been proposed and tested (Dengler, 2009(Dengler, , 2010Tjørve, 2003Tjørve, , 2009). Consequently, several comprehensive studies have been conducted on the fit of different functions and parameters for island SARs as well as other broadscale SARs. Triantis et al. (2012) compared 20 different models for 601 true island datasets around the world and found strong support for the power function overall. Matthews, Guilhaumon, Triantis, Borregaard, and Whittaker (2016) extended this study to 182 habitat islands, with a similar finding. In a further step, they tested how ecological context affects the slope parameter of the power function, and they found systematic differences between island types and spatial scales, but not between major taxa.
While knowledge of functions and parameters of island SARs has been broadly consolidated during the last decade, comparable empirical evidence on small-grain SARs in continuous habitats is still lacking (for theory see Storch, 2016;Williamson, 2003). With continuous habitat or vegetation, we refer to situations where the sampling units do not have a natural border such as islands or habitat islands, but are delimited by the researcher. The influential study of Crawley and Harral (2001) on how biodiversity depends on scale in continuous vegetation, for example, a priori only considered the power function. Some regional studies have found a prevalence of the power function using multimodel inference, but were restricted to less than 20 datasets (e.g. Dengler, 2009;Dengler & Boch, 2008). Furthermore, Rosindell and Cornell (2007) obtained power function SARs from a spatially explicit ecological drift model (Hubbell, 2001) within a homogeneous grid model assuming skewed dispersal kernels. By contrast, there is a belief that the logarithmic function should be more suitable at small spatial scales (Gleason, 1922;van der Maarel, 1997).
Saturated functions (i.e. functions with a horizontal upper asymptote) are also often assumed to represent SARs in continuous vegetation weöö, inspired by the still widespread, but flawed (see Barkman, 1989) concept of so-called 'minimal areas ' (e.g. Mueller-Dombois & Ellenberg, 1974), which was assumed to be the scale at which species richness is saturated for a given community. Additional confusion comprehensive studies of SARs at different grain sizes and for different taxa, thus supporting the general appropriateness of the power function for modelling species diversity over a wide range of grain sizes. The poor performance of the Michaelis-Menten function demonstrates that richness within plant communities generally does not approach any saturation, thus calling into question the concept of minimal area.

K E Y W O R D S
logarithmic function, Michaelis-Menten function, minimal area, nested-plot sampling, nonlinear regression, Palaearctic grassland, plant biodiversity, power law, scaling law, speciesarea relationship (SAR) Handling Editor: Holger Kreft around small-grain SARs was caused when confounding different sampling schemes with SARs in the strict sense (i.e. those originally considered by Arrhenius, 1921, or Preston, 1962. For example, Stiles and Scheiner (2007) and DeMalach et al. (2019) reported that the logistic function (a saturated function) and the logarithmic function, respectively, performed much better than the power function.
However, they had analysed species accumulation curves, merging non-contiguous sample units (also called species-sampling relationships, SSRs; see Dengler, 2009;Fridley, Peet, van der Maarel, & Willems, 2006), and not SARs in the strict sense. In conclusion, this situation calls for a comprehensive, multimodel inference analysis of small-grain SARs in continuous vegetation, comparable to the analyses of Triantis et al. (2012) and Matthews et al. (2016) for island SARs.
As the Palaearctic biogeographic realm comprises more than one third of the worlds' ice-free terrestrial surface and spans a wide range of climatic and topographic gradients, it harbours a high number of vegetation types and considerable biodiversity (Rounsevell, Fischer, Torre-Marin Rando, & Mader, 2018). Around 22% of the Palaearctic is composed of various grassland types (Török & Dengler, 2018), some of them being the world record holders of small-grain vascular plant diversity (Wilson et al., 2012). A large proportion of the Palaearctic grasslands are primary grasslands such as steppes and arctic-alpine grasslands. Even in regions where the potential vegetation is forest, natural grasslands occur in azonal and extrazonal conditions. Moreover, agricultural activities and pastoralism long present in the Palaearctic has resulted in the creation of secondary grasslands dependent on human land use that prevents succession towards shrublands or forests (Török & Dengler, 2018). The coverage of major ecological gradients and the high diversity of vegetation types across several biogeographic regions highlight Palaearctic grasslands as an excellent model system to study small-grain SARs and how they are affected by different factors.
Here, we used more than 2,000 nested-plot series from the GrassPlot database , from a wide range of grassland types across six biomes, to perform a comprehensive analysis of small-grain (1 cm 2 -1,024 m 2 ) SARs in continuous vegetation for vascular plants, bryophytes and lichens. Specifically, we aimed to address the following questions using the Palaearctic grasslands as an example: 1. Which function is most appropriate to describe small-grain SARs?
2. Does the performance of the different functions depend on factors such as sampling method, taxonomic group, biogeographic setting or vegetation type?

| Vegetation-plot data
We used plot data from the collaborative vegetation-plot database GrassPlot  http://b.link/grass plot), which is registered in the Global Index of Vegetation-Plot Databases (GIVD; Dengler et al., 2011) as EU-00-003. GrassPlot collects vegetationplot data (both richness and composition) together with methodological, environmental and structural information from grasslands as well as other plant communities dominated by herbs, dwarf-shrubs or cryptogams from the Palaearctic biogeographic realm (for delimitation see Figure S1.1). Requirements for inclusion are that the plots (sampling units) were precisely delimited in the field and carefully sampled with the aim of achieving complete species lists. One strength of GrassPlot is the numerous multi-scale datasets derived from a diversity of nested-plot sampling schemes (e.g. Dengler et al., 2016) of areas from 1 cm 2 to 1,024 m 2 (schemes of the three main types of sampling designs in Figure S2.1).
We retrieved all nested-plot series contained in GrassPlot (v.1.27 on 4 January 2019) that comprised at least seven different grain sizes (see overview of the 69 datasets with these data in Table S1.4). In total, there were 2,057 series with vascular plant information (Figure 1), of which 757 also contained bryophyte data and 780 lichen data ( Figure   S1.2). The plots were distributed over 26 different countries from 34.9° to 68.9°N, from 9.1°W to 161.8°E and covered an altitudinal gradient from 0 to 4,387 m a.s.l. (Figure 1, Figure S1.2). In total, the nested-plot series consisted of 139,265 individual subplots with richness data, often with several replicates per grain size. Further characteristics of the used datasets are provided in Appendix S1.
For those nested-plots series with more than one subplot for a certain grain size, we averaged richness values across subplots and stored the information on how many subplots the average was based on. Thus, we obtained one single richness value per each grain size within each nested-plot series, if possible, for four different taxonomic groups (1 -complete terricolous macroscopic vegetation; 2 -vascular plants; 3 -terricolous bryophytes; 4 -terricolous lichens). We also recorded whether plots were sampled with the shoot presence or with the rooted presence method (for terminology, see Dengler, 2008).

| SAR modelling
From the numerous different functions proposed for modelling SARs (Dengler, 2009;Tjørve, 2003), we selected three main functions that have specifically been suggested and used for SAR modelling in continuous vegetation (DeMalach et al., 2019;Dengler & Boch, 2008): the power function, the logarithmic function (often erroneously termed the exponential function) and finally the Michaelis-Menten function as a simple two-parameter example of a SAR with saturation (i.e. an upper threshold of richness). To account for possible 'scale dependence' of the SAR, we added two variants of the power function that allow for exponents to change with area: the 'quadratic power function' with a continuous change of the exponent, and the 'breakpoint power function' with an abrupt change of the exponent at a certain grain (e.g. Dengler, 2010). The five functions were selected to represent fundamentally different shapes (Dengler, 2008; see also Figure S3.1) as well as different complexities (number of fitted parameters; Table 1).
We fitted all five functions for both species richness S (S-space; 'linear space') and for log S (log S-space; 'logarithmic space') as dependent variables using nonlinear regression (Table 1). Both approaches are valid, have been used in the literature, and have different strengths and limitations (see Dengler, 2009). Due to the different treatment of the error structure, the parameter estimates in the two spaces usually slightly deviate. Generally, fitting in S-space gives more weight to good fit at large grain sizes, whereas fitting in log S-space gives more weight to good fit at small-grain sizes. Moreover, fitting in log S typically reduces heteroscedasticity in the residuals.
As fitting in log S-space is not possible if some subplots have S = 0 (excluding such cases is not recommended; Dengler, 2010;Williams, 1996), we addressed this issue as follows. Fitting nestedplot series in the optimal case means that the richness value for the smaller grain sizes is representative for the whole area of the largest plot, which could be achieved by full tessellation of its area and averaging the richness values of all resulting subplots. In such optimal sampling, evidently the mean richness value of any smaller F I G U R E 1 Density and spatial distribution of the 2,057 nested-plot series in the Palaearctic biogeographic realm that were analysed in this study. The colour scale indicates the number of available series per 10,000-km 2 grid cell. The map uses the Europe Lambert Conformal Conic projection [Colour figure can be viewed at wileyonlinelibrary.com] The five function types used in this study to model species-area relationships (SARs). All functions were fitted both in S-space and in log S-space. The following notations are used: S = mean species richness; A = area/m 2 ; log = log 10 . The k fitted parameters (except the variance) are termed c, z, z 1 , z 2 , b 0 , b 1 and T

Meaning of parameters
Power powSAR 2 S = c A^z log S = log c + z log A c = richness at unit area (1 m 2 ); z = steepness parameter (exponent in S-space or slope in log S-space) c = richness at unit area (1 m 2 ); z 1 = steepness parameter; z 2 = change of steepness with increasing area c = richness at unit area (1 m 2 ); T = breakpoint (area at which the steepness changes); z 1 = steepness parameter for A > T; z 2 = steepness parameter for A ≥ T adj. values for the model fit optimized using these starting parameter values were then calculated. We consider AICc and R 2 adj. as adequate measures for the relative appropriateness/superiority of the compared SAR functions despite the non-independence of the data points in our nested-plot data. In Appendix S4 (R codes in Appendix S5 and S6), we sampled from virtual landscapes where the shape of the SARs is known to test whether nested-plot sampling introduces biases in the model selection using AICc and R 2 adj . The 'true shape' of the SARs in these virtual landscapes was determined by averaging the results of several random non-nested plot series of different grain sizes. The results show that the model ranking obtained by nested-plot sampling is close to the true pattern and actually depicts, on average, the true pattern better than SARs constructed from a series of single non-nested plots in the same landscape would do. Accordingly, we consider our approach as valid.
In a small number of cases, there were multiple optimized model fits (i.e. with different parameter estimates) with identical (maximum) likelihood values; here, we simply selected one set of parameter values at random. Following standard statistical convention, the variance was always considered as an additional parameter when calculating AICc. Thus, for example, the power model was considered to have three parameters when calculating AICc. For the power breakpoint model, a further model-fitting step was implemented. In certain cases, the best-selected power breakpoint model using the aforementioned approach contained a z-value that was greater than 1 or less than 0. This z-value was then fixed at either 0 or 1 (depending on which of these values it was initially closest to) and the model fitting process repeated. If both original z-values were out of bounds, this additional step was not undertaken. For a given model and plot series, the above model fitting process was repeated across all four taxonomic groups.
For the log S-space analyses, the logarithmic, Michaelis-Menten and breakpoint power functions were fitted using the nonlinear fitting procedure outlined above, whereas the power model and the quadratic power function were fitted using linear regression and the standard 'lm' function in R. The overall model fitting process was relatively computationally demanding and took approximately 48 hr on a 24-core computer cluster (100 GB RAM). Due to the brute-force approach, we achieved convergence of all models for all taxa in all datasets in the log S-space, and a negligible amount of non-convergence in the S-space (maximum 4% for lichens, but 0% for complete vegetation). The R code used to run the analyses is available as Appendix S7.

| Ranking and comparison of the SAR functions
We ranked model performance in five ways. First, we counted for how many nested-plot series a certain function performed best among all compared functions, using model selection based on AICc (Burnham & Anderson, 2002). Second, for each function we calculated the Akaike weights based on AICc in each nested-plot series. Akaike weights can be interpreted as the probability that the function i is the best model for the observed data, given the set of five candidate models (Johnson & Omland, 2004). Third, for each function by nested-plot series combination we calculated Δ i , that is, the difference in AICc of the particular function compared to the respective best performing function ('delta AICc'). Fourth, we ranked models using R 2 adj. , which was calculated using the formula: 1 -(1 -R 2 ) (n -1)/k, where R 2 is the standard R 2 , n is the number of data points and k is the residual degrees of freedom. Fifth, we determined the best performing function based on the Bayesian Information Criterion (BIC) as there is no clear agreement on the superiority of AIC/AICc versus BIC (Burnham & Anderson, 2002;Johnson & Omland, 2004). The five comparisons were undertaken only in cases where our fitting procedure yielded a result for all five models. Note that model comparisons are restricted to each 'space', that is, measures of goodness of fit or information content (e.g. R 2 adj , AICc) cannot be compared between S-space and log S-space (Dengler, 2009).
As sampling methodology has been repeatedly suggested to influence the shape of SARs (Dengler, 2008;Williamson, 2003), we tested for an effect of some key sampling method aspects using  Table S1.1 and Figure S1.4). Furthermore, we tested whether the performance of the functions depended on (d) taxonomic group (vascular plants, bryophytes, lichens), (e) biome (Bruelheide et al., 2019; based on Schultz, 2005), (f) vegetation type or (g) richness in the largest plot of the series (see Figure S1.1, Tables S1.2 and S1.3).

| General suitability of the compared functions
Given the wide range of vegetation types studied, the species-area curves also varied widely ( Figures S8.1 and S8.2). For all taxonomic groups and irrespective of S-space versus log S-space, the power function was by far the best model when using AICc as a model selection criterion (Figure 2). For the richness of the complete veg-  (Figure S8.7). Owing to the one or two additional parameters, the quadratic and breakpoint power function, respectively, had a slightly better fit than the normal power function.
The resulting parameter estimates of all five models and their descriptive statistics are provided in Appendix S9. Here, we summarize only the results of the power function parameter estimates, as it was clearly the overall best model. In particular, we focus on a few parameters that are particularly relevant for interpretation. The slope parameter (z-value) of the overall best performing function (power function) in S-space was 0.20 ± 0.05 (mean ± standard deviation) for all taxa, with slight variation among the three taxonomic groups (vascular plants: 0.26 ± 0.11; bryophytes: 0.19 ± 0.12; lichens: 0.28 ± 0.14). In log S-space, the values showed a similar pattern with little deviation in absolute values from S-space (Table S9.

| Factors influencing function performance
The relative performance of the five models was strongly influ-

| The nature of the species-area relationship
We found strong support for the power function SAR at small-grain sizes in continuous vegetation using more than 2,000 nested-plots over large ecological, geographical and diversity gradients for three major taxa and when focusing on the complete vegetation. Using AICc and R 2 as measures, the 'normal' power function was on average the best model. Using BIC, the breakpoint power function prevailed, and the quadratic power function had a similar level of support to the normal power function. This difference is not astonishing as BIC penalizes complexity of a function differently than AICc, but actually less strongly for small sample sizes, which might lead to overfitting. If basing the conclusions on BIC, there might be some scale dependence of the SAR, that is, a minor change of the exponent z with grain size (see also Crawley & Harral, 2001). If all three variants of the power function are considered jointly, their prevalence as the best model increased from c. 60%-90% based on AICc to c. 90%-95% based on BIC. With our simulation (Appendix S4), we could further demonstrate that this result was not caused by the non-independence of the nested plots, but that this sampling approach, if at all, might even slightly underestimate the superiority of the normal power function.
The general superiority of the power function was largely unaffected by taxonomic group, biome or vegetation type. This finding is in line with previous regional studies analysing small subsets of the current database (Dengler, 2009;Dengler & Boch, 2008;Fridley et al., 2006). Although we restricted our comparisons for pragmatic reasons to a smaller set of functions, which still provides a good representation of the overall range of possible SARs, our findings are consistent with those of Triantis et al. (2012) and Matthews et al. (2016) for true islands and habitat This suggests that, in spite of the commonly accepted notion that contrasting factors influence species diversity at different spatial scales (Brown & Peet, 2003;Field et al., 2009;Shmida & Wilson, 1985;Siefert et al., 2012), the resulting SARs are astonishingly similar over many orders of magnitude (see also Wilson et al., 2012) and across taxa and ecological conditions. Although the power function has been repeatedly criticized (e.g. Pan, Zhang, Wang, & Zhu, 2016;Stiles & Scheiner, 2007), our study supports the idea that it is indeed the one model that, at least for plant communities, can be universally applied (note that even in those cases where it was not the best model, it performed very well; see Figure S8.7). In contrast, other models are suitable, at best, in only a few specific cases.
This poses the question of why a single function (but with varying parameters, see next subsection) can be suitable across so many different situations. In fact, power law SAR-like relationships are far from restricted to species diversity versus area, but can likewise be found in other natural phenomena, such as species frequency versus body size, or body size versus area (Southwood, May, & Sugihara, 2006), or even in completely different realms of science and everyday life (Nekola & Brown, 2007; but see Stumpf & Porter, 2012, for a critical view). A general finding from these different disciplines is that power functions most often result from non-equilibrium conditions (Mitzenmacher, 2012) or skewed underlying distributions (e.g. Rosindell & Cornell, 2007). Power law relationships are likely the consequence of complex dynamical systems, not necessarily of specific ecological mechanistic processes (Nekola & Brown, 2007), even if the slopes of the power law SARs might well be effected by such processes. In this respect, it is interesting to compare the relative performance of the power function across taxonomic groups. Performance was highest for all species groups combined, followed by vascular plants, bryophytes and lichens, which corresponds to the mean species richness of each group.
Moreover, in vascular plants (the groups with the biggest dataset), we found a strong increase in the superiority of the power function with the richness in the biggest plot. It seems that the more elements (here: species) with slightly different properties (e.g. frequencies, habitat preferences, sizes) are involved, the more closely power functions are approached.
In addition, our findings suggest that there is likely no saturation in SARs in continuous vegetation as our saturation function (Michaelis-Menten) performed on average much worse than the functions without saturation across a wide array of different ecological conditions. We believe that even the few individual datasets where the Michaelis-Menten function appeared to be superior are likely artefacts of insufficient replication at the smaller grain sizes.
As Dengler and Boch (2008) have shown, the relative performance of the power function versus saturation functions improves when the replication of smaller subplots is increased and thus the calculated average richness is closer to the true mean richness. This is in line with our finding of best fits for the Michaelis-Menten function for bryophytes and lichens in S-space. As these groups often have few species in grasslands, in many cases none of the smaller subplots (across several grain sizes) contained any species, resulting in a recorded richness of 0, despite the fact that the true average must be higher and increase with grain size (see Methods).
We thus recommend that the concept of 'minimal area' (which only has a meaning if saturation exists), that has been presented in numerous textbooks of vegetation science (e.g. Barbour, Burk, & Pitts, 1999;Kent, 2012;Mueller-Dombois & Ellenberg, 1974) for over a century, should be completely abandoned, as has already F I G U R E 5 Comparison of model performance of the power function expressed as AICc weights between the six biomes represented in the study. The displayed values are for the S-space (results in log S-space were consistent) [Colour figure can be viewed at wileyonlinelibrary. com] clearly been stated by previous studies (e.g. Barkman, 1989;van der Maarel, 1996).
The same holds for the logarithmic function (in the literature also termed 'exponential' or 'semi-log'). There is a widespread belief in vegetation science that this function is particularly appropriate at small-grain sizes (Gleason, 1922;He & Legendre, 1996;van der Maarel, 1997;Stohlgren, Falkner, & Schell, 1995), but the origin of this impression is unclear. For example, van der Maarel (1997) claims this despite the fact that the curvature of the SARs in his figure clearly suggests the better fit of a power law SAR. In addition, the logarithmic function cannot serve as an appropriate model for the SAR as it necessarily predicts negative richness values for small positive areas, often even within the fitted range (Dengler, 2008; see also Figure S3.1). The very poor performance of the logarithmic function across more than 2,000 grassland plots series throughout the Palaearctic matches the findings for multiple habitat types on the Curonian Spit, Russia (Dengler, 2009), and the south-eastern United States (Fridley, Peet, Wentworth, & White, 2005).
While we could rule out the logarithmic function and saturated functions as suitable models, at closer inspection, we found very small but consistent deviations from a power function with a uniform z-value. Increasing the number of replicates strongly increased the relative support of the power function variants with a varying exponent. Moreover, we found that the z 2 -value of the quadratic power function was significantly negative across all studied nested-plot series (e.g. -0.017 for all taxa), meaning that the actual slope is slightly decreasing towards larger grain sizes. While this pattern was to be expected for rooted sampling (Williamson, 2003), we found it also for shoot sampling (not shown), which would support the steep-flat-steep triphasic theory of SARs by Storch, Keil, and Jetz (2012). Our results show that such minor deviations could conveniently be accounted for in the power model by allowing z-values to change with grain size in a systematic manner (with the quadratic or the breakpoint variants of the power function).

| Methodological aspects
A few studies have found a much better performance of saturated and/or logarithmic functions compared to power functions at small spatial scales (DeMalach et al., 2019;Stiles & Scheiner, 2007). However, these authors analysed species accumulation curves and species-sampling relationships (SSRs; Dengler, 2009) rather than SARs in the strict sense (see the typology of Dengler, 2009), and thus these findings are not surprising. Even though their SSRs were also based on 'areas' (and thus many researchers continue calling them SARs in agreement with Scheiner, 2003), they have fundamentally different mathematical properties (Dengler, 2009). We illustrate this with our conceptual Figure 6 and Table 2. SSRs (whether based on individuals, samples or areas) increase sampling intensity within the same pre-defined focal area, while SARs in the strict sense actually increase the focal area. SSRs thus must be a saturation function by definition, as also shown by the simulations of Dengler & Oldeland (2010).
The fact that in such situations the logarithmic function also performs well (or even better than the rather inflexible saturated Michaelis-Menten function) as well as (or even better than) the rather inflexible saturated Michaelis-Menten function) has to do with the similar shapes of the two functions -at smaller grain sizes both become steeper or, in other words, shows an increasing negative deviation from the power function ( Figure S3.1).
Irrespective of whether area-based SSRs are called SARs or not, results (model superiority and parameter estimates) from analysing this type of curves are not comparable with those of SARs in the strict sense (whether these are nested-plot SARs in continuous habitats or island SARs).
Uncertain richness estimates, particularly underestimations might also mask fits of the power function and increase the relative performance of other models. For instance, Guilhaumon et al. (2008) reported relatively poor performance of power functions and large uncertainties in predictions of global hotspot species richness due to low or uncertain sample coverage. This coincides with our finding that the superiority of the power function was lowest for bryophytes and lichens, the two taxa with the lowest richness in most cases, because low absolute richness means that even a recording error of one species can be a substantial relative error.
Likewise, the superiority of the power function increased when the mean richness values at small-grain sizes were based on averages and thus more reliable than when they were based on single counts. This was also found by Dengler & Boch (2008), who argued that adding random noise to the true relationship by chance will lead to higher superiority of other functions in some cases. We found that other methodological aspects can have pronounced effects on model superiority even when focusing on SARs in the strict sense.
Specifically, we found that the power function performed much better for shoot presence sampling than for rooted presence sampling, which is in line with the predictions of Williamson (2003) and Dengler (2008). Theoretically, both of the widely applied ways to record plants in plots must theoretically lead to deviations from the shape of a perfect power function towards the smallest grain sizes, with z-values of the shoot presence method approaching 0 and those of the rooted presence method approaching 1. However, the deviation from a relatively constant z-value at larger grain sizes should appear at relatively larger grain sizes for rooted sampling and be more pronounced (Dengler, 2008), which evidently causes the much lower relative performance of the power function at the small-grain sizes studied here, in the case of rooted sampling in otherwise similar communities.

| Conclusions and outlook
While a perfect power function theoretically cannot hold across all grain sizes (Storch et al., 2012) from, for instance, 1 cm 2 to the terrestrial surface of the Earth (130 million km 2 , or 18 orders of magnitude), we found that it is a very good approximation for sessile organisms across the already large range of six orders of F I G U R E 6 Overview of the main types of species richness curves in terms of spatial arrangement of sampling units and their combination to 'areas' as well as the resulting function shapes. It is evident that species-area relationships (SARs) in the strict sense differ fundamentally from area-based species-sampling relationships (SSRs = species accumulation curves). The species richness curves are assigned to the typologies of Dengler (2009, bold) and Scheiner (2003, normal font) [Colour figure can be viewed at wileyonlinelibrary.com] magnitude in our study -despite the very wide ecological and floristical gradients included (e.g. 6 of the 10 global biomes, 18 major vegetation types). This is in line with the findings of the equally comprehensive studies of Triantis et al. (2012) and Matthews et al. (2016), who found a similar superiority of the power function across many orders of magnitude for multiple taxa in true and habitat islands at much larger grain sizes than in our study, but equally across many orders of magnitude. The superiority of the power function has also been shown at similar grain sizes and in continuous vegetation as well as habitats other than grasslands (e.g. forests and wetlands) (Dengler, 2009;Fridley et al., 2005). This leads us to conclude that the power function is a suitable (and mostly the best possible) model for SARs in nearly any situation, provided the areas from which the relationship is constructed are contiguous. For curves constructed from virtual areas consisting of non-contiguous subunits (as in the case of area-based species accumulation curves), a saturated function, rather than a power function, is to be expected (Dengler & Oldeland, 2010). As a consequence, power functions are usually not suitable for estimating species loss due to habitat loss, as the remaining habitat is typically highly fragmented (Hanski, Zurita, Bellocq, & Rybicki, 2013). However, in all cases with contiguous areas, be it islands (of any type) or areas in continuous habitats delimited by the researcher, power function SARs are suitable tools for interpolation and extrapolation of species richness, or for removing the area effect if other drivers of biodiversity are the focus. Moreover, power function SARs provide, with their exponent (z-value), a meaningful (and standardized) beta-diversity measure in continuous vegetation (Jurasinski et al., 2009;Polyakova et al., 2016), enabling the effective comparison of species turnover among taxa or between different ecological conditions.

ACK N OWLED G EM ENTS
We thank all vegetation scientists who carefully collected multi-

DATA AVA I L A B I L I T Y S TAT E M E N T
The GrassPlot data used in this study are collectively owned by the GrassPlot Consortium, consisting of several hundreds of contributors. The GrassPlot database  http://b.link/grass plot) is continuously growing, while specific older versions used for published analyses, such as version 1.27 analysed in this study, are permanently stored. They can be requested from the first author conditional on agreement with the GrassPlot Bylaws (https ://bit.ly/2IX2Svu), which ensure a fair involvement of data originators in emerging products. The data that support the findings of this study will become openly available in Dryad at http://doi.org/doi:10.5061/dryad.4j75656 after a 3-yr embargo period. The R code used in this study is provided in Appendix S5-S7. TA B L E 2 Major differences between species-area relationships (SARs) in the strict sense and species-sampling relationships (species-sampling relationships = species accumulation curves)