Inconsistency of Global Vegetation Dynamics Driven by Climate Change: Evidences from Spatial Regression

: Global greening over the past 30 years since 1980s has been conﬁrmed by numerous studies. However, a single-dimensional indicator and non-spatial modelling approaches might exacerbate uncertainties in our understanding of global change. Thus, comprehensive monitoring for vegetation’s various properties and spatially explicit models are required. In this study, we used the newest enhanced vegetation index (EVI) products of Moderate Resolution Imaging Spectroradiometer (MODIS) Collection 6 to detect the inconsistency trend of annual peak and average global vegetation growth using the Mann–Kendall test method. We explored the climatic factors that affect vegetation growth change from 2001 to 2018 using the spatial lag model (SLM), spatial error model (SEM) and geographically weighted regression model (GWR). The results showed that EVI max and EVI mean in global vegetated areas consistently showed linear increasing trends during 2001–2018, with the global averaged trend of 0.0022 yr − 1 ( p < 0.05) and 0.0030 yr − 1 ( p < 0.05). Greening mainly occurred in the croplands and forests of China, India, North America and Europe, while browning was almost in the grasslands of Brazil and Africa (18.16% vs. 3.08% and 40.73% vs. 2.45%). In addition, 32.47% of the global vegetated area experienced inconsistent trends in EVI max and EVI mean . Overall, precipitation and mean temperature had positive impacts on vegetation variation, while potential evapotranspiration and vapour pressure had negative impacts. The GWR revealed that the responses of EVI to climate change were inconsistent in an arid or humid area, in cropland or grassland. Climate change could affect vegetation characteristics by changing plant phenology, consequently rendering the inconsistency between peak and mean greening. In addition, anthropogenic activities, including land cover change and land use management, also could lead to the differences between annual peak and mean vegetation variations.


Introduction
Terrestrial vegetation controls the cycle of carbon, water and energy between land soil and atmosphere through biochemical processes such as photosynthesis and evapotranspiration and is strongly influenced by hydrothermal conditions and climate change and can affect the climate system in return [1][2][3]. In addition to directly providing ecosystem services such as food, raw materials, and landscape aesthetics, vegetation also has the potential of climate regulation, carbon sequestration and air purification, which are of great importance to maintaining the stability of a terrestrial ecosystem under global change [4,5]. Thus, systematic monitoring, detecting and quantifying vegetation dynamics and their response and feedback to global change have elicited a wide range of attention in multiple subjects.
It is currently unlikely to obtain vegetation properties and variations at the global scale using ground-based observations [5]. Vegetation indices (VIs) products, such as normalized difference vegetation index (NDVI) and enhanced vegetation index (EVI) derived from satellite observations, provide a possibility for long-term vegetation monitoring on a large scale [6,7]. VIs products from the Advanced Very-High-Resolution Radiometer (AVHRR) sensor, Système Pour l'Observation de la Terre VEGETATION (SPOT-VGT) and Moderate Resolution Imaging Spectroradiometer (MODIS) sensors on board different satellites have been widely used in vegetation growth monitoring and crop mapping [8][9][10]. Although providing the longest record of NDVI data (1981-present), AVHRR has the problems of orbit drifts and data inconsistency between various sensors. SPOT-VGT data discontinuity was also found in some areas due to sensor differences in spectral response functions [11]. MODIS data is considered to have a higher temporal consistency than AVHRR and SPOT-VGT data due to no orbit drift and sensor shifts problems and are usually used as reference data [12][13][14]. However, the negative influences of sensor degradation have been captured in MODIS-Terra Collection 5 VIs products [15]. Consequently, MODIS Collection 6 VIs products with the improved algorithm were released to address the sensor degradation problem, and the effects have been supported by some studies [16,17]. Compared to NDVI, EVI minimizes canopy-soil and atmospheric influences and improves sensitivity over dense vegetation conditions, and its reliability has been recognized by comparison analysis of multiple VIs products [10].
The overall greening trend of global vegetation has been widely reported and discussed since the early 1980s, supported by comparisons of multiple satellite observations [18], forest inventories [19] and process-based model simulations [20,21]. Pan et al. [22] explored the increasing global browning trend hidden in overall greening during 1982-2013 by using the ensemble empirical mode decomposition method and piecewise linear regression models; Gao et al. [23] verified the significant trend of global cultivated land greening during 1982-2015 by using two long-term satellite LAI datasets; Zhang et al. [24] found that the previous browning trend monitored by MODIS Terra-C5 needed to be reconsidered due to sensor degradation, while the trend from MODIS Terra-C6 was confirmed by AVHRR and enhanced land carbon sink data [25]. Although a lot of studies have re-examined global greening trend, whether there is a consistent trend between annual peak and average growth of vegetation remains unknown. The annual peak vegetation growth closely related to environmental change is critical in characterizing the capacity of terrestrial ecosystem productivity [20]. Vegetation growth is an ever-changing dynamic process. The annual average vegetation growth, which was selected in the most previous studies, can only reflect the overall state of a certain year but obscure details of vegetation growth. To comprehensively understand the changing ecosystems, it is necessary to track vegetation growth change from multi-dimensions.
Climate change characterized by global warming has been a dominant driver of greening over 28% of the global vegetation regions [21]. However, the global vegetation-climate relationship is complex and has firm spatial heterogeneity [26]. Warming has noticeable positive effects on vegetation growth in the temperate and arctic regions [27], while rising temperatures could limit vegetation growth in tropical regions [28] where the ambient temperature is the proximity of the physiological optimum [29]. In arid and semi-arid regions, vegetation is more sensitive to precipitation change than temperature due to water limitation. In contrast, in humid regions, the impacts of precipitation on vegetation variation are weaker than temperature [30]. In sum, the climate-vegetation relationship geographically varies with the climatic environment [5,23]. In addition, although temperature and precipitation have been the hot spots for a long time, the impacts of climate change on vegetation are more complex than that. For example, the impact of other climatic factors, such as vapour pressure deficit on vegetation growth, have been highlighted, which could offset the effect of rising CO 2 concentration, suggested to be considered in evaluating vegetation responses to climate change [31].
At present, most studies used the non-spatial regression methods to depict the relationship between vegetation variation and climate change, such as multiple linear regression [32] and partial correlation analysis [33]. However, the theoretical assumption of these methods is that the spatial data is statistically independent and uniformly distributed [34,35], which ignores influences of spatial autocorrelation. According to the first law of geography, spatial autocorrelation is ubiquitous, and the correlation between adjacent locations is usually stronger than that between distant locations [36,37]. Influenced by the spatial interaction and spatial diffusion of endogenous or exogenous factors, geographic data may no longer be independent of each other, but related to each other [38]. For example, the flow characteristics of the atmosphere and water supporting vegetation growth can suggest that vegetation growth are absolutely not spatially independent [37,39,40]. Previous empirical research within ecology have revealed that regression coefficients may radically shift between non-spatial and spatial (taking autocorrelation into account) regression modelling, resulting in erroneous conclusions [39,[41][42][43][44]. Given that spatial autocorrelation exists extensively in spatial data, it is important to use spatially explicit models to explore the climate-driving mechanism of vegetation growth. Therefore, the aims of this study are: (1) to examine inconsistency trends of global vegetation in peak and average growth during 2001-2018; (2) to reveal the overall relationship between climate change and vegetation growth using spatial regression models (at the global level); (3) to find out the spatial heterogeneity characteristic of climate driving using the geographically weighted regression model (at the local level).

Global Datasets and Pre-Processing
EVI, land cover and climatic datasets used in this study are displayed in Table 1. The raw MODIS data were mosaicked, re-projected and converted to Geo TIFF from HDF format using the MODIS Reprojection Tool. We averaged monthly EVI as the annual EVI mean and extracted the maximum monthly EVI as the annual EVI max . In the spatial regression analysis, we resampled EVI data into 0.5 • × 0.5 • spatial resolution by using the bilinear algorithm to match the climatic gridded data as well as the basic analytical unit. Global terrestrial land cover data are from the International Geosphere-Biosphere Programme (IGBP) classification layer of MCD12Q1 ( Figure 1). To eliminate the direct effects of transformation between vegetated and non-vegetated land during 2001-2018, we masked out non-vegetated areas including permanent wetlands, urban and built-up lands, permanent snow and ice, barren and water bodies by overlaying the IGBP land cover data during the period of 2001-2018. Global annual climatic gridded data were generated by extracting the maximum, minimum, mean or sum value from the monthly climatic gridded data. For example, annual precipitation data were generated by summating monthly precipitation, while annual maximum temperature data were obtained by calculating the maximum value of the monthly maximum temperature. All data pre-processing was accomplished in R x64 4.0.2 and RStudio (http://www.r-project.org/, accessed on 18 September 2020).

Methods
A methodological flowchart of this study is shown in Figure 2. Firstly, to find out the greening and browning areas and the area with inconsistent changes, we detected the trend of global vegetation growth variation during 2001-2018 from two dimensions of annual peak and mean using the Mann-Kendall test method recommended by the World Meteorological Organization in Section 2.2.1. Secondly, to prove our theoretical hypothesis that spatial autocorrelation exists in vegetation growth change, the global univariate Moran's I was first used to diagnose the spatial autocorrelation in the response variable in Section 2.  Thirdly, to incorporate spatial autocorrelation in the response of vegetation growth to climate change, we used multiple spatial regression models to analyse the spatial relationship between climate change and vegetation variation at the global level. The first step is to perform the ordinary least squares (OLS) model, which is under the assumption that there is no autocorrelation [45] because the spatial lag model (SLM) and the spatial error model (SEM) are developed from the OLS by incorporating spatial autocorrelation into the regression by means of a spatial weight matrix [34]. SLM was used when spatial autocorrelation exists in dependent variables. SEM was used when spatial autocorrelation exists in the residual. The second step is to select SLM or SEM based on the statistical significance of the Lagrange multiplier (LM)-lag and LM-error, or robust LM-lag and robust LM-error. The third step is to determine the better performance model, with a larger maximum likelihood logarithm (LIK) and a smaller Akaike information criterion (AIC) and Schwarz criterion (SC).
Finally, the relationship between variables varies with the change of geographical location due to the differences in the natural environment and human disturbance in different regions. This changing relationship also needs to be considered in the spatial analysis [46]. To reveal the differentiated local characteristics hidden in the overall correlation, we used the geographically weighted regression (GWR) model to measure the spatial heterogeneity of the relationship between vegetation growth and climate change by obtaining local regression results at each spatial unit in Section 2.2.4. The spatial autocorrelation analysis and the OLS, SLM and SEM were conducted in GeoDa software (http://geodacenter.github.io/, accessed on 8 February 2021). The GWR model was conducted in GWR4 software (https://gwr4.software.informer.com/, accessed on 11 March 2021).

Trend Detection of Vegetation Growth
The linear trend of global vegetation growth during 2001-2018 at each grid cell were estimated using Sen's slope method, also known as the Theil-Sen median method [47]. The method is a robust non-parametric statistical trend calculation method that has high computational efficiency and is insensitive to measurement error and outlier data [48]. Sen' slope was calculated by Equation (1): where the median refers to the mean value of all the slopes, and x i and x j represent the EVI values of years i and j.
Then the Mann-Kendall test method was used to test the significance of the Sen' slope [49,50]. The significant confidence level with p < 0.05 corresponds to the absolute value of the Z statistic >1.96. The Sen' slope and Mann-Kendall test for EVI data at each grid cell were accomplished in R x64 4.0.2.

Spatial Autocorrelation Analysis
Moran's I statistic is arguably the most commonly used indicator of global spatial autocorrelation. It was initially suggested by Moran [51], and popularized through the classic work on spatial autocorrelation by Cliff and Ord [35]. For an observation at location i, this is expressed as z i = x i − x, where x is the mean of variable x. Moran's I statistic is then: where n is the number of observations, w ij is the elements of the spatial weights matrix, x i and x j are the observed value of the location i and its surrounding location j, x is the mean of variable x. At a given significance level, when Moran's I > 0, it indicates a positive correlation between the observed values, and similar attributes cluster together. That is, the high value is adjacent to the high value, and the low value is adjacent to the low value; when I < 0, it indicates a negative correlation between the observed values, and the observations are dispersed; when I = 0, the observed value is randomly distributed.

Spatial Regression at the Global Level
Anselin [34] put forward the general form of spatial regression. When ρ = 0, β = 0, α = 0, the model is the ordinary least squares (OLS) model; When ρ = 0, β = 0, α = 0, the model is a spatial lag model (SLM), that is, the dependent variable of a location is not only related to the independent variable of the location, but also related to the dependent variable of the neighbourhood. When ρ = 0, β = 0, α = 0, the model is a spatial error model (SEM), that is, the dependent variable of a location is not only related to the independent variable of the location, but also related to other variables not considered at adjacent regions.
where y is the dependent variable, x is the independent variable, W 1 is the spatial weight matrix of the dependent variable, ρ is the coefficient of the spatial lag variable W 1 y, β is the coefficient of x, ε is the residual, W 2 is the spatial weight matrix of ε, α is the coefficient of ε, µ is the random error of normal distribution, σ is the variance of µ.

Spatial Regression at the Local Level
The GWR model is essentially the locally weighted least squares regression model, which is an extension of the OLS model [52]. In GWR, the weight of the observations is no longer constant during the regression process but is weighted by the degree of adjacency to location i [53]. The model structure is as follows:

Temporal Trend of Global Vegetation Growth
EVI max can reflect the potential productivity of terrestrial vegetation, and EVI mean represents the average vegetation growth state in a year, which depicts two aspects of vegetation change. Globally, peak (EVI max ) and average (EVI mean ) growth in global vegetation consistently showed linear increasing trends during 2001-2018, with the global averaged trend of 0.0030 yr −1 (p < 0.05) and 0.0022 yr −1 (p < 0.05). In terms of EVI max , 18.16% of the global vegetated areas showed a statistically significant (Mann-Kendall test, p < 0.05) greening during 2001-2018, and 3.08% of the global vegetated areas showed a statistically significant (Mann-Kendall test, p < 0.05) browning during 2001-2018 ( Figure 3a). By overlaying land cover, we found that the most dramatic greening occurred mainly in those areas with cropland agricultural activities, such as Northern China, India, Central-North America and Southeast Europe. The fastest degradation areas were mainly the grassland of Africa and South America, such as Tanzania, Nigeria and Brazil.
In terms of the EVI mean , 40.73% of the global vegetated areas showed a statistically significant (Mann-Kendall test, p < 0.05) greening during 2001-2018, and 2.45% of the global vegetated areas showed a statistically significant (Mann-Kendall test, p < 0.05) degradation trend (Figure 3b). The most obvious greening areas mainly occurred in China, India, Canada and Europe, covering a variety of land types including cropland, shrubland, forests and savannas. Vegetation degradation areas were mainly grassland and savannas in areas of South America and southern Africa.

Inconsistent Global Vegetation Growth in Terms of EVI max and EVI mean
It is worth noting that EVI max and EVI mean experienced consistent changes in 15.97% of the global vegetated areas (14.99% for greening and 0.98% for browning) (Figure 4). However, 32.47% of the global vegetated areas experienced inconsistent changes for EVI max and EVI mean . Specifically, 25.74% of the global vegetated areas that experienced significant greening in EVI mean , with no increase in EVI max simultaneously, occurred mainly in of Europe, Russia, Central Africa, North America and China. On the contrary, 3.17% of the global vegetated areas that experienced significant browning in EVI max , with no increase in EVI mean simultaneously, were mainly distributed in Northern Canada, Eastern Russia, Southern Australia, and were scattered in South America and Africa. There was also 2.10% of the global vegetated area showing a significant browning in EVI max, with no decrease in EVI mean simultaneously, scattered in Africa, South America and Canada, especially in Argentina and Madagascar. While 1.46% of the global vegetated area showed a decreased trend for EVI mean , with no decrease in EVI max simultaneously, such as Central Africa, Central Russia, Eastern Canada and Eastern Brazil. In addition, we found that EVI max and EVI mean had opposite trends in some areas; however, these areas together accounted for only 0.56% of the inconsistent area. EVI mean increased but EVI max decreased in 0.52% of the inconsistent area, while EVI max increased but EVI mean decreased in 0.04% of the inconsistent area.

Relationship between Climate Change and Vegetation Growth
In this study, the Moran's I values for the EVI max and EVI mean changes during 2001-2018 were 0.273 and 0.549 (p < 0.05), showing that a significant positive spatial autocorrelation existed in EVI changes in the past 18 years. The multicollinearity condition number was 4.74 (<30) from the OLS model estimation results, indicating no multicollinearity problems in the explanatory variables. The test statistics of LM-lag, LM-error, robust LM-lag and robust LM-error that form the OLS model were all significant (p < 0.0001); thus, SLM and SEM were both built to analyse spatial global correlation between vegetation and climate factors. By comparing the index of LIK, AIC and SC, SLM has the largest LIK value and the smallest AIC and SC value, showing that SLM is superior to SEM and OLS in this study (Table 2). R-squared represents the degree to which nine climatic factors explained the global vegetation change in the past 18 years. R 2 in SLM were 0.5370 and 0.2499 for EVI mean and EVI max , respectively, indicating climatic drivers and the spatial lags of EVI could explain more than half and the 24.99% changes of EVI mean and EVI max , respectively.   The regression results of the best fit model (SLM) in this study for EVI max and EVI mean were shown in Table 3. The spatial lag of EVI max in SLM passed the statistical significance test (p < 0.05) during 2001-2018, proving that changes in EVI were correlated not only to these climatic factors but also to EVI variation in its adjacent areas. The coefficient indicates that the vegetation in a certain location might change by 0.7472 units for every 1-unit change of vegetation in its adjacent areas. Precipitation and mean temperature had a statistically significant (p < 0.05) positive correlation with EVI max during 2001-2018, while potential evapotranspiration and vapour pressure had a statistically significant negative correlation with EVI max . Compared to EVI max , EVI mean had a higher R-squared than EVI max , 0.5370 vs 0.2499 (Table 2), indicating a higher climatic driving explanation of EVI mean . Except for precipitation and mean temperature, minimum temperature also had a statistically significant (p < 0.05) positive correlation with EVI mean change. The coefficient of spatial lag for EVI mean was 0.8721, suggesting EVI mean had a higher correlation with its adjacent EVI mean and, therefore, a stronger spatial dependence than EVI max .

Spatial Heterogeneity of the Climatic Driving
After examining the relationship at the global level by SLM, we then focused on the representative areas with significant greening or browning identified in the previous trend analysis. We further explored how climate change affects vegetation in these areas with the GWR model results. Five highly representative regions were selected for further analysis, located in China, India, North America, Brazil and Africa. The areas in China, India and North America were characterized by vegetation greening, while the areas in Africa and Brazil were characterized by vegetation browning. Local coefficients of climatic drivers for EVI max and EVI mean were mapped in Figures 5 and 6, respectively.
As shown in Figure 5, in Northern China, EVI max was positively affected by precipitation, potential evapotranspiration, minimum temperature and humid days, and was limited by maximum temperature and vapour pressure. However, in Southern China, maximum temperature and vapour pressure had positive effects, as well as potential evapotranspiration and humid days. The reason was that there are semi-arid areas in the north where precipitation is low and potential evapotranspiration is far greater than precipitation and the high temperature would lead to insufficient rainwater irrigation for crops and grass. It indicated that, due to the difference in climatic conditions (an arid or humid area) and vegetation types (cropland or grassland), the responses of EVI max to climate change in South and North China are inconsistent. In India, except minimum temperature, all other factors had positive impacts on the EVI max . In North America, EVI max was positively influenced by potential evapotranspiration and humid day but negative by precipitation and vapour pressure. In Brazil, the degradation of the EVI max might be due to the decrease in maximum temperature, potential evapotranspiration and humid days and the increase in mean temperature and diurnal temperature range. In Africa, the degradation of EVI max might be due to the joint decrease in mean temperature, humid days and diurnal temperature range and the increase in evapotranspiration and maximum temperature. As shown in Figure 6, we found that maximum temperature, minimum temperature and humid days also had opposite effects on EVI mean in Northern and Southern China. In addition, precipitation and vapour pressure had consistently positive effects on both the south and north. In India, potential evapotranspiration, maximum temperature and humid days had positive effects on EVI mean . In North America, temperature mainly had a positive effect in these greening areas, but potential evapotranspiration and diurnal temperature range had negative effects. However, precipitation and vapour pressure had opposite impacts on Canada and the United States, with negative and positive, respectively. In Brazil, all browning areas of EVI mean occurred in grassland, which might be affected by the decrease in precipitation, mean temperature, potential evapotranspiration and wet days. In Africa, degradation of vegetation also occurred in grassland and was caused by the decrease in precipitation, mean temperature, vapour pressure and wet days.
We defined the climate factor having the largest absolute value of the local coefficient with a statistical significance (p < 0.05) as the dominant climatic factor influencing vegetation change (Figure 7a,b). In the meantime, the local coefficients of the dominant climatic factors of EVI max and EVI mean were mapped in Figure 7c,d. We found that the EVI max change was strongly influenced by precipitation in 14.02% of the global vegetated areas during 2001-2018, followed by vapour pressure (12.56%), minimum temperature (9.87%), mean temperature (9.10%) and potential evapotranspiration (6.96%) (Figure 7a). 21.06% of the global vegetated areas where peak vegetation growth had no significant correlation (p > 0.05) with any climatic factor. In terms of the EVI mean , the global vegetated areas were strongly affected by mean temperature (17.36%) and precipitation (16.97%) (Figure 7b). There were 11.86% of the global vegetated areas where EVI mean change had no significant correlation (p > 0.05) with any climatic factor.

Comparison of Global Vegetation Trend Results and Uncertainties
The trends of EVI mean and EVI max during 2001-2018 were detected in this study and were compared with previous studies (Table 4). Zhang et al. [24] found that the global vegetation has an increasing trend of 0.0028 yr −1 and 0.0023 yr −1 at a significant level (p < 0.0001) in EVI mean and EVI max during 2001-2015, respectively. The results were consistent with our results of 0.0022 yr −1 and 0.0030 yr −1 in EVI mean and EVI max , respectively. Moreover, the results from NDVI also showed the greening trend, with 0.0022 yr −1 and 0.0015 yr −1 in NDVI mean and NDVI max during 2001-2015, respectively [24], 0.0013 yr −1 in NDVI max during 1982-2011 [20], and 0.0012 yr −1 in the growing season NDVI during 1982-2013 [54]. Overall, the global greening trends were found in both our and previous studies.
However, there were some differences in global vegetation greening or browning areas due to different satellite products, vegetation greenness indicators and time range. For example, Ding [57] found that 18.9% of the global vegetated area had a greening trend for EVI mean during 2000-2015, while~3% of browning (vs. 40.73% and 2.45% in this study); Chen [56] found that 34.1% of the global vegetated area showed a greening and 4.85% of browning for LAI from MODIS during 2000-2017, while 22.42% of greening and 13.54% of browning for LAI from AVHRR during 2000-2016. Furthermore, although the same satellite-derived data (MODIS Terra-C6 EVI) were used, there were still differences in area ratios of greening and browning between our study and the relevant studies, which might be due to the following reasons: (1) different vegetation monitoring time range; (2) different spatial resolution of EVI data; (3) different land cover data and starting reference year; (4) different trend analysis methods. Notes: NDVI max and EVI max refer to the annual maximum NDVI and EVI, respectively, and unsubscripted ones represent the annual mean values; NDVI gs and LAI gs refer to NDVI and LAI for growing season, respectively; **** p < 0.0001, *** p < 0.01, ** p < 0.05, * p < 0.1.

Potential Causes Inducing Inconsistencies in Vegetation Change
Our study found that only 15.97% of the global vegetated area experienced a consistent change (enhancement or degradation) in peak and mean growth of vegetation, and 32.47% of the global vegetated area experienced inconsistent changes in peak and mean vegetation growth. The first potential cause could be changes in plant phenology or growing season induced by climate change in these areas. In some areas, the annual EVI mean showed an increasing trend while the annual EVI max remained unchanged and even decreased, indicating that with time going on, the monthly EVI, which are larger than the original annual EVI mean , increased during 2001-2018, and in this case the growth circle of vegetation might be prolonged. In contrast, for the areas with increased EVI max but no significant change or decrease in EVI mean , the growth circle of vegetation might be shortened. This speculation can be supported by the consistent evidence for plant phenological changes from in situ and satellite observations [59]. The earlier beginning of the growing season and autumn postponement has been widely reported in Europe, North America and China since the 1980s [60][61][62][63]. Moreover, with the stagnation of warming, spring green-up advancement's trend might have slowed down or even reversed since the 2000s [64,65]. However, this speculation still is of great uncertainty due to the unclear vegetation phenological information and its variations in recent 20 years, which needs to be investigated in detail in further studies.
Plant phenology changes are determined by various biological and environmental factors such as nutrient and water availability, temperature and photoperiod. Temperature is generally regarded as one of the most critical controls of plant phenology through multiple processes, such as inducing the plant endodormancy by cold temperature [66] and breaking the ecodormancy by the accumulated warming [59]. For vegetation in pasture regions, precipitation variation had a significant limiting effect on its productivity [67]; this might result in the inconsistency of peak and mean greenness due to variations of annual maximum and minimum rainfall. What is more, the impact of climatic conditions on vegetation growth is more complex due to the interactions with other environmental and climatic factors [68].
In addition, anthropogenic activities, including land cover change and land use management, also could lead to the differences between annual peak and mean vegetation variations. On the one hand, the conversion between different vegetation types might directly result in the observed inconsistency between the peak and mean greenness, due to the different responses to environmental changes determined by different vegetation properties, such as thermal adaptability and photosynthetic efficiency [69]. For example, crops have higher photosynthetic efficiency than other non-crops [20]. On the other hand, land-use intensity and management could also explain the inconsistency largely. For example, anthropic seasonal irrigation and fertilization could improve the peak greenness and productivity of croplands [20,56,70], but the annual mean greenness might not enhance simultaneously.

Spatial Heterogeneity of Vegetation Growth Driven by Climate Change
Consistent with previous findings, we found that rainfall and temperature significantly impact vegetation both in peak and mean growth [20,21]. However, except for precipitation and temperature, vegetation growth was also found to be significantly correlated with other climatic factors such as potential evapotranspiration and vapour pressure in our study (Table 2). Furthermore, previous studies had neglected the spatial heterogeneity of this response. Using the GWR model, we revealed the spatial pattern of each climatic factor influencing vegetation variation. We found that not all regions had the strongest correlation between vegetation change and precipitation and temperature, which were generally considered to be the main climatic factors. For example, potential evapotranspiration contributed dominantly to the trend of EVI max and EVI mean over 6.96% and 7.54% of the global vegetated area, respectively. Significant positive effects of potential evapotranspiration occurred mainly in Brazil, Northern Europe and Northeast China, while adverse effects were mainly found in Africa (Figure 7). Vapour pressure had dominant contribution to the trend of EVI max and EVI mean in 12.56% and 8.95% of the global vegetation. The greening induced by climate change in the Tibetan Plateau was mainly attributed to increasing vapour pressure and temperature in this study (Figure 7), whereas that was rising in a previous study [21].

Limitations and Further Directions
Data continuity and accuracy of the VIs products are critical to accurately detect subtle changes in vegetation, which is important to the assessment of vegetation dynamics [12]. VIs are a kind of spectral vegetation index derived from satellite remote sensing, which is calculated by spectral band reflectance [14]. Thus, unlike LAI, VIs cannot be directly validated and calibrated continuously in time series by in situ measurements. A number of comparisons have been made between different satellite-based VIs datasets [12,14,16,71,72] and MODIS VIs data with improved technology, and no shifts of sensor is considered to be superior to other products in terms of data temporal consistency and has been widely used for reference data [10]. The latest MODIS Collection 6 (C6) VIs data providing several algorithmic improvements and calibration adjustments were released for correcting the negative influence of sensor degradation found in MODIS Collection 5 (C5) VIs data [17]. After release of MODIS C6 products, some studies began to evaluate differences between C5 and C6 VIs both on a local and global scale to verify that C6 products had eliminated effects of sensor degradation, and highlighted the need of re-analysing some previous results based on MODIS C5 VIs products [17,24,73]. Therefore, we selected MODIS C6 data to examine global vegetation growth trend in this study. However, the differences between MODIS EVI and NDVI were not sufficiently considered to reduce the uncertainty in detecting trend analysis in this study. Although EVI can improve reflectance sensitivity in dense vegetation areas, NDVI and EVI are generally regarded as two complementary datasets for providing more effective assessment of global vegetation dynamics [24]. Thus, it is necessary to evaluate the differences between the two datasets for monitoring vegetation dynamics using both EVI and NDVI data in the future. Moreover, research on the climatic driving mechanism of vegetation dynamic should be combined with VIs and other vegetation indicators, such as LAI and net primary productivity, which can be simulated in the process-based models, because either non-spatial or spatial regression analysis cannot explain the climate-driving mechanism from the ecological processes of vegetation growth but can only provide a hint of correlation.

Conclusions
This study detected the trend of global peak and average vegetation growth during 2001-2018 and mapped the inconsistency in vegetation growth, and the climatic factors that affected the inconsistency of vegetation growth were explored. The results showed that in terms of EVI max , 18.16% of the global vegetated areas are greening and 3.08% are browning, and in terms of EVI mean , there are 40.73% and 2.45%, respectively. The most dramatic greening of EVI max occurred mainly in those areas with cropland agricultural activities, and the fastest degradation areas of EVI max were mainly grassland and savannas of Africa and South America. Through mapping the consistency of global vegetation growth, it was found that from 2001 to 2018, 32.47% of the global vegetated area experienced inconsistent trends in EVI max and EVI mean .
The SLM was proved to be more suitable than SEM and OLS in this study for spatial regression at the global level. We found that precipitation and mean temperature had a statistically significant (p < 0.05) positive correlation with EVI max and EVI mean during 2001-2018, while potential evapotranspiration and vapour pressure had a statistically significant negative correlation. The results of SLM indicated that there was spatial autocorrelation in both EVI max and EVI mean change, which means the changes in EVI were correlated not only to these climatic factors but also to EVI variation in its adjacent areas.
The SLM results only indicated a correlation between EVI changes and climatic drivers on the whole but failed to reveal the spatial heterogeneity of climatic drivers.
The GWR model was used to explore the spatial heterogeneity of climatic drivers influencing vegetation change by obtaining regression results for each spatial unit. The results showed that the EVI max change was strongly influenced by precipitation in 14.02% of the global vegetated areas during 2001-2018, followed by vapour pressure (12.56%), minimum temperature (9.87%), mean temperature (9.10%) and potential evapotranspiration (6.96%). In terms of the EVI mean , the global vegetated areas were strongly affected by mean temperature (17.36%) and precipitation (16.97%). In China, maximum temperature and vapour pressure had opposite effects on EVI max in the north and south, and maximum temperature and humid days also had opposite effects on EVI mean .Due to the difference in climatic conditions (arid or humid area) and vegetation types (cropland or grassland), the responses of EVI to climate change were inconsistent. Climate change could affect vegetation characteristics by changing plant phenology, consequently rendering the inconsistency between peak and mean greening. In addition, anthropogenic activities, including land cover change and land use management, also could lead to the differences between annual peak and mean vegetation variations.

Data Availability Statement:
The data that support the findings of this study are available from the corresponding author upon reasonable request.