Biodiversity explains maximum variation in productivity under experimental warming, nitrogen addition, and grazing in mountain grasslands

Abstract Anthropogenic global warming, nitrogen addition, and overgrazing alter plant communities and threaten plant biodiversity, potentially impacting community productivity, especially in sensitive mountain grassland ecosystems. However, it still remains unknown whether the relationship between plant biodiversity and community productivity varies across different anthropogenic influences, and especially how changes in multiple biodiversity facets drive these impacts on productivity. Here, we measured different facets of biodiversity including functional and phylogenetic richness and evenness in mountain grasslands along an environmental gradient of elevation in Yulong Mountain, Yunnan, China. We combined biodiversity metrics in a series of linear mixed‐effect models to determine the most parsimonious predictors for productivity, which was estimated by aboveground biomass in community. We examined how biodiversity–productivity relationships were affected by experimental warming, nitrogen addition, and livestock‐grazing. Species richness, phylogenetic diversity, and single functional traits (leaf nitrogen content, mg/g) represented the most parsimonious combination in these scenarios, supporting a consensus that single‐biodiversity metrics alone cannot fully explain ecosystem function. The biodiversity–productivity relationships were positive and strong, but the effects of treatment on biodiversity–productivity relationship were negligible. Our findings indicate that the strong biodiversity–productivity relationships are consistent in various anthropogenic drivers of environmental change.


| INTRODUC TI ON
Anthropogenic impacts such as increasing temperature, higher nitrogen addition, and overgrazing all conspire to cause rapid declines in plant biodiversity worldwide, especially in mountain grassland ecosystems, which naturally elicits concern about the consequences for the maintenance of ecosystem functioning (Chapin et al., 2000;Cingolani, Noy-Meir, & Diaz, 2005;Roth, Kohli, Rihm, & Achermann, 2013;Urban, 2015). The relationship between plant biodiversity and ecosystem function has been a major research topic in ecology for several decades, and while there is general empirical support for a positive effect of biodiversity on function from manipulative experiments (Balvanera et al., 2006;Cardinale et al., 2006;Tilman, Isbell, & Cowles, 2014), there is a lack of clarity about how anthropogenic changes in plant biodiversity might affect biomass production in more natural systems (Zavaleta & Hulvey, 2007). Inconsistent biodiversity effects on productivity could result from how biodiversity is measured, the confounding effect of environmental heterogeneity, and the nature of the anthropogenic impacts.
One reason might be that traditional biodiversity measures, like species richness, do not sufficiently capture the critical processes such as resource complementarity and interspecific interactions that are responsible for ecosystem function, which might be better reflected in relevant traits or evolutionary histories of species in a community (Lavorel & Garnier, 2002;Partel, Laanisto, & Zobel, 2007). Recently, a number of studies have shown that measures based on phylogenetic or single or multiple functional traits appear to be superior to species richness in explaining variation in productivity of plant communities (Cadotte, 2013;Cadotte, Cavender-Bares, Tilman, & Oakley, 2009;Flynn, Mirotchnick, Jain, Palmer, & Naeem, 2011;Liu, Zhang et al., 2015), and further supply direct links to the mechanisms controlling productivity (Cadotte, 2017). In addition to this, some studies (Liu, Zhang et al., 2015) found that statistical models that combined different biodiversity facets maximally explained the effects of biodiversity loss on ecosystem functioning or services.
For example, Liu, Zhang et al. (2015) found that multivariate functional diversity was the single predictor that consistently outperformed other single-biodiversity measures in explaining variation in productivity, but phylogenetic diversity and community-level plant height combined to explain maximum variation. However, beyond biodiversity facets that represent species-level differences, intraspecific variation is critically important to fully capture the diversity of plant communities (Albert et al., 2012). Ali and Mattsson (2017) evaluated the relative power of intraspecific and interspecific tree size variation and found that intraspecific variation better explained variation in aboveground biomass.
Although biodiversity is a major determinant of ecosystem productivity, the estimation of the biodiversity effect might be confounded by environmental factors and potential drivers of environmental change such as elevated temperature, nitrogen addition, and herbivory (Fridley, 2002;Hooper et al., 2005;Seabloom et al., 2017;Steudel et al., 2012;Tilman, Reich, & Isbell, 2012;Tilman et al., 2014). Thus, it is important to disentangle the relative importance of biodiversity relative to other drivers along an environmental gradient for inferring the consistent effects of biodiversity on the primary productivity of ecosystems. The majority of research on biodiversity effects on ecosystem function has been in experimentally assembled communities, and these studies generally support a positive relationship between biodiversity and ecosystem functioning (Hector et al., 1999;Liu, Zhang et al., 2015;Tilman et al., 2001). In contrast, biodiversity levels produced by an environmental gradient such as elevation might reveal different response of ecosystem productivity (Gough, Grace, & Taylor, 1994). Hence, the direct relevance of these experiments for estimating the impacts of realistic biodiversity loss due to environmental changes on ecosystem functioning remains controversial (Hector et al., 2007;Jiang, Wan, & Li, 2009).
To address the biodiversity-productivity relationships of natural communities under different anthropogenic impacts, we developed a fenced warming-fertilizing experiment in mountain wetlands along an elevation gradient on Yulong Mountain, Yunnan, China. We employed open-topped, passive warming chambers and urea fertilizer to simulate the projected global warming and nitrogen addition, respectively. We used a multimodel comparative approach to assess the relative contribution of single and various combinations of multivariate biodiversity indices, both with and without intraspecific variation, to predict the variance in biomass production after accounting for potential confounding factors including local environmental heterogeneity, warming, fertilizing, and grazing. We aimed to answer the following questions: (a) Does phylogenetic and functional diversity outperform traditional richness and evenness regardless of environmental heterogeneity and anthropogenic impacts? (b) Does incorporating intraspecific trait variability enhance the explanatory power of functional diversity? (c) Are biodiversity-productivity relationships comparable in experimental warming, nitrogen addition, and grazing along environmental gradient of elevation in mountain grasslands?

| Study sites and experimental design
We established eighteen study sites in south-facing wetlands of regular topology of Yulong Mountain (100°10′E, 27°00′N) along an environmental gradient of elevation (2,700, 3,200, and 3,400m) within the Lijiang Alpine Botanical Garden of the Kunming Institute of Botany, Chinese Academy of Sciences in Lijiang, Yunnan Province, China. Yulong Mountain has the mean annual temperature of 12.8°C and the annual rainfall is 935 mm, which is mainly distributed from July to October with distinct dry and rainy seasons (Luo et al., 2016).
Plant communities of wetlands have obvious species turnover along the elevation gradient with the dominance of the genera Isachne, Juncus at the lowest elevation, the genera Ligularia, Agrostis at the middle elevation, and the genus Agrostis at the highest elevation. All three wetlands have long livestock-grazing histories, and each supports different types of livestock where sheep and horses graze at the lowest elevation, scalpers, and yaks graze at the middle and highest elevation, respectively.
We established six 12 × 12 m permanent fenced sites randomly distributed in wetlands within each elevation in May 2015 ( Figure 1).
Within each permanent site, we conducted a complete randomized block factorial experiment with each block of size 5 × 5 m. There were two factors of both experimental warming and nitrogen addition in each block and two levels for each factor. In both fertilized blocks, we applied urea fertilizer annually at the beginning of the rainy season approximately the end of May at a rate of 5 g m −2 year −1 . In both warmed blocks, we applied open top chambers (OTCs), commonly employed devices to study the effects of climate warming on ecosystems (Marion et al., 1997). Here, our open top chambers were octahedral frames made of angle iron, 1.5 m maximum diameter, and 45 cm height. Six sides of each open top chamber were fastened to transparent 1.5-mm-thick hard plastic with adjacent edges of two plastic pieces attached with adhesive. We regularly arranged two open top chambers and two corresponding plots of the similar area in four blocks with at least 3 m between the nearest edges of adjacent plots. Furthermore, we randomly positioned 3-4 plots around each permanent site with total of 20 grazed plots. Hence, there were five treatments (T C = control, T W = warming, T N = nitrogen addition, T WN = combination of warming and nitrogen addition, and T G = livestock-grazing; Figure 1).
We recorded species richness and their abundance in a rectangular subplot of 0.5 × 0.5 m from the center of each plot at the peak of the growing season in August 2016 ( Figure 1). We then harvested all the stems of each species in each subplot at ground level, dried, and weighed them to 0.1 mg to estimate biomass production (productivity).

| Environmental data
After cutting the stems to ground level, we collected soil core samples from three random locations in each subplot with a cylindrical soil auger (5 cm inner diameter, 15 cm length). We combined the three replicates from the same depth for each subplot as a single composite sample, dried it in the shade, and filtered it using a 2-mm sieve for stoichiometric analysis. We measured soil pH, concentration of nitrogen (N), phosphorus (P), and carbon (C) following the standard protocols (Sparks et al., 1996). Besides soil resources, we also collected climatic data for rainfall, air temperature, and air moisture using HOBO RG3-M, HOBO Pro v2, respectively (Onset Computer Corporation, Bourne, MA, USA) from July to October in 2016. For each elevation, we placed one HOBO RG3-M and two HOBO Pro v2, of which one was positioned inside an open top chamber and the other one was positioned in a control plot. We showed the detailed F I G U R E 1 Map of the study sites on Yulong Mountain, Lijiang, Yunnan Province, China, and the plot design. Shown are treatments: T C = control, T W = warming, T N = nitrogen addition, T WN = warming and nitrogen addition, and T G = livestock-grazing distributions of temperature and moisture during the experimental interval in Figure A1.
These plant traits might reflect fundamental resource complementarity and interactions among co-occurring species (Weiher et al., 1999;Wright et al., 2004). We recorded plant height of maximum five randomly selected individuals from each species in each subplot. We calculated the maximum of plant height for each species per plot for intraspecific variability among plots. We scanned at least 1 mature leaf of randomly selected five individuals per species in each subplot using an Epson-V200 scanner. We then measured leaf area with image analysis software (ImageJ; http://rsb.info.nih.gov/ ij). We weighed the leaves after dried to a constant weight at 60°C to 0.1 mg and calculated the specific leaf area as the ratio of leaf dry mass to leaf area. We pooled the leaves from different individuals of the same species and measured leaf carbon, nitrogen, and phosphorus content. For the missing traits data due to rare species, we substitute the average of the same traits of the same species or same genus within the same treatment.
We constructed the phylogeny for the 105 species recorded in our study using rbcL + matK regions of the chloroplast genome. The detailed descriptions of DNA extraction, amplification, and sequencing are provided in Liu, Yan et al. (2015). Here, we briefly described the inference method of phylogenetic reconstruction. We aligned the rbcL and matK sequences using MAFFT (Katoh & Standley, 2013) and concatenated matK to the rbcL to form a super matrix. We used the sequences from the same genus in BOLD as the substitutes for the missing sequences in 27 of the species. For each gene, we selected top-ranked maximum-likelihood model of nucleotide substitution using Akaike's information criterion, as implemented in the function modelTest in the phangorn library (Schliep, 2011) in R (R Core Team, 2016). Then, we estimated a maximum-likelihood phylogeny using PhyML 3.0 with the starting-tree estimated from the BioNJ (Guindon et al., 2010). We chose one representative of early diverging angiosperm lineage Amborella trichopoda as the root of phylogeny and then used a semiparametric rate-smoothing method to transform the phylogeny to an ultrametric tree using the chronopl function with parameter value 1,000 in the R ape library (Paradis, Claude, & Strimmer, 2004).

| Measures of biodiversity
Using species composition and number of individuals, we calculated traditional species richness (S) and Shannon's evenness index (H′) for each subplot. We also calculated a suite of single and multivariate functional diversity metrics based on plant traits, as well as phylogenetic diversity metrics using the maximum-likelihood phylogeny.
We listed the detailed descriptions of the measures of biodiversity in Table A1. Here, we give a brief description of important functional and phylogenetic metrics. To assess the potential effect of intraspecific trait variability, we averaged the traits for each species  (1996) Note. The order from top to bottom for the measures of biodiversity represents their relative ranking using Akaike's information criterion weights.
TA B L E 1 Measures of biodiversity for general multivariate linear mixed-effect models across all subplots in the study as its "fixed" traits and averaged the traits for each species in a given subplot as its "specific" traits. We then calculated a number of functional diversity metrics including single community-level plant traits and multivariate functional diversity metrics for each subplot using both "fixed" and "specific" traits (Leps, de Bello, Smilauer, & Dolezal, 2011). Here, multivariate functional diversity metrics included Rao's quadratic entropy

| General linear mixed-effect models
We constructed a series of general linear mixed-effect models to determine the most parsimonious relationships between productivity and the various measures of biodiversity, treatment, and local environmental factors including soil resources. We assumed that various measures of biodiversity, experimental treatments, and soil resources as fixed factors, whereas elevation, treatment, and plot were treated as hierarchical random factors. Here, the use of a normal distribution of model residuals was validated based on the normalized scores of standardized residual deviance (Q-Q plots). To evaluate model support, we used Akaike's information criterion corrected for small sample sizes (AIC c ; Burnham & Anderson, 2002. We also used the marginal R 2 values of the models (R m 2 ) as a measure of the model's goodness of fit (Nakagawa & Schielzeth, 2013).
To search for the most parsimonious models explaining patterns of biomass production, we firstly removed redundant predictors associated with phylogenetic, functional diversity metrics. We selected the relatively better-ranked single-biodiversity metric models in both phylogenetic and functional diversity metrics. Meanwhile, to testify whether experimental treatments affect biodiversity-productivity relationships, we regressed biomass production against each biodiversity metric with the addition or multiplication of treatment and compared the explanatory ability of these models using Akaike's information criterion weights. The detailed single-biodiversity model ranking is listed in Table A2 and the biodiversity metrics we used in the following model construction are listed in Table 1.
Because of the strong correlation between most biodiversity indices (Spearman's ρ > 0.3; Table A3) and because multivariate functional indices are derived from the same trait data, we avoided including more than one of these like indices in any one model. Notes. Fixed factors are number of species (S), Shannon's evenness (H'), and phylogenetic diversity (IAC, imbalance of abundance at the clade; MNND, mean nearest-neighbor distance), and communitylevel mean of single functional traits (H max , plot-specific maximum plant height; LN, mean leaf nitrogen content value for individual species used for all plots where the species is found) or multivariate functional trait indices (RaoQ, Quadratic entropy; FDis, Functional dispersion: weighted distances from a weighted centroid in multitrait space), and experimental treatments (T: T C = control, T W = warming, T N = nitrogen addition, T WN = warming and nitrogen addition and T G = livestockgrazing), and soil resources (C, soil carbon content; N, soil total nitrogen content; P, soil total phosphorus content). Hierarchical random factor is elevation (2,700, 3,200, and 3,400 m), treatment, and plot. Values are shown for the estimated number of model parameters (k), maximum log-likelihood (LL), and the information-theoretic Akaike's information criterion corrected for small samples (AIC c ), change in AIC c relative to the top-ranked model (ΔAIC c ), AIC c weight (wAIC c , model probability), and the marginal and total variance explained (R m 2 , R c 2 ) as a measure of the model's goodness of fit. The top 10 models are listed; the full table is shown in Appendix: Table A3.
TA B L E 2 General linear mixed-effect model (GLMM) results for biomass production as a function of several fixed factors and a hierarchical random factor following the same constraint of correlation among soil resources and between selected biodiversity metrics and soil resources.

| Comparisons between biodiversity metrics
As expected, phylogenetic and functional diversity indices alone outperformed traditional species richness and Shannon's evenness to explain the variation of biomass production when simultaneously considering elevation and treatment (Table 1, Table 2, Tables A2   and A4). Although only several functional diversity indices (H max , RaoQ, FDis, FDiv, detailed information see in Table A1) considering intraspecific variability attained greater model support than corresponding indices using species mean traits (Table 1, Table A2), most of these indices were selected as relatively better-ranked singlebiodiversity metrics (Table 1). Of all functional diversity indices, the community-level mean of "specific" maximum plant height (H max ) on average accounted for the most explained variation in biomass production (R m 2 > 50%; Table 1). Phylogenetic diversity (IAC) was the top-ranked single-biodiversity metric of all considered biodiversity metrics here and explained the most variation in biomass production (R m 2 > 66%; Table A2).

| Biodiversity effects
Of the 166 multivariate linear mixed-effect models, the most parsimonious model included species richness (S), phylogenetic diversity (IAC), the community-level mean of "fixed" leaf nitrogen content (LN f ), soil carbon content (C), and treatment (T) accounting for >91% of the deviance explained in productivity (Table 2). After accounting for confounding effects of environmental factors and experimental treatment, biomass production generally increased with increasing species richness, phylogenetic diversity, and the community-level mean of "fixed" leaf nitrogen content (Figure 2a-c).

| Environmental and treatment effects
We found relatively weaker environment and treatment effects on biomass production compared to those of selected biodiversity F I G U R E 2 Scatter plots of the best-supported variables combined in the general linear mixed-effect models to predict variation in biomass production: (a) species richness (S), (b) imbalance of abundance at the clade (IAC) based on a maximum-likelihood phylogeny, (c) community-level mean of mean leaf nitrogen content value for individual species used for all plots where the species is found (LN), (d) soil carbon content (C), and (e) experimental treatments (T C = control, T W = warming, T N = nitrogen addition, T WN = warming and nitrogen addition, and T G = livestock-grazing). Dashed lines are linear regression lines, gray ribbon are their confidence intervals, and points and error bar in (e) are predicted values and their confidence intervals using general linear mixed-effect model metrics, but few treatment effects on biodiversity-production relationship ( Table 2, Table A2). The top-ranked model showed that grazing strongly reduced the biomass production compared with nitrogen addition; however, nitrogen addition and experimental warming showed no impact on biomass production ( Figure 2e). We also found evidence for a weak negative relationship between biomass production and soil carbon content (Figure 2d).

| D ISCUSS I ON
Our results show that phylogenetic and functional diversity alone outperformed traditional biodiversity measures, species richness, and Shannon's evenness, for explaining variation in productivity.
This corroborates observational and experimental evidence that phylogenetic and functional measures better align with the mechanisms controlling community assembly and ecosystem function than taxonomic measures (Cadotte et al., 2009;Flynn et al., 2011;Liu, Zhang et al., 2015). Of all considered functional biodiversity indices, a single functional trait was the single best predictor of productivity patterns. This is not surprising since single functional trait might explain a larger amount of variation in productivity than multivariate functional indices likely due to functional trade-offs and coordinated variation of functional traits (Cingolani, Cabido, Gurvich, Renison, & Diaz, 2007;Roscher et al., 2012).
Meanwhile, our study revealed that transitioning from using species mean (e.g., "fixed") traits to plot level (e.g., "specific") traits enhanced the explanatory power of functional diversity irrespective of plant traits in isolation or combination. Including specific traits allows us to detect subtle differences in functional diversity that respond to environmental variation that does not involve species turnover (Luo et al., 2016). Indeed, Jung et al. (2014) reported that the response of subalpine grassland communities to short extreme drought events was more mediated by intraspecific trait variability than species turnover. Intraspecific trait variability, through phenotypic plasticity, can promote species coexistence through providing fitness advantages and acting as a buffer against rapid climate change (Aspinwall et al., 2015;Nicotra et al., 2010;Valladares, Gianoli, & Gomez, 2007). This might lead to the shift in plant strategies in association with resource capture and use efficiencies at the local scale, which in turn are more related to plot-specific aboveground biomass production. Furthermore, phenotypic plasticity, especially associating with maximum plant height, might ameliorate light competition, which is assumed to be an important mechanism explaining species loss and biodiversity effects (Borer et al., 2014;Cadotte, 2017;Fridley, 2003;Hautier, Niklaus, & Hector, 2009;Zhou et al., 2017). Our results generally supported these assumptions and highlighted the critical role of intraspecific trait variability in more precisely predicting the ecosystem functioning in the face of global climate change.
Although functional diversity could explain a substantial proportion of variation in productivity, the combination of phylogenetic diversity and a functional trait (leaf nitrogen) attained more model support and greater explanatory power. This implies that functional diversity and phylogenetic diversity could complement each other in the perspective of ecosystem functioning because of their own limitations. Functional diversity was limited by the absence of potential key functional traits, for example, belowground root traits in our study (Cadotte et al., 2009). Linkage between phylogenetic diversity and real ecological differences remains unclear (Cadotte, Davies, & Peres-Neto, 2017). Thus, the influence of unmeasured plant traits might be compensated by metrics that capture phylogenetic information, such as the distribution of abundances at the clades or the equitability of abundance-weighted entropic measure of the distribution of evolutionary distinctiveness in an assemblage (Cadotte et al., 2010). Such a combination of functional and phylogenetic information for explaining biodiversity-productivity relationships has received support from both biodiversity manipulation experiments and natural ecosystems (Liu, Zhang et al., 2015;Zhou et al., 2017). For example, Liu, Zhang et al. (2015) found that phylogenetic diversity and plant height represented the most parsimonious combination to predict aboveground biomass production in a removal experiment where species richness and functional diversity were manipulated in alpine meadows of the Tibetan Plateau.
In this study, we found strong and positive effects of species richness on productivity in natural ecosystems after accounting for potential confounding factors. This was consistent with a review by Our results revealed that the drivers of environmental change had negligible effects on the relationship between biodiversity and aboveground biomass production. Our finding showed that the relationship between IAC and biomass production was consistently strongest for all considered biodiversity metrics in various treatments. Cadotte (2013) showed that biomass production was strongly predicted by phylogenetic diversity and that this finding might result from species complementarity, and ultimately species coexistence mechanisms (Chesson & Warner, 1981;Hodapp, Hillebrand, Blasius, & Ryabov, 2016;Horn & Macarthur, 1972;Levins & Culver, 1971). IAC that quantifies the relative deviation in the abundance distribution of a local community from a null distribution where individuals are evenly partitioned between clade splits can be used to infer the relative importance of competition and environmental filtering for local assembly. IAC would tend toward 0 if the strength of competition was proportional to phylogenetic relatedness, while IAC would be far greater than 0 if environmental filtering was key to community structure (Cadotte et al., 2010). Meanwhile, Cadotte (2017) showed that multidimensional trait measures might drive complementarity effect through niche complementarity, while few, singular traits (mainly height) might drive selection effect through interspecific competition. Our results were generally in line with these studies, because on the one hand, maximum plant height outperformed the multivariate functional indices alone in the perspective of ecosystem productivity, implying the importance of selection effect in biomass production in natural mountain grassland ecosystems; on the other hand, we observed IAC values far greater than 0, implying the dominance of environmental filtering in local community assembly, which might contribute to the role of selection effect in our system.
Our results point to the importance of both complementarity effects and selection effects for aboveground biomass production in natural mountain grassland ecosystems.

CO N FLI C T O F I NTE R E S T
None declared.

AUTH O R CO NTR I B UTI O N S
JJL conceived the idea and designed the experiment, JJL collected the data, JJL and MWC analyzed the data, JJL and MWC led the writing of the manuscript, and all authors contributed critically to the drafts and gave final approval for publication.  (5)  Notes. Fixed factors are number of species (S), Shannon's evenness (H'), and phylogenetic diversity (PD, sum of branch lengths; IAC, imbalance of abundance at the clade; EAED, equitability of abundance-weighted entropic measure of the distribution of evolutionary distinctiveness; MPD, mean pairwise distance; MNND, mean nearest-neighbor distance; MPD ed , weighted mean pairwise distance; MNND ed , weighted mean nearest neighbor distance); and community-level mean of single functional traits using interspecific and intraspecific functional traits indicated by subscript "f" and "s," respectively (H max , maximum plant height; LC, leaf carbon content; LN, leaf nitrogen content; LP, leaf phosphorus content; SLA, specific leaf area) or multivariate functional trait indices using interspecific and intraspecific functional traits indicated by subscript "f" and "s," respectively (FDis, functional distribution; FRic, functional richness; FEve, functional evenness; FDiv, functional divergence; RaoQ, quadratic entropy), and experimental treatments (T: T C = control, T W = warming, T N = nitrogen addition, T WN = warming and nitrogen addition, and T G = livestock-grazing). Hierarchical random factor is elevation (  Notes. Shown are species richness (S), Shannon's evenness (H'), phylogenetic diversity (PD, sum of branch lengths; IAC, imbalance of abundance at the clade; EAED, equitability of abundance-weighted entropic measure of the distribution of evolutionary distinctiveness; MPD, mean pairwise distance; MNND, mean nearest neighbor distance; MPD ed , weighted mean pairwise distance; MNND ed , weighted mean nearest neighbor distance); multivariate functional trait indices based on interspecific and intraspecific functional traits indicated by subscript "f" and "s," respectively (FDis, functional distribution; FRic, functional richness; FEve, functional evenness; FDiv, functional divergence; RaoQ, quadratic entropy), community-level mean of single functional traits also using interspecific and intraspecific functional traits indicated by subscript "f" and "s," respectively (H max , maximum plant height; LC, leaf carbon content; LN, leaf nitrogen content; LP, leaf phosphorus content; SLA, specific leaf area). Notes. Fixed factors are number of species (S), Shannon's evenness (H'), and phylogenetic diversity (IAC, imbalance of abundance at the clade; MNND, mean nearest neighbor distance), and community-level mean of single functional traits (H max , plot-specific maximum plant height; LN, mean leaf nitrogen content value for individual species used for all plots where the species is found) or multivariate functional trait indices (RaoQ, Quadratic entropy; FDis, Functional dispersion: weighted distances from a weighted centroid in multitrait space), and experimental treatments (T: T C = control, T W = warming, T N = nitrogen addition, T WN = warming and nitrogen addition, and T G = livestock-grazing), and soil resources (C, soil carbon content; N, soil total nitrogen content; P, soil total phosphorus content). Hierarchical random factor is elevation (2700 m, 3200 m and 3400 m), treatment, and plot. Values are shown for the estimated number of model parameters (k), maximum log-likelihood (LL), and the information-theoretic Akaike's information criterion corrected for small samples (AIC c ), change in AIC c relative to the top-ranked model (ΔAIC c ), AIC c weight (wAIC c , model probability), and the marginal and total variance explained (R m 2 , R c 2 ) as a measure of the model's goodness of fit.