key: cord-0800125-4xwxuebk authors: Hua, Ting; Zhao, Wenwu; Cherubini, Francesco; Hu, Xiangping; Pereira, Paulo title: Sensitivity and future exposure of ecosystem services to climate change on the Tibetan Plateau of China date: 2021-08-24 journal: Landsc Ecol DOI: 10.1007/s10980-021-01320-9 sha: 1bcefec8309556aa3403a9034da075ab664680dd doc_id: 800125 cord_uid: 4xwxuebk CONTEXT: Climate change has imposed tremendous impacts on ecosystem services. Recent attempts to quantify such impacts mainly focused on a basin or larger scale, or used limited time periods that largely ignore observations of long-term trends at a fine resolution, thereby affecting the recognition of climate change’s effect on ecosystem services. OBJECTIVES: This study conducts a detailed and spatially explicit recognition of climate change’s effect on ecosystem services and provides an intuitive map for decision-making and climate change adaptation planning. METHODS: We used long-term time series of ecosystem service assessments and various future climate scenarios to quantify the sensitivity and future exposure of ecosystem services to climate change on the Tibetan Plateau. RESULTS: Carbon sequestration (CS) and habitat quality experience significant growth, while water retention did not show any trend. Sensitivity patterns of these ecosystem services vary largely. For CS, more than half of the pixels showed a positive sensitivity to climate change, even though the degree of sensitivity is not high. There is substantial spatial heterogeneity in the exposure of ecosystem services to future climate changes, and high levels of future climate change increase the intensity of exposure. CONCLUSIONS: This study illustrates the complex spatial association between ecosystem services and climatic drivers, and these findings can help optimize local response strategies in the context of global warming. For example, the existing protected areas have notable conservation gaps for disturbance of future climate change on ecosystem services, especially in the southeastern part of the study area. SUPPLEMENTARY INFORMATION: The online version contains supplementary material available at 10.1007/s10980-021-01320-9. decision-making and climate change adaptation planning. Methods We used long-term time series of ecosystem service assessments and various future climate scenarios to quantify the sensitivity and future exposure of ecosystem services to climate change on the Tibetan Plateau. Results Carbon sequestration (CS) and habitat quality experience significant growth, while water retention did not show any trend. Sensitivity patterns of these ecosystem services vary largely. For CS, more than half of the pixels showed a positive sensitivity to climate change, even though the degree of sensitivity is not high. There is substantial spatial heterogeneity in the exposure of ecosystem services to future climate changes, and high levels of future climate change increase the intensity of exposure. Conclusions This study illustrates the complex spatial association between ecosystem services and climatic drivers, and these findings can help optimize local response strategies in the context of global warming. For example, the existing protected areas have notable conservation gaps for disturbance of future climate change on ecosystem services, especially in the southeastern part of the study area. Keywords Ecosystem services Á Climate change Á Tibetan Plateau Á Sensitivity Á Exposure Ecosystem services are defined as benefits that humans obtain from the ecosystems (Costanza et al. 1997) and are key to achieve the Sustainable Development Goals (McElwee et al. 2020; Yang et al. 2020) . Climate change affects ecosystems' structure, composition, function, and the supply of ecosystem services (Holmgren and Hirota 2013; Grimm et al. 2013) . Several works explored climate change's influence on ecosystem services and put forward adaptive management strategies (Tanner-McAllister et al. 2017; Langhans et al. 2019 ). There are numerous pieces of evidence that climate change induces great disturbances in regulating and provisioning services (IPCC et al. 2019) . The global mean temperature is expected to rise by 1-5°C by the end of the twenty-first century, depending on future emission trajectories (IPCC 2018; Hoegh-Guldberg et al. 2019 ) and the warming trend are generally expected to intensify the negative impacts on ecosystem services, especially in vulnerable areas (Westerling et al. 2011; Seidl et al. 2014) . Therefore, it is vital to increase our understanding of the spatial association between ecosystem services and climate change. This is key to identify ecosystem services' exposure to future climate changes and to prevent the decline of ecosystem services under a changing climate. Previous works highlighted that the effects of climate change on ecosystem services are servicedependent and spatially explicit. Climate change impacts can have diverse effects in different regions due to different ecological, hydrological, and biophysical processes (Qiu et al. 2019) . For instance, warming may benefit vegetation distribution and productivity by boosting photosynthesis and plant growth rates (Sullivan et al. 2008 ). However, this increases evaporation and soil organic matter decline (Schroter et al. 2005) . Therefore, the increasing temperature may positively affect forest productivity and bioenergy production; nevertheless, it may have adverse impacts on water supply and soil fertility (Schroter et al. 2005; Weiskopf et al. 2020) . Several works specifically focused on the effects of climate change on ecosystem services. However, most of the works focused on the impacts of climate change on ecosystem services, just for a few snapshots (Qiu et al. 2020) . Considering that climate fluctuations take a long period to be stabilized, the study with short-term periods may limit their value. Besides, most of the studies had a coarse scale (e.g., basin or larger), and only a few had a finer resolution (e.g., pixel scale) (Hao et al. 2017a, b; Hou et al. 2020; Shen et al. 2020) . Coarse-scale studies overlook important information that may be essential for management. Thus, it is necessary to conduct detailed and high-resolution studies, essential to disentangle service-dependent and identify climate change impacts (Cord et al. 2017; Mao et al. 2019; Xu et al. 2020) . Increasing the understanding of ecosystem services' relationships contributes to achieving better ecosystem management that increases synergies and mitigates the trade-offs (Qiu et al. 2018; Zheng et al. 2019) . Understanding the trade-offs between ecosystem climate drivers is crucial for effective management as well. Multiple driving forces often inhibit each other's impact on ecosystem services. At the same time, ecosystem services' drivers may influence a single ecosystem service, but have limited effects on other services, or dramatically change multiple services (Bennett et al. 2009; Nelson et al. 2013) . Therefore, identifying the trade-offs and synergies among the ecosystem services' climatic drivers is essential. The Tibetan Plateau is an area that supplies a wide range of important ecosystem services (Sun et al. 2012; Liu 2020) , and it is one of the most sensitive areas to global climate change (Song et al. 2014; Kuang and Jiao 2016) . Long-term overgrazing, combined with very harsh climate conditions, increased the vulnerability of ecosystem services supply to climate change (Newbold et al. 2015; Lehnert et al. 2016) . The future value of ecosystem services is highly uncertain under different climate scenarios (Kubiszewski et al. 2017 ) and depends on policy choices. In this context, it is of utmost importance to assess the impacts of future climate change on ecosystem services and develop ecological-oriented regional management policies. This work aims to analyze the sensitivity of ecosystem services to historical climate change and infer future exposure of ecosystem services to climate change on the Tibetan Plateau. In particular, the objectives of the research are: (1) to measure the spatial patterns and trends of primary ecosystem services [i.e., carbon sequestration (CS), water retention (WR), and habitat quality (HQ)]; (2) to quantify the sensitivity of these ecosystem services to climate change (precipitation and temperature) at pixel scale; (3) to identify the spatial association between the ecosystem services and their climatic drivers, and (4) to investigate the exposure of ecosystem services under multiple future climate scenarios. Tibetan Plateau is located in western China, and covers Xizang, Qinghai, and part of the Gansu, Xinjiang, Yunnan, and Sichuan provinces of China (Fig. 1) . The mean elevation of the Tibetan Plateau is more than 4000 m above sea level, and it is known as ''the Roof of the World''. It is also the largest ice reservoir outside the polar regions and a source of several major Asian rivers such as the Yellow River, the Yangtze River, and the Mekong River. These rivers supply water for more than a billion people downstream (Immerzeel et al. 2020) . Given that it is known as ''the Third Pole of the Earth'' (Qiu 2008) and ''the Water Tower of Asia'' (Barnett et al. 2005; Immerzeel et al. 2010) , it acts as an important ecological security barrier including climate regulation and biodiversity conservation. Tibetan Plateau is characterized by low temperature, intense radiation, and low precipitation (see Fig. S1 for the spatial distribution of temperature and precipitation). From southeast to northwest, the vegetation types are broadleaf forest, needleleaf forest, shrub, grassland, and desert, respectively (Fig. S2 ). These ecosystems are fragile and sensitive to climate change ). We focused on three ecosystem services (i.e., CS, WR, and HQ) given the critical role of the Tibetan Plateau in climate regulation, water supply, and biodiversity conservation (Barnett et al. 2005; Immerzeel et al. 2010 Immerzeel et al. , 2020 . MODIS-based satellite and meteorological datasets were used to evaluate these ecosystem services (Table 1) . To investigate the exposure of ecosystem services to climate change, both historical baseline and future climate datasets from the World-Clim database are used. Baseline climate data were acquired from WorldClim version 2.1 (Fick et al. 2017) . Four climate scenarios considering representative concentration pathways (RCPs) and shared socioeconomic pathways (SSPs) were used for possible average future climate variables (annual mean temperature and annual precipitation) generated from 8 CMIP6 global climate models (BCC-CSM2-MR, CNRM-CM6-1, CNRM-ESM2-1, CanESM5, IPSL-CM6A-LR, MIROC-ES2L, MIROC6, and MRI-ESM2-0): SSP12.6, SSP2-4.5, SSP3-7.0, and SSP5-8.5, listed in the order of magnitude of temperature changes. Future climatic conditions are Fig. 1 Location of the Tibetan Plateau considered for a period of 20 years, between 2040 and 2060. For example, these future scenarios (2040-2060) resulted in an average temperature increase in the Tibetan Plateau of 2.5°C in SSP1-2.6 and 3.6°C in SSP5-8.5 (Fig. S3 ). The indicators and data used are described in Table 1 . The resolution of all raster data was resampled at 1 km 9 1 km. This analysis is centered on the ''sensitivity'' and ''exposure'' of ecosystem services to climate change. These concepts reflect the response of a system to stress (Turner et al. 2003) . ''Sensitivity'' is defined as the degree to which a system is modified or affected by perturbations, and ''exposure'' is the nature and degrees to which a system experiences environmental or social stress (Adger 2006) . According to the characteristics of the research object, the meaning of ''sensitivity'' and ''exposure'' has been also constantly enriched (Adger 2006; Smit et al. 2006; Smith et al. 2014; Seddon et al. 2016; Li et al. 2018a ). Adger (2006) thought vulnerability is a function of climate variation to which a system is exposed, its sensitivity, and its adaptive capacity. Smit et al. (2006) highlighted the importance of quantifying sensitivity and exposure in risk management and adaptation to climate change. Remote sensing technology has being increasingly applied and it uses vegetation indices (e.g., NDVI) as a proxy to quantify the ecosystem (Smith et al. 2014; Seddon et al. 2016; Li et al. 2018a) . In our analysis, we directly used the ecosystem services instead of vegetation indices, and evaluated the effect of climate change on the services. Therefore, following Li et al. (2018a) , sensitivity is thus understood as the degree of disturbance of climate change on ecosystem services or the ecosystem services' resistance to climate change. Exposure is to be interpreted as the gains or losses of ecosystem services driven by future climate. CS is an important regulation service that controls the atmospheric carbon concentration (Meersmans et al. 2016; Bai et al. 2018) . We evaluated the amount of CS based on net primary productivity (NPP). The basic assumption is that vegetation fixes 1.63 units of carbon when accumulating per unit NPP according to the photosynthesis equation (Li and Zhou 2016) . The formula is as follows: NPP data (units: kgÁC/m 2 ) is MOD17A3 from numerical terradynamic simulation group (Zhao et al. 2010) . WR refers to the water retained by ecosystems. We estimated WR based on the following formula (2, 3). WR is the amount of water retention services, P is precipitation (units: mm), R is storm runoff (units: mm), ET is evapotranspiration (mm), and a is the runoff coefficient based on Ouyang et al. (2016) , which was estimated from more than 300 publications. Precipitation data was from the climate hazards group infrared precipitation with station data (CHIRPS) dataset (Funk et al. 2015) , and evapotranspiration data is from MODIS global evapotranspiration project (MOD16) (Mu et al. 2007 ). HQ refers to the ability of an ecosystem to provide a suitable environment for individual and population persistence (Hall et al. 1997) . It is considered as a robust proxy to indicate the effectiveness and efficiency of the biodiversity conservation status (Maes et al. 2012; Moreira et al. 2018 ). Habitat quality module in the integrated valuation of ecosystem services and trade-offs (InVEST) model version 3.8.4 was applied to evaluate HQ by combining land cover data and multiple threat factors (Sharp et al. 2015) . The model's input requirements included land use and land cover data, threat factors, habitat types, and their sensitivities to each treat factor (Table S1) . We set urban, cropland, and crop/natural vegetation mosaics as threat factors and the decay laws of these threat factors are linear. Relevant parameters referred to the Tibetan Plateau research of Hou et al. (2020) and were slightly adjusted according to the classification system of MCD12Q1 (see Table S1 for details). For instance, according to the classification system of MCD12Q1, there are some mixed land cover types, such as crop/natural vegetation mosaics and mixed forests. Therefore, we set habitat suitability index and other parameters following the composition of the mixed land cover types and parameters of individual land types. We judged the robustness of the evaluation by comparing all the ecosystem services with other studies using Spearman correlation coefficient for all ecosystem services. Significant differences were considered at a p \ 0.05. HQ model was validated using above-ground biomass carbon and tree density coverage as used in previous works (Brandt et al. 2018; Gomes et al. 2021 ). A moving window was used to detect the sensitivity of three types of ecosystem services to historical climate change (temperature and precipitation) following Abel et al. (2019 Abel et al. ( , 2021 . A multi-year assessment of ecosystem services, temperature, and precipitation data was used (Fig. 2a) . A linear regression slope between ecosystem services and temperature or precipitation was calculated for all pixels within the moving window and considered time period (Fig. 2b , c). The regression slope was assigned to the centre pixel of the moving window as the sensitivity of ecosystem services to climate change (Fig. 2d ). The regression is not fitted, and the regression slope was not assigned to the centre pixel when more than twothirds of the pixels have no data. As the window moved, we can identify the pixel scale sensitivity of the whole study area (Fig. 2e ). The sensitivity of ecosystem services to climate change in the centre pixel is high when the regression slope's absolute value is high. A slope greater than 0 indicates positive sensitivity, while a slope lower than 0 shows a negative sensitivity. Positive sensitivity indicates the change direction of ecosystem services is consistent with climate factors and vice versa. The selection of the moving window instead of single-pixel ensured a high probability of statistically significant regression results by including a larger number of pixels in the analysis and considering the neighbour environmental effect (Abel et al. 2019) . Furthermore, according to Abel et al. (2019) , when a step of time of a temporal-spatial moving window is 14 or 15 years, the optimal size of the moving window is 7 9 7 according to the significance of the regression slope and the quality of the regression fit. Additionally, we conducted a sensitivity test on the moving window size to determine a suitable size (Fig. S4 ). For example, we took the relationship between CS and precipitation and tested different sizes from 3 9 3 to 19 9 19. Subsequently, we generated nearly 1500 random points to extract the p-value of regression analysis, and the proportion of points with a p-value less than 0.05 was counted. When the size is 7, there is a high proportion of statistical significance. Furthermore, a size of 7 can capture better-detailed features compared to the larger size. Therefore, the size of the moving window was set to 7. The annual precipitation data comes from CHIRPS dataset. Temperature data is from the day temperature and night temperature of the MOD11A2 and MYD11A2 products (see Table 1 for details). We used the average method to get the annual temperature product based on the above datasets. Further, we assumed that the effects of climate change on ecosystem services are linear, even though the relationship between environmental variables and ecosystem services are usually non-linear like logarithmic or exponential (Koch et al. 2009; Schlenker et al. 2009; Sutherland et al. 2016; Peng et al. 2017; Qiu et al. 2019 ) in a larger observation interval. This is based on the consideration that our analysis was conducted to explore the impact of climate change on ecosystem services for single or neighbour pixels, and the degree of climate change of a single pixel is not so large. Therefore, we can approximate the non-linear relationship as a linear relationship. Spatial association of ecosystem services' sensitivity to climate change We explored: (1) whether the change direction of the ecosystem service (ES1) is consistent with the change direction of precipitation and temperature (Fig. 3a) , and (2) whether the change direction of two ecosystem services (ES1 and ES2) is consistent with the change direction of precipitation/temperature (Fig. 3b) . This process provides a spatially explicit reference for the relationship between the effects of two climatic Fig. 2 Flow chart of the calculation of the sensitivity of ecosystem services to climate change drivers and the relationship between the effects of a single climatic driver on two ecosystem services. Specifically, each pixel in the study area was divided into four states (positive-negative, positive-positive, negative-positive, negative-negative) according to the direction of sensitivity to precipitation and temperature (positive or negative). Positive-negative means that the pixel's sensitivity to precipitation is positive and that to temperature is negative (blue colour in Fig. 3a) . Also, we calculated the proportion of the four states (e.g., perc1) in the four states. Spearman correlation coefficient (e.g., r1) was applied to assess the relationship between the sensitivity to precipitation and/or to temperature, which indicates the similarity of the effect intensity of two climatic factors. ''r1-4'' refers to the Spearman correlation coefficient of the sensitivity of ecosystem services to temperature and/or precipitation. Notice that in Fig. 3a , ''r1'' is the Spearman correlation coefficient between the sensitivity of the pixels with positive sensitivity to temperature and negative sensitivity to precipitation. In Fig. 3b , ''r1'' is the Spearman correlation coefficient between the sensitivity of the pixels with positive sensitivity of ES1 to precipitation (temperature) and negative sensitivity of ES2 to precipitation (temperature). For example, in Fig. 3a , if the coefficients of ''r2'' and ''r3'' are closer to 1 (-1), and the coefficients of ''r1'' and ''r4'' are closer to -1 (1), high sensitivity to temperature is always accompanied by high (low) sensitivity to precipitation. Similarly, we follow the instruction in Fig. 3b to detect whether the two ecosystem services (ES1 and ES2) are consistent with the trend of precipitation/temperature. Significant correlations were considered at a p \ 0.05. The calculation of exposure only considered extrinsic climate conditions. In combination with the sensitivity of ecosystem services to climate change, future climate velocity is applied to estimate the exposure. The future climate change velocity [V, see formula (4)] was obtained by dividing the average long-term temporal climate gradient by the average spatial Flow chart of detection of the association of sensitivity of ecosystem services to climate change. a Different colours represent different states, describing the direction of sensitivity of single ecosystem service to precipitation and temperature (positive or negative). For example, blue colour (positivenegative) represents the pixel's sensitivity to precipitation is positive and that to temperature is negative. b Illustrates the direction of sensitivity of two types of ecosystem service (ES1 and ES2) to temperature/precipitation. b The blue colour represents the sensitivity of ES1 to temperature (precipitation) is positive, and the sensitivity of ES2 to temperature (precipitation) is negative in this pixel. ''perc1-4'' and the length of the bar means the proportion (%) of the pixel area in the four states, respectively. ''r1-4'' means the Spearman correlation coefficient between the sensitivity of ecosystem services to temperature and/or precipitation. See method for more explanation climate gradient (Loarie et al. 2009 ). The spatial climate gradient represented the temperature or precipitation difference of adjacent grid cells, measured in°C km -1 and mm km -1 . The calculation of climate change velocity was based on climate data between the contemporary and future (2041-2060) periods following the research of Li et al. (2018a) . Four future pathway scenarios (SSP1-2.6, SSP2-4.5, SSP3-7.0, and SSP5-8.5) were selected from eight global climate models. We used the average output of eight global climate models to characterize the future climate to reduce a single model's uncertainty, according to Li et al. (2018a) . Additionally, we conducted a standard error estimate to map the future velocity uncertainty of climate change from 8 global climate models for each pixel and each future scenario. Historical baseline climate data from WorldClim database was used to measure the spatial gradient with 3 9 3 moving windows following the average maximum technique (Loarie et al. 2009; Burrows et al. 2011 ). The specific procedure was the following: for each moving window such as formula (5), spatial gradients of north-south direction were calculated as the difference for all northern and southern pairs (a-d, b-e, c-f, d-g, e-h, and f-i) within a moving window divided by the distance between them (111.325 km per degree). Same as north-south direction, spatial gradients of west-east direction were calculated as the difference for all eastern and western pairs (a-b, b-c, d-e, e-f, g-h, h-i) within a moving window divided by the distance between them [Distance WÀE , see formula (6)]. where lat is the latitude of the pixel to be processed. Average spatial gradients of north-south and westeast direction were assigned to the center pixel e. Then, overall spatial gradients were calculated as the vector sum of spatial gradients of these two types of directions as the formula (7). where S, S NS , and S EW is the overall spatial gradient, north-south spatial gradient, and east-west spatial gradient, respectively. And h is 0°as North and 180°as South. After the above calculations, we can get the temperature and precipitation change velocity separately. To evaluate the sensitivity differences observed in the different regions, the sensitivity result (see Methods for details) was introduced to evaluate the exposure. Ultimately, the formula of exposure of ecosystem services to future climate change is the following: where E is the exposure of ecosystem services to future climate change; Sensitivity is the sensitivity of ecosystem services to climate change (temperature or precipitation); V is the velocity of climate change (temperature or precipitation) for four future scenarios (see formula (4)). The temporal gradient was calculated based on the difference between baseline and future climate data in each pixel for the four climate scenarios. The formula of exposure and enhanced vegetation index (EVI), according to Li et al. (2018a) , was adjusted to ecosystem services in our framework. A zonal average of negative and positive exposure per 50 m altitude was applied to investigate the altitude discrepancies in the ecosystem services exposure to climate change. Also, we further explored the exposure intensity of ecosystem services to climate change along the altitude gradient under four climate scenarios. The low-emission pathway (SSP1-2.6) was regarded as the baseline. Exposure intensity is the ratio of exposure of SSP2-4.5, SSP3-7.0, SSP5-8.5 to the baseline (SSP1-2.6). Protected areas data was introduced to identify the gap between the current distribution of protected areas and adapt future exposure. According to the sign (negative/positive) of the exposure value to climate change, all pixels were reclassified into eight classes: ''positive-negative'', ''positive-positive'', ''negativenegative'', ''negative-positive'', ''positive-0'', ''negative-0'', ''0-positive'' and ''0-negative''. Among them, ''positive-negative'' means that the exposure to precipitation is positive and that to temperature is negative, and ''0'' means that exposure to precipitation/temperature is 0. ''Negative-negative'', ''negative-0'' and ''0-negative'' mean that future climate change may lead to losses of ecosystem services. Subsequently, we compared the gap between the distribution of these three states and protected areas based on spatial overlay analysis. In our framework, ''positive-negative'' and ''negative-positive'' were not considered, since we cannot yet determine the final direction of influence of future climate change. All the spatial and statistical analyses were conducted based on ArcGIS 10.5 and Matlab R2014a. The mean CS between 2000 and2015 was 303.992 t/km 2 (Fig. 4a) . The linear regression trend exhibited a significant increase with an average rate of 2.654 t/km 2 year (p \ 0.01) (Fig. S5) . High CS services occurred in the southeast and eastern regions, including Sannan, Nyingchi, Liangshan, Diqing, Aba, and Gannan. CS was close to zero in most areas, mainly located in the central and western regions of the Tibetan Plateau. The results suggested that the area where CS increased was larger than the area where it decreased (55.24% Vs. 9.91%). This contrast was small in the areas with statistically significant trends (13.60% Vs. 0.26%, p \ 0.05, Fig. S6 ). The annual average amount of WR was approximately 115.423 mm/km 2 from 2001 to 2014 (Fig. 4c) . It was higher in the southeast, including Lhoka and Nyingchi. The annual amount of water retained had a non-significant decreasing trend of -0.158 mm/km 2 (p = 0.876). Approximately 27.5% of the Tibetan Plateau exhibited an increasing trend of annual trend (3.82% with a significant increase, p \ 0.05). This is especially evident in the areas between 30 and 35°N, whereas decreasing WR was concentrated in the southeast region. This is consistent with the highvalue area of WR service that is sporadically seen in Haibei and Hainan prefecture. HQ changed slightly between 2001 and 2015 (Fig. 4e) . The annual average HQ was 0.3049/km 2year, while the annual variation trend was approximately 5.942e-05/year (p \ 0.01). The high values were mainly distributed in the southeast part of the studied area (forest area) or scattered (near the lakes). Overall, HQ is very stable, and only 3.50% (2.02%) of the Tibetan Plateau exhibited a significant increased (decreased) trend. The cross-validation combined with Spearman correlation coefficient showed acceptable robustness for CS and HQ evaluation (Fig. S7) . And the regression analysis also showed moderate but acceptable robustness of the HQ model: R 2 = 0.40 for above-ground biomass carbon and R 2 = 0.30 for tree density (Fig. S8) . See Figs. S7 and S8 for more detailed validation. The sensitivity varies notably geographically for all ecosystem services, whereas the sensitivity of ecosystem services to the climate in most areas is not high (Fig. 5) . As for CS, the average sensitivity to precipitation is minimal (-0.0015 t/mmÁkm 2 ). 56.94% of the pixels had a positive sensitivity and about half a negative sensitivity (26.72%, Fig. 5a ). High-negative sensitivity to precipitation regions was located in Qamdo and Diqing, whereas regions with highpositive sensitivity to precipitation are around Qinghai Lake and Sannan. As for the relationship between CS and temperature, 47.49 and 36.21% of study regions exhibit a positive and negative sensitivity, respectively. The negative sensitivity areas are mainly distributed in regions with high topographic relief, like the Hengduan Mountains. There is an apparent positive sensitivity between WR and precipitation with 0.614 mm/(mmÁkm 2 ) average sensitivity (Fig. 5c) . 41.56% of the pixels showed a positive sensitivity, and only 0.43% of the pixels showed a negative sensitivity. Compared to the sensitivity of WR to precipitation, the sensitivity to temperature was generally negative (average: -1.172 mm/°CÁkm 2 ) (Fig. 5d) . The high negative sensitivity was concentrated in the South of the study region like Nyingchi. There was a relatively low sensitivity between climate change and HQ (Fig. 5e, f) . The average sensitivity of HQ to precipitation and temperature was -4.85e-05/mm km 2 and -1.172°C km 2 , respectively. The regions sensitive to climate change were mainly concentrated in the southern and southeastern parts of the Tibetan Plateau. Spatial relation of ecosystem services' sensitivity to climate change The first row of Fig. 6 reflects whether the three ecosystem services' sensitivity to precipitation and temperature changes is consistent. In CS, the change of direction in temperature, precipitation, and ecosystem services is consistent in 27.78% of the area (Fig. 6a) . In 57.13% (27.78 ? 29.35%) of the pixels, we observed a CS change of direction consistent with precipitation. In 45.88% (29.35 ? 18.10%) of the area, the sensitivity of CS to temperature has an opposite trend to that identified in precipitation. There was a significantly high negative correlation when sensitivity to precipitation was opposite to sensitivity to temperature (-0.52 and -0.54, p \ 0.01, yellow and blue in Fig. 6a ). This shows that high positive (negative) sensitivity to precipitation is accompanied by a high negative (positive) sensitivity to temperature. In this area, temperature and precipitation have antagonistic effects. Among them, the precipitation effect on ecosystem services is larger than the opposite effect of temperature in 29.35% of the pixels. Also, the CS change of direction matched the observed precipitation. However, in 18.10% of the area, precipitation cannot offset the effects of temperature, and the direction of CS changes is consistent with the temperature. Concerning WR, over 40% of the area showed that the change direction was consistent with the precipitation observed (Fig. 6b) . Among them, precipitation and temperature were consistent with the direction of WR in 21.54% of the regions. In 18.86% of the pixels, the temperature direction was not consistent with WR and precipitation observed. If the area with zero sensitivity is not considered, in near 100% of the regions, the change direction of WR is the same direction as precipitation. This shows that the change of direction in WR matched with the precipitation in most cases. However, for HQ, the four states' proportion was relatively close (5-7%, Fig. 6c ). This showed that the change of direction in HQ was not correlated with climate change. We also explored whether the direction of sensitivity of two different ecosystem services to a single climatic driver (precipitation or temperature) is the same (the second and third rows of Fig. 6 ). In 29.06% of the regions, the change of direction in precipitation, CS and WR were the same, and the correlation was 0.43 (p \ 0.01, Fig. 6d ). If the area with zero sensitivity is not considered, this proportion will be close to 70%. This shows that the increase in precipitation has a high probability of increasing CS and WR simultaneously. In 14.00% of the area, the change of direction in WR and CS was consistent with the change of direction in temperature (Fig. 6g) . However, a larger area (16.72%) showed that CS and WR trends are opposite after experiencing temperature changes in the same direction. In CS and HQ (Fig. 6e) , there was no apparent correlation between precipitation and ecosystem services change of direction, and the areas of the four states were relatively close. Exposure of ecosystem services to future climate scenarios Due to spatial heterogeneity of climate velocity (Fig. S9) , the pattern of exposure exhibited a heterogeneous pattern, and a high exposure was founded in the southwest part of the studied area (Nyingchi and Sannan), and some lakes like Qinghai Lake surrounding area (Figs. 7 and S10 ). Under the SSP1-2.6 scenario, the mean exposure of CS to precipitation was negative (-14.87 t/km 2 ), while to temperature was positive (5.71 t/km 2 ). Also, climate change will boost WR's effect and slightly inhibit the stable supply of HQ (mean exposure of 0.417 mm/km 2 for precipitation and 2.05 mm/km 2 for temperature). The distribution pattern of exposure was relatively similar to the different scenarios. However, climate change under the high-emission pathway (SSP5-8.5) tends to have greater gains or losses on ecosystem services (Fig. 7) . In CS, high exposure areas are mainly distributed in the East, South of the Tibetan Plateau, and surrounding areas of Qinghai Lake, while in WR, the high exposure pixels were found around the scattered lakes in the southwest and east of Tibetan Plateau. High exposure was also observed in the southeast corner (forested area) of the studied area in HQ. Except for CS, the other ecosystem services' exposure pattern is similar to temperature and precipitation. This matched with the distribution of ecosystem services sensitivity. To investigate the altitude discrepancies in the exposure of ecosystem services to climate change, a zonal average of negative and positive exposure per 50 m altitude was applied (Fig. 7) . We found that exposure features have a unimodal or multimodal variation with altitude, and there was a remarkable symmetry between positive and negative exposure. The largest exposure to temperature was observed at low altitudes (0-2000 m), while sporadic peaks were identified in higher elevation ranges ([ 2000 m). Compared to temperature, peak exposure to precipitation was observed at higher altitudes ([ 3000 m). In other altitude ranges, the peak intensity was not evident. The exposure was compared along the altitude gradient with the low-emission pathway (SSP1-2.6) as the baseline (Fig. S11) to identify the exposure differences under the various scenarios. A high-emission pathway (i.e., SSP5-8.5) indicates a higher intensity of exposure, whether positive or negative. As expected, the exposure under SSP5-8.5 is the highest, and the exposure intensity is close to 1.5 times higher than the baseline (SSP1-2.6) in most altitude ranges. The trajectories of exposure in SSP2-4.5 and SSP3-7.0 scenarios are relatively consistent, and the exposure intensity of these two scenarios is slightly higher than the baseline. Even in some altitudes (e.g., exposure of WR from 2000 to 4000 m altitude), the exposure intensity under these SSP2-4.5 and SSP3-7.0 is lower than the baseline intensity. In general, low-altitude areas face greater temperature exposure, while mid-high-altitude areas face greater precipitation exposure. The higher emission pathways increased exposure in most of the cases. Fig. 6 Relationship of ecosystem services' sensitivity to climate change. The first row is the relationship between the single ecosystem service (a: CS; b: WR; c: HQ) and two climate drivers (precipitation and temperature). The second and third row is the relationship between the two ecosystem services and a single climate driver (precipitation or temperature). The numbers inside and outside the box represent the area proportion (%) and the correlation coefficient of pixels in this state. All correlation coefficients pass the significance test (p \ 0.05). See Fig. 3 and Section ''Spatial association of ecosystem services' sensitivity to climate change'' for the legend description We also explored the conservation gaps between the current distribution of protected areas and the need to adapt to future exposure. According to the sign (negative/positive) of the exposure value to precipitation/temperature, all pixels were reclassified into eight states (Fig. S12) . We then compared the distribution of protected areas and pixels where future climate change may cause loss of ecosystem services (Fig. S13) . Protected areas are mainly distributed in the northwest part of the study area, while the pixels with negative exposure are mainly distributed in the southeast of the study area. For CS, under these four climate scenarios, protected areas can cover approximately 29% of the negative exposure area. However, for WR and HQ, only approximately 20 and 17% of the negative exposure area are covered by existing protected areas, respectively. The protected areas cover 38.5% of the Tibetan Plateau. The above ratios are all far less than Fig. 7 Exposure of ecosystem services to climate change. The first column is the exposure of ecosystem services (CS, WR, and HQ) to climate change (precipitation and temperature) under SSP1-2.6. The second, third column is the positive and negative exposure of ecosystems to climate change along the altitude gradient under four climate scenarios (SSP1-2.6, SSP2-4.5, SSP3-7.0, and SSP5-8.5), respectively. Different colours represent different climate scenarios, and the rectangular box in the figure is a close-up of local features of exposure. And Figs. S9 and S14 showed climate velocity and standard error estimate of climate velocity under four climatic scenarios 38.5%, indicating large conservation gaps for future exposure. Interpretation of the sensitivity of ecosystem services to climate change Climate change may have positive effects on some ecosystem services and negative in others. The magnitude of influence is very heterogeneous and dependent on the region and ecosystem services types (Mina et al. 2017) . Understanding the direction and magnitude of climate change allows human communities to anticipate the impacts of climate change better and adapt if necessary (Weiskopf et al. 2020) . Precipitation and temperature are fundamental driving forces for ecosystem processes change and vegetation distribution and productivity, affecting CS substantially (Zhao and Running 2010; Chen et al. 2014; Li et al. 2014) . Increasing precipitation will reduce water deficit stress, promote vegetation photosynthesis, increase vegetation growth (Piao et al. 2003; Zhang et al. 2016) , and therefore enhance CS. Additionally, warming can increase NPP by boosting photosynthesis rates and plant growth (Sullivan et al. 2008) . Our results showed that CS in Guoluo, Haibei, and Yushu had a positive sensitivity to precipitation, where CS and precipitation had a growth trend (Figs. 4, 5 and S1). However, in the large area of the Hengduan Mountains, the change of CS and precipitation had the opposite trend. This is because the decline of precipitation in this area is not significant. Therefore, the effects on CS may be limited. On the other hand, implementing various ecological restoration projects such as the Grain for Green Project (initiated in 2000) improved vegetation conditions in the Hengduan Mountains (Qu et al. 2018 ). This may have played an essential role in improving CS (Tao et al. 2016; Qu et al. 2018) . Yin et al. (2020) also confirmed that climate change had a limited impact on NPP improvement in this region. Overgrazing, sloping cropland reclamation, and forest harvest, due to population growth and urbanization, reduced CS (Yin et al. 2020) . Furthermore, soil erosion events were identified in this area due to cropland reclamation in slopes, affecting vegetation growth (Liu et al. 2010; Zhao et al. 2015) . The change in the direction of WR is primarily determined by the precipitation, as confirmed by the results (Fig. 4) . Evapotranspiration, which is highly dependent on temperature, is likely to increase with global warming, including in the Tibetan Plateau (Dou et al. 2015) . WR of Xigaze, Sannan, and Nyingchi showed a negative sensitivity to temperature, whereas precipitation and WR exhibited a decreasing trend and temperature an increasing one. WR showed a positive sensitivity to Golog, Aba, and Yushu's temperature, where WR, precipitation, and temperature all exhibit an increasing trend. On the other hand, the temperatures will not necessarily lead to an increase in evapotranspiration, according to the ''evaporation paradox'' (Brutsaert and Parlange 1998 ). An evapotranspiration declining trend was reported in Qilian Mountain and Qinghai provinces (Deng et al. 2017) . Therefore, in terms of WR, in the Tibetan Plateau, the temperature cannot fully offset the impact from the precipitation. Climate change influences HQ by shifting plant distribution, phenology, biological events, food availability timing, and landscape connectivity (Chen et al. 2011 (Chen et al. , 2018 Guo et al. 2015) . The impact of climate change is particularly evident in mountainous areas (Kelly and Goulden 2008) . This is consistent with the results observed in the Hengduan Mountains. Frequent human activities and the ecological engineering measures mentioned can weaken or strengthen resistance to climate change, making the sensitivity's direction of ecosystem services to climate change heterogeneous. Antagonism and synergies between ecosystem services' climatic drivers Temperature and precipitation may mutually weaken or jointly promote the same ecosystem services supply. In 27.78% of the regions, precipitation and temperature have a coordinated effect on CS, where the trends of CS, temperature, and precipitation are similar. According to the impact of climate change on CS, the influencing direction of temperature and precipitation on CS is the same. However, in most regions (47.45%), temperature and precipitation trends have an antagonist effect on CS. Among them, 29.95% of the regions showed that the trend of CS is consistent with the precipitation trend but different from the observed temperature trend, mainly in Nagqu. In this region, the gains from warming cannot offset the loss caused by the precipitation decrease. When temperature and precipitation have an antagonist relation, the CS direction in several areas is mainly influenced by the precipitation (29.95% vs. 18.10%). This is because the warming trend in the Tibetan Plateau is mainly observed in the winter (Xu et al. 2011; Liang et al. 2013) , and winter temperature has a limited impact on vegetation development (Harris et al. 2010) . Regarding WR, the trend observed matched with that of precipitation, occupying near 99% of the area of non-zero value (sensitivity or value of WR service). Among the non-zero value area, near half (47.17%) of pixels showed an antagonism between temperature and precipitation, especially in the east of the study area. The reduced correlation (0.11, p \ 0.05) showed an antagonist effect. Also, no clustering pattern was identified. We also explored the effect of the same climatic driving force on multiple services. For CS and WR, these two services' trends matched precipitation in 29.06% of the area. In 12.63% of the studied area, the WR and precipitation trend is the same, contrary to that observed in CS. This result is in line with the influence of precipitation. Increased precipitation can create conditions boosting CS and WR. Moreover, the correlation coefficient (0.43, p \ 0.05) between the two sensitivities indicated that a large increase in WR often accompanies the high growth of CS caused by precipitation. Compared to precipitation, relationship between temperature and the studied ecosystem services is complex. This might be due to that precipitation can offset the opposite effects of the temperature. Also, excessive warming can cause abnormal conditions, such as a decline in NPP due to enhancing the stress of water deficit and increasing physiological drought (Wu et al. 2011; Gu et al. 2017 ) and a decline of evapotranspiration (Peterson et al. 1995; Brutsaert and Parlange 1998; Roderick and Farquhar 2002) . The impact of climate change on HQ is mainly distributed in the southeast region. This may be related presence of protected areas that increase the resistance of HQ to climate stresses. There are several limitations in the analysis. Data with different resolutions were used, especially the resolution of future scenario climate (2.5 min) is quite different from other datasets (500/1000 m). These problems and limitations were identified in previous works (e.g., Gomes et al., 2021) . This study quantified ecosystem services' sensitivity to historical climate change and inferred the possible future exposure combined with future scenario data. Climate models are highly uncertain (e.g., Wootten et al. 2017 ) and the projection of their impacts on ecosystem services may inherit these of uncertainties. To minimize this, we used an average of 8 different future climate models to reduce the errors and increase the robustness of the work. Although we assessed the potential impact of future climate change on ecosystem services, we have not performed future ecosystem service simulations. For this, we consider that the simulation of ecosystem services often requires robust data and simulations. Some recent studies also confirmed our concerns (Shi et al. 2018; Ploton et al. 2020) . Furthermore, there is a slight inconsistency in the period involved in the ecosystem services' evaluation and sensitivity calculation due to the data availability (e.g., CS and WR are 2000 and 2001 . The above process may cause some uncertainty. We assumed that the effects of climate change are linear and instantaneous, although we explained it in Section ''Spatial association of ecosystem services' sensitivity to climate change'', lacked consideration of the non-linear effects and lag effects of climate change. The analysis of non-linear relationships and lag effects depends on the development of processbased models and services evaluation with high temporal resolution, which is out the scope of this study. We used ecosystem services indicators (e.g., NPP for CS) that explain the ecosystem service partially. Therefore, more models are needed to properly assess ecosystem services . For instance, NPP does not consider soil carbon sequestration, one of the largest carbon sinks (Guenet et al. 2018) . Also, WR was based on a runoff coefficient from literature review conducted by Ouyang et al. (2016) and not on field data from the study site. Even though the application of these coefficients was due to the data limitations that may produce uncertainties in the estimations and the studies considered in Ouyang et al. (2016) were extensive, field data is necessary to have better ecosystem services assessment. HQ model was validated using tree cover density (%) and above ground-biomass carbon. The best results were identified using above ground-biomass carbon, showing that the health of the ecosystems is more related to biomass carbon than with tree cover. The lack of a better validation can be attributed to the different base data resolution and years used (Table 1 ) and the same issues were identified in other studies (e.g., Gomes et al. 2021) . Other uncertainties may arise as a consequence of future national and global policies. It is clear that substantial efforts are being taken at the national (e.g., Yu et al. 2021 ) and international level (e.g., Cortes et al. 2021) to restore the ecosystems increase carbon sequestration. Different global strategies are trying to pave the road for a more sustainable world (e.g., UN Sustainable Development Goals, UN Decade on Restoration). Nevertheless, it is not certain that global crises will not occur in the future, such as the current COVID-19 pandemic that will affect global goals and shift political options to solve urgent issues (Barbier and Burgess 2020; Miksa et al. 2020; Tonne 2021) . Although our work has these limitations, we developed a spatially explicit analysis framework for the sensitivity and exposure of ecosystem services to climate change. We compared the exposure risks under different climate scenarios on the Tibetan Plateau, which provides an important reference for regional ecological security. Potential implications for exposure of ecosystem services to future climate change and management According to the definition, exposure is a simplified indicator to quantify the response of ecosystem services to future climate change. It captures potential climate disturbance intuitively and indicate the specific need for climate change adaptation planning. We found a certain similarity between the exposure pattern and the sensitivity pattern (Figs. 5, 7) . Complex terrain increases the fragmentation of exposure, especially the direction (negative or positive) that is very irregular. Combined with patterns of sensitivity and future exposure, adaptive strategies may arise. Therefore, from a management point of view, the areas with high ecosystem services value, high exposure, and high sensitivity such as the Nyingchi, Sannan, Hengduan Mountains, the surrounding area of Qinghai Lake should have special protection, since ecosystem services in these regions experience higher climate disturbance and deviate from the equilibrium state more easily. The rapid urbanization also may weaken resistance to future climate change. In this context, it would be relevant to limit urban expansion. The establishment of protected areas or the construction of ecological corridors is considered as an effective mitigation measure (Li et al. 2018a, b) . However, according to the current distribution of protected areas, their effect on ecosystem services is likely reduced. The current distribution has a big gap for ecosystem service conservation, and only 15% of priority areas with high ecosystem service and biodiversity value are covered by protected areas . Also, the results of exposure revealed that low proportion of the negative exposure are covered by existing protected areas for CS, WR, and HQ, respectively (Fig. S12 ). This shows that the existing protected areas have important conservation gaps for ecosystem services and future exposure, especially in the southeastern of the study area . The distribution of the protected areas needs to be optimized. Given interlaced distributions of positive exposure and negative exposure (Fig. S9) , some future management options to improve the ecosystems could be established such as grazing exclusion and ecosystem restoration (e.g., grasslands, forests) that can increase the resilience to future climate change (Li et al. 2018b; Wang et al. 2018) . In regions with low exposure but high sensitivity, the future climate may be stable. However, ecosystem services here may be sensitive to small disturbances. Therefore, management strategies, including planting trees with particular traits (e.g., vegetative reproduction and disease tolerance), can be adopted (Millar et al. 2007) . For the regions with high intensity in climate change velocity (e.g., Nagqu and Golmud), ecosystems will be subjected to different climate conditions than they experienced today (Fig. S9 ). In these regions, the delimitations of areas with a reduced human footprint that could serve as climate refugees to plants would be essential (Li et al. 2018a ). Considering the natural system's regulate ability (Pielke et al. 2011) , high input-output ratio, and relatively stable climate, no additional adaptation measures are required in low exposure and low-value areas such as Bayingolin. Different emission scenarios will also cause varying degrees of climate disturbance to ecosystem services (Fig. 7) . We found that a high-emission scenario (i.e., SSP5-8.5) indicates a higher intensity of exposure, regardless of positive or negative (Fig. 7) . The intensity of exposure under SSP5-8.5 scenarios is close to 1.5 times higher than the baseline (SSP1-2.6) in most altitude ranges. The increase of carbon dioxide emissions, temperature, and precipitation in highemission scenarios tend to promote part ecosystem services such as CS and WR. Therefore, high-emission scenarios do not necessarily lead to severe degradation of ecosystem services. For example, Fan et al. (2020) found that soil conservation service and carbon fixation service on the Tibetan Plateau will decline under low-emission scenarios, while high-emission scenarios can bring about the highest ecosystem services' sustainability. Our analysis also confirmed that, for CS and HQ, the number of pixels negatively affected by the climate is slightly reduced under highemission scenarios (Fig. S13) . However, considering other potential negative effects (e.g., ecosystems disruption) under high carbon emission scenarios, determining the optimal climatic pathway still needs to be treated with caution. Overall, understanding the impacts of climate change on ecosystem services is vital to achieve sustainable development goals (Pereira 2020). The present study contributes especially to the achievement of Goal 13 (Climate Action) and Goal 15 (Life on Land). Climate change induces positive or negative impacts on ecosystem services. The direction and intensity of the effect of climatic drivers (temperature and precipitation) on ecosystem services needs be spatialized. This will help to understand the geographic patterns and adapt better strategies to climate change. As demonstrated in this study, under the interaction of climate, human activities, altitude, and other factors, the sensitivity distribution patterns of these ecosystem services studied to climate change vary obviously. For CS, more than half of the regions showed a positive sensitivity to climate change, although the degree of sensitivity was not high. The regions where HQ was sensitive to climate change were mainly concentrated in the southern and southeastern parts. Increasing precipitation has a high probability of improving CS and WR as well. The existing protected areas have important conservation gaps regarding the disturbance of future climate change on ecosystem services, especially in the southeastern part of the study area. In practice, management priorities should be given to these high exposure areas and conservation gaps during the ongoing establishment of the national park group. Given the high universality of the sensitivityexposure framework, it can be applied to other regions to identify the vulnerability area for ecosystem service conservation practices. Furthermore, exploring the spatial association between the ecosystem services and their climatic drivers can provide a comprehensive and spatially explicit understanding, which helps optimize the response strategy in the context of global warming. Towards improved remote sensing based monitoring of dryland ecosystem functioning using sequential linear regression slopes (SeRGS) The human-environment nexus and vegetation-rainfall sensitivity in tropical drylands Developing China's ecological redline policy using ecosystem services assessments for land use planning Sustainability and development after COVID-19 Potential impacts of a warming climate on water availability in snow-dominated regions Understanding relationships among multiple ecosystem services Satellite-observed major greening and biomass increase in South China Karst during recent decade Hydrologic cycle explains the evaporation paradox The Pace of Shifting Climate in Marine and Terrestrial Ecosystems Rapid range shifts of species associated with high levels of climate warming The impact of climate change and anthropogenic activities on alpine grassland over the Qinghai-Tibet Plateau Long-term changes in the impacts of global warming on leaf phenology of four temperate tree Towards systematic analysis of ecosystem service trade-offs and synergies: main concepts, methods and the road ahead Where are global vegetation greening and browning trends significant? The value of the world's ecosystem services and natural capital Temporal-spatial dynamic change characteristics of evapotranspiration in arid region of Northwest China Spatiotemporal distribution of temperature in Gansu Province under global climate change during the period from 1956 to 2012 Scenario-based ecological security patterns to indicate landscape sustainability: a case study on the Qinghai-Tibet Plateau WorldClim 2: new 1-km spatial resolution climate surfaces for global land areas The climate hazards infrared precipitation with stations -A new environmental record for monitoring extremes Future scenarios impact on land use change and habitat quality in Lithuania The impacts of climate change on ecosystem structure and function Effects of climate warming on net primary productivity in China during 1961-2010 Impact of priming on global soil carbon stocks Responses of spring phenology in temperate zone trees to climate warming: a case study of apricot flowering in China The habitat concept and a plea for standard terminology Relationship between paired ecosystem services in the grassland and agro-pastoral transitional zone of China using the constraint line method Impacts of changes in climate and landscape pattern on ecosystem services Rangeland degradation on the Qinghai-Tibetan plateau: a review of the evidence of its magnitude and causes The human imperative of stabilizing global climate change at 1.5°C Effects of interannual climate variability on tropical tree cover Relationships of multiple landscape services and their influencing factors on the Qinghai-Tibet Plateau Climate change will affect the asian water towers Importance and vulnerability of the world's water towers Ecosystem services of the Baltic Sea: an assessment and mapping perspective An IPCC special report on the impacts of global warming of 1.5°C above pre-industrial levels and related global greenhouse gas emission pathways, in the context of strengthening the global response to the threat of climate change, sustainable development, and efforts to eradicate poverty. World Meteorological Organization, Geneva IPCC (2019) Summary for policymakers. Climate change and land: an IPCC special report on climate change, desertification, land degradation, sustainable land management, food security, and greenhouse gas fluxes in terrestrial ecosystems Non-linearity in ecosystem services: temporal and spatial variability in coastal protection Review on climate change on the Tibetan Plateau during the last half century The future value of ecosystem services: global scenarios and national implications The potential of ecosystem-based management to integrate biodiversity conservation and ecosystem service provision in aquatic ecosystems Climate variability rather than overstocking causes recent large scale cover changes of Tibetan pastures Natural and human impacts on ecosystem services in Guanzhong-Tianshui economic region of China Vulnerability of the global terrestrial ecosystems to climate change Effect of degradation and rebuilding of artificial grasslands on soil respiration and carbon and nitrogen pools on an alpine meadow of the Qinghai-Tibetan Plateau Increasing sensitivity of alpine grasslands to climate variability along an elevational gradient on the Qinghai-Tibet Plateau Regional vegetation dynamics and its response to climate-a case study in the Tao River Basin in Northwestern China Enhancing protected areas for biodiversity and ecosystem services in the Qinghai-Tibet Plateau Climate change in the Tibetan Plateau three rivers source region: 1960-2009 The willingness to pay for ecosystem services on the Tibetan Plateau of China Spatial patterns and driving forces of land use change in China during the early 21st century The velocity of climate change Synergies and trade-offs between ecosystem service supply, biodiversity, and habitat conservation status in Europe Diverse policies leading to contrasting impacts on land cover and ecosystem services in Northeast China The impact of interventions in the global land and agri-food sectors on nature's contributions to people and the UN sustainable development goals Future C loss in mid-latitude soils: climate change exceeds land use mitigation potential Ecosystem services and legal protection of private property. Problem or solution? Climate change and forests of the future: managing in the face of uncertainty Future ecosystem services from European mountain forests under climate change Spatial assessment of habitat conservation status in a Macaronesian island based on the InVEST model: a case study of Pico Island Development of a global evapotranspiration algorithm based on MODIS and global meteorology data Climate change's impact on key ecosystem services and the human well-being they support in the US Improvements in ecosystem services from investments in natural capital Ecosystem services response to urbanization in metropolitan areas: thresholds identification Ecosystem services in a changing environment Evaporation losing its strength Seasonal dynamics of terrestrial net primary production in response to climate changes in China Land use/land cover changes and climate: Modeling analysis and observational evidence Spatial validation reveals poor predictive performance of largescale ecological mapping models The third pole Non-linear groundwater influence on biophysical indicators of ecosystem services Spatial and temporal variability of future ecosystem services in an agricultural landscape Understanding relationships among ecosystem services across spatial scales and over time What drives the vegetation restoration in Yangtze River basin, China: climate change or anthropogenic factors? The cause of decreased pan evaporation over the past 50 years Non-linear temperature effects indicate severe damages to US crop yields under climate change Ecosystem service supply and vulnerability to global change in Europe Sensitivity of global terrestrial ecosystems to climate variability Increasing forest disturbances in Europe and their impact on carbon storage InVEST 3.2.0 user's guide. The Natural Capital Project Uncovering the relationships between ecosystem services and social-ecological drivers at different spatial scales in the Beijing-Tianjin-Hebei region Model structures amplify uncertainty in predicted soil carbon responses to climate change Adaptation, adaptive capacity and vulnerability Remote sensing the vulnerability of vegetation in natural terrestrial ecosystems The role of changing multiscale temperature variability in extreme temperature events on the eastern and central Tibetan Plateau during 1960-2008 Temperature and microtopography interact to control carbon cycling in a high arctic fen Protection and construction of the national ecological security shelter zone on Tibetan Plateau Recovery trends for multiple ecosystem services reveal non-linear responses and long-term trade-offs from temperate forest harvesting Managing for climate change on protected areas: an adaptive management decision making framework Assessment of ecological effect of the natural forest protection project in Southwest China Tonne C (2021) Lessons from the COVID-19 pandemic for accelerating sustainable development A framework for vulnerability analysis in sustainability science Grazing management options for restoration of alpine grasslands on the Qinghai-Tibet Plateau Climate change effects on biodiversity, ecosystems, ecosystem services, and natural resource management in the United States Continued warming could transform greater yellowstone fire regimes by mid-21st century Characterizing sources of uncertainty from global climate models and downscaling techniques Responses of terrestrial ecosystems to temperature and precipitation change: a meta-analysis of experimental manipulation High positive correlation between soil temperature and NDVI from 1982 to 2006 in alpine meadow of the three-river source region on the Qinghai-Tibetan Plateau Identification of the geographical factors influencing the relationships between ecosystem services in the belt and road region from 2010 to 2030 Prioritizing sustainable development goals and linking them to ecosystem services: a global expert's knowledge evaluation What drives the vegetation dynamics in the Hengduan Mountain region, southwest China: climate change or human activity? Loess Plateau: from degradation to restoration Estimating net primary production of natural grassland and its spatialtemporal distribution in China Drought-induced reduction in global terrestrial net primary production from The current situation of slope farmland and comprehensive management of soil erosion in Yunnan Province China Realizing the values of natural capital for inclusive, sustainable development: informing China's new ecological development strategy Publisher's Note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations Author contributions All authors contributed to the study. Idea for the article was originated by HT, ZW and PP. Study conception and design was realized by HT and ZW. Data analysis was performed by HT and PP. The first draft of the document was written by HT and all authors (HT, ZW, PP, HX, and CF) commented on previous versions of the manuscript. All authors read and approved the final manuscript. Data availability All datasets were generated from publicly available sources but are available from the corresponding author on reasonable request.