Mounting evidence is starting to demonstrate the ability of ecosystems to simultaneously provide multiple functions and services, i.e. multifunctionality. The concept of multifunctionality has been defined in several ways, for instance, “the overall functioning of an ecosystem” (Hector and Bagchi, 2007), “the simultaneous provision of several ecosystem processes” (Byrnes et al., 2014), the “provision of multiple ecosystem functions and services as high or desired levels” (Allan et al., 2015), and “the potential of landscapes to supply multiple benefits to society” (Mastrangelo et al., 2014). However, such definitions are still unable to give a clear and general connotation, denotation and measurement of multifunctionality. Considering the differences among fundamental and applied research questions, a general framework was proposed which defines multifunctionality at two levels: ecosystem function multifunctionality and ecosystem service multifunctionality (Manning et al., 2018). The first one focuses on biological, geochemical and physical processes within an ecosystem, which is most relevant to fundamental research (Manning et al., 2018). The second one focuses on the simultaneous supply of multiple ecosystem services relative to human demands, and emphasizes applied research under the considerations of definable management objectives of stakeholders. Holistic understanding of “whole ecosystem” properties and measures is timely and feasible based on large datasets and model outputs with comprehensive knowledge from multidisciplinary. Knowledge of multifunctionality can avoid one-sided evaluation of ecosystems using only individual functions. Moreover, multifunctionality sets up linkages between the natural capacity of ecosystems with the requirements of stakeholders. Identifying thresholds of multifunctionality, along with climatic conditions and anthropogenic perturbation, can provide us with signals for degradation and catastrophic transition of ecosystems.

Different methods have been developed to quantify multifunctionality. The simplest is the averaging approach (Hooper and Vitousek, 1998), which aims to collapse multifunctionality into a single metric that estimates the average value of multiple functions achieved in a given assemblage or plot. The single threshold method (Gamfeldt et al., 2008; Zavaleta et al., 2010) tallies the number of functions that quantitatively exceed some pre-defined threshold of “functionality” in certain ecosystems. Recently a multiple threshold method was introduced (Byrnes et al., 2014), which has been widely used in studies of ecosystem multifunctionality. The various thresholds were defined based on the contributions of species diversity to different functions in the ecosystems. For example, minimum threshold (T_min): The lowest threshold where diversity begins to have an effect; maximum threshold (T_max): The value of the threshold beyond which the slope first declines such that it is not significantly different from zero; Threshold of maximum diversity effect (T_mde): The value of the threshold where diversity has its strongest positive or negative effect; And realized maximum effect of diversity (R_mde): The strength of the relationship where diversity has its strongest positive and/or negative effects, which indicates the maximum observed effect size of species richness on the number of functions surpassing the threshold. Such methods provide useful insights into holistic ecosystem functions. However, each has advantages and disadvantages.

As the third pole, the Tibetan Plateau is characterized as the highest mountain in the world with an average altitude above 4000 m a.s.l. (Piao et al., 2011; Li et al., 2017). Grasslands are the predominant vegetation types on the plateau, which cover over 35% of the total mountain area. The grassland ecosystems on the plateau play an important role in maintaining the livelihoods of humans and wildlife. Climatic conditions and geomorphic features on the plateau are highly heterogeneous along altitudinal gradients. Various grasslands types, including alpine meadow, alpine steppe, temperate meadow, temperate steppe, and temperate desert steppe, have developed on the mountains under the control of diverse climatic conditions (Editorial Committee for Vegetation Map of China, 2001). Temperate steppe and temperate meadow are the dominant grassland types when altitude is lower than 3300 m, while alpine steppe and alpine meadow dominate when altitude is above 3300 m. Grassland degradation is a common phenomenon on the high-elevation and low-temperature Tibetan Plateau, and the degraded grasslands amount to 5.0 × 10⁵ km² (Harris, 2010). The primary factors driving the degradation succession are herbivorous overgrazing together with the severe-climatic conditions (Harris, 2010). Livestock overgrazing generally suppresses edible grasses, for example Stipa spp. and Festuca spp., which induces a shift to dominance by sedge species, i.e. K. humilis or K. pygmaea. Meanwhile pikas, feeding on plant roots, prefer to establish their caves in K. pygmaea grassland. Thus, overgrazing from livestock and pikas together with the freeze-thaw cycles of soil are the main factors leading to degradation of K. pygmaea grasslands.

Numerous studies have been conducted which quantify individual functions of the alpine grassland, construct the relationships of individual functions with diversity, and test the effects of climate and grazing on individual functions. However, little is known about the spatial patterns of multifunctionality of the alpine grasslands on the Tibetan Plateau, or how climate and livestock grazing drive such patterns of multifunctionality. Recently, multifunctionality was quantified along a whole degradation succession series in an alpine meadow using three key individual functions out of a total of 14 functions (Song et al., 2020). Moreover, the effects of species richness and fragmentation intensity on multifunctionality were examined, and a threshold was identified to indicate the shift of multifunctionality from slow decrease to rapid collapse (Song et al., 2020). A critical next step is to assess multifunctionality of the alpine grassland on the Tibetan Plateau. Here we quantify the multifunctionality of the alpine grassland using data from remote sensing NDVI combined with data collected from transects in Qinghai Province and Tibet autonomous region which cover the major types of grassland on the Tibetan Plateau. We aim to: 1) map the spatial patterns of multifunctionality using the four key individual functions, i.e. above and belowground biomass, soil organic carbon storage (SOC), and species pool size (species richness); 2) test the relationships of the multifunctionality with climate and livestock grazing intensity; and 3) identify the threshold of the multifunctionality transition from slow decrease to rapid collapse.

The Tibetan Plateau lies in the west of China, and it serves as one of the most important ecological security barriers in China and plays a significance role in the conservation of global biodiversity, as well as the regulation of climate and water cycles (Yao et al., 2012). This region is characterized as an arid alpine climate with an average elevation over 4000 m above sea level. The annual average temperature ranges from -15 to 10 ℃ (You et al., 2013). The precipitation presents high spatial heterogeneity. Annual average precipitation is about 50-150 mm in the northwest, and 300-450 mm in the southeast of the Tibetan Plateau. Grassland is the most widely distributed vegetation type, occupying about 75% of the total area. Grassland types, such as humid alpine meadow, semi-arid alpine steppe, and arid alpine desert-steppe are distributed along the humid gradient from the southeast to the northwest of the Tibetan Plateau.

An extensive grassland investigation was carried out in Qinghai Province and Tibet autonomous region, which cover the major types of grassland on the Tibetan Plateau (Fig. 1). In Qinghai Province, the 35 experimental sites were mainly distributed in typical alpine grassland. We investigated 255 plots from 35 counties during the plant growing seasons from 2013 to 2017. In the Tibet autonomous region, a 1200-km transect was established in spring 2009 on the northern Tibetan Plateau which includes mainly four types of alpine grassland, i.e. alpine meadow, alpine steppe, temperate steppe, and temperate desert steppe. Nineteen sites were arranged along the transect. Dominant plant species vary from Kobresia pygmaea in meadows, and Stipa purpurea in steppes, to S. purpurea and S. glareosa in desert-steppes along the transect from east to west (Wu et al., 2014). In-situ observations were carried out annually and data were collected from 252 plots during 2013 to 2017.

At each site, three or five 50 cm × 50 cm plots were located at 20 m intervals within each of the two flat areas of 200 m × 200 m. Within each plot, plant species were recorded and their cover proportions were visually estimated. Plant materials were harvested with scissors at the soil surface and stored in separate envelopes by species. Within each plot, one soil block of 25 cm (length) × 25 cm (width) × 20 cm (depth) was sampled for root biomass estimation, because 85% of living roots are concentrated in that soil layer in alpine grasslands on the Tibetan Plateau (Deléglise et al., 2011). Soil samples were collected using a soil corer of 7-cm diameter.

In the laboratory, root samples were washed and cleaned using running water. The living roots were sorted by hand according to the color, elasticity, and attachment of fine roots. Finally, all the plant samples, above and below ground, were oven-dried at 65 ℃ for 48 h and weighed.

Soil moisture was measured gravimetrically and soil bulk density was determined. Subsamples used for C and N analysis were air-dried, sieved (2 mm mesh), handpicked to remove fine roots, and then ground in a ball mill. The concentration of total SOC was analyzed using the modified Walkley-Black method by boiling a soil-dichromate-sulfuric acid mixture, with a 1:20 ratio of soil to extractant, at 175 ℃ for 5 min (Walkley and Black, 1934; Lettens et al., 2005). A conversion factor of 1.1 was applied to correct the SOC measurements obtained from the modified Walkley-Black method (Lettens et al., 2005). Organic matter content values were corrected for water content and expressed on a dry weight basis.

Four individual functions were selected from the 14 individual functions in our previous study in an alpine meadow after considering correlations among the multiple functions and independence (Song et al., 2020). The four individual functions, i.e. above and belowground biomass, soil organic carbon storage (SOC), and species pool size (species richness), play important roles in grassland ecosystems and can be used to represent key aspects of their functions. In this study, we used these four individual functions to quantify the multifunctionality of the grassland on the Tibetan Plateau.

For aboveground biomass of grassland on the Tibetan Plateau, we adopted a model which had been established and used in a previous study (Cao et al., 2019) that is based on the relationship between remote sensing vegetation index and measured aboveground biomass in the field (AGB) in different types of grassland on the Tibetan Plateau. Data of normalized difference vegetation index (NDVI) version 6 were downloaded from the monthly MOD13A3 product of the Moderate Resolution Imaging Spectroradiometer with 1 km spatial resolution (MODIS) (https://lpdaac.usgs.gov/get_data/data_pool). The monthly NDVI data were developed using the Maximum Value Composition method (MVC), and had been calibrated for geometrical effects, atmospheric effects, and cloud contamination. Annual NDVI_max was calculated using the maximum value composite, which is widely used to estimate the annual maximum biomass of grasslands (Holben, 1986).

Soil organic carbon storage data (SOC) were available from the Land-Atmosphere Interaction Research Group (http://globalchange.bnu.edu.cn/research/soil2). Prior to the analyses, SOC data were resampled to the same spatial resolution of 1 km × 1 km. SOC data were further validated using measured soil organic carbon storage levels from our field survey.

For belowground biomass of grassland on the Tibetan Plateau, we used models which relate measured belowground biomass in the field (BGB) to the remote sensing vegetation index. Considering the significant differences of photosynthetic carbon allocation patterns between above and belowground organs among different types of grassland, we established separate models for meadow, steppe, and desert steppe.

where BGB_am is belowground biomass in alpine meadow; BGB_as is belowground biomass in alpine steppe; BGB_ad is belowground biomass in desert steppe.

For species richness of grassland on the Tibetan Plateau, we also adopted a model from our previous study for the estimation of species richness across the grassland (Li et al., 2020). The model was developed to include 20 environmental variables relevant to species settlement, development, and survival in order to predict the spatial pattern of plant species richness. This model was validated by our measured species richness data from the field, and it is the first attempt to map the species richness patterns across the entire grassland of the Tibetan Plateau (Please see more detailed information in Li et al. (2020) ).

The numbers of livestock in the northern Tibetan Autonomous Region were acquired from the yearly “Statistical yearbook”. Sheep and large livestock were both included in the original livestock data. Livestock were standardized into sheep units using the criteria that one sheep is equal to one sheep unit and one large livestock is equal to four sheep units (Fan et al., 2010). Grazing intensity (GI) was calculated as the number of livestock in a certain county divided by productivity in August in fenced grassland per square meter times grassland area within the county.

Daily meteorological records were obtained from China Meteorological Data Service Center (CMDC, http://data. cma.cn) for the 200 meteorological observation stations that are located within and around the Qinghai-Tibetan Plateau. Chen et al. (2014) generated monthly raster surfaces of temperature and precipitation with spatial resolutions of 1 km × 1 km by using ANUSPLIN 4.3 (Hutchinson and Xu, 2004), and found that the interpolated climatic values were very consistent with field observations. We integrated the monthly raster surfaces of temperature and precipitation which were interpolated by Chen et al. (2014) into growing season temperature (GST) and precipitation (GSP), respectively. Finally, we extracted annual GSP and GST values for each site according to the geographic coordinates in ArcGIS 10.2 for further analyses.

Data for each of the four functions were standardized to a mean of zero and standard deviation of one to allow for the easier comparison of effects. An averaging approach was used to calculate the multifunctionality of grassland by summing up the standardized values of each of the individual functions, then dividing by the number of functions, 4. In this calculation, we weighted the four individual functions equally. Density distributions of multifunctionality were analyzed individually in the alpine meadow, alpine steppe, temperate meadow, temperate steppe, and temperate desert steppe. A fragmentation intensity index was quantified using the relative cover of bare soil, calculated as the difference between 1 and the relative cover of vegetation patches within the plots.

Multifunctionality data were extracted and calculated according to the county boundaries in the northern Tibetan autonomous region, which includes alpine meadow, alpine steppe, temperate meadow, temperate steppe and temperate desert steppe. Regression analyses were performed to test the relationships of multifunctionality with climatic factors and grazing intensity across the grassland types. We also used bivariate regression analysis to test for relationships of multifunctionality with climatic factors. After confirming that precipitation is the primary factor which determines multifunctionality across the grassland, we further assessed the relationships between residuals of multifunctionality and grazing intensity in each type of grassland to reveal the net effect of grazing intensity on multifunctionality. We excluded temperate meadow from the bivariate analysis due to its small area in the transect. We employed piecewise linear regression with the segmented package (Muggeo, 2003) in R 3.5.2 (http://www.r-project.org/) to identify the threshold values of multifunctionality and individual functions along the fragmentation intensity (bare ground cover). Maps were generated using ArcGIS 10.2 (Environmental Systems Research Institute, ESRI).

Multifunctionality showed a decreasing trend from the alpine meadow concentrated in the eastern part to the desert steppe in the western part of the Tibetan Plateau (Fig. 2a). Among the different types of grassland, alpine meadow and temperate steppe showed high similarity in their patterns of multifunctionality, while temperate meadow and temperate desert steppe presented high similarity in their patterns of multifunctionality (Fig. 2b). Specifically, the highest value of multifunctionality in the alpine meadow was 0.91, and the mean value was 0.42. The highest and mean values of multifunctionality in the temperate steppe were 0.76 and 0.30, respectively. The highest value of multifunctionality in the alpine steppe was 0.53, and the mean value was 0.16. Most of the values of multifunctionality were concentrated around 0.13 in temperate meadow and temperate desert steppe, and the mean values of multifunctionality were 0.13 and 0.10 in these two types of grassland, respectively (Fig. 2). The lowest values of multifunctionality of around 0.05 were comparable in the different types of the grassland.

Both climate and livestock grazing significantly influenced the multifunctionality of grassland across the Tibetan Plateau. Annual precipitation showed a positive effect on multifunctionality. However, the effect of annual mean temperature was not significant, and therefore it was not included in the following equation. Grazing intensity showed a negative effect on multifunctionality. The effect of the interaction between annual precipitation and grazing intensity on multifunctionality was also significant.

A bivariate regression analysis showed that the relationship between multifunctionality and annual precipitation presented similar trends in alpine meadow, alpine steppe, temperate steppe and temperate desert steppe (Fig. 3). Meanwhile, the variation of multifunctionality residuals with grazing intensity showed shifts from positive to negative trends in each grassland type, suggesting that high grazing intensity has a negative effect on grassland multifunctionality.

Both above and belowground biomass showed similar decreasing trends with an increase of bare ground cover (Fig. 4). The difference between them was that aboveground biomass showed a smooth decrease with the decrease of vegetation cover, while belowground biomass showed a transition from smooth decrease to rapid decline as the bare ground cover exceeded 75% (Fig. 4). Species richness and SOC decreased slowly with increasing bare ground cover, but they shifted to quick declines after the bare ground cover surpassed 75% (Fig. 4). As a consequence, multifunctionality decreased gradually with an increase of bare ground cover, then switched to a rapid collapse after the bare ground cover amounted to 75% or more (Fig. 4). The threshold value of multifunctionality was 0.233.

Our research results showed a clear spatial pattern of multifunctionality in the grassland on the Tibetan Plateau by integrating the four important individual functions. Generally, multifunctionality demonstrated the highest value in alpine meadow. Multifunctionality values corresponding to the four individual functions that passed through the 3/4 quantile all occurred in the alpine meadow. In alpine steppe, only a few values of the above and belowground biomass passed through the 3/4 quantile. None of the individual functions passed through the 3/4 quantile in the other types of grassland. The lower values of multifunctionality were comparable in the different types of the grassland. Annual precipitation explained a large amount of the variation of multifunctionality across four primary types of grassland, which showed significantly positive effects on multifunctionality. Grazing intensity further explained the variation of multifunctionality, which showed a significantly negative effect on multifunctionality across the four types of grassland. The consistently rapid declines of belowground biomass, SOC, and species richness resulted in the collapse of multifunctionality as bare ground cover exceeded 75%, which is corresponded to a value of 0.233 for the multifunctionality.

The spatial pattern of grassland multifunctionality is consistent with the gradient of annual precipitation from the eastern to western parts of the Tibetan Plateau. Our specific analysis demonstrated that annual precipitation explained a large portion of the variation of multifunctionality using the extracted multifunctionality values across alpine meadow, alpine steppe, temperate steppe, and temperate desert steppe along the transect. Moreover, the bivariate regression analy-sis showed that the slopes of the relationships between multifunctionality and annual precipitation are steeper in alpine steppe and temperate steppe than in alpine meadow, which indicates that annual precipitation played a more important role in determining the multifunctionality in steppe. The alpine grassland on the northern Tibetan Plateau is mainly distributed in arid and semi-arid areas, where precipitation is the key factor limiting plant growth and development (Wu et al., 2014). Previous studies showed that individual functions, such as above and belowground biomass, species richness and soil organic carbon storage, were closely associated with annual precipitation across the northern Tibetan Plateau (Wu et al., 2014). Grassland is the most widely distributed vegetation type on the Tibetan Plateau, occupying about 75% of the total area (Liu et al. 2018). Alpine meadow, temperate meadow, alpine steppe, temperate steppe, and temperate desert steppe have devel-oped due to differences in temperature and precipitation. Although temperature is also an important factor controlling grassland multifunctionality (Wu et al., 2014), our results did not differentiate the significant effect of temperature. One potential reason is that the transect used in this study is characterized by a typical water gradient, so the strong effect of precipitation would be expected to override the effect of temperature. In addition, the interaction between temperature and precipitation makes it difficult to differentiate the effects of each of them from the other (Wu et al., 2014).

The results of our regression analysis demonstrated that grazing intensity negatively contributed to the variation of multifunctionality across alpine meadow, alpine steppe, and desert steppe in the northern Tibet transect. The result suggests that grassland located in regions with low annual precipitation is more sensitive to grazing intensity. Previous studies have shown that the alpine meadow dominated by Kobresia spp. is an equilibrium grazing system, with a coefficient of variance in the interannual precipitation variability of well below 30% (Ellis and Swift, 1988; Wehrden et al., 2012). In contrast to more variable semi-arid and arid non-equilibrium systems, equilibrium systems may degrade if livestock numbers increase to the point that the carrying capacity is exceeded (Liu et al., 2018b). Estimates of grazing-induced degradation vary, and the most frequently quoted value for the Tibetan highlands is that 30% of the grasslands are degraded (Harris, 2010; Wang and Wesche, 2016; Liu et al., 2018a). In the Ruoergai Plateau (Zoige Plateau), the ecosystem services value expressed as a multiple of the gross domestic product (GDP) has decreased by about 84% between 1990 and 2005 (Li et al., 2010). The intermediate disturbance hypothesis (IDH) posits that species richness (and diversity) at the local scale peak either at intermediate values of disturbance frequency or at intermediate spatial extents of disturbance (Grime, 1973). Moreover, studies have also found positive effects of intermediate disturbance on plant growth and productivity. In our study, precipitation explains the large variations of the multifunctionality, and the residuals showed a shift from positive to negative with an increase of grazing intensity. These results indicate a positive effect on multifunctionality at low grazing intensity, and a negative effect at high grazing intensity in all four types of grassland. It should be noted that temperate desert steppe possesses low multifunctionality, and even the low grazing intensity could cause a negative effect on its multifunctionality. Stakeholders should pay more attention to minimizing grazing in the temperate desert steppe to avoid catastrophic transition to desert.

We found consistent rapid decreases of species richness, SOC, and belowground biomass after vegetation cover declined to around 75%, which signify an impending collapse of the multifunctionality. That is to say, the alpine grassland would lose the multiple functions as the multifunctionality decreases and approaches a level around 0.233. A previous study has shown that degradation since the year 2000 is proceeding only in the less productive Kobresia pastures in the western part of the highlands, where it is largely driven by slight decreases in precipitation in combination with rising temperatures (Lehnert et al., 2015). In contradiction to the widely assumed high importance of human influence on the degradation process, stocking numbers were not strongly correlated with larger-scale plant cover changes (Lehnert et al., 2015). Thus, pasture degradation is a consequence driven by multiple factors, such as regional grazing, climate, soils and rangeland management (Wang and Wesche, 2016).

Grasslands on the Tibetan Plateau are fundamental for household welfare, species conservation, and environmental protection. Assessing ecosystem functions and the services they provide are urgently needed from multiple dimensions. Our results showed that high multifunctionality occurs in the alpine meadow, suggesting high carrying capacity. However, the lower multifunctionality values obtained are comparable among the major types of grassland. The significant negative effect of grazing intensity on the multifunctionality contributed to the large decrease of the multifunctionality within each type of grassland. We also found that the collapse of the multifunctionality occurred as vegetation cover became less than 25%, corresponding to the rapid decline of individual functions. Our results are expected to provide evidence and direction for sustainable development of alpine grassland and restoration management.