## Data and Methods
### Data Collection  
Time series of observed daily maximum, minimum, and mean air temperature (TX, TN, and Ta, respectively) and daily total precipitation (P) from meteorological stations were obtained from the China Meteorological Administration (CMA) [39]. To cover the whole study area, a 60-km buffer from the edge of the province boundary was used to determine the meteorological stations used in this study. Then, according to the criteria that the series length should cover 1961-2015 and the missing data should be no more than 1% in every year at each station, data from 52 stations were ultimately selected for analyses, which included 30 stations within the study area and 22 around the study area (Figure 1). To clarify the relationship between indices of climatic extremes and large-scale circulation modes of climate variability, time series of monthly atmospheric and oceanic climate indices, including Multivariate ENSO Index (MEI), Pacific Decadal Oscillation (PDO) and Atlantic Multidecadal Oscillation (AMO), were collected from Earth System Research Laboratory, National Oceanic and Atmospheric Administration (NOAA) [40].  
Vegetation variabilities in Hubei Province were evaluated with long-term global mapping leaf area index (LAI) product (GLOBMAP LAI Version 3) [41] and third-generation normalized  
Figure 1: Location map of Hubei Province with elevation. The dark triangles denote the 52 meteorological stations used in this study.  
difference vegetation index product (NDVI 3g.v1) [42]. The GLOBMAP LAI product was generated by quantitative fusion of the Moderate Resolution Imaging Spectroradiometer (MODIS) and the historical Advanced Very High Resolution Radiometer (AVHRR) data [43]. The pixel-level relationships between the AVHRR and the MODIS LAI were first established during the overlapping period 2000-2006 and used to estimate the AVHRR LAI back to 1981. Consequently, the long-term LAI series was constructed by a combination of the AVHRR LAI (1981-2000) and the MODIS LAI (2000-2016), which provided a global LAI dataset at a 15-day temporal resolution and an 8-km spatial resolution from 1981 to 2016. The NDVI data was derived from the new Global Inventory Modeling and Mapping Studies (GIMMS) [44], which can be found at the NASA Earth Exchange (NEX). This NDVI data was produced from AVHRR instruments that extended from 1981 to the present. In addition, this data has been widely used in previous land surface studies as a spatio-temporal indicator for vegetation activity [45; 46; 47].  
A new gross primary productivity (GPP) data, with 0.1\({}^{\circ}\) spatial resolution and monthly temporal frequency over China, was obtained from [PERSON] et al. [35], where a machine learning algorithm with eddy flux measurements was applied to 40 sites in China and its surrounding countries. The high-resolution GPP data can successfully capture the spatio-temporal variations of the GPP observed at the flux sites, including validation flux site data that were not put in training of model tree ensemble (MTE). All of the gridded datasets were first composited to annual temporal resolution and then aggregated to 0.05\({}^{\circ}\)\(\times\) 0.05\({}^{\circ}\) spatial resolution using bilinear interpolation.  
Brief information for all data used in this research is listed in Table 1.