The QB, which runs through Central China between 102°54′~112°40′E and 30°50′~34°59′N, is an important climatic and geographical boundary between northern and southern China ( Figure 1 ). The entire region consists of three parts, namely, the Qinling Mountains, the Daba Mountains and the Jianghan Valley (Liu et al., 2016). The terrain clearly undulates, with an average annual precipitation of 700-1500 mm and an average annual temperature of 12-16 °C (Zhang and Liang, 2020). As the area is located in the transitional zone between the warm temperate zone and northern subtropical zone, vegetation types are clearly differentiated along EG, and their response to climate change is more sensitive (Deng et al., 2019; Qi et al., 2021).

The SOS in the QB was extracted from the MOD13Q1 EVI dataset provided by the official NASA website (https://ladsweb.modaps.eosdis.nasa.gov/search/) for the period 2001-2018 with a spatial resolution of 250 m and a temporal resolution of 16 days. Due to the influence of the sensors themselves and other external factors, there were missing data or outliers, resulting in differences between the EVI time series curves and the real vegetation growth patterns (Qi et al., 2021). Therefore, we used Savitzky-Golay (S-G) filtering based on Timesat 3.3 software to eliminate noise that deviated from the normal growth trend line. By comparing the results of multiple calculations, we set the window size to 5, the envelope iterations to 2, and the adaptation strength to 8 to reconstruct the EVI time series of the QB from 2001 to 2018. Then, we defined the SOS as the date when the EVI increased to 20% of the seasonal amplitude (Wu et al., 2021). To reduce the uncertainty from the outliers, we excluded pixels with SOS earlier than the 30th day of the year or later the 180th day of the year.

The Chinese meteorological forcing dataset (CMFD) is widely used for its continuous time coverage and consistent quality (http://data.tpdc.ac.cn/zh-hans/), with an accuracy between meteorological observation data and satellite remote sensing data, and a spatial resolution of 0.1° × 0.1° (Yang et al., 2010; He et al., 2020). Therefore, we selected temperature and precipitation datasets from the CMFD during 2000-2018 and analysed the temporal trends in winter temperature (WT, from the previous November to the beginning of greening in March), spring temperature (ST, from the beginning of greening in March to the end of greening in June), and spring precipitation (SP, from the beginning of greening in March to the end of greening in June) along EG in the QB at 100 m altitude intervals.

Considering the low spatial resolution of the CMFD dataset, the number of pixels in the QB is relatively small, and there may be errors based on pixel statistical analysis. To this end, we used a geographically weighted regression (GWR) model to downscale the WT, ST, and SP data to a spatial resolution of 0.01° (Lu et al., 2018; Lu et al., 2019). We selected the EVI and elevation as influencing factors and downscaled the WT, ST, and SP data by three levels with an intermediate resolution of 0.05° ( Figure S1 ) (Zhang et al., 2020; Zhang and Cheng, 2020; Xu and Cheng, 2021). The values were calculated as follows:

Where (u _i ,v _i ) . denotes the spatial coordinates at a spatial position i, x _i1 and x _i2 . present the values of elevation and EVI, respectively. β ₀, β ₁, and β ₂ . present the constant terms, the regression coefficients of elevation and the regression coefficient of EVI at raster i, espectively.

The digital elevation model (DEM) data were obtained from the 90 m resolution SRTM product provided by the Geospatial Data Cloud (https://www.gscloud.cn/). Land Use Data, derived from Globe Land Cover in 2010 (http://www.globallandcover.com/), were used to extract the distribution of the forest vegetation in the QB. The DEM and Land Use Data were resampled to 250 m to match the resolution of the EVI by using ArcMap 10.3.

The slope of SOS along EG for each year from 2001 to 2018 was used to compare the difference in SOS along EG among different years ( Table S1 ). The smaller the slope was, the smaller the difference in SOS along EG. The quadratic curve fitting method based on Origin 2018 was used to detect abrupt changes in the slope of SOS along EG. The year corresponding to the point where the slope line of the curve intersected with the 0 slope line was considered an abrupt change point. Then, Theil-Sen trend analysis was used to analyse the temporal trends of SOS, WT, ST, and SP in the QB from 2001 to 2018 and before and after the abrupt years, and the significance of trends was tested by the Mann-Kendall (M-K) statistical test by using MATLAB R2016a (Sen, 2012; Liu et al., 2016). This value was calculated as follows:

where β is the temporal trend of the SOS; j and i denote the time series; and x_j and x_i denote the SOS at times j and i, respectively. β > 0 indicates that the temporal trend of the SOS has a delayed trend, and β< 0 indicates an advanced trend.

Partial correlation analysis based on MATLAB R2016a was used to analyse the relationship between SOS and WT, ST, and SP. The degree of association between the two variables was measured by the partial correlation coefficient after excluding the effects of other control variables (Shen et al., 2020). The climate factors (WT, ST, and SP) were used as independent variables, and the SOS was the dependent variable. Statistical significance was determined at a level of P< 0.05 based on a two-tailed t-test. Additionally, the temporal trend and partial correlation coefficient patterns in SOS and the climate factors patterns along EG were analysed at 100 m altitude intervals. Considering the rarity of forest pixels at high elevations, regions with fewer than 100 pixels along EG were excluded, and only the regions below 3100 m were analysed in the QB.

The mean SOS from 2001 to 2018 was early in the valleys and late on the mountains. Among them, 89.71% of the mean SOS were concentrated around 73-105 days. Along EG, the mean SOS was significantly (P< 0.001) delayed by 1.7 days per 100 m increase ( Figure 2A ). The temporal trends of delayed SOS (48.96% of total pixels) were mainly distributed in the valley and at marginal low elevations, and the magnitude of delay was concentrated around 0-0.8 days·a^-1 ( Figure 2B ). Significant delayed (P< 0.05) were mainly distributed in the northeastern low elevation regions of the QB ( Figure S2A ). The advanced SOS (51.04%) was mainly distributed at high elevations, and the magnitude of advance was concentrated in 0-0.8 days·a^-1 ( Figure 2B ). Significant advanced (P< 0.05) were mainly distributed in the western, northern, and southern alpine regions ( Figure S2A ). Along EG, the temporal trends in delayed SOS gradually decreased with increasing elevation (below 1000 m), and the advanced SOS gradually increased (above 1000 m). These divergent trends in SOS between high and low elevations showed a more uniform trend of SOS along EG in the QB.

The change in slope of SOS along EG from 2001 to 2018 decreased by 0.26 ± 0.01 days 100 m^-1 per decade (P = 0.06), the differences in SOS along EG decreased continuously, and the SOS showed a more uniform trend along EG ( Figure 3B ). Furthermore, 2011 was identified as a year of abrupt change according to the change in the SOS slope along EG from 2001 to 2018 based on the quadratic curve fitting method ( Figure 3A ). The change in the SOS slope along EG significantly (P = 0.001) decreased by 0.85 ± 0.02 days 100 m^-1 per decade between 2001 and 2011. Since 2011, the change was reversed by 0.54 ± 0.02 days 100 m^-1 per decade.

The temporal trends in SOS before 2011 mainly showed a delayed trend (67.54%) that was distributed in the regions below 1500 m, and the magnitude decreased with increasing elevation ( Figures 4A–D ). Significant delayed (P< 0.05) were mainly distributed at low elevations in the eastern and southern regions ( Figure S2B ). While the advanced SOS (32.46%) was mainly distributed at high elevations above 1500 m, and the magnitude was increased ( Figures 4A–D ). Significant advanced (P< 0.05) were mainly distributed in the western, northern, and southern alpine regions ( Figure S2B ). These divergent trends showed a more uniform trend of SOS along EG in the QB.

However, the temporal trends in SOS after 2011 mainly showed an advanced trend (73.38%) at all elevations, and the magnitude decreased by 0.06 days·a ^-1 for each 100 m increase ( Figures 4B–D ). Significant advanced (P< 0.05) were mainly distributed in the eastern and southern regions ( Figure S2C ). The advanced SOS (53.97%), which changed from a delayed SOS before 2011, was mainly distributed at low elevations ( Figure S3 ), and the magnitude was stronger than that at high elevations. The difference in SOS gradually widened along EG. Shifts in the temporal trend of SOS at low elevations determined the direction of the uniform trend in SOS along EG.

The results of the trend analysis showed that the temporal trends in WT before 2011 mainly showed a cooling trend (79.82%), which was mainly distributed in marginal low elevation regions ( Figure 5A, C ). Significant decreased (P< 0.05) were mainly distributed in the eastern high elevation regions ( Figure S4A ). Along EG, the temporal trends in WT decreased by 0.02°C·(10a) ^-1 for each 100 m increase ( Figure 5D ). After 2011, the temporal trends in WT mainly showed a warming trend (87.47%), which was mainly distributed in the marginal low elevation areas ( Figures 5B, C ). Significant increased (P< 0.05) were mainly distributed in the eastern and southern high elevation regions ( Figure S4D ). Along EG, the temporal trends in WT increased by 0.05°C·(10a) ^-1 for each 100 m increase ( Figure 5D ).

The temporal trends in ST before 2011 mainly showed a cooling trend (60.26%), which was mainly distributed in high elevation regions ( Figures 6A, C ). Significant decreased (P< 0.05) were mainly distributed in the northern and eastern high elevation regions ( Figure S4B ). While the warming ST (39.74%) was concentrated in the marginal low elevation regions and the central Daba Mountains ( Figure 6A ). Along EG, the temporal trends in ST decreased by 0.02°C·(10a) ^-1 for each 100 m increase, and the temporal trends in ST showed a slight increase below 800 m and then cooled ( Figure 6D ). After 2011, the temporal trends in ST mainly showed a warming trend (73.85%), which was mainly distributed in the central and western regions. While the cooling ST (26.15%) was concentrated in the eastern low elevation regions ( Figures 6B, C ), with a significant decreased (P< 0.05) in the northeastern regions ( Figure S4E ). Along EG, the temporal trends in ST increased by 0.09°C·(10a) ^-1 for each 100 m increase, and the temporal trends in ST showed a decrease below 500 m and then warmed ( Figure 6D ).

The temporal trends in SP before 2011 increased on average by 0.45 mm·a^-1. Among them, the increased regions (59.4%) were mainly distributed at high elevations, while the decreased regions (40.6%) were distributed at marginally low elevations. The temporal trends in SP were of low amplitudes and mainly concentrated around 0-1 mm·a^-1 ( Figures 7A, C ). Along EG, the temporal trends in SP increased by 0.18 mm·(10a)^-1 for each 100 m increase ( Figure 7D ). After 2011, the temporal trends in SP increased by 4.41 mm·a^-1 on average, which was approximately 10 times higher than that before 2011. Among them, the increased regions (96.45%) were mainly concentrated around 2-8 mm·a^-1 ( Figures 7B, C ), with a significant increased (P< 0.05) in the eastern low elevation regions ( Figure S4F ). Along EG, the temporal trends in SP decreased by 0.29 mm·(10a)^-1 for each 100 m increase ( Figure 7D ).

Based on the results from the partial correlation analysis, we found that the SOS and WT were mainly negatively partially correlated (76.18%) in the QB before 2011, with a significant negative partial correlation (P< 0.05) in the central and eastern low elevation regions ( Figure S5A ), while the regions with a positive partial correlation (20.18%) were concentrated in the western and southern high elevation areas ( Figure 8A ). Along EG, the negative partial correlation coefficient gradually weakened with increasing elevation (below 2100 m), and the positive correlation coefficient continued to increase (above 2100 m) ( Figure 8D ). The SOS before 2011 showed a partial negative correlation with ST (75.36%), with a significant negative partial correlation (P< 0.05) in the central and eastern low elevation regions ( Figure S5B ), while the regions with a positive partial correlation coefficient (24.64%) were concentrated in the western and southern high elevation areas ( Figure 8B ). Along EG, the negative partial correlation coefficient significantly weakened (P< 0.001) with increasing elevation ( Figure 8E ). The SOS showed a partial negative correlation with SP (57.91%), with a significant negative partial correlation (P< 0.05) in the eastern low elevation regions and the western and southern high elevation regions ( Figure S5C ), while the regions with a positive partial correlation coefficient (42.09%) were distributed within all elevation gradients ( Figure 8C ). Along EG, the negative partial correlation coefficient did not change much (P = 0.28) with increasing elevation ( Figure 8F ).

For 2011-2018, the SOS and WT were mainly positively partially correlated (63.81%), with a significant positive partial correlation (P< 0.05) in the eastern low elevation regions and the central and western regions ( Figure S5D ), while the regions with a positive partial correlation (36.19%) were concentrated in the northeastern and southern high elevation areas ( Figure 9A ). Along EG, the positive partial correlation coefficient significantly increased (P< 0.001) with increasing elevation ( Figure 9D ). The SOS showed a partial negative correlation with ST (74%), with a significant negative partial correlation (P< 0.05) in the eastern low elevation regions and the central and western regions ( Figure S5E ), while the regions with a positive partial correlation coefficient (26%) were concentrated in the northeastern and southern high elevation areas ( Figure 9B ). Along EG, the negative partial correlation coefficient significantly increased (P< 0.001) with increasing elevation ( Figure 9E ). The SOS showed a partial negative correlation with SP (67.52%), with a significant negative partial correlation (P< 0.05) in the eastern low elevation regions and the western and southern high elevation regions ( Figure S5F ), while the regions with a positive partial correlation coefficient (32.48%) were sporadically distributed within all elevation gradients ( Figure 9C ). Along EG, the negative partial correlation coefficient significantly weakened (P = 0.03) with increasing elevation ( Figure 9F ).

In this study, we found a decreased WT and a slightly increased ST below 500 m in the QB from 2001 to 2011, which somewhat provides more sufficient accumulation of chilling and forcing for the onset of spring phenology (Fu et al., 2015a; Ettinger et al., 2020; Pan et al., 2022). However, the SOS in this region was delayed ( Figure 4D ). Temperature alone did not sufficiently explain the SOS patterns in this region. Moreover, the reduced SP in this region ( Figure 7D ) somewhat limited the utility of water and heat conditions and offset the advanced SOS that would have been caused by the decreased WT and increased ST (Jewaria et al., 2021; Zheng et al., 2021), resulting in a delay in SOS. Hence, the decreased SP was the predominant controlling factor for the delayed SOS in this region. Similarly, at elevations of 500-1500 m, the decreased WT and increased SP somewhat facilitated the accumulation of chilling and the water demand for vegetation growth (Shen et al., 2015; Ettinger et al., 2020), which would have advanced the SOS. However, this advance would have likely been counteracted by reduced forcing as a result of the reduced ST ( Figure 6D ), which better explained why the SOS in this region was delayed (Piao et al., 2019; Ganjurjav et al., 2020). The decreased ST may have a stronger impact on SOS relative to the decreased WT and the increased SP. Meanwhile, at elevations above 1500 m, the decreased ST reduced the accumulation of forcing and delayed the SOS. Moreover, the decreased WT and the increased SP in this region ( Figures 5D , 7D ) somewhat increased the accumulation of chilling and satisfied the water demand for the advanced SOS (Lin et al., 2022; Wang et al., 2022). The SOS response to reduced ST may have been less than that to the reduced WT and increased SP, resulting in the advanced SOS in this region. This showed that the joint control of the reduced WT and increased SP, rather than the ST, was the driving factor for the advanced SOS.

For 2011-2018, the increased WT and the reduced ST at elevations below 500 m were not conducive for more chilling and forcing accumulation (Yu et al., 2010; Wang et al., 2020a), which somewhat delayed the SOS. However, this delay was counteracted by the increased water demand as a result of the increased SP ( Figure 7D ) and led to the advance of SOS (Shen et al., 2015). The advanced SOS below 500 m was influenced by the increased SP rather than temperature. Similarly, at elevations above 500 m, the increased WT may reduce chilling accumulation and delay the SOS (Wang et al., 2020b; Jewaria et al., 2021). However, this delay was offset by the increased water demand and forcing accumulation ( Figures 6D , 7D ), resulting in the advanced SOS (Piao et al., 2019; Li et al., 2021a). The effects of increased ST and SP on SOS were stronger than those of increased WT. The present results suggest that the SP also plays a crucial role in controlling the patterns of SOS along EG, and analyses conducted from the perspective of temperature alone may not fully elucidate the intrinsic mechanism of these patterns in SOS.

Our research demonstrated that the magnitude of SOS response to SP gradually decreased and the response to temperature gradually increased with increasing elevation, which was consistent with previous studies in the QB (Liu et al., 2016; Chen et al., 2019). The SOS at high elevations, where the SP showed a consistent increase from 2001 to 2018 ( Figure 7D ), was less sensitive to SP. In other words, sufficient SP at high elevations resulted in less water stress on vegetation growth (Jennifer and Florence, 2018; Gupta et al., 2020), and larger amounts of SP would not advance SOS. Along with this was the potential for greater temperature sensitivity of SOS to maximise thermal benefits (Shen et al., 2015; Fu et al., 2021). This hypothesis was further supported by the weaker advance of SOS in the context of stronger precipitation after 2011 and the stronger partial correlation between SOS and temperature ( Figures 4D , 9 ). In contrast, at low elevations, where the SP showed a decrease from 2001 to 2011 ( Figure 7D ), maximised water usage led to the stronger sensitivity of SOS to SP (Shen et al., 2015; Chen et al., 2019). Stronger water stress led to a delayed SOS even under better temperature conditions ( Figures 5D , 6D ). In addition, we speculated that the SOS would advance if more rainfall occurred after stronger water stress, even if the temperatures were less optimal. This hypothesis was confirmed by the significant advance of SOS at low elevations under the significant increase in SP after 2011 ( Figure 4D , Figure 7D ). This result showed that the temporal trends in SP played a crucial role in controlling the patterns of SOS at low elevations, which determined the direction of the uniform trend in SOS along EG.

Recent research has predicted that by the end of this century, the climate in east-central China will have a continuous warm-dry trend (Ma et al., 2019; Zhang et al., 2022), which may enhance water evaporation at low elevations in the QB. In addition, the structural overshoot due to the warm winter may further increase the water stress during the growing season (Bastos et al., 2020; Zhang et al., 2021b), thus strengthening the uniform trend of SOS along EG in the QB. Therefore, the sensitivity of SOS to precipitation may be further enhanced, and future changes in the spatiotemporal distribution of precipitation, rather than temperature, are likely to have a stronger control on the direction of the uniform trend in SOS. Since a more uniform SOS along EG may compromise the stability and serviceability of mountain ecosystems (Inouye et al., 2000; Lenoir et al., 2008; Wolf et al., 2017), it is crucial to improve the water-use efficiency of vegetation at low elevations. Previous studies have shown that a higher species diversity could notably enhance drought resistance (Grossiord, 2020; Liu et al., 2022), and different tree species have different strategies and abilities to cope with water stress (Anderegg et al., 2018; Li et al., 2020b). In the new round of ecological restoration projects, conversion of the current monoculture to mixed-species tree plantations and the planting of resilient tree species with a high water-use efficiency could relieve potential water stress in the future.

In this study, we confirmed that the SOS extracted by the MOD13Q1 EVI can be used to accurately trace the spring phenology in the QB. This conclusion was predominantly based on the strong consistency with previous studies in related areas that were based on different data sources ( Table S2 ). In addition, our findings demonstrated that the SOS showed a more uniform trend along EG between 2001 and 2018 in the QB. These divergent trends of SOS between low and high elevations led to a uniform trend of SOS along EG in the QB, which was consistent with the patterns of SOS in the Alps (Vitasse et al., 2018). However, the driving mechanisms of the more uniform SOS in the QB were different from those reported in the Alps. Vitasse et al. (2018), based on the hypothesis that temperature plays the dominant role, found that the reduced chilling accumulation at low elevations caused by the warming WT moderated the magnitude of SOS advance compared to high elevations. The SOS was more advanced at higher altitudes than at lower altitudes, resulting in a more uniform trend of SOS along EG. In this study, we added the effects of SP and found that the SP, rather than temperature, determined the direction of the uniform trend in SOS by controlling the SOS patterns at low elevations. These differences may have been due to potential mechanisms playing diverse roles across different areas.

In addition, the results demonstrated negative relationships between ST, SP and SOS during 2001-2018 and a positive relationship between WT and SOS between 2011 and 2018, which was consistent with most previous phenological studies (Piao et al., 2006; Zheng et al., 2021). However, the WT was negatively correlated with SOS below 2100 m during 2001-2011 and the cooling WT delayed the SOS in the QB. The simultaneous reduction of ST between 500 and 2100 m and the SP below 500 m may have counteracted the advance of SOS caused by decreased WT (Cong et al., 2017; Asse et al., 2018), possibly explaining why the SOS was delayed under a cooling WT.

In contrast to our results, Gao et al. (2019) found that there were no prevalent trends of elevational homogenization of SOS in most regions worldwide over the last 30 years. This discrepancy was likely due to the difference in the main drivers of vegetation growth at different spatial scales. For example, the stronger water stress in arid areas results in a greater sensitivity of SOS to precipitation, but the sensitivity of SOS to temperature and precipitation may change at a larger scale (Cleverly et al., 2016; Cong et al., 2021). In addition, the magnitudes and even the directions of the response to climate change are largely different among vegetation types (Lesica and Kittelson, 2010; Fu et al., 2015b). The proportion of vegetation types at different spatial scales varies, which may lead to differences in the variation patterns of SOS.

First, using different abrupt change tests may lead to different results. We used the M-K test, sliding t-test, and quadratic curve fitting test to detect abrupt changes in the slope of SOS along EG from 2001 to 2018 ( Figure S6 ). We found that the abrupt change years were approximately 2003, 2008, and 2011 in the M-K test, sliding t test, and quadratic curve fitting test, respectively. Moreover, previous research has shown that the vegetation coverage reversed in 2010 in the QB during 2001-2014 (Liu et al., 2016). Chen et al. (2021) found that 2010 was the significant acceleration point for gross and net primary production variations in China during 2001-2018. Therefore, using the quadratic curve fitting test, we defined 2011 as an abrupt change year. Second, since the CMFD dataset was reanalysis data, the temperature and precipitation data already contained elevation information (He et al., 2020). Downscaling by selecting the EVI and elevation as influencing factors may overestimate the effect of elevation on temperature and precipitation, which somewhat impacts the downscaling results. Finally, different vegetation types have different responses to climate change. As the dividing line between north subtropical and temperate forests, vegetation types in the QB are diverse and have an obvious altitudinal differentiation (Deng et al., 2019; Qi et al., 2021). Therefore, the patterns of SOS along EG on vegetation types and the north and south slopes need to be further studied.