Species richness and beta diversity patterns of multiple taxa along an elevational gradient in pastured grasslands in the European Alps

To understand how diversity is distributed in space is a fundamental aim for optimizing future species and community conservation. We examined in parallel species richness and beta diversity components of nine taxonomic groups along a finite space, represented by pastured grasslands along an elevational gradient. Beta diversity, which is assumed to bridge local alpha diversity to regional gamma diversity was partitioned into the two components turnover and nestedness and analyzed at two levels: from the lowest elevation to all other elevations, and between neighboring elevations. Species richness of vascular plants, butterflies, beetles, spiders and earthworms showed a hump-shaped relationship with increasing elevation, while it decreased linearly for grasshoppers and ants, but increased for lichens and bryophytes. For most of the groups, turnover increased with increasing elevational distance along the gradient while nestedness decreased. With regard to step-wise beta diversity, rates of turnover or nestedness did not change notably between neighboring steps for the majority of groups. Our results support the assumption that species communities occupying the same habitat significantly change along elevation, however transition seems to happen continuously and is not detectable between neighboring steps. Our findings, rather than delineating levels of major diversity losses, indicate that conservation actions targeting at a preventive protection for species and their environment in mountainous regions require the consideration of entire spatial settings.


Methods
Along a SW exposed, complete elevational gradient spanning from the lower montane zone (valley bottom) to the alpine zone (highest peak), we selected three replicate sites at four elevations each (1,000 m, 1,500, 2,000, 2,500, n total = 12, Fig. 1, site coordinates in Supplementary Table S1) within a horizontal distance of 2 km beeline.
All 12 sites have comparable slope (5-15°), are non-intensively grazed by cattle (0.5-1.5 livestock units per ha), are without additional fertilization or irrigation and were not subjected to substantial land-use-changes over the last 160 years 34 . The sites at the two lower elevation steps belong to an area which is protected by the Habitats Directive of the EU (code 6,240*) 35 . field sampling. We recorded species presence and where possible abundance for nine taxonomic groups with standardized methods. For animal groups, we carefully met basic sampling conditions [cf. 36 ] and always looked for exhaustiveness of species recording. A central point was established for each site wherefrom each group was sampled according to taxon-specific methods.
Scientific RepoRtS | (2020) 10:12516 | https://doi.org/10.1038/s41598-020-69569-9 www.nature.com/scientificreports/ • Vascular plants were recorded from five random subplots (3 × 3 m) according to Elzinga et al. 37 within an area of 531 m 2 (radius from central point = 13 m). Species occurrence was recorded for each subplot. • Bryophytes and lichens were recorded from five random subplots (50 × 40 cm) according to Elzinga et al. 37 within an area of 531 m 2 . Species occurrence was recorded for each subplot. Specimens of species which were difficult to identify in the field were collected and identified in the laboratory. • Grasshoppers including crickets and locusts were surveyed three times during summer period to ensure that all individuals were adult. They were caught by hand and with a sweep net for 20 min within 500 m 238,39 , identified and released. Additionally, vocals were used to identify species. • Butterflies and burnets were caught three times during summer period with a sweep net (45 min per site) within 500 m 238,39 , identified and released. • Earthworms were heat extracted from three soil core samples (20 × 20 cm, 15 cm deep, taken randomly within a 100 m 2 area) with a modified Kempson apparatus 40 according to Steinwandter et al. 33 . After extraction animals were identified to species level using a dissecting microscope. • Beetles were taken from soil core samples (cf. earthworms) and caught with two pitfall traps (plastic cups with 8.5 cm diameter filled with 150 ml of propylene glycol) installed on opposing sides within a 5 m distance from the central point for two weeks [cf. 41 ]. Individuals from pitfall traps and soil core samples were transferred to 75% ethanol and only adult specimens were identified to species level using a dissecting microscope. • Spiders were collected by visual search, within soil core samples and pitfall traps (cf. beetles 42 ). Individuals were transferred to 75% ethanol and only adult specimens were identified to species level using a dissecting microscope. • Ants were recorded by counting nests within 4 m 2 subplots and additionally collected within pitfall traps (cf. beetles) according to Schlick-Steiner et al. 43 . Individuals were collected and identified using a dissecting microscope.

Data analyses.
We first verified our sampling completeness by estimating sampling coverage (Table 1) with the function iNEXT within the R package inext 44 . Datatype was set to "incidence_freq" for bryophytes, lichen and vascular plants (due to data structure) and to "abundance" for all other taxa. In the case of incidence data, sample size refers to the number of sampling units, whereas for abundance data sample size equals the number of individuals in a sample 45 . Only a consistent and high sampling coverage estimation allows a compositional comparison of different taxonomic groups 46 . Species richness was calculated by computing the cumulative number of species per site (n = 12) for each taxon separately. By applying the function spline.plot within the R package drsmooth 47 , we plotted a splineestimated dose-response function on the actual richness data along the elevational gradient with its upper and lower 95 percent confidence bounds.
In order to have the same data base for all taxa, beta diversity was calculated using species occurrence data (not including abundance), which was then partitioned into turnover and nestedness components by applying the function beta.pair within the R package betapart 17 . This results in three matrices based on pair-wise comparisons of each site: the Sørensen dissimilarity index (β sor ) expresses the total compositional variation with values ranging between 0 and 1, the Simpson dissimilarity index matrix (β sim ) compositional changes due to species turnover, and β sor minus β sim is the resultant nestedness component β sne . Study design with visualization of along and step-wise elevation beta diversity approach. Beta diversity components turnover and nestedness were calculated for all possible pairs and further distinguished between (i) along elevation beta diversity accounting for pairs starting from the lowest elevation to all other elevations (1,000-1,500, 1,000-2,000, 1,000-2,500 m a.s.l.), and (ii) step-wise beta diversity including pairs from one elevation to the next neighboring one (1,000-1,500, 1,500-2,000, 2,000-2,500 m). Each elevation was surveyed with three replicate sites (n total = 12). www.nature.com/scientificreports/ We further processed beta diversity data (β sor, β sim, β sne ) to distinguish between (i) along elevation beta diversity ( Fig. 1) accounting for pairs starting from the lowest elevation to all other elevations (1,000-1,500, 1,000-2,000, 1,000-2,500 m), and (ii) step-wise beta diversity (Fig. 1) including pairs from one elevation to the next neighboring one (1,000-1,500, 1,500-2,000, 2,000-2,500 m).
The effect of elevation on both step-wise and along beta diversity was tested with a linear mixed effects model using the R package lme4 48 . As fixed effects, we tested the elevational distance, as random effects we included the elevational replicate sites. We also tested for geographical autocorrelations in our dataset. To do so, we decomposed the total Euclidian distance between two sites into elevational distance (the variable to be tested) and the remaining distance (considered as potential spatial autocorrelation). We calculated this remaining distance by applying the Pythagorean theorem with elevation as one leg of a hypothetical triangle and the total distance as the hypothenuse 49 , and used it as fixed factor in the above explained model. Visual inspection of residual plots did not reveal any obvious deviations from homoscedasticity or normality. Significance of fixed factors was confirmed if Akaike Information Criterion 50 was lowered by at least 2 points (AIC < 2) in comparison to the null model. Further, the R 2 partitioned into marginal and conditional R 2 was computed following Nakagawa and Schielzeth 51 . All analyses were conducted, and corresponding figures were produced using the open-source statistical programming language R (version 3.6.2, R Core Team 52 in R Studio, version 1.1.383, RStudio Team 53 ).  46 ] was very high for all elevations and taxa with a mean SC value of 0.95 (SD = 0.07). In few cases and only for singular elevations, sample-size-based R/E sampling curves, which compute diversity estimates for rarefied and extrapolated samples and plot the diversity estimates with respect to sample size 44,45 , did not flatten completely.

Species richness.
A hump-shaped relationship between elevation and species richness was detected for vascular plants, butterflies, beetles, spiders and earthworms, while species richness decreased linearly with elevation for grasshoppers and ants (Fig. 2). Species richness of lichens and bryophytes increased with increasing elevation. After this result, earthworms were excluded from all further analyses, since aside of very low species richness (gamma diversity equaled to only 4 species), on 1,000 and 2,500 m very few individuals were found. Beta diversity. Considering all possible pairs of comparisons of all sampled sites, the mean total beta diversity varied from 0.61 (grasshoppers) to 0.79 (lichens and ants). Partitioning total beta diversity and calculating portions, we found the highest percentage of turnover for vascular plants and ants (93.4 and 90.3, respectively), while nestedness was found to be highest for lichens (31.9%) and bryophytes (28.3%) for the investigated gradient.
Along elevational turnover (all comparisons starting from the lowest elevation) increased significantly for all groups (increase of β sim between 0.09 and 0.23), except for lichens (Table 1, Fig. 3).
Regarding step-wise beta diversity only bryophytes and lichens had a significantly lower turnover and overall beta diversity for higher elevated steps, while all other groups showed no significant differences (Table 1, Fig. 5). Nestedness increased significantly with higher elevated steps with respect to grasshoppers and lichens, while results on vascular plants are opposed, featuring a significant lower nestedness at higher elevated steps (Table 1, Fig. 6). The other groups showed no significant differences. Out of the eight analyzed taxonomic groups, ants were the only ones showing a geographical autocorrelation for step-wise pair comparisons.

Discussion.
Our results provide a comprehensive overview of how diversity is organized in a vertical space for eight taxonomic groups (earthworms were excluded after first data analyses) comprising animals, plants, bryophytes, and lichens sharing the same habitat, namely pastured grasslands.
For the majority of the analyzed groups we detected common patterns, when looking at multi-taxonomic species richness and along elevational beta diversity. On the contrary, we found dissimilarities among taxonomic groups in the step-wise beta diversity.
With regard to species richness, it has been shown, that scale of spatial extent strongly impacts the shape of relationship between elevation and species richness: analyses of complete gradients most often lead to humpshaped richness patterns, while omission of parts of the gradient at the upper and/or lower ends tends to favor monotonic increase or decrease with elevation 54 . According to that, with our study we captured the whole gradient for five out of our nine taxa, and the documented hump-shaped pattern is in line with previous work [e.g. 24,25,55 ]. Apparently, the upper limit of the gradient may be missing for lichens and bryophytes and the lower limit for grasshoppers and ants (Fig. 2). Some studies on bryophytes performed in warmer climate zones found a hump-shaped relation of richness to elevation 56,57 . However, in the temperate zone, an increase of richness with elevation for lichens and bryophytes was found also by Nascimbene and Marini 27 and Spitale 28 who explained their results with a negative relationship of richness with temperature and associated conditions such as slower evaporation rates. Regarding grasshoppers and ants, the decrease of species richness with elevation was detected also by Peters et al. 29 , who attributed positive effects of available land areas on grasshoppers' richness. Overall, Peters et al. 29 argue that such unimodal species richness patterns hold only for single taxa analyses [cf. 58  www.nature.com/scientificreports/ elevation. We cannot confirm this finding derived from a study on Mt. Kilimanjaro, since our joint analysis of nine taxonomic groups still shows a hump-shape pattern (Fig. 2, in line with Viterbi et al. 24 ). However, our elevational gradient is influenced by grazing and located in a completely different biogeographic region with other (temperate) climatic conditions, which may impact the shape of this relationship. We are also aware that particularly for animal groups, seasonality of sampling year or period may influence species recordings and hence total species pools of the study site might be larger. For our aim, which is not a species inventory of the sites, a snapshot embracing one season is still meaningful. Beside species richness, which is one of the most frequently used indicators for biodiversity assessments for example when determining locally delimited spaces deserving protection, beta diversity is the key to scale up to the regional extent. The comparison of compositional differences between sites, further partitioned into turnover and nestedness according to Baselga 13 , allows the understanding of how biodiversity is assembled across elevations. As already found for some insect taxa [cf. 18 ], turnover is the dominant component of beta diversity for all our studied taxa, which are vascular plants, lichen, mosses, ants, beetles, spiders, butterflies, and grasshoppers. According to our hypothesis that turnover will increase with increasing elevational distance, due to a diversification of communities, along elevational turnover (and total beta diversity i.e. ß Sorensen) increased significantly with increasing elevational distance for seven out of eight analyzed taxa. These results confirm the findings of studies focusing on single insect groups (e.g. 30 for ants, 20 for dung-beetles, 59 for bees and wasps) and support the assumption that a considerable portion of species is restricted to certain zones and does not colonize entire gradients 30 . For along elevation nestedness, we assumed a decrease with elevation, since communities usually adapt to given conditions and probabilities for subsets between different elevations are lower. Five out of 8 analyzed support our hypothesis, consequently suggesting a differentiation of communities at some point(s), with low occurrences of subsets, as found for dung-beetles 20,55 . Our results underpin the assumption that mechanisms such as environmental filtering and dispersal limitations (reflected by high turnover rates) seem to prevail colonization and extinction patterns in community assembly of pastured grassland along temperate elevational gradients [cf. 59 ].
When looking at similarities of neighboring steps (step-wise diversity), we could not confirm our second hypothesis assuming a decrease of total beta diversity and turnover when comparing sample units of higher www.nature.com/scientificreports/  Step-wise turnover (pair-wise comparisons) of neighboring steps. 1,015… comparison of steps 1,000 and 1,500, 1,520… comparison of steps 1,500 and 2,000, 2,025… comparison of steps 2,000 and 2,500. Figure was produced using R software (version 3.6.2, R Core Team, https ://www.r-proje ct.org/).

Figure 6.
Step-wise nestedness (pair-wise comparisons) of neighboring steps. 1,015… comparison of steps 1,000 and 1,500, 1,520… comparison of steps 1,500 and 2,000, 2,025… comparison of steps 2,000 and 2,500. Figure was produced using R software (version 3.6.2, R Core Team, https ://www.r-proje ct.org/). www.nature.com/scientificreports/ located steps (due to abiotic constraints at mountain-tops promoting community speciation). The 8 analyzed groups do not show a congruent directional pattern -so far only two groups significantly decreased turnover (lichens and bryophytes) or increased nestedness (lichen and grasshopper) on higher elevated neighboring steps. The few other studies which analyzed beta diversity according to Baselga 13 among elevational belts came to inconsistent findings. Paknia and Sh 31 did not find a uniform pattern of beta diversity for moths along elevation, while da Silva et al. 20 found a clustering of dung beetles into lowland (200-800 m) and highland (1,000-1,300 m) communities, driven by the turnover component. For vascular plants turnover was shown to be uneven along elevation and to peak between 1,800 and 2,200 m 60 . We assume that niche-breadth variation of some species 61 , barrier effects of treeline or ecotones 60 , or mobility of non-sessile groups like butterflies, beetles but also grasshoppers might influence turnover of species among neighboring elevational steps. Possibly, the range of our steps (500 m) is too large for detection of specifically adapted and hence similar communities at the upper limits of our gradient (2,000-2,500 m). Summarizing along and step-wise diversity results, we conclude there is a significant change of species communities along elevation; however this transition seems not to be detectable between neighboring steps. This result could either be a hint for a rather slow transition of communities across the analyzed gradient, or for the existence of an elevational threshold that separates communities, which we were not able to identify with our approach (such as 59 ).
conclusion. When unravelling spatial patterns of beta diversity, we found several studies focusing on arthropods in tropical to subtropical latitudes, while our study disentangles to our knowledge for the first time diversity of organisms belonging to different systematical classes in a temperate zone. The manifold possibilities to compare pairs of dissimilarities (e.g. within, along, step-wise, or complete gradient) and the frequent lack of specification of compared pairs present in literature, significantly hamper comparisons across studies and collective findings. However, our results confirm turnover to be the dominant component of beta diversity for all investigated taxa along the elevational gradient. Further, we detected a change of species communities for seven out of 8 analyzed taxa with increasing elevational distance but were not able to identify thresholds or delimitate zonal limits or boundaries of communities along the gradient. From a conservational point of view, this knowledge rather than delineating levels of major diversity losses, supports the idea that action sets or monitoring programs targeting at a capacious and preventive protection for species and their environment in mountainous regions require the consideration of entire gradients. Including information about specialist, rare or endemic species, niche breadth, or species traits into beta diversity assessments would significantly increase the future uncovering of spatial patterns of community assemblage.

Data availability
All species recorded during this study are included in Supplementary Information files of this article.