Impacts of climatic and varietal changes on phenology and yield components in rice production in Shonai region of Yamagata Prefecture, Northeast Japan for 36 years

ABSTRACT We investigated temporal changes in paddy rice phenology and yield under climatic and varietal changes in Shonai region, Yamagata Prefecture, northeast Japan across a period of 1982–2017. In the study region, the dominant variety shifted from Sasanishiki to Haenuki in mid 1990s. Temporal trends across the study period were significant in mean air temperatures for the months from May to July, and the temperature rise has accelerated heading, which was delayed by the varietal shift. We built models of the effects of climatic and varietal changes on phenological events and yield components: number of grains per land area, percentage ripened grains, and 1000-grain weight. Combining the yield component models, we built a yield model, with which we estimated the effect of the temperature rise from May to July on the yield. The temperature rise in May and June has reduced the yield by accelerating phenology, while that in July has increased the yield by raising 1000-grain weight. The net effect was no change in the yield across the study period, which, however, does not indicate an insensitivity of the yield to the temperature rise. Rather, the finding points to the dependency of the net change in the yield on seasonal pattern of the climate change. Climate change impacts on rice should therefore be quantified on the changes at monthly or shorter time scale rather than a seasonal mean. Our finding also points to the need to pay attention to the changes in agronomic practices such as varieties.


Introduction
Intergovernmental Panel on Climate Change has projected rise of global mean air-temperature in the range from 1.4 to 5.8°C relative to 1986-2005 period (Pachauri et al., 2014). The rising air temperature has already accelerated plant phenology in long-term observations (Estrella, Sparks & Menzel, 2007;Hu, Weiss, Feng & Baenziger, 2005;Menzel et al. (2006); Williams & Abberton, 2004). Menzel et al. (2006) summarized that 78% of 542 plant species in 21 European countries have had earlier leafing, flowering, and fruiting, while only 3% showed significant delay in the spring phenology. In Japan also, Fujisawa and Kobayashi (2010) reported accelerated bud break and flowering in Apple (Malus pumila var. domestica) as driven by the rising air temperature.
With the accelerated phenology in spring, crop species could reach maturity in a shorter period of crop growth as shown, for example, by Zhang, Huang and Yang (2013), who reported the shortening of growth duration by 4.1-4.4 days°C −1 for rice cultivation in China since 1980s. Rising temperature has also led to shortening of growth duration for winter wheat in North China (Wang, Wang, Feng, Yin & Yu, 2013). Similar report in Australian pome fruit growing regions ranging from 24 to 43 years also mentioned that expected advancement of phenophase ranged from 4.1 to 7.7 days°C −1 (Darbyshire, Webb, Goodwin & Barlow, 2013). Across the shorter growth duration, the crop plants can accumulate less amount of intercepted solar radiation, which could be translated into less biomass accumulation and reduced yield at harvest following the scheme of crop growth based on the radiation-use-efficiency concept (Monteith, 1977;Sinclair & Muchow, 1999). In the prediction of climatic change impacts on crop yield also, accelerated phenology and the resultant shortening of growth duration are among the causes of the yield reduction due to the temperature rise (Horie, Yajima & Nakagawa, 1992;Kang, Khan & Ma, 2009;Kim, Ko, Kang & Tenhunen, 2013;Tubiello, Donatelli, Rosenzweig & Stockle, 2000).
Supplemental data for this article can be accessed here.
We think, however, that the link between the accelerated phenology and reduced crop yield due to temperature rise must be tested rigorously against observations in the real world, where cultivation methods, varieties, and other agronomic practices change over time intervening into the link between climatic variables and crop phenology and yield. For example, Tao, Zhang, Zhang and Rotter (2014) reported growth duration of maize in China has prolonged over time due to combined effects of three main factors: temperature rise, changes in agronomic management, and varietal shift. In Japan, a similar study is yet to be reported, but the varietal shifts for rice in the Northeast Japan (Motoki, 1999) should have changed the relationships between climatic variables and the crop phenology and yield over time. Indeed, the earlier report by Nakagawa & Horie (1996) for 10 major varieties across Japan showed that phenology development was controlled by both climate changes and varietal differences.
In this study, we investigated the responses of paddy rice phenology and yield to the change in weather variables under the co-occurring shift of varieties in a major rice producing region of Northeast Japan across a period of 36 years from 1982. More specifically, we tested the following hypotheses across the study period, 1. climate has changed beyond the normal range of fluctuations, 2. rice plant phenology has changed in response to the climate change, 3. crop yield components have responded to the changes in phenology and climate, and 4. the crop yield has changed due to the changes in the crop yield components.
In testing the above hypotheses, we included the shift of varieties as a co-occurring change during the study period. Because our interest is in the effects of climate change beyond the normal fluctuation, we focused on the climatic trends and their effects on the crop phenology, the yield components, and eventually the yield.

Materials and methods
2.1. The study region, crop performance, and major varieties We studied rice production in Shonai region of Yamagata Prefecture, which is the 4 th largest in rice harvest among the prefectures of Japan as of 2017, when Shonai produced 156, 800 tonnes or 41% of the rice harvest in Yamagata Prefecture. Shonai is located in northeast Japan between latitudes about 38.5°N and 39.0°N along the coast of the Sea of Japan, and consists of 5 townships: Sakata, Tsuruoka, Shonai, Mikawa, and Yuza, as a rice-producing district of the crop statistics.
We obtained the monitoring data on rice yield, yield components, and phenological events: dates of transplanting (DTR), heading (DHD), and maturity (DMT), in Shonai for a 36 years period from 1982 to 2017 from the statistics office in Yamagata Prefecture of the Ministry of Agriculture, Forestry and Fisheries of Japan (MAFF). Since the crop performance records are based on an average of survey data at about 600 locations across the study region, the effects of varietal shifts, if any, are implicit in the temporal change of the crop performance.
To detect the varietal shift in the study region, we retrieved the areas planted to each of the major varieties from the data archive of MAFF (http://www.library. maff.go.jp/archive/Search?m=0&p=1&s=0&lpp=0&cId= 153) for the period from 1982 to 1995, and those from 1996 to 2017 downloaded at the web site for information on rice of Yamagata Prefecture (http://agrin.jp/hp/ kome/library/). Effects of other agronomic practices than the date of transplanting and the variety are implicit in the crop phenology and yield components in this study. Confounding effects of such implicit drivers could be detected, however, by our method as described below to analyze the effects of climatic change on crop phenology and yield components.

Weather variables
We used the observation records of monthly mean air temperature (MAT) and monthly total sunshine hours (TSH) at Sakata Regional Weather Observatory (latitude: 38°54.5ʹ N, longitude: 139°50.6ʹE) of Japan Meteorology Agency (JMA). The data were retrieved from the JMA website (http://www.data.jma.go.jp/obd/stats/etrn/). Recognizing the bias in TSH due to the change of the sensor types in 1986 as noted on the JMA website, we estimated TSH for the period from 1982 to 1985 from incident shortwave radiation (ISR) at the same observatory with linear regression of TSH on ISR for the period from 1986 through to 2009. While we could have used ISR in place of TSH in this study, we used the latter because the observation of ISR at Sakata was discontinued in October 2009. We also estimated a missing data of TSH at Sakata for September of 2010 from TSH at the closest AMeDAS (automated meteorological data acquisition system) Station in Tsuruoka (latitude: 38°44.1ʹ N, longitude: 139°49.7ʹ E) in Yamagata Prefecture using linear regression. Statistical summaries of MAT and TSH are shown in Fig. S1 for each month.

Analyses of the crop and weather data
To describe the shift of varieties, we defined the Top variety as the variety to which the largest area was planted in the study region. Name of the Top variety, e.g. Sasanishiki, was assigned to a nominal variable: NTV. In Shonai, as described later in Results, the shift of the Top variety occurred only once during the 36 years. Accordingly, NTV, in effect, split the entire duration into two periods before and after the varietal shift.
We tested our first hypothesis by regressing the weather variables on year. For robustness of the trend detection, Kendall's rank correlation with year was also calculated. The rank correlation detects monotonic trend without the linearity assumption. For the variables showing significant trend, the observation record was partitioned into two components: Trend across years and Fluctuation from year to year around the Trend. The Trend was estimated by a smoothing spline and the residual from the Trend was regarded as the Fluctuation (Fujisawa & Kobayashi, 2010). For the weather variables with no significant trend, the time series as a whole was regarded as the Fluctuation. We also analyzed the records of crop phenology, yield, and yield components, for their temporal trends in a similar manner as for the weather variables.
To test the second hypothesis that crop phenology has been changed by the climate change and the varietal shift, we fit the linear models for each phenological event (transplanting, heading and maturity) on MAT (mean air temperature), TSH (total sunshine hours) and NTV (name of the Top variety). We chose MAT and TSH for the months preceding the phenological event, e.g. from May to July for heading. Significance of the explanatory variables was judged with step-wise regression, where we chose the model that gives smallest AICc (the Akaike Information Criteria with the correction for the sample size). In some cases, where multiple models gave similar AICc, we chose the simplest one. Taking advantage of the temporal sequence of the phenological events, we included DTR (transplanting date) as an explanatory variable for the model of DHD (heading date), and DHD as an explanatory variable for the model of DMT (maturity date). In fitting the linear models as noted above and below, we also calculated VIF (variance inflation factor) to ensure that the problem of multicollinearity was avoided.
We tested the second hypothesis when a weather variable was chosen in the model of the phenological event, and the same variable showed a significant temporal trend. We maintained the hypothesis when the sensitivity of the phenological event to the Trend of the weather variable is not statistically different from that to the Fluctuation. Otherwise, we judged that changes in other drivers than the weather variable caused the apparent phenological trend. This method assumes that a phenological response to the Fluctuation shows the true sensitivity, but that the response to the Trend could include the effects of other variables that changes across the study period. For rationales and details of this methodology, the reader should refer to Fujisawa and Kobayashi (2010).
The third hypothesis that yield components responded to the changes in phenology and climate was tested in the same way as for the second hypothesis. The only difference was that the yield components were log-transformed in the linear models on MAT, TSH, NTV, and the phenological timings. This is to describe the log-transformed yield as the simple summation of the log-transformed yield components, viz.
ln YLD = ln NGL + ln PRG + ln TGW + ln 10 −7 , [1] where YLD (t ha −1 ) is yield at harvest, NGL (m −2 ) is number of grains per unit land area, PRG (%) is percentage ripened grains, and TGW (g) is 1000-grain weigh. The original record had number of heads per land area and number of grains per head separately, but they were combined as NGL in this study. This is to avoid the often-reported negative correlation between them (e.g. Nemoto, Hamasaki, Matsuba, Hayashi & Yanagihara, 2016), which discredits the idea that they are two separate components. The last term of eqn. 1 is to adjust for the differences among the units of the yield components. Eqn. 1 indicates that the fourth hypothesis can be tested by summing up the third hypotheses for each of the yield components.
We compared the model-estimated yield and yield components with the reported ones according to Kobayashi and Salam (2000). In their methodology, the model performance is measured by the mean squared deviation (MSD) between the model estimates and reported values, and MSD is partitioned into three components: bias of the model estimation (SB), difference in the data dispersion (SDSD), and lack of correlation (LCS).
We used JMP Pro 13.0.0 (SAS Institute, Cary, USA) for the data analyses.

Temporal change in crop varieties
Sasanishiki dominated the planted area in the study region from 1982 to 1992, but thereafter the area sharply declined, while the area of Haenuki started to increase ( Figure 1). Varieties were in transition for three years from 1995 to 1997, when no single variety took more than 50% of paddy area in this region. In 1998, the area of Haenuki rose to 52%, while that for other varieties diminished ( Figure 1). Haenuki stayed at the top of varieties through to 2017. The other major varieties included Domannaka, Hitomebore, Koshihikari, and Tsuyahime, none of which had the largest planted area in this region. We therefore set the nominal variable: NTV (name of Top variety) at 'Sasanishiki' for the period from 1982 to 1994, and 'Haenuki' from 1998 to 2017. We did not use the crop data for the transition period from 1995 to 1997 in the analyses below.

Temporal trends in weather variables
We investigated the temporal changes in weather variables from April to September considering the distribution of the phenological timings in the study region ( Figure 2): transplanting has been done in early-to mid-May, heading occurred in early August, and maturity has been reached between late September and early October. Monthly mean air temperatures from May to July and September showed statistically significant increase with similar slopes in the range from 0.04 to 0.05°C y −1 across the study period (Table 1). The warming trends were significant even without the linearity assumption. In contrast, no significant trend was found in mean air temperature in April and August, and monthly total sunshine hours for most months but September (Table 1).

Temporal trends in crop phenology, yield, and yield components
Among the phenological timings, only transplanting showed significant trend: it has become later ( Table 2). With the date of heading showing no clear trend despite the delaying transplanting, the period from transplanting to heading has become shorter significantly. The period from heading to maturity has not changed across the years (Table 2).
Yield has shown no significant temporal trend, but the yield components have shown the trends. NGL (number of grains per unit land area) has become less, while TGW (1000-grain weight) and PRG (percentage ripened grains) have increased (Table 3).

Effects of climate and variety on crop phenology
DTR (date of transplanting) was best explained by MAT in April and year (Table 4). When the varietal shift was considered by including NTV in place of year, its effect was significant (p = 0.012) but R 2 of that model was much lower (0.466) than that shown in Table 4. Because MAT in April showed no temporal trend (Table 1), the delaying trend in transplanting (Table 2) remained unexplained by the variables considered in this study.
DHD (date of heading) was best explained by NTV, DTR and MAT_May-July (mean air temperature from May to July) (Table 4). When estimating DHD, the coefficient (1.28) for NTV is added for Haenuki and subtracted for Sasanishiki. Haenuki is therefore estimated to reach heading later by 2.6 (2 x 1.28) d than Sasanishiki. It is noteworthy that a delay in transplanting by 1 d delayed DHD by only 0.37 d. As monthly mean air temperatures for the period from May to July increased significantly (Table 1), the mean across the 3 months (MAT_May-Jul) also showed a significant increase (r = 0.641, p < 0.001; τ = 0.448, p < 0.001). To test the hypothesis 3 on DHD: heading has been accelerated by the increasing air temperature, we partitioned MAT_May-July into its Trend and Fluctuation. We then fit a linear model with the two components to DHD minus effects of the other terms: Constant, NTV, and DTR. The close temperature sensitivity of DHD between the Trend and the Fluctuation (Table 4) supported the hypothesis. It can therefore be judged that the apparent lack of temporal trend in DHD (Table 2) could be due to the counteracting effects of the delaying DTR and the rising air temperature.
DMT (date of maturity) was strongly affected by DHD, whose change by 1 d would change DMT by 1 d (Table 4). DMT was also affected by TSH_Aug-Sep (total sunshine hours from August to September), which, however, showed no temporal trend (r = 0.084, p = 0.628).

Effects of climate, variety, and crop phenology on yield components
Log-transformed NGL (number of grains per unit land area) was best explained by the varietal shift (NTV) and sunshine   hours from May to June (TSH_May-Jun) ( Table 5). With no trend in the latter variable (r = 0.207, p = 0.242), the apparent decline in this yield component (Table 3) was accounted for by the varietal shift, which resulted in less grains per land area in the new Top variety, Haenuki, than the former one, Sasanishiki. Log-transformed PRG (percentage ripened grain) was also affected by the varietal shift (NTV), which, however, increased this yield component as opposed to NGL (Table 4). While TSH_Aug (sunshine hours in August) had significant impact on PRG, it showed no temporal trend (Table 1). The apparent increase in PRG (Table 2) can therefore be explained by the varietal shift (NTV). The effect of the varietal shift on this yield component was dependent on sunshine hours in August as shown by the significant interaction between the two variables (Table 5). Indeed, PRG was reduced more in Sasanishiki by the decline of sunshine hours than that in Haenuki (Figure 3).
Log-transformed TGW (1000-grain weight) was best explained by the model with NTV, DMT, and MAT_Jul-Aug (Table 5). This yield component has increased due to the varietal shift, but has also been influenced by the rise of MAT_Jul (Table 1) via MAT_Jul-Aug. To test the hypothesis 3 on the effect of the increasing trend in MAT_Jul on log TGW, the time series of MAT_Jul was partitioned into the Trend and the Fluctuation, and a linear model with these two terms only (no Constant term) was fit to log TGW after all the other effects including that of MAT_Aug being subtracted. The very close sensitivity of log TGW to the two components supports the hypothesis that the rise in MAT_Jul has actually increased TGW. On the other hand, as noted for DHD before, the rise of MAT_Jul, together with that of MAT_May and MAT_June, has accelerated heading, i.e. less DHD, which has led to less DMT and eventually smaller TGW. The rising temperature from May to July had thus counteracting effects on TGW. Net effect of the counteracting effects is evaluated by combining the models of phenology and yield components into the model of the yield as shown below.
of the log yield ( Figure 4). Moreover, the phenological timings, e.g. DHD, are among the explanatory variables in the model of a yield component (TGW) (Table 5 and Figure  4). Since the phenological events are also described by linear models of weather variables and varietal shift (Table   4), they are simply incorporated into the linear models of log yield components (Table 5) and eventually the log yield ( Table 6). The log yield model includes all the monthly mean air temperatures from April to August and total sunshine hours from May to September except that in July in addition to the effect of the varietal change. The significant interaction between NTV and TSH_Aug on PRG is accounted for by assigning different values to the coefficient of TSH_Aug for the two varieties ( Table 6).
As noted above, MAT_Jul has counteracting effects on 1000 grain weight and accordingly the yield (Figure 4). The coefficient for MAT_Jul is positive (0.00532) in the yield model (Table 6), which means that the net effect of the counteracting effects of this weather variable is to increase the yield. This yield increase could, however, be overridden by the yield decrease due to the increase of MAT_May and MAT_Jun, both of which have negative coefficients for the yield ( Table 6). As noted earlier, the hypothesis 3 is supported for each of the counteracting effects: accelerated heading and increased 1000grain weight due to the temperature rise. Since the yield is simple sum of the yield components, the hypothesis 4 should be supported: the yield has changed in response to the climate change particularly the rise of temperature during the May-July period. The direction and extent of the climateinduced yield change, however, hinge on the counteracting effects of the temperature rise as addressed in Discussion.

Testing the models of yield components and yield against the reported values
Since the models of the yield components were fitted to the reported values, there is no wonder that MSD (mean squared deviation) is dominated by LCS (lack of correlation) with negligible SB (bias of estimation) in NGL, TGW, and PRG (columns NGL and PRG in Table 7). The yield components can thus be estimated by the models with ca. 2% relative root mean square errors (columns NGL, PRG, and TGW in Table 7).
Unlike the yield components, the reported yield was not used to make the log yield model, which was simple summation of the yield component models. The comparison between the reported values and model-estimates of yield showed non-negligible SB (column YLD in Table 7), whose contribution to MSD was nevertheless minor. When the comparison between the reported and modeled yields was visualized ( Figure 5(a)), it became obvious that the reported  Tables 4 and 5 are divided by 2 (3) to get the coefficient for the individual months, whereas coefficients for TSH across more than one months were set at the same values in Tables 4 and 5. c. Coefficients for ln TGW are the sum of the direct effects (Table 5) and the indirect effects via phenological timings, e.g. DHD that appears in both Tables 4 and 5, for each explanatory variables. Note that the models for phenology are nested following the sequence of the phenological events, i.e. DMT includes DHD, which includes DTR, among the explanatory variables. d. The coefficient for ln YLD is sum of the coefficients for the yield components except for the Constant, for which see caption g below. e. The effect of TSH_Aug on ln PRG has significant interaction with the effect of NTV (Table 5), and, hence, coefficient for TSH_Aug is shown for the two varieties separately. f. Coefficients for Haenuki are shown, while those for Sasanishiki can be obtained by inverting the sign. It must be noted that the effect of NTV on PRG varies by TSH_Aug due to the interaction between them. It follows that the effect of the varietal change on the yield can only be evaluated with the coefficients for both NTV and TSH_Aug. g. Constant for ln YLD includes the last term (ln 10 −7 ) of Eqn. 1 to account for the discrepancy between the units. yield in 2004 was by far lower than that estimated by the model. The reported yield was below the 99.99% confidence limit of the regression line fitted with no intercept to the values of the other 32 years. The probability of observing such a deviation for a 33 years period is less than 0.5%. When the anomalous year was omitted from the model-report comparison in the yield, MSD was reduced by 34%, and the relative root mean square error shrank from 3.6 to 2.8% (column YLD-2004 in Table 7).
The same anomaly for 2004 was also evident in the comparison between the reported yield and the product of the reported yield components (Figure 5(b)). In theory, the latter 'yield' should be equal to the reported yield, but this is not the case due to the difference in what they actually represent. The reported yield components are essentially based on the field survey conducted by the statistics office, whereas the reported yield is the ratio of the crop harvest amount to harvested area in the region. The reported yield and the yield calculated from the reported yield components are close to each other in most years, but never be identical ( Figure 5(b)). The anomaly in 2004 is also evident in Figure 5(b), which partly accounts for the large deviation of the modeled yield from the reported one ( Figure 5(a)). The cause of the anomalously low yield in 2004 is sought in Discussion.

Overall effect of rising air temperature on rice yield
As noted above, the increase in mean air temperature in July (Table 1) has had counteracting effects on 1000grain weight: it reduced TGW via accelerated maturity, but it also increased TGW directly (Figure 4). The coefficient (0.00554) for net effect of MAT_Jul is given by the difference between the coefficient for MAT_Aug (0.00937) and that for MAT_May (as well as MAT_Jun)   (−0.00383) ( Table 6). It indicates that the positive effects overridden the negative one, and that the temporal increase in air temperature in July had increased the yield. On the other hand, the increase in air temperatures in May and June (Table 1) should have reduced the yield with the negative coefficient (−0.00383) on the yield (Table 6). We can calculate the overall effect of the temperature rise for the period from May to July on the yield as follows.
Let us first estimate the negative impact of increasing trend of temperature from May to June. Assuming the linear trends from May to June (Table 1) for 36 years multiplied by the coefficient (−0.00411) for ln YLD (Table 6) Sum of the above changes to opposite directions represents net effect of the temperature increase for the period from May to July: a decrease of yield by only 0.2% across the 36 years. This is practically naught, and the hypothesis 4: rice yield has changed due to the climate change, can be rejected in a quantitative sense. This does not mean insensitivity of the yield to temperature rise, however. The nearzero change is a result of the seasonal variability of the temperature rise. The temperature rise was not detected for August, but if it were, the positive effect should have prevailed. On the contrary, if the temperature rise was not significant in July as reported for many locations in northeast Japan (Shimono, 2008), the negative effect should have prevailed. We may not therefore expect that the climate change impacts on rice yield can be quantified on a seasonal mean basis not just because of the climatic injuries on a shorter time scale, e.g. chilling damages to developing panicles, but also due to the counteracting effects as found in this study.
In the yield model (Table 6), the effect of year is also included to account for the temporal trend in transplanting time (Table 4). The effect of the delaying trend in transplanting on the yield across the 36 years is estimated to be only 0.3%, however, viz.
The effect of the unknown driver of the delaying transplanting should therefore be negligible.

Impacts of the varietal shift
As we analyzed the data reported by the statistics office across 36 years, the change of varieties was inevitable but presented us a challenge in the analyzing the date for climatic change impacts. Motoki (1999) has previously reported the shift of varieties from Kiyonishiki to Sasanishiki and to Haenuki in Yamagata Prefecture for the period from 1970 to 1995. The varietal shift took place mainly for the sake of better quality (Motoki, 1999) or greater field size (Saito, 2007). Nihei (2010) mentioned the repealing of the Food Control Law in 1995 as a driver of the changes in rice varieties across northeast Japan. The law was enacted due to the widespread crop failures caused by the cold summer of 1993. In Shonai region, Haenuki was widely adopted for its good grain quality along with higher resistance against diseases and lodging due to its shorter stature. Our analysis on the effects of the varietal shift estimated that Haenuki had heading date by about 2.5 d later (Table 4), smaller number of grains per land area, but higher percentage ripened grains and greater 1000-grain weight (Table 5) as compared with its predecessor, Sasanishiki. These varietal features are very close to those presented by Satoh et al. (1992) in their Figure 3 and Table 8 for Yamagata No. 45, which was later named as Haenuki. The percentage ripened grains in Haenuki is higher than Sasanishiki particularly under lower sunshine hours in August (Figure 3).

Climatic effect on crop phenology and yield components
The impact of climate change on rice phenology in this study is in line with the past studies (Estrella et al., 2007;Hu et al., 2005;Menzel et al., 2006;Williams & Abberton, 2004). The temperature sensitivity (5.11 d°C −1 ) of heading date (Table 4) is only slightly higher than that for rice in China (Zhang et al., 2013) and within the range of interspecies variability (e.g. Menzel et al., 2006). In most plant species, vegetative development usually has a higher optimum temperature than that for reproductive development (Hatfield & Prueger, 2015), and the duration from transplanting to maturity is shortened (Table 2). Our result thus confirmed the findings of Zhang et al. (2013), who showed that increasing temperature shorten rice growth duration in China in the last three decades (1980s-2000s).
The effects of weather variables on yield components has been studied for many years, but its use to detect the signal of the climate change impacts on crop production is relatively recent. The study by Nemoto et al. (2016) identified the effects of weather variables on yield components in rice variety Kirara 397 in the fields mostly in Hokkaido in northern Japan for the period from 1989 to 2013. Their findings have some commonalities with those in this study (Figure 4), such as the increase of 1000-grain weight due to higher temperature for the period from panicle initiation to heading, and the increase of percentage ripened grains as affected by solar radiation during the grain filling period. The positive effect of temperature on 1000-grain weight is of particular interest, since it is responsible for the counteracting effects of rising air temperature on the yield. Fukushima, Ohta, Kaji and Tsuda (2015) also found a positive correlation between monthly mean temperature in July and 1000-grain weight in an analysis of field performance of rice varieties from 2000 to 2013 at an experimental station in Akita Prefecture, which adjoins the study region. It might appear counter-intuitive that temperature before heading affects 1000-grain weight at harvest, but this effect has been demonstrated by Matsushima, Tanaka and Hoshino (1964), who subjected potted rice plants to different water and air temperatures at different growth stages. They reported that 1000-grain weight increased by raising water temperature from 21 to 31°C at late panicle development stage. They attributed this effect to enhancement of spikelet growth and increased hull size by the higher temperatures during the panicle development stage. It has indeed been shown that hull size is a determinant of individual grain mass (Satoh, 1968).
The rising trend in 1000-grain weight has been reported by Shimono (2008), who analyzed the temporal trends in phenology, yield components, and yield of rice variety Sasanishiki at agricultural experiment stations across northern Japan for the period from 1980 to 2005. The significant trend of increasing 1000grain weight was the most salient among the temporal trends detected in his study, and indeed mostly responsible to the significant increase in the grain yield. No trend was detected for the number of panicles per land, and the increasing trend in the length of a panicle is much weaker than that for 1000-grain weight (Shimono, 2008).
To our knowledge, the climate change impact studies on rice have not paid much attention to the effect of temperature on 1000-grain weight, and accordingly this effect has not been addressed in crop growth models for the impact prediction. One of the requirements for the impact estimation by using a process-based model is that the model must describe the individual processes explicitly. This is indeed strength of the approach, but could be a weakness as well, because it cannot address an unknown process without sacrificing the strength.
The yield component approach could augment the process-based approach by offering a complete coverage of the yield via its components. The components include various growth processes only implicitly, but they could identify a need to look at growth processes that have been studied less than they should have been.

The anomaly in 2004
In 2004, for which we detected the anomaly in reported rice yield ( Figure 5), they had a record low crop index of 87 in Shonai region due to damages by a Typhoon, which traveled over the Sea of Japan from the southwest to the northeast and came close to the coast on August 20 th . The rice plants were from 1 to 3 weeks after the heading stage, and vulnerable to the saline wind brought by the Typhoon (Fujii et al., 2006). It was later reported that the main cause of the damage was the reduction of percentage ripened grains caused by the damages by the saline water on rachis branches (Mori, Matuda, Shibata & Fujii, 2008). The reported percentage ripened grains showed a dip to 86% in that year from c. 90% in other years, but it underestimated the region-wide damages on the yield component in the range from 80 to 40% for Haenuki (Mori et al., 2008).
The Typhoon damage in 2004 presents us a challenge that a climatic anomaly could damage the regional rice production to the unprecedented and unexpected degree. Unpredictable nature of such a climatic event would make the quick and accurate evaluation of the damages across the region one of the strategies for greater resilience against the climatic hazards.

Limitations of this study and required developments
This study found the effects of the varietal change on a wide range of traits in rice phenology and yield components. The interaction of the varietal shift and sunshine hours in August (Figure 3) was particularly revealing. While the weather variable did not show a temporal trend, the shift of varieties could introduce temporal trend into the sensitivity of rice production to that weather variable. Admittedly, our way of accommodating the varietal shift is too simplistic and requires elaborations in future studies. It would be challenging, however, to work on a region (or regions) with more frequent shifts of varieties or more major varieties planted in overlapping periods of time.
Further elaborations could also be made on the temporal resolution in aggregating the weather variables. Our models of crop phenology, yield components, and yield are based on monthly mean and monthly sum of the weather variables (Table 4-6, and Figure 4). This is because we aimed in this study to quantify the impacts of changes in monthly climate across the study period on rice yield. It must be noted, however, that some climatic anomalies, e.g. chilling injury to pollen development, exert impacts at a shorter time period, and that quantifying such impacts under the changing climate require models with a finer temporal resolution.
Despite these challenges, we believe, improved understandings of the climate change impacts in the real-world situations is crucial to better adapt rice production to the changes in climate and other drivers in the decades to come. Regionally-based studies like this could make their own contributions to the science of climate change.