Spatiotemporal Variation of Vegetation Productivity and Its Feedback to Climate Change in Northeast China over the Last 30 Years

Gross primary productivity (GPP) represents total vegetation productivity and is crucial in regional or global carbon balance. The Northeast China (NEC), abundant in vegetation resources, has a relatively large vegetation productivity; however, under obvious climate change (especially warming), whether and how will the vegetation productivity and ecosystem function of this region changed in a long time period needs to be revealed. With the help of GPP products provided by the Global LAnd Surface Satellite (GLASS) program, this paper gives an overview of the regional feedback of vegetation productivity to the changing climate (including temperature, precipitation, and solar radiation) across the NEC from 1982 to 2015. Analyzing results show a slight positive response of vegetation productivities to warming across the NEC with an overall increasing trend of GPPGS (accumulated GPP within the growing season of each year) at 4.95 g C/m2. yr−2 over the last three decades. More specifically, the growth of crops, rather than forests, contributes more to the total increasing productivity, which is mainly induced by the agricultural technological progress as well as warming. As for GPP in forested area in the NEC, the slight increment of GPPGS in northern, high-latitude forested region of the NEC was caused by warming, while non-significant variation of GPPGS was found in southern, low-latitude forested region. In addition, an obvious greening trend, as reported in other regions, was also found in the NEC, but GPPGS of forests in southern NEC did not have significant variations, which indicated that vegetation productivity is not bound to increase simultaneously with greening, except for these high-latitude forested areas in the NEC. The regional feedback of vegetation productivity to climate change in the NEC can be an indicator for vegetations growing in higher latitudes in the future under continued climate change.


Introduction
Gross Primary Productivity (GPP) of a terrestrial ecosystem refers to the amounts of organic compounds accumulated by green vegetations through assimilating carbon dioxide from the atmosphere via photosynthesis and a series of internal physiological processes [1][2][3]. GPP can partially offset anthropogenic CO 2 emissions [3][4][5][6], so that vegetation plays a critical role in the global carbon balance. Statically, vegetation can store about 500 Gt of carbon, which accounts for 20% of all carbon in the terrestrial biosphere [7].
To date, studies [8,9], using leaf area index (LAI) as a reference, found a significant global vegetation greening; however, the question of whether the vegetation productivity and ecosystem function varies under the obvious vegetation greening has yet to be unraveled. Compared with widely discussed leaf area (LAI or some vegetation indexes) primarily by the Amur, Argun, and Ussuri Rivers, from North Korea to the south by the Yalu and Tumen Rivers, and from the Inner Mongolian Autonomous Region to the west by the Greater Khingan Range.  [38] in 2015).

GPP Datasets
GPP data used in this paper were provided by the GLASS program [39]. Generally, the GPP can be estimated using the LUE model as follows: where PAR is the incident photosynthetically active radiation (MJ m −2 ) per time period (e.g., day or month), fPAR is the fraction o PAR absorbed by the vegetation canopy, εmax is the potential LUE (g C m −2 MJ −1 APAR) without environment stress, and f is a scala varying from 0 to 1 and represents the effects of temperature, moisture, and other environmental conditions on LUE. Across from mid to high latitude (approximately 39-53 • N), the NEC has a large variation in annual average temperature from 10 • C in south to −5.5 • C in north. Moreover, located at the northernmost region of China, the NEC contains the coolest region in China.
The accumulated temperature greater than 0 • C is around 2000~4200 • C, and the average temperature during summer is 20-25 • C. Influenced by the Eastern Monsoon of China, the NEC has a great variation in annual precipitation, decreasing from wet southeast (900 mm) to semi-arid northwest (400 mm), and a majority of precipitation occurs in summer (July to September).
There are three mountains: the Greater Khingan Mountain (GKM), the Lesser Khingan Mountain (LKM), and the ChangBai Mountain (CBM) (Figure 1a) that form important natural barriers for the NEC, so that the NEC can be divided into four zones (one plain and three mountains). The heartland of the NEC is the Northeast China Plain (Figure 1a, including Sanjiang Plain, Songnen Plain, and Liaohe Plain), which mainly located at the western NEC and grows a lot of crops ( Figure 1b) for its rich, deep, black soil. Owning these plains, the NEC is one of the most important centers of commercial agriculture in China for both food grains (wheat, rice, and maize) and economic crops (soybean and sugar beets). The rest three mountainous zones also play a crucial role in ecosystem function. There are a great amount of forests growing in these three mountains, e.g., the LKM and the CBM are widely covered with broadleaved deciduous forests, and the GKM owns the special bright coniferous forest which is the southern boundary of the Boreal forests of Eurasia.

GPP Datasets
GPP data used in this paper were provided by the GLASS program [39]. Generally, the GPP can be estimated using the LUE model as follows: where PAR is the incident photosynthetically active radiation (MJ m −2 ) per time period (e.g., day or month), fPAR is the fraction of PAR absorbed by the vegetation canopy, ε max is the potential LUE (g C m −2 MJ −1 APAR) without environment stress, and f is a scalar varying from 0 to 1 and represents the effects of temperature, moisture, and other environmental conditions on LUE. GLASS-GPP products were estimated using an improved EC-LUE (Eddy Covariance Light Use Efficiency) model, which is derived by satellite data and eddy covariance measurements [16,17]. There are four kinds of input parameters: NDVI, PAR (photosynthetically active radiation), air temperature, and the Bowen ratio of sensible-to-latent heat flux. NDVI datasets are provided by Advanced Very High Resolution Radiometer (AVHRR), which cover a long period starting from 1979. Net radiation (Rn), air temperature (T), relative humidity (Rh), and photosynthetically active radiation (PAR) from the MERRA (Modern Era Retrospective-Analysis for Research and Applications) archive [40] were used in GPP retrieval. The MERRA dataset was provided by the NASA. Its theory uses a fixed global atmospheric model and data assimilation system to analyze the historical satellite and conventional data records into a continuous global gridded dataset. The spatial resolution was 0.5 • latitude by 0.6 • longitude.
The GLASS-GPP datasets provides datasets covering from 1982 to 2015 with a spatial resolution of 0.05 • and a temporal resolution of 8 days. We downloaded 1564 scenes with 46 scenes per year, and then used the Savitzky-Golay [41] filter to correct errors and fill gaps. After a series of re-projection, extraction, we acquired the available time series of GPP of the NEC. Aiming to analyzing inter-annual variation of GPP, we calculated accumulated GPP during the growing season (April to October) of each year (denoted as GPP GS , hereafter) and, finally, acquired the 34-year GPP GS time series for further analysis.

Auxiliary Datasets
For the attribution analysis of the spatiotemporal GPP variation, there are two kinds of datasets: the meteorological datasets and human intervention indices.
The meteorological datasets used in this study came from the China Meteorological Forcing Datasets [42], which can be acquired from the Cold and Arid Regions Science Data Center [43]. Specifically, the air temperature dataset was produced by merging observations from 740 meteorological stations of China Meteorological Administration (CMA) into the corresponding Global Land Data Assimilation System version 1 (GLDAS-1) dataset [44]. The precipitation dataset was produced by combining the observations from the 740 CMA stations, the Tropical Rainfall Measuring Mission (TRMM) 3B42 precipitation products [45], and the GLDAS-1 rainfall rate data. The downward solar radiation was constructed by correcting the Global Energy and Water Cycle Experiment-Surface Radiation Budget shortwave radiation dataset [46] with radiation estimates from meteorological station observations. We downloaded reanalyzed meteorological data (including daily average temperature, daily precipitation rate, and daily solar radiation) with monthly temporal resolution and spatial resolution of 0.1 • × 0.1 • . First, daily average temperature and daily radiation data are averaged within the growing season of each year, and the daily precipitation rate were translated into accumulated precipitation within the growing season of each year to obtain yearly accumulated precipitation from 1982 to 2015. In the next step, the meteorological datasets were downscaled to 0.05 • to match the spatial resolution of GPP datasets for correlation analysis. Finally, time series of yearly average temperature, yearly accumulated precipitation, and yearly average radiation within the growing season from 1982 to 2015 were acquired for further analysis.
Human intervention indices are also considered in this study since crops are heavily influenced by human activities, so the assessment of the impacts of human intervention on crops is necessary. This study considered three indices, which can quantify the impacts of human interventions, i.e., Cropland Area (CA), Crop Yield per Area (CYA), and Irrigated Area (IA). All these indices were provided by the National Agriculture Database (Statistics Bureau of China) [47]. The detailed impacts of human intervention on GPP are shown in Section 4.1 as an additional analysis.

Mann-Kendall Trend Test
The Mann-Kendall (MK) test is a rank-based non-parametric test for detecting a monotonic trend in a time series. The method was originally used by Mann (1945) [48], and the test-statistic distribution was subsequently derived by Kendall (1975) [49]. It is frequently used for trend detection in hydrologic research [50]. In this study, we use the MK test to detect monotonic trends in the GPP data series as the test is more powerful than ordinary parametric trend tests for non-normally distributed series [51].
The MK test statistics (S) for a series x 1 , x 2 , . . . , x n can be given by where N is the length of the dataset and x i and x j denote the data values at times i and j. A negative value of S indicates a decreasing trend, while a positive value indicates an increasing trend [52]. Statistical significance of the trend is checked using the variance of the MK statistics, while N > 10, which is calculated as where P is the number of tied groups, the summary sign (Σ) denotes the summation over all tied groups, and t i is the number of data in the ith (tied) group. If the tied groups are not available, this summation process is excluded from the equation. The standard Z value of the test statistics can be estimated by the following formula.
If the computed standard Z value is greater than the critical value of the standard normal distribution at the specified significance level (α), it is considered to be significant trend, and the null hypothesis is rejected.

Pearson Correlation
Pearson correlation is a standard parametric method used to analyze the relationship between two variables. It can be expressed as follows.
where x i and y i (i = 1, 2, . . . , n) are time series for two variables, x and y are multi-year average values of two variables, and r is the correlation coefficient. A value of r > 0 indicates a positive correlation between two variables, whereas r < 0 suggests a negative correlation. The absolute value of r can represent the closeness of two variables, i.e., the larger the value of r, the closer the two variables are and vice versa. All the correlation coefficients were conducted using the 0.05 significance level.

Standardized Multivariate Linear Regression Model
Unlike the Pearson correlation method, the Standardized Multivariate Linear Regression (SMLR) method uses Gaussian normalized values in place of the original ones. Thus, the effects of changes in different factors can be compared. The SMLR model separates the relative contribution of each factor. The coefficient before each variable represents its importance in partitioning the SMLR model, and the fitting coefficients (R) of SMLR represent the combined effects of driven factors to GPP. The SMLR model is expressed as follows: where ED i (i = 1, 2, . . . , n) are influencing factors; GPP and ED 1 are multi-year mean values of GPP GS and influencing factors, respectively; σ GPP and σ ED are standard deviations of GPP GS and influencing factors, respectively, during the study period; and β 1 , β 2 , . . . , β n are regression coefficients that represent the importance to GPP GS . Factors with a larger coefficient are assumed to contribute more to the variation in GPP GS , and the one with largest coefficient was identified as the dominant factor to GPP GS in our study area during 1982-2015.
Before using the SMLR model, the Variance Inflation Factor (VIF) is calculated to exam multicollinearity among the influencing factors. The F-test is conducted on the model, while the t-test is conducted on the coefficients before each influencing factor. Significant tests were those below the 0.05 significance level.

Spatiotemporal Distribution and Variation of GPP
Figure 2a reveals an obvious increasing trend of multi-year average GPP GS from north to south and from west to east in the NEC, which is consistence with the spatial distribution of vegetation types and climate. The productivity of forests ranges from 700 to 1500 g C/m 2 . yr −1 , which is obvious larger than that of crops (50-700 g C/m 2 . yr −1 ), and the productivity of forests in warm-humid lower latitude regions is about 1000-1500 g C/m 2 . yr −1 , larger than that of forests in cold-dry regions (700-1000 g C/m 2 . yr −1 ).
Spatial distribution of inter-annual variation of GPP GS in the NEC is shown in Figure 2b which owns an overall increasing trend, though great spatial heterogeneity exists. Specifically, the western NEC growing with crops has the largest GPP GS increment (over 11 g C/m 2 . yr −2 ), followed by a slight increment of GPP GS in the GKM and LKM (around 3-7 g C/m 2 . yr −2 ) covered with forests, while the variation of GPP GS is not significant for the majority regions of CBM with sporadic pixels showing increases. Statistically, around 66.4% of the regions in the NEC has an increasing trend in GPP GS , among which the area of croplands accounts for 67.3%.
We further identified the inter-annual variation of zonal GPP GS at five regions: the whole NEC, plain, GKM, LKM, and CBM. There is an overall increasing trend of the whole NEC (the change ratio (denoted as CR, hereafter) = 4.95 g C/m 2 . yr −2 ; p < 0.01) from 1982 to 2015 as shown in Figure 3a, which is mostly attributed to the growth of crops in the plain (Figure 3b). Crops growing in the NEC have a stable increase with the largest amplitude (CR = 7.47 g C/m 2 . yr −2 ; p < 0.01), while, in the other three forested regions, GPP GS shows a much smaller increasing trend (LKM: CR = 5.35 g C/m 2 . yr −2 , p < 0.01; GKM: CR = 3.28 g C/m 2 . yr −2 , p < 0.05; CBM: CR = 0.75, p > 0.05) with quite a bit of inter-annual fluctuation over the past three decades.   Figure 4 shows a warming, wetting, and dimming trend across the whole NEC, and the impacts of this climate change on vegetation GPP GS are analyzed in the following sections.  Table 1 indicates that GPP GS of the NEC is related with the three meteorological factors (daily average temperature, precipitation, and daily solar radiation), especially the obviously warming; nevertheless, this correlation has a great spatial heterogeneity. In plains, warming (Figure 4a) can improve the productivity of crops, while precipitation ( Figure 4b) and solar radiation (Figure 4c) do not have significant impacts. In mountainous areas, positive impacts of warming on forest productivity only exist in high-latitude GKM and LKM, and solar radiation comes into positive effect on GPP GS in GKM and CBM. Precipitation has a negative effect on GPP GS of forests according to the correlation results, nevertheless, those negative effects were not significant when eliminating the coupling effects between temperature and radiation with precipitation as partial correlation coefficients display. The difference in correlation and partial correlation results manifest that in the NEC, precipitation indirectly impacts vegetation growth by affecting temperature and solar radiation.  The contribution of each meteorological factors on GPP GS can be quantified using the SMLR method, and thus the dominant meteorological factors for GPP GS variation over the last three decades can be identified as shown in Figure 5. In mountainous areas, the combination of the three meteorological factors has great impacts on forest productivity. Figure 5 shows that 70.5% of the GPP GS variation in LKM is caused by these three factors among which temperature is the dominant factor. Meteorological factors can explain 61.0% of the variation in GPP GS of CBM where solar radiation is the dominant factor. In GKM, climate effects accounts for 84.0% of the variation in GPP GS , and both temperature and solar radiation are equally important for vegetation productivity. Contrary to productivity of natural forests in mountainous areas, crops growing in plain are less related with the three meteorological factors, although warming still has positive effect on vegetation productivity in some areas. Furthermore, the linear combination of temperature, precipitation, and solar radiation can only explain 40% of the significant GPP GS increase in crops, indicating that some other facilitating factors exist for crop productivity increments in the NEC.

Relationship between GPP and LAI
The comparison between GPP and LAI is noteworthy. We download GLASS-LAI products with an 8-day temporal resolution and 0.05 • spatial resolution from 1982 to 2015 [39,53,54] and then pre-processed it to a 34-year time series of mean value of LAI within the growing season (denoted as LAI GS ) across the whole NEC (Figure 6a). The detection of LAI trend across the NEC shows a dramatic increase in LAI GS across the whole NEC over the last 30 years (Figure 6b). It is the same as the widely reported worldwide greening [8,[55][56][57]. Different from the significant increase in LAI GS across the whole NEC as Figure 6b shows, there is a relatively smaller increment of GPP GS in forested region across the NEC (the slight increment in the GKM and LKM, and even no increments of GPP GS in the CBM); quantitatively, the absolute change ratio (ACR) and relative change ratio (RCR, defined as ACR * study period multi−year mean value ) of GPP GS and LAI GS in each part of the NEC are displayed in Table 2. In forested regions, the increments of LAI GS are larger than that of GPP GS , especially in CBM, while, in croplands, the RCR of GPP GS is twice that of LAI GS . The relationship between GPP GS and LAI GS seems to vary in different regions. According to correlation analysis (Figure 7), GPP GS is correlated with LAI GS in relatively high-latitude forests (GKM: R = 0.652, p < 0.05; LKM: R = 0.691, p < 0.05) and parts of cropland; however, in mid-latitude CBM, there is non-significant correlation (R = 0.334; p > 0.05), indicating that GPP is not always changing simultaneously with LAI.

Differences of Crops and Forests in GPP GS Trend
In both pixel or regional scales, the inter-annual variation trend of GPP GS in the NEC (Figures 2 and 3) indicates that the crops growing in plains provide the majority interannual increments of terrestrial ecosystem productivity while forests in mountainous area contribute less in productivity increasing over the last 30 years. Furthermore, meteorological factors are not much important for crops in the NEC ( Figure 5). Previous studies [58][59][60] identify two main reasons for the variation of crop productivity, i.e., climate change and agricultural technological progress. Fang et al. (2000) [60] concluded that agricultural technological progress plays a much more important role than climate for rice yield in Heilongjiang, which may be suited for the whole NEC.
The relationship between GPP GS of crops and the three human intervention indices (IA, CA, and CYA) are displayed in Figure 8. All these indices are highly correlated with GPP GS where the IA (R = 0.897, p < 0.05) owns the largest correlation, followed by the CA (R = 0.896, p < 0.05), and finally, the CYA (R = 0.784, p < 0.05). The increments of vegetation productivity of crops in the NEC benefits greatly from human intervention (e.g., fertilization, irrigation, and cropland area expansion), instead of meteorological factors identified in Table 1, which is in line with Yuan's findings [61].  Song et al. (2017) [62] defined vegetation productivity as a combination of net assimilation rate, LAI, and photosynthesis duration. That is to say, apart from the LAI, net assimilation rate and photosynthesis duration are also decisive for GPP.

Different Trends of GPP and LAI
Firstly, the net assimilation rate is related with the vegetation type as well as vegetation age. There is a great spatial difference of vegetation type across the whole NEC as Figure 1b and Qian et al.'s study [63] display, so that the relationship between GPP and LAI has great spatial heterogeneity. In addition, the vegetation age is different with time went by which the correlation between GPP and LAI is also influenced. Second, combined with the correlation analysis between GPP and meteorological factors ( Figure 5), we found that the CBM is a radiation-dominant region. so, the decreasing trend in solar radiation (or sunshine duration) (Figure 4c) may offset the effect of increasing LAIGS (Figure 6b), leading to non-significant variation in GPPGS over the 30 years.
Actually, LAI and GPP describes vegetation growth from two different aspects. LAI is defined as one half of the total green leaf area per unit horizontal ground surface area [64], which aims to quantify the amount of leaf area in an ecosystem [65], so in the NEC, the increasing LAIGS reflects that the leaf quantity and vegetation coverage of the canopy in the NEC increases, while the small variation in GPPGS indicates that vegetation's contribution to regional carbon balance has not changed during the three decades of interaction with its environment across the NEC.

Uncertainties and Future Works
This paper has revealed the slight increasing of GPPGS in the NEC from 1982 to 2015, and simultaneously identified the effects of climate change on GPPGS. Nevertheless, some uncertainties still exist for this analysis.
First, we used the GLASS-GPP products which can provide vegetation productivity from 1982; however, the spatial resolution of GLASS-GPP (around 5 km) is rather rough, and thus some variations of GPP may be overlooked, especially in heterogeneous regions, leading to overestimation or underestimation of the trend of vegetation productivity. All in all, human intervention contributes a lot to vegetation productivity of crops, which greatly facilitates vegetation growth of crops directly leading to the differences of crops and forests in GPP trend. Song et al. (2017) [62] defined vegetation productivity as a combination of net assimilation rate, LAI, and photosynthesis duration. That is to say, apart from the LAI, net assimilation rate and photosynthesis duration are also decisive for GPP.

Different Trends of GPP and LAI
Firstly, the net assimilation rate is related with the vegetation type as well as vegetation age. There is a great spatial difference of vegetation type across the whole NEC as Figure 1b and Qian et al.'s study [63] display, so that the relationship between GPP and LAI has great spatial heterogeneity. In addition, the vegetation age is different with time went by which the correlation between GPP and LAI is also influenced. Second, combined with the correlation analysis between GPP and meteorological factors ( Figure 5), we found that the CBM is a radiation-dominant region. so, the decreasing trend in solar radiation (or sunshine duration) (Figure 4c) may offset the effect of increasing LAI GS (Figure 6b), leading to non-significant variation in GPP GS over the 30 years.
Actually, LAI and GPP describes vegetation growth from two different aspects. LAI is defined as one half of the total green leaf area per unit horizontal ground surface area [64], which aims to quantify the amount of leaf area in an ecosystem [65], so in the NEC, the increasing LAI GS reflects that the leaf quantity and vegetation coverage of the canopy in the NEC increases, while the small variation in GPP GS indicates that vegetation's contribution to regional carbon balance has not changed during the three decades of interaction with its environment across the NEC.
In summary, although LAI GS indicates a greening trend, ecological productivity may not have corresponding variation. GPP GS may be a more suitable parameter to reflect variation in vegetation growth and vegetation's impacts on carbon cycle [1][2][3][4][5][6][7].

Uncertainties and Future Works
This paper has revealed the slight increasing of GPP GS in the NEC from 1982 to 2015, and simultaneously identified the effects of climate change on GPP GS . Nevertheless, some uncertainties still exist for this analysis.
First, we used the GLASS-GPP products which can provide vegetation productivity from 1982; however, the spatial resolution of GLASS-GPP (around 5 km) is rather rough, and thus some variations of GPP may be overlooked, especially in heterogeneous regions, leading to overestimation or underestimation of the trend of vegetation productivity.
Another uncertainty comes from correlation analysis between GPP GS and three meteorological factors. In this study, we chose the fixed growing season (April to October) of each year as our study period, and then considers the impacts of the three meteorological factors (daily average temperature, precipitation, and daily radiation) within the growing season, nevertheless, these three factors may influence vegetation growth during the whole year. Take temperature as an example, winter temperature also influences spring phenological events of the next year [66]. Therefore, using fixed growing season period as well as considering annual average temperature during growing season to analyze the effects of warming on vegetation productivity may bring some uncertainties. In the next step, we should consider a variable growing season length, such as using temperature thresholds (e.g., 0 and 5 • C) as a sign for growing season, so that a more comprehensive impact of meteorological factors on GPP can be taken into account.
Finally, uncertainty may exist when quantified the effects of each meteorological factors on GPP using SMLR models. The SMLR model assumed that the combined effect of meteorological factors is linear, actually, there is a more complex relationship between GPP GS and meteorology, in addition, there is a coupling effect between vegetation and atmospheric dynamics, that is to say, apart from the effect brought by climate change on vegetation, the variation of vegetation would impact climate simultaneously as a lot of studies [67][68][69] have revealed. The relationship between weather and vegetation is too complex to quantify using a statistic-based model. So, a land process models based on biochemical mechanisms, such as CABLE [70], ORCHIDEE [71], and LPJ [72], is needed for more accurate analysis. Apart from reducing uncertainties mentioned above, one important thing that needs further exploration is the relationship between GPP GS and LAI GS . This paper found that the relationship between GPP GS and LAI GS across the NEC has large spatial heterogeneity, indicating that the widely reported greening is not bound to have great contribution to improve carbon sink capability. Whether this relationship between GPP GS and LAI GS found in this paper is specific to the NEC or a common phenomenon in the terrestrial ecosystem is still to be seen.

Conclusions
This study focuses on the regional feedback of vegetation in the Northeast China (NEC) to the obvious climate change using GLASS GPP products as indicators. Results shows that, from 1982 to 2015, GPP GS across the NEC has an overall increasing trend (4.95 g C/m 2 . yr −2 , p < 0.01) within which great spatial heterogeneity exists. Vegetation productivity of crops distributing in western NEC is the smallest but has the largest increment (7.47 g C/m 2 . yr −2 , p < 0.01). On the contrary, vegetation productivity of forests is obvious larger than that of crops but only has slight increments in northern NEC (GKM: CR = 3.28 g C/m 2 . yr −2 , p < 0.05; LKM: CR = 5.35 g C/m 2 . yr −2 , p < 0.01) and non-significant variation in southern NEC (CBM).
Correlation analysis between GPP GS and meteorological factors indicates that warming can improve vegetation productivity of crops as well as forests located in temperaturelimited high latitude (GKM and LKM), while this positive effect vanishes in mid-latitude CBM. In addition, for crops of the NEC, apart from the facilitating effect of warming, human intervention plays a decisive role in vegetation productivity increasing over the last 30 years, which directly leads to the difference of crops and forest in GPP GS trend.
Contrary to the slight or even non-significant variation of vegetation productivity in forested region, an obvious vegetation greening is found across the whole NEC. Correlation analysis between GPP GS and LAI GS over the last three decades illustrates that the simultaneously increment of GPP GS and LAI GS only exists in high-latitude NEC (GKM and LKM), while in other region, facilitating effect of greening on GPP GS may be offset by other factors, such as the shortage of solar radiation in CBM. The different variation trend of GPP GS and LAI GS reminds us that apart from focusing on vegetation growth dynamics (greening or browning), more attention should be paid on vegetation's feedback to climate as well as its role in carbon cycle.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author or from the reference website.