Response of rice phenology to climate warming weakened across China during 1981–2018: did climatic or anthropogenic factors play a role?

Climate warming has substantially shifted plant phenology, which alters the length of growing season and consequently affects plant productivity. Recent studies showed a stalled or reversed impact of climate change on vegetation phenology since 1998, as well as an asymmetric warming effect. However, how field crop phenology responded to the recent climate warming and the asymmetric warming remains unknown. In addition, the relative roles of climate change, sowing date and cultivars shifts in the spatiotemporal changes of crop phenology at different regions need to be better understood. Here, using the latest 9,393 phenological records at 249 agro-meteorological stations across China over 1981–2018, we critically investigated the spatiotemporal dynamics of rice phenology and disentangled the effects of different drivers by exploiting the physiological relationship between crop phenology and thermal accumulation. The results showed that length of growing period (GP) increased by 3.24 ± 0.15 days/decade for single rice, 1.90 ± 0.22 days/decade for early rice and 0.47 ± 0.14 days/decade for late rice. Although climate warming during rice GP did not slow down, the trends in rice GP and the correlations between GP and temperature decreased generally from 1981–1999 to 2000–2018. The weakened phenological response to climate change was mainly caused by agronomic managements, especially cultivar shifts. Climate warming shortened GP by 0.84 ± 1.80, 1.23 ± 0.77, and 1.29 ± 1.24 days/decade for single rice, early rice and late rice, respectively. However, cultivar shifts prolonged it respectively by 3.28 ± 3.68, 2.15 ± 2.38, and 2.31 ± 3.36 days/decade, totally offsetting the negative effects of climate warming. Rice responded to daytime and night-time warming differently with night-time temperature affecting GPs more. Our study provided new insights that rice phenology responded to night-time warming more than daytime warming across China however the response to climate warming weakened, and cultivar shifts outweighed climate change in affecting rice phenology.


Introduction
Phenological stages, which are strongly controlled by short-and long-term variability in climate, have been the first-observed biological footprint of climate change impact , Ettinger et al 2020, Wu et al 2021. Global warming has greatly shifted the timing of plant phenological events, altered the length of the growing season, and consequently affected ecosystem structure, function and their climatic feedbacks (Liu et al 2016, Piao et al 2017, Vitasse et al 2018.
Crop phenology has great ecological and economic implications, considering its great biophysical influences on land surface heat and water fluxes, the dominant harbor of biodiversity and the close relation to crop productivity (Challinor et al 2014, Zhao et al 2017, Asseng et al 2018. Crop systems are valuable natural laboratories for exploring ecosystem's response and adaptation to climate change, especially for such complicated and highly heterogeneous environments in China. Comparing with more concerns on forest or grass, and consistent results on natural ecosystems (Vitasse et al 2018, Wu et al 2018, 2021, Meng et al 2020, Ettinger et al 2020, Rosbakh et al 2021, crop phenology is more complex because it is subject to both natural and anthropogenic factors (Tao et al 2012, Abbas et al 2017, Hu et al 2017. For example, climate warming can advance crop phenology and shorten growth duration, but agronomic management such as sowing date and cultivar shifts may offset the negative impacts to some extent (Tao et al 2012. In addition, crop phenology is also complicatedly affected by many other factors such as photoperiod and vernalization (Tao et al 2012, Zhang et al 2014a, He et al 2015, Rezaei et al 2017, Wang et al 2017a, Ortiz-Bobea et al 2021. Attributing the changes in the phenology is thus exceptionally challenging. Nevertheless, understanding the phenological dynamic and its underlying drivers is essential for crop management, crop growth monitoring and production prediction, as well as developing climate change adaptation options (Tao et al 2006, Siebert and Ewert 2012, Sánchez et al 2014, Shew et al 2020.
Recently, some studies showed a stalled or reversed impact of climate warming on plant phenology since 1998 (Fu et al 2015, Ballantyne et al 2017, Güsewell et al 2017, Piao et al 2017, Wang et al 2019c. For example, Fu et al (2015) showed that the sensitivity of spring leaf unfolding to climate warming decreased by 40% from 1981-1994 to 1999-2013 using long-term observations for seven dominant European tree species. Wang et al (2019c) showed no trends in spring and autumn phenology globally during 1998-2012. Empirical evidence for this, however, is limited to natural vegetation. Whether climate warming during crop GP slows down, and how field crop phenology responds to such warming remains unclear. In addition, many studies have evaluated crop phenological response to temperature, but few of them explicitly have investigated which temperature variable, i.e. daily mean (T mean ), night-time (T min ) or daytime (T max ) temperature, play a dominant role (Tao et al 2013, Piao et al 2017, Wu et al 2018, 2021. The increase rate of T min over the past five decades was 1.4 times that of T max . This asymmetric warming was reported to have contrasting effects on the timing of autumn plant phenology in the Northern Hemisphere (Peng et al 2013, Wu et al 2018.
Yet how field crops phenology responded to the asymmetric warming remains unclear. Moreover, most of previous studies have focused on estimating the effects of climate variables (He et al 2015, Gao et al 2018, but few studies have quantified the contributions of agronomic practices to the phenological shifts in specific regions (Abbas et al 2017, Hu et al 2017, Rezaei et al 2017, Wang et al 2017a. The few available analyses have rarely taken the mechanical processes into account and not thoroughly investigated the effects of different drivers such as climate, sowing date, cultivar thermal requirement and photoperiod sensitivity. The analyses focused on the influences of anthropogenic managements have mostly conducted on the observations before 2013 (Tao et al 2006, Siebert and Ewert 2012, Zhang et al 2014a, Parent et al 2018. Crop phenological responses to recent climate warming and its regional disparity have not yet been well documented, especially in China.
Previous studies have mainly used statistical or process-based models to characterize sensitivity of phenology to climate warming and to isolate the contribution of different drivers , Wang et al 2017a. The results on the impacts of climate change on crop phenology were inconsistent, which showed climate change shortened (Tao et al 2012, Abbas et al 2017 or prolonged (Zhang et al 2014a, Hu et al 2017 crop growth period (GP) or had insignificant effects on it (Wang et al 2017a, Tariq et al 2018, Ye et al 2019. Such divergent results may be ascribed to the limited ground observations with coarse spatiotemporal resolution or different analysis methods. With few ground observations, it is difficult to separate the impacts of climate change and agronomic practices on crop phenology by statistical methods , de Los Campos et al 2020. Alternatively, crop model can be applied to investigate the mechanisms on and quantify the contributions of different drivers to crop phenology changes (Jones et al 2003, Müller et al 2017, Xin and Tao 2019, Zhang et al 2020. Yet calibrating crop model with limited field trials and then applying them for the studies over long-term or broad regions will inevitably cause large uncertainties (Asseng et al 2013, Li et al 2015, Rezaei et al 2018. In complement, integrating the two methods on the basis of accumulated thermal development unit, a well-known horticultural principle, provides a promising way to obtain more robust estimations (Lobell and Asseng 2017, Laskin et al 2019, Peng et al 2020. Along this line, in this study, we collected the latest rice phenology dataset including 9,393 phenological records at 249 agro-meteorological experiment stations across China over 1981-2018. Based on the largest dataset on rice phenology so far, we integrated the statistical methods and crop model simulations to (a) investigate the spatiotemporal dynamics of rice phenology over the past four decades, as well as the differences between 1981-1999 and 2000-2018; (b) identify the key temperature variables controlling rice phenology using the updated climate observations at 2459 meteorological stations; (c) disentangle the roles of climatic and anthropogenic factors in rice phenology change.

Rice phenology and climate observation data
Data on the observation of transplanting, heading and maturity date from 1981 to 2018 were obtained from the national agro-meteorological experiment stations across the entire rice cultivation areas in China. The phenological dates were observed and recorded by well-trained agricultural technicians in the experimental field and then checked and managed by the Chinese Agricultural Meteorological Monitoring System. It was the latest dataset with the longest time series and the most observational stations for rice phenological events in China. Only the stations with more than 15 years of records were used for analysis, with 9393 phenological records at 249 observational stations. Rice cultivation areas were classified into five agricultural ecosystem zones (AEZs), according to climate, topography, hydrography, soil, geology and cropping systems (figure 1). Details on the cropping system, number of stations and mean date in each AEZ are shown in table 1 and supplementary figure S1 (available online at stacks. iop.org/ERL/17/064029/mmedia). The daily weather data on T mean , T max and T min were obtained from 2459 meteorological stations across China, which were obtained from the National Meteorological Data Center (http://data.cma.cn/).

Calculating temperature and ADTU during rice GP
Rice GP is defined as the duration from transplanting to maturity, which is divided into the vegetative growing period (VGP, from transplanting to heading) and the reproductive growing period (RGP, from heading to maturity). Average temperature during VGP, RGP and GP was calculated at each station in a year.
Accumulated thermal development unit (ATDU) refers to the photo-thermal requirement of a certain cultivar to reach a specific development stage, which represents the comprehensive characteristic of a cultivar including thermal requirements and photoperiod sensitivity. Trend in ATDU thus represents the effects of cultivar shifts on the phenology. ATDU is calculated as the methods used by cropping system model-CERES-Rice in Decision Support System for Agrotechnology Transfer 4.7 (Jones et al 2003, Hoogenboom et al 2019. ATDU is the sum of daily thermal development unit during VGP and the sum of daily thermal time (DTT) during RGP where T min and T max are daily minimum and maximum temperature, T base is the lower limit and T high is the upper limit at which rice suffered tissue injuries. T opt is the optimal temperature giving the highest rate of rice physiological process. T h is hourly temperature and T he is effective temperature in h hour in a day. DF is a photoperiod factor, P2R is photoperiod sensitivity parameter, P2O is the critical day length. DL is mean day length and estimated as a function of latitude and day of the year (DOY) at a station using the R package geosphere. The parameters in the above formulas (table 1) were obtained from Zhang et al (2014a).

Disentangling the effects of climate change and agricultural managements on rice phenology
We first defined the average of ATDU and transplanting date from 1981 to 1983 as referenced cultivar (ATDU ref ) and transplant date (transplanting ref ), respectively. We then estimated heading and maturity  The residuals between the trends in the observed GP and GP E1 , GP E2 and GP E3 were ascribed to the role of other factors such as fertilization and irrigation. In addition, the effects of all agricultural managements together on rice phenology were estimated as the differences between the trends in the observed GP and GP E1 .

Statistical analysis
We first performed piecewise regression with two linear segments to detect the potential turning point in the time series of temperature (T mean , T max and T min ), phenological events (transplanting, heading and maturity) and GPs (VGP, RGP and GP).  An independent-sample t-test was tested the difference between the two time periods. Trends in the phenological events, GPs and temperatures were examined using linear regression. Partial correlation analyses were conducted to determine the relationship between rice GPs and three temperature variables. Two-tailed t-test was used to test the significant level. Statistical analyses were implemented using R 4.0.3 and Stata 13.

Spatiotemporal changes of rice phenology
In the past 38 years, rice phenological dates advanced at more than half of the stations, although there were remarkable differences among rice types, AEZs and time periods (figure 3). Transplanting date advanced significantly at 50% of the stations for early rice, followed by single rice (43%) and late rice ( 3(g)). The trend for late rice was similar to that for early rice, but with relatively great degree (0.81 ± 0.28 days/decade) ( figure 3(h)). In contrast, maturity dates of all the three rice types consistently delayed across AEZs during the two periods (figures 3(i)-(l)), especially for early rice (the average of 2.23 days/decade), followed by single rice (1.54 days/decade) and late rice (1.18 days/decade).

Changes in length of rice GP and their relations to different temperature variables
Rice GP lengthened at most stations however depending on rice types and studied periods. The trends decreased generally from 1981-1999 to 2000-2018 (figure 4). Rice GP increased significantly at 56%, 46% and 24% stations for single, early and late rice, respectively (figure 4(a)). Single rice GP increased across AEZs during the two time periods (figure 4(b)), especially for zone II (>4 days/decade). Early rice GP behaved similarly as single rice, but with a smaller magnitude (figure 4(c)). In contrast, late rice GP shortened during 1981-1999 whereas lengthened during 2000-2018 ( figure 4(d)). Rice VGP significantly lengthened at more than 30% stations (figure 4(e)). Single rice VGP lengthened in the three AEZs in southern China during the two time periods (figure 4(f)). The trends reversed totally from 1981-1999 to 2000-2018 for early and late rice (figures 4(g) and (h)). Rice RGP significantly lengthened at 43% stations for single rice, followed by early (26%) and late rice (12%) (figure 4(i)). The RGP lengthened in almost all AEZs during the two periods for all the three rice types (figures 4(j)-(l)). Under the intensified climate warming in the past four decades, T max increased more than T mean and T min for single rice whereas T min increased more for early and late rice (supplementary figure S3). Rice GP was generally affected more by T min than T max and T mean for all the rice types (figure 5). The GP was negatively correlated with T min at 89% stations (figure 5(a)), with the highest correlation coefficient for single rice, followed by early and late rice (figures 5(b)-(d)). T max was also negatively correlated with the GP in most AEZs (figures 5(e)-(h)). However, it was positively correlated with single rice GP in zone II and III during 1981-1999 and zone IV during 2000-2018. In contrast, T mean affected rice GP less independently of rice types, AEZs and time periods (figures 5(i)-(l)). Collectively, climate warming during rice GP did not slow down significantly from 1981-1999 to 2000-2018 (supplementary figures S2 and S3); however both the trends in duration and the correlation coefficients between duration and temperature were generally decreased from 1981-1999 to 2000-2018, indicating that responses of rice phenology to climate warming weakened.

Effects of climate warming and agricultural managements on rice GPs
Once disentangling the contributions of different drivers, we found that climate warming negatively, cultivar shifts positively, and transplanting date changes and the other factors such as fertilization and irrigation weakly affected the VGP, RGP and GP, although the impacts varied by rice types and studied periods (figure 6). For single rice, increased temperature reduced GP by ∼0.57 days/decade whereas cultivar shifts increased GP by more than 2.53 days/decade during 1981-2018 (figures 6(a), (d) and (g)). Transplanting date changed GP by −0.15 and 0.10 days/decade during 1981-1999 and 2000-2018, respectively. For early rice, temperature decreased VGP by 1.51 days/decade whereas increased RGP by 0.57 days/decade, leading to a decrease of GP by 0.94 days/decade during 1981-1999 (figures 6(b), (e) and (h)). During 2000-2018, increased temperature negatively affected both VGP and RGP, leading to a decrease in GP by 0.63 days/decade. Similarly, for single rice, cultivar shifts increased GP by more than 1.5 days/decade during the two periods, and transplanting date and other factors had small impacts on GP by <0.5 days/decade. For late rice, GP decreased by 0.58 days/decade during 1981-1999 because of the integrated impacts of cultivar shifts which shortened VGP by −2.24 days/decade whereas lengthened RGP by 1.23 days/decade, and climate warming which lengthened RGP by 0.45 days/decade (figures 6(c), (f) and (i)). During 2000-2018, however, the observed GP was increased by 1.31 days/decade owing to the integrated impacts of cultivar shifts (1.96 days/decade), climate warming (−0.71 days/decade), and shifts in transplanting dates and other agronomic managements. Overall, the growth cycle was prolonged for all the three rice types during both the periods, despite of the negative effects of warming. This was mainly because agronomic managements, especially cultivar shifts, outweigh climate change in affecting the GP changes ( figure 7).

The declined rice phenological response to climate change
We showed that climate warming did not slow down significantly during rice GP (supplementary figure  S3) but the phenological responses declined generally from 1981-1999 to 2000-2018 (figures 3-5). The decline was largely because agronomic managements, especially cultivar shifts, dominantly controlled the GP trends (Tao et al 2013, Zhang et al 2014a, Wang et al 2017a, Ye et al 2019. Climate warming generally advanced the phenology and shortened the GP, whereas cultivar shifts compensated or reversed the negative effects ( figure 6). Yet this compensation effect was reduced from 1981-1999 to 2000-2018 because the rate of cultivar Besides anthropogenic factors, other mechanisms may also play a role, such as the 'photoperiod limitation' mechanisms (Fu et al 2015, Meng et al 2020. This was demonstrated by the negative correlation between day length and ATDU requirements for all the rice types (figures 8(a)-(c)), suggesting that the decreased day length (supplementary figure S4) may decelerate rice development rate thus limit the phenological advancement when the events occur too early in the season. This study reveals that rice phenological response to climate change declined across China, suggesting that rice phenology may continue to advance with climate warming but the rate may slow down. The results show that anthropogenic managements outweigh climate warming in affecting rice phenology, highlighting that climate change impact studies should adequately account for anthropogenic adaptation.

The different effects of daily mean, daytime, and night-time warming on rice phenology
Our analyses showed that rice phenology responded differently to T min , T max and T mean , and T min affected more significantly the GPs for all the three rice types, which was well supported by chambers and field-based studies in cereals across major cropping regions of the world (Peng et al 2004, Impa et al 2021, Schaarschmidt et al 2021, Sakai et al 2022. The stronger negative effects of T min may mostly result from its greater increasing rate (supplementary figure S3). According to equation (1), the low  temperature could become above the base temperature with climate warming and then start to affect rice development and phenology. In addition, T min is negatively correlated with standardized precipitation evapotranspiration index (SPEI) in dry temperate regions (Peng et al 2013, Wu et al 2018, a stronger drought (i.e. a low SPEI) may affect rice growth negatively and consequently lead to an earlier phenology event. A lower pollen viability and pollen germination percentage caused by night-time warming was noticed in rice (Jain et al 2007, Mohammed andTarpley 2009, Zhang et al 2018). An additional mechanism to interpret the negative effects of T min may be that T min is generally positively correlated with autotrophic respiration and evapotranspiration (Wang et al 2017b, Su et al 2018, Lian et al 2020. That is, elevated T min could enhance autotrophic respiration, reduce soil moisture availability and indirectly limit the duration of photosynthesis during the following daytime. Moreover, increased nightrespiration and reduced photosynthesis by higher night-temperature is significantly higher during the post-flowering stage comparing with pre-flowering stage in rice, which will negatively affect grain-filling duration, post-flowering senescence and grain protein composition (Peraudeau et al 2015, Bahuguna et al 2017, Impa et al 2021. That is why RGP is more sensitive to night-time warming than VGP (supplementary figures S4 and S5).
For the advanced phenology and shortened GPs associated with T max , it is mainly because rising T max could reduce photosynthetic activity through enhancing evaporation and reducing soil water content (Schlaepfer et al 2017, Reich et al 2018, Samaniego et al 2018. There is observational evidence that T max is negatively correlated with soil water content retrieved from multiple microwave satellite sensors in dry temperate regions (Owe et al 2008). Furthermore, a higher T max is always associated with stronger radiation and potentially a higher chance of water stress (Piao et al 2017, Wu et al 2018, 2021. By contrast, T mean affected the phenology less and even had a slightly positive correlation with the GPs (figure 5), implying rice growth is more vulnerable to extreme temperatures, especially night-time warming. Therefore, the effects of T min on crop phenology need to be paid more attention to.

The roles of climatic or anthropogenic factors in affecting rice phenology
After disentangling the effects of climatic or anthropogenic factors, we found the factors controlling crop phenology varied by growth stages and studied periods. In 1981-1999, changes in temperature shortened GP for early rice, prolonged GP for late rice, and impacted less on GP for single rice, which were consistent with many previous studies using statistical approaches (Tao et al 2013, Zhang et al 2014a, Hu et al 2017, Ye et al 2019 and crop model simulations (Zhang et al 2016a, Wang et al 2017a. In 2000-2018, however, increase in temperature consistently shortened rice GPs, implying a totally different adaptation strategy (figure 6). During 1981-1999, the cultivars with a longer GP (higher ATDU requirements) were adopted for single rice (supplementary figure S5) to take advantage of the ameliorated photo-thermal resources for higher yields (Tao et al 2013, Zhang et al 2014a. The cultivars with a longer GP by shortening VGP and lengthening RGP were adopted for early rice to improve yield meanwhile leave more days for transplanting late rice earlier (Wang et al 2017a, Ye et al 2019. The cultivars with a shorter GP were adopted for late rice to escape the cold stress during heading and grain-filling stages (Zhang et al 2016b, Wang et al 2016, 2019a. Note that, the cultivars with a longer GP by lengthening both VGP and RGP were adopted for all the three rice types during 2000-2018 except single rice in zone III. The change in VGP could be ascribed to the combined effects of temperature and photoperiod. We showed that ATDU requirements were positively correlated with T min for most rice types and AEZs (figures 8(d)-(f)). Rice yield has been limited by insufficient photo-thermal resources for decades, especially in southern China. With climate change, increase in solar radiation and temperature advanced transplanting date and delayed heading date, lengthening VGP (increasing ATDU requirements) to make full use of the improved heat resources. Interestingly, a negative correlation between ATDU and temperature was found for late rice in zone VI (figure 8(f)), the warmest zone in southern China. Late rice in low latitudes in China was often damaged by extreme high-temperature events (Zhang et al 2014b, Wang et al 2019b. Temperature, above T opt , contributed only partly even negatively to ATDU according to equation (2). The negative correlation suggested that intensive warming surpassed current cultivars' photo-thermal requirements and hence threatened rice development in the region and probably exacerbates with climate warming. In addition, day length affected negatively ATDU (figures 8(d)-(f)), suggesting that decreased day length (supplementary figure S4) may decelerate rice development rate, increase ATDU requirements, and thus prolong VGP. The increased RGP can be largely attributed to cultivar turnover (Tao et al 2013, Zhang et al 2014a, Wang et al 2017a, Ye et al 2019. The cultivars with a high ATDU requirement delayed maturity date and consequently lengthened the growth duration (supplementary figure S5). Yet single rice RGP in zone III was strongly decreased because heading date delayed significantly. In the zone, rice heading dates were around 15 July-10 August , coinciding with the occurrence of extreme high temperature events caused by subtropical high. Cultivars with a later heading date were adopted by farmers to avoid the heat stress during the heading-flowering stage (Tao et al 2013, Wang et al 2014, Zhang et al 2014a.
Transplanting is a field management activity, which is directly decided by farmers, but it could be largely influenced by many factors, including local climate, suggestions from agricultural technicians for a specific cultivar and others such as soil physical problems with droughts or water-logging (Waha et al 2013, Ding et al 2020. There was a notion that crop establishment was not affected much more with sowing change by one or two weeks because of increased temperature and good soil-moisture condition (Tao et al 2012, Zhang et al 2014a, Gao et al 2021. Our study showed, however, delaying transplanting decreased the length of RGP for late rice by 2 days/decade during 2000-2018, which was comparable with climate impacts and exhibited greater variability (figure 6(f)). Field experiments also proved that adjusting sowing date greatly influenced crop growth and development rate, water use efficiency, source-sink relationship as well as crop resilience to climate (Bonelli et

Uncertainties in the study
This study, like many others, also had some uncertainties. First, the regional comparisons could cause some uncertainties due to the uneven distribution of agro-ecological stations. Second, a single temperature threshold over rice development phases might overestimate the contribution of cultivar shifts to rice phenology change. A more detailed model instead of the phenology model used in our study could improve the reliability and robustness of estimations. Third, statistical analysis had inherent uncertainties because the output was affected by sample size and method. Despite the limitations, our study provided an observational evidence for a weakened rice phenological response to climate change using the latest phenological and meteorological observations at a larger number of field stations with a longer period than any previous study on crop phenology (Tao et al 2012

Conclusions
In this study, we used 9393 phenological records at 249 agro-meteorological stations and climate observations at 2400 meteorological stations across China to investigate the spatiotemporal changes of rice phenology during 1981-2018, as well as their relations to temperature change, transplanting date and cultivar shifts. The larger and valuable dataset in both space and time allows us to investigate the response of rice phenology to climate change more critically and reliably. We found rice phenological response to climate warming declined from 1981-1999 to 2000-2018, but climate warming during rice GP did not slow down. Rice GP was generally affected more by T min than T max and T mean for all the rice types. Increased temperature significantly shorten the GP, especially for early rice, followed by single rice and late rice. Cultivar shifts offset totally the effects of climate warming on GP change for all the three rice types. Changes in transplanting date and other agronomic managements had relatively small effect on rice GPs.

Data availability statement
All data that support the findings of this study are included within the article (and any supplementary files).