Rainy season onset mainly drives the spatiotemporal variability of spring vegetation green-up across alpine dry ecosystems on the Tibetan Plateau

It is still debatable whether temperature or precipitation mainly triggers spring vegetation green-up (SVG) in alpine dry ecosystems on the Tibetan Plateau. As phenological sensitivity to the arrival of monsoon-season rainfall would allow plants to simultaneously avoid drought and frost damages in the early growing season, we hypothesize that rainy season onset (RSO) rather than temperature mainly drives the spatiotemporal variability of SVG across alpine dry ecosystems over the Tibetan Plateau. Dates of RSO and SVG across 67 target areas nearby 67 weather stations over the Tibetan Plateau were calculated from time-series data of daily mean temperature and precipitation (1974–2013) and of the Normalized Difference Vegetation Index from the Moderate Resolution Imaging Spectroradiometer (2001–2013), respectively. Satellite-derived SVG was validated by 7-year observations (2007–2013) for leaf emergence of dominant species in alpine meadows along elevations (4400–5200 m) in Damxung of Tibet. We found that SVG generally synchronized with or was somewhat later than RSO although seasonal air temperatures were already continuously above 0 °C in 1 month before SVG dates. In pooled data across sites and years, the analysis of linear mixed model indicated that RSO (F = 42.109) and its interactions with pre-SVG precipitation (F = 6.767) and temperature (F = 4.449) mainly explained the spatio-temporal variability of SVG, while pre-SVG temperature and its interaction with precipitation did not have significant effects on SVG. Our data supported the hypothesis, suggesting that synchronization of SVG and RSO is a general spring phenological strategy across alpine dry ecosystems under influence of monsoon climate.


Scientific Reports
| (2020) 10:18797 | https://doi.org/10.1038/s41598-020-75991-w www.nature.com/scientificreports/ rainfall 22 , and desert plants in Africa green up ahead of rainy season 23 . Because phenological responses to climate change may be involved in complicated evolutionary adaptations to environmental stresses 3,5,18,19,22,23 , it is challenging to understand climatic controls on spring vegetation phenology in dry regions, especially at high elevations where variable precipitation regimes and the interaction of temperature and precipitation may further confound the issue. In alpine dry ecosystems on the Tibetan Plateau, it is still debatable whether temperature or precipitation mainly drives the spatiotemporal variability of spring phenology though numerous studies have endeavored to search the possible cues 14,[20][21][22] . Plant phenology is a key component of ecological fitness, serving as an adaptive trait in shaping species distribution and ecosystem function 24 . On the Tibetan Plateau, the typical Indian monsoon circulation is established during late May and June when the plateau air temperature is generally above 0 °C, and the interaction of Indian monsoon and the westerlies generally creates a cold and dry climate in the pre-monsoon season 25 . The vast alpine meadows and steppes on the Tibetan Plateau are the zonal climax vegetation types adapted to the cold and dry climate, and their fan-shaped distribution pattern is mainly controlled by the warm and moist air of Indian summer monsoon through the atmospheric water transports along the high mountain valleys in eastern Himalayas 26 . The monsoon-westerly interaction may have selected genotypes that respond better to the pre-monsoon cold and dry climate. Our observed data from the Kobresia meadow in Damxung of Tibet suggest that phenological sensitivity to the arrival of monsoon-season rainfall would allow plants to simultaneously avoid drought and frost damages in the early growing season 22 . However, little attention has been paid to the effect of monsoon rainfall onset. Because of lack of long-term observed data in remote high mountains, few studies have investigated whether such a spring phenological strategy found in the Kobresia meadow may generally occur in different alpine vegetation types (alpine meadows, steppes, dry shrubs, etc.) over the Tibetan Plateau. To further test the generality of our previous findings in the Kobresia meadow 22 , we hypothesize that rainy season onset (RSO) rather than temperature mainly drives the spatiotemporal variability of spring vegetation green-up (SVG) across alpine dry ecosystems over the Tibetan Plateau. Such knowledge is important to understand the mechanisms underlying variable and complex responses of spring vegetation phenology to climate warming in the alpine dry region.
To test the hypothesis, dates of SVG and RSO across 67 target areas locating nearby 67 weather stations were calculated from the time-series data of daily mean temperature and precipitation  and the Normalized Difference Vegetation Index (NDVI) from the Moderate Resolution Imaging Spectroradiometer (MODIS, 2001(MODIS, -2013, respectively. First, we tested whether the satellite-derived SVG dates compare well to our 7-year observations (2007-2013) for leaf emergence of two dominant species in alpine meadows along elevations (4400-5200 m) in Damxung of Tibet. Second, we investigated whether SVG dates of 2001-2013 across vegetation types (alpine dry shrub, alpine meadow, and alpine steppe) generally synchronize with RSO dates although seasonal air temperatures are already continuously above 0 °C before SVG dates. Third, we examined whether there are similar patterns in interannual change trends of SVG and RSO across the 67 sites and among the three vegetation types. As earlier onset of Indian summer monsoon may result in earlier RSO and higher pre-monsoon precipitation 22,25 , we also examined whether SVG generally shows a positive correlation with RSO and a negative correlation with pre-SVG precipitation. Finally, we identified the key factor driving SVG by quantifying the relative effects of RSO and pre-monsoon temperature and precipitation and their interactions on the interannual variation of SVG across the 67 sites.

Data and methods
Target areas. By means of Google Earth and the Vegetation Atlas of China at Scale 1:1 Million 27 , 67 target areas (> 1000 pixels of 250 × 250 m MODIS image for each target area) locating nearby 67 weather stations (with a distance < 5 km) of China Meteorological Administration were selected to represent the typical vegetation types in arid and semi-arid zones of the Tibetan Plateau. The 67 target areas included 18 for alpine steppes, 33 for alpine meadows, and 16 for alpine dry shrubs. Each target area was selected on the same slope with an elevation difference of < 300 m relative to the weather station. Deserts and evergreen shrubs/woodlands were excluded in this study because the resolution of MODIS NDVI data is not enough to detect the small variation in seasonality of leaf area in desert plants and evergreen canopies 28 . To reduce the possible effects of irrigation and water level, croplands and wetlands were also excluded. Detailed information on geographical locations, vegetation types and related climate factors of the 67 target areas is found in Table S1.  Table S2). Using the pixel-specific QA data, we conducted a data preprocessing procedure to remove the noises of snow cover and clouds according to the commonly accepted methods in the literature 15,[29][30][31][32] . In alpine meadows and steppes, there is short-time snow cover in the non-growing season (November to March), which often reduces NDVI values and then leads to errors in retrievals of spring phenology 30,33 . To remove the snow-cover effect on NDVI, we first calculated the uncontaminated winter NDVI for each pixel as the mean value of 75-95% quantile of the uncontaminated winter NDVI values (flagged as 'Good Data') during November to March 15,31,32 . The contaminated winter NDVI value (flagged as 'Snow/Ice' or 'Cloudy') that is lower than the uncontaminated winter NDVI was substituted by the latter 30 .

Scientific Reports
| (2020) 10:18797 | https://doi.org/10.1038/s41598-020-75991-w www.nature.com/scientificreports/ Clouds and poor atmospheric conditions generally depress NDVI values, which may result in an abrupt drop of NDVI in the growing season 29 . When the NDVI data were not flagged as 'Good Data' , the abrupt drops of NDVI values during April and August (prior to the occurrence of maximum NDVI) were reconstructed by using the method of Savitzky-Golay filter 29 . To remove the noises of bare land and rock fields with sparse vegetation, each pixel for calculation should meet the following 3 criteria 15,16,32 : (1) the average NDVI during June and September should be higher than 0.10; (2) the annual maximum NDVI should be higher than 0.15 and occur in the growing season of July to September; (3) the average NDVI of July-September should be 1.2 times higher than that of November-March.
The preprocessed data of time-series NDVI were used for calculating SVG dates according to the relative threshold method of NDVI ratio developed by White, Thornton and Running 28 . The NDVI ratio (i.e., the ratio of NDVI change to annual maximum amplitude) was calculated as: NDVI ratio = (NDVI t − NDVI min )/(NDVImax − NDVI min ) × 100%. NDVI t is the NDVI value at a certain time t; NDVI max and NDVI min are the maximum and minimum NDVI values, respectively, during the annual NDVI cycle. The SVG date for each year and each site was defined as the date when a NDVI ratio reached a threshold of 20%, as indicated in previous research on the Tibetan Plateau 8 .
To validate the methods used in this study, we tested the extent to which the satellite-derived SVG dates are consistent with our 7-year observations (2007-2013) for leaf emergence of two dominant species (Kobresia pygmaea sedge at 4650-5200 m; Stipa capillacea grass at 4400-4500 m) in alpine meadows along a large elevation gradient in Damxung of Tibet. The observed dataset was obtained from our recent report 22 .

RSO date and other climate factors.
Time-series data of daily mean temperature and precipitation across the 67 weather stations during 1974-2013 were obtained from China Meteorological Administration (https ://cdc.nmic.cn/home.do). In each of the 67 weather stations, the RSO dates during 1974-2013 were calculated from the time-series data of daily mean temperature and precipitation. There is still no consensus on the most appropriate definition of rainy season onset 23,34 . In previous studies, rainy season onset is usually determined by using a variety of empirical precipitation thresholds which depend on the growth requirements of plants under different climatic conditions 23,[35][36][37] . In this study, the RSO date was defined as the first day when a 5-days moving-average daily rainfall was higher than a constant threshold value, in which case the daily mean temperature was continuously above 0 °C 22 . By using different rainfall thresholds (ranging from 1 to 5 mm by a step of 0.5 mm) within a weather station, we calculated the correlation coefficients between satellite-derived SVG dates and the estimates of rainfall onset during 2001-2013. According to the highest correlation coefficient, the optimum rainfall threshold was determined for a weather station and then was used for calculating its RSO dates during 1974-2013. In this way, the interannual variation of RSO dates was independent of the constant, optimum rainfall threshold used for a weather station. Because of the geographical variations in rainfall intensity and plant growth requirements, the estimated optimum rainfall thresholds varied with study sites and vegetation types (Fig. S2a), ranging from 3.2 to 3.5 mm in steppes and meadows to 4.0 mm in shrubs (Fig. S2b).
To further explore ecologically meaningful, interannually comparable climate factors driving the change of SVG, the mean temperature and precipitation of 30 or 60 days before a multi-year mean SVG date (T -30d and Pr -30d or T -60d and Pr -30d ; so called pre-SVG temperature and precipitation) were calculated as in Li et al. 22 .
Data analysis. Simple linear regression was used for testing the interannual change trends in SVG, RSO and pre-SVG temperature (T -30d , T -60d ) and precipitation (Pr -30d , Pr -60d ) as well as their relationships to each other. At site level and in grouped data for shrubs (16 sites), meadows (33 sites) and steppes (18 sites), partial correlation analysis of multiple linear regressions was used for assessing the relative importance of RSO, T -30d and Pr -30d in determining the interannual variation of SVG during 2001-2013. In pooled data across the 67 sites during 2001-2013, the analysis of variance with Type III sum of squares under the framework of linear mixed model was used for quantifying the relative effects of RSO, T -30d , Pr -30d and their interactions on the variation of SVG, in which the station ID was taken as the subject; RSO, T -30d , Pr -30d and their interactions were taken as fixed factors, while random effects of RSO, T -30d and Pr -30d were considered.
All statistical analyses were performed using SPSS 18.0 (SPSS Inc., Chicago, Illinois, USA) and OriginPro 9.0 for Windows (OriginLab Corporation, Northampton, USA). All the maps were drawn by ArcGIS 10.3 for Desktop (Environmental Systems Research Institute, Inc., RedLands, California, USA). All significant correlations were taken at P < 0.05.

Results
Comparison between satellite-derived SVG and ground observations. The satellite-derived SVG dates synchronized well with the 7-year observations of leaf emergence dates (2007-2013) for two dominant species of alpine meadows along elevations from 4400 to 4800 m (Fig. 1a-d; R = 0.90, P < 0.05), but not at higher elevations close to the upper limit of Kobresia meadows (4950-5200 m, Fig. 1e-g). As the elevations of all the 67 target areas were less than 4800 m (Table S1), our observed data validated the methods for satellite-derived SVG dates at elevations below 4800 m on the Tibetan Plateau.
During 2001-2013, T -30d had an increasing trend in 28 of the 67 sites (P < 0.05 for 17 sites, P < 0.10 for 11 sites) but varied little in remaining other sites (Fig. S3a). Pre-SVG precipitation (Pr -30d ) varied little in most cases (except 5 sites with a significant decreasing or increasing trend; Fig. S3c). RSO dates were negatively correlated with pre-SVG precipitation (Pr -30d , Pr -60d ) (Fig. 3a,b, P < 0.05 for 32/37 sites, P < 0.10 for 11/8 sites), but varied little with pre-SVG mean temperatures in most cases (Fig. 3c,d, except for 8-11 of the 67 sites where there were positive or negative correlations between RSO and T -30d /T -60d ). Thus, SVG dates of 2001-2013 showed a general negative correlation to pre-SVG precipitation in almost half of the 67 sites (Fig. 4c,d). In comparison, SVG dates were not correlated with pre-SVG mean temperatures (T -30d , T -60d ) in most cases, only in 9-10 of the 67 sites where there were positive or negative correlations between SVG and T -30d (T -60d ) (Fig. 4a,b).
SVG dates of 2001-2013 showed an advancing trend in 19 of the 67 sites (P < 0.05 for 16 sites, P < 0.10 for 3 sites) but typically varied little in remaining other sites (except one site with a significant delaying trend; Fig. 5a). RSO showed an advancing trend in 10 of the 67 sites during 2001-2013 (Fig. 5b, P < 0.05 for 5 sites, P < 0.10 for 5 sites) and in 29 of the 67 sites during 1974-2013 (Fig. 5c, P < 0.05 for 26 sites, P < 0.10 for 3 sites), but varied little in remaining other sites (Fig. 5b,c). In grouped data for shrubs, meadows and steppes, mean RSO dates generally showed an advancing trend during 1974-2013 (P < 0.001, Fig. 5d-f). During 2001-2013, mean RSO and SVG dates both advanced significantly in meadows and steppes but varied little in shrubs (Fig. 5d-f), in which SVG dates synchronized with RSO dates in meadows and shrubs (Fig. 5d,e) and were somewhat later than RSO dates in steppes (Fig. 5f).

Relative effects of RSO and other climate factors on interannual variation of SVG.
Partial correlation analysis indicated that RSO had the highest significant effect on SVG in grouped data for each of the 3 vegetation types (R = 0.52-0.73, P < 0.001), compared to lower or insignificant effects of T -30d and Pr -30d (Table 1; the highest effect of RSO on SVG was also found in site-level data, Table S3).

Discussion
The Tibetan Plateau, which is known as the roof of the world, is one of the regions with strongest warming 38 . In recent decades, the whole plateau has become warmer and wetter 39 , which is expected to have a positive effect on the world's highest and largest alpine vegetation (meadows, steppes, dry shrubs, etc.) [40][41][42] . In response to the recent warming and wetting, however, there are variable trends (advanced, delayed, unchanged) in satellitederived green-up dates 8,11,12,[15][16][17] and long-term observed leaf emergence of dominant species 9,10,14 . In this study, our data indicated that SVG dates across alpine dry ecosystems generally synchronized with or were somewhat later than RSO dates although seasonal air temperatures were already continuously above 0 °C in 1 month before SVG dates (Figs. 2 and 5). Earlier RSO and higher pre-monsoon precipitation are generally linked to earlier onset of Indian summer monsoon 22,25 , which explains why earlier RSO dates typically resulted in higher pre-SVG precipitation (Fig. 3a,b) regardless of pre-SVG temperature (Fig. 3c,d). Accordingly, SVG dates showed a general negative correlation with pre-SVG precipitation (Fig. 4c,d) but did not correlate with pre-SVG temperature in www.nature.com/scientificreports/ most cases (Fig. 4a,b). There were similar patterns in change trends of RSO and SVG across the 67 sites and among the three vegetation types (Fig. 5). The analysis of linear mixed model further indicated that RSO and its interactions with pre-SVG temperature and precipitation mainly explained the spatio-temporal variability of SVG dates across arid and semi-arid ecosystems on the Tibetan Plateau, while pre-SVG temperature and its interaction with pre-SVG precipitation did not have significant effects on SVG ( Table 2). The findings suggest that the cold and dry climate in the pre-monsoon season is one of main filters for the survival and distribution of dominant species in alpine dry ecosystems, and phenological sensitivity to rainy season onset allows plants to simultaneously avoid drought and frost damages in the early growing season. Such a mechanism provides a general explanation for the spatially variable greening responses to climate change in the vast alpine grasslands and shrublands.
Since the arrival of monsoon rainfall mainly depends on the onset and intensity of Indian summer monsoon circulation and the complex effects of topography 43 , large spatial and temporal variations in RSO and SVG dates are expected (Figs. 3, 4, 5). In some cases, there were spatially positive or negative correlations between SVG and pre-SVG temperature (Fig. 4a,b) as reported in previous studies 14,32 , which would be mainly caused by the positive/negative correlations between RSO and pre-SVG temperature (Fig. 3c,d). It has been suggested that Indian summer monsoon is weakening and the westerlies are reinforcing since the 1970s, resulting in decreased precipitation in the Himalayas but increased precipitation in northern Tibetan Plateau 44 . In southern Tibetan Plateau, the weakened monsoon may lead to later onset of RSO and thus less precipitation and higher temperature (due to decreased cloud cover and increased solar radiation) in the pre-monsoon season, which explains the positive correlation between RSO and pre-SGV temperature (Fig. 3c,d) and the negative correlation between pre-SVG precipitation and temperature (Fig. S4a,b) at some sites in the southern part. In northern Tibetan Plateau, the enhanced westerlies may lead to increased spring snowfall and the spring warming may indirectly induce earlier    13 ) and/or the ratio of rainfall to snowfall 22 , which explains the negative correlations of RSO and SVG to pre-SGV temperature at some sites in the northern part (Fig. 3c,d).
It should be noted that at high elevations close to the upper limit of the Kobresia meadows (4950-5200 m, Fig. 1e-g), satellite-derived SVG dates (DOY 110-170) were much earlier than the observed leaf emergence dates (DOY 170-180), compared to the significant consistency between satellite-derived and observed values at lower elevations (4400-4800 m, Fig. 1a-d). In the early growing season (April to June) when precipitation often occurs at night time, it often snows and sleets at high elevations above 4800 m but rains at lower elevations. While the snowfall has no effect on plant growth, the early-season snow melting may increase the NDVI values and then result in incorrect estimation of SVG dates derived from satellite data at elevations above 4800 m. In this study, such uncertainty should be minor because the elevations of all the 67 target areas were less than 4800 m (Table S1). In the future, field experimental data are needed to verify the optimum rainfall threshold for RSO and the NDVI ratio threshold for SVG.
Our observed data in Damxung of Tibet further indicated that leaf emergence of dominant species in the alpine meadows was highly sensitive to the arrival of monsoon rainfall but insensitive to that of snowfall 22 , resulting in later green-up dates above 4800 m than at lower elevations ( Fig. 1). At the upper limit of the Kobresia meadows (5100-5200 m, Fig. 1f,g), the observed green-up dates varied little during 2007-2013 and typically occurred at the end of June (DOY 180) when the daily mean air-temperature was continuously above 0 °C, setting a threshold of minimum growing days that allows plants to have enough time to safely complete their life cycle. This is consistent with our previous reports that the number of days with daily mean soil temperature ≥ 5 °C (at − 10 cm) during June and August is about 80-90 days at 5100-5200 m 45 , with a seasonal mean soil temperature of 6.8-7.3 ℃ 46 being close to the low temperature threshold that significantly limits treeline growth and photosynthesis [47][48][49][50] .
The ratio of rainfall to snowfall usually decreases with increasing elevation. Spring warming may increase the rainfall fraction and then indirectly advance the green-up dates of alpine plants at high elevations. This is typically consistent with some previous observations on the spring phenological response to precipitation change in this region. Using observed data from 23 phenological stations on the Tibetan Plateau, Chen et al. 14 indicated that the different interactions between snowfall and temperature during late winter and early spring likely determine the spatiotemporal variations of green-up dates. In the central Himalaya, Pangtey et al. 51 found that the growth initiation of 184 alpine plant species generally synchronizes with the beginning of snow melt period.
In conclusion, our data supported the hypothesis, indicating that monsoon rainfall onset and its interactions with pre-monsoon temperature and precipitation mainly drives the spatio-temporal variability of spring green-up date across alpine dry ecosystems over the Tibetan Plateau. Our data further suggest that synchronization of spring vegetation green-up and rainy season onset is a general spring phenological adaptation to the pre-monsoon cold and dry climate. The threshold of minimum growing days determining the upper limit of alpine meadows/steppes would be linked to the sensitivity of spring phenology to the arrival of monsoon-season rainfall. Long-term observations and experiments, which are very scarce in remote high mountains, are required to verify the findings based on satellite-derived data. Table 2. Summary of Type III tests of linear mixed model for fixed effects of RSO, T -30d , Pr -30d and their interactions on the variation of SVG in pooled data across the 67 sites during 2001-2013. The fixed effects were tested by using analysis of variance with Type III sum of squares under the framework of linear mixed model. Denominator df is equal to within-group sum of squares.