*2.2. Data Source*

The GIMMS NDVI3g dataset provided by NASA was used to estimate the EOS on the Qinghai-Tibetan Plateau. The dataset was available from 1982 to 2015 with an 8-km spatial resolution and 15-day temporal resolution [32]. Some previous processes (e.g., calibration, noise removal) were performed for this version to better detect the vegetation dynamics [32]. This dataset has been widely used to detect long-term vegetation dynamics [33–35]. Due to the normalized difference vegetation index (NDVI) data might be misrepresented by snow [36]; we used the average temperature of a sequence of five days less than 0 ◦C to screen out the pixels that might be covered by snow. Temperature, precipitation, and insolation data from 1982–2015 were extracted from the China meteorological forcing dataset (1979–2015) downloaded from the Big Earth Data Platform for Three Poles with a spatial resolution of 0.1◦ and temporal resolution of 3 h (http://poles.tpdc.ac.cn/, Accessed on 15 August 2021) [37].

Human activities, including grazing density and economic development, were quantified using economic statistic data. For example, the grazing density were represented by the number of large animals (one large animal equal to five sheep unit) and sheep, and uniformly converted into sheep units. The economic development levels were quantified as the production of primary, secondary, and tertiary sectors. These data come from the statistical yearbooks of Qinghai and Tibet from 1982 to 2015.

## *2.3. Retrieval of EOS*

Numerous methods have been used to fit the NDVI changes from seasonal vegetation cycles. After comparing the fitting results of HANTS [38], Polyfit [39], and Double logistic [40] in the nine subregions, we found that the RMES of HANTS (1.26 ± 0.24 × <sup>10</sup>−5) and Polyfit (1.28 ± 0.24 × <sup>10</sup>−5) were similar and smaller than the Double logistic results (1.93 ± 0.43 × <sup>10</sup>−5) (Figure A1). HANTS and Polyfit, were therefore selected as the two most simple and effective methods to fit the NDVI change curves. Dynamic thresholds were adopted to determine the EOS. Further details of these two fitting methods are described below.

The HANTS method involves the harmonic analysis of a time series, is adapted from the fast Fourier transform, and eliminates cloud noise using the least square method [38,40]. The HANTS methods can quickly smooth the data, remove outliers, and fill gaps of missing data. The following Equation (1) was used to fit the NDVI seasonal fluctuation curve:

$$NDVI(t) = a\_0 + \sum\_{i=1}^{n} a\_i \cos(2\pi t - \varphi\_i) \tag{1}$$

where *t* is the Julian date, *a*0 is the average of all *NDVI* observations, and *ϕi* and *ai* are the phase and amplitude of the curve, respectively.

The Polyfit method adopts a polynomial function to fit the *NDVI* records [39]. The following sixth order Equation (2) is used to describe the *NDVI* curve:

$$NDVI(t) = a\_0 + a\_1t + a\_2t^2 + \dots + a\_6t^6 \tag{2}$$

where *a*0–*a*6 are regression coefficients determined using the Levenberg–Marquardt method.

The EOS was determined by the day when the smoothed curve of the 34-year mean passed a designated threshold. We first fitted the *NDVI* changes with HANTS and Polyfit methods and then calculated the NDVIratio (described in Equation (3)) for 365 days with multi-year mean *NDVI* values, next detected the time *t* with the minimum *NDVI*ratio and used the corresponding *NDVI*(*t* + 1) at time (*t* + 1) as the *NDVI* threshold for the EOS. Finally, we obtained the EOS for 34 years using the threshold:

$$NDVI\_{ratio}(t) = \frac{NDVI(t+1) - NDVI(t)}{NDVI(t)}\tag{3}$$

#### *2.4. Quantification of the EOS Trends, Turning Points, and Controls*

After extracting the EOS at the pixel scale, we first quantified the tendency of the EOS changes using greenness changes methods, and then detected the turning points using the piecewise regression method. The mean EOS values and EOS trends within the subregions were calculated as the EOS and EOS trends at the subregion level. The turning points at a subregion and province level were calculated by the majority values.

The EOS trends were calculated using the greenness rate of change [41]. The EOS was considered delayed if the slope was a positive value; otherwise, the EOS advanced.

$$slope = \frac{n \times \sum\_{i=1}^{n} (i \times NDVI) - \sum\_{i=1}^{n} i \sum\_{i=1}^{n} NDVI\_i}{n \times \sum\_{i=1}^{n} i^2 - \left(\sum\_{i=1}^{n} i\right)^2} \tag{4}$$

where *i* is the order of the year, *n* is the number of years, *NDVIi* is the *NDVI* in the *i*th year, and the slope is the vegetation change rate. Alternatively, we can use the unary linear regression, in which the *P* values and confidence levels can be calculated.

Turning points were identified by piecewise regression [42] analysis, as defined in Equation (5), which can be used to detect sudden and sharp changes in directionality. This method has been widely applied for analyzing vegetation dynamics [19,43,44].

$$y = \begin{cases} \begin{array}{c} \beta\_0 + \beta\_1 t + \varepsilon \ t \le a \\ \beta\_0 + \beta\_1 t + \beta\_2 (t - a) + \varepsilon \ t > a \end{array} \end{cases} \tag{5}$$

where *t* is the order of the year, *α* is the estimated turning point of the vegetation change trend determined using the least square error method, *β*1 and (*β*1 + *β*2) are the change rates before and after the turning points, respectively, and *ε* is the residual error. We performed *t*-tests to check the significance of the piecewise regressions.

Redundancy analysis (RDA) is a powerful analysis technique that could be applied in separating the contributions of climate, human activities, and their intersections to the EOS changes. RDA is a method to extract and summarize the variation in a set of response variables that can be explained by a set of explanatory variables [45]. In this study, RDA was performed with the vegan package in R language [46]. In RDA, climatic variables or human activity variables were chosen as predictors to maximize the extent of their correlation with the EOS changes as the response variable. RDA had been widely used in ecology-related studies [47,48]. The turning points of human activities were also calculated with Equation (5). The relationships between the turning points of the EOS and climatic variables were quantified using partial regression analysis or the correlation coefficient.
