Species richness captures plant functional and phylogenetic diversity variations along different ecosystems on the hillsides of Damavand Mountain (Iran)

The phylogenetic diversity (PD) can be used as a proxy for measures of functional diversity (FD); this relationship between PD and FD is premised on the reasonable assumption that evolutionary diversification has generated trait diversification, which in turn may result in greater niche complementarity. Functional traits have been used in various applications, from assemblage ecology looking at the relationship between functional diversity, environmental filtering, and assemblage structure. We used a trait-based approach, quantifying variations in the leaf-height-seed (LHS) scheme and functional and phylogenetic indices in semi-arid communities from three different habitats (grasslands, shrubland, and mixed habitats). To do so, three functional traits including specific leaf area (SLA), vegetation height (VH), and seed mass (SM) were measured in three habitats. We calculated at the plot scale, community weighted means (CWM) for each trait, species richness, Faith’s phylogenetic diversity index (PD), functional richness (Frich), functional and phylogenetic mean pairwise dissimilarities (MPD), and nutrient concentration. The results showed that the first Soil_PCA axis (explaining 37.1% of the total variance) showed a high loading of soil nutrients (N, P, K, OM) and pH, while the second axis (which explained 22.1% 224 of the variance) exhibited a high loading of pH at the top and K at the bottom. Also, the higher soil nutrient concentration and pH were significant and positively related with species and functional richness and Faith’s phylogenetic diversity.


Introduction
The importance of soil and vegetation management to achieve the United Nations Goals for sustainability (Keesstra et al. 2016) has been highlighted in the last decade due to the central role that both elements play in ecosystem services. Plant-soil interactions provide mechanistic information on the ecosystem's adaptation to the environment (Berdugo et al. 2017), especially in dryland and montane ecosystems, where the soil resources and land use are key factors that explain vegetation structure evolution (Eldridge et al. 2016;Salgado-Luarte et al. 2019;Vaieretti et al. 2018). In this regard, the functional trait-based approach has emerged as a powerful tool to evaluate the evolutionary history and adaptations of plants to a wide range of abiotic factors (Lavorel and Garnier 2002;Jafarian et al. 2019). In addition, many studies have provided information on the potential impact of the global change on vegetation using functional trait analysis of plant communities along environmental gradients across different biomes (e.g., Pérez-Ramos et al. (2012) and Jung et al. (2014)) or at global scale (e.g., Yang et al. (2015) and Wieczynski et al. (2019)). Two assemblages could have identical species richness but could differ dramatically with respect to functional and/or phylogenetic diversity-facets that, respectively, address the range of ecological functions and unique evolutionary histories of assemblages. This multi-faceted approach to quantifying biodiversity has become increasingly common in biodiversity conservation research due to the identification of links between functional traits and ecosystem functioning. Phylogenetic and functional diversities and their relationship are important for understanding community assembly, which relates to forest sustainability. Thus, both diversities have been used in ecological studies evaluating community responses to environmental changes.
However, notwithstanding the relevant information obtained by these studies, additional empirical information on functional trait variation in plant communities is needed to accurately infer their sensitivity to ongoing global change (Pérez-Ramos et al. 2017).
Arid and semi-arid regions will be among the areas on the planet most affected by global change (Delgado-Baquerizo et al. 2013). For instance, recent studies from Iran have shown that global change caused a shift in vegetation composition and soil biogeochemical cycles (Ashrafzadeh and Erfanzadeh 2014). But the impact from a functional and phylogenetic point of view has been largely ignored (Jafarian et al. 2019). The species diversity and the soil characteristics seem to be a relevant driver of community structure and highlight the importance of considering the effects of ecological legacy on management plans (Souahi et al. 2022). In this regard, favoring the expansion of shrubs patches on grasslands may promote communities with a higher nutrient concentration and species diversity, which could attenuate the potential impact of global change on these rangeland communities.
Investigating the functional and phylogenetic structure of communities can thus provide useful insight to understand the outcome of both environmental and historical factors involved in the species assembly. Thus, while the phylogenetic structure provides tools to explore the role of the different evolutionary histories of the species in the community assembly, the functional structure is closely related to the different ecological strategies used to cope with environmental filters (Tanaka and Sato 2015;Xu et al. 2017). For instance, phylogenetic and/or functional over dispersion can occur among communities situated in similar habitats as a direct consequence of competitive exclusion between similar species, while strong abiotic filtering may promote phylogenetic and/or functional clustering by the selection of certain traits or clades imposed by the environmental constraints (Cavender-Bares et al. 2009;Xu et al. 2017). According to Macheroum et al. (2021), the variation of phylogenetic diversity, estimated using taxonomic structures (i.e., genus/ species and family/species ratios), was related to local environmental conditions. Phylogenetic diversity increased with the increase in total vegetation cover.
However, the leading dimensions of ecological variation among plants include not only traits related with resource uptake but also those associated with competition for light and reproductive ability (Díaz et al. 2016). To this end, Westoby (1998) proposed the leaf-height-seed (LHS) strategy scheme, which captures the strategy of plants by three trait-related independent dimensions of ecological variation. Hence, the LHS scheme assumes that leaf economics is essentially independent of the dimensions of plant height and seed production (Laughlin et al. 2010).
Moreover, the community assembly patterns cannot be assessed by a single measure of functional composition (i.e., the dominant trait values expressed by the community weighted means; Garnier et al. 2004) but rather needs a multi-faceted approach (de Bello et al. 2013). Other descriptors of functional and phylogenetic structure may allow inferring relevant processes and can help predict community responses to global change (Salgado-Luarte et al. 2019; Chang and HilleRisLambers 2019). The understanding of how environmental variables and biotic factors influence phylogenetic and functional diversity of different plant communities and also ecosystem functioning remains very poorly investigated, especially in drylands. Thus, carrying out an extensive study based on complementary diversity indices and approaches (i.e., taxonomic, phylogenetic, or functional, e.g., Le Bagousse-Pinguet et al. (2019)), estimations of sampling effort and data analytical tools would help to better understand plant diversity of semi-arid steppe rangelands.
In this study, we tested these predictions: firstly, by examining the responses of the functional LHS trait composition at community level between habitats along the soil nutrient gradient; secondly, by comparing the species richness and the functional and phylogenetic diversity indices with the above mentioned factors.

Study site
Of the total land area of Iran, about 52.4% are rangelands, 8.6% are forests, and 19.5% are deserts including bare salty lands. There are three main climatic zones in Iran, including arid and semi-arid regions of the interior and far south, Mediterranean climate (mainly in the western Zagros Mountains, the high plateau of Azerbaijan, and the Alborz mountains), and humid and semi-humid regions (mainly in the Caspian but also in West Azerbaijan and the southwest Zagros). The study area is located in plant geography zone of Europe-Siberia (the Hyrcanian area). The study area, the southern hillsides of Damavand Mountain, is located in the north part of Iran (52°11′-52°5′ N; 35°46′-35°49′ E). It covers ca. 6000 ha and has an elevation stretching from 2500 to 34,600 m asl. The climate is semi-arid cold. The minimum monthly mean is − 0.3 °C, and the maximum monthly mean is 19.2 °C. The mean annual temperature is 12.7 °C, and the mean annual precipitation is 652 mm (Fig. 1).
The soil type of the study area is sandy-loamy and silty soil. The bedrock is primarily calcareous. The soils are variable, ranging from deep brown soils to shallow lithological soils, the latter exerting a greater influence on vegetation. Basic cations and carbonates promote an alkaline pedogenesis that can be followed by an acidic pedogenesis after carbonate lixiviation.

Experimental design
The experiment was carried out in three different habitats: grassland, dominated by Agropyron repens, Festuca ovina, Bromus tectorum, and Bromus tomentellus; shrubland, dominated by Artemisia aucheri, Astragalous aegobromus, Eryngium billardieri, and Thymus kotschyanus; and a transitional fringe between both habitats, which we called "mixed area," dominated by Astragalus ochrodeucus, Festuca ovina, Onobrychis cornuta, and Bromus tomentellus. These habitats are adjacent to each other in our study area, and they have different community compositions, structures, and dominant species. From mid spring and summer of 2016, sampling was down in a 1 × 2-m plots within each habitat. In each habitat, in total, 150 sampling plots were selected randomly with the constraint that it was at least 1 m from the margin to avoid edge effects. On the field, each plot was separated at least by 100 m from the others,

Plant functional trait measurements
We select the 44 most abundant species. Between four and 20 individuals per species were measured in sites (Dubuis 2013). One well-developed entire leaf was collected per individual for trait measurement. The leaves were scanned and analyzed in the lab using the leaf area matter software to measure their surface. We also measured the leaf dry mass after drying and weighting the leaves. The vegetation height (VH, m) was measured in the field as the distance between the top of the photosynthetic tissue and the ground. Specific leaf area (SLA; mm 2 mg −1 ) was calculated as the ratio of the leaf surface to its dry mass (Rossier 2011;Dubuis 2013). Seed mass (SM; mg) was compiled from seed databases (Goethe University Frankfurt, http:// www. seed-dispe rsal. info, and the Royal Botanic Gardens, http:// data. kew. org/ sid/). For the few species that we did not find seed mass values, we estimated the seed mass of each of these species from the average value of the species of the same genus that were registered (de la Riva et al. 2019), because seed mass is strongly conserved through the phylogeny (Lord et al. 1995). Given the multidimensionality of plant functions, we chose these three traits because of their importance in providing information about different independent axes of ecological strategies (see Díaz et al. 2016;Bruelheide et al. 2018). Specifically, we selected SLA, VH, and SM as a subrogate of the plant strategies related with resource acquisition, competition for light, and reproduction, respectively (Laughlin et al. 2010).

Soil collection and processing
Soil samples of the top 30 cm of depth (where nutrient uptake mostly occurs; Jobbágy and Jackson 2001) were collected from each sampling plot. Soil samples were randomly collected from three points using a bucket auger and mixed into a single soil sample and were brought into the laboratory in airtight plastic bags. All of the soil samples were air-dried and then filtered through a 0.2-mm sieve, discarding the visible roots and other plant debris. Soil pH was measured using a pH meter with a glass electrode (soil/KCl ratio is 1 V 2:5). Soil organic carbon was determined using the Walkley-Black method, and a factor of 1.3 was applied to adjust the organic C recovery (Nelson and Sommers 1982). The soil available nitrogen, phosphorus, and potassium were measured using the methods of Miller and Keeney (1982).

Statistical analysis
Following Garnier et al. (2004), the traits were weighted by the relative abundance of their constitutive species to calculate the community weighted mean (CWM) in each plot. The phylogenetic tree was obtained with the comprehensive angiosperm species-level phylogeny from Zanne et al. (2014), as updated by Qian and Jin (2015), which is included in the R package "S. PhyloMaker" (Qian and Jin 2015). The distances of the few species that were not found in the PhytoPhylo were supplanted by the distance of the closest species of the same genus found in the megaphylogeny tree (de la Riva et al. 2019). For each plot, we calculated the species richness (Pielou 1969;Marcon 2016), which represented by the total number of species per plot, the Faith's phylogenetic diversity index (PD), the functional richness (Frich), and the mean pairwise dissimilarity (MPD) index for both functional and phylogenetic diversity. PD was calculated with the R package "PICANTE" (Kembel et al. (2010)). To quantify the Frich index, we built a functional space through a principal coordinate analysis (PCoA), using a Euclidean trait matrix after standardizing traits (mean = 0, SD = 1) (Maire et al. 2015). We calculated the functional and phylogenetic MPD indices (FMPD and PMPD, respectively) based on existing algorithms (de Bello et al. 2016). Faith's phylogenetic diversity (PD; Faith 1992) was calculated for each community using the R package Picante. PD measured community phylogenetic richness and was calculated as the sum of the lengths of all those branches that are members of the corresponding minimum spanning path (Faith 1992). As with functional diversity, other methods that have been developed to represent phylogenetic diversity usually fall under the categories of pairwise distance methods or nearest-neighbor distance methods (Swenson, 2014, Macheroum et al. 2021. Mean pairwise distance, or MPD, is the most commonly used pairwise distance metric. As with Faith's PD index, MPD can be weighted by the abundances of the species in the assemblage. These four indices integrate the main components of community structure enclosed by other related indices. In order to account for trait coordination and to identify the major axes of trait variation, we run a principal component analysis with the three CWM trait values (VH, SM, and SLA) (Biotic_PCA). In addition, to reduce the number of variables characterizing the soil environment and their collinearity, other PCA was performed with the soil variables measured in this study (Soil_PCA). When necessary, variables were previously standardized and log-transformed to fulfill assumptions of normality (based on Kolmogorov-Smirnov) and homoscedasticity. We used an ANOVA analysis to check for differences on the distribution of the three types of habitats (grassland, shrubland, and mixed) along the principal component analysis (PCA) in both biotic and soil PCAs. Post hoc Tukey tests were performed to check the significance of differences between the least square means of each type of habitats. We also performed a one-way ANOVA to check for the variation of the diversity indices among habitats.
In order to evaluate shifts in the functional and phylogenetic structure distributions in response to edaphic factors, we used the generalized least square model (GLS) controlling for the potential spatial correlation between soil variables associated with the type of habitat. We included a spatial autocorrelation structure using the cor-CAR1 function. We built one generalized least square regression model per community property; for each CWM trait and diversity index (respond variable) with the first and second Soil_PCA axes as predictors, using the gls function of the nlme package (Pinheiro et al. 2016). The variance explained by the GLS model was estimated by extracting the marginal R 2 value, which was calculated with the r squared function in the "piecewise SEM" package (Lefcheck 2016).
In addition, we used null models to check if the functional and phylogenetic diversity patterns were a trivial consequence of species richness variation. We randomized the trait combinations across species, while fixing the number of species of each plot (999 runs) and estimating the same functional and phylogenetic indices. For each index, we estimated the standardized effect size (SES, = (observed difference − mean null difference)/SD null differences), in order to assess if the relationship observed between the functional and phylogenetic indices and the soil variables was significantly different from those obtained between randomized communities (null differences).
All the statistical analyses were performed using the R 3.1.3 software (R Development Core Team 2014).

Results
Trait space distribution at community level was relatively well captured by the first two axes of the Biotic_PCA (Fig. 2). The first PCA axis accounted for 39.4% of overall variation and was related mostly with vegetation height. The second principal component, which explained 33.8% of the overall variance, was represented mainly by variations in SLA and SM. In addition, the first Soil_PCA axis (explaining 37.1% of the total variance) showed a high loading of soil nutrients (N, P, K, OM) and pH, while the second axis (which explained 22.1% of the variance) exhibited a high loading of pH at the top and K at the bottom (Fig. 3). The positions of the plots along the first and second Soil_PCA axes reflect the segregation of the different types of habitats ( Fig. 3 and Table 1), while no significant differences were found among habitats for the Biotic_PCA (data not shown). Thus, shrublands were significantly segregated from the other habitats, showing the highest nutrient concentration. In contrast, the grasslands show significantly the lowest nutrient concentration ( Fig. 3 and Table 1). Significant differences were also observed in the second axis, where the mixed ecosystem was segregated by its higher values of pH, while the shrubland was in the opposite end with higher values of K.
The PCA of soil environmental factors showed a first main axis associated with soil fertility, while the second axis represents a trade-off between pH and soil K, and both axes, from which plots segregate, were the most important factors associated with diversity variation. Appendix S1 shows that soil properties did not significantly capture the communitylevel traits (Fig. 4 and appendix S1). Also, box plot shows the difference of species richness and phylogenetic indices among habitats. No significant difference was observed in species richness between grassland and mixed site (Fig. 5).
In this regard, the higher soil nutrient concentration (first Soil_PCA axis) and pH (second Soil_PCA axis) were significant and positively related with species and functional richness and Faith's phylogenetic distance, while no significant results were obtained for any functional trait at community level (SLA CWM , VH CWM , and SM CWM , data not shown) or for the functional and phylogenetic mean pairwise distance (appendix S1). Nevertheless, the results of this significant relationship were not consistent when applying the null model (SES). Therefore, the null model suggests that changes in functional richness and Faith's phylogenetic distance were the reflection of the variations in species richness (Appendix S2).

Discussion
Our results report several insights on the relationship between vegetation structure, soil nutrient concentration, and land use in Damavand Mountain. The LHS strategy scheme assumes a certain degree of independence between these functional traits along three main dimensions (Laughlin et al. 2010;Reich 2014). Yet, our approach to this scheme focused on the community level does not uphold it. Our results show that trait distribution was relatively well capture by the first two axes of a principal component analysis (PCA) on which community trait values converge: one axis shows the positive covariation between SLA and plant height, while the positive covariation between SLA and seed mass could be loaded on separate axes. Surprisingly, community-level trait composition was poorly captured by the habitat or soil variables. This result suggests that species with opposing trait values co-exist in the same habitat, which is unexpected. Similar results under different grazing regimes were obtained in arid environments (Vesk et al. 2004;Kouba et al. 2021;Macheroum et al. 2021) or fertility gradients (Rusch et al. 2009), contrasting with others (Laliberté et al. 2012;Vaieretti et al. 2018). It is widely recognized that increased and intense grazing leads to vegetation homogenization at large scale and that grazing exclusion may increase the spatial variation in species assemblages.
The discrepancy with our expectations could be due to the fact that we deliberately chose similar climatic conditions, to focus on grazing disturbance and soil nutrient availability. The study of biodiversity and its variables is important to deepen our understanding of arid and semi-arid ecosystems and then assist rehabilitation activities through grazing exclusion supported by plantations (Macheroum et al. 2021). Arid conditions promote specific adaptations to the imposed climatic constraints, leading to relatively small differences in terms of functional syndromes (Reu et al. 2011;Bruelheide et al. 2018;de la Riva et al. 2018), and thus, despite some diversification of strategies exist, they are limited within the range of viable traits that allow plants to persist in this arid environment (de la Riva et al. 2016). Hence, aridity may limit the success of species that are not physiologically able to tolerate such water scarcity, reducing the range of functional traits. We acknowledge that our study is observational, and therefore, we cannot demonstrate if the results are a consequence of water deficit or of pre-existing patterns in the species assemblage, such as land use history. However, the fact that there were no different community trait syndromes among habitats suggests that the climatic filtering acts as dominant structuring processes of community trait composition, promoting strong trait convergence independently of the habitat (Gross et al. 2013;Bruelheide et al. 2018), while the local variation in environmental filters tends to alter species composition rather than functional composition,  which has been previously observed in areas with relatively long disturbance history (Cingolani et al. 2005;Conti and Diaz 2013).
Our results showed that soils from shrublands have higher nutrient concentration than grasslands. Trapping by livestock can lead to congestion and changes in levels of infiltration, bulk density, and reduced destructive activities (Li et al. 2011;Yu and Jia 2014). Moreover, under longterm pressure of grazing, some energy and nutrients are transmitted to the diet (Li et al. 2011;Lu et al. 2015), while Box plot shows the variation of species richness and phylogenetic indices among habitats. The line inside the box represents the median value, the box limits are the 25th and 75th percentiles, error bars show 10th and 90th percentiles, and filled symbols show outliers. Different letters denote significant differences between habitats (Tukey's test). Faith's phylogenetic diversity (PD), functional richness (FR), and species richness (SR) shrublands can collect a large amount of soil organic matter due to the increase of litter on the soil surface by the lowering of bed decomposition (Nieder and Benbi 2008). This could be the consequence of the reduction segregation of the nutrient concentration among habitats. Overall, it seems that the phylogenetic and functional patterns responded more to the taxonomic richness variations than to trait-based assembly processes. These results are consistent with other studies in grazing areas where variations of functional diversity were related with species diversity (Niu et al. 2015;Simova et al. 2015;Salgado-Luarte et al. 2019). However, we observed an overall low covariation between species richness and functional and phylogenetic mean pairwise dissimilarity. This implies that the mechanisms that support the coexistence of the species not necessarily support the functional or phylogenetic differentiation among those species (de Bello et al. 2006), and thus, while aridity could more strongly constrain the overall functional and/ or phylogenetic dissimilarity among species, competition for available resources could induce a niche differentiation among plants reaching the higher species richness in more productive plots (Turner 2008;Ceulemans et al. 2017). This pattern is especially common on resource-poor environment, where an increased availability of limiting resource may provide niche opportunities for coexisting plant species (Willems et al. 1993;Turner 2008;Souahi et al. 2022).
Overall, our research shows how some deterministic processes influenced functional and phylogenetic changes related with species richness in the semi-arid ecosystem from Damavand Mountain. However, the different traits and dimensions used show some correlations, and this, together with the lack of certain traits, which may be relevant to estimate a wide range of the functional niche occupation (e.g., belowground traits or leaf nutrients), could be limiting our assessments. Therefore, further studies based on meaningful syndromes or strategies should provide stronger tests of the underlying mechanisms causing community changes under different grazing regimes, which could be considerable importance in the effort to conserve taxonomic, functional, and phylogenetic richness of these arid communities from west Asia.

Conclusions
In this study, we tried to discern the assemblage patterns of the semi-arid plant communities from Iran along natural nutrient gradient on three habitats. Our results show that the community-level trait composition was poorly captured by the habitat or soil variables. Soils from shrublands have higher nutrient concentration than grasslands, but this pattern has little effect over the functional or phylogenetic structure which tended to increase with nutrient content.
Also, this suggests that the arid conditions of the region promote specific adaptations to the imposed climatic constraints, leading to relatively small differences in terms of functional syndromes, while the local variation in nutrient content tends to alter species composition rather than functional composition. The results of this research can be used in major projects or implementation projects, because the functional and phylogenetic changes may help us to understand ecological processes, such as succession patterns and sustainability of communities.