Climatic and Anthropogenic Contributions to Vegetation Changes in Guangdong Province of South China

Abstract

How to distinguish the relative role of climate change and human activities in vegetation dynamics has attracted increasing attention. However, most of the current studies concentrate on arid and semiarid regions, while the relative contributions of climate change and human activities to vegetation changes remain unclear in warm-humid regions. Based on the normalized difference vegetation index (NDVI) and climatic variables (temperature, precipitation, radiation) during 2001–2020, this study used the Theil–Sen median trend analysis, partial correlation analysis, and residual trend analysis to analyze the spatiotemporal pattern of vegetation trends, the response of vegetation to climate variations, and the climatic and anthropogenic contributions to vegetation dynamics in the warm and humid Guangdong Province of south China. Results showed that the NDVI in most areas exhibited an increasing trend. Changes in climatic variables displayed different spatial variations which, however, were not significant in most areas. Vegetation responded diversely to climate change with temperature as the most important climatic factor for vegetation improvement in most areas, while precipitation was the dominant climatic factor in the southern edge region and radiation was the dominant climatic factor in the central and western regions. Vegetation in most areas was influenced by both climate change and human activities, but the contribution rate of human activities was commonly much higher than climate change. The findings of this study are expected to enhance our understanding of the relative climatic and anthropogenic contributions to vegetation changes in warm-humid regions and provide a scientific basis for future ecological policies and ecosystem management in highly urbanized regions.

1. Introduction

Terrestrial ecosystems are significant global carbon stocks and carbon sinks, and vegetation variations are highly related to carbon budgets [1,2,3,4]. Vegetation change is a complex process affected by various natural and anthropogenic factors [5,6,7]. Hence, disentangling the driving forces of vegetation dynamics is essential for understanding carbon balance and developing adaptive strategies.

As one of the few countries where vegetation has significantly improved since 2000, China leads the global greening [8]. Specifically, forests in the southern part and croplands in the central part of China display the most obvious greening trend. Contrary to naturally vegetated areas, the urbanization process usually causes vegetation degradation in many regions [9,10,11]. However, recent studies indicate that urban vegetation growth enhancement is also widely observed in Chinese and worldwide cities [12,13]. Thus, vegetation changes in urban areas are much more complicated than in undisturbed areas. Previous work found that vegetation improvement happened in urban cores, while vegetation degradation occurred in newly developed peripheries in the Guangdong–Hong Kong–Macao Greater Bay Area urban agglomeration [14]. Vegetation dynamics are attributed to numerous driving forces, which are generally concluded as two categories: climate change and human activity [15,16].

Climate variation is extensively believed to have widespread impacts on vegetation change. Across the globe, water availability, temperature, and radiation constrained vegetation growth for 40%, 33%, and 27% of the vegetated areas on Earth, respectively [17]. The increasingly intense dry seasons have caused global vegetation productivity to decrease, and it is projected to potentially reduce by 10–30% in 2100 [18]. Effects of climate change on vegetation vary in different regions and ecosystems [19,20]. For example, worldwide plains are believed to be more vulnerable to climate variability than alpine regions [21]. Temperate broadleaf forest and temperate grassland have a higher sensitivity to extremely low precipitation than other biomes [22]. Crops are more sensitive to water and meadows are more sensitive to temperature in the Heihe River Basin of northern China [23]. Precipitation is the limiting factor of vegetation changes in the northeast and southwest of the Tibetan Plateau, while temperature and the combination of radiation and temperature are the prominent climatic drivers for the mid-east and the south, respectively [24]. Water-limited regions of vegetation productivity are primarily situated in the southern and eastern basins and piedmont plains of the Taihang Mountains of north China, while the northern and central regions with high elevations mainly suffer from low temperatures [25].

The effects of anthropogenic activities on vegetation change have attracted increasing concern in recent decades, particularly for rapidly developing regions [8,9,26,27]. On the one hand, industrialization and urbanization usually cause vegetation degradation and vegetation loss. For instance, forest vegetation decline has lasted for 30 years since the construction of industrial complexes in southern Korea [28]. Urban expansion widely resulted in vegetation degradation around urbanized regions, such as the Beijing–Tianjin–Hebei Region [27], the Yangtze River Delta [5], and the Pearl River Delta and Shantou city of Guangdong Province [9]. On the other hand, ecological projects and land-use management greatly contribute to vegetation greening. During the late 1990s and early 2000s, a series of national ecological restoration programs including The Grain for Green Program, The Natural Forest Protection Program, and The Shelterbelt Development Programs were carried out in China, and ecological restoration measures such as forest planation and mountain closure played a vital role in vegetation improvement [29]. Agricultural intensification and urban green space management such as fertilization and irrigation also dramatically enhanced vegetation greening [8,26].

Since both climate change and human activities have important impacts on vegetation dynamics, how to disentangle their relative contributions has become the focus of concern, especially in recent years. A large number of studies have made remarkable progress on this task and at least two kinds of approaches are proposed and widely used so far. First, statistical techniques and models such as the multiple linear regression model [30], ridge regression [31], structural equation modeling [32], random forest algorithms [33], redundancy analysis [34], geographically weighted regression models [35], and geographical detectors [36,37] are often used to explore the relationship between vegetation indices and their climatic and anthropogenic drivers to quantity the relative role of their influences. These models need both climatic and human-related variables as input data, but the indicators of human activities are usually difficult to obtain and often inconsistent with the spatial resolution of climatic indicators. Hence, the residual trend analysis has been proposed and then widely adopted to distinguish the impacts of climate variations and human activities on vegetation changes without the requirement of human-related data [38,39]. For this residual approach, a regression model is usually established between vegetation indices and climatic variables to obtain the climate-induced vegetation change, and then the residuals of the model are calculated to represent the human-induced vegetation change. Considering the difficulty in obtaining comprehensive data of anthropogenic variables, the residual trend approach has become increasingly popular and has been proven to be an effective method to separate the relative effects of climate variability and human disturbance on vegetation greenness or productivity in many studies [15,16,40,41,42,43]. However, most of these studies focus on arid and semiarid regions where vegetation generally has distinct phenological dynamics and plant growth is obviously limited by climatic factors. Few studies focus on warm-humid regions with vastly different climatic conditions to arid and semiarid ecosystems, so attributing vegetation dynamics in warm-humid regions is beneficial to enrich our understanding of the characteristics and mechanisms of vegetation changes [9].

Guangdong Province, situated in South China with a warm and humid climate, is an important carbon sink with vegetation cover at the forefront of China. Climate extremes in Guangdong are much higher than in most other areas of China [44], and Guangdong Province is one of the most vulnerable regions to climate change [45]. During the past decades, Guangdong Province has experienced a soaring economy and population boom, and land use has changed dramatically, thereby affecting ecosystem services and habitat quality [46,47]. Current studies indicate that the vegetation in Guangdong Province changed significantly during the urbanization process [33]. As such, Guangdong Province is a typical area for studying the combined effects of climate change and human disturbance on vegetation. By taking land use/cover change (LUCC) to represent human activities, previous studies emphasized human activities as the main driving force of vegetation changes and explored how the interaction of climate variations and human activities affected vegetation dynamics in Guangdong Province [9,33]. LUCC is one of the most direct manifestations of human activities, but it cannot reflect all aspects of human activities. Many other human-related activities without LUCC, including the increase in population intensity, improvement of low-quality forests, and management of green space, also have tremendous impacts on vegetation variations. Thus, considering only LUCC may underestimate the influence of human activities on vegetation changes.

This study took Guangdong Province as a case study to (1) investigate the spatiotemporal trends of vegetation and climate, (2) explore the relationships between vegetation and climate, and (3) quantify the relative contributions of climate change and anthropogenic activities to vegetation dynamics. To achieve the above aims, Moderate Resolution Imaging Spectroradiometer (MODIS) Normalized Difference Vegetation Index (NDVI) data during 2001–2020 were collected, and multiple analytical methods and techniques including Theil–Sen median trend analysis with a Mann–Kendall significance test, partial correlation analysis, and residual trend analysis were used to detect vegetation variations and their attribution. The results of this study are expected to separate the climate-induced and human-induced vegetation changes in the warm-humid Guangdong Province, which is helpful to evaluate the effectiveness of ecological conservation and restoration measures over the last 20 years, thereby providing guidance for future policy making.

2. Materials and Methods

2.1. Study Area

Guangdong (GD) Province (20°09′–25°31′N, 109°45′–117°20′E) is located in the south of China and covers a land area of ~179,700 km² (Figure 1a). Guangdong Province contains 21 cities and can be divided into 4 eco-regions according to ecosystem types and geographical characteristics (www.ecosystem.csdb.cn, accessed on 5 May 2020): the northern Guangdong (NG) eco-region with mid-subtropical mountainous and hilly evergreen broad-leaved forest, the middle Guangdong (MG) eco-region with south-subtropical mountainous and hilly evergreen broad-leaved forest, the southern Guangdong (SG) eco-region with tropical monsoon forest and rainforest, and the Pearl River Delta (PRD) eco-region with south-subtropical urban and suburban agriculture (Figure 1b). The MG eco-region can be further divided into the western part (West MG) and the eastern part (East MG). Vegetation in these eco-regions is different, thus the averaged NDVI over 2001–2020 displays distinct spatial patterns (Figure 1c). The warm-humid subtropical climate dominates GD Province with abundant light, heat, and water resources. According to the spatially interpolated meteorological records as described in Section 2.2.2 below, the mean annual temperature (MAT), precipitation (MAP), and radiation (MAR) ranges during 2001–2020 are 11.77–24.65 °C, 1546–2688 mm, and 4147–5308 MJ/m², respectively (Figure 1d–f).

2.2. Data Preparation

2.2.1. NDVI Data

The MODIS NDVI data (MOD13A1) for the period of 2001–2020, with the temporal resolution of 16 days and the spatial resolution of 500 m, were collected from the Land Processes Distributed Active Archive Center (LP DAAC) within the U.S. National Aeronautics and Space Administration (NASA) Earth Observing System Data and Information System (EOSDIS) (https://lpdaac.usgs.gov, accessed on 15 April 2020). Since the data might be contaminated by clouds and haze, we used the Savitzky–Golay filter in the TIMESET software package of MATLAB R2020b to smooth the annual cycle of NDVI data and improve the data quality of the NDVI time series [48]. Then, the maximum value composite (MVC) method was used to generate monthly and yearly NDVI datasets.

2.2.2. Climate Data

Daily temperature (°C), precipitation (mm), and sunshine duration (hour) records of 128 meteorological observation stations within and nearby Guangdong Province during 2001–2020 were provided by the China Meteorological Data Service Center (http://data.cma.cn, accessed on 13 December 2021). Incorrect and suspicious data were first removed via data cleaning and quality control and the percentage of missing data was kept within 10% for every meteorological station [49]. The missing data were imputed via the “mice” package in RStudio [50]. Then, these daily climatic records were composited to monthly datasets for spatial interpolation. The ANUSPLIN 4.2 software was used to interpolate meteorological point data into raster surfaces with the spatial resolution in line with NDVI pixels [51]. The sunshine duration (hour) was converted to radiation (MJ/m²) via the method recommended by the Food and Agriculture Organization [52]. Finally, the yearly average temperature and total precipitation and radiation were calculated.

2.2.3. Eco-Region Data

The eco-region data were obtained from the datasets of ecosystem assessment and ecological security in China (www.ecosystem.csdb.cn, accessed on 5 May 2020). Based on the eco-environmental problems, ecosystem sensitivity, and ecosystem services in China, eco-regions were delineated according to natural conditions such as landforms, hydrothermal conditions, vegetation characteristics, ecosystem types, and geographical features via multiple methods including spatial overlay technique, correlation analysis, and expert consultation. There are 50 eco-regions in China, and Guangdong Province involves 4 of them (Figure 1b).

2.3. Methods

To distinguish the relative contributions of climate change and human activities to vegetation dynamics via residual trend analysis, it is essential to understand the characteristics of vegetation and climate trends and the relationships between vegetation and climatic variables in advance. Using the time series of NDVI and climatic variables (temperature-TEM, precipitation-PRE, radiation-RAD) during 2001–2020, the Theil–Sen median trend analysis and linear regression analysis were first adopted to reveal the features of vegetation and climate changes at both pixel and regional scales, and then partial correlation analysis was used to explore the relationships between vegetation and climatic factors, and finally residual trend analysis was applied to separate the climatic and anthropogenic contributions to vegetation changes in this study (Figure 2).

2.3.1. Theil–Sen Median Trend Analysis with a Mann–Kendall Significance Test

The Theil–Sen median trend analysis was used to analyze the yearly trend of NDVI or climatic variables (temperature, precipitation, radiation) at pixel scale, and the Mann–Kendall test was used to detect the significance level of the changing trend. As a robust non-parametric technique, this method is insensitive to outliers and missing noise values [53,54,55], so it has been widely employed for long time-series trend analyses, especially in vegetation dynamic and climate change studies [14,15,56,57].

First, the slope of the Theil–Sen median was calculated by the Theil–Sen estimator:

$$S_x=Median\left(\frac{x_j-x_i}{j-i}\right), 2001\le i<j\le 2020$$

(1)

where S_x is the median from the set of slopes of the NDVI or climatic variables time series; x_i and x_j are the NDVI values or climatic variables of year i and year j, respectively. S_x > 0 indicates an increasing trend, and S_x < 0 indicates a decreasing trend for the period of 2001–2020.

Then, the significance of the trend was determined by the Mann–Kendall test:

$$S={\sum }_{j=1}^{n-1}{\sum }_{i=j+1}^{n}sgn\left(x_j-x_i\right)$$

(2)

$$sgn\left(x_j-x_i\right)=\left\{\begin{array}{cc}1,& x_j-x_i>0\\ 0,& x_j-x_i=0\\ -1,& x_j-x_i<0\end{array}\right.$$

(3)

$$\delta \left(S\right)=\frac{n\left(n-1\right)\left(2n+5\right)}{18}$$

(4)

$$Z=\left\{\begin{array}{c}\frac{S-1}{\sqrt{\delta \left(S\right)}}, S>0\\ 0, S=0\\ \frac{S+1}{\sqrt{\delta \left(S\right)}}, S<0\end{array}\right.$$

(5)

where n is the length of years; sgn is a sign function. Z value, calculated by S statistics, was used to test the significance of the changing trend.

For a given significance level α, |Z| > Z_1−α/2 denotes the changing trend of the time series is significant [15,57]. This study selected α = 0.05 as the significance level to classify the trend magnitude (

Table 1. Trend magnitude classification of the Theil–Sen median trend analysis.

Sx	Z	Trend Magnitude
Sx > 0	|Z| > 1.96	Significant increase
Sx > 0	|Z| ≤ 1.96	Slight increase
Sx < 0	|Z| > 1.96	Significant decrease
Sx < 0	|Z| ≤ 1.96	Slight decrease

).

2.3.2. Linear Regression Analysis

The linear regression analysis based on the ordinary least-squares method was applied to detect the overall interannual trend of NDVI during 2001–2020 at regional scale. The slope coefficient was calculated as below:

$$slope=\frac{n\times {\sum }_{t=1}^{n}t\times {NDVI}_t-\left({\sum }_{t=1}^{n}t\right)\left({\sum }_{t=1}^n{NDVI}_t\right)}{n\times {\sum }_{t=1}^n{t}^2-{\left({\sum }_{t=1}^{n}t\right)}^2}$$

(6)

where slope is the changing trend of NDVI; n is the number of years; NDVI_t is the averaged NDVI value of year t for the study area or eco-regions. The positive slope value indicates an overall increasing trend and the negative slope value implies an overall decreasing trend in the study period.

2.3.3. Partial Correlation Analysis

Partial correlation analysis was performed to explore the degree of association between every individual climatic factor (temperature, precipitation, or radiation) and monthly NDVI while controlling the potential effects of the other two climatic variables. The partial correlation coefficient can reflect a more intrinsic correlation compared with the simple correlation coefficient by eliminating the interactive influences of climatic factors on vegetation, so the partial correlation analysis has been widely used to investigate climatic driving forces on vegetation growth [11,15,24,36]. We used the “ppcor.test” function in the “ppcor” package in software RStudio version 2023.03.0 to explore the partial correlations between NDVI and each climatic variable. The calculation of the partial correlation coefficient was as follows:

$$r_{ab,c~p}^2=\frac{R_{a(b,c,\dots ,p)}^2-R_{a(c,\dots ,p)}^2}{1-R_{a(c,\dots ,p)}^2}$$

(7)

where r_ab,c~p is the partial correlation coefficient between variables a and b without the influence of variables (c~p); R_a(b,c,…,p) is the Spearman correlation coefficient of regression analysis between variable a and (b~p); R_a(c,…,p) is the Spearman correlation coefficient of regression analysis between variable a and (c~p). The non-parametric Spearman’s rank correlation has no assumptions of the distribution of variables, so the Spearman correlation coefficient rather than the conventional Pearson correlation coefficient was used in this study because the climatic variables may not be normally distributed. According to the significance test, the partial correlation between NDVI and each climatic variable was classified into four types: significant negative (SN), non-significant negative (NN), non-significant positive (NP), and significant positive (SP).

In order to screen the dominant climatic factor driving vegetation change, the absolute partial correlation coefficient values of the three climatic variables (temperature, precipitation, radiation) were compared, and the climatic factor with the maximum absolute partial correlation coefficient value was considered to have the greatest importance [58,59].

2.3.4. Residual Trend Analysis

Since both climate variations and human activities have tremendous impacts on vegetation in Guangdong Province [14,33,45], residual trend analysis was used to distinguish the relative impacts of climate change and human activities on vegetation changes [38]. This method assumes that the unexplained variation in the model between vegetation and climatic variables is attributed to anthropogenic activities [60]. In other words, this assumption regards human activities as the only remaining cause except for climate change in driving vegetation dynamics. Hence, the impact of human activities on vegetation dynamics can be indirectly quantified by establishing a model between climatic variables and vegetation indices and extracting the residual of the model. Heat, water, and light are widely regarded as the three main climatic limitations for vegetation growth [17], and vegetation responds sensitively to monthly climate extremes in Guangdong Province [45]. Thus, we selected temperature, precipitation, and radiation as climatic factors and established a multiple linear regression model between monthly NDVI and climatic factors to predict the climate-driven NDVI in this study. Then, the human-induced NDVI was calculated by the model residual (NDVI_res), which was the difference between the observed NDVI (NDVI_obs) and predicted NDVI (NDVI_pre). The model was expressed as follows:

$${NDVI}_{pre}=\alpha \times TEM+\beta \times PRE+\gamma \times RAD+\epsilon$$

(8)

$${NDVI}_{res}={NDVI}_{obs}-{NDVI}_{pre}$$

(9)

where TEM, PRE, and RAD refer to the climatic factors of temperature, precipitation, and radiation, respectively; α, β, γ, and ε refer to the model coefficients.

2.3.5. Relative Contribution under Various Scenarios

Based on the residual trend analysis as mentioned above, the impacts of climate change and human activities on vegetation variations were separated and measured by NDVI_pre and NDVI_res, respectively. The trends of both NDVI_pre and NDVI_res were then analyzed by linear regression analysis. An increasing/decreasing trend of NDVI_res indicates positive/negative effects of human activities on NDVI improvement and so does the NDVI_pre. Thus, various scenarios based on the slopes of NDVI_obs (slope_obs), NDVI_pre (slope_pre), and NDVI_res (slope_res) were developed to determine the relative contributions of climate change and human activities in vegetation change [39] (

Table 2. Relative contributions of climate change and human activities to vegetation change [39].

Vegetation Trend	Scenario		Relative Contribution (%)		Contribution Classification
	slopepre	sloperes	Climate Change	Human Activity
Increase	>0	<0	100	0	Climate change induced vegetation improvement (CI).
(slopeobs > 0)	<0	>0	0	100	Human activities induced vegetation improvement (HI).
	>0	>0	$\frac{{slope}_{pre}}{{slope}_{obs}}\times 100$	$\frac{{slope}_{res}}{{slope}_{obs}}\times 100$	Both climate change and human activities induced vegetation improvement (BI).
Decrease	<0	>0	100	0	Climate change induced vegetation degradation (CD).
(slopeobs < 0)	>0	<0	0	100	Human activities induced vegetation degradation (HD)
	<0	<0	$\frac{{slope}_{pre}}{{slope}_{obs}}\times 100$	$\frac{{slope}_{res}}{{slope}_{obs}}\times 100$	Both climate change and human activities induced vegetation degradation (BD)

).

3. Results

3.1. Spatiotemporal Changes in the NDVI and Climatic Variables

The spatial pattern of yearly NDVI trends analyzed by the Theil–Sen median trend analysis showed that an increasing trend was observed in most land areas (61.87% for significant increase and 25.13% for slight increase) of Guangdong Province from 2001 to 2020 (Figure 3a). Areas with a decreasing NDVI trend were mainly concentrated in the PRD, SG, and East MG eco-regions, in which the significant decrease trend (4.03%) primarily happened in the PRD and East MG eco-regions, and the slight decrease trend (8.97%) chiefly occurred in the SG eco-region.

Area statistics illustrated that NDVI trends varied greatly among eco-regions (Figure 3b). Overall, the total area percentage of increased NDVI ranked as follows: NG > MG > SG > PRD. The NDVI in approximately two-thirds of the MG (67.72%) and NG (66.00%) eco-regions exhibited a significant increasing trend, and the proportion for the SG and PRD eco-regions was 49.71% and 37.43%, respectively. In contrast, many more areas (14.30%) in the PRD eco-region experienced a significant NDVI decrease than in the other eco-regions, while this percentage was only 1.08% for the NG eco-region. The significant decrease areas for the MG and SG regions were comparable, accounting for 4.30% and 3.76%, respectively.

From the result of linear regression analysis, the regional area-averaged NDVI generally presented a significant increasing trend (p < 0.01, n = 20) during 2001–2020 for the whole study area but with different changing rates and processes for the four eco-regions (Figure 3c). The changing rate of interannual NDVI was highest for the MG and NG eco-regions with a similar rate of 0.0030–0.0031 per year and the minimum NDVI value appeared in 2004. The SG eco-regions had a slightly lower changing rate at 0.0027 per year, and the NDVI increased sharply from 2003 but fluctuated between 2005 and 2014 and increased again after 2014. The increasing rate for the PRD eco-region was lowest at 0.0015 per year with the NDVI increasing between 2006–2016 but decreasing before 2006 and after 2016.

During 2001–2020, temperature, precipitation, and radiation displayed diverse fluctuations but non-significant yearly trends in most areas (Figure 4). A warming trend was observed in most areas especially for several regions in the southern PRD, northern NG and eastern East MG eco-regions with more evident warming, but only a small region exhibited a significant changing trend (Figure 4a). Wetting mainly happened in most of the NG and PRD eco-regions, while the majority of the MG and SG eco-regions became drier, particularly for some regions at an annual decreasing rate of more than 16 mm (Figure 4b). Radiation decreased in the whole study area except for a few areas in the SG eco-region with a slightly increasing trend, and the significant changing trend was merely observed in a small part of the central NG eco-region (Figure 4c).

3.2. Relationships between Climatic Variables and the NDVI

The partial correlation between the monthly NDVI and climatic factors (temperature, precipitation, and radiation) was mapped to determine the impacts of climate change on vegetation variations (Figure 5). A significant positive correlation between the NDVI and temperature was observed in most of the study area (80.73%), while the area showed a significant negative correlation only accounting for 0.32%. Regions with non-significant correlations between the NDVI and temperature were mainly scattered in the PRD and West MG eco-regions (Figure 5a). Contrary to temperature, the impact of precipitation on the NDVI was not significant in most areas (71.63%). The NDVI and precipitation showed significant positive correlations in the SG eco-region and parts of the West MG and northern NG eco-regions, whereas the NDVI and precipitation showed significant negative correlations in the East MG and eastern NG eco-regions (Figure 5b). Similar to PRE, regions showing a significant negative correlation between the NDVI and radiation were mainly distributed in the East MG eco-region. There was a noticeable tendency of transition from significant negative correlation to significant positive correlation between the NDVI and radiation along the east-west gradient, except for the negative correlation at the southern edge of the SG eco-region. The proportion of regions presenting a significant positive correlation between the NDVI and radiation was 30.35%, which was widely situated in the PRD, West MG, southern NG, and northern SG eco-regions (Figure 5c).

The absolute values of the partial correlation coefficients between temperature, precipitation, radiation and the NDVI were compared to determine the dominant climatic factor in driving vegetation changes (Figure 6). The area percentage of the three dominant climatic factors was calculated for each eco-region and the whole study area. Results indicated that vegetation was predominantly influenced by temperature (73.42%), especially for the East MG and northern NG eco-regions. Areas with precipitation as the dominant factor were mainly distributed in the southern SG eco-region, accounting for 32.98% of the SG eco-region, while the impact of precipitation on the other three eco-regions was relatively small. Vegetation dynamics in the northern PRD, West MG, southern NG, and northern SG eco-regions were primarily dominated by radiation.

3.3. Contributions of Climate Variations and Human Activities to NDVI Change

Vegetation responded differently to climate change and human activities based on their relative contribution rates calculated by the residual trend analysis (Figure 7). For climate change, regions with a negative contribution to the NDVI were mainly distributed in the southern NG, southern and northern West MG, eastern East MG, massive areas of the PRD, and a few areas of the SG eco-regions, accounting for 31.18% of the study area, against 68.82% for the regions with a positive contribution (Figure 7a). The contribution of climate change for most areas was small (−25% < contribution rate < 25%), while regions with a relatively large contribution (contribution rate < −50% or contribution rate > 50%) were scattered sparsely throughout all the eco-regions with an area proportion of only 1.57%. For human activities, regions with a negative contribution to the NDVI increase merely occupied 8.32%, and most (87.74%) of these regions were concentrated in the PRD and eastern East MG eco-regions with a contribution rate of less than −75% (Figure 7b). Conversely, human activities in most areas (91.68%) exhibited a positive contribution, in which 96.98% presented a high contribution rate (>75%). Overall, human activities dominated (contribution rate < −50% or contribution rate > 50%) 98.43% of the vegetation changes in the whole study area.

The contribution type of climate change and human activities to vegetation dynamics was classified under various scenarios (Table 2, Figure 8). Vegetation increase or decrease caused by climate change only happened in a few regions, accounting for 0.56% and 0.14% of the study area, respectively. Human-induced vegetation improvement (23.42%) mainly occurred in the NG and MG eco-regions, while human-induced vegetation degradation (4.71%) was concentrated in the PRD and eastern East MG eco-regions. Vegetation variations in most (68.26% for NDVI increase and 2.92% for NDVI decrease) of the study area were influenced by both climate change and human activities. The area proportion of vegetation improvement caused by both climatic and anthropogenic contributions was highest in all four eco-regions, occupying 64.71%, 69.62%, 79.19%, and 60.04% of the NG, MG, SG, and PRD eco-regions, respectively. In contrast, regions of vegetation degradation induced by the combination of climate change and human activities had a relatively small area and were distributed dispersedly, except for the concentrated distribution in the PRD eco-region.

4. Discussion

4.1. Vegetation Trends and Their Climatic Drivers

This study found that vegetation growth in most areas of Guangdong Province displayed an increasing trend during 2001–2020, and the decreasing trend was mainly concentrated in the urban areas and peripheries in the PRD and eastern East MD ecoregions (Figure 3), in agreement with previous studies [9,33]. However, these studies demonstrated a non-significant increasing or even decreasing trend of vegetation greenness in most area of northern Guangdong [9,33], contrasting with our results that vegetation improved significantly in the NG eco-region. This is partially due to the more extended study period of our study (2001–2020) compared with the prior work (2001–2013). Our results showed that the NDVI in the NG eco-region fluctuated highly before 2011, and has sharply increased in recent years (Figure 3c). Another possible reason is the difference in zonation. The eco-regions in this study were delineated based on natural conditions, while previous studies used the administrative division map. Thus, the boundary definition of northern Guangdong (NG eco-region in this study) was different. The vegetation trend in many regions of the northern NG eco-region in this study was not significant either, similar to these previous studies. The results of this study indicated that trends of climatic variables (temperature, precipitation, radiation) and their effects on vegetation changes were spatially different (Figure 4 and Figure 5). First, significant positive correlations between the NDVI and temperature were observed across most areas (Figure 5a) and temperature was the dominant climatic factor for 73.42% of the study area (Figure 6), especially the northern NG and East MG eco-regions, in line with previous studies [9,33]. These regions are mostly mountainous areas with higher latitudes and elevations, and the temperature is relatively low in comparison with other eco-regions (Figure 1d), so vegetation growth is more sensitive to low temperatures [61]. The increasing temperature is widely believed to prolong the growing season with an earlier spring and a delayed autumn [62,63], thus promoting vegetation greenness in these heat-limited regions. Second, radiation had a predominant influence on vegetation in most areas of the PRD, West MG, and southern NG eco-regions (Figure 6), where the NDVI and radiation were significantly and positively correlated (Figure 5c). The deficiency of sunshine time may be the reason why vegetation tends to be limited by light in these regions (Figure 1f). Conversely, negative effects of radiation on vegetation were extensively detected in the East MG eco-region and the southern edge of the SG eco-region. This is possibly because of photoinhibition as sufficient radiation exists in these regions and strong light inhibits vegetation photosynthesis [64]. Last, compared with the widespread impacts of temperature and radiation on vegetation variations, regions with precipitation as the dominant climatic factor only accounted for 6.42% of the study area (Figure 6). Since precipitation in Guangdong Province is usually abundant, drought is generally not the stressor of vegetation growth, which was also concluded in a global study showing that tropical forests are highly sensitive to temperature and cloudiness rather than water availability [19]. The influence of precipitation on vegetation showed distinct north-south patterns with negative impacts in most areas except for the SG eco-region (Figure 5b), which was also observed in previous studies [9,45]. Specifically, significant negative correlations between the NDVI and precipitation were mainly distributed in the eastern NG and East MG eco-regions. This is presumably because these regions have high susceptibility to geohazards [65], and the excessive precipitation is likely to reduce photosynthetic activity via aggravating soil erosion and even triggering floods [66]. In contrast, the significant positive effect of precipitation on vegetation greenness was primarily situated in the SG eco-region. The precipitation in the SG eco-region is relatively low with a decreasing trend (Figure 1e and Figure 4b), thus the drier environment may cause stress on vegetation growth. In addition, the SG eco-region has the best heat and light resources in the study area, and the persistent high temperature and radiation may accelerate evapotranspiration and exhaust soil moisture, thereby causing adverse impacts on vegetation [67,68,69]. Precipitation is believed to alleviate these unfavorable effects and promote vegetation growth [45,70].

4.2. The Dominant Role of Anthropogenic Activities in Vegetation Change

As Guangdong Province belongs to the warm-humid subtropical climate region with relatively high temperature, abundant precipitation, and sufficient radiation (Figure 1) and these climatic variables did not change significantly in most areas over 2001–2020 (Figure 4), vegetation dynamics are largely independent of climate change. Instead, anthropogenic activities play a dominant role in vegetation variations (Figure 7 and Figure 8), also supported by other studies [7,9,33,36], suggesting that the residual trend analysis is an effective method in the study area.

This study indicated that 98.43% of vegetation changes were predominantly induced by human activities with the absolute contribution rate >50% (Figure 7), which is much higher than the value of 79.40% in a previous work [33]. One possible reason is that our results included the contribution of both climate change and human activities just with a higher contribution rate of human activities than climate change. In contrast, the previous work only referred to the explanation rate of human activities. Additionally, the previous work took LUCC to represent human activities, which might neglect human-related activities without land use conversions and thus underestimate the anthropogenic contribution to vegetation changes. By comparison, we used the residual trend analysis in this study which assumed that vegetation variations were induced by human activities after excluding the effect of climate change. Human activities have both negative and positive effects on vegetation greenness. Land cover changes caused by urbanization are widely regarded as the main reason for vegetation degradation in Guangdong province [9,14,33,71]. Our results also found that the NDVI decreased around the urban areas, particularly for the PRD eco-region and the eastern part of the East MG eco-region (Figure 3a), where urban expansion mainly occurred [47]. Apart from urban regions, vegetation in most other areas was improved by human activities, which is principally due to a series of ecological protection and restoration programs since 1999 [72]. In this study, we found that NDVI had the largest increasing trend in the MG eco-region (Figure 3c). This is probably because conversions from croplands and grasslands to forests via afforestation mainly happened in eastern and western Guangdong [9], and fast-growing forest species were selected in these large-scale reforestation programs [73], which significantly contributed to vegetation greenness. In addition, urban greening policies [26], the urban heat island effect [74], the fertilization effect of carbon dioxide [75], and nitrogen deposition [76] may also promote vegetation enhancement in built-up areas. thus, vegetation changes were more complicated in the PRD eco-region (Figure 3c).

4.3. Uncertainties and Challenges

Our findings highlighted the predominant role of human activities in vegetation change in Guangdong Province, but some limitations and uncertainties remain. The residual trend analysis used in this study assumes that vegetation changes are driven by either climate change or human activities. However, some situations may conflict with this assumption including (1) the difficulty in accurately quantifying the relationship between climatic variables and vegetation indices, and (2) where other factors rather than climate change and human activities are the controlling factors of vegetation variations. Firstly, the residual trend analysis based on the multiple linear regression model has been well applied to disentangle the relative contributions of climate change and human activities to vegetation dynamics in arid and semiarid regions [39,40,42,43], where vegetation growth is usually water-limited and linearly correlated to precipitation change. However, vegetation growth in warm and humid regions such as Guangdong Province is generally not inhibited by climatic conditions and the response of vegetation to climate change may be nonlinear [77,78]. The multiple linear regression model may not be able to accurately reflect the complicated relationship between vegetation responses and climate variations. A recent study attempted to combine a nonlinear model with residual trend analysis to better reveal the nonlinear vegetation responses to climatic variables [41]. Thus, the method applicability and selection need more exploration based on a deeper understanding of vegetation–climate relationships. Moreover, other factors such as atmospheric carbon dioxide concentration [62], nitrogen deposition [76], soil condition [25], agricultural intensification [26], and exotic plant invasion [79] may also largely affect vegetation greenness and productivity. For example, the invasion of the exotic grass red brome in the northeastern region of the Mojave Desert played a vital role in vegetation greening on post-fire landscapes [79]. The carbon dioxide fertilization effects explained most of the greening trends in the tropics [62]. Hence, the contribution of human activities to vegetation changes separated by the residual trend analysis is likely to be overestimated when there are other important factors besides climate change and human activities. Additionally, the lag effect of climate change on vegetation growth was not considered in this study but has been considered in many other studies [9,58,80]. However, these studies only analyzed the lag time of climatic influences on vegetation, but how to integrate the lag effect into the attribution analysis of vegetation dynamics is still a challenging work.

5. Conclusions

Based on data of the NDVI time series and climatic variables (temperature, precipitation, radiation) during 2001–2020, this study analyzed the spatiotemporal patterns of vegetation trends, explored the response of vegetation to climate variations, and distinguished the relative climatic and anthropogenic contribution to vegetation dynamics in Guangdong Province of South China. Vegetation enhancement happened in most areas, while vegetation degradation was mainly concentrated around urban areas. Vegetation trends varied largely among the four eco-regions in Guangdong Province with the highest rate of change in the MG eco-region and lowest rate of change in the PRD eco-region. Changes in temperature, precipitation, and radiation presented distinct spatial patterns over the past 20 years, but their changing trends were not statistically significant in most areas. Vegetation responded differently to climate change with significant positive correlations to temperature in most areas, and negative correlations in the east but positive correlations in the west were generally observed between vegetation and both precipitation and radiation. Temperature was the dominant climatic factor in most areas, while regions with radiation or precipitation as the dominant factor were mainly concentrated in the center-west or the southern edge, respectively. The contribution of climate change to vegetation variations was much smaller compared with human activities. Human-induced vegetation greening widely happened in the NG and MG eco-regions, while human-induced vegetation browning was concentrated in the PRD and eastern East MG eco-regions. Vegetation dynamics in most areas were influenced by both climate change and human activities, although the contribution rate of human activities was much higher than climate change. The results of this study highlight the effectiveness of ecological policies in vegetation restoration, which is of great significance for future ecosystem management and conservation.