Seasonal Divergent Tree Growth Trends and Growth Variability along Drought Gradient over Northeastern China

With the increasing temperature and intensified drought, global climate change has profound impacts on tree growth in temperate regions, which consequently regulates terrestrial-atmosphere biogeochemical processes and biophysical feedbacks. Thus, increasing numbers of studies have addressed the long-term annual trends in tree growth and their response to climate change at diverse spatial scales. However, the potential divergence in tree growth trends and growth variability (represented by coefficient of variance) in different seasons across large-scale climate gradients remains poorly understood. Here, we investigated the tree growth trends and growth variability in different seasons across diverse drought conditions in forested regions over northeastern China during the period 1982–2015, using both remote sensing observations and in situ tree-ring measurements. We found clear seasonal divergence in tree growth trends during 1982–2015, and the apparent increase was mainly observed in spring and autumn, attributed mainly to the increase in spring temperature and autumn solar radiation, respectively, but not in summer. The magnitudes of increasing trends in tree growth decrease with the increase of the multi-year average dryness index (MAI) in semi-arid areas (1.5 < MAI < 4.0) in all seasons. We further revealed that the interannual variability in tree growth was much larger in the semi-arid regions than in the humid and semi-humid regions in all seasons, and tree growth variability was significantly and negatively correlated with the variations in temperature and water deficit. Our findings improve our understanding of seasonal divergence in tree growth trends and provide new insights into spatial patterns in forest vulnerability in a warmer and drier climate.


Introduction
Forest ecosystems cover a large percentage (circa 26%) of Earth's land surface [1], and play a crucial role in regulating the global carbon cycle, and biophysical land-atmosphere feedbacks [2][3][4][5]. Along with global warming, the spatiotemporal variations in climate conditions have changed dramatically, particularly in mid-latitudinal regions [6], which exerts remarkable effects on both the structure and functioning of forest ecosystems [7]. It is widely reported that temperate forests are increasingly susceptible to intensified drought stress, and global-change-type drought has already triggered widespread forest growth decline, and even forest mortality, across diverse bioclimatic regions [8][9][10][11]. However, on the other hand, forest ecosystems tend to be resilient to and to some extent have already adapted to the changing climate [12]. Trees may adapt to the increasing drought by adjusting their structural and functional attributes, and physiological processes (e.g., changing water use strategies and regulating the distribution of non-structural carbon) [13,14]. Nonetheless, the general picture is still lacking regarding the spatial pattern of tree growth variability at the regional level during past decades. A comprehensive understanding of this feature could provide crucial insight into the vulnerability and stability of temperate forests in response to the changing climate.
Seasonal variations in vegetation growth over temperate regions are exemplified as a direct manifestation of the short-term response/adaptation of vegetation to climate change, and a few studies have revealed marked differences in vegetation growth in different seasons [15,16]. Changes in vegetation growth in spring and autumn were even more pronounced and were attributed to interactive shifts in plant phenology (i.e., the timing of the beginning or end of the vegetation growing season: BOS/EOS, respectively) and asymmetric seasonal climate change [17,18]. The seasonal asymmetry of climate change is characterized by a higher warming rate in spring and winter than in other seasons [19]. And the increasing temperature in winter and spring may accelerate the melting of snow and frozen soil in many parts of temperate regions, which could benefit vegetation growth in spring [20]. The delay of EOS could prolong the forest growing season, and this may increase the growth of vegetation in autumn and alleviate strong drought stress in the peak growing season in water-limited regions [17]. Nevertheless, the spatial pattern of seasonal tree growth across large-scale climate gradients during past decades remains inconclusive.
Northeastern China is the main distribution region of forests over northern China, and is one of the central parts of the Three North Shelterbelt Development Program (http://tnsf.forestry.gov.cn/) across China. Changes in forest growth have crucial impacts on ecosystem services and regional terrestrial-atmosphere carbon and water exchanges. Therefore, increasing efforts have been devoted to the understanding of forest growth and its response to the changing climate [21][22][23]. Studies based on both satellite observations and in situ forest investigation have revealed that the temperate forest is particularly sensitive to climate change in northeastern China [17,22,24,25]. Past studies mainly focused on the annual tree growth trend and its drivers [8,[26][27][28][29]. However, the seasonal patterns in, and underlying drivers for, the tree growth trend and growth variability across large-scale climate gradients remain poorly understood. In this paper, we attempt to address these issues using a satellite-derived vegetation index, regional tree ring chronologies, and gridded climate data. Specifically, we try to: (1) explore the potential divergence in tree growth trends and their response to climate in different seasons; and (2) investigate the spatial pattern of tree growth variability and its possible drivers across large-scale climate gradients.

Study Region
Our study region (33 • 15 N~53 • 32 N, 102 • 44 E~135 • 1 E) is located in the eastern part of northern China. This region is the main area of forest distribution over northern China (Figure 1). Our study region is characterized by a continental monsoon climate, with~60.2% of annual precipitation occurring in summer (June-August) and regular snow in winter. There are distinct seasons with marked differences in hydrothermal conditions. The mean annual precipitation (MAP) varies from 282 mm to 1150 mm and gradually decreases from the southeast to the northwest. The mean annual temperature (MAT) shows great spatial variations and ranges from −9.0 • C to 1.8 • C. Our study region is dominated by semi-humid and semi-arid regions, and can be divided into seven climatic regions, including: the cold temperate humid zone, middle temperate humid zone, middle temperate semi-humid zone, middle temperate semi-arid zone, warm temperate humid zone, warm temperate semi-humid zone and warm temperate semi-arid zone. Forest types distributed in this region mainly include evergreen coniferous forest, evergreen broad-leaved forest, deciduous coniferous forest, deciduous broad-leaved forest and mixed forest. Forests in semi-arid and semi-humid regions over northern China are generally characterized by patchy distribution, with smaller sized patches in drier regions ( Figure 1).

Tree Growth Proxies and Land Cover Data
Interannual variations in tree growth were evaluated using two independent datasets: the latest Normalized Difference Vegetation Index (NDVI) derived from advanced very high resolution radiometer observations, produced by the Global Inventory Modeling and Mapping Studies (GIMMS) group (i.e., GIMMS NDVI3g, https://ecocast.arc.nasa.gov/data/pub/gimms/3g.v1/); and site-level standard tree-ring index (TRI) chronologies over northeastern China. The GIMMS NDVI3g dataset has a spatial resolution of 8 km and biweekly temporal resolution, which covers 1982 to 2015. This dataset has been comprehensively corrected for atmospheric effects, sensor degradation, solar zenith angle, cloud contamination, and intersensor difference [30]. The GIMMS NDVI3g dataset has been widely used to investigate and evaluate spatio-temporal patterns in vegetation growth/productivity across diverse bioclimatic regions [30][31][32][33]. In this paper, we limited our analyses to forested regions, and resampled the NDVI of these regions into a spatial resolution of 0.1° to match the gridded precipitation and temperature data. Regions with multi-year average NDVI < 0.1 during the growing season (April-October) were further excluded from our final analyses. Eight stand-level standard TRI chronologies across different hydrothermal conditions were used to verify the interannual variation in tree growth (Table 1). We sampled two Pinus tabuliformis Carr. forests and six Larix gmelinii (Rupr.) Kuzen. forests along a large climate gradient. In each forest stand, we established a plot of 25 m × 25 m avoiding wet microhabitat or rock outcrops. We measured the tree height and diameter at breast height (DBH) ( Table 1) and we cored all living trees with DBH > 5 cm [34]. At least two cores (through the pith of the tree or as close as possible to the pith) per tree were collected. All cores were dried, mounted and sanded with successively finer grades of sandpaper until annual rings could be distinguished easily.

Tree Growth Proxies and Land Cover Data
Interannual variations in tree growth were evaluated using two independent datasets: the latest Normalized Difference Vegetation Index (NDVI) derived from advanced very high resolution radiometer observations, produced by the Global Inventory Modeling and Mapping Studies (GIMMS) group (i.e., GIMMS NDVI3g, https://ecocast.arc.nasa.gov/data/pub/gimms/3g.v1/); and site-level standard tree-ring index (TRI) chronologies over northeastern China. The GIMMS NDVI3g dataset has a spatial resolution of 8 km and biweekly temporal resolution, which covers 1982 to 2015. This dataset has been comprehensively corrected for atmospheric effects, sensor degradation, solar zenith angle, cloud contamination, and intersensor difference [30]. The GIMMS NDVI3g dataset has been widely used to investigate and evaluate spatio-temporal patterns in vegetation growth/productivity across diverse bioclimatic regions [30][31][32][33]. In this paper, we limited our analyses to forested regions, and resampled the NDVI of these regions into a spatial resolution of 0.1 • to match the gridded precipitation and temperature data. Regions with multi-year average NDVI < 0.1 during the growing season (April-October) were further excluded from our final analyses. Eight stand-level standard TRI chronologies across different hydrothermal conditions were used to verify the interannual variation in tree growth (Table 1). We sampled two Pinus tabuliformis Carr. forests and six Larix gmelinii (Rupr.) Kuzen. forests along a large climate gradient. In each forest stand, we established a plot of 25 m × 25 m avoiding wet microhabitat or rock outcrops. We measured the tree height and diameter at breast height (DBH) ( Table 1) and we cored all living trees with DBH > 5 cm [34]. At least two cores (through the pith of the tree or as close as possible to the pith) per tree were collected. All cores were dried, mounted and sanded with successively finer grades of sandpaper until annual rings could be distinguished easily.  The distribution of forests over our study region was identified from the Moderate Resolution Imaging Spectroradiometer dataset (MCD12Q1). This dataset has a spatial resolution of 500 m and a yearly temporal resolution, which provides gridded land use data from five different global vegetation classification systems (https://e4ftl01.cr.usgs.gov/MOTA/MCD12Q1.051/). We selected the MCD12Q1 map in 2013 to identify forests, and the land use types were defined by the International Geosphere-Biosphere Programme (IGBP) (Table S1). We categorized evergreen coniferous forest, evergreen broad-leaved forest, deciduous coniferous forest, deciduous broad-leaved forest and mixed forest as forests. We resampled the land cover into a spatial resolution of 0.1 • by a moving window of 24 × 24 (pixels), and forest was consequently defined when the forest pixels with a spatial resolution of 500 m accounted for more than 60% in the moving window.

Climate Data
Gridded climate data used in this paper, including precipitation (PRE), temperature (TEM), and solar radiation (SR), were obtained from the China Regional High-temporal-resolution Surface Meteorological Elements Driven Data Set (CMFD) (http://westdc.westgis.ac.cn/data/7a35329c-c53f-4267-aa07-e0037d913a21). The CMFD dataset includes seven climate variables, with a spatial resolution of 0.1 • and temporal resolution of 3-h, daily, and monthly. The monthly precipitation, temperature and net solar radiation data for the period 1982-2015 were used in this paper. We calculated the total precipitation, mean temperature and mean net solar radiation in each season for the subsequent seasonal analyses.
The multi-year average dryness index (MAI) and the water deficit index (WD) were used to evaluate the drought conditions. MAI was used for identifying the regional mean drought condition, which averaged the annual dryness index from 1982 to 2015 following Equation (1), and the annual dryness index was calculated as the ratio of annual potential evapotranspiration (PET) to total annual precipitation following Equation (2). We approximately divided the study region into three drought gradients, with MAI < 1.0 standing for humid areas, 1.0 < MAI < 1.5 for semi-humid areas, and 1.5 < MAI < 4 for semi-arid areas [35] (Figure 1).
where, n is the number of years from 1982 to 2015, and AI i is the aridity index in ith year. The PET was calculated with the adjusted Thornthwaite equation [36,37], which has been widely used in large-scale research, where detailed parameters for the Penman-Monteith equation were sometimes hard to obtain [37,38]. The annual total PET was calculated by summing monthly PET (PETm), which was calculated using Equation (3): where Ta is the monthly mean temperature ( • C). We assigned Ta as zero when Ta was below zero [37,39]. N is the number of days of the corresponding month. L is the mean daily hours of sunshine. I is the caloric index and α is a parameter, and they are calculated using Equation (4) and Equation (5), respectively.
The PET showed a significant decrease from the southeast to the northwest, with the highest values of PET ranging from 1200 mm to 1600 mm in Huang-Huai-Hai Plain ( Figure S1) [40], and the lowest values of PET ranging from 400 mm to 600 mm in the northern part of our study region ( Figure S1) [41].
WD was calculated by total precipitation minus total potential evapotranspiration in a period of time following Equation (6) [42]: where t is the length of time. Thus, a more negative WD indicates more severe water deficit.

TRI Chronology
Ring widths were measured to the nearest 0.001 mm using a LINTAB system (version 6, Frank Rinn Company, Heidelberg, Germany). All measured ring-width series were plotted and cross-dated by comparing the patterns of wide and narrow rings to identify possible false rings, missing rings or measurement mistakes. The quality of visual cross-dating was further examined using the software COFECHA [43]. A few cores were discarded because of missing rings or rot in the pith. Although the number of trees varied widely across the eight forest stands (Table 1), at least 40 cores were used to construct the TRI chronologies for each of the eight forest stands. The standard TRI chronologies were generated using the dplR package in R [44]. Negative exponential curve fitting or the linear regression method were used to remove the long-term tree growth trend, which is commonly linked to the species-specific age-dependent growth trend, stand disturbance history, and gradual effects of other external forces (e.g., CO 2 fertilization). The eight TRI chronologies were first divided into two groups based on the MAI for subsequent analyses regarding tree growth response to climate variations. The first group with MAI <1.5 had two TRI chronologies, and the other group with MAI >1.5 had six TRI chronologies.

Statistical Analysis
Principal component analysis (PCA), linear regression analysis, partial correlations analysis, the generalized linear model (GLM) and one-way ANOVA were used in this paper. First, PCA was applied to identify the main features among the eight TRI chronologies, and the first principal component (PCA1) was used to characterize the regional tree growth pattern. Then, we applied linear regression analysis to assess the linear trends in the seasonal TEM, WD, SR and tree growth proxies (NDVI and TRI) at both regional and pixel levels, and the magnitudes of linear trends were quantified by the slopes of the linear regression. Third, partial correlation analysis was applied to explore the interannual response of tree growth to variations of TEM, WD and SR, and the GLM was further used to identify the relative contribution (CB) of multiple independent variables (TEM, WD and SR) to tree growth (NDVI) in each pixel. We averaged all CB in pixels for whole study region and the three categories of MAI [38,42]. We calculated the coefficient of variance (CV) to analyze the spatiotemporal patterns of tree growth (NDVI) variability. Finally, one-way ANOVA was used to assess the differences in CV between the three regions (humid region, semi-humid region and semi-arid region). All variables were linearly detrended prior to partial correlation and GLM analyses. In our study, we consistently defined the growing season for tree growth as April to October, with spring (SP), summer (SU) and autumn (AU) correspondingly defined as April to May, June to August, and September to October, respectively [33].

Seasonal Trends in Tree Growth and Climate Response
We first investigated the linear trends in regional mean vegetation growth during the period 1982−2015 over northeastern China. We observed a marked increasing trend (p < 0.01) in regional mean  (Figure 2a), as well as in six (three of which show a significant increasing trend) of eight stand level TRI chronologies ( Figure 3). Among different seasons, the increasing trends in regional mean NDVI were more prominent in spring and autumn (Figure 2b,d), with significant increasing rates of 0.011 decadal −1 and 0.015 decadal −1 , respectively ( Table 2). In contrast, no significant increasing trend was shown in summer (Table 2; Figure 2c). In addition, large spatial variations were found in linear trends of seasonal mean NDVI in different seasons in the period 1982-2015, with prominent increasing trends in the northern part of our study region in spring and autumn ( Figure 4). There did not appear to be a clear pattern in linear trends of NDVI along the drought gradient in humid and semi-humid regions in summer and autumn, in contrast to a general increase in linear trends in spring. The linear trends in seasonal NDVI showed greater spatial variations in semi-arid regions (1.5 < MAI < 4). The positive linear trends approximately decreased, and even became negative, along with increasing MAI in spring and autumn (Figure 4e-h). We first investigated the linear trends in regional mean vegetation growth during the period 1982−2015 over northeastern China. We observed a marked increasing trend (p < 0.01) in regional mean growing season vegetation growth in the forested region of northeastern China from the perspective of NDVI (Figure 2a), as well as in six (three of which show a significant increasing trend) of eight stand level TRI chronologies ( Figure 3). Among different seasons, the increasing trends in regional mean NDVI were more prominent in spring and autumn (Figure 2b, d), with significant increasing rates of 0.011 decadal −1 and 0.015 decadal −1 , respectively ( Table 2). In contrast, no significant increasing trend was shown in summer (Table 2; Figure 2c). In addition, large spatial variations were found in linear trends of seasonal mean NDVI in different seasons in the period 1982-2015, with prominent increasing trends in the northern part of our study region in spring and autumn ( Figure 4). There did not appear to be a clear pattern in linear trends of NDVI along the drought gradient in humid and semi-humid regions in summer and autumn, in contrast to a general increase in linear trends in spring. The linear trends in seasonal NDVI showed greater spatial variations in semi-arid regions (1.5 < MAI < 4). The positive linear trends approximately decreased, and even became negative, along with increasing MAI in spring and autumn (Figure 4e-h).  (a-d) is for growing season, spring, summer and autumn, respectively. Shaded areas represent the 95% confidence intervals of the linear trends, and *, ** and *** represent the significance of linear fitting, which indicate statistical significance at p < 0.1, p < 0.05 and p < 0.01, respectively, with Student's t-test.        During past decades, the air temperature, moisture conditions (here represented by WD), and solar radiation have undergone tremendous changes in our study region, with great differences in temporal trends (Figures S2 and S3; Table 2). In spring, temperature exerted significantly positive effects on NDVI in humid (MAI < 1.0) and semi-humid (1.0 < MAI < 1.5) regions, with partial correlation coefficients (PR) of 0.66 (p < 0.05) and 0.63 (p < 0.05) (Figure 5), and contribution rates (CB) of 38.26% and 33.95%, respectively ( Figure 6); whereas, NDVI exhibited significantly positive partial correlations with WD and SR in semi-arid regions in spring, with a PR of 0.46 (p < 0.05) and 0.38 (p < 0.05), respectively (Figure 5b). In summer, NDVI was significantly and positively correlated with SR in both humid and semi-humid regions, whereas temperature was negatively correlated with tree growth in semi-arid regions, although not at a significant level. In autumn, NDVI was significantly and positively correlated with SR in semi-humid and semi-arid regions, with CB ranging from 18% to 22% (Figures 5 and 6; Figures S2 and S3). In addition, water deficit acted as the main limiting factor for annual tree radial growth from the perspective of TRI chronologies in both semi-humid regions and semi-arid regions (Figure 7). Compared with the semi-humid region, tree growth showed a much stronger water limitation in semi-arid regions, with PR of 0.38 (p < 0.05) and 0.52 (p < 0.05), respectively (Figure 7). Note: GS, SP, SU, AU stand for growing season, spring, summer and autumn, respectively. And *, ** and *** represent statistical significance at p < 0.1, p < 0.05 and p < 0.01 levels, respectively, with Student's t-test.
During past decades, the air temperature, moisture conditions (here represented by WD), and solar radiation have undergone tremendous changes in our study region, with great differences in temporal trends (Figures S2 and S3; Table 2). In spring, temperature exerted significantly positive effects on NDVI in humid (MAI < 1.0) and semi-humid (1.0 < MAI < 1.5) regions, with partial correlation coefficients (PR) of 0.66 (p < 0.05) and 0.63 (p < 0.05) (Figure 5), and contribution rates (CB) of 38.26% and 33.95%, respectively ( Figure 6); whereas, NDVI exhibited significantly positive partial correlations with WD and SR in semi-arid regions in spring, with a PR of 0.46 (p < 0.05) and 0.38 (p < 0.05), respectively (Figure 5b). In summer, NDVI was significantly and positively correlated with SR in both humid and semi-humid regions, whereas temperature was negatively correlated with tree growth in semi-arid regions, although not at a significant level. In autumn, NDVI was significantly and positively correlated with SR in semi-humid and semi-arid regions, with CB ranging from 18% to 22% (Figures 5 and 6; Figures S2 and S3). In addition, water deficit acted as the main limiting factor for annual tree radial growth from the perspective of TRI chronologies in both semi-humid regions and semi-arid regions (Figure 7). Compared with the semi-humid region, tree growth showed a much stronger water limitation in semi-arid regions, with PR of 0.38 (p < 0.05) and 0.52 (p < 0.05), respectively (Figure 7).

The Spatiotemporal Pattern in Tree Growth Variability and the Potential Drivers
We investigated and compared the coefficient of variation (CV) of seasonal mean NDVI along with the drought gradient in the period 1982-2015 (Figure 8a-d). Significantly (p < 0.05) higher CV of mean growing season NDVI was found in the semi-arid regions than in the humid and semi-humid regions. Among the seasons, a higher CV of seasonal mean NDVI was observed in spring than in summer and autumn in all climate regions; whereas, in summer and autumn, a higher CV of NDVI . The values on the bars in (b) are the corresponding partial correlation coefficients, and *, ** and *** indicate statistical significance at p < 0.1, p < 0.05 and p < 0.01 levels, respectively, with Student's t-test.

The Spatiotemporal Pattern in Tree Growth Variability and the Potential Drivers
We investigated and compared the coefficient of variation (CV) of seasonal mean NDVI along with the drought gradient in the period 1982-2015 (Figure 8a-d). Significantly (p < 0.05) higher CV of mean growing season NDVI was found in the semi-arid regions than in the humid and semi-humid regions. Among the seasons, a higher CV of seasonal mean NDVI was observed in spring than in summer and autumn in all climate regions; whereas, in summer and autumn, a higher CV of NDVI mainly occurred in semi-arid regions (Figure 8a-d). In contrast, a higher CV in both TEM and WD was found in humid and semi-humid regions than in semi-arid regions in all seasons (Figure 8e-l). The CV of SR showed higher values in semi-arid regions than in other climate regions in autumn, whereas there was no marked difference in the CV of SR among climate regions in spring and summer (Figure 8m-p). mainly occurred in semi-arid regions (Figure 8a-d). In contrast, a higher CV in both TEM and WD was found in humid and semi-humid regions than in semi-arid regions in all seasons (Figure 8e-l). The CV of SR showed higher values in semi-arid regions than in other climate regions in autumn, whereas there was no marked difference in the CV of SR among climate regions in spring and summer (Figure 8m-p). Partial correlation analyses were further conducted to explore the potential drivers for variability of NDVI in forested regions in different seasons (Figure 9). The CV of seasonal mean NDVI showed consistently negative correlations with the CV of both TEM and WD in all seasons, especially in summer (p < 0.05) for TEM and in summer and autumn (p < 0.05) for WD. However, the relationship between the CV of NDVI and the CV of SR among different seasons was inconsistent, showing a significant positive correlation in autumn and a significant negative correlation in spring, but no significant correlation in summer ( Figure 9).  Partial correlation analyses were further conducted to explore the potential drivers for variability of NDVI in forested regions in different seasons (Figure 9). The CV of seasonal mean NDVI showed consistently negative correlations with the CV of both TEM and WD in all seasons, especially in summer (p < 0.05) for TEM and in summer and autumn (p < 0.05) for WD. However, the relationship between the CV of NDVI and the CV of SR among different seasons was inconsistent, showing a significant positive correlation in autumn and a significant negative correlation in spring, but no significant correlation in summer (Figure 9).

Divergent Tree Growth Trends in Different Seasons
The changing climate has profound effects on the vegetation growth [4,18,45]. It is vital to explore the spatiotemporal patterns in tree growth trends and the underlying potential drivers in different seasons. Here, we discovered significant divergence in tree growth trends among different seasons, with a marked increase in NDVI occurring mainly in spring and autumn, but not in summer, over temperate Northeastern China. Previous studies reported significant changes in the beginning and end of growing seasons for trees in most regions over the Northern Hemisphere in response to the changing climate [18]. Increased temperatures and shifting precipitation patterns mainly contributed to an earlier spring and delayed growing season, which extends the length of the growing season for tree growth and could increase tree growth in spring and autumn [7,17]. We consistently observed that increasing spring temperatures tend to enhance NDVI in spring. In contrast, we identified that variation in solar radiation act as the dominating factor in regulating tree growth in autumn in our study region [7,46]. Further, due to the spatial heterogeneity in hydrothermal conditions and the differences in tree species (e.g., coniferous forests vs. broad-leaf forest) in northeastern China, the tree growth trends are significantly different across diverse regions [47][48][49]. Our results show the linear trends of NDVI declined with the increase of MAI in semi-arid regions in spring and autumn but show a vague pattern in humid and semi-humid regions (Figure 4). In humid regions, water availability was likely not the main limiting factor for NDVI; instead, variations in temperature, solar radiation as well as CO2 concentration interactively controlled tree growth [46,50] (Figure 5; Figure S2). Tree growth in drier regions tends to be increasingly affected by early growing season water availability and winter snow exerts important effects on tree growth, with contributions ranging from 2% to 23% across diverse climate zones over northeastern China [25]. More importantly, the water-limited regions are becoming drier accompanied due to rapid climate warming [51], which exerts intensified drought stress on tree growth and could further trigger tree growth decline ( Figure 5).
Regarding the TRI chronologies, there were no consistently significant increasing trends among the eight TRI chronologies, but these trends exhibited a similar spatial pattern as the NDVI results (Figures 3 and 4). However, the small samples with only eight TRI chronologies may not fully capture Figure 9. The partial correlation coefficient (PR) between the CV of seasonal mean NDVI and the CV of temperature (TEM), water deficit index (WD) and net solar radiation (SR) in our study region. GS, SP, SU and AU are growing season, spring, summer and autumn, respectively.

Divergent Tree Growth Trends in Different Seasons
The changing climate has profound effects on the vegetation growth [4,18,45]. It is vital to explore the spatiotemporal patterns in tree growth trends and the underlying potential drivers in different seasons. Here, we discovered significant divergence in tree growth trends among different seasons, with a marked increase in NDVI occurring mainly in spring and autumn, but not in summer, over temperate Northeastern China. Previous studies reported significant changes in the beginning and end of growing seasons for trees in most regions over the Northern Hemisphere in response to the changing climate [18]. Increased temperatures and shifting precipitation patterns mainly contributed to an earlier spring and delayed growing season, which extends the length of the growing season for tree growth and could increase tree growth in spring and autumn [7,17]. We consistently observed that increasing spring temperatures tend to enhance NDVI in spring. In contrast, we identified that variation in solar radiation act as the dominating factor in regulating tree growth in autumn in our study region [7,46]. Further, due to the spatial heterogeneity in hydrothermal conditions and the differences in tree species (e.g., coniferous forests vs. broad-leaf forest) in northeastern China, the tree growth trends are significantly different across diverse regions [47][48][49]. Our results show the linear trends of NDVI declined with the increase of MAI in semi-arid regions in spring and autumn but show a vague pattern in humid and semi-humid regions (Figure 4). In humid regions, water availability was likely not the main limiting factor for NDVI; instead, variations in temperature, solar radiation as well as CO 2 concentration interactively controlled tree growth [46,50] (Figure 5; Figure S2). Tree growth in drier regions tends to be increasingly affected by early growing season water availability and winter snow exerts important effects on tree growth, with contributions ranging from 2% to 23% across diverse climate zones over northeastern China [25]. More importantly, the water-limited regions are becoming drier accompanied due to rapid climate warming [51], which exerts intensified drought stress on tree growth and could further trigger tree growth decline ( Figure 5).
Regarding the TRI chronologies, there were no consistently significant increasing trends among the eight TRI chronologies, but these trends exhibited a similar spatial pattern as the NDVI results (Figures 3 and 4). However, the small samples with only eight TRI chronologies may not fully capture the spatial pattern in tree radial growth over the regional level and more samples are needed in future studies. We found water deficit acts as the main limiting factor for annual tree radial growth, which is different from the results of NDVI. In fact, NDVI is an index of vegetation canopy cover or the greenness of vegetation, and TRI reflects the radial growth of trees [52,53]. Although some previous studies show tree radial growth and forest canopy activity may be coupled owing to their consistent sensitivity to seasonal temperature/water conditions in northern high latitudinal forests [54,55], the relationship between tree radial growth and forest canopy activity varies among regions and tree species in arid and semi-arid regions because of the potential different sensitivity to temperature and drought [56]. Thus, the different features of these two proxies for representing vegetation growth could also partly explain the difference in the response patterns [52].

Tree Growth Variability and Drought Adaptation of Temperate Forests
Higher variability of tree growth in response to changing climate might imply a higher climate vulnerability. Drought stress tends to be further intensified in temperate regions in a warmer climate. How forests cope with increasing drought stress will exert great effects on both the functioning and stability of forest ecosystems. The adaptation of trees to drought manifests in many ways. At the individual level, trees adapt to drought by adjusting their growth rates and carbon allocation strategy [57,58]. At the ecosystem level, forests in drier conditions may have stronger resilience in response to drought [12,58]. Further, the sensitivity of vegetation to drought varies with bioclimatic regions and forest structures [49,59,60]. Our results showed that the interannual variability (CV) in NDVI was much larger in the semi-arid regions than in humid and semi-humid regions, which suggests tree growth in the semi-arid regions might be more vulnerable to climate variations, particularly intensified drought stress [60]. Interestingly, we revealed a negative correlation between the interannual variability in NDVI and the variability in TEM and WD in all seasons, further validating the notion that forests in a more variable climate tends to be more resilient to climate variations [42,61], and owing partly to the ability of tree growth to recover rapidly after drought events, benefiting from their substantial xylem plasticity for hydraulic transport in dry conditions [12,62,63]. Besides this, the interannual variability of tree growth is a comprehensive result and is regulated by multitudinous factors (e.g., SR and vapor pressure deficit) [64][65][66][67]. In this study, we observed a close positive relationship between the interannual variability of tree growth and SR variations, and such a relationship is season-dependent ( Figure 9); however, the underlying mechanisms are still unclear and need more comprehensive understanding. In addition, different tree species have developed specific adaptive strategies in ecological function [58] and physiological structure [68] for acclimatization, and the potential divergent vulnerability/adaptation of tree growth to drought among different tree species could partly reshape the spatiotemporal pattern in tree growth variability [52,69]. Although we identified clear patterns in both tree growth variability among different seasons and climate response at the regional level, further efforts are urgently needed to better understand the contribution of species-specific drought response to the stability of regional forest ecosystems in a warmer and drier climate.

Conclusions
In summary, we discovered clearly marked divergence in seasonal trajectories in tree growth and asymmetric tree growth variability along with the drought gradient over temperate Northeastern China, using both remote sensing retrievals and tree ring measurements. We observed a much higher increasing rate of regional tree growth in spring and autumn by analyzing NDVI, attributed to increasing temperature and solar radiation, respectively. In contrast, we did not observe a significant increasing trend in regional NDVI in summer. Interestingly, we found much higher interannual variability in tree growth in the semi-arid regions than in humid and semi-humid regions using NDVI, and the variability of NDVI was negatively correlated with the variability in hydrothermal conditions. Our findings provide crucial insights into a comprehensive understanding of spatiotemporal pattern in forest growth and the vulnerability of forests to the changing climate.
Supplementary Materials: The following are available online at http://www.mdpi.com/1999-4907/10/1/39/s1, Figure S1: The spatial distribution (a) and interannual variability of regional mean PET (b) in northeastern China from 1982 to 2015, Figure S2: The linear trends of regional mean TEM (a-d), WD (e-h) and SR (i-l) in different seasons in northeastern China from 1982 to 2015. Shaded areas indicate the 95% confidence interval of the linear trends and k in y-axis is the slope of the linear fitting. The green points are the values of TEM, WD and SR, Figure S3: The distribution of linear trends of seasonal mean TEM (a-d), WD (e-h) and SR (i-l) along with MAI in different seasons over northeastern China. The red dashed lines is the line of MAI = 1.0 and MAI = 1.5, respectively, which means the thresholds for different drought regions, Table S1: The information of land cover classification in MCD12Q1.