Heat and Drought Stress Advanced Global Wheat Harvest Timing from 1981–2014

Studying wheat phenology can greatly enhance our understanding of how wheat growth responds to climate change, and guide us to reasonably confront its influence. However, comprehensive global-scale wheat phenology–climate analysis is still lacking. In this study, we extracted the wheat harvest date (WHD) from 1981–2014 from satellite data using threshold-, logistic-, and shape-based methods. Then, we analyzed the effects of heat and drought stress on WHD based on gridded daily temperature and monthly drought data (the Palmer drought severity index (PDSI) and the standardized precipitation evapotranspiration index (SPEI)) over global wheat-growing areas. The results show that WHD was generally delayed from the low to mid latitudes. With respect to variation trends, we detected a significant advancement of WHD in 32.1% of the world’s wheat-growing areas since 1981, with an average changing rate of −0.25 days/yr. A significant negative correlation was identified between WHD and the prior three months’ normal-growing-degree-days across 50.4% of the study region, which implies that greater preseason effective temperature accumulation may cause WHD to occur earlier. Meanwhile, WHD was also found to be significantly and negatively correlated with the prior three months’ extreme-growing-degree-days across only 9.6% of the study region (mainly located in northern South Asia and north Central-West Asia). The effects of extreme heat stress were weaker than those of normal thermal conditions. When extreme drought (measured by PDSI/SPEI) occurred in the current month, in the month prior to WHD, and in the second month prior to WHD, it forced WHD to advance by about 9.0/8.1 days, 13.8/12.2 days, and 10.8/5.3 days compared to normal conditions, respectively. In conclusion, we highlight the effects that heat and drought stress have on advancing wheat harvest timing, which should be a research focus under future climate change.


Introduction
Vegetation phenology has attracted a lot of attention in the face of global climate change over the past several decades. Due to strong dependence on climate conditions, vegetation phenology dynamics are an effective tool to assess the magnitude of regional and global climate change [1][2][3]. Meanwhile, vegetation phenology spatiotemporal variability determines, to some extent, the terrestrial carbon budget and ecosystem productivity, and even impacts local climate [3,4]. Therefore, studying vegetation phenology can enhance our knowledge of how vegetation growth responds to climate change and how to confront the influence of global change. However, most studies have focused on the phenological variability of natural vegetation and only a few have investigated the phenological changes of agricultural plants [5][6][7][8].
Given that wheat is the most widely grown crop in the world, changes in its phenology could greatly affect regional land surface material and energy exchange, as well as yields [9,10]. Using long-term monitoring, many studies have detected a significant advance in wheat phenology events (such as the timing of shooting, heading, and yellow ripeness) over the past several decades in Germany [11], the U.S. Great Plains [6], and China [8,12,13], with rates of change from 0.1 days/yr to 0.31 days/yr. All of these shifts were primarily attributed to the significant increase in average preseason or seasonal temperatures. Remote sensing data provide an excellent opportunity to study the spatiotemporal patterns of wheat phenology at regional and global scales [14,15]. Using Moderate Resolution Imaging Spectroradiometer (MODIS) data from 2000-2009, Lobell et al. [14] found that extremely high temperatures (>34 • C) contributed more to the advancement of wheat senescence timing than average temperature increases in India. In addition to air temperature, water supply also impacts wheat growth by altering photosynthesis, photochemical activities, and enzyme activities [16,17]. However, the role of water conditions on wheat growth is still debated. Through statistical analysis and controlled experiments, some studies have found that increased precipitation or irrigation can slow wheat development and lengthen the growing period [14,18], while another study found that precipitation had no significant effects on wheat heading and flowering dates [6]. Therefore, studying the impacts of heat and moisture stress on wheat phenology is of crucial importance to wheat production and variety breeding.
Additionally, most studies on wheat phenology have focused on specific sites or covered a small regional scale. Global-scale wheat phenology-climate analysis is still lacking. In the current study, we used three methods to extract wheat harvest date (WHD) from satellite data covering 1981-2014 over global wheat-growing areas. Then we analyzed the effects of heat and drought stress on WHD and evaluated spatial differences at the global scale using gridded temperature and drought data. The goals of this study were to (1) evaluate the impacts of normal-growing-degree-days and extreme-growing-degree-days on WHD, and (2) assess how WHD responds to drought stress.

Global Wheat Planting Map
In the present study, global wheat distribution was identified using the spatial production allocation model 2005 (SPAM2005, http://mapspam.info/global-data), which was released by the International Food Policy Research Institute [19]. The SPAM2005 relies on a collection of spatially explicit input data covering three years around 2005, including crop production statistics, cropland data, biophysical crop "suitability" assessments, population density, and any prior knowledge about the spatial distribution of specific crops or crop systems. The spatial resolution is 5 arc minutes (0.083 • , 10 km). To reduce the effects of other crops and natural vegetation types as much as possible, we only focused on grids containing greater than 30% wheat. This wheat-planting extent was used for all years. Here, we assume there is no change in wheat-planting area during 1981-2014. In addition, to clarify the spatial differences in wheat phonological response to climate change, we divided the study into nine major wheat-planting regions, i.e., North America, South America, Europe, North Africa, Central-West Asia, East Africa, South Asia, China, and Australia ( Figure 1).

Figure 1.
Global wheat (green dots) distribution map. (i)-(ix) represent the nine major wheat-planting regions, i.e., North America, South America, Europe, North Africa, Central-West Asia, East Africa, South Asia, China, and Australia.

Remote Sensing Data and Processing
We estimated WHD from two-band enhanced vegetation index data (EVI2, 1981(EVI2, -2014 [20]. This product has a spatial resolution of 0.05° and a temporal resolution of 7 days. We first chronologically joined the EVI2 time series across all years and removed poor-quality EVI2 observations (contaminated by cloud, snow, and ice, or affected by the viewing geometry (solar zenith angle >75°)) for each pixel based on the associated quality data. If the poor-quality EVI2 observations accounted for more than 30% of the total annual observations, the data from that year were discarded. Next, we smoothed the EVI2 time series and interpolated it to daily resolution using a cubic smooth spline function, which could simulate the continuous growth process of plants and create a curve with high smoothness. Using the growing calendar for global wheat summarized by the United States Department of Agriculture [21], the corresponding EVI2 wheat profile was further truncated from the full-year EVI2 time series. To eliminate the effects of other land cover types that experience small seasonal changes, we removed pixels with a multiyear average EVI2 seasonal amplitude time series less than 0.1. Finally, three methods, i.e., threshold-, logistic-, and shape-based methods, were used to estimate WHD from the wheat EVI2 time series for each pixel, respectively.
(1) Threshold-based method. This method defined the WHD as the time point corresponding to 20% of the wheat EVI2 time series amplitude each year, which has been widely used to estimate phenological parameters from remote sensing data [14,15,22].
(2) Logistic-based method. A double logistic function (eq 1) was used to fit the annual wheat EVI2 time series and calculate WHD (eq 2) based on the maximum curvatures of the fitted curve: where base denotes the annual minimum EVI2 value; frange represents the range from spring minimum EVI2 to annual maximum EVI2; srange represents the range from autumn minimum EVI2

Remote Sensing Data and Processing
We estimated WHD from two-band enhanced vegetation index data (EVI2, 1981-2014), released by the University of Arizona Vegetation Index & Phenology Laboratory (https://vip.arizona.edu/ viplab_data_explorer.php). The data combine Advanced Very High Resolution Radiometer (AVHRR, 1981(AVHRR, -1999 and MODIS (2000-2014) datasets by standardizing AVHRR to the same level of MODIS [20]. This product has a spatial resolution of 0.05 • and a temporal resolution of 7 days.
We first chronologically joined the EVI2 time series across all years and removed poor-quality EVI2 observations (contaminated by cloud, snow, and ice, or affected by the viewing geometry (solar zenith angle >75 • )) for each pixel based on the associated quality data. If the poor-quality EVI2 observations accounted for more than 30% of the total annual observations, the data from that year were discarded. Next, we smoothed the EVI2 time series and interpolated it to daily resolution using a cubic smooth spline function, which could simulate the continuous growth process of plants and create a curve with high smoothness. Using the growing calendar for global wheat summarized by the United States Department of Agriculture [21], the corresponding EVI2 wheat profile was further truncated from the full-year EVI2 time series. To eliminate the effects of other land cover types that experience small seasonal changes, we removed pixels with a multiyear average EVI2 seasonal amplitude time series less than 0.1. Finally, three methods, i.e., threshold-, logistic-, and shape-based methods, were used to estimate WHD from the wheat EVI2 time series for each pixel, respectively.
(1) Threshold-based method. This method defined the WHD as the time point corresponding to 20% of the wheat EVI2 time series amplitude each year, which has been widely used to estimate phenological parameters from remote sensing data [14,15,22].
(2) Logistic-based method. A double logistic function (Equation (1)) was used to fit the annual wheat EVI2 time series and calculate WHD (Equation (2)) based on the maximum curvatures of the fitted curve: where base denotes the annual minimum EVI2 value; frange represents the range from spring minimum EVI2 to annual maximum EVI2; srange represents the range from autumn minimum EVI2 to annual maximum EVI2; mS and mA are the maximum slope values of the fitted curve in spring and autumn, respectively; S and A are the corresponding dates of mS and mA, respectively. The logistic-based method has been used extensively to retrieve vegetation phenology from remote sensing data and has had better performance in monitoring phenology for different land cover types [23][24][25][26].
(3) Shape-based method. This method assumes that the shape of the vegetation index (VI) growth curve (named "shape" curve) is relatively stable for a specific crop in a given location [27], which was defined as the multiyear mean EVI2 time series in this study. Geometric scaling differences exist between the annual VI curves and the shape curve due to local weather conditions and management types. The shifted shape curve for each year was obtained by fitting the annual EVI2 time series f(x) to the shape curve h(x) (Equation (3)): where x is the day of year; xscale, tshift, and yscale are scaling parameters that were determined based on the smallest root mean square error (RMSE, Equation (4)) between the shifted shape curve and the raw EVI2 time series for each year using the Markov chain Monte Carlo method.
where, P EVI2 is the estimated EVI2 value from Equation (1); R EVI2 is the observed EVI2 value; n is the number of observations. A predefined WHD was extracted from the shape curve using the threshold-based method. Annual WHD values were then acquired by entering the estimated parameters and the predefined WHD into Equation (3). For the logistic-based and shape-based methods, we also computed the goodness-of-fit (R 2 ) and RMSE between the fitted and the raw EVI2 values. To ensure the accuracy of the results, estimated WHD values were discarded if R 2 was less than 0.8, or if RMSE was greater than 0.05.

Climate Data
Daily mean air temperatures (1981-2014) from the gridded WATCH-Forcing-Data-ERA-Interim (WFDEI) dataset (http://www.eu-watch.org) were provided by the European Union Water and Global Change (WATCH) project [28], and were used to analyze the thermal impacts on temporal WHD patterns. The dataset has a spatial resolution of 0.5 • . When comparing with 2 m temperature observed with flux tower, the average error of WFDEI air temperature is about 0.79 • C. It has been previously used by Stacke and Hagemann [29] to drive the Max Planck Institute-Hydrology Model for global wetland mapping.

Drought Data
Two gridded drought index products, the Palmer drought severity index (PDSI) and the standardized precipitation evapotranspiration index (SPEI), were selected to study the effects of drought stress on global WHD. The PDSI was originally designed by Palmer [30] to characterize drought conditions of the United States and further developed by Wells et al. [31] to make it suitable for worldwide drought monitoring through self-calibration based on local conditions. The PDSI value is derived by a physical water-balance model, which uses both precipitation and air temperature at 2 m as inputs and considers preceding conditions. The PDSI data product used in this study was acquired from the Climatic Research Unit (CRU, https://crudata.uea.ac.uk/cru/data/drought/) and driven by the CRU TS3.26 monthly climate dataset [31,32].
SPEI is also extensively used for global drought monitoring but has the advantage of being multiscalar (1−48 months) over PDSI (9−12 months) [33]. It is based on monthly precipitation and potential evapotranspiration from the CRU TS 3.23 dataset. In this study, three-month SPEI data were downloaded from the Global SPEI database (https://digital.csic.es/handle/10261/153475). Both PDSI and SPEI have a spatial resolution of 0.5 • × 0.5 • and a temporal resolution of 1 month. According to PDSI and SPEI, drought was divided into five levels: near normal (NN), mild drought (MID), moderate drought (MOD), severe drought (SD), and extreme drought (ED) (

Analysis
Estimated WHD values from the three methods were first compared and averaged for each pixel. To match the spatial resolution of the climate data, WHD (0.05 • × 0.05 • ) values were then resampled up to 0.5 • × 0.5 • by averaging all the records within a given climate grid. To explore the effects of normal and extreme thermal conditions on WHD, we calculated the normal-growing-degree-days (NDD, Equation (5)) and the extreme-growing-degree-days (EDD, Equation (6)) within 90 days prior to the WHD based on daily mean air temperatures for each climate grid. Considering the uncertainties of maximum temperature threshold for wheat growth [35] and such a large spatial range, we took 30 • C as the threshold to distinguish NDD and EDD.
where T is daily air temperature. Pearson's correlation coefficients were calculated for WHD and NDD and WHD and EDD at each grid.
To assess the impacts of drought conditions on wheat harvest, we first grouped WHD values according to drought levels for the current month, the previous month before WHD, and the second month prior to WHD in various years at each grid. WHD values at the same drought level were then averaged. Finally, WHDs were compared among different drought levels across the whole study region. In addition, long-term trends of WHD, NDD, EDD, PDSI, and SPEI were obtained using a Theil-Sen estimator, which selects the median of the slopes of all lines through pairs of two-dimensional sample points [36,37]. The Theil-Sen estimator is a robust method for detecting temporal series trends and is not affected by extreme outliers. The nonparametric Mann-Kendall test was then used to conduct significance testing for phenological trends at the level of 0.05.

Spatial Pattern of WHD
The results show that WHD values calculated using the logistic-based method occurred slightly earlier by an average of 6.0 days than those derived by the threshold-based method, while WHD values from the shape-based method occurred slightly later by about 6.5 days than those from the threshold-based method. Spatially, WHDs estimated by all three methods were generally delayed from the low to mid latitudes (Figure S1a-c). By averaging WHD values from the three methods over all the years in each grid, we found that except for Chile, where wheat is harvested in late December or early January, the earliest WHD values occurred primarily in southern South Asia (from late February to March) ( Figure 2). In northern South Asia, as well as Egypt and Morocco in North Africa, farmers usually harvest wheat from April to mid-May. In China, West Asia, Europe, as well as southern Central Asia and North America, wheat normally matures from mid-May to mid-August. For the northernmost wheat-growing areas in Canada, Russia, and Kazakhstan, wheat is usually harvested in late August and September. By comparison, WHD in Ethiopia and the Southern Hemisphere generally range from October to the end of the year.

Variation Trends in WHD and Climate
At the pixel scale, a significant advancement in WHD was detected in 32.1% of the world's wheat-growing areas since 1981, with an average rate of change of −0.25 days/yr ( Figure 3). However, the direction and magnitude of WHD variability were different across the nine wheat-growing areas. During the past 34 years, WHD was significantly delayed in North America (0.15 days/yr), but significantly advanced in Europe, East Africa, and Australia by 6.5, 2.4, and 12.9 days, respectively ( Figure 4). For China, a shift in WHD variability occurred around 2000, namely, gradually advancing between 1985 to 2000 and becoming delayed after 2000.

Variation Trends in WHD and Climate
At the pixel scale, a significant advancement in WHD was detected in 32.1% of the world's wheat-growing areas since 1981, with an average rate of change of −0.25 days/yr ( Figure 3). However, the direction and magnitude of WHD variability were different across the nine wheat-growing areas. During the past 34 years, WHD was significantly delayed in North America (0.15 days/yr), but significantly advanced in Europe, East Africa, and Australia by 6.5, 2.4, and 12.9 days, respectively (

Variation Trends in WHD and Climate
At the pixel scale, a significant advancement in WHD was detected in 32.1% of the world's wheat-growing areas since 1981, with an average rate of change of −0.25 days/yr ( Figure 3). However, the direction and magnitude of WHD variability were different across the nine wheat-growing areas. During the past 34 years, WHD was significantly delayed in North America (0.15 days/yr), but significantly advanced in Europe, East Africa, and Australia by 6.5, 2.4, and 12.9 days, respectively ( Figure 4). For China, a shift in WHD variability occurred around 2000, namely, gradually advancing between 1985 to 2000 and becoming delayed after 2000.

Interannual WHD Fluctuation in Relation to Thermal Conditions
A significant negative correlation was identified between WHD and the prior three months' NDD in 50.4% of the study region, which means greater preharvest effective temperature accumulation would cause WHD to appear earlier. The stronger sensitivity of WHD to NDD was primarily detected in South Asia, Central-West Asia, Europe, and Australia ( Figure 5a). Meanwhile, WHD was also significantly and negatively correlated with the prior three months' EDD in 9.6% of the study region (mainly northern South Asia and northern Central-West Asia) (Figure 5b). In other words, higher preharvest EDD values led to earlier WHD.

Interannual WHD Fluctuation in Relation to Thermal Conditions
A significant negative correlation was identified between WHD and the prior three months' NDD in 50.4% of the study region, which means greater preharvest effective temperature accumulation would cause WHD to appear earlier. The stronger sensitivity of WHD to NDD was primarily detected in South Asia, Central-West Asia, Europe, and Australia ( Figure 5a). Meanwhile, WHD was also significantly and negatively correlated with the prior three months' EDD in 9.6% of the study region (mainly northern South Asia and northern Central-West Asia) (Figure 5b). In other words, higher preharvest EDD values led to earlier WHD.

Effect of Drought on WHD
According to the PDSI, the frequency of serious drought events (from moderate to extreme drought) which occurred in the current month, in the month prior to WHD, and in the second month prior to WHD was 25%, 25%, and 23.8% during 1981-2014, respectively ( Figure S3a-c). Spatially, places that were easily affected by serious drought events were primarily distributed in northwestern North America, southern India, Central-West Asia, and western Australia. For serious drought events measured based on SPEI, the frequency of serious drought events occurred in the current month, in the month prior to WHD, and in the second month prior to WHD was 15%, 15.6%, and 15.6%, respectively ( Figure S3d-f). After comparing the spatial patterns (across the whole study region) of SPEI-and PDSI-derived serious drought events based on correlation analysis ( Figure S4), we found the spatial distribution of SPEI-derived serious drought events for the month when WHD occurred was similar with that based on PDSI. For the second month before WHD, however, all South Asia and China were often impacted by serious drought events, while serious drought events scarcely took place in North America.
By averaging WHD values for each drought level, we found that WHD gradually occurred earlier as drought level increased. Extreme drought occurring in the current month forced wheat harvest timing to advance by about 9 days for PDSI and 8.1 days for SPEI compared to normal conditions (Figure 6a,d). Simultaneously, extreme drought occurred in the last/second month before

Effect of Drought on WHD
According to the PDSI, the frequency of serious drought events (from moderate to extreme drought) which occurred in the current month, in the month prior to WHD, and in the second month prior to WHD was 25%, 25%, and 23.8% during 1981-2014, respectively ( Figure S3a-c). Spatially, places that were easily affected by serious drought events were primarily distributed in northwestern North America, southern India, Central-West Asia, and western Australia. For serious drought events measured based on SPEI, the frequency of serious drought events occurred in the current month, in the month prior to WHD, and in the second month prior to WHD was 15%, 15.6%, and 15.6%, respectively ( Figure S3d-f). After comparing the spatial patterns (across the whole study region) of SPEI-and PDSI-derived serious drought events based on correlation analysis ( Figure S4), we found the spatial distribution of SPEI-derived serious drought events for the month when WHD occurred was similar with that based on PDSI. For the second month before WHD, however, all South Asia and China were often impacted by serious drought events, while serious drought events scarcely took place in North America.
By averaging WHD values for each drought level, we found that WHD gradually occurred earlier as drought level increased. Extreme drought occurring in the current month forced wheat harvest timing to advance by about 9 days for PDSI and 8.1 days for SPEI compared to normal conditions (Figure 6a,d). Simultaneously, extreme drought occurred in the last/second month before WHD caused it to advance by about 13.8/12.2 days for PDSI and 10.8/5.3 days for SPEI compared to normal years (Figure 6b,c,e,f). Similar results were also detected in most of the nine wheat-growing areas (Figure 7). Except for a slightly delayed WHD (1.5 days) triggered by severe drought (SPEI) in the last month before WHD in Europe, severe drought occurring in the current month when WHD occurred to the second month before WHD caused it to advance by 0.6-20 days compared to normal years in South America, Europe, North Africa, Central-West Asia, East Africa, China, and Australia. In North America, delayed WHD values (0.5-4.2 days) were observed under an extremely arid climate measured by PDSI, whereas advanced WHD (7.3-17.6 days) occurred under severe drought conditions measured by SPEI. On the contrary, in South Asia we identified advanced WHD values (7.4-10.6 days) in extremely dry years calculated by PDSI but delayed WHD (2.2-7.6 days) in extremely dry years calculated by SPEI.  (Figure 6b,c,e,f). Similar results were also detected in most of the nine wheat-growing areas (Figure 7). Except for a slightly delayed WHD (1.5 days) triggered by severe drought (SPEI) in the last month before WHD in Europe, severe drought occurring in the current month when WHD occurred to the second month before WHD caused it to advance by 0.6-20 days compared to normal years in South America, Europe, North Africa, Central-West Asia, East Africa, China, and Australia. In North America, delayed WHD values (0.5-4.2 days) were observed under an extremely arid climate measured by PDSI, whereas advanced WHD (7.3-17.6 days) occurred under severe drought conditions measured by SPEI. On the contrary, in South Asia we identified advanced WHD values (7.4-10.6 days) in extremely dry years calculated by PDSI but delayed WHD (2.2-7.6 days) in extremely dry years calculated by SPEI.

Discussion
This study revealed a dominant advancing trend of global WHD during 1981-2014, as observed in other studies based on ground observations and modelling results [6,11,12,38]. Unlike the gradual delay in wheat harvest timing detected in North America from 1981 to 2014, WHDs gradually occurred earlier in other wheat-growing regions, which were closely connected with spatial differences in heat and drought conditions ( Figure S2). The delayed wheat harvest timing in North America may have been caused by a significant increase in moisture ( Figure S2c,d). Many studies based on field experiments and statistical analysis have demonstrated the delaying effects of more water availability on wheat phenology occurrence [18]. Meanwhile, declining NDD in central North America may have strengthened the delay in WHD (Figure S2a), as a larger time span was needed to satisfy the heat requirement for wheat to mature [39]. In most of the eight other wheat-planting areas, however, the increased NDD was primarily responsible for the advance in WHD, especially in Eurasian wheat-planting areas ( Figure S2a). The aggravated drought from 1981-2014 may have further contributed to advanced wheat development in West Asia, central India, and Australia ( Figure S2c,d). Specifically, extreme high temperature events characterized by EDD were more frequent from 1981 to 2014 in northern South Asia and southern Central-West Asia, which may have further forced wheat into earlier senescence beyond the effects of NDD ( Figure S2b).
A significant negative correlation between interannual WHD and NDD variability was observed in almost half of the study region, which indicates that wheat matured earlier in years with higher NDD. Heat condition is the primary factor controlling wheat growth and is usually correlated positively with wheat growth rate. Without considering the effects of other environmental factors, more heat in a previous period facilitates wheat photosynthesis and reduces the time span needed to accumulate enough heat to move wheat growth into the next stage [35]. On the contrary, extreme heat exposure can inhibit wheat growth, increase leaf senescence, and shorten the duration of wheat growth by slowing down photosynthesis (closing the stomata and reducing CO2 input) [40]. However, excessive exposure to extremely high temperatures was only detected in specific regions

Discussion
This study revealed a dominant advancing trend of global WHD during 1981-2014, as observed in other studies based on ground observations and modelling results [6,11,12,38]. Unlike the gradual delay in wheat harvest timing detected in North America from 1981 to 2014, WHDs gradually occurred earlier in other wheat-growing regions, which were closely connected with spatial differences in heat and drought conditions ( Figure S2). The delayed wheat harvest timing in North America may have been caused by a significant increase in moisture ( Figure S2c,d). Many studies based on field experiments and statistical analysis have demonstrated the delaying effects of more water availability on wheat phenology occurrence [18]. Meanwhile, declining NDD in central North America may have strengthened the delay in WHD (Figure S2a), as a larger time span was needed to satisfy the heat requirement for wheat to mature [39]. In most of the eight other wheat-planting areas, however, the increased NDD was primarily responsible for the advance in WHD, especially in Eurasian wheat-planting areas ( Figure S2a). The aggravated drought from 1981-2014 may have further contributed to advanced wheat development in West Asia, central India, and Australia ( Figure S2c,d). Specifically, extreme high temperature events characterized by EDD were more frequent from 1981 to 2014 in northern South Asia and southern Central-West Asia, which may have further forced wheat into earlier senescence beyond the effects of NDD ( Figure S2b). A significant negative correlation between interannual WHD and NDD variability was observed in almost half of the study region, which indicates that wheat matured earlier in years with higher NDD. Heat condition is the primary factor controlling wheat growth and is usually correlated positively with wheat growth rate. Without considering the effects of other environmental factors, more heat in a previous period facilitates wheat photosynthesis and reduces the time span needed to accumulate enough heat to move wheat growth into the next stage [35]. On the contrary, extreme heat exposure can inhibit wheat growth, increase leaf senescence, and shorten the duration of wheat growth by slowing down photosynthesis (closing the stomata and reducing CO 2 input) [40]. However, excessive exposure to extremely high temperatures was only detected in specific regions (mainly South Asia and northern Kazakhstan), so most regions did not suffer severe heatwaves. Nevertheless, the effects of extremely high temperatures on wheat growth were weaker than those of normal temperatures, which differed from the stronger impacts of extremely high temperature stress found in India by Lobell et al. [14]. Increased temperatures, including overall temperature rise and frequent high-temperature events, may also intensify plant and soil evapotranspiration, and thus indirectly influence wheat growth by regulating the water supply. According to the projected results, air temperatures will continue to rise and simultaneously increase the frequency of high temperature events in the future [41]. Therefore, wheat growth will face serious challenges under future climate warming.
Water supply is another critical limitation on plant growth in agricultural systems. The severity of drought in previous periods had significant impacts on wheat harvest timing. When drought severity increased over the prior three months, the wheat maturity date gradually advanced (5.3-13.8 days) in most wheat-growing areas. This was in agreement with results from field experiments carried out in Syria, in which mid-and late-season moisture stress shortened the grain-filling period by 10-11 days [42]. Drought can reduce leaf expansion and longevity, stomatal conductance, and the light-saturated rate of photosynthesis [43]. In addition, transpiration is also significantly reduced under drought stress; this in turn slows leaf heat loss and increases leaf temperature [44]. As a result, increased CO 2 concentrations and photosynthesis can promote wheat development. It is worth noting that the effects of drought stress on wheat harvesting varied across the different developmental stages. In addition, the impacts of extreme drought on WHD were different across the nine wheat-growing areas, which may be closely linked to the difference in drought tolerance of wheat cultivars in different regions [38,45]. As drought is projected to be more frequent and serious in 2080-2099 than 1981-1999 for global wheat-growing areas (especially in North America, Europe, and Australia), more attention should be paid to drought impacts on wheat growth [46].
It should be noted that there are still some uncertainties in this study. First of all, because there were no published annual wheat distribution maps available, we were forced to use the 2005 wheat extent map to determine the range of wheat distribution for all the years. However, in fact, the crop planting system is very complicated. In some places, farmers may grow wheat this year, but switch to another crop next year. For a particular pixel, the percentage of wheat may vary from 1981 to 2014. It is urgent to produce annual wheat distribution maps for a more accurate temporal pattern analysis of wheat phenology. Second, we defined the percentage of wheat grown >30% as wheat pixels. Thus, the estimated wheat harvest timing was inevitably affected by other crops or land covers. Given the relatively coarse spatial resolution of the pixels in this study, only a few pixels have a wheat-planting percentage >50%. Future research is best done on pixels with a higher percentage of wheat grown or even on pure wheat pixels.

Conclusions
In this study, we assessed the impacts of thermal conditions (including normal and extreme conditions) and drought stress on wheat harvest timing with remote sensing data and gridded temperature, and drought data at a global scale. The results demonstrated that increases in normal-growing-degree-days and extreme-growing-degree-days can both significantly advance wheat harvest date. Normal thermal conditions had stronger impacts on wheat harvest timing worldwide, while extremely high temperatures only showed a weaker influence on wheat growth in limited areas. Serious drought could also notably advance the timing of global wheat harvest. Overall, we highlight the advancing effects of heat and drought stress on wheat development across most wheat-planting areas and suggest to focus on strategies to address future climate change impacts on wheat development.

Supplementary Materials:
The following are available online at http://www.mdpi.com/2072-4292/11/8/971/s1, Figure S1. Spatial distribution of multi-year average WHD values estimated by the (a) threshold-based method, (b) logistic-based method, and (c) shape-based method. doy in the color bar denotes the day of year. Figure S2.
Variation trends in (a) NDD and (b) EDD for the three months before WHD; (c) PDSI and (d) SPEI from the second month before WHD to the current month. Figure S3. Frequency of serious drought events. (a)-(c): PDSI; (d)-(f): SPEI. Figure S4. Correlation of spatial pattern of serious drought events measured by PDSI and SPEI for (a) the month when WHD occurs, (b) the month before WHD, and (c) the second month before WHD.
Funding: This study was supported by the National High-Resolution Earth Observation Project (11-Y20A16-9001-17/18) and the National Natural Science Foundation (41771371).