Elevational gradients in species richness and species-area relationships are among the most interesting patterns in ecology and biogeography. Both patterns can be characteristic of the same system; however, current knowledge of how these patterns co-exist and how we can disentangle their contributions to biodiversity structure is insufficient. In this article, we tested the effect of elevation and area on the formation of plant species diversity patterns in the forest-free Javakheti Highlands (Georgia). Samples (170 plots) were collected within elevations of 1400-3100 m, and the diversity distribution was examined in relation to altitude, available band area, and sampling. In total, 564 species from 67 families were recorded. Plant species richness was highest at mid-elevations (1900–2200 m), irrespective of area and sampling effort. This was in line with other studies from the Caucasus indicating the generality of plant elevational diversity patterns in the region. Area was not an important predictor of species richness; however, this may be considered a result of insufficient sampling. Our study shows that more research is needed to understand the effect of area on patterns of elevational biodiversity distribution.

Caucasus, species diversity distribution, plants, species-area relationship

Understanding the patterns and processes of elevational gradients in species diversity (EGSD) is one of the major focuses for ecologists (Lomolino et al. 2006; McCain and Grytnes 2010). Reviews show that the response of species richness to elevation is predominantly hump-shaped (Guo et al. 2013; Rahbek 1995). That is, species diversity approaches its maximum at intermediate altitudes. Although the mid-elevational maximum seems rather general, the variation is also remarkable, displaying the two basic patterns. The first is related to geography, according to which the location of an elevational diversity peak varies greatly (Guo et al. 2013); and the second is that different taxonomic groups exhibit different patterns even in the same geographical areas (e.g. Bhattarai and Vetaas 2003; Lee and Chun 2015; Mumladze et al. 2017a). The possible sources of these variations are supposed to be numerous and related to spatial, climatic, biotic, historical, anthropogenic, and other factors (McCain and Grytnes 2010). These factors are considered potential drivers of EGSD and are routinely studied when the explicit mechanism underlying the pattern needs to be understood. At this end, the necessity of making a distinction between pattern and cause is important. In particular, available land area as a cause of EGSD is infrequently tested (Rahbek 1995; Odland and Birks 1999; Sanders 2002; Bhattarai and Vetaas 2003; Bachman et al. 2004; Kluge et al. 2006; Xu et al. 2017) due to the well-known species-area relationship (SAR) (Arrhenius 1921). On the other hand, SAR is an empirically better supported biodiversity pattern than EGSD itself, and they are not mutually exclusive by definition. Although the same factors might be responsible for both patterns (at least partly), no one can be thought of as a mechanism under the other. Hence, when studying the EGSD, one ought to contrast it with SAR rather than simply search for a causal relationship between these two (Rosenzweig 1995; Jones et al. 2003; Bachman et al. 2004; Mumladze et al. 2017b).

From a methodological point of view, there are two main approaches used in identifying species diversity patterns, along with an elevational gradient. The first approach is based on standardized species distribution data collected on elevational transects. In such cases, sampling effort (SE) is the same at all elevations, and no sampling bias due to unequal sampling exists. A second approach is based on regional data, which mostly comes from herbarium collections and/or literature with usually spatially biased SE and frequently inaccurate spatial metadata. The limitations of both approaches are related to the question we aim to answer. For instance, when studding the EGSD in general, neither approach independently is able to find an optimal solution since transect data could not confidently be extrapolated to a regional scale and vice versa. This is due to the differences in extent and strength of factors affecting species diversity and distributions at different scales (Willig et al. 2003; Chase and Knight 2013). The same is true for species-area relationships (SAR) as well as transect diversity data, which is limited to revealing area effects while most regional diversity data is subject to sampling bias. Although there are some means to statistically correct the biases, this is not possible in most cases until the data are collected using a standardized methodology (i.e. plot surveys). A clear example is Zhang et al. (2015), who, based on herbarium data, showed that specimen records are perfectly correlated with species diversity when studying elevational diversity patterns.

In this article, we report the results of a plot based field survey of mountain grasslands in Javakheti highland (Georgia) in order to evaluate regional plant species diversity in an elevational gradient spanning from 1400-3100 m a.s.l. In particular, by studding the structurally homogenous and continuous mountain grassland ecosystems, we aimed to evaluate the contribution of EGSD and SAR to the overall pattern of grassland community diversity while explicitly considering the SE.

In total, we found 564 vascular plant species belonging to 287 genera and 67 families. The most diverse family was Asteraceae with 97 species (17%), while 13 families were represented with single species in the study area. The full species/site presence/absence data set is represented in Suppl. material 1, table S4. Average plot species richness was 24 with a declining pattern towards higher elevations (R²=0.13, p(F_1,168)<0,001). On average, species frequency of occurrence is relatively low (7.2 (±6.8-SD)). The most widespread species are Trifolium ambiguum and Cephalaria gigantea, with 44 and 41 occurrences, respectively, while 60 (11%) singletons and 71 (13%) doubletons are represented in the data set. The species inventory (calculated at 100 m elevational bends) was quite incomplete according to asymptotic estimation, as on average 111 (±107 sd) species are to be expected in addition to the observed 130 (±94) species for each elevational band.

Regression analyses of log-transformed species richness vs. altitude revealed a hump-shaped pattern of plant diversity irrespective of the grain size or richness measure used (Table 1; Suppl. material 2, figure S1). However, with increasing grain size, the absolute altitude of the observed richness peak varied from 1800–2000 m at the finest grain size to 2200–2300 m at coarser grain sizes. Only at a grain size of 100 m, the response of rarefied species richness to altitude decline linearly. However, this is due to the absence of data points from lower elevations (a result of the small sample size; Suppl. material 1, table S1).

Asymptotic estimator and rarefaction did not reveal a significant effect of incomplete sampling on species richness patterns. However, after accounting for the effect of SE, richness declined linearly with increasing elevation, as shown by residual analyses (Suppl. material 2, figure S1). Likewise, species richness (in log-transformed space) showed a strong linear relationship with available area at all grain sizes. At coarser grains, this relationship was getting stronger (Table 1; Suppl. material 2, figure S2).

The ESB approach showed that the effect of SE is negligible when considering the general pattern of species richness. In particular, the SE-standardized elevational diversity peak did not shift at all while higher elevations proved relatively species poorer than lower elevations in a studied gradient (Fig. 2). The distribution of both raw and estimated species richness was best explained by elevation (second-order linear regression) after removing the effect of SE. The area did not expose any explanatory power as expected (Table 2). When diversity was standardized on area (EAB approach), the general shape of the elevational richness curve was not different from that observed in the EEB approach. Indeed, the diversity peak maintained its elevational position irrespective of the richness measure used (Table 3; Suppl. material 2, figure S3). Only the relative diversity at higher elevations tended to be higher than at lower elevations (comparing the left and right tails of the response curve in Supplementary Material 2, figure S3).

Although the SAR is one of the best-supported diversity patterns, its relevance to elevational species richness distribution is questionable and hardly tested. While some studies have shown that the area can play an important role in shaping the elevation pattern of plant species richness, others do not support it (see Romdal and Grytnes 2007 for a review). However, the application of elevational SAR (eSAR) is still poorly developed (but see Bachman et al. 2004; Xu et al. 2017). Indeed, most of the data sets represent just the species richness-sampling effort relationship rather than actual SAR (Dengler 2009). On the other hand, incomplete sampling can strongly affect pattern perception and obscure the true pattern (Zhang et al. 2015; Yamaura et al. 2016; Mumladze et al. 2017b). The complete diversity data for mainland elevational bands (at the regional scale) is hardly available, which may be a reason why eSAR is so controversial. Indeed, although the elevational band area has been found to be strongly correlated to mite species richness in Georgia (even better than elevation), this was in fact a false interpretation until the SE was considered (Mumladze et al. 2017b). In contrast, we did not find a strong effect of sampling on species elevational richness curves, as its removal did not result in changing the observed pattern. So did the available area, as its effect was negligible beyond the sampling and elevation. On the other hand, standardization of sampling means that we are artificially discarding area per-se effects (Hill et al. 1994; Rosenzweig 1995). Indeed, if there are a large number of rare and patchily distributed species, then area proportional sampling is needed to properly account for the full effect of area (Mumladze et al. 2017b; Schoereder et al. 2004). In our data, 24% of species are represented by singles and doubletons, suggesting a high rate of rarity. Hence, the observed weak effect of area may be a result of insufficient sampling rather than the absence of its effect. Here, we did not manually rarify the species richness proportional to area (as has been done by Mumladze et al. (2017b)) since SE was in strong correlation with the area of EEBs (0.68<r<0.87, p<0.05). In this case, area alone was a better linear predictor of diversity than altitude (richness vs. area (first order): R²=0.65, p<0.001); richness vs. elevation (second order): R²=0.5, p<0.001). However, since the species inventory for each elevational band (at all grain sizes) was quite incomplete, we could not confidently consider area as an important factor beyond the poor statistical sampling. Clearly, more data is needed for an unambiguous assessment of the area effect.

Recent studies of plant elevational diversity in the Caucasus region have demonstrated that plant species richness is highest near the tree line at 2100–2400 m (Erschbamer et al. 2010; Mumladze et al. 2017a,c). As the tree line represents the forest-meadow transition zone, this was thought to be the cause of the highest plant diversity at respective elevations (Mumladze et al. 2017c). That is, Caucasian deciduous temperate forests are comparatively species poor at a site level (compared to subalpine meadows (Dolukhanov 2010); see also case studies given in Mumladze et al. (2017a,c)), resulting in a rapid increase in species richness per area unit near the tree line and then a linear decline with further increases in altitude. On the other hand, all the studies referenced above employed transects or local-scale biodiversity data instead of regional inventories, and thus the observed richness pattern of plant species seems trivial. However, could we extrapolate these results on a regional scale? In our study, we have shown that Javakheti highland, which lacks the high mountain forests (1700–2200 m), is still characterized by the highest species richness at 1900–2200 m (a potential treeline, given the general treeline pattern in the Caucasus), and this pattern is independent of sampling. Although the available area at these elevations is the largest (Suppl. material 1, table S1), the species richness pattern is robust without the effect of area at a given sampling intensity. Presumably, other factors related to elevation (for instance, temperature, precipitation), land use, etc. rather than factors related to area (for instance, habitat diversity) may be responsible for the plant species richness in the Javakheti highland meadows.

The unimodal EGSD of plants is well pronounced in Javakheti Highland, independent of the effect of available area and SE. Species richness peaks at 1900–2200 m in the absence of treeline ecotones. The observed pattern fits perfectly with other cases from the Caucasus region (e.g. Nakhutsrishvili and Batsatsashvili 2017) and suggests a climate-dependent response of plant diversity along an elevation gradient instead of dependence on area-related factors. While eSAR is also observable in the Javakheti Highlands, it seems to be a result of pure statistical sampling rather than the true area effect. Although we cannot completely reject the effect of area on elevational plant species richness at a local scale, apparently eSAR at a regional scale is irrelevant (at least in the Javakheti highlands) as the available area is largest below 1800 m.

In summary, the unimodal EGSD of plants is well pronounced in Javakheti Highland, independent of the effect of available area and SE. Species richness peaks at 1900–2200 m in the absence of treeline ecotones. The observed pattern fits perfectly with other cases from the Caucasus region (e.g. Nakhutsrishvili and Batsatsashvili 2017) and suggests a climate-dependent response of plant diversity along an elevation gradient instead of dependence on area-related factors. While eSAR is also observable in the Javakheti Highlands, it seems to be a result of pure statistical sampling rather than the true area effect. Although we cannot completely reject the effect of area on elevational plant species richness at a local scale, apparently eSAR at a regional scale is irrelevant (at least in the Javakheti highlands) as the available area is largest below 1800 m.

This work is dedicated to the memory of our colleague and friend, the first author – Shamil Shetekauri, who passed away while we were trying to complete the manuscript. Field work was supported by the Shota Rustaveli National Science Foundation under the research grant "Diversity of the Flora of Samtskhe-Javakheti Region (Lesser Caucasus)" (FR/37/7-120/13).