Multi-spatiotemporal heterogeneous legacy effects of climate on terrestrial vegetation dynamics in China

ABSTRACT Investigating vegetation–climate interactions is critical for understanding behavioral patterns of terrestrial ecosystems and formulating food security strategies. However, the multi-spatiotemporal legacy effects of climate on terrestrial vegetation remain unclear. In this study, we examined the dynamic trends of vegetation distribution and climatic factors at multiple temporal and spatial scales. Moreover, using cross-wavelet transform, wavelet coherence transform, and partial correlation analysis, a paradigm framework was established to determine the multi-spatiotemporal legacy effect of climate on vegetation in China in 2000–2019, as well as the response of different vegetation types. The results indicate a significant greening trend in China, accompanied by a warming and wetting pattern over the past 20 years. The phase difference of the wavelet coherence transform in the time-frequency domain revealed remarkable legacy effects and regional variations, indicating a complicated relationship between vegetation and climate. Meanwhile, different vegetation types exhibited heterogeneous responses of legacy effect of precipitation and temperature on multi-spatiotemporal scales; moreover, the lag time in spring was shorter than that in summer and autumn. The average legacy effect on different vegetation types was approximately 1–2 months. Therefore, the heterogeneity of the legacy effects is a complicated process of dynamic variation, which can be summarized as the comprehensive characteristics of vegetation response to climate with regional discrepancy, ecosystem category, and multi-temporal scale. These findings advance our understanding regarding the preference of vegetation to hydrothermal conditions across biomes and ecosystems and provide a future framework for elucidating the dynamic response of vegetation to other more complex factors in this warmer world. Furthermore, our results emphasize that multi-spatiotemporal legacy effects should be incorporated into the vegetation–climate interaction model and the formulation of macro-environmental management policies.


Introduction
Vegetation, an indispensable component of terrestrial ecosystems, plays a fundamental role in the exchange of global carbon, nitrogen, water, and energy fluxes (Gao et al. 2020;Shi et al. 2020b;Sun et al. 2015;Wang et al. 2011); moreover, vegetation is recognized as a crucial indicator of biological responses to environmental factors (Menzel and Fabian 1999). Generally, the dynamic variations of vegetation are closely related to climate-induced factors (Pearson et al. 2013), which provide the specific hydrothermal conditions for vegetation growth and result in significant effects on the structure and function of ecosystems (Doherty et al. 2010). Indeed, this phenomenon has been revealed at regional, national, and global scales (Chen et al. 2019;Detsch et al. 2016). Furthermore, terrestrial ecosystems are likely to be more profoundly affected by ongoing global climate change (Garcia et al. 2014;Seddon et al. 2016). Therefore, it is essential to comprehensively quantify the impacts of climatic constraints on terrestrial vegetation and improve our understanding of the response mechanism of vegetation dynamics for accurate projection of future vegetation trends, ecosystem evolution, global change, and effective implementation of macroenvironmental management.
Precipitation and temperature have been widely recognized as the dominant climatic parameters driving variations in vegetation dynamics (Bégué et al. 2011;Gao et al. 2020;Richardson et al. 2018;Wu et al. 2015). Meanwhile, the responses of vegetation growth to climatic variations are often complicated and vary greatly in different regions (Bianchi, Villalba, and Solarte 2019). For example, precipitation has a more significant influence in arid and semi-arid regions, whereas the dominant factor at high northern latitudes is temperature (Cleland et al. 2007;Jeong et al. 2011;Xu et al. 2013). In light of global climate change, warming can lead to greening trends in ecosystems with sufficient water supply, but it also produces moisture stress on the vegetation in arid areas (Emmett, Renwick, and Poulter 2018;Sulla-Menashe, Woodcock, and Friedl 2018). In addition to the spatial heterogeneity of vegetation responses, quantifying the diversity characteristics of vegetation dynamics on temporal scales is also particularly complex in the comprehensive understanding of vegetation-climate interactions. Furthermore, the generalization of knowledge information on hierarchical spatial scales could be understood as morphological abstraction at different perspective levels. Meanwhile, the response of vegetation to environmental variations on a finer scale can provide a wealth of physiological and ecological information. It is well documented that terrestrial vegetation in China has been exhibiting a statistically significant greening trend (Piao et al. 2006Zhu et al. 2016). Under such a background, it is particularly urgent to study the variation trends of terrestrial vegetation growth and interaction relationships with climatic constraints on multi-spatiotemporal scales.
Legacy effects are caused by asynchrony between the vegetation phenological cycle and climatic change (Bianchi, Villalba, and Solarte 2019;Ehleringer et al. 1991;Zhou et al. 2020a). Incorporating the legacy effect of climate on vegetation and its potential mechanisms into the climate model can not only improve the prediction performance, but also facilitate the development of food security and environmental management strategies (Zhang et al. 2015a). Previous studies have recognized the importance of the legacy effect of climate on vegetation growth and recommended its use as a benchmark for coupled vegetation-climate models (Chen et al. 2014;Vicente-Serrano et al. 2013;Wen et al. 2019;Wu et al. 2015). Nevertheless, existing studies still adopt synchronous meteorological and vegetative parameters in climate-vegetation interactions (Jong et al. 2013), which affect the performance of models, thereby hindering future projections. Some studies have examined the legacy effect of climate on vegetation growth on a regional scale (Anderson et al. 2010;Zhao et al. 2020), but the conclusions are not always consistent. Because of the spatial heterogeneity of vegetation growth and meteorological factors, the paradigm framework based on a regional scale provides only the tip of the iceberg, and cannot provide a comprehensive perception of the legacy effect. Meanwhile, considering the diversity of ecosystem behavior, the temporal heterogeneity of vegetation growth should be understood as the discrepancy in its preferences for hydrothermal conditions, which is a dynamic and complicated process. The response mechanism of vegetation growth to hydrothermal conditions in specific periods is often obscure when the spatial heterogeneity of the legacy effect is described alone. Therefore, a paradigm framework for summarizing the multi-spatiotemporal heterogeneity of legacy effects in the context of complex terrain is still lacking.
It is well documented that wavelet transform, including continuous wavelet transform, crosswavelet transform and wavelet coherence transform, can extract localized intermittent periodicities by extending time series to time-frequency space (Grinsted, Moore, and Jevrejeva 2004). Thus, the method has been extensively employed to analyze the characteristics of geophysical and ecological time series (Jia et al. 2018;Vargas et al. 2012;Zhou et al. 2020b). For instance, previous studies have reported that the cross-wavelet transform is effective in exploring the interaction relationship (i.e. their common power and relative phase in timefrequency space) between time series (Grinsted, Moore, and Jevrejeva 2004;Sun et al. 2020). Meanwhile, some studies have analyzed the multiscale characteristics of climatic factors using wavelet transform (Chang et al. 2018;Tan, Gan, and Shao 2016). Furthermore, the wavelet coherence transform can be useful for determining the degree of association between two non-stationary series (Gallegati 2018), and can capture the quasi-periodic component (Grinsted, Moore, and Jevrejeva 2004). Hence, the advantage of the wavelet coherence transform is that it can investigate the dynamic characteristics and phase difference between two time series in the time and frequency domains. Recently, existing studies have utilized the wavelet coherence transform to examine the coherence and lags between climate and phenology (Detto et al. 2018). Given that crosswavelet transform and wavelet coherent transform can acquire the common power region, coherence, and phase angle of time series across different scales, introducing them to analyze the multi-spatiotemporal heterogeneity of legacy effects would be promising for revealing the complicated responses of vegetation to climate change.
The primary objective of our study is to examine vegetation-climate interactions ascertaining the multi-spatiotemporal heterogeneity of the legacy effect of climate on vegetation growth in China, and to quantify the response of different vegetation types to the legacy effect. Our study is based on the normalized difference vegetation index (NDVI) dataset from the Terra Moderate Resolution Imaging Spectroradiometer (MODIS) NDVI datasets (MOD13A3) and meteorological data from the Climate Research Unit (CRU) Time-series (TS) data. Specifically, our research focused on the following questions: (1) How do terrestrial vegetation dynamics and climatic constraints (temperature and precipitation) change on multi-temporal (monthly, seasonal, and inter-annual) scales in China? (2) Can crosswavelet transform and wavelet coherence transform be used to explain the legacy effect of climate on vegetation growth in the time-frequency domain?
(3) Do multi-spatiotemporal characteristics exist in the legacy effect of climate on vegetation? (4) Does the legacy effect of climate vary with vegetation discrepancy?

NDVI data
Vegetation dynamics are described using satellitebased NDVI datasets, which have been reported in previous studies (Fensholt et al. 2012;Hmimina et al. 2013;Liu et al. 2018). The MODIS NDVI (MOD13A3) Version 6 datasets were used to detect and characterize vegetation trends for the period from 2000 to 2019 (Didan 2015). The datasets can be accessed from https://ladsweb.modaps.eosdis. nasa.gov/. The original 1-km spatial resolution datasets representing the monthly gridded product in the sinusoidal projection have been applied to evaluate the global greenness in previous studies (Asri et al. 2020). Moreover, previous studies utilized NDVI as a proxy for vegetation growth (Jiao et al. 2021;Wang et al. 2011;Zhou et al. 2020a), so this terminology (i.e. vegetation growth) was adopted in this study.

Land cover dataset
The MODIS Land Cover Type (MCD12Q1) Version 6 datasets were obtained in this study (Friedl and Sulla-Menashe 2019). The datasets were primarily used to establish the legacy effects of climate on different vegetation types. The MCD12Q1 products have five land-cover classification schemes. We used the dataset of International Geosphere-Biosphere Programme (IGBP) classification scheme containing 17 land-cover types with a spatial resolution of 500 m for the period from 2001 to 2019 (https:// ladsweb.modaps.eosdis.nasa.gov/). Because landcover data are not available in this product for the year 2000, it is replaced by land-cover data for the year 2001. In this study, we adopted the central pixel approach to match land-cover data and NDVI (Jong et al. 2013). To examine the response of different vegetation types to the legacy effect, unchanged regions that were defined as pixels and always presenting the same state during 2000-2019, were mainly considered ( Figure 1). In addition, the spatial distributions of six sub-regions were obtained from the Resource and Environment Science and Data Center (http://www.resdc.cn/Default.aspx).

Meteorological data
The monthly precipitation and temperature data were obtained from the CRU TS dataset version 4.04 with a spatial resolution of 0.5°×0.5° (Harris et al. 2014) (TS dataset are monthly gridded products [http://badc.nerc.ac.uk/data/cru/]). For the analysis to be consistent with vegetation dynamics, the datasets were restricted to the period of 2000-2019. This dataset was calibrated and updated using as earlier product. Previous versions of the dataset have been thoroughly assessed and used in global research Wang et al. 2011;Wu et al. 2015). To maintain consistency with the spatial resolution of the NDVI dataset, precipitation and temperature data were resampled to a resolution of 1 km using the bilinear interpolation method (Bond-Lamberty et al. 2007;Ershadi et al. 2013;Peng et al. 2019).

Methods
Our study provided a framework to examine vegetation-climate interactions ascertaining the multispatiotemporal heterogeneity of the legacy effect of climate on vegetation growth in China, as well as the response of different vegetation types. The overall procedures implemented in this study are shown in (Figure 2).

Trend analysis of vegetation dynamics and climatic factors
The least squares method was used to determine the linear trend of the vegetation and meteorological parameters at the pixel scale in 2000-2019 (Fu et al. 2019). The trend function is: where P i refers to the monthly vegetation or meteorological parameters of the year i and Tr is the linear regression slope. The vegetation and meteorological parameters exhibit a negative trend if Tr< 0 and a positive trend if Tr> 0.
To further illustrate the significance of the linear trend, we used the Mann-Kendall test as a nonparametric approach recommended by the World Meteorological Organization , to determine the presence of a significant trend (Dong et al. 2020). The Mann-Kendall test statistic (Z c ) is defined as (Geng et al. 2021): S ¼ where Z c < 0 indicates a decreasing tendency and Z c > 0 indicates an increasing trend (Zhou et al. 2020b).
When Z c j j > Z 1À α = 2 (i.e. Z 1À α = 2 is the standard normal deviation), the time series of vegetation or meteoro-logical parameters have a significant trend from 2000 to 2019. The significance level (α) of 0.05 was used in this study.

Partial correlation analysis
In probability theory and statistics, partial correlation measures the direction and strength between two random variables while removing the effect of the controlling covariates Z ). Thus, we performed pixel-level partial correlation analysis using NDVI, temperature, and precipitation during 2000-2019. Previous studies have suggested that the legacy effect of climate on vegetation growth is less than three months (Chen et al. 2014). Therefore, we used the NDVI of the current month and the meteorological parameters of the current and previous 1, 2, and 3 months for partial correlation analysis. The partial correlation coefficient is expressed as follows (Xie et al. 2019): ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffiffi where r yx�z refers to the partial correlation coefficient between y and x, when z is the control variable, and, r yx , r yz , and r xz represent the correlation coefficients among y, x and z. In this study, y represents NDVI, x denotes one of the climate parameters, and z represents another meteorological factor. The partial correlation coefficients ranged from −1 to 1. The higher the absolute value, the stronger is the correlation. A partial correlation analysis was performed using MATLAB R2017b for the present study.

Ascertaining legacy effect
The partial correlation coefficient PCC n (n ¼ 0; 1; 2; 3, denotes current and previous months) was first calculated between NDVI and climate factors (precipitation and temperature) based on Equation (5). Then, when PCC n corresponding to the maximum determination coefficient PCC 2 n was selected as the optimal PCC n (OPCC), and n was chosen as the optimal legacy effect (OLE) for vegetation responses to climate (Wen et al. 2019).
Permanent wetlands, urban areas, snow and ice, barren or sparsely vegetated and water bodies (i.e. IGBP classification scheme) were excluded when determining the legacy effects. Meanwhile, the determined duration of legacy effects was at a significance level of 0.05. Additionally, in winter (i.e. December, January, and February), vegetation coverage may be relatively low in many regions of China (e.g. Northern China and the Qinghai−Tibet Plateau). Analyzing the legacy effect of climate on vegetation is not significant at the pixel scale, and even contains some uncertainties. Therefore, the NDVI time period covers March-November, while the time span of meteorological parameters is from January to December (i.e. climate parameters in winter may affect vegetation in March) in the present study.

Wavelet power spectrums, cross-wavelet transform, and wavelet coherence transform
Continuous wavelet transform (CWT) measures the similarity between a signal and an analyzing function. The CWT is defined as follows (Martínez and Gilabert 2009): where ψ represents the mother wavelet containing a trade-off between time and frequency, ψ � refers to the complex conjugate of the mother wavelet, and s and u denote the scaling and translation parameters, respectively. Therefore, CWT reveals the variation characteristics (i.e. dominant oscillation modes) of NDVI and meteorological parameters from three dimensions (i.e. time, periodic components, and color). The wavelet power spectrum (W PS ), a useful tool for determining the distribution of energy, of the CWT can be expressed as follows (Detto et al. 2018): The cross-wavelet transform is applied to analyze the co-relationship (i.e. the power distribution of two signals) of two time series in the time-frequency domain (Lisiecki 2010;Markonis et al. 2018). For two series X and Y, the cross-wavelet transform can be expressed as follows (Shi et al. 2020a): where W X α; τ ð Þ represents the wavelet coefficient of X, and � denotes the complex conjugation. Not only does the cross-wavelet transform reflect the common highenergy region, but it also describes the phase structure and local characteristics of two time series. The crosswavelet transform is thus used to exhibit areas with high common power and to reveal the phase difference between NDVI and meteorological factors in this study.
Wavelet coherence reveals the coherence (i.e. a measure of the correlation between two signals) of the cross-wavelet transform in time-frequency space (Peng et al. 2020). The wavelet coherence (R 2 n ) is defined as (Grinsted, Moore, and Jevrejeva 2004): where S denotes the smoothing operator, in which the wavelet transform is first performed in time and then processed along the axis of the wavelet scale. The phase angle difference represented by arrows in the wavelet coherence transform reflects the hysteresis characteristics (i.e. legacy effect) of the two sequences, such as NDVI and climate parameters. Specifically, right-oriented arrows denote synchronous status and left-oriented arrows represent the anti-phase status. The upward arrows show the delay of the first factor by 90°, while the downward arrows indicate the lead of the first factor by 90° (Jia et al. 2018).

Multi-spatiotemporal trends of NDVI and climate factors
Our goal was to establish spatial variation trends of NDVI and climate parameters on multiple time scales (i.e. inter-annual, seasonal, and monthly scales). Vegetation, precipitation and temperature exhibited significant cycle characteristics, with periods of 8-16 months ( Figure 3). Trends in NDVI and meteorological factors showed remarkable spatial heterogeneity on the monthly time scale during 2000-2019 ( Figures S1-S3). The annual change rate of NDVI was concentrated in the range of −0.06/yr to 0.06/yr. Considering the regional discrepancy, the study area was divided into six sub-regions: South Central China (SC), Northwest China (NW), East China (EC), North China (NC), Northeast China (NE), and Southwest China (SW). Greening trends were located in the south of NC, east of NW and SW, and north of EC and SC in spring (Figure 4), and this trend increased significantly in summer and autumn in the north of China. Precipitation substantially increased in the south of SC in spring (Figure 4). In summer and autumn, areas with increased precipitation were concentrated in southern and northeastern China (slope>1.36 mm/ yr). The proportion of the area with increased precipitation was 51.77% in spring, 64.17% in summer, and 57.38% in autumn (Table S1). In addition, the seasonal mean temperature showed a noticeable increase. The proportion of the area under warming was 90.90% in spring, 67.09% in summer, and 60.85% in autumn (Table S2). The proportion of regions with no change in precipitation and temperature was less than 1%. The regions in northern China showed the most extensive increase in spring, which was statistically significant. In summer and autumn, the temperature decreased in NW (i.e. junction areas of Qinghai and Xinjiang) and NE, while it increased in other regions, especially the southeastern coastal areas and the Qinghai-Tibet Plateau (see significance level in Figure S4). Consistent with previous studies (Chen et al. 2019;Piao et al. 2015), a significant greening trend was observed in China, accompanied by a warming and wetting pattern over the past 20 years.

Legacy effect of climate on vegetation growth in China
NDVI, precipitation and temperature exhibited a synchronous increasing trend before July ( Figure  S5). Interestingly, precipitation and temperature peaked in July, while the NDVI peaked in August. Meanwhile, the drastic reduction in precipitation from August to September resulted in a significant decrease in NDVI after September. Hence, these findings preliminarily demonstrate that the legacy effect of climate on vegetation does exist. To further reveal the relationships between vegetation dynamics and meteorological parameters, the cross-wavelet transform and wavelet coherence transform reflecting the hysteresis characteristics of two sequences in different time domains were adopted (Grinsted, Moore, and Jevrejeva 2004). The results of the cross-wavelet transform indicated that the significantly positive phase relationships between NDVI and meteorological factors (i.e. precipitation and temperature) occurred in the period of 8-16 months (i.e. regions with high common power were observed) from 2000 to 2019 (Figure 5(a)). NDVI and climate parameters were not completely synchronized (i.e. incompletely horizontal right arrows). Moreover, the results of wavelet coherence transform showed that the major resonance periods between NDVI and precipitation were approximately 8-16 months during 2000-2019 (WTC: NDVI-Pre; Figure  5(b)). A significantly similar pattern in the WTC was also found between NDVI and temperature at the 8-16month band (WTC: NDVI-Tem; Figure 5(b)). Meanwhile, NDVI exhibited strong coherence (i.e. >0.8) with precipitation and temperature in the wavelet coherence transform at the major resonance periods. The legacy effect of climate on vegetation was further determined according to sloping arrows. It is worth noting that there are some high coherence peaks at time scales between 0 and 2 months. The diverse arrows also illustrate the complexity of the vegetation response to climate across different time scales.

Multi-spatiotemporal heterogeneity of legacy effect
NDVI reached its maximum in August. Nevertheless, precipitation and temperature reached the maximum in July, except for SC and EC (the largest precipitation in June, monthly average temperature >0°C; see Figure S6). We performed the same analysis including wavelet power spectrum, cross-wavelet transform, and wavelet coherence transform at the sub-regions as we did at the national scale (see Figure S7-8; Figure 6). Our findings demonstrated that the major resonance periods between NDVI and meteorological factors were approximately 8-16 months in all sub-regions during 2000-2019. In addition, there were some high peaks of coherence between NDVI and climatic factors at timescales of 0-2 months ( Figure 6). The legacy effect of climate on vegetation is still obvious in subregions, which is consistent with results at the national scale. Nevertheless, some coherence peaks exist at time-scales of 2-4 months in SC and EC. Hence, the legacy effect of climate on vegetation growth has multiple spatiotemporal characteristics.
Heterogeneity of the multi-spatiotemporal legacy effects was further quantitatively described at the pixel level and monthly time scale (Figure 7). In March and April, vegetation growth was significantly dependent on the temperature requirement in northern China and the Qinghai-Tibet Plateau. Thus, vegetation growth produced the greatest response to the temperature of the same month (i.e. no legacy effect was observed at The cross-wavelet transforms between monthly NDVI and climatic factors (i.e. Pre denotes precipitation and Tem denotes temperature); Horizontal axis corresponds to the original time series, the vertical axis represents the period (i.e. months), and the color represents the intensity of the change period. The color bar denotes the energy density. b. The wavelet coherence transforms between NDVI and climatic factors. The thick contour designates the 95% confidence interval against red noise. The thin black line is the cone boundary influenced by wavelet. The wavelet coherence ranges from 0 to 1 (i.e. 0 denotes the perfectly uncorrelated series, whereas 1 denotes the perfectly correlated series). The arrows represent the relative phases of the two sequences. Right-oriented arrows denote synchronous status, and the left-oriented arrows represent the anti-phase status. The upward arrows show the delay of the first factor by 90°, while the downward arrows indicate the lead of the first factor by 90°. The input parameters in XWT and WTC are the monthly average NDVI, monthly average precipitation, and monthly average temperature. Further details are described in .sections 2.1 and 2.3 lower temperatures). However, precipitation did not contribute positively to vegetation growth under low temperature conditions in early spring (negative correlation with precipitation in the previous month; see Figure S10). When the temperature meets the requirement of vegetation growth, water supply plays a dominant role in arid and semi-arid regions (see Figures S9-10). The legacy effect of precipitation is approximately 1-2 months in those regions (Figure 7). The heterogeneity of spatial distribution of legacy effects on monthly time-scales also describes the significance of precipitation to vegetation growth at important time nodes (e.g. March and June). In the eastern part of NW and the southern part of NC, which were mainly covered by grasslands and croplands, vegetation growth in April presented the greatest response to precipitation in the same month, and the legacy effect lasted until July. In the response of vegetation growth to precipitation, the proportion of pixels with three-month lagged duration was 17.58% (i.e. in July). Interestingly, in the eastern part of NW, the spatiotemporal pattern of the legacy effect of temperature on vegetation growth was similar to that produced by precipitation from May to July (Figures 7 and 8). Pixels with 1-month legacy effect accounted for 27.45% of the total pixels, followed by a 3-month delay, accounting for 26.57% (Figure 8, Temperature in July). This also indicates that vegetation growth in these areas is equally dependent on precipitation and temperature. In the NE sub-region, the legacy effect of precipitation in May continued until August and showed a 1-month time-lag in October. A large area of croplands is distributed in NE, and the legacy effect of climate may also The horizontal axis corresponds to the original time series, the vertical axis represents the period (i.e. months). The thick contour designates the 95% confidence interval against red noise. The thin black line is the cone boundary influenced by wavelet. The wavelet coherence ranges from 0 to 1 (i.e. 0 denotes the perfectly uncorrelated series, whereas 1 denotes the perfectly correlated series). The arrows represent the relative phases of the two sequences. Right-oriented arrows denote synchronous status, and the left-oriented arrows represent the anti-phase status. The upward arrows show the delay of the first factor by 90°, while the downward arrows indicate the lead of the first factor by 90°. Pre denotes precipitation and Tem denotes temperature. Input parameters in WTC are the monthly average NDVI, monthly average precipitation, and monthly average temperature in each sub-region during 2000-2019. Further details are described in .sections 2.1 and 2.3 be affected by the crop planting system. A later spring arrival and a shorter average growing season in NE makes this region adapted for only one cropping season per year (Liu et al. 2013), including one region for thermophilic crops and warm crops. Because of the requirements for a more thermal environment, it is sensitive to temperature and does not exhibit longer lag time. In southern China and NC, the legacy effect of precipitation on vegetation growth began to be more obvious after August. In comparison, the legacy effect of temperature on vegetation growth is obviously less remarkable than that of precipitation in NC. It should be emphasized that high temperatures can inhibit vegetation growth through the depletion of soil moisture in arid and semi-arid regions (July; see Figure S11). In the eastern part of SW, the temperature in August has a profound effect on vegetation growth during the following three months. Therefore, our results demonstrate that the multispatiotemporal heterogeneity of the legacy effect of climate on vegetation can explain the preference of vegetation growth to hydrothermal conditions.

Responses of vegetation types to legacy effects
The legacy effect of climate varies with vegetation type and temporal scale, and the lag time in spring is shorter than that in summer and autumn (Figure 9). After May, the legacy effect of climate on vegetation is more than one month at the national scale (based on mean value). Longer legacy effects of temperature are observed in some vegetation types, such as mixed forests, open shrublands, woody savannas and savannas (i.e. in November), with average lag times of 1.8 ± 1.2, 1.9 ± 0.5, 2.2 ± 1.1, and 2.1 ± 1.2 months, respectively (Figure 9(a)). This indicates the dependence of vegetation growth on temperature in the previous months. Croplands (0.8 ± 1.1 months) have periodicity and exhibit a short lag time for temperature in November. Closed shrublands are more sensitive to temperature in the early spring. Nevertheless, the legacy effect of precipitation is more common and has obvious patterns (based on the mean value at the national scale) when the temperature requirement for vegetation growth is optimal. (Figure 9(b)). For  instance, deciduous forests (needleleaf and broadleaf forest) show a remarkable response to precipitation during May-September. A longer lag time for precipitation is observed in closed shrublands. Hence, different vegetation types exhibit a remarkable inhomogeneity response of the legacy effect to precipitation and temperature on temporal scales as well.

Multi-spatiotemporal heterogeneity of legacy effect
The legacy effect of climate on vegetation is essential for understanding the vulnerability of ecological environments and should be considered when exploring the mechanisms of climate-vegetation interactions (Wu et al. 2015). Based on the monthly NDVI dataset, our findings demonstrated that the legacy effect of climate on vegetation growth had significant heterogeneity on a multi-spatiotemporal scale. Previous studies have shown that the legacy effect varies with vegetation type and climatic factors (Wu et al. 2015). Our study not only quantified the spatial heterogeneity of the legacy effect of climate on vegetation growth, but described the multi-temporal variations of lag response, providing an indispensable mechanism within the vegetation-climate interactions. The legacy effects of climate on vegetation should be understood as the temporal integration of environmental variability (Fabricante, Oesterheld, and Paruelo 2009). The advantage of the wavelet coherence transform is that abundant scale information and coherence coefficients with phase difference can be obtained in the time-frequency domain (Peng et al. 2020). Although similar coherence patterns for all subregions can be observed between NDVI and meteorological factors (i.e. precipitation and temperature) at the 8-16-month band, the right arrow (i.e. denoting the phase information) in some subregions is more slanted, such as SC, EC, and SW ( Figure 6). Interestingly, for sub-regions SC and EC, a more slanted right arrow was also observed in the NDVI-pre (i.e. coherence relationship between NDVI and precipitation) than in NDVI-tem. That is, the legacy effect of precipitation on vegetation may have a more significant relationship than temperature during the long period (8-16 months) in those sub-regions. The coherence relationship between vegetation and precipitation and temperature in SW may be more complicated at longer periods because a more slanted right-pointing arrow is presented. Meanwhile, discontinuous oblique phase relationships between NDVI and temperature were observed at periods of 2-4 months in SC and EC. Peaks of coherence at lower timescales (i.e. periods <2 months) indicate highly unstable relationships between NDVI and meteorological parameters, which include anti-phase, inphase, and leading to each other.
Vegetation growth relies on water from soil conditions that can be replenished by precipitation. Therefore, the mechanism of the legacy effect is the time gap between the occurrence of precipitation and water absorption by plant roots (Jobbágy, Sala, and Paruelo 2002). Nevertheless, soil temperature, organic matter decomposition, and nutrient availability determine the legacy effect of temperature on vegetation growth (Braswell et al. 1997). Some studies have indicated that precipitation is the dominant factor in arid and semi-arid regions (Nemani et al. 2003), and others have demonstrated that the lag time was approximately 1 month (Wu et al. 2015). However, our results showed that there was a significant negative correlation between vegetation growth in March and precipitation in previous months at high latitudes and high altitudes (i.e. Qinghai-Tibet Plateau, QTP; see Figure S10). This led us to ask whether precipitation contributed little to vegetation growth under low temperature conditions in these regions. Evidence has demonstrated that rising temperatures can promote vegetation growth in the plateau (Wang et al. 2012), while low temperatures limit vegetation growth (Piao et al. 2015). Vegetation growth requires more energy and temperature to cope with colder environmental conditions in early spring. Therefore, they may produce a mechanism to compensate for the negative effects by responding to temperature rather than precipitation (Slot and Kitajima 2015;Wan et al. 2009;Wen et al. 2019). Additionally, over the past 50 years, the temperature rise on the QTP has been significantly higher than the magnitude of the global average (Chen et al. 2015). Rising temperature is one of the most important driving factors for glacier melting in the QTP (Zhang et al. 2021). Previous studies have also reported that the snow melting effect may be the reason for higher NDVI values over the QTP during January-April You et al. 2020). Hence, the effect of snow and glacier melting caused by climatic warming on vegetation should not be ignored.
The spring temperature is the main driver of vegetation growth in northern China. Previous studies have reported that a temperature rise in spring can alter the growing season length and affect biomass production and vegetation growth (Mulder, Iles, and Rockwell 2017;Vitasse et al. 2009;Wang et al. 2019). Our findings also indicate that precipitation does not exhibit a legacy effect on vegetation growth until the temperature requirement is met. The monthly responses of terrestrial vegetation reflect the sensitivity and dependence on hydrothermal conditions at important time nodes. When the ideal temperature is reached, the water supply plays a dominant role in arid and semi-arid regions, in agreement with earlier results (Guo et al. 2018;Zhou et al. 2020a). Note that the ideal temperature is defined as the corresponding temperature when the relationship between vegetation and precipitation changes from negative to positive; precipitation has a legacy effect on vegetation in subsequent months. Nevertheless, a hightemperature environment may increase evaporation, reduce soil moisture, and ultimately inhibit vegetation growth Yan et al. 2013). This further demonstrates that the legacy effect should be determined by considering temporal variations. Furthermore, legacy effects should be understood as a dynamic process that varies with meteorological parameters and multiple temporal and spatial scales.

Discrepancy response of vegetation types to legacy effects
Given the interaction relationship between various ecological and physiological characteristics and bioclimatic states, legacy effects vary greatly across biomes and plant functional types . Our findings confirmed that the legacy effect of climate on different vegetation types presented remarkable heterogeneity. Consistent with previous studies (Wu et al. 2015), precipitation and temperature show inconsistent legacy effects on different vegetation types, and the same vegetation type also appears to have discrepant responses to precipitation and temperature. In this study, unchanged regions that were defined as pixels and always presenting the same state were mainly considered to examine the legacy effect of climate on different vegetation types during 2000-2019. Previous studies have reported that legacy effects of precipitation on vegetation may be related to structural changes in the ecosystem (Reichmann, Sala, and Peters 2013;Wen et al. 2019;Wu et al. 2015). The heterogeneity of the legacy effect elaborates the preference and sensitivity of vegetation types to hydrothermal conditions. Meanwhile, the length of the legacy effect duration reflects the dependence of vegetation growth on precipitation and temperature in the previous months. Nevertheless, the discrepancy of vegetation types in physiology, structure, and adaptability may lead to significant differences in responses to precipitation and temperature (Pan et al. 2013;Wen et al. 2019;Xu et al. 2015).
Climate is one of the factors that causes spatial heterogeneity in light use efficiency, which is a momentous plant attribute (Hilker et al. 2010). For example, precipitation has a significant legacy effect on the light use efficiency (LUE) of most biomes, suggesting that precipitation is an indirect mechanism affecting vegetation growth during the long-term period (Zhang et al. 2015b). Previous studies have reported that LUE exhibits significant differences in different forest cover types (Ahl et al. 2004). We found that deciduous forests showed the longest response to legacy effects of precipitation during July-September, while evergreen forests exhibited shorter lags. LUE may be reduced for evergreen needleleaf forests compared with deciduous forests because of the energy used for maintenance of foliage and biomass (Still, Randerson, and Fung 2004). Furthermore, the maximum light use efficiency (LUE m ) is defined as the light energy utilization efficiency when all other factors (e.g. temperature and water) are optimal . It is well documented that LUE m decreases under suboptimal temperature and water deficit, and varies with vegetation type and environmental conditions (Madani et al. 2014;Madani, Kimball, and Running 2017). The LUE m of some typical vegetation types in China has been simulated in existing studies (Zhu et al. 2006). The simulation results show that broadleaf forests have a larger LUE m than shrublands and grasslands. Our analysis indicated that grasslands were more sensitive to precipitation variations and presented smaller lags than forestland and shrubland during May-July. This is because grasses draw water from the upper layers of the soil, while shrubs and forests get water from deeper layers (Bianchi, Villalba, and Solarte 2019;Sepúlveda et al. 2018). Meanwhile, water stress on vegetation is always related to temperature (Brzostek et al. 2014). Additionally, the legacy effect of precipitation on closed shrublands appears to be a periodic characteristic. However, an obvious discrepancy in the legacy effect of precipitation was detected between closed shrublands and open shrublands. The discrepancy in LUE m between closed shrublands and open shrublands has also been reported in existing research ).

Limitations
This study has some limitations. It aimed to analyze the multi-spatiotemporal heterogeneity of the legacy effect of precipitation and temperature on vegetation growth. However, the response of vegetation to other factors was not considered. For instance, previous studies have reported that biogeochemical drivers (e.g. CO 2 fertilization effects, nitrogen deposition) (Ukkola et al. 2015;Zhu et al. 2016), human activities (e.g. irrigation, pesticides, and fertilizers) and natural disturbances (e.g. solar activity, wildfire, and climate oscillations) (Li, Fan, and Xu 2015;Verón et al. 2012;Zhao et al. 2020) also play a vital role in vegetation growth. Meanwhile, using the vegetation index based on remote sensing to detect more complicated vegetation-climate interactions, such as the legacy effect of photosynthesis and photorespiration on shorter scales, is a daunting task. In addition, the response of vegetation growth to the cumulative effects of meteorological factors was not revealed in the present study. Nevertheless, our findings demonstrate that vegetation growth may be more sensitive to climate at important time nodes. Thus, we speculate that the cumulative effect of meteorological parameters may weaken the significance of the legacy effect. In contrast, these findings also indicate that the dynamic process of the legacy effect on a monthly scale can well describe its variation mechanism.

Conclusions
In summary, investigating the legacy effect is of great significance for exploring the mechanism of climatevegetation interactions, understanding the balance of terrestrial ecosystems, and formulating strategies for food security and environmental management. In this study, we determined the multi-spatiotemporal heterogeneity of the legacy effect of climate on vegetation growth, as well as the response of different vegetation types. Our results show that the heterogeneity of the legacy effects is a complicated process of dynamic variation, which can be summarized as the comprehensive characteristics of vegetation response to climate with regional discrepancy, ecosystem category, and multi-temporal scale. The multispatiotemporal heterogeneity of the legacy effect of climate on vegetation can explain the preference of vegetation growth to hydrothermal conditions. Moreover, different vegetation types exhibited heterogeneous responses of legacy effect to precipitation and temperature on multi-spatiotemporal scales, and the lag time in spring was shorter than that in summer and autumn. The average legacy effect for different vegetation types was approximately 1-2 months. Hence, our results provide insights into the interrelationships between vegetation and climate change, which advance the understanding of the preference of vegetation to hydrothermal conditions across biomes and ecosystems. Meanwhile, this study has established a paradigm framework for revealing the dynamic response of vegetation to other more complex factors in this warmer world. Furthermore, our results emphasize that multi-spatiotemporal legacy effects should be incorporated into the vegetationclimate interaction model and the formulation of macro-environmental management policies.

Highlights
• A paradigm framework summarizes the multispatiotemporal legacy effects • Cross-wavelet and wavelet coherence transform reveal the legacy effects • Legacy effects of climate vary with vegetation types and temporal scales • Preference of vegetation to hydrothermal conditions is explained • Legacy effects should be incorporated into the vegetationclimate model

Disclosure statement
No potential conflict of interest was reported by the author(s).

Funding
This study was supported by the National Natural Science Foundation of China (Grant No. 41971368).