Species elevational richness gradient and species-area relationship in mountain vegetation of Javakheti highland (Georgia)

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


Introduction
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 lo-cation 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.

Study area
The Javakheti region is located in the southern part of Georgia, on the Lesser Caucasus Mountains (Fig. 1).It is a highland of volcanic origin with a total area of 6421 km 2 within the borders of Georgia.Javakheti Highland is represented by lower valleys (1400-1700 m), a highland plateau (1700-2200 m), and two longitudinally arranged mountain ranges -Samsari and Javakheti, with a maximum absolute elevation of 3300m a.s.l.The region is very rich with lentic and lotic water bodies and supports the largest mountain wetland ecosystems in the Caucasus (Maruashvili 1964).
The climate of Javakheti is dry and continental.In the altitudes between 1700 and 2500 m a.s.l., the mean annual temperature (mean: 7 C, min: 17.8 C, max: 20.4 C) and annual averaged precipitation (550mm) are low in comparison with other mountain areas of the Caucasus at the same elevations.Annual total precipitation increases towards higher elevations, up to 1400 mm, and the mean temperature decreases to roughly 1.4 C per 100m elevation (Zanina 1961;Maruashvili 1964).
Javakheti highland represents a transitional province between relic Kolchic, western Asian, and partly Holarctic and Mediterranean flora (Nakhutsrishvili 2013).A great deal of the highlands is represented by mountain meadows and wetlands, while the natural forest cover is extremely small (only a few hectares) and fragmented.Forest ecosystems on Javakheti highland are thought to have been much more widespread during the Holocene, but before the last 2-3 millennia, forests disappeared as a result of human impact (Arabuli et al. 2008;Messager et al. 2013;Connor and Kvavadze 2014).Wetland vegetation is widely spread on the Javakheti plateau.Apart from the lake shores and river banks, they appear in numerous mountain bogs represented in the region (Zedelmaier 1929(Zedelmaier , 1933;;Nakhutsrishvili 1966;Tedoradze et al. 2023).

Data collection
Field data collection was conducted during the three-year (2014-2017) period during the whole vegetation season.Field sampling was focused only on the natural grassland plant communities, including mountain steppes and wetlands.In total, 170 plots (Fig. 1) with a size of 100 m 2 (10 m x 10 m) were investigated using maximum variation sampling, according to Meier and Hofer (2016).Samples were taken haphazardly in the studied region, with a minimum distance between plots of 1 km.Plant species were identified using regional identification keys either in the field or in the laboratory after the preservation of specimens (Dimitrieva 1959; Ketskhoveli et al. 1971Ketskhoveli et al. -2013;;Czerepanov 1994;Gagnidze 2005;Fischer et al. 2018).

Data analyses
To obtain the planimetric surface area of the study region, we used a 30 m resolution Shuttle Radar Topography Mission Digital Elevation Model and software ArcGIS v.9.3 (ESRI, Redlands, CA).The total area was then split into 100-meter elevational bands (equal elevational bands -EEB), and the areas were calculated for each of them (Suppl.material 1, table S1).We also calculated three dimensional surface area with the slope method using the approach of Berry (2002) and Jenness (2013).Since the correlation between the planimetric and three-dimensional areas was very high (Pearson's r > 0.9) for all grain sizes (see below), only the planimetric area was used in the follow up analyses.
The general elevational trend of plant species diversity was examined using first-and second-order simple regression techniques.Analyses were applied to the log-transformed raw species richness.The same regression techniques were also used for rarefied, and estimated species richness (Chao 1984;Gotelli and Colwell 2001;Chao et al. 2016) to reduce the bias sourced from unequal sampling.We also examined the residual (derived after regressing raw richness vs. SE) distribution along with elevations to see the direct effect of sampling on the perception of EGSD.These procedures were repeated for equal elevational bands (EEB) of different grain sizes (100, 200, 300, and 400 m) (Suppl.material 1, table S1).
The methodology used by Mumladze et al. (2017b) was followed to explicitly disentangle the effect of altitude, area, and sampling effort on observed species richness.In particular, when no complete diversity data is available for a study region, sampling may play a pivotal role in determining the richness pattern.To see the effect of sampling effort, we first divided the total elevational profile into equal sampling bands (ESB).We set up 10 consecutive bands of varying area, each with 17 samples, and recalculated species richness for each band, band mean elevation, and area (Suppl.material 1, table S2).In this way, we directly filtered out the effect of SE and examined the effects of altitude and area using GLM (Mumladze et al. 2017b).However, this approach only partly reflects the influence of area on species diversity as the area per-se effect (Rosenzweig 1995;Schoereder et al. 2004) is a priori excluded.In this regard, we also employed the equal area approach, in which case we split a total elevational gradient into equal area bands (EAB) and recalculated species richness (Suppl.material 1, table S3).In this way, we directly removed the effect of area (Bachman et al. 2004) and examined the changes in the response of the species richness curve compared to the EEB approach.Here too, unequal SE may also hinder the inference, so we also modeled rarefied, estimated species richness, and residuals as described above.
R software (R Core Team 2022) was used to fit the regression models.In the case of GLM, we test the Poisson error distribution for overdispersion (Quinn and Keough 2002;Kleiber and Zeileis 2008;Zuur et al. 2009).As the overdispersion was always strong (around 15), we used negative binomial GLMs to correct the variance (Quinn and Keough, 2002;Zuur et al. 2009).Competing models were selected based on a definite sample-corrected AICc and model weights using a multimodel selection framework (Burnham and Anderson 2002;Barton 2016).

Results
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 Supplementary Material 1, Table S4.Average plot species richness was 24 with a declining pattern towards higher elevations (R 2 =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).

Discussion
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 2 =0.65, p<0.001); richness vs. elevation (second order): R 2 =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.

Figure 1 .
Figure 1.A map of the sampling area (Javakheti highland) showing 170 sampling localities.

Table 1 .
Summary table of the best regression models of log-transformed species richness (various measures) vs. elevation and area at different elevational grain sizes.

Table 2 .
Model selection results based on negative binomial regression applied to ESB data (species richness vs. elevation, and area).nnumber of observations; k -model parameters; θ -distribution parameter of the NB function; AICc -value of small sample corrected Akaike information criteria; Wi -model weight.

Table 3 .
Model selection result based on negative binomial regression applied to EAB data (species richness vs. elevation, area, and SE).n -number of observations; k -model parameters; θ -distribution parameter of the NB function; AICc -value of small sample corrected Akaike information criteria; Wi -model weight.