The Greenness of Major Shrublands in China Increased from 2001 to 2013

Shrubs have been reported to expand into grassland and polar regions in the world, which causes complex changes in ecosystem carbon, nutrients, and resilience. Given the projected global drying trend, shrubs with their superior drought resistance and tolerance may play more important roles in global ecosystem function. Shrubland exists in all of the climate zones in China, from subtropical to temperate and high cold regions, and they occupy more than 20% of the land area. In this paper, we analyzed the spatiotemporal trend of MODIS (Moderate Resolution Imaging Spectroradiometer) EVI (Enhanced Vegetation Index) for six shrubland types in China from 2001 to 2013 and its relationship to intraand inter-annual regional climate dynamics. Existing literature reported that the vegetation index did not change significantly in China during 2000–2012. However, we found that the shrubland EVI in China increased significantly at a rate of 1.01 ˆ 10 ́3 EVI ̈ a ́1 from 2001 to 2013. Two major shrubland types (subtropical evergreen and temperate deciduous) and two desert types (high-cold desert and temperate desert) increased significantly, whereas subalpine evergreen shrubland decreased at a rate of  ́0.64 ˆ 10 ́3 EVI ̈ a ́1. We also detected a significantly lengthened growing season of temperate deciduous shrubland. The growing season length contributed significantly to the annual averaged EVI for temperate deciduous, subalpine deciduous and subtropical evergreen shrublands. Furthermore, the precipitation variation contributed more to the annual averaged EVI than the temperature. The year-round decrease in rainfall and the increase in temperature led to a significant reduction in the subalpine evergreen shrubland EVI. The enhancement of countrywide shrubland EVI may promote its contribution to the regional ecosystem function and its potential to invade grasslands.


Introduction
The prevailing analysis and synthesis of climate models [1,2] project that a persistently increasing trend of aridity dominates the global climate change in the 21st century despite a short period of cooling over the past 17 years.Frequent and severe droughts would put many forest ecosystems under significant moisture stress, with severe consequences for the ecosystem structure and function [3][4][5].Because of the increased droughts and heat waves, the global net primary productivity (NPP) has been found to decrease in the new century [6], and tree mortality has been found to increase [7] in many places around the world.Ecosystems with drought resistance and tolerance, such as shrubs [8], may become increasingly prominent in the future.
Shrubs have been reported to expand to arid, semiarid, and sub-humid regions [9][10][11][12] and to invade arid and semiarid grassland [8,13].The expansion of shrubs in grassland has complicated the consequences of ecosystem structure and function.For instance, the invasion of shrubs in arid grassland reduced ecosystem production and carbon sequestration [14].However, recent studies have indicated that shrub encroachment has positive effects on the ecosystem functions of herbaceous biomass [15], microbe diversity, and nitrogen mineralization [16].Based on the experiment in southeast Spain [17], shrub encroachment has been found to increase soil organic carbon, total nitrogen, and biological soil crusts.Instead of leading to desertification in previous studies, shrub encroachment is considered to reverse desertification in the Mediterranean grasslands [17].
The superior tolerance of shrubs to moisture stress can be attributed to their small leaf-to-sapwood-area ratio [18,19] and their deep rooting profile [20,21].Because of the lack of a main trunk and their multiple branches, the leaf-to-sapwood-area ratio of shrubs is on the order of 10 3 , whereas the ratio for arbores is within the order of 10 5 to 10 6 [22,23].The stomata of shrubs is more resistant to moisture stresses in both soil and air compared to that of other functional types [24][25][26][27].Therefore, shrubs should perform better than other functional types in the global drying trend.
The seasonality of the climate has been considered to be an important modifier of grass-shrub competition.Shrubs sensed warming in the spring earlier than grasses due to the shrubs higher buds in the air, and warm spring favored shrubs during leaf out [8].Warming tends to shorten the growing season of mono-cultivated crops; however, it has been found to prolong the growing season of multispecies ecosystems [40].The rainfall seasonality might be the most important driver, and shifting rainfall from summer to other seasons favors shrubs rather than grasses [41].
Shrublands occupy more than 20% of the land in China [42] and can be classified into six functional types: subtropical evergreen (EVGNST) in the south; temperate deciduous (DCDSTP) in the north and northeast; subalpine evergreen (EVGNMT) and subalpine deciduous (DCDSMT) in the east slopes of the Tibet; temperate desert (DSRTTP) in the north and northwest; and high-cold desert (DSRTHC) shrublands in the west of Tibet.EVGNST is dominated by subtropical/tropical evergreen broad-leaved shrubs mixed with a small deciduous component.DCDSTP and DCDSMT are dominated by deciduous broad-leaved shrubs.EVGNMT is composed of a majority of sclerophyllous evergreen broad-leaved shrubs and a small part of evergreen needle shrubs.DSRTTP is dominated by temperate subshrubs and temperate desert shrubs, while DSRTHC is dominated by high-cold cushion subshrubs.EVGNST (the second largest in area) and DCDSTP have the best climate conditions and thus the largest average EVI of 0.429 and 0.321, respectively.On the other hand, the two desert shrublands (DSRTTP and DSRTHC) have the smallest average EVI of 0.088 and 0.071, respectively, despite the greatest distribution area of 1,068,218 km 2 of the former.The two sub-alpine shrublands have relatively small distribution areas and moderate average EVI values.Thus, the countrywide average EVI and its temporal patterns are mostly contributed by EVGNST and DCDSTP.The contribution of shrublands to the NPP in China has been found to be comparable to the arbores forest [43].Shrublands in China expanded in this century because of the policy-guided shrub planting after abandonment of crops in northern China [44,45].While EVI in China has been found to increase in a non-significant manner [46], we expected that the shrubland EVI should have a stronger increasing trend than the overall EVI in China.
The objective of this study is to test three hypotheses: (1) the annual average EVI in China Shrublands significantly increased from 2001 to 2013; (2) intra-and inter-annual changes in temperature and precipitation greatly altered the growing season of shrublands; and (3) the significant increase in shrub EVI can be explained by the extension of the growing season and the favorable climate.We analyzed the annual and seasonal EVI at both pixel-based and shrubland type scales to derive the trends in annual/seasonal EVI and climate variables, as well as the trends of the growing season length, and quantified the relationship between the trend in EVI and the climate variations and the growing season length.

Data Acquisition and Preprocessing
We obtained the MODIS EVI that covers continental China over the period from 2001 to 2013, with a spatial resolution of 1 km and temporal resolution of 16 days, from the NASA Data Center [47].Firstly, we applied the approach suggested by Samanta et al. [48] to screen for good data.Secondly, we temporally filled the missing or unreliable EVI based on the quality assurance flags using a simple linear interpolation [49].Then, we aggregated the monthly EVI by applying the maximum value composite (MVC) method to the two images from each month.Finally, based on the 1:1 million vegetation map [42], we extracted the EVI for the subtropical evergreen mixed with a small part of deciduous (EVGNST), temperate deciduous (DCDSTP), subalpine evergreen (EVGNMT), subalpine deciduous (DCDSMT), temperate desert (DSRTTP), and high cold desert shrublands (DSRTHC).Furthermore, pixels with an annual mean value less than 0.05 were excluded to reduce the impact of bare soil and too sparsely vegetated pixels [50].
To explain the changes in EVI, we obtained the monthly temperature and precipitation records of 659 meteorological stations across China from the Meteorological Data Center for the same period of time [51].The monthly temperature and precipitation were then interpolated to all of the shrublands using Kriging.The monthly EVI, temperature, and precipitation maps were spatially averaged over the six shrubland types and then averaged seasonally and growing seasonally to obtain the times series of the EVI and the climate variables.Our study focused on the growing season (from April to October), which can be divided into three seasons: spring (April and May), summer (June, July and August) and fall (September and October) [50,52].To explore the effects of winter climate change on the starting days of the growing season, we defined winter as November of the prior year to March of the following year [53].

Growing Season Length Detection
To explore the temporal trend of the growing season for each shrubland type, the TIMESAT program was used to obtain the profile of the 16-day EVI time series using asymmetric Gaussian functions [54].TIMESAT program is a common tool for phenology detection by setting a threshold value to determine the start of the season and end of the season [30,55,56].The start of the season was defined as the point at which the EVI increased to a value that was set to 30% of the distance between the minimum and maximum of the left edge.Similarly, the end of the season was defined as the point at which the EVI decreased to 30% of the distance between the minimum and maximum of the right edge.The growing season length was extracted by calculating the days between the starting day and the ending day.

Pixel-Based Non-Parametric Trend Analyses of EVI
Mann-Kendall's test and Theil-Sen median slope were used to obtain the possible monotonic trend in the time series of the annual average EVI at the pixel level, which were widely used to determine vegetation variation [57][58][59].The Mann-Kendall's test [60] is a robust nonparametric significance test that calculates Z values based on the signs of changes in the time series.Z values greater than 1.96 indicates a significantly increasing trend, and a Z value smaller than ´1.96 is regarded as a significantly decreasing trend.
Mann-Kendall's test does not count the magnitude of changes between any two consecutive pairs of numbers in the time series.The Theil-Sen median estimator, β TS , is a robust simple linear regression that chooses the median of slopes calculated between any pair of data points in the time series, and hence compensates for the Mann-Kendall's test in the quantitative aspect.This method is appropriate for assessing the rate of change in short or noisy time series because it is robust to outliers.

Regression with Serial-Correlation for Trends in EVI, Climate Variables and Growing Season
A more quantitative approach to time series analysis is the generalized least squares regression, which is capable of deriving the quantitative trend together with the serial correlation in the time series.Generalized least squares regression is applied to annual/seasonal EVI, temperature and precipitation, and EVI-derived growing season length in this study.To avoid spatial autocorrelation, we averaged these variables over the six shrubland types and applied the generalized least squares regression to the lumped data.We tested whether there was a significant temporal trend and quantified the strength of the trends in EVI, temperature and precipitation, as well as the growing season length of the growing seasons.
The seasonally and annually averaged EVI, climate, and growing season length variables were regressed on time using the following equation: where y is any of the above variables, x stands for year, and a is an intercept.Slope k stands for the trend.The analysis of EVI, seasonal climate variables and growing season length will allow us to quantify the changes in seasonal profiles of these variables.

Regression of Annual Average EVI on Growing Season Length and Climate Variables
The Kriging interpolated precipitation and temperature may not reflect the true water supply and temperature of the scattered shrublands.However, the temporal trends of the interpolated climate variables should be close to those of the shrublands.
To quantify the effects of climate variations on the growing season, we regressed the growing season length on the seasonal precipitation and temperature.The annual averaged EVI (April to October) were regressed on the growing season length and selected seasonal climate variables to establish dependence of EVI on the growing season length and climate.The regression equation can be expressed as follows: y " a `k1 where y is the annual averaged EVI, x represents the seasonal climate variables and growing season length.Regression coefficients are analyzed, interpreted, and discussed in the context of vegetation responses to regional climate dynamics.We used the Generalized Least Squares regression included in the "nlme" package of the free statistical software R [61] for all of the regression analyses.

Spatial Variation in the Interannual Trend of EVI at the Pixel Level Using Non-Parametric Analyses
The "significantly improved" and "insignificantly improved" classes (Figure 1, Table 1) occupy 14.67% and 50.23% of the total shrubland area of China, respectively, and are mostly distributed in the subtropical evergreen and temperate deciduous shrublands.The "significantly degraded" and "insignificantly degraded" classes occupy 3.36% and 30.88%, respectively, and are mostly located in the two subalpine shrublands, and the northwest of temperate desert shrubland.Moreover, 0.86% of the area is "unchanged".

Spatiotemporal Trends of Annual and Seasonal EVI at the Shrubland-Type Level Using Regression
The regression-obtained trend coefficient, k (Table 2), indicates that the two major shrubland types, subtropical evergreen (EVGNST) and temperate deciduous shrublands (DCDSTP), showed significant and similar increasing trends with k = 1.88 × 10 −3 EVI•a −1 and 1.89 × 10 −3 EVI•a −1 , respectively (p < 0.01), over the 13-year period.The high cold desert (DSRTHC) and temperate desert shrublands (DSRTTP) also exhibited an increasing trend coefficients with k = 0.36 × 10 −3 EVI•a −1 and 0.55 × 10 −3 EVI•a −1 , respectively (p < 0.05).Although DSRTTP EVI increased on the whole, we should note that a large area decreased in the northwest of DSRTTP.In contrast to the increasing trends of the above shrubland types, the two subalpine shrubland types showed a decreasing trend.Countrywide, shrublands exhibit an increasing trend with k = 1.01 × 10 −3 EVI•a −1 (p < 0.01).The finding is approximately consistent with those of the pixel-based analysis (Figure 1).

Spatiotemporal Trends of Annual and Seasonal EVI at the Shrubland-Type Level Using Regression
The regression-obtained trend coefficient, k (Table 2), indicates that the two major shrubland types, subtropical evergreen (EVGNST) and temperate deciduous shrublands (DCDSTP), showed significant and similar increasing trends with k = 1.88 ˆ10 ´3 EVI¨a ´1 and 1.89 ˆ10 ´3 EVI¨a ´1, respectively (p < 0.01), over the 13-year period.The high cold desert (DSRTHC) and temperate desert shrublands (DSRTTP) also exhibited an increasing trend coefficients with k = 0.36 ˆ10 ´3 EVI¨a ´1 and 0.55 ˆ10 ´3 EVI¨a ´1, respectively (p < 0.05).Although DSRTTP EVI increased on the whole, we should note that a large area decreased in the northwest of DSRTTP.In contrast to the increasing trends of the above shrubland types, the two subalpine shrubland types showed a decreasing trend.Countrywide, shrublands exhibit an increasing trend with k = 1.01 ˆ10 ´3 EVI¨a ´1 (p < 0.01).The finding is approximately consistent with those of the pixel-based analysis (Figure 1).
We calculated the relative changes (∆R in Table 2) as 13k/EVI 2001 , where EVI 2001 is the mean EVI of the growing season in 2001.The temperate desert (DSRTTP) shrubland showed the maximum relative increase of 8.86%, whereas the subalpine evergreen (EVGNMT) shrubland showed the maximum relative decrease of ´2.97%.The significant relative changes of the other types are between 5.90% and 8.10%, and the resultant countrywide relative change is 5.93% (Figure 2a).We calculated the relative changes (ΔR in Table 2) as 13k/EVI2001, where EVI2001 is the mean EVI of the growing season in 2001.The temperate desert (DSRTTP) shrubland showed the maximum relative increase of 8.86%, whereas the subalpine evergreen (EVGNMT) shrubland showed the maximum relative decrease of −2.97%.The significant relative changes of the other types are between 5.90% and 8.10%, and the resultant countrywide relative change is 5.93% (Figure 2a).Seasonal EVI showed more diverse trends among the six types (Table 2).The EVI of the temperate deciduous and subtropical evergreen shrublands showed a significantly increasing trend in all of the three seasons.Similarly, the high-cold desert and temperate desert shrublands showed a growing trend in all of the three seasons but was only significant in the spring and summer.However, the two subalpine shrublands showed a decreasing trend in almost all of the seasons except the spring of the deciduous shrublands, but only the trend in the fall of the subalpine evergreen shrubland was significant.Furthermore, the countrywide shrublands had significant increasing trends in all three seasons (Table 2).

Spatiotemporal Trends of Annual and Seasonal Precipitation and Mean Temperature
The annual precipitation of the four shrubland types showed increasing trends (Table 3).The trends of DCDSTP and DCDSMT were significant at the rate of 6.43 mm•a −1 (p < 0.05) and 2.80 mm•a −1 (p < 0.1), respectively.Moreover, DSRTHC and DSRTTP exhibited increasing trends of 0.83 and 0.57 mm•a −1 , respectively, but these were not significant.In contrast, EVGNST and EVGNMT had decreasing trends in annual precipitation at −3.80 and −3.70 mm•a −1 , respectively.Five shrubland types, DSRTHC, DSRTTP, EVGNMT, DCDSMT and EVGNST, showed increased annual mean temperature over the 13-year period, in which DSRTHC, DSRTTP, EVGNMT and DCDSMT had significant trends of 4.72, 3.36, 5.50, and 5.56 × 10 −2 °C•a −1 , respectively.In contrast, the annual mean temperature of DCDSTP exhibited a decreasing trend, but this was not significant.The countrywide Seasonal EVI showed more diverse trends among the six types (Table 2).The EVI of the temperate deciduous and subtropical evergreen shrublands showed a significantly increasing trend in all of the three seasons.Similarly, the high-cold desert and temperate desert shrublands showed a growing trend in all of the three seasons but was only significant in the spring and summer.However, the two subalpine shrublands showed a decreasing trend in almost all of the seasons except the spring of the deciduous shrublands, but only the trend in the fall of the subalpine evergreen shrubland was significant.Furthermore, the countrywide shrublands had significant increasing trends in all three seasons (Table 2).

Spatiotemporal Trends of Annual and Seasonal Precipitation and Mean Temperature
The annual precipitation of the four shrubland types showed increasing trends (Table 3).The trends of DCDSTP and DCDSMT were significant at the rate of 6.43 mm¨a ´1 (p < 0.05) and 2.80 mm¨a ´1 (p < 0.1), respectively.Moreover, DSRTHC and DSRTTP exhibited increasing trends of 0.83 and 0.57 mm¨a ´1, respectively, but these were not significant.In contrast, EVGNST and EVGNMT had decreasing trends in annual precipitation at ´3.80 and ´3.70 mm¨a ´1, respectively.Five shrubland types, DSRTHC, DSRTTP, EVGNMT, DCDSMT and EVGNST, showed increased annual mean temperature over the 13-year period, in which DSRTHC, DSRTTP, EVGNMT and DCDSMT had significant trends of 4.72, 3.36, 5.50, and 5.56 ˆ10 ´2 ˝C¨a ´1, respectively.In contrast, the annual mean temperature of DCDSTP exhibited a decreasing trend, but this was not significant.The countrywide shrublands exhibited an increasing trend for the mean temperature at a rate of 2.87 ˆ10 ´2 ˝C¨a ´1 (p < 0.1) but a decreasing trend for annual precipitation (Table 3, Figure 2b,c).The annual precipitation of EVGNST exhibited a trend of shifting to the fall with increased fall precipitation but decreased precipitation in other seasons.It also exhibited a significantly increased summer temperature.The precipitation of DCDSTP increased from the spring to fall, which is consistent with the significantly enhanced EVI.Moreover, the mean temperature exhibited a decreasing trend in the spring fall and winter but increased in the summer.
The precipitation of EVGNMT decreased for all of the seasons, and the temperature increased from the spring to winter.The decreasing rainfall may have caused the significantly decreased EVI.For DCDSMT, precipitation and temperature increased for all of the seasons except winter temperature.The increased spring temperature and summer precipitation and the decreased summer temperature may be responsible for the significantly increased spring and summer DSRTTP EVI values.Similarly, the increased rainfall and temperature for DSRTHC may have been responsible for the significantly increased spring and summer EVI.

Spatiotemporal Trend of the Growing Season Length and its Relationship with Climate
The application of the TIMESAT successively detected the growing season length for four shrubland types (Table 4).The effort failed for the two desert shrubland types, largely because of the small EVI and too noisy EVI values.Subtracting the starting day from the ending day yielded the growing season length for the four shrubland types.The temporal trend of the growing season length was also obtained based on generalized least squares regression.Only the DCDSTP growing season was significantly increased, by 0.46 day¨a ´1 (Table 4).The growing season was shortened for EVGNMT but was extended for DCDSMT and EVGNST.The altered growing season length may have contributed to the changes in annual average EVI.
The seasonal climate significantly affected the growing season length (Table 5).The DCDSTP growing season length was significantly increased with the spring temperature and precipitation, the summer temperature and the fall precipitation.Temperature in the summer and fall and precipitation in the winter, spring, and summer prolonged the EVGNMT growing season.However, the DCDSMT growing season length was significantly increased with the spring temperature and the fall precipitation.The growing season length of EVGNST, which is a warm and wet area, was prolonged by the winter temperature but was shortened by the spring precipitation and the fall temperature.

Regression of Annual Averaged EVI on Growing Season Length and Climate Factors
The model-predicted annual average EVI versus observed ones (Figure 3) indicates that the data points are approximately evenly scattered on the two sides of the 1:1 diagonal lines for the six shrubland types.Therefore, the model performed reasonably well in predicting the shrubland EVI.

Regression of Annual Averaged EVI on Growing Season Length and Climate Factors
The model-predicted annual average EVI versus observed ones (Figure 3) indicates that the data points are approximately evenly scattered on the two sides of the 1:1 diagonal lines for the six shrubland types.Therefore, the model performed reasonably well in predicting the shrubland EVI.Climate factors showed different effects on annual average EVI for the different shrubland types (Table 6).The annual average EVI of EVGNST, DCDSTP, and DCDSMT significantly increased with the growing season length.The annual mean EVI of two desert shrubs increased with precipitation in the summer, and we also detected a positive relationship between DSRTTP and its spring precipitation.The DCDSTP EVI significantly increased with the fall temperature and the summer precipitation but decreased with the spring precipitation.EVGNMT EVI responded to the summer temperature negatively but positively to fall precipitation, whereas DCDSMT EVI was the reverse.Furthermore, EVGNMT EVI decreased with the summer precipitation and DCDSMT EVI decreased with the spring, fall temperature and the fall precipitation.The EVGNST EVI significantly increased with the spring, summer temperature and the spring, fall precipitation.Table 6.Regression coefficients of annual average EVI on seasonal mean temperature and precipitation.P values were based on one-tailed hypothesis, and "*", "**", and "***" are for p-values less than 0.1, 0.05, and 0.01, respectively.Climate factors showed different effects on annual average EVI for the different shrubland types (Table 6).The annual average EVI of EVGNST, DCDSTP, and DCDSMT significantly increased with the growing season length.The annual mean EVI of two desert shrubs increased with precipitation in the summer, and we also detected a positive relationship between DSRTTP and its spring precipitation.The DCDSTP EVI significantly increased with the fall temperature and the summer precipitation but decreased with the spring precipitation.EVGNMT EVI responded to the summer temperature negatively but positively to fall precipitation, whereas DCDSMT EVI was the reverse.Furthermore, EVGNMT EVI decreased with the summer precipitation and DCDSMT EVI decreased with the spring, fall temperature and the fall precipitation.The EVGNST EVI significantly increased with the spring, summer temperature and the spring, fall precipitation.Vegetation indices in China increased significantly until the mid-1990s, after which the trend was weakened or even reversed according to previous reports [50,52,62].Similar changes were observed in North America [63].Specifically, Peng et al. [50] found that the June to August NOAA AVHRR (Advanced Very High Resolution NDVI (Normalized Difference Vegetation Index) over china decreased at a rate of ´0.2 ˆ10 ´3 a ´1, the April to May NDVI increased at 0.1 ˆ10 ´3 a ´1 and the growing season (April to October) NDVI increased at 0.4 ˆ10 ´3 a ´1 since 1990s.Although NOAA AVHRR NDVI and MODIS EVI are not directly comparable, the relationship between MODIS EVI and NOAA AVHRR NDVI can be described as EVI = 0.92 NDVI + 0.001 with R 2 of 97% for EVI in the range from 0.1 to 0.85 [64].The equation indicates that the increasing/decreasing EVI trend would be translated into a stronger corresponding NDVI trend.Therefore, we can translate the trends of summer EVI of the two major shrubland types (temperate deciduous and subtropical evergreen) and the countrywide shrublands in this study into 2.27 ˆ10 ´3, 1.21 ˆ10 ´3, and 1.01 ˆ10 ´3 a ´1 of NDVI, respectively.Furthermore, the countrywide EVI trend (April to October) can be translated to NDVI as 1.09 ˆ10 ´3 a ´1.Hence, the trend of the countrywide shrublands is not only positive and significant but also much stronger than that of the countrywide vegetation reported in the literature.
The significance of the regression coefficients also depends on the methodology used.The generalized least squares method produces greater p-values than ordinary linear regression if the serial correlation is positive, which is the case in this study [65].Hence, if the result of the generalized least squares regression is statistically significant, the corresponding result of ordinary linear regression, which has been used in the abovementioned literature, should also be significant.The recognized drought resistance of shrubs may have led to the differences between our findings for shrublands and the results for countrywide vegetation in the literature.
Rainfall shifts from the summer to the spring and fall have been hypothesized to favor shrubs in the competition with summer-active grasses [41,66].In this study, we found a regional case for the rainfall shift from the summer to fall in the subtropical evergreen shrublands, associated with a strong increase in EVI.
The significantly increasing trend of annual average EVI of the subtropical evergreen and temperate deciduous shrublands is mostly because of the significantly enhanced EVI in the summer and fall (Table 2).Furthermore, the significantly increased summer and fall EVI also coincided with the increased summer temperature and the rainfall shift from the summer to the fall for the former, and the increased rainfall in both the summer and the fall for the latter (Tables 3 and 6).
Unlike the increasing trend of two major shrublands and two desert types, the annual average EVI of two subalpine shrublands decreased over the 13-year period (only EVGNMT was significant, p < 0.1).The continuously reduced precipitation and the increased temperature in the four seasons (Table 3) might lead to a drier condition, and finally resulted in the EVI reduction of the two subalpine shrublands.Similarly, the decreasing trend (although insignificant) of the subalpine deciduous shrublands may have been caused by a combination of the spring and summer warming and a marginally increased precipitation during the growing season (Tables 3, 5 and 6).These shrublands are located in the steep slopes of the eastern Tibet Plateau, and the harsh environment is extremely sensitive to soil erosion and water retention.
The two desert shrublands are located in the arid region and is mainly limited by water.Hence, the EVI of the two desert shrublands was sensitive to precipitation, and the small increase in precipitation in the spring and summer led to a significant increase in EVI in the spring and summer (Tables 2, 3  and 6).Although a significantly increasing trend in the temperature in DSRTHC was found, the EVI showed no relationship with the temperature.

Growing Seasonal Length and Climate Explained the Changes in Annual Average EVI
The changes in the growing season length are responsible for the changes in the annual average EVI (Table 6) of the subtropical evergreen, temperate deciduous, and subalpine deciduous shrublands, conforming the current understanding [67].The annual average EVI increased significantly with the growing season length for the three shrubland types.
The growing season length is not the only explanation for the changes in the shrubland EVI (Table 6).Seasonal temperature in the spring and summer contribute to the additional increase in annual EVGNST EVI, as well as precipitation in the spring and fall.Increased fall temperature and summer precipitation positively affected the annual DCDSTP EVI, while the spring precipitation reduced it.The summer precipitation contributed significantly to the EVI of two desert shrublands largely because their active growth mostly occurs in the summer.We failed to extract the phenology of two desert shrublands due to the noisy EVI time series caused by sparse vegetation coverage.Nevertheless, the rainfall increase in the winter and spring may advance the starting day of DSRTHC, which was suggested by former study that in the water limited areas increasing precipitation in preceding months can advance the growing season [68].In addition, the increased spring temperature may advance the starting day of DSRTHC and eventually lead to the EVI increase in the spring (Tables 2 and 3) [62,69].Similarly, a weaker increase in the winter rainfall and spring temperature led to a small increase in the spring DSRTTP EVI.The continuous warming and precipitation reduction throughout the year led to a drier condition, and eventually resulted in a reduction in the annual average EVGNMT EVI (Tables 2, 5 and 6).However, Table 5 shows all positive coefficients for climate variables for EVGNMT, and the significantly increased temperature and non-significantly changed precipitation indicates that the growing season length of EVGNMT should have extended, but the growing season length of EVGNMT was decreased (Tables 3 and 4).Hence, the EVGNMT EVI variation was not only derived by regional climate dynamics but also other factors, such as topography (EVGNMT located in the southeast of Tibet-plateau), which is complicated in high mountainous regions and can affect regional climate through altering spatial patterns of local hydrothermal conditions and finally have effect on vegetation variation [70] and human activity is an important factor of vegetation variation which was ignored in our study.The Tibet plateau is a sensitive region; human activity may affect the vegetation variation through land use change, soil erosion and so on.In summary, increasing temperature and marginally decreased precipitation can partly explain the decrease of EVGNMT EVI, and the other factors such as land use change and deforestation need further research in more detail.

Conclusions
In this paper, we analyzed the spatiotemporal trend of MODIS (Moderate Resolution Imaging Spectroradiometer) EVI (Enhanced Vegetation Index) of six shrubland types in China from 2001 to 2013 and its relationship to intra-and inter-annual regional climate dynamics, which confirmed our hypotheses.The shrublands EVI of China increased significantly from 2001 to 2013 at a rate of 1.01 ˆ10 ´3 EVI¨a ´1.The two main shrubland types (subtropical evergreen and temperate deciduous shrublands) increased at a very high speed of 1.88 and 1.89 ˆ10 ´3 EVI¨a ´1, respectively.In addition, the EVI of two desert shrub types also increased significantly, whereas the subalpine evergreen shrubland EVI decreased significantly.We also detected a lengthened growing season of temperate deciduous shrubland, which partly followed the second hypothesis.The last hypothesis, that regional climate dynamics and growing season length significantly affected the annual averaged shrublands EVI, was also verified.The precipitation variation played a more important role in the EVI variation than temperature, which can be attributed to the water-limited habitat of shrublands.
The increased annual average EVI of the shrublands in China may have significantly enhanced their contribution to the regional ecosystem production and carbon sequestration.The strong increasing trend of EVI for the temperate deciduous shrubland in northern China may point on increasing capacity for their expansion into grasslands, so as the extended growing season to fall that may favor shrublands in competition with the summer-active grassland.

Figure 1 .
Figure 1.Pixel-based nonparametric trend analysis from 2001 to 2013.The circles indicate the distribution of different shrubland types.The first number in bracket indicates averaged EVI during the 13 years, and the second number indicates the percentage of total shrubland area.DCDSMT and EVGNMT are in one circle because there is no clear boundary between them.DCDSMT mainly distributed in the upper half of the circle and EVGNMT distributed in the bottom half.

Figure 1 .
Figure 1.Pixel-based nonparametric trend analysis from 2001 to 2013.The circles indicate the distribution of different shrubland types.The first number in bracket indicates averaged EVI during the 13 years, and the second number indicates the percentage of total shrubland area.DCDSMT and EVGNMT are in one circle because there is no clear boundary between them.DCDSMT mainly distributed in the upper half of the circle and EVGNMT distributed in the bottom half.

Figure 2 .
Figure 2. Annual mean EVI (a); temperature (b); and annual precipitation (c) of the China shrublands.Solid lines represent the regression with Generalized Least Squares method to derive trend from the time series.

Figure 2 .
Figure 2. Annual mean EVI (a); temperature (b); and annual precipitation (c) of the China shrublands.Solid lines represent the regression with Generalized Least Squares method to derive trend from the time series.

Figure 3 .
Figure 3. Regression of annual averaged EVI of the six shrubland types on pertinent climate variables: Regression predicted vs. observed values.Observed EVI indicates EVI obtained from MODIS and predicted EVI indicates values predicted by generalized least square based regression model.

Figure 3 .
Figure 3. Regression of annual averaged EVI of the six shrubland types on pertinent climate variables: Regression predicted vs. observed values.Observed EVI indicates EVI obtained from MODIS and predicted EVI indicates values predicted by generalized least square based regression model.

Table 1 .
Pixel-based nonparametric trend analysis.TS stands for Theil-Sen median slope estimator.

Table 2 .
Temporal trends in annual and seasonal EVI based on the Generalized Least Squares Regression.k denotes Linear coefficient for year.ΔR is the relative change computed as the slope k in the second column times 13 then divided by the annual average EVI in 2001.CONTRY denotes the pooled EVI of the countrywide shrublands; p values were based on one-tailed hypothesis, and "*", "**", and "***" are for p-values less than 0.1, 0.05, and 0.01, respectively.

Table 1 .
Pixel -based nonparametric trend analysis.TS stands for Theil-Sen median slope estimator.

Table 2 .
Temporal trends in annual and seasonal EVI based on the Generalized Least Squares Regression.k denotes Linear coefficient for year.∆R is the relative change computed as the slope k in the second column times 13 then divided by the annual average EVI in 2001.CONTRY denotes the pooled EVI of the countrywide shrublands; p values were based on one-tailed hypothesis, and "*", "**", and "***" are for p-values less than 0.1, 0.05, and 0.01, respectively.

Table 3 .
Temporal trends in annual and seasonal climate based on the Generalized Least Squares Regression.p values were based on one-tailed hypothesis, and "*", "**", and "***" are for p-values less than 0.1, 0.05, and 0.01, respectively.

Table 4 .
Temporal trends in the length of growing season in 2001-2013.GSL stands for growing season length and ∆GSL denotes the changes in the growing season length between 2001 and 2013.p values were based on one-tailed hypothesis, and "***" is for p-values less than 0.01.

Table 5 .
Regressed coefficients of the growing season length on seasonal precipitation (P, day¨mm ´1) and mean temperature (T, day¨˝C ´1).SPR, SMR, FAL, and WNT indicate spring, summer, fall, and winter, respectively.p values were based on one-tailed hypothesis, and "**" and "***" are for p-values less than 0.05 and 0.01, respectively.

Table 5 .
Regressed coefficients of the growing season length on seasonal precipitation (P, day•mm −1 ) and mean temperature (T, day•°C −1 ).SPR, SMR, FAL, and WNT indicate spring, summer, fall, and winter, respectively.P values were based on one-tailed hypothesis, and "**" and "***" are for p-values less than 0.05 and 0.01, respectively.