-
PDF
- Split View
-
Views
-
Cite
Cite
Jia-Jia Liu, Yong Xu, You-Xia Shan, Kevin S Burgess, Xue-Jun Ge, Biotic and abiotic factors determine species diversity–productivity relationships in mountain meadows, Journal of Plant Ecology, Volume 14, Issue 6, December 2021, Pages 1175–1188, https://doi.org/10.1093/jpe/rtab064
- Share Icon Share
Abstract
Species diversity–productivity relationships in natural ecosystems have been well documented in the literature. However, biotic and abiotic factors that determine their relationships are still poorly understood, especially under future climate change scenarios.
Randomized block factorial experiments were performed in three meadows along an elevational gradient on Yulong Mountain, China, where open-top chambers and urea fertilizer manipulations were used to simulate warming and nitrogen addition, respectively. Besides species diversity, we measured functional diversity based on five traits: plant height, specific leaf area and leaf carbon, nitrogen and phosphorus contents. Several abiotic factors relating to climate (air temperature and precipitation) and soil chemistry (pH, organic carbon concentration, total nitrogen concentration and phosphorus concentration) were also measured. Generalized linear mixed-effect models were used to investigate the responses of species diversity and productivity to elevation, warming, nitrogen addition and their interactions. The effects of biotic and abiotic factors on the direction and magnitude of their relationship were also assessed.
Species diversity decreased with increasing elevation and declined under warming at mid-elevation, while productivity decreased with increasing elevation. Functional richness, maximum air temperature, soil pH and their interactions showed strong but negative influences on the species diversity–productivity relationship; the relationship shifted from positive to neutral and then to slightly negative as these sources of variation increased. Our study highlights the negative effects of short-term warming on species diversity and emphasizes the importance of both biotic and abiotic drivers of species diversity–productivity relationships in mountain meadow communities.
摘要
自然条件下物种多样性-生产力相互关系取决于生物和非生物因素,但其相对重要性及相互作用仍不清晰,特别是在未来的气候变化情景下。为此,我们在中国玉龙雪山3处不同海拔的高山草甸开展了模拟气候变暖和大气氮沉降的完全随机组块析因试验。除物种多样性外,我们根据株高、比叶面积、叶片碳、氮、磷含量计算了实验处理下草甸植物群落的功能多样性,并将其作为关键生物因素。此外,我们测量了气温、降雨以及土壤的化学属性作为潜在重要的非生物因素。我们利用广义线性混合模型研究了物种多样性和植物生产力对海拔、增温、施肥及其可能的交互作用的响应,同时评估了上述生物和非生物因素对物种多样性-生产力相互关系的影响。研究结果表明,物种多样性随海拔升高而降低并且在增温处理下有下降趋势且在中间海拔最为强烈。相对而言,植物生产力仅随海拔升高而下降。功能丰富度、最高气温、土壤pH对物种多样性-生产力相互关系表现出强烈的负交互作用,即物种多样性-生产力相互关系随着这些因素的增加从正相互关系变为中性关系,然后变为轻微的负相互关系。我们的研究指出短期增温对高山草甸物种多样性的负面影响,并强调生物和非生物因素决定了自然条件下物种多样性-生产力相互关系。
INTRODUCTION
The plant species diversity–productivity relationship (DPR) is known to play an important role in ecosystem functioning (Hector et al. 1999; Tilman et al. 2001; Whittaker and Niering 1975). Given that the rates of species loss are expected to rise with future increases in global temperatures and atmospheric nitrogen deposition (Nogues-Bravo et al. 2007; Payne et al. 2017; Sunday 2020), future decreases in ecosystem primary productivity also seem likely (Root et al. 2003; Stevens et al. 2004). Most studies on the DPR utilize artificial assemblages, where species diversity is manipulated either by filtering species from a particular pool or by removing them from natural assemblages (Hector et al. 1999; Huang et al. 2018; Liu et al. 2015; Pfisterer and Schmid 2002; Tilman et al. 2001), and the relationship is often a positive one. For studies in natural ecosystems, however, results are conflicting, where the DPR is positive (Flombaum and Sala 2008), neutral (Adler et al. 2011; Assaf et al. 2011) or even negative (Rose and Leuschner 2012; Zhang et al. 2016). Biotic and abiotic factors that determine the DPR in natural ecosystems, especially under future global change scenarios, remain poorly understood.
Biotic factors such as functional diversity, which measures the mean or variance of single or multiple trait values of species in a community, may mediate the relationship between species diversity and productivity in natural ecosystems (Roscher et al. 2012). A high level of functional diversity may suggest high trait dissimilarity among species in a community (Cadotte et al. 2009), which in turn, may be associated with a strong positive DPR when the diversity effect on productivity is driven by resource acquisition due to niche partitioning or facilitation among species (i.e. complementarity effect hypothesis; Loreau and Hector 2001). However, a high level of functional diversity may also be associated with a neutral or even negative DPR when productivity is mainly determined by the functional traits of dominant species in a community (i.e. selection effect or competitive exclusion hypothesis; Loreau and Hector 2001). Furthermore, the influence of functional diversity on the direction and magnitude of DPR might depend on how functional diversity is measured as well as which traits are included in analyses (Petchey et al. 2004). For example, community-weighted means of trait values (CWM) and trait dissimilarity among species in a community are often assumed to be more related to selection and complementarity effects, respectively (Cadotte 2017; Roscher et al. 2012). Although many studies have indicated that functional diversity may be more important than species diversity in predicting variation in ecosystem productivity (Cadotte et al. 2009; Flynn et al. 2011; Liu et al. 2015, 2018), relatively few studies have investigated the mediating effect of functional diversity on the relationship between species diversity and productivity in natural ecosystems (Roscher et al. 2013).
The DPR tends to be context dependent (Waide et al. 1999), and elevation, as one of the major natural gradients (Körner 2007), has the potential to regulate the relationship between species diversity and ecosystem productivity (Li et al. 2012). Communities at a higher elevation, in particular, can experience increased levels of environmental stress, such as low temperatures and soil fertility, coupled with high levels of solar radiation (Sanders and Rahbek 2012). It is plausible that plant species diversity in high elevations might have a stronger influence on ecosystem productivity than those in low elevations, due in part, to reduced competition and increased facilitation among species experiencing high levels of environmental stress (Brooker et al. 2008; Callaway et al. 2002; Mulder et al. 2001). A warming global climate and increasing atmospheric nitrogen deposition might also have effects on the DPR in plant communities (Cowles et al. 2016; Dong et al. 2019; Fornara and Tilman 2009; Jucker et al. 2016; Ma et al. 2010; Paquette et al. 2018; Yue et al. 2020). For example, warming experiments show that increased temperature can drive species loss through light competition or reduction in soil moisture (Klein et al. 2004), although the total species biomass may increase due to enhanced ecosystem CO2 uptake, even in the short term (Li et al. 2017; Lin et al. 2010). In addition, similar patterns occur with nitrogen addition, especially in long-term experiments (Dickson and Foster 2011; Niu et al. 2014; Wang et al. 2010). However, which abiotic factors (i.e. air temperature, precipitation, soil properties) drive variation in the DPR along the elevational gradient, or under global warming and nitrogen deposition, in natural plant communities has yet to be investigated (Kimmel et al. 2020; van der Sande et al. 2018).
The collective effects of biotic and abiotic factors on the DPR in natural ecosystems remain largely unknown, although a few experimental studies suggest that their interaction is important (Xu et al. 2018; Zhang et al. 2019). For example, there is some evidence that abiotic factors, including minimum temperature and light availability, are closely associated with plant community composition and functional traits along elevational gradients (Jiang and Ma 2015). Further, experimental warming often decreases soil water availability but increases soil available nitrogen concentration, which might drive the aggregation of the functional traits of drought-tolerant species (e.g. warm-season C4 grasses) and legumes (Cowles et al. 2016). Nitrogen addition usually leads to soil acidification and eutrophication, which might strengthen the dominance of faster-growing or taller species (e.g. gramineous) in a community (Hautier et al. 2009; Suding et al. 2005; Zhao et al. 2019b). Nonetheless, the relative importance, and potential interactions, of biotic and abiotic factors that underlie the direction and the magnitude of the DPR in natural ecosystems requires further experimental investigation.
In this study, we evaluated the DPR in natural alpine meadow ecosystems that may be impacted by future global climate change (Gottfried et al. 2012). Particularly, we performed a series of randomized block factorial experiments along an elevational gradient to simulate two main divers of global climate change (i.e. increases in temperature and atmospheric nitrogen deposition) that are likely to have far-reaching consequences on ecosystem structures and functions (Bloor et al. 2010; Day et al. 2008; De Boeck et al. 2007; Hutchison and Henry 2010; Lin et al. 2010; Sebastia 2007). We aim to address one overarching question: What are the effects of functional diversity, climatic factors, soil chemical properties and their interactions, on the DPR under climate warming and nitrogen addition along an elevational gradient in natural alpine meadow ecosystems? Here, we make the following predictions: (i) plant species diversity will decrease with increasing elevation and under experimental warming and nitrogen addition (Humbert et al. 2016; Klein et al. 2004; Körner 2007), while plant productivity will decrease with increasing elevation, but increase under experimental warming and nitrogen addition (Baldwin et al. 2014; Peng et al. 2020; Stevens et al. 2015), and (ii) functional diversity, climatic factors and soil properties will affect the DPR (Cowles et al. 2016; Li et al. 2012; Pires et al. 2018; Roscher et al. 2016). Specifically, decreasing temperature and soil fertility with increasing elevation will decrease functional diversity by increasing the ecological constraints on trait expression (Ottaviani et al. 2018). These might strengthen the positive relationship between plant species diversity and productivity by excluding inferior species through biotic interactions and environmental filtering (Cheng et al. 2020).
MATERIALS AND METHODS
Study site
We established our study sites in three different south-facing meadows along a gradient of elevation (2700, 3200 and 3400 m) on Yulong Snow Mountain (100°10′ E, 27°00′ N, Fig. 1a), China. The climate is oceanic monsoon with a mean annual rainfall of 935 mm and a mean annual temperature of 12.8 °C (~5.9 °C in January and 17.9 °C in July). The soil type is classified as wet iron bauxite (Gong 1999). The plant community is dominated by a few graminoids (Agrostis, Festuca and Kobresia) as well as dicots such as Primula, Artemisia, Trifolium and Ligusticum (Liu et al. 2018). The dominant animals observed in the region include livestock (e.g. yaks, horses), and various insect species (Zhao et al. 2019a).

Schematic representation of the randomized block factorial experiments along an elevational gradient on Yulong Mountain (a). Two factors of experimental warming and nitrogen addition and two levels for each factor (b). For experimental warming, two OTCs were used to raise the air temperature by 1.0–2.0 °C consistently along the elevational gradient. For nitrogen addition, half of the area of each block was fertilized using urea fertilizer at the beginning of the rainy season at a rate of 5 g/m2/year.
We established eighteen 12 m × 12 m permanent, fenced blocks (six blocks at each of three elevations, Fig. 1b) in May 2015; each block included four plots of 5 m × 5 m. At the center of two randomly selected plots, we arranged two open-top chambers (OTCs), a common device used for assessing passive warming effects on terrestrial biota (Marion et al. 1997). The air temperature in each of the OTCs reached 1.0–2.0 °C higher than the ambient temperature along the elevational gradient (Supplementary Fig. S1). We also applied urea fertilizer to half the area of each block at the beginning of the rainy season (~end of May) at a rate of 5 g/m2/year according to the regional distribution of atmospheric N deposition estimated in China (Xu et al. 2015). In total, we had four treatments: TC = control, TW = warming, TN = nitrogen addition, TWN = combination of warming and nitrogen addition. At the peak of the growing season in August 2016, we recorded species composition and abundance (i.e. the number of stems) in a rectangular quadrat of 0.5 m × 0.5 m measured from the center of each plot. We then harvested at ground level all the stems for each species in each quadrat; these were dried and weighed to estimate aboveground biomass production (i.e. productivity).
Trait data
We measured five plant traits that are important for plant resource acquisition and interactions (Liu et al. 2015): maximum plant height (Hmax, cm), specific leaf area (SLA, cm2/g), leaf carbon content (LC, mg/g), leaf nitrogen content (LN, mg/g) and leaf phosphorus content (LP, mg/g). We recorded the plant heights of five randomly selected individuals from each species in each quadrat. To estimate leaf area for the selected individuals, we measured 30 mature and undamaged leaves using an Epson-V200 scanner (Epson, Beijing, China) and image analysis software (ImageJ; http://rsb.info.nih.gov/ij). We then measured leaf dry mass for each of the leaves collected from the same individuals, and then pooled the leaves for each species to measure LC, LN and LP content (expressed as mg/g) using a potassium dichromate external heating method, an H2SO4–H2O2 digestion method and a Mo–Sb Antiluminosity method, respectively (Institute of Soil Sciences 1978).
Soil and climate variables
In August 2016, we used a cylindrical soil auger to collect five soil cores (5 cm diameter; 0–15 cm depth) along diagonal lines in the quadrat. We combined all soil cores per quadrat as a single composite sample. We air-dried the soil samples and filtered them using a 2-mm sieve for the analysis of soil properties. We measured soil pH in a 1:5 soil/water suspension, soil organic carbon (C, mg/g), soil total nitrogen (N, mg/g) and total phosphorus (P, mg/g) following the Walkley–Black wet oxidation procedure, the semi-micro Kelvin procedure and the sulphuric acid–perchloric acid digestion procedure, respectively (Institute of Soil Sciences 1978). Climatic data for air temperature and rainfall were also collected using HOBO Pro v2 and HOBO RG3-M, respectively (Onset Computer Corporation, Bourne, MA, USA) for each meadow from July to October in 2016.
Plant species diversity and functional diversity
We measured plant species diversity as the number of species per quadrat, a traditional measure of local species richness. We also calculated a number of functional diversity metrics: single community-level plant traits, multivariate metrics for each quadrat using the particular values of Hmax and SLA measured within a quadrat, and the averaged values of LC, LN and LP across all quadrats due to the practical difficulties to measure those leaf chemical traits by quadrat. We used this approach because although the species-level mean of functional traits is essential for plant biomass production and interspecific interactions (Cadotte 2017; Flynn et al. 2011; Roscher et al. 2012), recent work suggests that intraspecific trait variation, which reflects heritable genetic variation or phenotypic plasticity, is also important for ecosystem functioning and species coexistence (Jung et al. 2010; Mao et al. 2017). For single community-level plant traits, we calculated CWM for each trait by weighting the single trait values of species in a quadrat based on their abundance (Garnier et al. 2004). For multivariate functional diversity metrics, we calculated (i) functional richness (FRic), which measures the volume of the functional space occupied by the community, (ii) functional evenness (FEve), which measures the abundance-weighted pairwise functional distances, (iii) functional dispersion (FDis), which measures the weighted distances from a weighted centroid in multitrait space (Villéger et al. 2008) and (iv) Rao’s Q, which is the sum of pairwise functional distances between species weighted by their relative abundances (Rao 1982). Given the influence of functional diversity on DPR might be dependent on which traits are included in analyses, we also calculated each multivariate functional diversity metric using all possible combinations of the five functional traits measured in this study.
Data analysis
To explore how plant species diversity and productivity respond to experimental warming and nitrogen addition along the elevational gradient, we first determined whether experimental warming, nitrogen addition and elevation showed a strong influence on plant species diversity and productivity. To do this, we constructed a series of generalized linear mixed-effect models, where each variable (i.e. experimental warming, nitrogen addition and elevation) was set as fixed effects and blocks were set as random effects (i.e. Diversity/Productivity ~ Warming/Nitrogen addition/Elevation + [1|block]). Here, we constructed the corresponding null models assuming that plant species diversity and productivity showed no responses to experimental warming, nitrogen addition or elevation (i.e. an intercept model). The ‘glmer’ function in the lme4 library in R (Bates et al. 2015) was used to conduct this analysis. For plant species diversity, the use of a Poisson distribution of model residuals with a log link was validated based on the normalized scores of standardized residual deviances (Q–Q plots). For productivity, we validated the use of a Gamma distribution of model residuals with a log link. The model support of the single variable model against their corresponding null model was evaluated using Akaike’s information criterion corrected for small sample sizes (AICc) (Burnham and Anderson 2002, 2004). Once the strong predictors were selected for plant species diversity and productivity (ΔAICc = AICc of the null model − AICc of the model with predictors > 2), we constructed a series of multivariate models, including all possible combinations and interactions, and then used the top-ranked models with the minimum AICc as the final models. We then ran pairwise comparisons and multiple testing adjustments to determine whether two different treatments (e.g. 2700 vs. 3200 m a.s.l.) show significant differences in means using the function ‘glht’ in R package multcomp (Hothorn et al. 2008).
To determine whether there was a shift in the relationship between plant species diversity and productivity in our systems, we modeled plant productivity as a function of plant species diversity and its quadratic term (i.e. Productivity ~ Richness + Richness2) using a generalized linear mixed-effects model similar as above. To explore how FRic, climatic factors (i.e. mean, minimum and maximum air temperature and total precipitation), soil chemical properties and their interactions influence DPR, we first determined whether each single biotic and abiotic factor showed a strong effect on DPR. To do this, we first constructed a base model (i.e. a null model) assuming that plant productivity was a function of plant species diversity and a single biotic and abiotic factor (i.e. Productivity ~ Richness + Predictori; subscript i represents different biotic and abiotic factors). Then, we added an interactive term of plant species diversity and every single predictor into the above addition model to construct a corresponding interaction model (i.e. Productivity ~ Richness + Predictori + Richness:Predictori). Using both models for each single biotic and abiotic factor, we calculated the information-theoretic evidence ratio (ER), which explicitly evaluates the parameter bias-corrected likelihood of the null hypothesis against the alternative, as the AICc weight of the interaction model divided by the AICc weight of the addition model. Higher ER values provide stronger support for the alternative hypothesis versus the null hypothesis (Kass and Raftery 1995), which translates as a stronger influence of the single biotic and abiotic factor on DPR. We also tested whether experimental warming, nitrogen addition and elevation influenced the DPR in our systems following the same procedure as above. With the strong biotic and abiotic factors selected (ER >3), we considered all their possible interactive effects with species diversity (e.g. Productivity ~ Predictori + Predictorj + Richness + Predictori:Richness + Predictorj:Richness for two-way interactive effects; Productivity ~ Predictori + Predictorj + Richness + Predictori:Predictorj:Richness for higher-order interactive effects; subscripts i and j represent different biotic and abiotic factors). The generalized linear mixed-effect model was the same as above for productivity with blocks nested within each elevational gradient as random effects. To avoid serious collinearity, we removed the models constructed with any of the variance inflation factors of all parameters in a model greater than three (Cade 2015). Similarly, the remaining models were ranked using AICc, and the top-ranked model with the minimum AICc was chosen as the final model. For all models constructed, we used the marginal R2 values of the model () as the measure of the model’s goodness-of-fit (Nakagawa and Schielzeth 2013). All of the data analyses were performed in R 3.5.3 (R Core Team 2015).
RESULTS
Responses of plant species diversity and productivity to experimental warming and nitrogen addition along the elevational gradient
For plant species diversity, our results showed that elevation outperformed warming and nitrogen addition as the strongest influence, followed by warming; the effect of nitrogen addition was negligible (Table 1). The top-ranked model included the interaction between elevation and experimental warming (wAICc = 0.995), accounting for >67% of the variance explained in plant species diversity. The effect plot of the top-ranked model showed that species richness generally declined with increasing elevation (Fig. 2a). Further, experimental warming tended to decrease plant species diversity over the entire elevational gradient, however, the effect was more pronounced at the middle elevation. By contrast, plant productivity was only subject to the influence of elevation (Table 2), showing a moderately decreasing tendency with increasing elevation (, Fig. 2b).
Generalized linear mixed-effect models (GLMMs) results for species richness as a function of elevation, warming and nitrogen addition
Model . | k . | LL . | AICc . | ΔAICc . | wAICc . | . |
---|---|---|---|---|---|---|
E + W + E:W | 8 | −179.072 | 376.429 | 0.000 | 0.995 | 67.9% |
E + W | 6 | −186.875 | 387.043 | 10.614 | 0.005 | 61.2% |
E | 5 | −199.053 | 409.015 | 32.585 | <0.001 | 47.0% |
W | 4 | −216.098 | 440.792 | 64.363 | <0.001 | 15.2% |
1 | 3 | −221.525 | 449.403 | 72.973 | <0.001 | 0.0% |
N | 4 | −221.501 | 451.598 | 75.169 | <0.001 | 0.1% |
Model . | k . | LL . | AICc . | ΔAICc . | wAICc . | . |
---|---|---|---|---|---|---|
E + W + E:W | 8 | −179.072 | 376.429 | 0.000 | 0.995 | 67.9% |
E + W | 6 | −186.875 | 387.043 | 10.614 | 0.005 | 61.2% |
E | 5 | −199.053 | 409.015 | 32.585 | <0.001 | 47.0% |
W | 4 | −216.098 | 440.792 | 64.363 | <0.001 | 15.2% |
1 | 3 | −221.525 | 449.403 | 72.973 | <0.001 | 0.0% |
N | 4 | −221.501 | 451.598 | 75.169 | <0.001 | 0.1% |
Notes: Fixed effects are elevation (E: 2700, 3200 and 3400 m), experimental warming (W: Control and Warming), nitrogen addition (N: Control and Nitrogen addition) and their interactions. Random effects are experimental blocks. Shown are maximum log-likelihood (LL), the estimated number of model parameters (k), the information-theoretic Akaike’s information criterion corrected for small samples (AICc), the change in AICc relative to the top-ranked model (ΔAICc), AICc weighted (wAICc = model probability) and the marginal variance explained () as a measure of the model’s goodness-of-fit.
Generalized linear mixed-effect models (GLMMs) results for species richness as a function of elevation, warming and nitrogen addition
Model . | k . | LL . | AICc . | ΔAICc . | wAICc . | . |
---|---|---|---|---|---|---|
E + W + E:W | 8 | −179.072 | 376.429 | 0.000 | 0.995 | 67.9% |
E + W | 6 | −186.875 | 387.043 | 10.614 | 0.005 | 61.2% |
E | 5 | −199.053 | 409.015 | 32.585 | <0.001 | 47.0% |
W | 4 | −216.098 | 440.792 | 64.363 | <0.001 | 15.2% |
1 | 3 | −221.525 | 449.403 | 72.973 | <0.001 | 0.0% |
N | 4 | −221.501 | 451.598 | 75.169 | <0.001 | 0.1% |
Model . | k . | LL . | AICc . | ΔAICc . | wAICc . | . |
---|---|---|---|---|---|---|
E + W + E:W | 8 | −179.072 | 376.429 | 0.000 | 0.995 | 67.9% |
E + W | 6 | −186.875 | 387.043 | 10.614 | 0.005 | 61.2% |
E | 5 | −199.053 | 409.015 | 32.585 | <0.001 | 47.0% |
W | 4 | −216.098 | 440.792 | 64.363 | <0.001 | 15.2% |
1 | 3 | −221.525 | 449.403 | 72.973 | <0.001 | 0.0% |
N | 4 | −221.501 | 451.598 | 75.169 | <0.001 | 0.1% |
Notes: Fixed effects are elevation (E: 2700, 3200 and 3400 m), experimental warming (W: Control and Warming), nitrogen addition (N: Control and Nitrogen addition) and their interactions. Random effects are experimental blocks. Shown are maximum log-likelihood (LL), the estimated number of model parameters (k), the information-theoretic Akaike’s information criterion corrected for small samples (AICc), the change in AICc relative to the top-ranked model (ΔAICc), AICc weighted (wAICc = model probability) and the marginal variance explained () as a measure of the model’s goodness-of-fit.
Generalized linear mixed-effect models (GLMMs) results for biomass production as a function of elevation, warming and nitrogen addition
Model . | k . | LL . | AICc . | ΔAICc . | wAICc . | . |
---|---|---|---|---|---|---|
E | 5 | −386.050 | 783.009 | 0.000 | 0.843 | 13.64% |
1 | 3 | −390.654 | 787.660 | 4.651 | 0.082 | 0.00% |
N | 4 | −390.088 | 788.773 | 5.764 | 0.047 | 2.02% |
W | 4 | −390.644 | 789.885 | 6.876 | 0.027 | 0.04% |
Model . | k . | LL . | AICc . | ΔAICc . | wAICc . | . |
---|---|---|---|---|---|---|
E | 5 | −386.050 | 783.009 | 0.000 | 0.843 | 13.64% |
1 | 3 | −390.654 | 787.660 | 4.651 | 0.082 | 0.00% |
N | 4 | −390.088 | 788.773 | 5.764 | 0.047 | 2.02% |
W | 4 | −390.644 | 789.885 | 6.876 | 0.027 | 0.04% |
Notes: Fixed effects are elevation (E: 2700, 3200 and 3400 m), experimental warming (W: Control and Warming), nitrogen addition (N: Control and Nitrogen addition) and their interactions. Random effects are experimental blocks. Shown are maximum log-likelihood (LL), the estimated number of model parameters (k), the information-theoretic Akaike’s information criterion corrected for small samples (AICc), the change in AICc relative to the top-ranked model (ΔAICc), AICc weighted (wAICc = model probability) and the marginal variance explained () as a measure of the model’s goodness-of-fit.
Generalized linear mixed-effect models (GLMMs) results for biomass production as a function of elevation, warming and nitrogen addition
Model . | k . | LL . | AICc . | ΔAICc . | wAICc . | . |
---|---|---|---|---|---|---|
E | 5 | −386.050 | 783.009 | 0.000 | 0.843 | 13.64% |
1 | 3 | −390.654 | 787.660 | 4.651 | 0.082 | 0.00% |
N | 4 | −390.088 | 788.773 | 5.764 | 0.047 | 2.02% |
W | 4 | −390.644 | 789.885 | 6.876 | 0.027 | 0.04% |
Model . | k . | LL . | AICc . | ΔAICc . | wAICc . | . |
---|---|---|---|---|---|---|
E | 5 | −386.050 | 783.009 | 0.000 | 0.843 | 13.64% |
1 | 3 | −390.654 | 787.660 | 4.651 | 0.082 | 0.00% |
N | 4 | −390.088 | 788.773 | 5.764 | 0.047 | 2.02% |
W | 4 | −390.644 | 789.885 | 6.876 | 0.027 | 0.04% |
Notes: Fixed effects are elevation (E: 2700, 3200 and 3400 m), experimental warming (W: Control and Warming), nitrogen addition (N: Control and Nitrogen addition) and their interactions. Random effects are experimental blocks. Shown are maximum log-likelihood (LL), the estimated number of model parameters (k), the information-theoretic Akaike’s information criterion corrected for small samples (AICc), the change in AICc relative to the top-ranked model (ΔAICc), AICc weighted (wAICc = model probability) and the marginal variance explained () as a measure of the model’s goodness-of-fit.

Effect plots for the generalized linear mixed-effect models (GLMMs) of plant species diversity (a) and biomass production (b). Fixed effects are the interaction between elevation and experimental warming (Control and Warming) for plant species diversity and only elevation for biomass production. Random effects are experimental blocks nested in each elevational gradient. Hollow points represent the observed values for plant species diversity and biomass production associated with each experimental treatment and elevational gradient. Solid points represent their estimated means and error bars represent their 95% confidence intervals. Shown in the text are the marginal variance explained () as a measure of the model’s goodness-of-fit and group classification (e.g. a, b, c, etc.) based on the multiple comparisons of group means using Tukey contrasts. Two groups with no common alphabet between them (e.g. a vs. b) means a significant difference at a critical confidence level (P = 0.05).
Effects of functional diversity, climatic factors and soil chemical properties on the relationship between plant diversity and productivity
Elevation, rather than experimental warming and nitrogen addition, accounted for the variation in DPR in our systems (Fig. 3; Supplementary Fig. S2), which showed an evident nonlinear relationship (Fig. 4a). For all functional diversity metrics considered, only FRic and FEve showed strong influences on DPR (Supplementary Fig. S3); FRic based on SLA attained the highest explanatory power () and thereafter was treated as the focal biotic factor. The climatic and soil chemical characteristics along the elevational gradient are shown in Table 3. Among the abiotic factors, only soil pH, maximum air temperature (Tmax) and mean air temperature (Tmean) showed considerable effects on DPR (Supplementary Fig. S4). Thus, we treated Tmax and soil pH as the focal climatic factor and soil property, respectively. The top-ranked model identified for the interaction of the focal biotic and abiotic factors on DPR included all the interactions among species richness, FRic, soil pH and Tmax (i.e. S:FRic:pH:Tmax, wAICc = 0.796, Table 4; Supplementary Table S1), accounting for >48% of the variance explained in productivity. Although species richness, FRic, soil pH and Tmax had strong positive relationships with each other (Supplementary Fig. S5), the standard coefficients of the top-ranked model showed that the interaction term had the strongest but negative influence on plant productivity, i.e. the DPR shifted from positive to neutral and then to slightly negative with the increasing values of FRic, pH and Tmax. Only FRic had a significant effect on productivity after accounting for their interactive effects (Fig. 4b).
Elevation (m) . | Temperature (°C) . | Precipitation . | pH . | C (mg/g) . | N (mg/g) . | P (mg/g) . |
---|---|---|---|---|---|---|
2700 | 17.0 [7.56, 40.1] | 4140 | 5.75 ± 0.09 | 14.1 ± 3.09 | 1.06 ± 0.17 | 0.07 ± 0.01 |
3200 | 13.7 [2.34, 35.3] | 4367 | 5.00 ± 0.08 | 18.6 ± 3.37 | 1.47 ± 0.24 | 0.17 ± 0.02 |
3400 | 12.1 [0.98, 25.3] | 4067 | 4.89 ± 0.10 | 14.9 ± 2.16 | 1.32 ± 0.16 | 0.31 ± 0.02 |
Elevation (m) . | Temperature (°C) . | Precipitation . | pH . | C (mg/g) . | N (mg/g) . | P (mg/g) . |
---|---|---|---|---|---|---|
2700 | 17.0 [7.56, 40.1] | 4140 | 5.75 ± 0.09 | 14.1 ± 3.09 | 1.06 ± 0.17 | 0.07 ± 0.01 |
3200 | 13.7 [2.34, 35.3] | 4367 | 5.00 ± 0.08 | 18.6 ± 3.37 | 1.47 ± 0.24 | 0.17 ± 0.02 |
3400 | 12.1 [0.98, 25.3] | 4067 | 4.89 ± 0.10 | 14.9 ± 2.16 | 1.32 ± 0.16 | 0.31 ± 0.02 |
Notes: Climatic factors included temperature (°C) and precipitation and soil chemical properties included soil pH, organic carbon concentration (C, mg/g), total nitrogen concentration (N, mg/g) and total phosphorus concentration (P, mg/g). Shown values are mean, minimum and maximum for temperature, total precipitation events and mean ± standard error for soil chemical properties. Climatic factors were recorded during experimentation (June–October in 2015).
Elevation (m) . | Temperature (°C) . | Precipitation . | pH . | C (mg/g) . | N (mg/g) . | P (mg/g) . |
---|---|---|---|---|---|---|
2700 | 17.0 [7.56, 40.1] | 4140 | 5.75 ± 0.09 | 14.1 ± 3.09 | 1.06 ± 0.17 | 0.07 ± 0.01 |
3200 | 13.7 [2.34, 35.3] | 4367 | 5.00 ± 0.08 | 18.6 ± 3.37 | 1.47 ± 0.24 | 0.17 ± 0.02 |
3400 | 12.1 [0.98, 25.3] | 4067 | 4.89 ± 0.10 | 14.9 ± 2.16 | 1.32 ± 0.16 | 0.31 ± 0.02 |
Elevation (m) . | Temperature (°C) . | Precipitation . | pH . | C (mg/g) . | N (mg/g) . | P (mg/g) . |
---|---|---|---|---|---|---|
2700 | 17.0 [7.56, 40.1] | 4140 | 5.75 ± 0.09 | 14.1 ± 3.09 | 1.06 ± 0.17 | 0.07 ± 0.01 |
3200 | 13.7 [2.34, 35.3] | 4367 | 5.00 ± 0.08 | 18.6 ± 3.37 | 1.47 ± 0.24 | 0.17 ± 0.02 |
3400 | 12.1 [0.98, 25.3] | 4067 | 4.89 ± 0.10 | 14.9 ± 2.16 | 1.32 ± 0.16 | 0.31 ± 0.02 |
Notes: Climatic factors included temperature (°C) and precipitation and soil chemical properties included soil pH, organic carbon concentration (C, mg/g), total nitrogen concentration (N, mg/g) and total phosphorus concentration (P, mg/g). Shown values are mean, minimum and maximum for temperature, total precipitation events and mean ± standard error for soil chemical properties. Climatic factors were recorded during experimentation (June–October in 2015).
Ten top-ranked generalized linear mixed-effect models (GLMMs) results for biomass production as a function of biotic and abiotic factors and their interactions with plant species richness
Model . | k . | LL . | AICc . | ΔAICc . | wAICc . | . |
---|---|---|---|---|---|---|
S + FRic + pH + Tmax + FRic:pH:Tmax:S | 9 | −371.780 | 764.464 | 0.000 | 0.416 | 48.3% |
S + FRic + Tmax + FRic:S + FRic:Tmax:S | 9 | −372.169 | 765.241 | 0.777 | 0.282 | 48.0% |
S + FRic + pH + Tmax + pH:S + FRic:pH:Tmax:S | 10 | −371.775 | 767.156 | 2.692 | 0.108 | 48.2% |
S + FRic + Tmax + FRic:Tmax:S | 8 | −375.358 | 769.002 | 4.538 | 0.043 | 43.6% |
S + FRic + pH + Tmax + E + FRic:pH:S + FRic:E:S + pH:Tmax:S + FRic:pH:Tmax:S | 13 | −368.900 | 770.075 | 5.611 | 0.025 | 54.8% |
S + FRic + pH + Tmax + E + FRic:pH:S + pH:Tmax:S + Tmax:E:S + FRic:pH:Tmax:S | 13 | −368.903 | 770.081 | 5.617 | 0.025 | 54.7% |
S + FRic + pH + FRic:pH:S | 8 | −375.967 | 770.219 | 5.755 | 0.023 | 45.4% |
S + FRic + FRic:S | 7 | −377.334 | 770.419 | 5.955 | 0.021 | 38.8% |
S + FRic + pH + Tmax + E + FRic:pH:S + FRic:Tmax:S + pH:Tmax:S + Tmax:E:S + FRic:pH:Tmax:S | 14 | −368.903 | 773.174 | 8.710 | 0.005 | 54.7% |
S + FRic + pH + Tmax + E + FRic:pH:S + FRic:Tmax:S + Tmax:E:S + FRic:pH:Tmax:S + FRic:pH:E:S | 14 | −368.920 | 773.208 | 8.744 | 0.005 | 54.7% |
Model . | k . | LL . | AICc . | ΔAICc . | wAICc . | . |
---|---|---|---|---|---|---|
S + FRic + pH + Tmax + FRic:pH:Tmax:S | 9 | −371.780 | 764.464 | 0.000 | 0.416 | 48.3% |
S + FRic + Tmax + FRic:S + FRic:Tmax:S | 9 | −372.169 | 765.241 | 0.777 | 0.282 | 48.0% |
S + FRic + pH + Tmax + pH:S + FRic:pH:Tmax:S | 10 | −371.775 | 767.156 | 2.692 | 0.108 | 48.2% |
S + FRic + Tmax + FRic:Tmax:S | 8 | −375.358 | 769.002 | 4.538 | 0.043 | 43.6% |
S + FRic + pH + Tmax + E + FRic:pH:S + FRic:E:S + pH:Tmax:S + FRic:pH:Tmax:S | 13 | −368.900 | 770.075 | 5.611 | 0.025 | 54.8% |
S + FRic + pH + Tmax + E + FRic:pH:S + pH:Tmax:S + Tmax:E:S + FRic:pH:Tmax:S | 13 | −368.903 | 770.081 | 5.617 | 0.025 | 54.7% |
S + FRic + pH + FRic:pH:S | 8 | −375.967 | 770.219 | 5.755 | 0.023 | 45.4% |
S + FRic + FRic:S | 7 | −377.334 | 770.419 | 5.955 | 0.021 | 38.8% |
S + FRic + pH + Tmax + E + FRic:pH:S + FRic:Tmax:S + pH:Tmax:S + Tmax:E:S + FRic:pH:Tmax:S | 14 | −368.903 | 773.174 | 8.710 | 0.005 | 54.7% |
S + FRic + pH + Tmax + E + FRic:pH:S + FRic:Tmax:S + Tmax:E:S + FRic:pH:Tmax:S + FRic:pH:E:S | 14 | −368.920 | 773.208 | 8.744 | 0.005 | 54.7% |
Notes: Fixed effects are biotic factors (FRic, i.e. FRic calculated based on SLA) and abiotic factors including elevation (E), soil pH and maximum air temperature (Tmax) and their interactions with species richness (S). Random effects are experimental blocks nested in each elevational gradient. Shown are maximum log-likelihood (LL), the estimated number of model parameters (k), the information-theoretic Akaike’s information criterion corrected for small samples (AICc), the change in AICc relative to the top-ranked model (ΔAICc), AICc weighted (wAICc = model probability) and the marginal variance explained () as a measure of the model’s goodness-of-fit. The full model table is shown in Supplementary Table S1.
Ten top-ranked generalized linear mixed-effect models (GLMMs) results for biomass production as a function of biotic and abiotic factors and their interactions with plant species richness
Model . | k . | LL . | AICc . | ΔAICc . | wAICc . | . |
---|---|---|---|---|---|---|
S + FRic + pH + Tmax + FRic:pH:Tmax:S | 9 | −371.780 | 764.464 | 0.000 | 0.416 | 48.3% |
S + FRic + Tmax + FRic:S + FRic:Tmax:S | 9 | −372.169 | 765.241 | 0.777 | 0.282 | 48.0% |
S + FRic + pH + Tmax + pH:S + FRic:pH:Tmax:S | 10 | −371.775 | 767.156 | 2.692 | 0.108 | 48.2% |
S + FRic + Tmax + FRic:Tmax:S | 8 | −375.358 | 769.002 | 4.538 | 0.043 | 43.6% |
S + FRic + pH + Tmax + E + FRic:pH:S + FRic:E:S + pH:Tmax:S + FRic:pH:Tmax:S | 13 | −368.900 | 770.075 | 5.611 | 0.025 | 54.8% |
S + FRic + pH + Tmax + E + FRic:pH:S + pH:Tmax:S + Tmax:E:S + FRic:pH:Tmax:S | 13 | −368.903 | 770.081 | 5.617 | 0.025 | 54.7% |
S + FRic + pH + FRic:pH:S | 8 | −375.967 | 770.219 | 5.755 | 0.023 | 45.4% |
S + FRic + FRic:S | 7 | −377.334 | 770.419 | 5.955 | 0.021 | 38.8% |
S + FRic + pH + Tmax + E + FRic:pH:S + FRic:Tmax:S + pH:Tmax:S + Tmax:E:S + FRic:pH:Tmax:S | 14 | −368.903 | 773.174 | 8.710 | 0.005 | 54.7% |
S + FRic + pH + Tmax + E + FRic:pH:S + FRic:Tmax:S + Tmax:E:S + FRic:pH:Tmax:S + FRic:pH:E:S | 14 | −368.920 | 773.208 | 8.744 | 0.005 | 54.7% |
Model . | k . | LL . | AICc . | ΔAICc . | wAICc . | . |
---|---|---|---|---|---|---|
S + FRic + pH + Tmax + FRic:pH:Tmax:S | 9 | −371.780 | 764.464 | 0.000 | 0.416 | 48.3% |
S + FRic + Tmax + FRic:S + FRic:Tmax:S | 9 | −372.169 | 765.241 | 0.777 | 0.282 | 48.0% |
S + FRic + pH + Tmax + pH:S + FRic:pH:Tmax:S | 10 | −371.775 | 767.156 | 2.692 | 0.108 | 48.2% |
S + FRic + Tmax + FRic:Tmax:S | 8 | −375.358 | 769.002 | 4.538 | 0.043 | 43.6% |
S + FRic + pH + Tmax + E + FRic:pH:S + FRic:E:S + pH:Tmax:S + FRic:pH:Tmax:S | 13 | −368.900 | 770.075 | 5.611 | 0.025 | 54.8% |
S + FRic + pH + Tmax + E + FRic:pH:S + pH:Tmax:S + Tmax:E:S + FRic:pH:Tmax:S | 13 | −368.903 | 770.081 | 5.617 | 0.025 | 54.7% |
S + FRic + pH + FRic:pH:S | 8 | −375.967 | 770.219 | 5.755 | 0.023 | 45.4% |
S + FRic + FRic:S | 7 | −377.334 | 770.419 | 5.955 | 0.021 | 38.8% |
S + FRic + pH + Tmax + E + FRic:pH:S + FRic:Tmax:S + pH:Tmax:S + Tmax:E:S + FRic:pH:Tmax:S | 14 | −368.903 | 773.174 | 8.710 | 0.005 | 54.7% |
S + FRic + pH + Tmax + E + FRic:pH:S + FRic:Tmax:S + Tmax:E:S + FRic:pH:Tmax:S + FRic:pH:E:S | 14 | −368.920 | 773.208 | 8.744 | 0.005 | 54.7% |
Notes: Fixed effects are biotic factors (FRic, i.e. FRic calculated based on SLA) and abiotic factors including elevation (E), soil pH and maximum air temperature (Tmax) and their interactions with species richness (S). Random effects are experimental blocks nested in each elevational gradient. Shown are maximum log-likelihood (LL), the estimated number of model parameters (k), the information-theoretic Akaike’s information criterion corrected for small samples (AICc), the change in AICc relative to the top-ranked model (ΔAICc), AICc weighted (wAICc = model probability) and the marginal variance explained () as a measure of the model’s goodness-of-fit. The full model table is shown in Supplementary Table S1.

Effects of elevation, experimental warming and nitrogen addition on the relationship between species richness and productivity. The dashed lines represent the linear regression of productivity against species richness. The ribbons represent the confidence intervals of the linear regressions. Shown texts are the estimated intercepts, slopes, R squared as a measure of the goodness-of-fit and P values of the linear regressions.

Relationship between plant species diversity and biomass production (a) and the standard coefficients of the biotic and abiotic factors in the generalized linear mixed-effect models (GLMMs) of biomass production (b). Different point size represents the different values of FRic calculated based on SLA (cm2/g), different colors represent the different values of soil pH, and different shapes represent different maximum temperature (Tmax, °C). The dashed line represents the estimated mean value of biomass production as a function of species richness (S) and its quadratic term (i.e. Productivity ~ S + S2). The gray ribbon represents the confidence interval (95%) for the estimation. The standard coefficient plots show the estimated standard coefficients (i.e. solid points) and their confidence intervals (95%) in the generalized linear mixed-effect models of biomass production (detailed information see Materials and methods). Random effects are experimental blocks nested in each elevational gradient. Shown in the text is the marginal variance explained () as a measure of the model’s goodness-of-fit.
DISCUSSION
Compared with experimental warming and nitrogen addition, elevation showed the strongest influence on plant species diversity, suggesting that the elevation effect might dominate the distribution of plant species diversity under short-term climate change in mountain meadows. Experimental warming tended to decrease plant species diversity along the entire elevational gradient, suggesting that even short-term rises in temperature could have a detrimental influence on the plant species diversity of mountain meadows. This is consistent with previous studies showing that experimental warming can cause large and rapid species loss, especially for tundra and mountain grasslands ecosystems (e.g. Klein et al. 2004; Walker et al. 2006). Given that our top-ranked model for the interaction between elevation and warming was also significant, it seems plausible the effects of warming on plant species diversity can vary, even within the same type of ecosystem. Gruner et al. (2017) in a quantitative meta-analysis suggested that the negative effects of experimental warming on species richness in terrestrial ecosystems were often attributable to pronounced shifts in the functional groups of a community and were highly contingent on local climatic conditions and species composition. Our results partially agree with theirs because not only did experimental warming show little influence on functional diversity and soil chemical properties in our systems (results not shown) but the effect of experimental warming was only evident at the middle elevation.
It is not surprising that experimental warming showed little influence on aboveground biomass production. While previous empirical studies and meta-analyses indicate that warming generally increases aboveground plant biomass (Day et al. 2008; Lin et al. 2010), the response for plant productivity with high species richness may also be negative, due in part, to intense interspecific competition (De Boeck et al. 2007). This latter finding seems to be relevant to our results, where on the one hand, experimental warming decreased plant species diversity in our quadrats, while on the other hand, it increased biomass production for those plants remaining in the plots. Although nitrogen addition showed little effect on productivity in our study—a result consistent with its null effect on plant species diversity—there is increasing evidence showing that the impacts of nitrogen addition on plant biodiversity in mountain grasslands depend on the dose, application duration and climate (Humbert et al. 2016).
We found a nonlinear relationship between plant species diversity and biomass production in our systems, implying a shift in the direction of the relationship from positive to neutral, or even slightly negative, with increasing plant species diversity. This agrees with variable, and even conflicting, DPRs observed in natural ecosystems (Assaf et al. 2011; Flombaum and Sala 2008; Zhang et al. 2016). Partially in line with our expectations, elevation exhibited a strong influence on DPR in our systems and might be due to increasing elevation, which often represents an increasing level of physical stress that can enhance DPR through increasing facilitative interactions among species in a community (Mulder et al. 2001). While these results generally agree with previous studies, our study underscores the collective effects of biotic and abiotic factors in regulating DPR in natural ecosystems.
In our systems, functional diversity showed a substantial effect on DPR, which is consistent with the assumption of the important role of biotic factors in regulating DPR in natural ecosystems. However, FRic outperformed the CWM of single traits and trait dissimilarity among species in a community (e.g. Rao’s Q and FDis), which are often assumed to be more related to the selection and complementarity effects, respectively (Roscher et al. 2012), in explaining the variation in the DPR in our systems. Further, our results showed a substantial variation in the effect of FRic on DPR due to the different traits included in the analyses. Interestingly, SLA appeared to be the most pertinent trait to regulate DPR from the five traits measured in our systems. This is not surprising since SLA has been shown to be quite informative for local environmental variation probably due to its high plasticity (Dwyer et al. 2014; Wright et al. 2004). Moreover, our results showed that DPR tended to be positive when FRic is low, but negative when FRic is high. Since FRic was calculated using a single trait, lower FRic indicates the smaller absolute range of the character, implying species in a community were more similar to each other (Mason et al. 2005). Thus, our results suggest that competition exclusion might predominate the relationship between plant species diversity and productivity in our systems and that plant productivity was mainly determined by species with the particular values of SLA.
Results also revealed the strong influences of soil pH and maximum air temperature on DPR in our systems. High acid soil and lower maximum temperatures suggest severe physical stress on the community, which in turn, can cause more constraints on the expression of functional traits (Ottaviani et al. 2018) and increase the probability of positive interactions among species (Brooker et al. 2008; Callaway et al. 2002). More importantly, we found a strong interactive effect of FRic, soil pH and maximum air temperature on DPR in our systems. This result underscores the strong influence of elevation on the DPR in our systems, where elevation affected our focal biotic and abiotic factors, and their strong positive relationships. Interestingly, after accounting for their interactive effects on plant productivity, only FRic showed a positive relationship with plant productivity. This supports the robust relationship between functional diversity and plant productivity, which has been extensively reported in experimental studies (Cadotte 2017; Flynn et al. 2011; Liu et al. 2015; Roscher et al. 2012).
CONCLUSIONS
Climate change is projected to be detrimental to the maintenance of biodiversity on a global scale (Bellard et al. 2012). Our study underscores the vulnerability of local plant species diversity in alpine mountain meadows to climate warming, even in a short term. Despite the predominant role of elevation in regulating plant species diversity, productivity and their relationships in our systems, our results suggest that the interactions among biotic factors (i.e. FRic) and abiotic factors (i.e. soil pH and maximum air temperature) drive the shift in the relationship between plant species diversity and productivity from slightly negative to neutral and then to positive with elevation. We suggest that future analyses of biotic factors that are coupled with abiotic sources of variation can provide important insights into the DPRs in natural ecosystems.
Supplementary Material
Supplementary material is available at Journal of Plant Ecology online.
Table S1: Generalized linear mixed-effect models (GLMMs) results for biomass production as a function of biotic and abiotic factors and their interactions with plant species richness.
Figure S1: Boxplot for the air temperature recorded in the control and open-top chambers (OTC) along the elevational gradient.
Figure S2: Barplot of information-theoretic evidence ratio (ER) comparing two generalized linear mixed-effect models (GLMMs) for experimental treatments.
Figure S3: Barplot of information-theoretic evidence ratio (ER) comparing two generalized linear mixed-effect models (GLMMs) for functional diversity.
Figure S4: Barplot of information-theoretic evidence ratio (ER) comparing two generalized linear mixed-effect models (GLMMs) for abiotic factors.
Figure S5: Correlation matrix plot for biotic and abiotic factors.
Funding
This study was financially supported by the Strategic Priority Research Program of the Chinese Academy of Sciences (XDB31000000) and the National Natural Science Foundation of China (31500335).
Acknowledgements
We thank Associate Editor and two anonymous reviewers for their helpful comments and suggestions. This work was done in the Lijiang Alpine Botanical Garden of the Kunming Institute of Botany, Chinese Academy of Sciences. We thank Lian-Ming Gao, Kun Xu and De-Tuan Liu for the help with the experiment.
Conflict of interest statement. The authors declare that they have no conflict of interest.
Authors’ Contributions
J.-J.L. and X.-J.G. conceived and designed the experiments; J.-J.L. performed the experiments; J.-J.L., Y.X. and Y.-X.S. analyzed the data; J.-J.L. and K.S.B. wrote the manuscript, with significant contributions from X.-J.G.
Data Availability
The essential data and R script for reproducing the data analyses of this paper are deposited in https://doi.org/10.6084/m9.figshare.14707575.v1.