Impacts of Future Climate Changes on Spatio ‐ Temporal Distribution of Terrestrial Ecosystems over China

: Understanding the response of terrestrial ecosystems to future climate changes would substantially contribute to the scientific assessment of vegetation–climate interactions. Here, the spatiotemporal distribution and dynamics of vegetation in China were projected and compared based on comprehensive sequential classification system (CSCS) model under representative con ‐ centration pathway (RCP) RCP2.6, RCP4.5, and RCP8.5 scenarios, and five sensitivity levels were proposed. The results show that the CSCS model performs well in simulating vegetation distribu ‐ tion. The number of vegetation types would increase from 36 to 40. Frigid–perhumid rain tundra and alpine meadow are the most distributed vegetation types, with an area of more than 78.45 × 10 4 km 2 , whereas there are no climate conditions suitable for tropical–extra ‐ arid tropical desert in China. Some plants would benefit from climate changes to a certain extent. Warm temperate–arid warm temperate zone semidesert would expand by more than 1.82% by the 2080s. A continuous expansion of more than 18.81 × 10 4 km 2 and northward shift of more than 124.93 km in tropical forest would occur across all three scenarios. However, some ecosystems would experience inevi ‐ table changes. More than 1.33% of cool temperate–extra ‐ arid temperate zone desert would contin ‐ uously shrink. Five sensitivity levels present an interphase distribution. More extreme scenarios would result in wider ecosystem responses. The evolutionary trend from cold–arid vegetation to warm–wet vegetation is a prominent feature despite the variability in ecosystem responses to cli ‐ mate changes.


Introduction
The interactions between terrestrial ecosystems and climate changes at the regional scale have long received widespread attention from ecologists [1][2][3][4][5][6][7][8][9][10]. The International Geosphere and Biosphere Programme (IGBP) has recognized the accurate prediction of global and regional climate changes and their possible impact on ecosystems as a major research objective [11]. Regional and global climate changes driven by anthropogenic emissions of greenhouse gas (GHG) have placed heavy pressure on terrestrial ecosystems through the processes of rising average air temperatures and interannual variability of rainfall in recent decades [12][13][14][15]. These anticipated effects pose a great threat to vegetation growth and survival in susceptible regions, especially climate-sensitive and climate-vulnerable vegetation [16][17][18][19]. Basically, mature vegetation is a synonym for the final successional stage of vegetation [20]. Potential natural vegetation (PNV) is defined as the expected state of mature vegetation in the current environmental conditions and the absence of human intervention [20,21], which has been widely explored to reflect the climate-vegetation interactions [3]. The study of PNV change is essential to predict the potential impacts of climate changes on terrestrial vegetation dynamics and their response to future climate changes at a large scale, which sets the starting point for vegetation-climate relationship research [1,22].
There are many methods to modeling PNV response to climate changes at the regional and global scales [23][24][25]. The representative models, such as the ʺHoldridge life zoneʺ (HLZ) model and the BIOME series models (BIOM1, BIOME2, BIOME3, and BI-OME4), have been widely applied to explore the PNV on various spatial scales, which enhance our understanding of the relationship between climate and vegetation patterns [5,[26][27][28][29][30][31]. Based on the biotemperature (BT), the annual precipitation (P), and the potential evapotranspiration rate (PER), the HLZ model sets up a coordinate system, which divides the world into 38 life zones and more than 100 types [32,33]. However, this system does not consider the difference in the quantitative standard between horizontal and mountain zones, and it is still controversial for taking 30 °C as the upper limit of biological temperature [3]. The BIOME series models, only used in rough biome simulations at a large scale, are equilibrium vegetation models, while the selection of vegetation function models and the evaluations of the performance of model-based simulations are subject to several limitations, both inherent in the models themselves and in the input data utilized to run the models [5,34]. The Lund-Potsdam-Jena dynamic global vegetation model (LPJ-DGVM), a vegetation-atmosphere model, represents large-scale terrestrial vegetation dynamics and provides a way to comprehensively examine the responses of vegetation ecosystems to climate changes, land processes, and atmospheric CO2 increases [29,30,35]. Nevertheless, these benefits have been largely constrained by the complex input parameters required and a lack of data for validation, and further improvement of the DGVM is therefore necessary [12,36]. The comprehensive sequential classification system (CSCS), based on the relationships of climate, soil, and vegetation, was proposed by Ren et al. [35,37,38] in accordance with moisture and heat conditions of the current environment. The model has been successfully applied to simulate biomes at multiple spatial scales, which could contain a greater amount of information. For example, it can reflect the genetic relationship and zonal rules and unify a standard classification [1,13,38,39].
The fifth report of the Intergovernmental Panel on Climate Change (IPCC-AR5) has predicted that the resilience and adaptability of many climate-sensitive ecosystems, by the mid-21st century and beyond, could suffer and experience an unparalleled challenge posed by climate changes coupled with associated disturbances [40]. There are various thermal and humidity zones in China, and the diversity in climate environments has bred rich vegetation types [40,41]. In the context of global climate changes, the vegetation in China would also undergo adaptive changes in distribution and growth [3]. Therefore, the projection of the PNV change is of great significance to alleviate and effectively adapt to the impacts of climate changes on terrestrial ecosystems in China. The representative concentration pathways (RCPs) implemented by the IPCC-AR5 provided a long time series, in the range of historical to the future, for global land areas at a very high spatial resolution [42]. Combined with the CSCS method, the data have been successfully applied to simulate climate changes and their effects on terrestrial ecosystems. Numerous studies have been conducted on PNV in China based on the CSCS model. For instance, based on the CSCS model, Xiu et al. [43] carried out a study on the distribution dynamics of PNV in China using meteorological station data , coupled with future climate data (the 2030s, 2050s, and 2080s) under the RCP scenarios issued by the IPCC-AR5. Du et al. [11] explored the spatio-temporal distribution of PNV by using meteorological station data  and climate data (the 2030s, 2050s, and 2080s) under RCP4.5 scenarios as inputs into the CSCS model, and they divided its sensitivity to climate changes into four levels. Nevertheless, it is with high confidence that the number and uniform distribution of meteorological stations are expected to both directly and indirectly affect the interpolation effect of data, with consequences of biased results. Recently, based on the CSCS model, in combination of climatic datasets in the period of last interglacial period, etc., and projected data in the 2050s and 2070s under RCP2.6, the PNV and its response to climate changes were analyzed by Ren et al. [44]. However, these previous efforts have been mostly based on a single future climate scenario or biased data provided by the interpolation method; few have tried to conduct future intercomparison studies on the spatio-temporal succession of PNV in China and its sensitivity to climate changes under different RCPs. In addition, detailed comparisons between the CSCS model, the HLZ model, and a map of actual vegetation at a national scale have only been implemented minimally. Therefore, considering climate changes and the possible extent of their impacts on vegetation ecosystems in China, covering long time series under different climate scenarios, is essential to fill the gaps in this field.
This study aims to simulate terrestrial vegetation under projected climate changes across China at a 30 arc-second resolution (~1 km resolution), in particular with the dynamics of the spatio-temporal distribution of PNV and its succession as a consequence of climate changes, coupled with its sensitivity under different future climate scenarios. The objectives of this study are (1) to validate the accuracy of the CSCS model in simulating PNV of China; (2) to produce long time series maps of PNV in China with a 30 arc-second resolution (~1 km resolution) under current conditions  and future periods, covering the 2030s (2020-2049), 2050s (2040-2069), 2070s (2060-2089), and 2080s (2070-2099) under RCP2.6, RCP4.5, and RCP8.5 scenarios based on the CSCS model; (3) to compare the spatio-temporal distribution of PNV in China under three future climate scenarios; (4) to provide great potential in providing insights into how vegetation ecosystems in China will respond to climate changes in the upcoming decades; (5) to propose several sensitivity levels of PNV in China to climate changes.

Study Area
There are various thermal and humidity zones in China, and extra-arid zones, arid zones, semiarid zones, subhumid zones, humid zones, and perhumid zones dominate different regions in China (Figure 1), and more humid characteristics in China are found successively from northeast to southwest. The diversity in climate results in various vegetation types, and has a remarkable influence on the growth and richness of vegetation [45]. Interactions between vegetation and climate changes would vary in each vegetation type under different climate scenarios.

Data
The four datasets used in this study are current climate data, CPIM5 scenario projected data, basic geographic data, and a map of actual vegetation, as shown in Table 1.  [42]. The data are employed to drive the CSCS model to simulate PNV under the current scenario. Additionally, the modeling outcomes are used to validate the effect of the CSCS model in simulating PNV at a national scale.
The administrative zoning map of China was utilized as a mask to extract the study area.
The vegetation map of China at a scale of 1:1,000,000 was used to validate the CSCS model. It reflected in detail the actual vegetation distribution status in China, and 43 original vegetation types were then aggregated into five broad vegetation types. The map in raster form at a 1 km resolution is an essential scientific material and an important basis for studying climate changes, biodiversity, and vegetation protection and monitoring.

The Potential Natural Vegetation Model
The CSCS model, proposed by Ren et al. [37], is a comprehensive sequential classification system of climate-soil-vegetation [11]. The system is composed of three levels: class, subclass, and type, which are established through zonal bioclimatic characteristics, namely, the grouping or clustering of units with similar moisture and temperature properties [1,38]. The theory is that vegetation distribution is affected by the bioclimatic conditions, edaphic conditions, and vegetation characteristics, while the bioclimatic conditions are the most stable [11]. Therefore, the class level, the basic unit, is relatively stable and mainly determined by bioclimatic conditions. Classification in terms of the same class is supposed to have a consistent interpretation despite the difference in geographic locations [12]. The subclass level has intermediate stability, and is determined by the edaphic conditions [38]. Subclasses are integrated into classes based on the indexes of moisture and temperature which capture the natural occurrence of vegetation ecosystems [38]. The type level, the "less stable" hierarchy, is based on vegetation characteristics [12,38].
In the class level, the two key inputs employed in the CSCS model are the index of moisture based on humidity (K) and the index of heat based on cumulative annual temperature above 0 °C (θ) (i.e., growing degree-days on a 0 °C base, GDD0) [55]. The formula is as follows: where MAP represents the mean annual precipitation (mm); θ is cumulative annual temperature above 0 °C; 0.1 is an empirical parameter, and K is the humidity. By considering the relationship between vegetation type, soil type, natural zone, agricultural zone, animal husbandry production type, and contour line of the K value, the humidity is divided into 6 levels [1] (Figure 2). According to the heat zone distribution, the heat is divided into 7 levels ( Figure 1). Fifty-six PNV types can be theoretically classified by 7 thermal and 8 humidity levels (6 humidity levels, 1 meadow level, and 1 swamp level) [1]. To further explore the dynamics of distribution of PNV in China, 56 classes were regrouped into 10 super-classes based on the zonal characteristics. (Table 2) [1]. The overall approach is illustrated in Figure 3.

The Holdridge Life Zone Model
The HLZ model divides the world into over 38 life zones to determine the distribution of potential terrestrial ecosystems based on biotemperature (BT), average total annual precipitation (P), and potential evapotranspiration ratio (PER) [28,33]. The specific calculation formula of each HLZ model factor is as follows: .
where Ti is monthly mean temperature (°C), with values below 0°C substituted with 0 °C and above 30 °C substituted with 30 °C; BT is mean annual biotemperature (°C); PET is potential evapotranspiration (mm); P is average total annual precipitation (mm); PER is potential evapotranspiration ratio; BT(i,j), P(i,j), and PER(i,j), respectively, represent the average annual biological temperature value, annual precipitation value, and possible evapotranspiration rate value of the i th row and the j th column grid in the spatial raster data; BTk, Pk, and PERk are the central climate indexes of the k th hexagon in the life zone; Dk(i,j) is the distance from the center of the hexagon in the k th Holdridge model to i th row and the j th column grid.

The CSCS Model Validation
To make direct comparisons between the CSCS model, the HLZ model, and a vegetation map in China at a scale of 1:1,000,000, the following PNV maps in China from 1970-2000 were employed. Firstly, we generated PNV maps with 1 km spatial resolution based on the CSCS and HLZ model to validate the accuracy of the CSCS model and make appropriate comparisons over the same period. To further assess the capability of the CSCS model and HLZ model in simulating PNV distribution on a national scale, we compared two PNV maps with a vegetation map of China at a scale of 1:1,000,000.
Considering that each PNV type in the HLZ model may not necessarily correspond to one or more PNV types in the CSCS model because of the significant differences between the two models in terms of classification indices and methods, it is crucial for direct comparisons (i.e., CSCS vs. HLZ, CSCS vs. a vegetation map, HLZ vs. a vegetation map) to merge classes. Therefore, after trial and error, we proposed a reclassification where the 10 types of the CSCS model, 38 categories of the HLZ model, and 43 vegetation types (Table A1) were aggregated into the same 5 broad categories based on previous work on accuracy verification of the CSCS model (Table 3) [43,56]. Considering that validation samples would involve the majority of the study area and validate the reliability and accuracy of simulation results based on the CSCS and HLZ models well, we randomly created 750 points across China in ArcGIS software ( Figure 1).   Table 2; (b) the class names are found in Table A1. The Kappa coefficient has been widely used in assessing model-simulated vegetation types and spatial distributions due to its advantages. It takes chance agreement into account, regardless of the number of types being compared in the maps [5,57]. We employed the Kappa statistic to evaluate the similarities and validate the accuracy between the two kinds of PNV maps on a national scale. The Kappa statistic is expressed by [58]:

Class Code (b) in the HLZ model Vegetation code (b) in the Vegetation Map
where N is the total number of samples; n denotes the number of rows and columns in the error matrix; Pii represents the individual entry for the row and column on the main diagonal of the constructed error matrix; Pi+ denotes the row total for each category i; P+i is the column total for each category i. In general, the Kappa statistic (k) ranges from 0.0 to 1.0, with 0.0 representing totally different patterns and 1.0 indicating complete agreement. The evaluation thresholds are: 0.0-0.2: no to poor agreement; 0.2-0.4: poor to fair agreement; 0.4-0.55: fair to good agreement; 0.55-0.70: good to very good agreement; 0.7-1.0: very good to perfect agreement [35].

Accuracy Verification Based on the CSCS Model
To verify that the CSCS model we employed performs well in predicting the vegetation distribution, the differences in the simulation results between the CSCS model and HLZ model are compared and studied. Figure 4a Table A1 in Appendix A.

Spatial Distribution of PNV in China under the Current Scenario
There were 36 types of PNV in China during the period of 1970-2000 (T0) (Figure 4a and Table A2 in Appendix A), indicating that six types of vegetation, such as the tropicalextra-arid tropical desert (VIIA), did not appear at all in China in T0. The three PNV types with a wide distribution area are frigid-perhumid rain tundra, alpine meadow (IF), cool temperate-humid forest steppe, deciduous broad-leaved forest (IIIE), and warm temperate-extra-arid warm temperate zone desert (IVA), with the areas of 149.43 × 10 4 km 2 , 92.32 × 10 4 km 2 , and 90.22 × 10 4 km 2 , respectively, covering nearly 15.6%, 9.61%, and 9.39% of the total national area. IF, the most extensively distributed vegetation, is extensively distributed in the cold and wet mountainous areas of the Qinghai-Tibet Plateau, the Tianshan Mountains, and the Altai Mountains. IIIE is centrally distributed on the south side of Xiaoxingʹanling, Daxingʹanling, and Qinghai Nanshan. IVA is mainly distributed in the Tarim Basin, Inner Mongolia Plateau, and Zhungeer Basin.

Spatial Variations of PNV in China under Future Climate Scenarios
The spatial distribution and areas of PNV in China in the 2030s, 2050s, 2070s, and 2080s (T1, T2, T3, and T4) under three RCPs are shown in Figures 5-7, and Table A2. Our maps indicate that there are 39 (T1, T2 of RCP2.6, and T1 of RCP4.5) or 40 (others) classes of PNV in China. It is worth noting that VIIA does not appear in any periods of the three RCPs, which is consistent with the result of Du et al. [11]. Combined with the outcome of the current scenario, we found that VIIA never appeared in China from T0 to T4 in the four RCPs. Therefore, there is no suitable habitat for VIIA in China. The modeled results showed consistency with Che et al. [3]. Meanwhile, frigid-extra-arid frigid desert, alpine desert (IA) only appears in T4 of the RCP2.6 scenario with 0.012 × 10 4 km 2 , and subtropical-extra-arid subtropical desert (VIA) is not observed in T1, T2, and T4 of RCP2.6 and T1 of RCP4.5. All the other 39 vegetation types are projected to be distributed in the five periods of the three RCPs; nevertheless, the area changes for each vegetation type show obvious differences.  Table A1 in Appendix A.  Table A1 in Appendix A.  Table A1 in Appendix A.
Specifically, three possible change trends of areas were explored across all periods under the three RCPs, including the expansion, reduction, and fluctuation types (Table 4). Under the three scenarios, warm temperate-arid warm temperate zone semidesert (IVB) would benefit from the future climate, and is predicted to increase continuously. By contrast, the area of cold temperate-extra-arid montane desert (IIA) and cool temperate-extra-arid temperate zone desert (IIIA) would decrease consistently, in which IIA would cover less than 0.01% of the total land area of China in T3 and T4 under RCP8.5. These vegetation types would suffer and experience irreversible changes caused by climate changes if the radiative forcing level reflects the higher-level emission scenario (RCP8.5) by the end of the century.  , VIIE  IID, IIF, IVB, IVD, VA, VC,  VIA, VID, VIIF   IID, IIF, IVB, IVC, IVD,  VIA, VIB, VIC, VIIC,  VIID, VIIE, VIIF   Reduction pattern  IIA, IIC, IIIA, VF  IIA, IIC, IIIA, IVE, VF  IA, IB, IC, ID, IE, IF, IIA,  IIIA, IIIE, VF   Fluctuation  pattern   Expansion of  fluctuation   IB, IE, IIB, IIE, IIF, IIIB, IIIF, IVA,  IVC Under the RCP2.6 scenario ( Figure 5 and Table A2), the greatest continuous decreases are projected for warm-humid deciduous-evergreen broad-leaved forest (VF), followed by cool temperate-extra-arid temperate zone desert (IIIA), of 0.78 × 10 4 km 2 and 0.48 × 10 4 km 2 per decade, respectively. By contrast, the distribution of warm temperate-arid warm temperate zone semidesert (IVB) and tropical-humid seasonal rain forest (VIIE) would expand continuously by 0.67 × 10 4 km 2 and 0.50 × 10 4 km 2 per decade, respectively. In addition, the area of tropical-semiarid savanna (VIIC) and tropical-subhumid tropical xerophytic forest (VIID) would increase, by 38% and 15.13%, respectively, relative to T0.
Under the RCP4.5 scenario ( Figure 6 and Table A2), among all vegetation types, the distribution of warm-humid deciduous-evergreen broad-leaved forest (VF) would consistently shrink the most, followed by warm temperate-humid deciduous broad-leaved forest (IVE), by 1.09 × 10 4 km 2 and 1.04 × 10 4 km 2 per decade, respectively. By contrast, the highest and the second highest increases are projected for warm-extra-arid subtropical desert (VA) and warm temperate-arid warm temperate zone semidesert (IVB) of 1.99 × 10 4 km 2 and 1.13 × 10 4 km 2 per decade. In addition, the highest increase in area would occur in VIIC by 170.23% relative to T0, and VIID is also predicted to expand by 111.48% relative to T0.
Under the RCP8.5 scenario (Figure 7 and Table A2), the largest continuous decrease in area would occur in frigid-perhumid rain tundra, alpine meadow (IF), followed by cool temperate-humid forest steppe, deciduous broad-leaved forest (IIIE) by 3.44 × 10 4 km 2 and 2.28 × 10 4 km 2 every 10 years, respectively. By contrast, the largest and the second largest consistent expansions are expected to happen in warm temperate-arid warm temperate zone semidesert (IVB) and the warm temperate-semiarid warm temperate typical steppe (IVC) by 2.28 × 10 4 km 2 and 1.36 × 10 4 km 2 per decade, respectively. It is worth noting that the areas of VIIC and VIID are estimated to increase by 341.93% and 277.31% relative to T0, respectively. Based on the analysis of RCP2.6 and RCP4.5, the results indicate that VIIC and VIID would be more sensitive to climate changes than the other vegetation types.

Geometrical Center Shift of Super-Classes for PNV in China
To further reflect the macroscopic nature of PNV in the spatial movement, the centers of super-classes for PNV in China are simulated. Considering the error from elevation, we calculated the ʺGEODESICʺ distance, which presented different changing trends in directions and distances under the three scenarios ( Figure 8 and Table 5).  Table 2. Table 5. Geometrical center shifting distance of super-classes for PNV in China in T0 (1970-2000), T1 (2030s), T2 (2050s), T3 (2070s), and T4 (2080s) under RCP2.6, RCP4.5, and RCP8.5.

Super-Class Code
RCP2.6 (km) RCP4.5 (km) RCP8.5 (km)  T0-T1 T1-T2 T2-T3 T3-T4 T0-T1 T1-T2 T2-T3 T3-T4 T0-T1 T1-T2 T2-T3 T3-T4   1 Tundra As for movement directions, under the RCP2.6 scenario, the gravity center of the semidesert would move toward the southeast from T0 to T4. This super-class would continuously expand toward the southeast; in particular, its increasing rate of area would accelerate in T0-T1 then slow down after T1. Additionally, tropical forest is estimated to significantly shift toward the north, with the distribution persistently expanding to the north. Among the shrunken super-classes for PNV, the movement trend towards the west would occur in tundra and alpine steppe, which is projected to gradually shrink toward the west. By contrast, warm desert is expected to move eastward, and its area would continuously decrease eastward.
Under the RCP4.5 scenario, both subtropical forest and tropical forest would move to the north, and their distributions are predicted to continuously expand toward the north. Meanwhile, both semidesert and warm desert would significantly shift to the south. Semidesert is estimated to move to the southeast from T0-T1 then southwest from T1-T4, and its distribution would expand from southeast to southwest successively. By contrast, warm desert presents a southwest trend, which would also persistently increase toward the southwest.
Under the RCP8.5 scenario, a movement trend towards the north is projected for tropical forest, and the super-class for PNV would continuously expand toward the north. Tropical desert grassland would move to the southwest, and is estimated to increase gradually toward the southwest. Specifically, savanna would be observed to move southeast from T0 to T1, then northeast from T1 to T4, and would expand southeast and northeast successively.
From the movement distances of each super-class from the PNV perspective, the largest movement distance of steppe in each period would occur in RCP2.6, followed by RCP4.5, reaching 1023.02 km and 733.52 km, respectively, whereas tundra and alpine steppe would move the least, only 79.5 km and 93.12 km, respectively. By contrast, the largest movement distance of savanna in each period is expected to happen in RCP8.5, reaching 1534.52 km, whereas the shortest movement distance is projected for temperate forest, at only 472.23 km.
In terms of the movement distance of each period, in RCP2.6, the longest movement distance from T1 to T2 would be observed in tropical desert at 269.42 km, whereas steppe would move the most in all the other periods, by 545.71 km, 152.58 km, and 201.3 km, respectively. By contrast, the shortest movement distances during all the periods are projected for warm desert, tundra and alpine steppe, cold desert, and semidesert, at only 52.25 km, 6.32 km, 3.73 km, and 7.5 km, respectively.
In RCP4.5, the longest and shortest shifting distances from T0 to T1 are projected for steppe and tundra and alpine steppe, at 493.02 and 493.02 km, respectively. Warm desert and tundra and alpine steppe from T1 to T2 would present the longest and shortest movement distances, at 256.58 km and 11.62 km, respectively. In the T2-T3 period, cold desert and tundra and alpine steppe are estimated to move the longest and the smallest distances, with moving distances of 175.92 km and 8.31 km, respectively. During the period of T3-T4, temperate humid grassland would have the longest moving distance, reaching 85.01 km, whereas warm desert would move the least, at only 5.94 km.
In RCP8.5, in the T0-T1 and T3-T4 periods, steppe and temperate humid grassland would present the largest and the smallest movement trends, at 496.96 km and 190.11 km, respectively, whereas the center of tundra and alpine steppe is projected to move the least during the corresponding period, with a movement distance of 67.14 km and 33.58 km, respectively. In the T1-T2 and T2-T3 periods, the longest movement distance would occur in savanna, reaching 814.22 km and 406.76 km, respectively, whereas the shifting distance of temperate forest and semidesert would be the smallest during the corresponding period, at only 200.98 km and 36.93 km, respectively.

Spatio-Temporal Pattern Evolution of Super-Classes for PNV in China
The spatio-temporal pattern evolution of super-classes for PNV across China in T4 relative to T0 is simulated under the RCP2.6, RCP4.5, and RCP8.5 scenarios (Figure 9). From each super-class from the PNV perspective, temperate forest would evolve the most across all scenarios, with areas of 45.05 × 10 4 km 2 , 96.38 × 10 4 km 2 , and 120.30 × 10 4 km 2 , respectively (Table A3), indicating that temperate forest would be more sensitive to climate changes than the other vegetation types. By contrast, the smallest evolution of the areas of both warm desert and savanna would occur in RCP4.5 and RCP8.5, and the absence of any succession is expected to happen in RCP2.6. Based on the results, these two super-classes for PNV are projected to not include climate-sensitive and climate-vulnerable vegetation. In addition, the total succession areas of super-classes for PNV in China are estimated across three scenarios. The largest succession of 441.80 × 10 4 km 2 would occur in RCP8.5, followed by RCP4.5, of 249.23 × 10 4 km 2 , whereas the smallest succession of 148.75 × 10 4 km 2 would be observed in RCP2.6, covering nearly 46.04%, 30.66%, and 15.50% of the total national area. The outcomes indicate that more intense climate changes would result in a more extensive response of terrestrial vegetation.
As for evolution types, super-classes for PNV would present different evolution trends under the three scenarios. Under RCP 2.6, the largest area, of 34.36 × 10 4 km 2 , of temperate forest would convert into subtropical forest, which amounts to 3.58% of the total national area, and is evaluated to be concentrated in the North China Plain, the middle and lower reaches of the Yangtze River Plain, the northern part of the Lingnan Mountains, the Yunnan-Guizhou Plateau, the southern foot of the Qinghai-Tibet Plateau, the Sichuan Basin, and the Loess Plateau. Under RCP4.5, cold desert is projected to evolve into warm desert in an area of 38.87 × 10 4 km 2 , accounting for 4.05% of the total national area, which would contribute the most succession among all ecosystems. The ecosystem would be mainly distributed in the Tarim Basin and Tianshan Mountains. Under RCP8.5, tundra and alpine steppe succession is observed to be highest in temperate forest, of 69.23 × 10 4 km 2 , accounting for 7.21% of the total land area of China, which would be distributed extensively in the Qinghai-Tibet Plateau, Hengduan Mountains, Kunlun Mountains, Himalayas, and Qilian Mountains.

Analysis of Sensitivity of PNV in China to Climate Changes
The concept of ecosystem sensitivity mainly refers to the possible extent of terrestrial ecosystems in response to impacts associated with future climate changes. Global and regional climate changes as a consequence of anthropogenic emissions and projections of the possible extent of the effects on terrestrial ecosystems of China are essential for CSCS ecosystems in China to adapt to climate changes. The climate-sensitive and climate-vulnerable regions are prone to generating ecological environment problems, which are the focus of ecological environment protection and restoration [44,60,61]. Therefore, the sensitivity of PNV in China is explored to serve as guidance in ecological protection and restoration work. The five sensitivity levels, non-sensitive area (0 changes), low-sensitivity area (one change), medium-sensitivity area (two changes), high-sensitivity area (three changes), and extreme sensitivity area (four changes), are proposed through comparing the change numbers for each grid of PNV maps in China from T0 to T4 based on previous work in the sensitivity analysis of PNV in China [3,11].
Generally, all the three RCPs agreed that the five sensitivity levels present an interphase distribution ( Figure 10). The results of merging the low-sensitivity, middle-sensitivity, high-sensitivity, and extreme sensitivity regions into sensitive regions suggest that sensitive regions are concentrated, whereas the non-sensitive regions are relatively dispersed, which is consistent with the outcomes of Che et al. [3]. Under RCP2.6, the extremely sensitive region has the smallest distribution of 15.81 × 10 4 km 2 , covering 1.65% of the total land area of China (Figure 11a and Table 6), and is projected to only be distributed in the Northeast Plain, North China Plain, and Inner Mongolia Plateau. By contrast, the most widely distributed area is the low-sensitivity zone, at 242.42 × 10 4 km 2 , accounting for 25.26% of the total national area, which would be distributed across the Yinshan Mountains, Loess Plateau, Junggar Basin, Tarim Basin, Kunlun Mountains, Qinghai-Tibet Plateau, Hengduan Mountains, and Yunnan-Guizhou Plateau. The corresponding sensitivity of PNV in China to climate changes presents the regional differences, which are mainly caused by the joint action of climate changes and regional topography. It is noteworthy that the total area of the sensitive areas would be 496.53 × 10 4 km 2 , covering 51.74% of the total national area, whereas the non-sensitive areas, at 463.09 × 10 4 km 2 , account for 48.26% of the total national area (Figure 11a,d, and Table 6).   (Figure 11b and Table 6), would occur in the extremely sensitive zones, which are expected to be distributed in the Northeast Plain, North China Plain, and southern Inner Mongolia. By contrast, a low-sensitivity area of 412.70 × 10 4 km 2 , amounting to 43.01% of the total national area, would be distributed the most extensively. The zone would be observed in the Tarim Basin, the Northeast Plain, the North China Plain and the Sichuan Basin, Qilian Mountains, Kunlun Mountains and Tianshan Mountains, Quasi Algor basin, the Yunnan-Guizhou Plateau, Taiwan Mountains, and the hilly mountainous area in southern China. In addition, the total area of the sensitive areas is estimated to be 598.10 × 10 4 km 2 , covering 62.33% of the total national area, whereas the non-sensitive areas would be 361.51 × 10 4 km 2 , amounting to 37.67% of the total land area of China. The total area of the sensitive area simulated in our study is similar to that in Che et al. [3], who assessed the sensitive area in China to be 64.10% of the total national area. It is worth noting that the total sensitive area would increase by 20.46%, whereas the non-sensitive areas are predicted to decrease by 21.94% under RCP4.5 relative to RCP2.6 ( Figures 10 and  11d, and Table 6).
Under RCP8.5, extremely sensitive zones, with the smallest distribution area of 7.63×10 4 km 2 , accounting for 0.79% of the total land area of China (Figure 11c,d, and Table  6), are projected to be distributed in the Northeast Plain and North China Plain. In contrast, the most widely distributed areas would occur in the low-sensitivity zone, at 460.00 × 10 4 km 2 , covering 47.94% of the total area of China. In addition, the total area of the sensitive areas would be 782.48 × 10 4 km 2 , covering 81.54% of the total area of China, and the non-sensitive area would be 177.14 × 10 4 km 2 , accounting for 18.46% of China. Moreover, the total areas of the sensitive area and the non-sensitive area under RCP8.5 are projected to increase by 57.59% and decrease by 61.75%, respectively, relative to RCP2.6 (Figure 11d and Table 6). Based on the above analyses, it is with high confidence that a higher sensitivity of the terrestrial ecosystem would occur as a consequence of more intense climate changes.
In this paper, the modeled sensitivity results showed consistency with the distribution of Chinaʹs ecologically vulnerable areas described by Du et al. [11]. The common basic characteristics of ecologically fragile areas are that the ecosystems have weak resistance to disturbance and readaptation to new habitats, a strong spatio-temporal volatility, significant edge effects, high environmental heterogeneity, and sensitivity to global climate changes. In ecologically fragile areas, both environmental and biological factors are in a critical state of phase transition, and the ecosystems are sensitive to climate changes.

Conclusions
Climate changes featuring consistently increasing air temperatures and altered rainfall patterns since the 1980s have become a consensus in the context of global climate changes. This study concludes that future climate impacts on terrestrial vegetation in China could occur by the end of the 21st century, and the extents of terrestrial ecosystems in response to future climate changes would vary greatly among different ecosystems.
The types of PNV in China are projected to increase from 36 (1970-2000) to 39 (2030s, 2050s in RCP2.6, and 2030s in RCP4.5) and 40 (others). The CSCS model can successfully simulate vegetation distribution at a national scale. Frigid-perhumid rain tundra, alpine meadow (IF) is the most widely distributed vegetation type by the 2080s across the three RCP scenarios, whereas there are no climate conditions suitable for tropical-extra-arid tropical desert (VIIA) in China. The simulations of the dynamics of terrestrial ecosystems indicate that warm temperate-arid warm temperate zone semidesert (IVB) would benefit from the warmer and wetter climate, with a continuous expansion by the 2080s across the three RCP scenarios at the alongside a consistent contraction of cold temperate-extra-arid montane desert (IIA) and cool temperate-extra-arid temperate zone desert (IIIA). As a consequence of a warming climate, tropical-semiarid savanna (VIIC) and tropical-subhumid tropical xerophytic forest (VIID) would respond to climate changes with rapid increases of an average of 183% and 134% in the 2080s relative to 1970-2000, respectively. Meanwhile, warm temperate-arid warm temperate zone semidesert (IVB) is expected to continuously increase the most in response to changes in future climate under the RCP2.6 and RCP8.5 scenarios, whereas under RCP4.5, the largest continuous expansion would occur in warm-extra-arid subtropical desert (VA). Among the three scenarios, RCP8.5 shows the highest change rate of PNV in China. A continuous expansion and northward shift in tropical forest would occur across all the three scenarios, whereas cold desert would suffer a consistent contraction by the 2080s. The largest and smallest shifts are projected for steppe and tundra and alpine steppe, respectively, under RCP2.6 and RCP4.5, whereas the largest and smallest movements would occur in savanna and temperate forest, respectively. In addition, the largest shift in distance is expected to happen in RCP8.5 and the shortest distance would occur in RCP2.6. Five sensitivity levels present an interphase distribution, and more intense climate changes would result in a higher sensitivity of the CSCS ecosystem in China.
All the evolution assessments of the CSCS ecosystems agree that future climate changes would benefit the succession and growth of some vegetation to a certain extent, and the most succession would occur in the higher emission scenario. However, some ecosystems in China would suffer inevitable changes. Although the greatest contribution to succession could be different under the three scenarios, the CSCS ecosystems would present an overall change trend from cold and arid to warm and wet. In general, this paper contributes to improving our understanding of the CSCS ecosystem processes in China to mitigate and effectively adapt to the impacts of future climate changes.
In future studies, it should be possible to further study the advantages and disadvantages of the CSCS model through detailed comparisons between the Holdridge life zone (HLZ) model and the dynamic global vegetation models (DGVMs), etc., and put forward optimization methods (e.g., modification of parameters). Acknowledgments: The authors appreciate the anonymous reviewers and the editors for their valuable comments for significantly improving this manuscript. We also thank the WorldClim database and the Climate Change, Agriculture and Food Security program for providing climatic data, the National Natural Resources Standard Geographic Service Network for providing the Chinese administrative division map, and the Resource and Environment Science and Data Center for providing the vegetation map of China at a scale of 1:1,000,000.

Conflicts of Interest:
The authors declare no conflict of interest.