Direct and Lagged E ﬀ ects of Spring Phenology on Net Primary Productivity in the Alpine Grasslands on the Tibetan Plateau

: As a key biotic factor, phenology exerts fundamental inﬂuences on ecosystem carbon sequestration. However, whether spring phenology a ﬀ ects the subsequent seasonal ecosystem productivity and the underlying resource limitation mechanism remains unclear for the alpine grasslands of the Tibetan Plateau (TP). In this study, we investigated the direct and lagged seasonal responses of net primary productivity (NPP) to the beginning of growing season (BGS) along a precipitation gradient by integrating ﬁeld observations, remote sensing monitoring and ecosystem model simulations. The results revealed distinct response patterns of seasonal NPP to BGS. Speciﬁcally, the BGS showed a signiﬁcant and negative correlation with spring NPP ( R = − 0.73, p < 0.01), as evidenced by the direct boosting e ﬀ ects of earlier BGS on spring NPP. Moreover, spring NPP was more responsive to BGS in areas with more annual precipitation. The boosting e ﬀ ects of earlier BGS on NPP tended to weaken in summer compared with that in spring. Sequentially, BGS exhibited stronger positive correlation with autumn NPP in areas with less annual precipitation, which suggested the enhanced lagged suppressing e ﬀ ects of earlier spring phenology on ecosystem carbon assimilation during the later growing season under aggravated water stress. Overall, the strengthened NPP in spring was o ﬀ set by its decrement in autumn, resulting in no obvious relationship between BGS and annual NPP ( R = − 0.34, p > 0.05) for the entire grasslands on the TP. The ﬁndings of this study imply that the lagged e ﬀ ects of phenology on the ecosystem productivity during the subsequent seasons should not be neglected in the future studies. ﬀ ects of earlier BGS on summer NPP was weaker than on spring NPP. Moreover, our study revealed the negative lagged e ﬀ ects of earlier BGS on autumn NPP, especially in dry parts of the TP. The counterbalance among each season resulted in neutral e ﬀ ects of earlier BGS on annual NPP. The lagged e ﬀ ects of spring phenology on the subsequent seasonal productivity may help provide suggestions on optimizing the timing of grazing in alpine grasslands under changing growing season conditions.


Introduction
Climate change affects ecosystem carbon balance directly through changing photosynthesis, and indirectly through adjusting biological processes [1]. Under certain circumstances, the indirect effects might contribute to a higher proportion of variation than the direct effects do [2,3]. As one of the key biogeochemical processes influenced by climate, phenology exerts tight controls on ecosystem carbon sequestration and productivity [4,5]. During recent decades, satellite and ground-based evidences show that phenology has experienced significant shifts in response to climate change [6][7][8][9], which unavoidably influences exchanges of carbon, water and energy between vegetation and the atmosphere [10,11].
As the 'third pole' of the planet, the Tibetan Plateau (TP) is the highest plateau in the world. It covers an area of approximately 2.5 million km 2 and has an average altitude higher than 4000 m [12]. The alpine grassland (mainly composed of alpine meadow and steppe) is the dominant vegetation type on the TP, which is extremely sensitive to climate change [13]. As a major constituent of the carbon pool in the alpine ecosystem, alpine grasslands play a crucial role in regulating the regional carbon balance [14]. In recent decades, the TP has experienced unprecedented climate warming [15], which has affected various aspects of the fragile ecosystems, including shifted phenology and fluctuated ecosystem productivity [16,17]. However, how phenology regulates ecosystem productivity is still a bottleneck for us to understand the responses of global carbon cycles to climate change.
In the TP, various land surface phenological metrics (e.g., the start, end and length of the growing season) have been investigated for their effects on ecosystem carbon assimilation [18][19][20][21]. However, these studies mainly concerned the relationship between phenological metrics and annual ecosystem productivity. How phenology affects ecosystem productivity of each individual season is still largely obscure for the TP. The lagged effects of spring phenology on ecosystem productivity during the subsequent summer and autumn in particular have not been adequately investigated. Earlier spring growth may cause mixed effects on the subsequent ecosystem productivity of each individual season in different environments [5,22]. Previous studies reported that warming-induced earlier onset of the growing season would result in a longer growing season in some forest ecosystems, which in turn leads to strengthened ecosystem productivity [23,24]. However, earlier spring phenology might inhibit carbon assimilation, due to aggravated water stress in summer or autumn for the water-limited grassland ecosystems [25,26]. Water availability is one of the critical environment factors that regulate the response of vegetation productivity to spring phenology. Previous studies had revealed earlier spring phenology could aggravate soil moisture deficits in later growth stages by increasing evapotranspiration [27,28], leading to suppressed vegetation activities. However, few studies have reported the linkage between water resource limitation and the lagged effects of spring phenology on the ecosystem productivity in the subsequent seasons on the TP. The TP is characterized with frigid and dry environments [13]. Spatially, the annual precipitation in this region decreases from southeast to northwest [29]. Under climate warming, start of growing season has advanced significantly in the past decades in the TP [16,17]. The combination of arid environments and earlier spring phenology provides a unique opportunity for investigating the direct and lagged effects of spring phenology on ecosystem productivity.
In this study, our objectives were to investigate how spring phenology and climatic variables co-regulate carbon assimilation in direct and lagged ways for the alpine grasslands of the TP during 2001-2015, based on remote sensing data. Utilizing regionally optimized ecosystem productivity model and in-situ validated land surface phenology, our specific objectives were (1) to identify how spring phenology influences seasonal net primary productivity (NPP) during spring and the subsequent summer and autumn, and (2) to reveal the water resource limitation mechanism underlying the linkage between spring phenology and ecosystem productivity by direct and lagged effects.

Datasets
The Normalized Difference Vegetation Index (NDVI) is commonly used for land surface phenology and ecosystem productivity estimation. In this study, we used the 16-day composite NDVI product (MOD13Q1, Collection 6, downloaded from: https://ladsweb.modaps.eosdis.nasa.gov) covering 2001-2015 at 250 m spatial resolution from the Moderate Resolution Imaging Spectroradiometer (MODIS) Terra satellite. Noises caused by cloud and snow contamination were further removed [29]. A time series of annual minimum and uncontaminated NDVI was generated for each pixel during 2001-2015. The median value in the minimum NDVI time series was then extracted as the background NDVI value, which was later used to replace the smaller values and snow flagged values [30]. Moreover, we removed the remaining cloud contamination in NDVI time series by using Savitzky-Golay filter [31]. The 16-day temporal resolution NDVI was used to extract beginning of growing season (BGS). In addition, the maximum value composite (MVC) method was used to construct the monthly NDVI data for the calculation of monthly NPP [32]. To avoid the uncertainty caused by land cover change, we only considered pixels where grassland cover type remained identical during 2001-2015, according to the criteria set in a previous study [33]: (1) the average NDVI for June-September should be greater than 0.1; (2) the annual maximum NDVI should exceed 0.15 and occur within July-September; (3) the average NDVI for July-September should be greater than 1.2 times the average NDVI for November-March; and, (4) the average NDVI in winter (December-February) should be lower than 0.4. In addition, the green-up (also termed as MCD12Q2-derived BGS) derived from the MODIS Land Cover Dynamics Collection 6 dataset (MCD12Q2) (https://lpdaac.usgs.gov/products/mcd12q2v006/) was compared in performance with BGS derived from MOD13Q1 NDVI in this study.
Field data during 2001-2012 collected by China Meteorological Administration (CMA) were used to validate satellite-derived BGS and model-simulated NPP (Figure 1) [17]. The field sites are composed of eight phenology observations and seven biomass observations, and they cover both alpine meadow and steppe. Phenology and biomass observations were carried out in fenced natural pastures with an area of 100 m × 100 m at each station. The species-specific phenology was observed every two days, in accordance with uniform observation criteria [34]. The green-up date is identified when 10% of individual grasses per species at each station display leaves and grow up to 1 cm in spring. The monitored species is listed in Table 1. The multispecies mean green-up date in each year at each station was derived to represent the phenological status of the local plant community [35]. According to the observation criteria [34], four random plots (1 m × 1 m) were selected to measure aboveground dry biomass in each month during the growth period. For grasslands, NPP can be estimated by peak biomass [36]. The belowground biomass was estimated according to the average ratio of belowground to aboveground net production (3.18) for alpine grasslands [29]. The total biomass [g m −2 ] was converted to NPP [g C m −2 ] with a factor of 0.45 [37].
Remote Sens. 2020, 12, x FOR PEER REVIEW 3 of 14 for each pixel during 2001-2015. The median value in the minimum NDVI time series was then extracted as the background NDVI value, which was later used to replace the smaller values and snow flagged values [30]. Moreover, we removed the remaining cloud contamination in NDVI time series by using Savitzky-Golay filter [31]. The 16-day temporal resolution NDVI was used to extract beginning of growing season (BGS). In addition, the maximum value composite (MVC) method was used to construct the monthly NDVI data for the calculation of monthly NPP [32]. To avoid the uncertainty caused by land cover change, we only considered pixels where grassland cover type remained identical during 2001-2015, according to the criteria set in a previous study [33]: (1) the average NDVI for June-September should be greater than 0.1; (2) the annual maximum NDVI should exceed 0.15 and occur within July-September; (3) the average NDVI for July-September should be greater than 1.2 times the average NDVI for November-March; and, (4) the average NDVI in winter (December-February) should be lower than 0.4. In addition, the green-up (also termed as MCD12Q2derived BGS) derived from the MODIS Land Cover Dynamics Collection 6 dataset (MCD12Q2) (https://lpdaac.usgs.gov/products/mcd12q2v006/) was compared in performance with BGS derived from MOD13Q1 NDVI in this study. Field data during 2001-2012 collected by China Meteorological Administration (CMA) were used to validate satellite-derived BGS and model-simulated NPP ( Figure 1) [17]. The field sites are composed of eight phenology observations and seven biomass observations, and they cover both alpine meadow and steppe. Phenology and biomass observations were carried out in fenced natural pastures with an area of 100 m × 100 m at each station. The species-specific phenology was observed every two days, in accordance with uniform observation criteria [34]. The green-up date is identified when 10% of individual grasses per species at each station display leaves and grow up to 1 cm in spring. The monitored species is listed in Table 1. The multispecies mean green-up date in each year at each station was derived to represent the phenological status of the local plant community [35]. According to the observation criteria [34], four random plots (1 m × 1 m) were selected to measure aboveground dry biomass in each month during the growth period. For grasslands, NPP can be estimated by peak biomass [36]. The belowground biomass was estimated according to the average ratio of belowground to aboveground net production (3.18) for alpine grasslands [29]. The total biomass [g m -2 ] was converted to NPP [g C m -2 ] with a factor of 0.45 [37].   The monthly meteorological data of 172 stations on the TP and surrounding areas during 2001-2015 were collected from the China Meteorological Data Service System (http://data.cma.cn/), including mean temperature, total precipitation, percentage of sunshine and total solar radiation. Since solar radiation data is only available in 26 stations, we used the Angstrom-Prescott model to calculate the monthly total solar radiation for all stations based on the percentage of sunshine [38,39]. Lastly, the meteorological data was interpolated into 250 m × 250 m raster grids using ANUSPLIN 4.3 software [40]. The digitalized China Vegetation Map with a scale of 1:1000000 was utilized to obtain the distribution of alpine meadow and steppe ( Figure 1) [41]. In this study, all data were re-projected to the Albers conic equal area projection.

Determination of Land Surface Phenology
The temporal variations in NDVI during the growth period were fitted by a piecewise logistic function as follow [42] y where t is the time (day of year, DOY); y(t) is the NDVI value at time t; a and b are fitting parameters; c + d is the maximum NDVI value; and d is the initial background NDVI value. The BSG is defined as the date when the curvature change rate reaches the first local maximum value.

Simulation of NPP
The Carnegie-Ames-Stanford Approach (CASA) model, a light use efficiency model, was adopted to simulate NPP [43,44]. In the CASA model, NPP is estimated from the absorbed photosynthetically active radiation (APAR) and the actual light use efficiency (ε) as follows: where x is the specific geographic coordinate and t is the time; SOL is the total solar radiation; FPAR is the fraction of photosynthetic active radiation (PAR) absorbed by vegetation canopy and is determined from NDVI; the coefficient of 0.5 is the fraction of active incoming solar radiation (wavelength range of 0.4-0.7 µm) utilized by vegetation; ε max refers to the maximum light use efficiency under ideal conditions; T ε1 and T ε2 are the temperature stress coefficients representing the restriction of low and high temperature on LUE; and W ε is the water stress coefficient. More detailed information on the parameters calculation can refer to our previous studies [29,45].
The ε max was set to be a constant value of 0.389 g C MJ −1 for all vegetation types in the original CASA model [44], but it actually differs with vegetation types [45]. Our prior study had optimized the parameter of ε max for the alpine grasslands (0.4812 g C MJ −1 ) using field measurements [29].

Statistical Analyses
To validate satellite-derived phenology and ecosystem productivity using field observations, a 3 × 3 window average encompassing each site was extracted in consideration of the geospatial inaccuracies of the observed sites. The performance of satellite-derived BGS and NPP were assessed by the coefficients of determination (R 2 ) and p-values. The mean error (ME), median absolute error (MAE) and root mean square error (RMSE) were adopted as indicators of model accuracy. The direct and lagged effects of spring phenology on ecosystem productivity were estimated utilizing the Pearson correlations between BGS and NPP during spring and subsequent seasons. The Pearson correlations between BGS and NPP were also generated at a pixel level to analyze the spatial variations in the response of NPP to BGS for spring (March-May), summer (June-August), autumn (September-November) and the entire year, respectively. To evaluate regulation effects of precipitation on the relationships between BGS and seasonal NPP, we also conducted the Pearson correlation analysis between BGS and NPP for different biomes and precipitation intervals on the TP. The annual total precipitation was divided into three classes: 0-400 mm (arid areas), 400-600 mm (semi-arid areas) and >600 mm (dry sub-humid or humid areas) [46]. To identify the predominant direction of the detected significant relationships, the ratios of area with negative correlations to that with positive correlations at both 1% and 5% significant levels were used [47]. The threshold values of > 2.0 (or < 0.5) indicates a strong asymmetry biasing towards significant negative (or positive) correlations [47].

Ground Validation of Satellite-Derived BGS and NPP
The satellite-derived BGS used in this study displayed close agreement with field observations, as evidenced by ME of 5.8, MAE of 7.0 and RMSE of 9.7 days (R 2 = 0.62, p < 0.001) (Figure 2a

Spatial Variations in Seasonal Responses of NPP to BGS
The responses of NPP to BGS exhibited high seasonal variations ( Figure 4, Table 2). Spring NPP showed predominant negative correlations with BGS, with 38.74% of the grasslands being significantly correlated at the P < 0.05 level and mostly distributed in the central and eastern parts of the TP (Figure 4a). Significant positive correlations (P < 0.05) between BGS and spring NPP were only observed for 0.57% of the grasslands. In contrast, summer NPP showed weaker negative responses to BGS but the negative relationship was still predominated (Table 2). Significant negative correlations (P < 0.05) accounted for 10.96% of the grasslands, mostly in the central and southwestern

Spatial Variations in Seasonal Responses of NPP to BGS
The responses of NPP to BGS exhibited high seasonal variations ( Figure 4, Table 2). Spring NPP showed predominant negative correlations with BGS, with 38.74% of the grasslands being significantly correlated at the P < 0.05 level and mostly distributed in the central and eastern parts of the TP (Figure 4a). Significant positive correlations (P < 0.05) between BGS and spring NPP were only observed for 0.57% of the grasslands. In contrast, summer NPP showed weaker negative responses to BGS but the negative relationship was still predominated (Table 2). Significant negative correlations (P < 0.05) accounted for 10.96% of the grasslands, mostly in the central and southwestern

Spatial Variations in Seasonal Responses of NPP to BGS
The responses of NPP to BGS exhibited high seasonal variations ( Figure 4, Table 2). Spring NPP showed predominant negative correlations with BGS, with 38.74% of the grasslands being significantly correlated at the p < 0.05 level and mostly distributed in the central and eastern parts of the TP (Figure 4a). Significant positive correlations (p < 0.05) between BGS and spring NPP were only Remote Sens. 2020, 12, 1223 7 of 14 observed for 0.57% of the grasslands. In contrast, summer NPP showed weaker negative responses to BGS but the negative relationship was still predominated (Table 2). Significant negative correlations (p < 0.05) accounted for 10.96% of the grasslands, mostly in the central and southwestern parts of the TP (Figure 4b). Significant positive correlations (p < 0.05) between summer NPP and BGS were mostly found in the eastern TP, covering 2.71% of the grasslands. In autumn, more widespread positive correlations between NPP and BGS were observed compared to spring and summer (Figure 4c). Autumn NPP exhibited predominant positive correlations with BGS (Table 2). Significant positive correlations (p < 0.05) were mainly observed in the southwestern part of the TP, covering 8.21% of the grasslands. Only in 1.27% of the grasslands were there statistically significant negative correlations (p < 0.05) between BGS and autumn NPP. In terms of the entire year, the negative correlation between BGS and annual NPP was found to be statistically significant (p < 0.05) in 9.28% of the grasslands, primary in the central and southwestern TP (Figure 4d). Meanwhile, significant positive correlations (p < 0.05) accounted for only 1.74% of the entire grasslands.
Remote Sens. 2020, 12, x FOR PEER REVIEW 7 of 14 parts of the TP (Figure 4b). Significant positive correlations (P < 0.05) between summer NPP and BGS were mostly found in the eastern TP, covering 2.71% of the grasslands. In autumn, more widespread positive correlations between NPP and BGS were observed compared to spring and summer ( Figure  4c). Autumn NPP exhibited predominant positive correlations with BGS (Table 2). Significant positive correlations (P < 0.05) were mainly observed in the southwestern part of the TP, covering 8.21% of the grasslands. Only in 1.27% of the grasslands were there statistically significant negative correlations (P < 0.05) between BGS and autumn NPP. In terms of the entire year, the negative correlation between BGS and annual NPP was found to be statistically significant (P < 0.05) in 9.28% of the grasslands, primary in the central and southwestern TP (Figure 4d). Meanwhile, significant positive correlations (P < 0.05) accounted for only 1.74% of the entire grasslands.

Response of NPP to BGS in Different Biomes
For the entire alpine grassland ecosystem on the TP, the relationships between BGS and NPP varied among seasons ( Figure 5). Overall, spring NPP was significantly and negatively correlated with BGS (R = −0.73, P < 0.05), suggesting the promoting effects of earlier BGS on spring NPP. In summer, a weaker negative correlation between BGS and NPP was observed (R = −0.24). It indicated that BGS exerted weaker impacts on summer NPP than on spring NPP. However, the positive, but non-significant, correlations (R = 0.27) between BGS and autumn NPP indicate the down-regulation effects of earlier BGS on autumn NPP.
The correlation analysis was further conducted to explore the relationships between BGS and NPP for the two major biomes (alpine meadow and steppe) separately. BGS exhibited predominant negative correlations with spring NPP and summer NPP, and a predominant positive correlation with autumn NPP in both alpine meadow and steppe (Table 2). Overall, the negative correlation magnitude between BGS and NPP weakened from spring to summer, and then the direction reversed in autumn for both alpine meadow and steppe ( Figure 5). A significant negative correlation between  Table 2. Areal percentage of significant correlations between beginning of growing season (BGS) and seasonal or annual net primary productivity (NPP), and the ratio of negative correlated area to positive correlated area in different grassland types. Ratio values in bold and in italics indicate ratio of negative correlated area to positive correlated area > 2.0 and < 0.5, respectively.

Response of NPP to BGS in Different Biomes
For the entire alpine grassland ecosystem on the TP, the relationships between BGS and NPP varied among seasons ( Figure 5). Overall, spring NPP was significantly and negatively correlated with BGS (R = −0.73, p < 0.05), suggesting the promoting effects of earlier BGS on spring NPP. In summer, a weaker negative correlation between BGS and NPP was observed (R = −0.24). It indicated that BGS exerted weaker impacts on summer NPP than on spring NPP. However, the positive, but non-significant, correlations (R = 0.27) between BGS and autumn NPP indicate the down-regulation effects of earlier BGS on autumn NPP.    The correlation analysis was further conducted to explore the relationships between BGS and NPP for the two major biomes (alpine meadow and steppe) separately. BGS exhibited predominant negative correlations with spring NPP and summer NPP, and a predominant positive correlation with autumn NPP in both alpine meadow and steppe (Table 2). Overall, the negative correlation magnitude between BGS and NPP weakened from spring to summer, and then the direction reversed in autumn for both alpine meadow and steppe ( Figure 5). A significant negative correlation between BGS and spring NPP Remote Sens. 2020, 12, 1223 9 of 14 was observed in the alpine meadow (R = −0.76, p < 0.05), while a relatively weaker negative correlation was found in the alpine steppe (R = −0.51). In summer, the negative correlation was weaker for the alpine meadow (R = −0.20) than for the alpine steppe (R = −0.46). The BGS also showed a weaker positive correlation with the autumn NPP for the alpine meadow (R = 0.21) than for the alpine steppe (R = 0.42).

Correlation between BGS and NPP in Different Precipitation Classes
We further explored the relationship between BGS and NPP in different moisture conditions. The relationship varied among different precipitation classes ( Figure 6, Table 3). In spring, BGS showed predominant negative correlation with NPP for all precipitation classes at both p < 0.01 and p < 0.05 levels ( Table 3). An overall significant negative correlation (p < 0.05) between BGS and spring NPP was observed for the precipitation classes of 400-600 and > 600 mm, and the correlation magnitude strengthened from low to high precipitation class (Figure 6a). In summer, BGS exhibited a predominantly negative correlation with NPP in precipitation classes below 600 mm (Table 3). However, the overall significant negative correlation (p < 0.05) was only observed when precipitation was below 400 mm ( Figure 6b). The relationship between BGS and summer NPP was extremely weak when precipitation was above 600 mm. In autumn, the positive relationship between BGS and NPP dominated in all precipitation classes at both p < 0.01 and p < 0.05 levels ( Table 3). Autumn NPP was overall positively correlated with BGS and was less responsive to BGS in more adequate precipitation classes (Figure 6c).
Remote Sens. 2020, 12, x FOR PEER REVIEW 9 of 14 was below 400 mm ( Figure 6b). The relationship between BGS and summer NPP was extremely weak when precipitation was above 600 mm. In autumn, the positive relationship between BGS and NPP dominated in all precipitation classes at both P < 0.01 and P < 0.05 levels ( Table 3). Autumn NPP was overall positively correlated with BGS and was less responsive to BGS in more adequate precipitation classes (Figure 6c). Figure 6. The correlation coefficients between beginning of growing season (BGS) and net primary productivity (NPP) in spring (a), summer (b) and autumn (c) for different mean annual precipitation classes. Circles filled with black indicate significant correlations at P < 0.05 level. Table 3. Areal percentage of significant correlations between beginning of growing season (BGS) and seasonal net primary productivity (NPP), and the ratio of negative correlated area to positive correlated area in different precipitation classes. Ratio values in bold and in italics indicate ratio of negative correlated area to positive correlated area > 2.0 and < 0.5, respectively.   Table 3. Areal percentage of significant correlations between beginning of growing season (BGS) and seasonal net primary productivity (NPP), and the ratio of negative correlated area to positive correlated area in different precipitation classes. Ratio values in bold and in italics indicate ratio of negative correlated area to positive correlated area > 2.0 and < 0.5, respectively.

Direct and Lagged Effects of Spring Phenology on NPP
To fully understand the interannual variability of carbon sequestration, the direct and lagged temporal effects of phenology on ecosystem productivity needs to be considered. Earlier spring phenology are expected to have direct effects on spring productivity, and lagged effects on the subsequent summer and autumn productivity. Our study revealed the different response patterns of each individual season NPP to spring phenology. According to the asymmetry analysis in correlation, BGS exhibited predominant negative correlations with spring NPP and summer NPP, and a predominant positive correlation with autumn NPP. At a regional scale, the boosting effects of BSG on spring NPP (i.e., increased NPP caused by earlier BGS), was obvious for the whole alpine grassland ecosystem (R = −0.73, p < 0.05). However, the boosting effects weakened dramatically from spring to summer, and the direction even reversed in autumn. It indicated that the boosting effects of earlier BGS on ecosystem NPP abated gradually in a year.
The boosting effect of earlier BGS on spring NPP was also observed in other ecosystems around the world [23,48]. Earlier, BGS could extend the period for photosynthetic CO 2 uptake in spring by terrestrial ecosystems. Moreover, greater display of leaf area could enhance light interception and photosynthetic potential of vegetation canopy [49,50]. In addition, earlier spring phenology related to warming might increase the foliar nitrogen concentrations by accelerating the rate of nitrogen mineralization, thereby further enhancing photosynthetic activity [23]. We also observed that spring NPP was less correlated with BGS in areas with less precipitation. Limited water provision in arid regions can constrain plant physiological activity and undermine the stimulating effects caused by earlier BGS [51,52].
Compared with spring NPP, summer NPP was overall less correlated with BGS, which is in accord with another study in temperate China [26]. This might be resulted from the declined water availability caused by earlier BGS [28]. However, the correlation between NPP and BGS weakened more from spring to summer for alpine meadow than for the alpine steppe, as well as for areas with relatively adequate precipitations than arid regions. The BGS occurred earlier in the alpine meadow than that in the alpine steppe [16]. Then the extra water consumption in spring due to earlier BGS would be larger for the alpine meadow than for the alpine steppe. This might cause the larger decline in correlation between NPP and BGS from spring to summer for the alpine meadow.
In autumn, NPP exhibited positive correlation with BGS. In particular, the correlation magnitude strengthened from high to low precipitation class. It indicates that earlier spring phenology could aggravate the suppression of carbon assimilation in autumn in drier areas. The negative lagged effect of earlier BGS on the autumn NPP was consistent with a field-measured result on the TP [53]. Previous studies also reported a similar linkage between spring phenology and ecosystem productivity in the subsequent seasons for subalpine forests and temperate grasslands [25,54]. Differently to those studies, the suppression effect on productivity occurred in autumn, rather than summer in our study. The negative lagged effect of spring phenology on the autumn productivity might be caused by several reasons. Earlier spring growth could deplete soil water resources by increasing transpiration [25,55,56], consequently suppressing the physiological activities of vegetation. In addition, a decline in water availability can decrease plant nutrient uptake by reducing nutrient supply through mineralization, nutrient diffusion and mass flow in the soil [57][58][59], which also results in declined NPP.
Previous studies reported the boosting effects of earlier spring phenology on annual carbon assimilation [23,60]. However, for the unique alpine grasslands of the TP, the increased NPP in spring was offset by decreased NPP in autumn on the whole, resulting in no obvious relationship between BGS and annual NPP. It indicated that earlier BSG might not result in increased annual NPP, which corroborates the finding of previous study [20]. Other studies also confirmed that BGS was not the dominant contributor to annual carbon assimilation in the TP [18,19]. However, the role of spring phenology in regulating seasonal carbon dynamics through direct and lagged ways should be paid mounting attention. For example, warming-induced earlier BGS might lead to declined grassland productivity in autumn, which would limit the grassland carrying capacity. Then the livestock stocking rates need to decrease to match the grassland productivity and to maintain the human-nature balance. Therefore, improved knowledge on the relationship between BGS and ecosystem NPP is beneficial for implementing a plethora of sustainable development policies, such as grazing management and vegetation restorations.

Uncertainties
In this study, we were aimed to explore the direct and lagged effects of spring phenology on ecosystem NPP by seasons and assess the related water resources availability caused by earlier start of growing season. Autumn phenology can affect autumn and annual NPP to a certain extent. Due to its low accuracies as captured by remote sensing [33,61,62], and its minor effects on annual water budget, autumn phenology was not considered in this study. In the future, effects of autumn phenology entail to be considered in combination with spring phenology. Field monitoring is also critical. The field phenology and NPP observation sites in this study are mostly located in the central and eastern parts of the TP. Additional field sites need to be set up in the remote western plateau to provide in-situ validation on remote sensing products. In addition, anthropogenic activities, including grazing, can affect phenology and ecosystem productivity [20], which also need to be further investigated in the future. Besides, snowmelt can influence the effective precipitation for vegetation growth, which should be further considered to obtain a comprehensive understanding of the relationship between water availability and the response of vegetation productivity to phenology.

Conclusions
In this study, we revealed the direct and lagged effects of BGS on NPP for the alpine grassland on the TP, and how they were regulated by water resource limitation. The results showed direct positive effects of earlier BGS on spring NPP, and the positive effects were more pronounced in moisturized areas. Though earlier BGS exerted positive lagged effects on summer NPP, the boosting effects of earlier BGS on summer NPP was weaker than on spring NPP. Moreover, our study revealed the negative lagged effects of earlier BGS on autumn NPP, especially in dry parts of the TP. The counterbalance among each season resulted in neutral effects of earlier BGS on annual NPP. The lagged effects of spring phenology on the subsequent seasonal productivity may help provide suggestions on optimizing the timing of grazing in alpine grasslands under changing growing season conditions.