This study was conducted along seven transects (four west-east and three north-south transects) is mainly alpine and subalpine grasslands located across the eastern and central part of the Qinghai-Tibetan Plateau. Due to the restrictions of terrain and available roads, the transects’ shape appears to be irregular and non-linear. The study area is 1,220 km long and 1,070 km wide, covering a geographic area within the latitude of 29.14–38.84° N and longitude of 88.15–101.51° E.

Field surveys were conducted during late July to early September in 2014 and 2017 and 118 sites were finally surveyed along the seven transects according to the criterion of every 50–60 km interval unless there were no grassland vegetations (e.g., forest, shrub, salt marsh, barren and stony land, lake, farmland, or residential area) (Figure 1 and Supplementary Appendix 1). Field time is the peak aboveground biomass and cover time for most species in the vegetation (Zhu et al., 2015, 2016). All sampling sites were used as yak and sheep winter grazing pasture (from late autumn to early spring) and no obvious summer grazing by livestock or rodents was observed. In total, 118 sites were selected which were as follows: 17 desert steppes, 19 steppes, 26 steppes meadow, 37 meadows, 12 wet meadows, and 7 cushion vegetation (altogether six grassland types). Moreover, seven plots were in the subglacial belt, 35 in the subalpine (including cold temperate), and 76 in the alpine belt. These sites represent an elevational range from 2,792 to 5,217 m, with the climatic range covering a gradient in mean annual temperature (MAT) of -9.04 to 6.07°C and mean annual precipitation (MAP) of 44–716 mm. For selected sites, the number proportion of each grassland vegetation or climate belt is roughly equal to its area proportion in the studied zone (Wu, 1980). In total, we recorded 719 angiosperm species which included almost all dominant species and more than 65% common species in the study region (Wu and Chen, 2004) and we used 719 as the regional species pool number.

For each site, we recorded the geographical coordinates, elevation, slope gradient, and aspect. Five 1 m × 1 m quadrats were laid at a distance of 30–40 m from each other within each 100 m × 100 m site and all vascular plant species were recorded. In the meadow, wet meadow and cushion vegetation individual plants have a very high abundance (up to 2,200/m²), thus, for aboveground biomass and abundance, we sampled a smaller0.5 m × 0.5 m quadrat located within the 1 m² quadrats and calculated each species 1 m² abundance and biomass. In total, we sampled 581 quadrats due to the loss of plant material or recording problems in 9 quadrats of 6 sites.

For the vegetation measurements, in each quadrat we recorded the canopy cover, species number, individual plant number, height, and the percent cover and aboveground biomass of each species. Each species aboveground biomass was clipped to the soil level and then weighed after removing any dead parts (standing dead and litter) and dried at 75°C for 48 h to a constant weight.

We selected four traits tightly associated with plant life strategy and functional tradeoff: specific leaf area (SLA), leaf size, plant height, and seed mass. SLA is correlated with other plant traits such as leaf nitrogen content, leaf lifespan, relative growth rate, representing, thus, a trade-off between the leaf carbon acquisition rate and the longevity of plant tissues (Schöb et al., 2013; Qi et al., 2014; Lemoine et al., 2015; Rolo et al., 2016). Plant height and leaf size indicate a trade-off between the competitive or interception ability for light and support structure (e.g., stem, twig, and petiole biomass construction) cost; while seed mass signifies the trade-off between dispersal probability and a seedling establishment’s ability (Cadotte et al., 2009; Niu et al., 2014; Qi et al., 2015; Wagg et al., 2017). Thus, community-level distribution patterns (e.g., clustering or divergence) for these traits can disentangle the importance of these trade-offs during community assembly and help predict how assembly processes respond to environmental changes (Lemoine et al., 2015; Qi et al., 2015; Schmitz et al., 2015; Rolo et al., 2016; Wagg et al., 2017). There were 7,905 species records (samplings) in 581 quadrats (average 13.61 species per quadrat). For each sampling, 3–5 mature individuals and 5–20 fully developed but not senescing leaves were randomly chosen to measure plant height and leaf traits, respectively. The seed mass of most species was from a large database based on our previous studies. The methods of trait measurement referred to Qi et al. (2014, 2015), and their details are described in Supplementary Appendix 2.

For each quarter, we first calculated each species’ important value index (IVI) as the average of its relative biomass, relative cover, and relative density. We used IVI instead of the frequently used relative density because the former (1) represents more clearly the status and performance of species in a community, and (2) standardizes the data and makes the variance independent of the mean, which is better for this study with obvious differences in species composition and community structure among sites. Species diversity in each quarter was estimated with both species’ richness and the Shannon–Wiener index [H = –Σ (p_i × ln p_i), where p_i is the IVI of species i].

Rao’s quadratic index (FD_Q) was used to estimate the community functional diversity. FD_Q is the sum of the pairwise trait distances between species weighted by their IVI. FD_Q for single and multiple traits was calculated by using R and FDiversity (Casanoves et al., 2011). We used Rao’s quadratic index (FD_Q) for two reasons. First, Carboni et al. (2013) has advocated that this measure is more comprehensive for describing the community trait diversity. Secondly, FD_Q can be used to examine both single and multiple trait space, therefore better to predict variations in ecosystem functions (Butterfield and Suding, 2013).

We calculated the plant community phylogenetic diversity, which was assessed on the basis of a published angiosperm supertree (Phylomatic tree R20120829)¹. Branch lengths were assigned using the BLADJ algorithm in Phylocom 4.2 software (Webb et al., 2008). To quantify community phylogenetic diversity, we calculated the mean pairwise phylogenetic distance (MPD) using the “construct” function in Phylocom 4.2 software (Webb et al., 2008). The calculation of MPD is mathematically similar to FD_Q, and thus, the combination of MPD and FD_Q can provide a uniform and effective assessment of functional and phylogenetic community patterns. Based on the phylogenetic tree, we also calculated the phylogenetic signal of each plant trait, and the methods and results of the calculation were referred to in Supplementary Appendix 3.

Species number and abundance varied dramatically among sampled quadrats in our data. To remove any effect of sample size, we calculated the standardized effect size (SES) of MPD and FD_Q against a null model by generating 10,000 random assemblages for each community, while preserving the species IVI and shuffling the taxon names. For MPD, we calculated SES_M as: SES_M = (MPD_OBS – MPD_RANDOM)/sd(MPD_RANDOM), where sd(MPD_RANDOM) is the standard deviation of the random MPD values. For FD_Q, we calculated SES_FD as: SES_FD = (FD_OBS – FD_RANDOM)/sd(FD_RANDOM). A positive/negative SES_FD or SES_M indicates that traits or phylogeny are dispersed/clustered within a community. The difference in species, functional, and phylogenetic diversity among six grassland types is shown in Supplementary Appendix 4, in which meadow had highest species richness, aboveground biomass, SES_M and SES_FD of leaf size, SLA, and multiple traits, steppe had highest SES_FD of plant height and seed mass, desert steppe had lowest aboveground biomass, species richness and SES_FD of SLA, the wet meadow had lowest SES_M and SES_FD of seed mass, and cushion vegetation had lowest SES_FD of leaf size, plant height, and multiple traits.

Data on functional traits, aboveground biomass, and species richness were log-transformed to improve normality. We used regression analyses to test for the relationship among aboveground biomass, species diversity (species richness and H), and various functional and phylogenetic diversity indexes. For all regressions, we examined the fit of linear, quadratic, and piecewise models using stat and segmented packages from Software R 4.2.1, and identified the best fit model based on the explained variance (higher R²) and parameter significance (lower P-values). We constrained the maximum number of breakpoints to one per model to avoid overfitting. In piecewise models, the significance of the difference in regression slope between either side of breakpoint was evaluated by using Davies’ test. In addition, we used the difference in AIC (Akaike information criterion) value (ΔAIC) to compare among models. If ΔAIC was < 5 between linear and non-linear (quadratic and piecewise) models, models were considered not different (Segura et al., 2015), and the simple linear model was selected as the best model.

Communities’ phylogenetic or functional diversity may be structured by multiple niche processes (see the “Introduction”), resulting in a potential high variance among communities even with similar species diversity. Thus, there may be more than a single slope (i.e., rate of change) describing the relationship between species and phylogenetic (or functional) diversity (de Carvalho and Tejerina-Garro, 2015). Accordingly, we then used quantile regression analysis with phylogenetic or each functional diversity index as the response variable. The analysis can identify the limits, boundaries, and shifting within our bivariate distributions by estimating slopes not only through the median but also through each quantile (or percentile) of the relationship. We examined the nature of this upper (lower) bound with quantile regressions with data points in the 0.95 (0.05) quantile using Quantreg package (Koenker, 2018) for R. The significance of the slopes of quantile regressions was assessed with bootstrapped standard errors (with 999 permutations; Silva et al., 2016).

Our analysis evidenced an overall pattern of strongly phylogenetic clustering, but weakly functional overdispersion in the individual and multivariate trait space. The strongly phylogenetic clustering appears to be a result of a stochastic process, such as dispersal limitation. Tibetan Plateau is characteristic of high elevation and a series of huge mountains, which restricts the spread of alpine plants, resulting in the coexistence of closely related plant species with similar evolutionary history. However, to maintain the coexistence and avoid inter-specific resource competition because of similar ecological adaptation, these related species may undergo trait shifts (“character displacement”; Schluter, 2000), thereby reducing niche overlap. Some studies have evidenced that, in many cases, closely related species coexist in close proximity, but are well-separated ecologically and functionally (Lemoine et al., 2015; Qi et al., 2015; Lyu et al., 2017). Moreover, a degree of random functional trait distribution within the phylogeny may also contribute to the inconsistency between functional and phylogenetic community patterns to some extend (Kraft et al., 2007). For instance, the diversity (i.e., mean SES_FD value) of leaf size and plant height, the traits whose variation are restricted less by phylogeny (based on traits’ phylogenetic signal, seen in Supplementary Appendix 3), have more significant differences to phylogenetic diversity (i.e., mean SES_M value) in this study.

Our study supports hypothesis 1, i.e., an overall positive relationship between aboveground biomass and species’, functional and phylogenetic diversity, indicating the important role of biodiversity at different aspects in maintaining ecosystem functioning for the QTP grasslands. But surprisingly, the optimum species diversity-productivity relationship shows a two-stage pattern, in which, with increasing species diversity, communities’ productivity increases non-significantly at low diversity levels but significantly at middle and high levels. The pattern is opposite to general relationships (e.g., positively linear, logarithmic, unimodal, or neutral), and thus it cannot be explained by common hypotheses such as rivet, compensatory/keystone species, redundant species, or null hypothesis (O’Connor and Crowe, 2005; Kang et al., 2015), but rather a unique mechanism potentially adapted to plant communities along environmental pressure gradients. The mechanism emphasizes a species diversity threshold, and below it, extreme environmental conditions (e.g., extreme low temperatures and intense radiation in some alpine zones, and low soil nutrients and drought stress in desert steppe zones) filter out most species, resulting in the coexistence of few species. The life history strategy of surviving species is to adapt to these environments, resulting in the lack of strong interaction and significant niche complementarity among them. Thus, the biomass production of these communities may be largely determined by restrictive resources or environmental factors, rather than by species composition and diversity. In other communities, no resources are significantly limited, which allows the coexistence of multiple species with different ways of resource utilization. Increasing species diversity will enhance plant use efficiency of different resources, and ultimately increase community productivity (Liu et al., 2015; Duffy et al., 2017). These results, meantime, suggest the common finding that facilitative plant interactions drive high-altitude plant diversity and every species contributes to community functioning, especially when species diversity is high (Maestre et al., 2009; Michalet et al., 2014). On the contrary, the best model for explaining the functional diversity-productivity relationship for most traits is positively linear, confirming the general prediction that functional differences among coexisting species are the basis of niche complementarity, and a great inter-specific functional difference will enhance community-level resource utilization efficiency (Wagg et al., 2017; Mensah et al., 2018).

Contrary to our predictions (hypothesis 2) and most previous studies (Cadotte et al., 2009; Chalcraft et al., 2009; Lemoine et al., 2015; Liu et al., 2015; Rolo et al., 2016; Zhu et al., 2016; Rosenfield and Müller, 2020), phylogenetic and functional diversity explain less biomass variation of the QTP grassland communities than species richness. Strong environmental filtering effect on functional traits may have caused this inconsistency (Kraft et al., 2007; Rolo et al., 2016). In general, only fewer evolutionarily related species with similar morphological features can survive under stressful environments. In this case, species’ functional traits may better predict their response to environmental stress but not their community status and relationships to other species (Schöb et al., 2013; He and Bertness, 2014). For most QTP grassland communities, plants live under rich soil nutrition and light resources but harsh climate conditions. As a result, species’ adaption to climate, rather than interspecific competition for resources, may be the primary force structuring these communities, resulting in the decoupling of species traits and community function to a certain extent (Schmitz et al., 2015; Cadotte, 2017). Meantime, our result implies the importance of species diversity in maintaining the function of some special ecosystems even though coexisting species show similar morphological traits.

Consistent with hypothesis 3, our best model for explaining the species’ and functional or phylogenetic diversity relationship for the QTP grasslands is non-linear, in which functional or phylogenetic diversity (SES_FD of muli-traits and SES_M) increases significantly with species richness after it exceeds a threshold. The result demonstrates that high functional and phylogenetic divergence, representing significant species’ niche differentiation or interspecific difference (dissimilarity) in resource acquisition or other life-history strategies (Violle et al., 2007; Bu et al., 2014; Rolo et al., 2016; Lyu et al., 2017), is important to maintain species coexistence in high species diversity communities. In low species diversity communities, however, interspecific competition for resource and (or) interspecific life history differences are not strong determinants of species coexistence, and as a result, species’ niche differentiation may be determined mainly by stochastic processes, leading to less dependence of average trait divergence (dissimilarity) to species diversity.

It has been frequently reported that two opposite niche processes, species interaction and habitat filtering, may not be mutually exclusive but instead operate simultaneously across grassland communities (Anthelme et al., 2014; Michalet et al., 2014; Zhu et al., 2016). Quantile regression for the functional species diversity relationship may help to disentangle their effects in structuring community, in which the upper (or lower) bound should represent the maximum competitive effect (or filtering effect) due to the higher (or lower) functional divergence among species. Our result, however, shows that the relationship changed across quantiles, shifting from a significantly and steeply positive slope in the lower quantile to a slightly positive or negative slope for a single trait and even a near-zero for multi-traits or phylogeny in the upper quantile. This suggests a filled triangle for the limit of the bivariate distribution, implying that functional diversity differs significantly among communities when their species diversity is low but finally converges to be a constant with increasing communities’ species diversity. The reason for big functional diversity differences for low-species diversity communities may be that they are dominated by strong species competition, which excludes low-competitive species and increases inter-specific trait dissimilarity, or alternatively, by strong habitat filtering, which excludes low-adaptive species and increases inter-specific trait similarity. By decreasing the effect of species competition or habitat filtering, more species can coexist, and these opposite niche processes may operate simultaneously and tend to balance each other (He and Bertness, 2014; Fichtner et al., 2017). Accordingly, a constant interspecific trait divergence (similar to limiting similarity theory; Kembel and Hubbell, 2006), representing the result of species’ response to the combined effect of opposite niche processes, enables the community to adapt to the fluctuations of environment, to sufficiently use local resources, and to keep community productivity at a high level.

According to the best model, we show significantly different patterns in the relationship between functional diversity of plant height, leaf and seed traits and community productivity and (or) species diversity, suggesting the different roles of these traits in structuring the QTP grasslands. Leaf traits (leaf size and SLA) significantly contribute to achieving a predictive framework for community species diversity and ecosystem functioning. Leaf traits are tightly related to light resource acquisition and plant growth strategy and identified as a key functional predictor of the plant capacity of adaptation to changeful and stressful environments. The result is in accordance with previous studies suggesting that the diversity pattern of leaf traits, especially SLA, can be used to predict communities’ productivity, carbon storage and assembly rule (Long et al., 2011; Lemoine et al., 2015). By contrast, the relationships between diversity patterns of plant height and seed mass and community species diversity and productivity are generally weak or non-significant, implying that plant’s light interception capability and seedling establishment’s ability may be less important in structuring QTP grassland communities with a rich light resource for under plants. Moreover, most Tibetan grassland species are typically rosette (or semi-rosette) plants or primarily reproduce asexually (Nagy and Grabherr, 2009; Qi et al., 2014; Chen et al., 2017). These unique plants’ architectural and reproductive traits would weaken the effects of plant height and seed mass on supporting leaf light interception and enhancing seedling establishment and survival, respectively, and thus diminish their roles in driving community assembly. In addition, the piecewise model shows that communities with lowest functional diversity of leaf traits usually have no lowest productivity (Table 1, Figure 3), implying that, in some special zones (e.g., desert steppe) where light competition is not severe, similar leaf traits may be more helpful for coexisting species to adapt to the environment and use various resources.