## Materials and Methods  
### Data Sources  
We used Moderate Resolution Imaging Spectroradiometer (MODIS) datasets, including global Gross Primary Productivity product (MOD17A2H V6) and land cover type products (MCD12Q1 V6), both of which have a spatial resolution of 500 m. MOD17A2H is a dataset containing cumulative 8-day composite and time-series GPP applicable to extracting vegetation phenology (SOS and EOS). This GPP dataset contains a quality control 8-bit band (Psn_QC) to flag poor quality pixels; the first bit with value zero demonstrates good quality. We selected the MODIS datasets from 2001 to 2020. The International Geosphere-Biosphere Programme (IGBP) from MCD12Q1 provides consistent vegetation cover classification schema with many global field observation programs. Particularly, FLUXNET consists of a global network of tower sites that use eddy covariance to measure the exchanges of carbon dioxide, water vapor, and energy between the biosphere and atmosphere [29].  
FLUXNET2015 is the latest freely downloadable FLUXNET data ([[https://fluxnet.org](https://fluxnet.org)]([https://fluxnet.org](https://fluxnet.org)), accessed on 10 August 2022). Compared with its predecessor, the 2007 LaThuile database, the FLUXNET2015 database has been improved in both data quality control and data processing procedures [8]. There are 212 stations worldwide [29]. We selected the data from station sites according to the following criteria. (1) The sites had complete and continuous observations for at least one year, (2) the vegetation type of the site was consistent with the MODIS IGBP type, and (3) quality control was performed (>80% of the sampled data show good quality) [8; 30]. The raw observations from the stations have a fine temporal resolution of a half hour/hour. The fine-temporal-scale datasets are aggregated into secondary products with daily, monthly, and annual resolutions. Sample points from the daytime (GPP_DT_VUT_REF) and nighttime (GPP_NT_VUT_REF) methods are available from the FLUXNET2015, but the phenology extracted from either of them was proved to be identical to each other, thus we only applied GPP_NT_VUT_REF to extract the site-based phenology to validate remote-sensing-derived phenological metrics. A detailed description of the process for preparing a good quality station-based dataset has been documented previously [8]. The IGBP land (vegetation) cover types from MCD12Q1 were further reclassified into five groups, including forest(s) (evergreen needleleaf forests, evergreen broadleaf forests, deciduous needleleaf forests, deciduous broadleaf forests, and mixed forests), shrub (open shrublands and closed shrublands), savannas (woody savannas and savannas), grass (grasslands), and crop (croplands and cropland/natural vegetation mosaics). The vegetation cover and the locations of FLUXNET stations are presented in Figure 1.  
In addition, monthly average temperature and total precipitation were extracted from the ERA5-Land with a spatial resolution of \(\sim\)10 km [31]. ERA5-Land was produced by replaying the land component of the ECMWF ERA5 climate reanalysis [32] and made available to facilitate many applications requiring easy and fast access to the data. The dataset is freely available in the Google Earth Engine (GEE) platform and has been applied in many ecological studies [33; 34; 35]. To make all data have an identical resolution, both temperature (unit \({}^{\circ}\)K) and precipitation (in centimeter/month) were resampled to 500 m by bi-cubic interpolation.