Evaluating irrigated rice yields in Japan within the Climate Zonation Scheme of the Global Yield Gap Atlas

Abstract The Global Yield Gap Atlas (GYGA) is an international project that addresses global food production capacity in the form of yield gaps (Yg). The GYGA project is unique in employing its original Climate Zonation Scheme (CZS) composed of three indexed factors, i.e. Growing Degree Days (GDD) related to temperature, Aridity Index (AI) related to available water and Temperature Seasonality (TS) related to annual temperature range, creating 300 Climate Zones (CZs) theoretically across the globe. In the present study, the GYGA CZs were identified for Japan on a municipality basis and analysis of variance (ANOVA) was performed on irrigated rice yield data sets, equating to actual yields (Ya) in the GYGA context, from long-term government statistics. The ANOVA was conducted for the data sets over two decades between 1994 and 2016 by assigning the GDD score of 6 levels and the TS score of 2 levels as fixed factors. Significant interactions with respect to Ya were observed between GDD score and TS score for 13 years out of 21 years implying the existence of favourable combinations of the GDD score and the TS score for rice cultivation. The implication was also supported by the observation with Yg. The lower values of coefficient of variance obtained from the CZs characterized by medium GDD scores indicated the stability over time of rice yields in these areas. These findings suggest a possibility that the GYGA-CZS can be recognized as a tool suitable to identify favourable CZs for growing crops.


Introduction
The year 2020 was experienced as an unprecedented year all over the globe. As a response to the Covid-19 pandemic declared by the World Health Organization on 11th March, food export restrictions were imposed by several nations (Hepburn et al., 2020). In addition, reports of flooding along the Yangtze river (Myers, 2020) and locust attacks in various parts of the world (FAO, 2020) led to concerns about food supply and attention was focused on food security in particular in nations like Japan that are heavily dependent on food imports. An international project named Global Yield Gap Atlas (GYGA) initiated by  has been helping to provide the international community with the basic knowledge to address the issue of food security by adopting the simulation approach using crop models. In the GYGA project, the possibility for achieving the world's food production capacity is expressed in the form of yield gaps (Yg), i.e. the difference between potential (Yp) and actual yields (Ya). The visualization has been aiding the next step to follow, i.e. closing and narrowing Yg (Hochman et al., 2016). What makes GYGA different from other yield gap studies (Sentelhas et al., 2015) is that its original Climate Zonation Scheme (CZS) was developed in the project (Van Wart et al., 2013) allowing it to be used on a global basis. The climate here is used in a conservative context in accordance with the definition 'the weather conditions prevailing in an area in general or over a long period' (Lexico Dictionaries, 2019).
The Köppen climate classification system (Köppen, 1936;Peel et al., 2007) developed by the Russian-German botanist-climatologist Wladimir Köppen (Arnfield, 2021) is a means to categorize climate zones over the world. Something special about the Köppen climate classification is that it depicted world climate defined by temperature and precipitation in relation to terrestrial vegetative biomes. In the Köppen's system, most of the areas in Japan fall into four classes, i.e. Dfa, Dfb and Dwb (humid continental climate) and Cfa (humid subtropical climate). These four climate categories are, however, too coarse to accommodate the climatic heterogeneity that is relevant to agriculture across a mountainous island landscape as found in Japan. The GYGA-CZS is composed of three factors, i.e. Growing Degree Days (GDD: an index to address temperature), Aridity Index (AI: an index to address precipitation and evapotranspiration) and Temperature Seasonality (TS: an index to address annual temperature range) (Van Wart et al., 2013; GYGA, 2021). The number of combinations of the three factors theoretically amounts to 300, approximately 20 of which have been observed in Japan. This would suggest that the CZS may be more suitable than the Köppen climate classification in explaining yield differences of crops between regions. A growing number of papers have been published using CZS (Edreira et al., 2018;Deng et al., 2019), although its validity especially in terms of TS, the index that is considered to differentiate temperate and tropical climates (Edreira et al., 2018), has not been tested explicitly. Hochman et al. (2016) raised the point that information contained within smaller climate zones can be missed out under the current GYGA protocols where a 100 km buffer zone around the selected Reference Weather Station (RWS) defines the smallest scale. They compared GYGA with a data-rich approach in a study performed on wheat in Australia. Actual annual wheat grain yields were obtained from the national agricultural data collated by the Australian Bureau of Statistics at the level of statistical division annually, and at the finer scale of statistical local area every 5 years when a census was performed. As Hochman et al. (2016) have shown, there is the potential to lose information under the current GYGA protocols; this would be a concern for a small and long-shaped country like Japan. The width of the main island is very often below 100 km, the size of buffer zones defined in the GYGA protocol. In addition, a chain of mountains runs through the main island as if it were its spine creating climate gradation caused by diverse altitudes. This is where an alternative approach needs to be explored so that Japan and other countries in similar positions could also participate in the GYGA project and seek a way to contribute to the problem of future global food security. Apart from the difficulties associated with its geography, rice production in Japan might provide a perfect test case to explore the potential of the CZS as a tool to explain yield differences for three reasons: (1) Rice cultivars of similar backgrounds have been produced from north to south (Yamamoto et al., 2010). (2) Statistical records of rice yield over the last 30 years on a municipality scale are easily accessible. (3) Most rice in Japan is irrigated, removing the necessity to include AI, an index related to available water to crops. Removal of AI makes the CZS a simplified scheme composed of only two factors, i.e. GDD and TS.
The prime objective of the present study was to find out using actual rice yield (Ya) if the GDD and TS scores defined in the GYGA-CZS could be utilized to identify rice production areas of favourable characteristics, especially those that are highyielding and stable-yielding. To do this, we developed an alternative approach to deal with CZS in order to overcome the relatively small geographical scale and mountainous topography of Japan as well as to utilize statistical information stored, mostly, on a municipality basis. This allowed the evaluation of the GYGA-CZS using estimated Yg, i.e. the gaps between simulated Yp by a crop growth model and Ya.

Identification of climate zones on a municipality basis
The original map of GYGA-CZS (GYGA, 2021) is composed of approximately 10 km × 10 km grid-cells of various Climate Zones (CZs) expressed as four or five-digit-integers. The integers are defined by the combination of three factors, i.e. GDD score (1000, 2000, 3000 … 10 000, where smaller numbers coincide with colder climate), AI score (000, 1000, 2000 … 900, with smaller numbers coinciding with drier climate) and TS score, in other words, annual temperature range (1, 2 or 3: smaller numbers coincide with smaller annual temperature range) (Fig. 1). So, for example, a CZ score of 4903 has a GDD score of 4000, an AI score of 900 and a TS score of 3, indicating a cool, wet climate with large variation in temperature over the course of the year. The actual ranges for each factor are listed in Supplementary information. However, it is on a municipality basis, not on a gridcell basis, where statistical information such as crop yields is accumulated and policies are implemented. In the present study, the CZs in the framework of GYGA-CZS were identified on a municipality scale to enable an analysis of rice yields recorded in governmental statistics in the framework of GYGA-CZS. CZs of 1718 municipalities that were present in 2016 were identified by visual observation using QGIS Desktop (v2.18.13). Municipalities covered by only a single CZ and those covered by more than two CZs where one of them covers the largest area were subjected to the analysis (Fig. 2). The same identification was conducted on the municipalities that were present in 1993, as many municipalities were merged around 2004 reflecting partly the falling and ageing population. The number of municipalities was 3232 and 1727 on 31 March in 1999 and on 31 March in 2010, respectively (e-Stat, 2021). In total, 23 CZs (comprising eight GDDs and two TSs) were identified in Japan in the original GYGA map ( Fig. 3  (a)). The CZs of 1803, 3802 and 4703 were omitted from the maps as well as the analysis due to their small shares (<0.1%). Similarly, the CZs of 1902, 7901, 7902, 8901 and 8902 were omitted as well due to the lack of significant rice production (MAFF, 2021). As a result, 15 CZs (comprising six GDDs and two TSs) were considered in the map as well as the analysis between 1993 and 2002 ( Fig. 3(b)) and between 2005 and 2016 ( Fig. 3 (c)). In Fig. 3(d ), the regions of Japan referred to in the present study are presented to facilitate reading.

Data collection
Rice yields recorded in each municipality between 1993 and 2016 were collected from the Statistical Survey on Crops of the Ministry were omitted from any analysis in the present study, as no rice yield was recorded in many municipalities possibly due to the merger of municipalities that occurred in that year. Weather data recorded at selected weather stations and locations by the Automated Meteorological Data Acquisition System (AMeDAS) were obtained from MeteoCrop DB (2021), the website run by the National Institute for Agro-Environmental Sciences (NIAES), currently known as the Institute for Agro-Environmental Sciences, National Agriculture and Food Research Organization (NARO) and the website of Japan Meteorological Agency (JMA, 2021). The target crop being the irrigated rice, collection of soil data was not conducted due to the GYGA protocols (GYGA, 2021).

Estimation of yield gaps
Yg was calculated by subtracting actual yield (Ya) from potential yield (Yp). Unlike Ya that was taken from the real world, Yp needed to be simulated by means of a crop growth model such as ORYZA v3 in the case of irrigated rice. Possible issues of using ORYZA to estimate rice yields in a temperate climate have been pointed out by Espe et al. (2016b). Therefore, we tested crop models developed in Japan, including SIMRIW (Horie, 1987;Horie et al., 1992) and CYGMA (Iizumi et al., 2017(Iizumi et al., , 2018 for their suitability to this study. SIMRIW is specifically designed for irrigated rice in Japan, while CYGMA is a generic crop model currently applicable to six major crops, including rice. The Yp values simulated by SIMRIW and those simulated by CYGMA were comparable in most of the areas employed for the test runs. In the present study, CYGMA was employed to calculate Yp for all locations of AMeDAS. CYGMA originated from a model that was designed to calculate agro-climatic indices to simulate potential crop growth under optimal conditions (Iizumi and Ramankutty, 2016). The model operates at a daily time step to resolve major eco-physiological processes to simulate growth and yields under rainfed and irrigated conditions separately at a given location. In the model, crop development is computed as a fraction of the accumulated GDD (with the base temperature of 10°C for rice) relative to the crop thermal requirements. Leaf growth and senescence are calculated according to the fraction of the growing season using the prescribed shape of the leaf area index curve. Maximum leaf area index usually reached around the time of heading was set from 4 to 7 m 2 /m 2 (Yoshida, 1972). The details of the CYGMA model can be found elsewhere (Iizumi et al., 2017). Phenology data, i.e. dates of transplanting, heading and maturity of major cultivars, i.e. the cultivars planted in the largest area in each prefecture, were taken from the governmental crop progress reports' characteristics tables of recommended cultivars of irrigated rice (MAFF, 1993(MAFF, , 1995(MAFF, , 1997(MAFF, , 1999(MAFF, , 2002(MAFF, , 2003(MAFF, , 2005(MAFF, , 2008(MAFF, , 2009(MAFF, , 2011(MAFF, , 1993. All cultivars relevant to the present study were japonica type. The crop thermal requirements were computed year by year and used for the simulation. Daily mean 2 m air temperature and global solar radiation data obtained from AMeDAS were used as the weather input.
Yg was estimated by subtracting Ya from Yp. It should be noted here that Ya and Yp come from different sources: Ya was taken from statistical records on a municipality scale, while Yp was simulated using weather data taken from AMeDAS locations. As a way of combining the two, all municipalities that have records of rice yields in the statistical information database were manually related to one of the nearest locations of AMeDAS considering primarily the CZs and secondarily topography of both the target municipality and candidate locations of AMeDAS. The simplest combination occurred when a given municipality belonging to a certain CZ could be related to an AMeDAS location situated in the same CZ as the target municipality. It happened, however, that the CZ of some municipalities differed from that of any relatable location of AMeDAS. Yg values were estimated in such cases in order to draw a map of Yg; however, they were not subjected to the statistical analysis described in the next section. As the number of locations of AMeDAS was
smaller than that of municipalities, all combinations of Yp and Ya involving the same location of AMeDAS were averaged to generate a single value (Fig. 4). Some Yg values were observed to be negative and considered to be partly due to the mountainous topography of Japan. The percentage of negative values was lower than 1% of the total number of Yg calculations in any year and the negative values were treated as missing values. Yg was not estimated in 1993 and 2003 due to severe cold damage that occurred especially in the Hokkaido and the Tohoku regions, while the lack of yield data, probably as a consequence of the merger of municipalities, did not allow the estimation of Yg in 2004. comparisons by Scheffé were performed setting GDD score and TS score of CZs as fixed factors using SPSS Ver. 24 (IBM). The numbers of data sets subjected to ANOVA are presented in Table 1. Statistical yields were analysed, while CV was calculated by summing the observations over 21 years to increase the number of observations as much as possible by selecting municipalities that were existent both before and after the merger. Estimated Yg values were subjected to ANOVA and multiple comparisons by Scheffé following the method presented in Fig. 4. Yg data sets for the 9 years of data sets between 1994 and 2002 and the 12 years between 2005 and 2016 were each bulked and then subjected to ANOVA with GDD and TS scores as fixed factors and the year as a random factor, respectively. QGIS Desktop 2. 18. 13 was used for the presentation of maps. Shapefiles of blank maps of Japan in 1993 and 2016 were sourced from date-specific administrative boundary data (Tokyo Map Research Inc., Tokyo) and were used to draw administrative boundary maps of Japan in Figs 3(a)-(d ), 7(a), (b) and 9.

Identification of climate zones on a municipality basis
The original map of GYGA CZs expressed on a grid-cell basis is shown in Fig. 3(a)

Rice yields
Following ANOVA and the multiple comparison analysis, significant interactions between GDD score and TS score for rice yield (Ya) were observed for 8 years out of 9 between 1994 and 2002 (Table 2 and Fig. 5). Overall rice yields (Ya) excelled in the CZs characterized by the TS score 3 (P < 0.001), which  characterizes a greater annual temperature range, compared to those CZs characterized by the TS score 2. The combination of TS score 3 and GDD score 3 tended to give higher rice yields (P < 0.001). Similar trends were observed for the data sets between 2005 and 2016, although interactions between GDD score and TS score were observed less frequently (Table 3 and Fig. 5). Interactions between GDD score and TS score were observed for CV (P = 0.009). The CZs characterized by the GDD score 3 had a lower CV when they were characterized by the TS score 3 than when they were by the TS score 2 (Fig. 6). Overall, the CZs characterized by the GDD score 3, 4 and 5 had a lower CV compared to those characterized by the GDD score 1, 2 or 6 (P < 0.001) (Fig. 6). The areas covered by the CZs characterized by the TS score 3 appeared to have coincided with mountainous areas in the eastern part of Japan. The Kanto plain and the western part of Japan were mostly covered by the CZs characterized by the TS score 2 (Figs 3(a)-(c)).    score of 2 irrespective of GDD score. The difference, however, in Yg between CZs with a TS score of 2 or 3 was reduced when the CZ had a GDD score of 4 implying that Yp was exploited better in these CZs. The lowest Yg were found from the CZs with a GDD score of 3 and a TS score of 3.

Identification of climate zones on a municipality basis
Identification of CZs in each municipality was, without doubt, a time-consuming task. Doing so, however, made it possible to utilize a larger number of statistical data sets prepared by MAFF (2021). Considering that research and the release of cultivars have been performed by each prefecture, it seems meaningful to identify the GYGA CZs on a municipality basis. Identification of the GYGA CZs on a municipality basis allowed representative rice-growing environments in Japan to be identified without missing any major producing area. The method presented in the present study would be applicable to countries and areas where taking 100 km buffer zones of homogenous CZ around selected RWS as in the current GYGA protocols is found to be difficult for mostly two reasons, i.e. small size and heterogeneous topography. The concern with the current GYGA protocols was also addressed by Hochman et al. (2016) who studied wheat in Australia. They found the major difference between the original GYGA protocol and a data-rich approach was in the coverage of the target crop area rather than in explanatory power, and they commented that the GYGA protocols were successful in an unbiased sampling of the target crop areas.

Rice yields and yield gaps under GYGA climate zonation scheme
High rice yields were obtained in the CZs characterized by the TS score 3 (Fig. 5). Akita prefecture (the Tohoku region), Yamagata prefecture (the Tohoku region) and Nagano prefecture (the Koshin region) known for high rice yields had large areas of the CZs of 3903 and 4903. The CZ of 5903 happened to coincide with the Niigata city, the rice production area of the largest share among all municipalities producing rice in the country. The CZ of 5903 is limited to Niigata prefecture (the Hokuriku region) and small areas in Ishikawa (the Hokuriku region), Fukui (the Hokuriku region), Gifu (the Tokai region) and Okayama (the Chugoku region) prefectures. TS scores in the GYGA-CZS were defined as the standard deviation of the mean temperature of 12 months reflecting the annual temperature range (GYGA, 2021). It was speculated, however, that rice yield was possibly influenced by the diurnal temperature range considering that annual temperature range is often related to diurnal temperature range (Scheitlin, 2013), both of which are associated with the degree of maritime influence, in other words, the distance to the sea (Scheitlin, 2013). Fields located in and behind mountains, however, are less likely to be influenced by the sea. Experienced Japanese rice growers know that rice yields well in a climate of cool night-time temperatures, which is an observation made overseas as well (Espe et al., 2016a). In a report published in the 1950s, rice grown in a basin in Fukushima prefecture (the Tohoku region) was found to give a greater number of panicles, a greater percentage of effective tillers and greater weights of both grains and straws per hill compared to the crop grown along the coastal area (Yamamoto, 1953). From more recent experiments conducted in growth cabinets, a rise in night temperature was reported to have reduced grain weight of rice, while day temperature did not have such an effect (Morita et al., 2002). The CV of Ya tended to be high in the CZ of 2903, the coldest area among rice-cultivating areas in Japan indicating that cold weather is one of the important factors influencing the stability of rice yields, a situation well recognized by rice researchers around the globe (Da Cruz et al., 2013;Shimono, 2018). It should be noted that rice cultivars achieving both cold tolerance and good taste became the target of breeding, especially in the Hokkaido and the Tohoku regions during the last quarter of the 20th century (Motoki, 1999), and with some success (Sasaki et al., 1990;Yoshimura et al., 2002;Mikami et al., 2007). The high CVs observed with CZs characterized by the GDD score 6 indicate that the current japonica rice cultivars are not necessarily suited to these areas of warm environment, a situation wellrecognized among rice-research communities (Morita, 2008;Kim et al., 2011;Samol et al., 2015;Yang et al., 2017). Partly to address the issue, a rice cultivar named 'Nikomaru' that shows resistance to high temperature during the ripening period (Maeda and Watanabe, 2013) has been bred and released by NARO (2021).
Both GDD score and TS score were found to be useful in giving an idea before planting of the level of rice yields likely in a given CZ. It is probable, however, that CZs will need to be reappraised regularly to account for changes caused by, for example, global warming, a concern previously raised for and worked on in relation to the Köppen climate classification (Peel et al., 2007;Beck et al., 2018). ANOVA and multiple comparisons performed on Yg suggested that the potential yielding ability of rice is well exploited in the CZs characterized by the GDD score 3 and 4 for japonica rice cultivars planted in Japan over the last two decades. Potential yielding ability appears to be utilized better in the CZs characterized by TS score 3 than those characterized by TS score 2. The interaction observed between GDD score and TS score with respect to Yg implied that CZs characterized by TS score 3, related to 'continentality', were well suited to exploit rice's yield potential as was the CZ characterized by the GDD score 4 and TS score 2. It might be, therefore, worthwhile taking the effects of TS score, in other words, the effects of annual temperature range and/or diurnal temperature range, into consideration when estimating Yp.

Possibility for future rice production
In general, Yg tends to be associated with insufficient inputs such as fertilizers and agrochemicals, although the situation with Japanese rice could be different. Nishio (2002) pointed out that between 1980 and 1999, the reduction in the rate of fertilizer applied to rice was accompanied by the fall in the government purchase price and that rice yields during the same period, however, stayed the same. His speculation as to the cause of this phenomena included the increased use of controlled-release fertilizer and expansion of paddy areas planted with Koshihikari, the symbol cultivar of good taste. It might be surprising for some readers that Koshihikari, the most commonly grown rice cultivar in Japan today has been cultivated for more than 60 years: it was bred in 1956 in Fukui prefectural agricultural experimental station (the Hokuriku region). Most of the cultivars grown in large areas often share similar genetic backgrounds with Koshihikari (Saito et al., 1989;Sasaki et al., 1990;Yagi et al., 1990). Japanese consumers appreciate the taste of rice and therefore tend to choose brand cultivars and the brand place where the particular rice was grown. In this type of environment, growers are not necessarily likely to focus on obtaining very high yields and minimising Yg over brand acceptability. Identifying prefectures where either Koshihikari or an apparently related cultivar has ever shared the largest area of rice cultivation between 1993 and 2016 revealed that the only potential exception appeared to be Gifu prefecture (the Tokai region) confirming the broad and deep involvement of Koshihikari in Japan's rice cultivation (Fig. 9). Even Gifu is not the exception, however, as Hatsushimo, the major rice cultivar grown there also shares the same ancestors such as Asahi and Ginbozu with Koshihikari (Institute of Crop Science, 2021). Yoshida (2001) argued that the genetic backgrounds of rice cultivars in Japan had been rather narrow even before the advent of Koshihikari. Similar issues with a genetic base of irrigated rice can be found in other countries as well (De Oliveira Rabelo et al., 2015).
From the viewpoint of both genetic diversity and possible expansion of potential rice production capacity, the comparison of rice yields in terms of the GYGA-CZS with other countries such as China and the USA would potentially be useful for policymakers in Japan. As rice yields as high as 6-10 t/ha have been
obtained from areas of comparable CZs in China and the USA (GYGA, 2021), it can be deduced that higher rice yields could be achieved in Japan by adopting updated cultivars with greater yielding ability. This has already been confirmed by rice cultivars newly bred for usages other than human consumption such as livestock feed (Hirabayashi et al., 2010;Fukushima et al., 2018). With a change in social environment, e.g. market acceptability of new varieties, one might be able to harvest a greater quantity of rice in Japan in the future, a piece of information useful for policymakers concerned with food supply. Areas characterized by the TS score 3, a greater annual temperature range, tended to give higher rice yields. These areas are often associated with higher altitudes compared to the areas characterized by the TS score 2. Mountainous areas in Japan have often been accused of lower efficiency in terms of mechanized farming, namely, scattered small patches of fields and unpopular slopes that require skilled operators, compared to the fields located in plains. In addition, the continental nature of the area characterized by the TS score 3 in the present study has sometimes been perceived as a negative element for rice production, especially in cold seasons, though this could at least partly be alleviated by cultivation techniques such as deep-water management (Watanabe et al., 2006). Admitting that we avoided using data sets from years of severe cold damage (1993 and 2003) in the present study, the analysis of Ya under the GYGA-CZS casts doubt on the view that these fields are inherently less suitable for rice production. Focusing rice production in fields on the plains that often happen to coincide with the CZs characterized by the TS score 2 may not necessarily be an optimal policy decision.
With the predicted expansion in global population and the increased need for food, decision-making that pursues 'apparent efficiency' when there might be 'hidden efficiency' that we are unaware of at present would not be advantageous. Quantifying the benefit of TS score 3 would be possible by multiplying the average rice yield of TS score 3 to its area, which will be addressed in a future study. The GYGA-CZS, when performed on a municipality scale and if properly updated in accordance with ongoing climate change, could serve as a strong tool to evaluate crop production in the climates in Japan and possibly in other parts of the world, especially the Korean Peninsula and China where japonica-type rice is produced in CZs with the TS score of 3.
Supplementary material. The supplementary material for this article can be found at https://doi.org/10.1017/S0021859621000186