Assessing the vulnerability of ecosystems to climate change based on climate exposure, vegetation stability and productivity

Global warming has brought many negative impacts on terrestrial ecosystems, which makes the vulnerability of ecosystems one of the hot issues in current ecological research. Here, we proposed an assessment method based on the IPCC definition of vulnerability. The exposure to future climate was characterized using a moisture index (MI) that integrates the effects of temperature and precipitation. Vegetation stability, defined as the proportion of intact natural vegetation that remains unchanged under changing climate, was used together with vegetation productivity trend to represent the sensitivity and adaptability of ecosystems. Using this method, we evaluated the vulnerability of ecosystems in Southwestern China under two future representative concentration pathways (RCP 4.5 and RCP 8.5) with MC2 dynamic global vegetation model. (1) Future (2017–2100) climate change will leave 7.4% (under RCP 4.5) and 57.4% of (under RCP 8.5) of areas under high or very high vulnerable climate exposure; (2) in terms of vegetation stability, nearly 45% of the study area will show high or very high vulnerability under both RCPs. Beside the impacts of human disturbance on natural vegetation coverage (vegetation intactness), climate change will cause obvious latitudinal movements in vegetation distribution, but the direction of movements under two RCPs were opposite due to the difference in water availability; (3) vegetation productivity in most areas will generally increase and remain a low vulnerability in the future; (4) an assessment based on the above three aspects together indicated that future climate change will generally have an adverse impact on all ecosystems in Southwestern China, with non-vulnerable areas account for only about 3% of the study area under both RCPs. However, compared with RCP 4.5, the areas with mid- and high-vulnerability under RCP 8.5 scenario increased by 13% and 16%, respectively. Analyses of future climate exposure and projected vegetation distribution indicate widespread vulnerability of ecosystems in Southwestern China, while vegetation productivity in most areas will show an increasing trend to the end of twenty-first century. Based on new climate indicators and improved vulnerability assessment rules, our method provides an extra option for a more comprehensive evaluation of ecosystem vulnerability, and should be further tested at larger spatial scales in order to provide references for regional, or even global, ecosystem conservation works.


Background
Recent climate changes have caused widespread negative impacts on global ecosystems (Allen et al. 2010;Allen et al. 2015). Continued emission of greenhouse gases will cause further warming and long-lasting changes in all components of the climate system, increasing the likelihood of severe, pervasive and irreversible impacts on various human and natural systems (IPCC 2014). Therefore, assessing the vulnerability of ecosystems in the context of global climate change is one of the major components of current global change study. According to IPCC's definition, vulnerability to climate change is the degree to which a system is susceptible to, or unable to cope with, the adverse effects of climate change, including climate variability and extremes. It is a function of the character, magnitude, and rate of climate variation to which a system is exposed (exposure), its sensitivity, and its adaptive capacity (IPCC 2001). The complexity of ecosystems makes the relevant vulnerability assessment studies have no uniform method (Turner et al. 2003;Füssel and Klein 2006), and researchers have worked on this issue from multiple perspectives and spatial scales. For instance, Wu et al. (2007) and Zhao & Wu et al. (2014) used the temporal variation of net primary productivity (NPP) as the indicator of ecosystem vulnerability. Gonzalez et al. (2010) used vegetation change to assess ecosystem vulnerability. Watson et al. (2013) proposed an assessment based on a spatial analysis of the ecoregion's natural integrity (vegetation intactness) with its relative exposure to future climate change, which was adopted by subsequent studies (e.g. Segan et al. 2016). Eigenbrod et al. (2015) further introduced relative species richness to improve the assessment of ecosystem vulnerability.
The main shortcoming in present ecosystem vulnerability assessment works is that researchers tend to focus on one specific aspect of the ecosystem such as biome shift , vegetation intactness (Watson et al. 2013) or vegetation productivity (Wu et al. 2007), while less attention has been paid on a method that can provide comprehensive assessment for these critical aspects of ecosystem. Here we recommend a vulnerability assessing methodology that integrates ecosystem's net primary productivity, vegetation stability (a new indicator that combines vegetation intactness and vegetation change) and its exposure to climate change. In our analysis, the degree of 'exposure' was indicated by climate changes, while NPP, vegetation intactness and its temporal change together represent the 'sensitivity' and 'adaptive capacity' of ecosystems. Our method is in part an inheritance and integration of existing methods, and improvements have also been made in order to make a better assessment of ecosystem vulnerability.
In the study of Gonzalez et al. (2010), ecosystem exposure was quantified as the probabilities of climate change (temperature and precipitation), and its value was actually determined by the one that has larger magnitude of change. However, ecosystem was not affected by temperature or precipitation alone but by the combined effects of the two, thus their method may not fully characterize the degree of exposure of ecosystems to climate change. In addition, their analysis did not consider whether the trend of climate change is positive or negative, which should have contrasting effects on ecosystems. It is well known that climate change has both negative and positive impacts on ecosystems. For instance, studies have found increased productivity as a result of increased temperature and atmosphere CO 2 , when water is not limiting . In this situation, ecosystem should not be considered as exposed to the adverse effects of climate change. Moisture and temperature work together in determining terrestrial ecosystem dynamics (Heimann and Reichstein 2008), thus a more suitable way to depict exposure is to use indicators that can simultaneously reflect the changes of both. In this study, a moisture index (MI) was used to estimate exposure, which has been shown by a variety of studies that can well reflect the integrate effects of temperature and precipitation on long-term vegetation dynamics (Yang et al. 2005;Grundstein 2009;Fensham et al. 2015a;Fensham et al. 2015b).
In the study of Watson et al. (2013), vegetation intactness was defined as the proportion of natural vegetation that was not subject to human disturbance in a region. However, as the authors stated themselves, this metric did not take into account the influence of vegetation degradation or change of vegetation type. This means that, the vegetation intactness cannot reflect the change of vegetation type or productivity under changing climate. Therefore, in our method, vegetation (or biome) shift due to climate and vegetation intactness were combined into one new indicator (hereafter referred to as 'vegetation stability'), and then overlapped with the dynamics of vegetation productivity as a measure of ecosystem sensitivity and adaptation ability.
In recent decades, varies of ecosystem process models have emerged, which has greatly facilitated researchers in various fields, especially in the study of global climate change (Rosenzweig et al. 2014;Pacifici et al. 2015). Among these models, 'MC' is a dynamic global vegetation model (DGVM) that was developed around the year 2000. The first version of MC model (MC1) has long been proven that can provide reasonable predictions for ecosystem carbon dynamics and potential vegetation cover Lenihan et al. 2008), and more importantly, the model also has been increasingly used in recent ecosystem vulnerability assessments Eigenbrod et al. 2015). In this study, an updated version of MC model 'MC2' was used, which is more efficient and user-friendly compared to its ancestor. We utilized the model to predict the spatial and temporal variation of vegetation productivity and potential vegetation distribution in Southwestern China, a region with the most complex topography and the most diverse ecosystem but also suffered from severe human disturbance in China. We analyzed the future (2017-2100) vulnerability of ecosystems for the study area under two representative concentration pathways (RCP 4.5 and RCP 8.5) with our assessing method, and hopefully provided references for regional, or even global, ecosystem conservation and assessment works.

Study site
Our study region was located in the Southwestern part of China, including four provinces (Sichuan, Guizhou, Yunnan, and Guangxi), Chongqing municipality and the southeastern part of the Tibetan plateau, with a total area of 2.33 million square kilometers (Fig. 1). The topography and climatic conditions of the region are complex and have a variety of vegetation and species, including two of the world's important biodiversity hotspots: the Hengduan Mountains (at the junction of Tibetan Plateau, Sichuan and Yunnan) and the southern part of Yunnan (Myers et al. 2000;Mittermeier et al. 2011). However, long-term human disturbances and dramatic climate changes in recent decades have severely threatened ecosystems in the study area. In addition, due to the combined influence of the East Asian, the South Asian and the plateau monsoon (Bao et al. 2018), the impact of future climate change on ecological security of this region is highly uncertain. However, systematic assessment of climate change vulnerability for ecosystems across Southwestern China has seldom been conducted.

MC2 model simulation and validation
MC2 (short for MC version 2) is a new implementation of the MC1 dynamic global vegetation model Bachelet et al. 2001a), which has been proved to provide reasonable vegetation dynamics for many studies from regional Kim et al. 2018), national (Bachelet et al. 2003;Mills et al. 2015) to global scales Kim et al. 2017). MC2 can simulate potential vegetation distributions, and the movement of carbon, nitrogen, and water through ecosystems (Lenihan et al. 2008). The model requires five climate variables (monthly mean, maximum and minimum temperature; precipitation; vapor pressure) and five soil variables (soil depth; bulk density; clay, sand and Fig. 1 Map of the study area and its location in China rock fractions) to run three interacting modules of biogeography, biogeochemistry and wildfire.
Our simulation included two phases. First, MC2 was run under historical climate for validation purpose. Historical climate inputs required by the model were extracted from the CRU TS 4.01 dataset (Harris et al. 2014 (Zhao et al. 2005), which has been used in many studies for validation of carbon-cycling models (Tang et al. 2010;Zhang et al. 2005;Buddenbaum et al. 2015;Slevin et al. 2017). Average annual NPP of the study area from 2000 to 2015 was extracted and calculated from the dataset and regridded to 0.5°resolution, which was then compared with the mean annual NPP of the same period simulated by MC2.
Second, future vegetation distribution and vegetation productivity were simulated for the year of 2017-2100 under the two RCPs (RCP 4.5 and RCP 8.5). Future climate inputs came from the BCSD-CMIP5 Climate Monthly Projections dataset produced by the "Downscaled CMIP3 and CMIP5 Climate and Hydrology Projections" (Brekke et al. 2013;Maurer et al. 2014). This dataset provided downscaled 0.5°global "Bias-Correction Spatial Disaggregation (BCSD)" climate projections for 37 major global circulation models (GCMs) from the IPCC Fifth Assessment Report, from which we selected the outputs of BCC_CSM 1.1(m). The BCC_CSM 1.1(m) model was developed by the Beijing Climate Center, and has been proven to be able to provide preferable projections for climate change study in China ). Since the model did not provide vapor pressure projections required by MC2, future monthly vapor pressure was estimated using the method of Harris et al. (2014), which was also adopted by the CRU dataset to calculate vapor pressure: where VAP is vapor pressure (hPa), T min is monthly average daily minimum temperature (°C). These 0.5°-grid climate data were used to simulate the future vegetation productivity of study area. However, for vegetation distribution, we first interpolated the data to a resolution of 0.05°and then performed the simulation. The reason is that the calculation of vegetation stability requires a finer scale vegetation map to calculate the proportion of changes in vegetation types within each 0.5°grid (see details below).

Future climate exposures
As mentioned above, we used MI to indicate the future exposure of ecosystems to climate change, which was defined as: where MI is the moisture index (mm), P is annual precipitation (mm) and PET is annual potential evapotranspiration (mm). A negative MI value indicated a dry climate where precipitation was not sufficient to meet climatic water demands while positive values indicated precipitation supplies were in excess of water demands (Grundstein 2009). By adding up monthly precipitation of each year (using the same data for MC2 inputs), we calculated 2017-2100 annual precipitation for the study area under RCP 4.5 and RCP 8.5 scenarios. Since the BCC_CSM 1.1(m) model did not provided projections for future PET, we used annual PET calculated by MC2 model.
We then calculated the probabilities of MI change for both RCP 4.5 and RCP 8.5 using a method modified from that of Gonzalez et al. (2010): where erf(x) is the error function, C MI is the projected MI change of 2017-2100, and σ is the 2017-2100 standard deviation, thus P MI is the probability that C MI falls within a calculated number of standard deviations of the mean. However, in the study of Gonzalez et al. (2010), climate exposure was considered to be vulnerable as long as temperature (or precipitation) changes, no matter the trend of climate change is positive or negative, thus the value of C was actually the absolute value of the climate change rate. In our calculation, the value of C MI was directly used without taking the absolute value, and thus the value of P MI ranged from − 1 to 1, with the sign (positive or negative) indicating the direction of MI change. The exposure of future climate was considered vulnerable only when MI showed a decreasing trend (see Introduction). We classified P MI into five vulnerability classes based on IPCC treatment of uncertainty (IPCC 2007 and very high (P MI ≤ − 0.95).

Future ecosystem sensitivity and adaptation ability
We used vegetation stability (S v ) and ecosystem productivity as the measure of 'sensitivity' and 'adaptation ability' of ecosystems to future climate change. Vegetation stability is defined here as the proportion of intact natural vegetation that remains unchanged in an area over a period. This index integrates two previous indices that are important in assessing ecosystem vulnerability: (1) vegetation change under climate change ), (2) vegetation intactness (Watson et al. 2013). As shown in Appendix S1, the calculation of vegetation stability had two steps. In each 0.5°grid, we first excluded the area that is not covered by natural vegetation due to human disturbances from the remaining area. We extracted the land cover map of the study area from the GlobCover 2009 version 2.3, which had a spatial resolution of 300 m (Arino et al. 2008), and regridded to 0.05°resolution to be consistent with the MC2 modeled map of vegetation in the next step. Then we excluded all sub-grids that were subject to human disturbances in each 0.5°grid, including Cultivated Terrestrial Areas and Managed Lands, and Artificial Surfaces and Associated Areas (Watson et al. 2013).
The second step required further excluded the area that the vegetation type has changed due to climate change for the remain area in each 0.5°grid. To achieve this, a more detailed vegetation map was required to identify the proportion of areas where vegetation type changed in each 0.5°grid. Therefore, future (2017-2100) potential vegetation map was simulated by MC2 model at a resolution of 0.05°(thus there were a hundred 0.05°s ub-grids in each 0.5°grid). For each 0.05°sub-grid, we used the majority of vegetation types during each of the two periods (2017-2027 and 2090-2100) as the begin and end state of the future vegetation distribution, respectively. By comparing the vegetation types of the two states, all sub-grids changed in vegetation type in each 0.5°grid were then excluded. Thus, the remaining part in each 0.5°grid will be the proportion of intact natural vegetation that remains unchanged in the future, namely vegetation stability. The value of vegetation stability (S v ) ranged from 0 to 100% and was negatively related with the vulnerability of ecosystems. We also classified S v into five vulnerability classes: very low (S v ≥ 80%), low (60% ≤ S v < 80%), medium (40% ≤ S v < 60%), high (20% ≤ S v < 40%) and very high (S v < 20%). By integrating present land use (vegetation intactness) and future potential vegetation change, this approach could provide a more comprehensive assessment of the vegetation dynamics.
To evaluate the vulnerability of vegetation productivity under future climate change, we used the annual NPP of 2017-2100 simulated by MC2 model to calculate the probability of NPP changes (P NPP ) for each 0.5°grid under the RCP 4.5 and RCP 8.5 scenarios, with the same method as for P MI (Eq. 3), and also classified into five vulnerability classes.

Assessment of future ecosystem vulnerability
For each of the two future climate scenarios (RCP 4.5 and RCP 8.5), we determined the overlap of the vulnerability classes that derived separately from future climate (P MI ), vegetation dynamics (S v ) and productivity (P NPP ). As the flowchart showed (Fig. 2), S v and P NPP were first overlapped following the rules as shown in Fig. 3, which could prevent from under-or overestimation of the vulnerability of certain grids of very high or very low S v or P NPP . Then the result was further overlapped with the vulnerability classes for P MI .

Model validation
According to the Vegetation Regionalization Map of China, there were five major vegetation types in our study region (Fig. 4a), including alpine grassland and shrubland, subalpine forest and tropical monsoon forest in the west, and subtropical evergreen broadleaf forest, subtropical and tropical monsoon forest in the east part of our study region. This geographic pattern of vegetation distribution was largely similar to the current vegetation map predicted by the MC2 model (Fig. 4b), suggesting that MC2 was qualified to model vegetation patterns in Southwestern China.
MC2 also well simulated the spatial pattern of ecosystem productivity for the study area, Fig. 5 showed that, modeled NPP was well correlated with MODIS NPP (Slope = 0.89, R 2 = 0.79). These tests suggest that MC2 can be well used to evaluate ecosystem vulnerability to future climate change in Southwestern China.

Future climate exposures
The exposure of ecosystems to the adverse effects of future climate change, as assessed by P MI (Eq. 1), was shown in Fig. 6a and b. Under the RCP 4.5 scenario, the exposure in the study area was relatively moderate by the end of the century, with only 7.4% of the total area classified as 'high vulnerability', and no area revealed 'very high' vulnerability (Table 1). Medially vulnerable areas accounted for the largest area (58.95%), and the remaining 35% were occupied by areas with low or very low vulnerability. While under the RCP 8.5 scenario, more than half of the study area (57.47%) revealed high or very high vulnerability, with less than 7% of the area showed low or very low vulnerability. In addition, future climate exposure under both RCP 4.5 and RCP 8.5 showed a clear latitudinal pattern, climate vulnerability generally increased with latitude, indicating that the impacts of climate change were more drastic at high latitude and altitude areas.

Future ecosystem sensitivity and adaptation ability
The MC2 simulations showed that, future climate change will have contrasting impacts on potential vegetation distribution in Southwestern China between the RCP 4.5 and RCP 8.5 scenarios. Under RCP 4.5, each major vegetation types showed a significant northward shift in the future (Figs. 4c), when compared with current vegetation distribution ( Fig. 4a and b). For instance, the south edge of alpine grassland on the Tibetan Plateau will move from the current 30°N to about 32°N , and the reduced distribution area will be occupied by subalpine forests. The subtropical evergreen broadleaf and monsoon forest zone will also move northwardly for about 1°. However, all vegetation types under RCP 8.5 will move southward significantly (Fig. 4d). Almost all area north of 29°N on the Tibetan Plateau will be covered by alpine grassland and shrubland. Current tropical monsoon forests in the southwest part of the study area will disappear and be replaced by the subtropical evergreen broadleaf forest zone in the future.
The vegetation stability in most areas of Southwestern China will decline to varying degrees ( Fig. 6e and f), with nearly half of the area classified into high or very high vulnerability under both pathways (45.61% and 42.76%, respectively), especially in low altitude areas. However, the northwestern part of the Tibetan Plateau will keep being a region with high natural vegetation intactness and moderate vegetation shift (very low vulnerability).
As for the vulnerability of vegetation productivity (2017-2100), about 92% of the study area showed very low or low vulnerability under the RCP 4.5 scenario, with the remaining area classified into medium vulnerability (Fig. 6e). Under RCP 8.5, about 16% of the study area revealed medium vulnerability, which were mostly located in the areas east of 105°E. The rest west regions generally showed low or very low vulnerability (Fig. 6f).

Future ecosystem vulnerability
By overlapping the vulnerability layers of future climate exposure, vegetation stability and productivity vulnerability, we obtained a comprehensive assessment of ecosystem vulnerability for the study area ( Fig. 6g and h). Under the RCP 4.5 scenario, the majority of our study region showed varying degrees of vulnerability, with 39.5% of areas classified as low vulnerability and 57% as medium vulnerability. Only 0.5% of areas were marked as highly vulnerable and no 'very high' vulnerable areas exist (Table 1). Under the RCP 8.5 scenario, low and medium vulnerable ecosystems still accounted for the vast majority of the study area (10.4% and 70%, respectively). There was still no area been classified as "very high" vulnerable, but the proportion of highly vulnerable areas has increased significantly to 16% when compared with RCP 4.5 scenario.

Future ecosystem vulnerability in southwestern China
In this study, we conducted a comprehensive analysis of ecosystem vulnerability of Southwestern China based on future climate exposure, vegetation stability, and vegetation productivity under two representative concentration pathways (RCP 4.5 and RCP 8.5). The results showed that, future climate change has different degrees of impact on water availability, vegetation distribution and  Gonzalez et al. (2010). Abbreviation: VL, very low; L, low; M, medium; H, high; VH, very high productivity of ecosystems, and increased climate warming will lead to more areas with high or very high vulnerability ( Fig. 6g and h).
Vegetation stability is an index that integrates vegetation intactness and climate change induced vegetation change. Since vegetation intactness is the same under both future scenarios (RCP 4.5 and RCP 8.5), then the differences in vegetation stability ( Fig. 6c and d) is actually determined by the changes of potential vegetation distribution. Future climate will drive significant latitudinal shift in vegetation distribution, but the directions of movements can be different under different climate change scenarios ( Fig. 4c and d). We found that all vegetation types will significantly move northward under RCP 4.5, which was consistent with many studies that global warming will lead species and biomes to expand their ranges towards high latitudes (Bachelet et al. 2001b;Hickling et al. 2005;Lucht et al. 2006) or high altitudes (Kelly and Goulden 2008;Moritz et al. 2008). In addition, studies have also shown that climate warming and its impact on ecosystems were greater at higher latitudes and altitudes (Deutsch et al. 2008;Qin et al. 2009;Mountain Research Initiative EDW Working Group 2015). Thus it can be predicted that the magnitude of biome shift at higher latitudes should be larger than that at lower latitudes. As the results showed, the alpine grassland and shrubland on the Tibetan Plateau moved northward by more than 2°(about 220 km) under RCP 4.5, but the subtropical evergreen broadleaf and monsoon forest at relatively lower latitudes only moved by about 1° (Fig. 4c). Under the RCP 8.5 scenario, however, the biomes in Southwestern China did not move northward but all shifted southward (Fig. 4d). In recent years, increasing studies have shown that warm-induced drought is a major limiting factor for the growth and distribution of natural vegetation during global warming, and the distribution of vegetation are closely related to water availability Huber et al. 2011;Xie et al. 2019). In our analysis, water availability was  well quantified by the moisture index (MI). As the results showed, for most of the study area, a relatively moderate warming (RCP 4.5) will not laid severe impact on water availability (Fig. 6a). Consequently, it is natural that the biomes shift northward as a result of climate warming, as commonly reported. However, a drastic and continuous warming (RCP 8.5) will lead more than half of the study area under serious water deficit (Fig. 6b). In this case, water will be the limiting factor for vegetation distribution and lead the vegetation to move to southern areas with relative better water availability.
Although the future vegetation distribution can differ significantly from current distribution, the productivity of future vegetation in most areas is not very vulnerable under both RCPs (Fig. 6e and f). This indicates that future vegetation productivity in the study area may not decrease significantly but will increase or remain stable in most areas (especially for alpine grassland and subalpine forest on the Tibetan Plateau). This is consistent with many reports that climate warming revealed positive effect on vegetation productivity (Beck et al. 2011;Hudson and Henry 2009;Miller and Smith 2012). Interestingly, areas with relatively vulnerable productivity were mostly under the control of the southeast monsoon (East Asia monsoon) (Jiang et al. 2017), and were also areas under the most severe drought stress ( Fig. 6a and  b). This may suggest that the regions controlled by the southeast monsoon (e.g. Central and South China) may suffer from drier and warmer climate in the future, and a large area of subtropical evergreen broadleaved forest in eastern China is likely to experience significant productivity decline.

Advantages of our method
By comparing with similar studies assessing ecosystem vulnerability to climate change, our analysis builds on existing efforts in three important aspects.
First, we used an indicator (moisture index, MI) that comprehensively reflects the hydrothermal conditions to calculate future climate exposure. In the study of Gonzalez et al. (2010), climate exposure was quantified by separately calculating the probability of temperature and precipitation change, and then taking the larger value of the two to evaluate climate exposure. We also applied this method to our study area in the preliminary analysis. As shown in Fig. S1, the vulnerability of future climate exposure (2017-2100) under RCP 4.5 scenario for all grids were marked as "very high" without exception (Fig. S1c), and this result was almost all determined by temperature alone (Fig. S1a) since the future temperature of the study area will increase significantly. However, climate exposures calculated with moisture index showed different results ( Fig. 6a and b). As an integration of temperature and precipitation, moisture index has long been used to describe vegetation distribution (Mather 1978) and has been proved to be aided in the interpretation of possible mechanisms controlling the distribution of vegetations (Stephenson 1990), while changes in temperature alone may be inadequate for understanding distributional shifts of plant species ). This has also been confirmed in our research that the MIbased climate exposure could well explain the difference of the vegetation shift in the study area under different future climate change scenarios ( Fig. 4c and d).
Second, we integrated biome shift due to climate change into the "vegetation intactness" proposed by Watson et al. (2013). Although "vegetation intactness" provided a practical way of identifying potential natural refugia (Watson et al. 2013;Eigenbrod et al. 2015), it is still a conservative measure of ecosystem adaptive capacity, because it is purely related with human-induced land use change but did not consider the climate change-induced vegetation change. Vegetation change has long been recognized as a significant consequence of climate change (Hutyra et al. 2005), and it has therefore been used in many studies to assess the vulnerability of ecosystems (Villers-Ruiz and Trejo-Vázquez 1997;Dlamini 2011;Iglesias et al. 2018). However, since most studies on vegetation changes use models to simulate potential distributions of vegetation types, it is generally assumed that all areas are covered by natural vegetation, making it difficult to predict human-induced land use. In this study, we proposed a method to evaluating biome shift and vegetation intactness simultaneously (i.e. vegetation stability), which provides a more robust way for assessing the vulnerability of vegetation to climate change. In addition, as one of the most important components of ecosystem function, vegetation productivity has been added to our assessment framework, which has long been used for indicating the temporal dynamics of ecosystems in response to climate change (Beck et al. 2011) but not often been used for assessing ecosystem vulnerability together with vegetation intactness, biome shift and climate exposure (e.g. Watson et al. 2013;Eigenbrod et al. 2015). Third, in our calculation, climate exposure and vegetation productivity were considered to be vulnerable only when they declined significantly, and otherwise the grids were classified as "very low" vulnerability. Assessment of vulnerability should not only focus on the magnitude of the change itself, but also should consider how ecosystems were affected by these changes. Therefore, when conducting similar climate vulnerability analysis (especially for future research), careful consideration should be given to the indicators used and the method for vulnerability classification to avoid over-or underestimation.
Limitations and suggestions for future research MC2 model (and its ancestor MC1) has been used in many ecosystem vulnerability studies in recent years from regional to global scales Choi et al. 2011;Zhao and Wu 2014;Eigenbrod et al. 2015).
Here we also showed that the model can provide a good simulation of vegetation distribution and productivity for Southwestern China (Figs. 4 and 5). However, the commonly used DGVMs such as LPJ , IBIS (Foley et al. 1996), SEIB (Sato et al. 2007) and CLM (Levis et al. 2004) simulate potential but not realized vegetation distributions, as well as MC2. These models generally use climate variables as main inputs for determining the potential vegetation type and generally are difficult to simulate the influences of harvesting, management, pests and diseases (Peng 2000), as well as competition at individual level (Fisher et al. 2010;Quillet et al. 2010). Therefore, the simulated results of DGVMs must be verified using reliable data. In our results, a small number of grids were classified with high (0.91%) or very high (0.23%) productivity vulnerability under RCP 8.5, and these grids was irregularly distributed on the Tibetan Plateau (Fig. 6f). We think this is likely to be related to the soil inputs of the MC2 model. The study of Peterman et al. (2014) has showed that MC2 model is very sensitive to soil depth in simulations of carbon and hydrological variables. The soil data we used indeed had a small number of grids on the Tibetan Plateau where the soil depth was very shallow or even zero (for example, bare land), which may cause the model to produce unpredictable results in simulating the productivity of these grids.
When calculating the vegetation stability, we use vegetation intactness to reflect the impact of human activities on ecosystems. However, so far there still lacks a reliable method to predict land use and cover change (LUCC) under future climate change. Consequently, we followed the method of Watson et al. (2013) and assumed that vegetation intactness will remain unchanged in the future. Although many studies have proposed LUCC prediction models, most of them focused on urbanization (Hepinstall et al. 2008;Halmy et al. 2015), and so far there have been few cases applied to natural ecosystems at a large scale. At present, it is still very difficult to make a reasonable prediction of land use change on dozens or even hundreds of years, further studies are still required (Verburg et al. 2019).

Conclusions
In this study, we built on recent works and carried out an ecosystem vulnerability assessment method based on climate exposure, vegetation stability and vegetation productivity. We applied the method for Southwestern China and assessed the vulnerability of ecosystems under two future climate change scenarios (RCP 4.5 and RCP 8.5). As the climate continues to warm, over 95% of the ecosystems in the study area will be vulnerable to varying degrees. Current vegetation distribution will show significant latitudinal shift, but the direction of movement will be different due to the difference in water availability under the two scenarios. While in the meantime, the productivity of future vegetation will generally show an increasing trend in most areas under both RCPs. Compared with similar studies, we used a moisture index (MI) that integrates the effects of temperature and precipitation to calculate future climate exposure, and put forward "vegetation stability" by integrating vegetation distribution shift into the "vegetation intactness" proposed by Watson et al. (2013). Based on these improvements, our research provided an extra option for a more comprehensive assessment of ecosystem vulnerability, and have the potential to be used on lager spatial scales such as national or even global scales.
Additional file 1: Figure S1. Schematic of the calculation of vegetation stability. Here the large square represents a 0.5°grid and is divided into four sub-grids for example. In the actual analysis, each 0.5°degree grid was divided into 100 sub-grids. Figure S2. Vulnerability of climate exposure for the study area under future RCP 4.5 scenario, calculated based on 2017-2100 (a) annual mean temperature, (b) annual precipitation, and (c) overlap of (a) and (b) with the method of Gonzalez et al. (2010).

Funding
This study was supported by the National Key Research and Development Program of China (No. 2016YFC0502104, No. 2017YFC0503901), and the National Natural Science Foundation of China (No. 31870430).

Availability of data and materials
The datasets used and/or analyzed in this study are available from the corresponding author on request.
Ethics approval and consent to participate Not applicable.

Consent for publication
Not applicable.