Introduction

Biodiversity conservation is one of the most important issues that closely link human welfare and ecosystem functions in a rapidly changing world1,2,3,4. The current studies on the relationship between biodiversity and ecosystem functions are increasingly moving their focus from species diversity indices to trait-based approaches5,6,7,8,9,10. Interactions between biodiversity and ecosystem functions are argued not necessary to be linear or unimodal across different scales. Therefore, the relationship between biodiversity and ecosystem functionality calls for causal network analyses rather than simply bivariate analyses to understand mechanisms. However, it is still unclear and under debates how biodiversity changes affect community productivity on alpine grasslands.

On the Tibetan Plateau, alpine grasslands are sensitive and vulnerable to both natural and anthropogenic disturbances11. Degradation of alpine grasslands is mainly attributed to climate warming and overgrazing there12,13,14,15,16. For example, Piao et al.17 found that the spatial patterns of vegetation green-up and its change in response to warming are altitude dependent. As these pastures are an important ecological security by protecting the headwaters of major rivers in Asia, their ongoing degradation likely threatens the livelihood of local residents and their distinctively nomadic culture11,18, but may also have wider impacts on water security in East China and South-Asia. Assumed to effective to recover degraded pastures, metal fences were built on this plateau to control grazing by livestock12,18. However, the magnitude and direction of fencing effects on plant communities and soil properties are still under debates19,20,21,22,23,24,25. The questions whether a general biodiversity-productivity relationship exists for Tibetan alpine grasslands and how this relationship is affected by climatic factors are still under discussion.

A positive linear relationship between species richness and productivity has been reported for Tibetan alpine grasslands26,27. However, no relation was found between the corresponding residuals of species richness and productivity after environmental and climatic influences were removed. Ma, et al.26 and Wang, et al.27 argued that plant species richness can weakly affect community productivity across Tibetan alpine grasslands. In fact, intrinsic community properties, for instance, community assembly of plant functional groups, can be as important as extrinsic abiotic divers in explaining spatial and temporal variabilities of grassland productivity across the Tibetan Plateau19,28,29,30. Current findings are mainly deducted from bivariate analyses, and insight studies on direct and indirect interactions between potential biotic and abiotic are still rare.

Plant functional trait diversity is useful to address questions of plant performance, community assembly and ecosystem functions along climate, resources and disturbance gradients10,31,32,33,34. For example, changes in plant trait composition in response to grazing vs. mowing disturbance can be distinguished and predicted by functional trait means and divergence35,36. Maintaining and enhancing functional trait diversity might be vital to buffer the negative effects of climate change and human activities on ecosystem multifunctionality34,37, especially for dry and alpine grasslands on our planet. He, et al.38 explored the general bivariate relationship between plant leaf traits of 74 species sampled from 49 sites on the Qinghai-Tibetan Plateau, but they did not link plant functional traits with ecosystem functionality.

Precipitation use efficiency (PUE) is a widely used proxy for the sensitivity of grassland productivity in response to precipitation39,40,41. Hu, et al.42 and Yang, et al.43 examined the spatial PUE pattern across different alpine grassland types, explained the potential influences of climatic variables and soil properties, but did not examine the direct and indirect influences of plant functional trait diversity. Current studies that follow trait-based approaches are limitedly conducted in the humid meadows in the eastern Tibetan plateau28,44,45,46 and rarely in the arid steppe and desert-steppe zones in the northwestern region47. At large spatial scales, from both field observations and remote sensing modelling, scientists have confirmed that Tibetan alpine grasslands are very sensitive to changes in precipitation28,29,48,49,50. However, gaps in understanding the interactions between biodiversity, functionality and precipitation changes still remain.

In this paper, we concretely aim (1) to examine differences in biodiversity and functionality among zonal alpine grassland types, especially focusing on for functional trait diversity indices; (2) to clarify the spatial patterns of biodiversity and functionality along geographical, edaphic and climatic gradients at a regional scale; (3) to figure out direct and indirect links between biodiversity and functionality in response to the regional water availability gradient, an indicator of habitat moisture. In addition, we also want to find out whether plant functional trait diversity indices are more effective than the traditional plant species diversity (PSD) indices, such as richness, dominance and evenness, in explaining the sensitivity of alpine grassland productivity to regional climatic gradients across the Tibetan Plateau.

Methods

Study area and sampling

Study area

We sampled 75 plots at 15 sites in total during the peak-growing season of 2014, with five sites for each grassland type, westwardly across humid alpine meadow (AM) dominated by Kobresia pygmaea, semi-arid alpine steppe (AS) dominated by Stipa purpurea, and arid alpine desert-steppe (ADS) co-dominated by S. purpurea and S. glareosa on the Northern Tibetan Plateau (29°53′–36°32′N; 78°41′–92°16′E) (Fig. 1 and Supplementary Table S1). Plant growing season (GSP) there generally starts in early May and ends in late September when up to 85% of annual precipitation happens and mean daily temperature is always over 5.0 °C51. To minimize disturbances by domestic livestock and wild herbivores, sampling was conducted at fenced pastures. Thus, the peaking aboveground biomass can be accepted as aboveground net primary productivity (ANPP).

Figure 1
figure 1

Study region, climate and sites sampled across the Northern Tibetan Plateau.

(a)-regional context of alpine grassland types (AGTs) showed site locations. (b)-mean annual temperature (MAT) and (c)-mean annual precipitation (MAP) over the Tibetan Autonomous Region, China, indicating the climatic conditions of our sampling sites. ArcGIS10.2 was used to create these maps. Climatic data ranges from 1979 to 2008 and daily weather records are available from the China Meteorological Data Sharing Service System.

Full size image

Sampling design

At each site, we established a sub-enclosure of 200 m × 200 m in size for measuring species composition, leaf traits and aboveground biomass. Five 50 cm × 50 cm plots were independently arranged at 20-m intervals along a random transect line. We identified plant species occurring, measured general leaf height (GLH, cm), and estimated percentage cover for each species within each plot. Aboveground biomass was harvested, sorted by species and stored in paper bags. Plant samples were dried at 65 °C for 48 h and weighted to 0.001 g so that rare species could also be identified. Species frequency was additionally sampled using thirty 0.1 m2 random circles at each site. We also measured specific leaf area (SLA, cm2g−1) and leaf mass fraction (LMF, %) as plant functional traits52.

Data management

Species and trait diversity indices

We firstly calculated species relative dominance based on percentage cover, leaf height and frequency measured at each site as did Wu, et al.22. Then we defined species richness (SR) as the total number of species at each quadrat (Plot_SR). Next, we calculated Shannon index of diversity (Shannon), Simpson index of diversity (Simpson), and Pielou evenness index (Pielou)53. Plot-SR, Shannon, Simpson and Pielou were grouped as plant species diversity (PSD) indices.

To compare community trait composition of GLH, SLA and LMF, we calculated community weighted mean (CWM) values for each quadrat following the equation of Violle, et al.10 and Ricotta and Moretti54.

where the Pij is the relative dominance (cover percent) of the species i in the quadrat j, Tij is the mean trait value of the species i in the quadrat j, and CWMj is the community-weighted trait of the quadrat j.

We also followed the distance-based framework of Laliberte and Legendre53 to quantify the variance of community trait distribution. Functional trait divergence (FTD) that reflects niche complementarity (overlaps/divergence) was calculated for plant GLH, SLA and LMF, respectively.

where the Pij is the relative dominance of the species i in the quadrat j, Tij is the mean trait value of the species i in the quadrat j, and FTDj is the community-weighted trait of the quadrat j.

Climatic variables and topsoil nutrients

We downloaded daily records of precipitation and temperature from the China Meteorological Data Sharing Service System (http://cdc.nmic.cn/home.do). As low-temperature stress is common over the Tibetan Plateau, plants generally start to grow when it is warm. In addition to GSP, accumulated temperature when daily mean is over 5 °C (AccT) also accumulated27. GSP and AccT in 2014 were firstly compiled from daily climatic raster surfaces that were interpolated using ANUSLPIN 4.355, and then extracted for each site in ArcGIS 10.2 (ERSI, Redlands, CA, USA). Finally, we accepted the ratio of GSP/AccT to describe the habitat moisture index (HMI = GSP/AccT) as did Wang, et al.27. Soil variables, including soil organic carbon (SOC) and total nitrogen (STN) were extracted from Wu, et al.51 and Li, et al.56.

Statistical analyses

Geospatial responses of alpine grasslands to regional abiotic gradients

As plant communities are not necessary to linearly respond to environmental changes57, generalized additive models (GAM) were induced to cope with problems such as collinearity and non-linearity between responsible and explanatory variables58 using the mgcv package in R59. Here, abiotic explanatory variables include site geographical locations (longitude, latitude, and altitude), climatic factors (GSP, AccT and HMI) and soil properties (SOC, STN and C: N ratio). The response variables were ANPP and PUE at the plot level. We followed a backward selection approach with Akaike Information Criterion (AIC) and Bayesian Information Criteria (BIC) to find the optimal models out58,60 (Supplementary Figures S1 and S2 and Tables S2 and S3).

Functional and structural differences between alpine grassland types

Functional and structural differences among zonal alpine grassland types were initially examined by one-way analysis of variance (ANOVA) with Tukey HSD test (P < 0.05) (Supplementary Figure S3). Next, we conducted principal component analyses (PCA) to examine functional and structural differences of communities surveyed, firstly did separately for PSD, CWM and FTD, and then pooled them together. We used the PCA coordinates with an eigenvalue >1 to identify the plant strategy along which trait or species diversity indices co-vary across communities.

Direct and indirect effects of water availability on ecosystem function

For the biodiversity-functionality relationship, the causal network analysis is likely more effective to clarify direct and indirect interrelations than bivariate regressions or analyses of variances59. Plant functional trait diversity may provide more mechanistic understandings than species diversity indices when relates to ecosystem functionality10. We conducted structural equation model (SEM) using the lavaan package61 in Rstudio59 to explore interrelations among climate, soil, community structure, and ecosystem function (sensitivity) as hypothesized in Fig. 2.

Figure 2
figure 2

Hypothetical interrelations among climate, soil, community structure, and ecosystem function (sensitivity) of alpine grasslands on the Northern Tibetan Plateau.

Red lines are for direct impacts, blue lines for indirect impacts, and green lines for covariations. Habitat moisture index (HMI) equals the ratio of growing season precipitation (GSP) to accumulated temperature when daily mean is over 5 °C (AccT). Soil C: N ratio (SCNR) was also considered in this study. Plant species diversity (PSD) includes the indices of richness, Shannon, Simpson, Pielou evenness at the plot level. Community weighted means (CWM) and functional trait divergence (FTD) were calculated from general leaf height (GLH, cm), specific leaf area (SLA, cm2 g−1) and leaf mass fraction (LMF, %).

Full size image

Results

Geospatial patterns of ANPP and PUE across studied communities

Due to collinearity between explanatory variables (Supplementary Figures S1 and S2), soil parameters were excluded from the optimal models (Supplementary Tables S2 and S3). The smoothers of HMI and site altitude were significant, entered the optimal models, and together explained 80.3% of the variation in ANPP and 68.8% of variation in PUE (Table 1), respectively. The patterns of ANPP and PUE are similar, showing clear increasing trends along the HMI gradient from approximately 0.18 mm °C−1 to the most humid sites 0.4 mm °C−1, and unimodal patterns at the dry end of the HMI gradient, lower than 0.18 mm °C−1 (Fig. 3a,c). Along the site altitudinal gradient, both ANPP and PUE show unimodal patterns with peaks at approximately 4800 m (Fig. 3b,d).

Table 1 The summary of the optimal generalized additive models (GAMs) with habitat moisture index (HMI) and site altitude (Alt.) as explanatory variables.
Full size table
Figure 3
figure 3

The estimated smoothers for ANPP (a,b, solid lines) and PUE (c,d, dashed lines), respectively, obtained by the optimal generalized additive models (GAMs) that include habitat moisture index (HMI, a,c) and site altitude (b,d) as explanatory variables. The grey areas show the 95% point-wise confidence bands.

Full size image

Grassland functional structure based on plant species and leaf traits

The biodiversity indices that were separately grouped into PSD, CWM and FTD, showed different explanatory powers for variance of community structure in this study. The PCA results indicated that the explanatory power of PSD is comparable to the FTD group, the former accounting for 89.68% (Fig. 4a) vs. the latter for 94.60% (Fig. 4c) of the total variance explained by the first component. The CWM group can account for 78.66% of the total variance, with the first component (43.04%) highly correlated with LMF and the second component (35.62%) highly correlated with GLH and SLA (Fig. 4b). When all community structural predictors were pooled together, the first PCA components (63.23% of the variance) separated communities according to the PSD and FTD values (Fig. 4d, with PSD being negatively correlated to FTD, to see Supplementary Figure S1) while the second PCA component (11.38% of the variance) discriminated communities according to the CWM indices.

Figure 4
figure 4

Principal component analyses (PCA).

(a)-based on plant species diversity (PSD) that includes species richness at the plot level (Plot_SR), Shannon diversity index (Shannon), Simpson dominance index (Simpson) and Pielou evenness index (Pielou); (b)-based on community weight means (CWM) of general leaf height (CWM_GLH), leaf mass fraction (CWM_LMF) and specific leaf area (CWM_SLA); (c)-based on functional trait divergence (FTD) of general leaf height (FTD_GLH), leaf mass fraction (FTD_LMF) and specific leaf area (FTD_SLA); and (d)-with all potential explanatory variables being pooled together. Yellow diamonds, blue triangles and green squares represent alpine meadow (AM), steppe (AS) and desert-steppe (ADS), respectively. For each component we indicate the percentage of variance explained. See Supplementary Figure S1 for correlations between PSD, CWM and FTD values, Figure S3 for site locations within each alpine grassland type (AGTs) and climatic background, and Figure S4 for variable comparisons among the three AGTs.

Full size image

The causal networks for ANPP and PUE across studied communities

The influence of HMI on ANPP was mainly mediated through PSD and CWM of plant functional traits. HMI directly (β = 0.51, standardized coefficient), and indirectly through CWM (β = −0.24, standardized coefficient) and PSD (β = 0.24, standardized coefficient), impacted ANPP (Fig. 5a). The indirect influences of FTD of plant functional traits and soil properties (SCNR) on ANPP were not significant (Supplementary Table S4). SCNR was controlled by HMI (β = 0.52, standardized coefficient) while FTD was mainly determined by PSD (β = −0.69, standardized coefficient). The influence of HMI on PUE was significantly mediated by PSD (β = 0.32, standardized coefficient) and CWM (β = −0.23, standardized coefficient) (Fig. 5b). The indirect influence of HMI through PSD on PUE was approximately equal in magnitude to the direct influence of HMI (β = 0.34, standardized coefficient). The FTD of plant traits and soil properties (SCNR) were no significant impacts on PUE (Supplementary Table S4). The strongest relationship observed in the SEM analyses was between PSD and FTD (β = −0.69, standardized coefficient).

Figure 5
figure 5

Structural equation models examining effects of environmental factors and biodiversity components on ANPP (a) and PUE (b), respectively, at the regional scale across zonal alpine grassland types on the Northern Tibetan Plateau. Green and red arrows indicate significant positive and negative effects, respectively. Line width illustrates path strength. Values associated with path arrows represent standardized path coefficients. Values near to the round arrows show the residuals. Non-significant paths and variables were eliminated from the final SEMs. Results of model fitting: (a) χ2 = 1.513, P = 0.679; (b) χ2 = 1.513, P = 0.679 (Note: higher P-values associated with χ2 tests indicate better fitting to data. See Supplementary Table S3 for summaries of the full SEMs with non-significant relationships included for ANPP and PUE, respectively.

Full size image

Discussion

Geospatial nonlinear patterns of productivity and collinearity between abiotic factors

Alpine grasslands are generally believed to be sensitive to climate change. Unfortunately, the question how community structure under ongoing climate change affects ecosystem functionality is still poorly understood. To mechanistically understand this question, we have examined the spatial patterns of productivity-biodiversity relationship within an effect-response framework across three zonal alpine grasslands types on the Northern Tibetan Plateau. Our results have shown that, at the regional scale, both ANPP and PUE (termed as the ratio of ANPP to GSP) nonlinearly respond to the increasing water availability (HMI) gradient at the regional scale (Fig. 3). Importantly, the unimodal ANPP-altitude relationship at such large geospatial scale was consistent with the finding of Wang, et al.27 at a local scale, along a mountain slope in the central Tibetan Plateau. Therefore, our results supported hypothesis that ecological processes are not necessary to linearly respond to climate change either62,63,64.

The responses of alpine grassland communities subject to abiotic and biotic factors have been explored by a few studies based on bi-variate analyses29,42,43. But the collinearity between explanatory variables and the nonlinear response patterns have rarely documented. For example, Yang, et al.43 reported a unimodal pattern of PUE along mean annual precipitation (MAP) across alpine steppe and meadow communities on the Tibetan Plateau, and attributed it to the linear bivariate relationships of PUE with species richness, soil organic carbon and soil silt content. As soil nutrients were highly correlated with HMI (Supplementary Figure S2), in this study soil explanatory variables were excluded from the optimal GAMs. Soil C: N was also confirmed to have no significant impacts on either ANPP or PUE in the full SEMs (Supplementary Figure S4 and Table S4). This may be attributed to the historical co-evolution of climate, vegetation and soils, because soil development in alpine biomes and mountain regions is mainly controlled by climatic condition there65. Therefore, our findings on the collinearity between soil nutrients and climatic factors implied that the dynamics of community structure and ecosystem functionality might be weakly influenced by soil nutrients. This can also explain why no relation exists between corresponding residuals of species richness and productivity after removing environmental influences26,27.

Vegetation at high altitude is generally believed to be more sensitive to global changes than other biomes because of the extreme habitat conditions65,66, for instance, stress due to low temperatures and poor soil nutrient availability affect plant performance and reproduction. In addition, we found the optimal GAMs that included HMI and site altitude (adj. R2 = 0.803, AIC = 559.77 for ANPP; adj. R2 = 0.688 AIC = −327.74 for PUE) were a bit better that the sub-optimal models that included GSP and site altitude (adj. R2 = 0.798, AIC = 561.71 for ANPP; adj. R2 = 0.676, AIC = −325.40 for PUE) (Supplementary Tables S2 and S3), although a few studies have argued that precipitation tends to be more important than temperature for alpine grassland productivity29,67. Our results support that alpine grassland dynamics closely related the wet/dry condition and influenced by the combination of temperature and precipitation27,50,68. The combination of temperature and precipitation varies along an increasing altitudinal gradient to affect community assembly and ecosystem productivity27. Therefore, the unimodal patterns of ANPP and PUE at the dry end of HMI gradient in this study might be modified by site altitude differences, because we have examined unimodal patterns of ANPP and PUE against site altitude at the regional scale. Therefore, our results indicate that elevation dependency of climate changes, especially in the combination of temperature and precipitation, should be also clarified for predicting vegetation dynamic in montane regions27,69.

Ecosystem causal networks and explanations for ANPP and PUE from biotic factors

Previous studies showed that drier places tend to have lower ANPP because of sparser community, higher evaporation potential, and higher tolerance to drought stress39,41,70,71. In general, we found that PUE showed an increasing trend along the HMI gradient from approximately 0.18 mm °C−1 to the most humid sites 0.4 mm °C−1 (Fig. 3c), being partly consistent with the positive linear relationship between PUE and mean annual precipitation (MAP) reported by Hu, et al.42. However, a unimodal pattern between PUE and HMI was found at places with HMI lower than 0.18 mm °C−1 (Fig. 3c), being partly consistent with findings of Yang, et al.43. In addition to the influence of site altitudes, the shifts of ANPP and PUE patterns, from unimodal to approximately linear, may also be attributed to changes in species assembly and trait divergence. Community assembly of plant functional groups have been confirmed by bivariate regressions to be as important as precipitation in controlling and explaining the spatial variation in ANPP28. In this study, we also found high correlations of ANPP (PUE) with PSD and FTD indices (Supplementary Figure S1). In addition, comparisons of FTD among alpine grassland types showed that plant functional traits were more convergent with lower averages and smaller variations in humid meadows than semi-arid steppes and arid desert-steppes (Supplementary Figure S3). This may also contribute to the unimodal patterns of ANPP and PUE at the end of HMI.

Interior changes of community assembly, including species composition, plant functional groups and traits, can consequently regulate ecosystem function and sensitivity in response to exterior driver for predicting potential responses to regional environmental changes40,72,73,74,75,76. Surprisingly, we found that PSD had similar importance as FTD in explaining the variance of alpine grassland communities from the PCA results (Fig. 4a,c) and that PSD had a significant impact negative impact on FTD with the highest standardized coefficient (β = −0.69) in the final SEMs (Fig. 5). This may be explained by the environmental filtering theory that community assembly is not a random process but based on the filtering of environmental factors on plant functional traits31,77,78. The indirect impact of HMI through PSD was thereby found as strong as the indirect impact through CWM (having same standardized coefficients, Fig. 5a). FTD had no significant effect on ANPP or PUE due to covariations with PSD in the SEM analyses (Fig. 5 and Table 2). These findings indicate that alpine plants and communities have been adapted to local abiotic conditions because of environmental filtering effects on plant trait evolution. It should be clarified that both CWM and FTD, in addition to plant trait values, are also determined by species relative dominance as shown in their definitions. The productivity sensitivity to climate changes at the community level can also be modified by changes in plant relative dominance. Therefore, we suggest that the flexibility of plant functional traits within and across species as well as their dominance changes in response to climate change and human disturbance should be seriously considered in long-term and large-scale studies on ecosystem functionality.

Table 2 Summaries of the optimal structure equation models (SEMs) for aboveground net primary productivity (ANPP) and precipitation use efficiency (PUE), respectively.
Full size table

In this study, we explored the patterns of alpine grassland productivity (ANPP) and its sensitivity (PUE) in response to regional climate changes, and examined both strength and direction of biotic and abiotic drivers within the effect-response framework. First, the nonlinearity of ANPP and PUE in response to abiotic variables was detected by general additive models, in which HMI and altitude entered and together accounted for 80.3% and 68.8% of the variation in ANPP and PUE respectively. This result highlights the combinational influences of climate warming and precipitaoitin fluctuation on alpine grassland dynamics. Because the combination of temperature and precipitation always varies along increasing altitudinal gradients, the elevation-dependency should also pay more attention in modelling alpine grassland responses to climate change in this plateau and other mountainous regions. Second, both principal component analyses and structural equation model confirmed that functional trait diversity, including community weighted means and functional trait divergence, is as important as plant species diversity for explaining the nonlinear relationship between climate changes and ecosystem productivity in alpine grasslands on the Northern Tibetan Plateau. These findings are likely due to the collinearitiy of abiotic explanatory factors with climate gradients. On the other hand, due to the environmental filtering effects, plants have been adapted to alpine habitat conditions with highly differential functional traits. Trait flexibility within and across species and changes in plant relative dominance should be examined for long-term studies on the biodiversity-ecosystem functionality relationship.

Additional Information

How to cite this article: Wu, J. et al. Plant functional trait diversity regulates the nonlinear response of productivity to regional climate change in Tibetan alpine grasslands. Sci. Rep. 6, 35649; doi: 10.1038/srep35649 (2016).