Model‐based evaluation of climate change impacts on rice grain quality in the main European rice district

Abstract

Crop simulation models are used to forecast the impacts of climate change on yield levels and to identify adaptation strategies. Nevertheless, crop quality has been almost neglected in available studies, despite its relevance on the economic and nutritional value of agricultural products. We present here a modeling study to evaluate the future trends of rice quality in the main European rice district, placed in Northern Italy. A rice growth model was coupled with a library of models of rice milling and cooking suitability, using current farmer management and baseline/future climatic scenarios as input for the simulations. Four general circulation models (NOResm, MIROC‐ESM, HadGEM2‐ES, and GISS‐ES) and two CO₂ representative concentration pathways (2.6, 8.5) were used to generate 20‐year future climatic data centered on 2030 and 2070. Spatially distributed simulations were run at 2 × 2 km spatial resolution considering a Japonica (Loto) and a Tropical Japonica (Gladio) rice cultivar. The results depicted an overall decline in rice quality, especially for Loto (−5% of milling suitability considering GISS‐ES‐2.6 in all time frames; −8% in 2030 and −20% in 2070 under HadGEM2‐ES‐8.5). While the revenues of millers are expected to decrease of about 50 € t⁻¹ in 2030 and 100 € t⁻¹ in 2070 for Loto, minor changes are projected on Gladio milling and cooking suitability, except in the worst scenario in 2070 (−10 to −5% of cooking suitability, corresponding to −30 to −72 € t⁻¹). Despite the need of reducing models uncertainty, this study provides variety‐specific indications on rice grain quality ready‐to‐use for crop specialists, farmers, and millers and that could raise the interest of different stakeholders of the agri‐food sector, including food scientists, geneticists, and policy makers.

INTRODUCTION

Crop quality is a key determinant of the economic and nutritional value of agricultural products, as it influences their purchase by consumers and their acceptability to buyers (Koutroubas et al., 2004). The definition of food security provided by the Food and Agriculture Organization of the United Nations (FAO) lists the quality of agricultural products as a principal pillar, intended as the production of nutritious food to allow people to meet dietary needs and food preferences for an active and healthy life. Savary et al. (2017) combined three different definitions of food security to provide a new framework consisting of six components of food security, two of which are especially relevant in European economically wealthy countries, where rice production must meet superior quality standards in line with national and European regulations while satisfying the consumers' preferences. These two components are as follows: (i) the stability of food availability, intended as the availability of adequate food irrespective of sudden shocks (e.g., climatic) at all time and (ii) the Utility-Safety-Quality-Nutritive value, which is based on the suitable intrinsic value of food as nourishment. The achievement of high quality standards is of paramount importance also for farmers, in order to gain competitive advantage in the internal market as well as toward low price/quality exporting countries (Griglione et al., 2015). Recent agricultural policies in Europe aimed at enhancing the qualitative aspects of processes and products and to safeguard the environment through a more responsible use of resources, besides merely focusing on crop productivity (European Commission, 2018).

Over the last 50 years, despite the yields of major crops have steadily increased due to genetic improvement and technological progress (Edgerton, 2009; Martre et al., 2011), the quality of productions has leveled off or even decreased for cereals (Oury et al., 2003), grain legumes (Graham & Vance, 2003), oilseeds (Triboi & Triboi-Blondel, 2002), fruits for processing (Grandillo et al., 1999), or fresh fruits (Causse et al., 2003). The few available studies on the impacts of climate change on crop quality are aligned to this trend, depicting an overall decline in many cereal crops (DaMatta et al., 2010), especially for rice (Nagarajan et al., 2010; Sreenivasulu et al., 2015). Among the causes of the worsening of rice quality, the increased frequency of heat stress during ripening is responsible for the alteration of sink strength and storage composition in the grains, with a direct impact on the nutritional, rheological, and milling attributes (Sreenivasulu et al., 2015). Many other physiological processes contribute to decrease the quality of rice productions, such as the shortening of the grain filling period (Ishimaru et al., 2009), the reduced nitrogen translocation from the leaves to the panicles (Cheng et al., 2010), the predominance of the nighttime respiration on photosynthesis (Sreenivasulu et al., 2015), the non-uniform moisture desorption from the grain surface before harvest (Lan et al., 1999) and the reduced biosynthesis of starch and its concurrent degradation via α-amylase (Yamakawa et al., 2007). Given the paramount importance of rice as staple food, with increasing consumption even in non-rice-eating countries (e.g., +3% year⁻¹ in northern Europe and +5% year⁻¹ in the United States; Suwannaporn & Linnemann, 2008), research efforts are needed to anticipate the future dynamics of qualitative aspects of productions, and to put in place effective adaptation strategies.

Process-based simulation models represent effective research tools to address these questions, due to their capability to reproduce nonlinear responses to boundary agrometeorological conditions (Porter & Semenov, 2005). Despite a variety of models has been developed to simulate pre-harvest grain quality (Cappelli et al., 2014; Martre et al., 2011), very few studies focused on the future trends, for example, of grain nitrogen concentration in wheat (Asseng et al., 2013, 2019; Erda et al., 2005; Porter & Semenov, 2005) and of chalkiness in rice (Masutomi et al., 2019; Okada et al., 2011; Takimoto et al., 2019). These studies agree in identifying the rationalization of nitrogen fertilizations, the optimization of sowing dates, and the use of heat-tolerant cultivars as appropriate adaptation strategies to counterbalance the quality decay in cereals. The projected decay in rice quality due to climate change strengthens the need of modeling studies to forecast its trend in the near future, considering the connections with the stability and utility-safety-quality-nutritive components of food security. Indeed, the expected increase in thermal regimes and in the frequency of weather extremes may cause higher incidence of damaged grains and reduced head rice yield even in rice top producing areas (Cappelli, Pagani, et al., 2018), causing a decreased stability of food production and supply and of the nutritive, eating and cooking value of rice grains. Decrease in cooking suitability could exacerbate grain quality decay and related decreases in market value.

We present here a simulation study aimed at evaluating the impacts of climate change on rice milling and cooking suitability performed on the main European rice district, placed in Italy and accounting for 50% of total harvested area in EU-27. National breeding programs aimed at improving rice milling and nutritional aspects as well as at promoting low-impact agronomic practices (i.e., organic rice varieties, Griglione et al., 2015) were recently launched (Russo & Callegarin, 1997) to support rice stakeholders in identifying priorities for future investments. The objective of this study is to support these programs by assessing the future trends of rice quality in the short (2030) and long (2070) term, via the spatially distributed application of process-based quality models.

MATERIALS AND METHODS Study workflow

The workflow of this study articulates in six steps. Firstly, we developed a simulation environment by coupling a rice growth model with ten simulation models of rice quality traits under non-limiting conditions for water and nutrients. The quality traits considered were amylose and protein content, viscosity properties, head rice yield, milky white, broken, and pecky grains. Secondly, we derived the model parameterizations from Cappelli, Pagani, et al. (2018) and Cappelli, Confalonieri, et al. (2018), where a multi-site, multi-year (2006–2014) evaluation of rice grain quality models was performed on Loto (japonica) and Gladio (tropical japonica) cultivars across 16 sites in the main European rice area (Lombardo Piemontese district, Northern Italy, study area). Thirdly, we performed spatially distributed simulations (2×2 km spatial resolution throughout the study area) of grain quality variables of Loto and Gladio cultivars in the baseline scenario (2003–2014), considering current farmer management and potential conditions for water and nutrients. The variability in quality traits, their internal relationships and the differences between cultivars were analyzed. Fourthly, we assessed the future trends in rice quality variables (2 × 2 km simulations, no adaptation), using 20-year data series centered in 2030 and 2070, generated by combining four general circulation models (NOResm, MIROC-ESM, HadGEM2-ES, and GISS-ES) and two CO₂ representative concentration pathways (2.6, 8.5). The impact of climate change was evaluated as percentage difference compared to the baseline for each quality variable. Internal relationships between quality variables were also inspected. Fifthly, we aggregated rice quality variables according to their effect on cooking (amylose, protein, and viscosity properties) or milling (i.e., head rice yield, milky white, pecky, and broken kernels) suitability of rice grains by considering preferences of the Italian consumers. Expected future trends in rice quality were explored, analyzing spatial and temporal patterns of milling and cooking quality. The combined effect of air temperature and of the length of crop cycle on the variability of rice quality was also investigated. Sixthly, we performed economic projections on the future profitability of rice using percentage changes in milling and cooking quality compared to the baseline and the price that milled rice is currently exchanged in the local market.

Rice cultivation in the study area

The Italian Lombardo-Piemontese district, representing the leading rice-growing area in Europe with 1,512,000 t produced on 230,000 ha in 2018 (Ente Nazionale Risi, ENR; www.enterisi.it) and an overall turnover of about one billion euros (National Institute of Agricultural Economics Research, CREA, 2017) was chosen as the study area (Figure 1).

FIGURE 1. Input information of the spatially distributed simulations performed in the Lombardo-Piemontese district. (a) Rice distribution map, (b) time windows of sowing dates (DOY = day of year), (c) average air temperature and (d) mean cumulative precipitation in the period 2003–2014. Mean temperature and precipitation changes expected in 2030 and 2070 are presented as the average ± one standard deviation of the projections considering all combinations of emission scenarios × general circulation models

The climate in the area is temperate humid, with large variability of pedo-climatic conditions (Fumagalli et al., 2011). Average annual air temperature is 13.5℃ (0–4℃ in winter months, with peaks over 30℃ in summer), and annual precipitation is around 750 mm. Main limiting factors during the rice-growing season (May–October) are low temperatures at sowing, cold-induced sterility due to air irruptions around heading, weeds infestation, and blast disease epidemics. The main irrigation strategies are continuous flooding of paddy fields, or dry sowing with delayed flooding at the third-fourth leaf stage. Nitrogen is mainly distributed in the form of urea in two or three events, one in pre-sowing, and the others at beginning of tillering and/or at panicle initiation. According to EU standards, Italian rice cultivars are categorized in four groups based on the grain shape and length/width ratio: “long A,” “long B,” “medium,” and “round” (European Parliament, 2013). Around 70% of the rice cultivars grown in the area belong to the Japonica type (average yield of 6.5 t ha⁻¹) and 30% to Tropical Japonica (average yield of 7 t ha⁻¹), with cycle length ranging from short (120–140 days) to medium-long (140–170 days). Among Japonica cultivars, long rice (long A) with chalky grains (i.e., with short spots of pearl in the ventral—white belly, dorsal—white back—or central—white core—part of the grains) and intermediate (i.e., 20%–25%) amylose content is predominant and mainly used for the parboiling process and risotto dishes. Round and medium rice grains are instead characterized by low amylose content (i.e., <20%) and are utilized for soups, timbales, dessert, and puffed rice preparation. Tropical Japonica rice (long B) includes both high and intermediate amylose cultivars that are used for salads, side dishes, or subjected to parboiling process. The export of rice production is limited to chalky grains to both EU and non-EU countries (above 60% of total milled rice), accounting for approximately € 500 million in 2016 (CREA, 2017). The sustainability and the innovation in rice sector are supervised by the ENR (National Rice Authority) and by the Council for Agricultural Research and Economics (CREA), which cooperate on research programs to achieve yield stability and superior grain quality by integrating agronomic management, genetic improvement, and cultivar selection.

The simulation environment Dynamic simulation of the rice cropping system

The WARM model (Confalonieri et al., 2009) was used in this study to simulate rice phenology, micrometeorology, and qualitative aspects (Cappelli et al., 2014). Phenology development was simulated based on hourly thermal time accumulation, considering a development stage code (DVS) ranging from 0 to 4, where 0 = sowing, 1 = emergence, 2 = flowering, 3 = maturity, and 4 = harvest. The nonlinear temperature response curve (Yan & Hunt, 1999) is driven by the minimum, optimum, and maximum crop cardinal temperatures. The floodwater effect was estimated via a micrometeorological model simulating hourly temperatures of water body and of air layers above the air-water interface according to an energy balance approach (Confalonieri et al., 2005); temperature estimated at the meristematic apex layer was used as input for the simulation of crop development. The simulation of rice grain quality during ripening was performed by models driven by weather data, that is, temperature, radiation, rainfall, wind speed, and relative air humidity, needing as input dynamic information on rice phenology (i.e., DVS). These models were presented in Cappelli, Pagani, et al. (2018) and Cappelli, Confalonieri, et al. (2018) and consist of broken linear, curvilinear, or logistic response functions driven by biological parameters which reduce the daily rate of change or the final value of a given quality variable. Considered physiological processes were the synthesis of chemical compounds (i.e., starch, amylose, and protein), the length of the sensitive period to weather conditions (e.g., number of days/growing degree days from the onset of heading/flowering) and the cultivar-specific susceptibility to heat stress damage. The quality variables considered in this study were the amylose (AC) and protein concentrations (PC), the indices describing the starch viscosity during cooking (i.e., the gelatinization temperature—GT, the breakdown—BDV, the setback—SV, and the peak viscosity—PV), the percentage incidence of severely chalky (Milky white—MWG, almost entirely chalky, opaque, and floury grains), broken (BK), pecky grains (PG), and the head rice yield (HRY). MWG, PC, and SV were dynamically simulated using a daily time step, whereas the other quality variables were estimated at physiological maturity. Simulations were performed at 2 × 2 km resolution on the whole Lombardo-Piemontese rice district under current weather (baseline) and future climate projections, assuming non-limiting conditions for water, nutrients, pests, and weeds.

Models selection and parameterization

The Italian cultivars Gladio and Loto were selected as representative of Tropical Japonica and Japonica type genotypes due to their high representativeness in the study area and the availability of calibrated parameters sets (Cappelli, Confalonieri, et al., 2018). Gladio is a short cycle long B cultivar released in 1998, presenting long slender crystalline grains recommended for rice salads, side dishes cooking, or parboiling, with an average yield of 6.3 t ha⁻¹. Loto is an early long A cultivar released in 1988, with medium slender grains suitable for parboiling and risotto preparation, with a slightly lower yield than Gladio (5.8 t ha⁻¹). The features of the two cultivars in terms of cooking and milling quality are summarized in Table 1.

TABLE 1 The quality characteristics of the Italian rice cultivars Loto and Gladio

QV	Description/properties	Unit	Average range values
			Loto	Gladio
Cooking suitability
AC	Determinant of grain texture and stickiness of boiled rice	% grain−1	14–18.7	23.3–27.2
PC	Determinant of nutritional value, taste, and cooking properties	% grain−1	4.9–7.5	5.0–7.6
GT	Mean thermal transition of the starch in presence of water and heat. Determinant of cooking duration and the texture	°C	65.5–71.3	73.4–79.4
BDV	Starch granules disruption. Measure of paste stability	BU	479–772.5	480–662.5
PV	Maximum viscosity during the heating phase. Measure of granules swelling during gelatinization	BU	1021–1257.0	934–1137.5
SV	Starch retrogradation tendency during cooling phase. Determinant of grain texture	BU	591–685.5	665–846
Milling suitability
BK	Amount of fissured kernels	%	5.3–23.6	3.3–8.9
HRY	Relative weight presence of entire kernels after milling	%	47.6–70.8	62.6–68.6
MWG	Amount of almost entirely chalky, opaque, and floury grains	%	0.8–6.1	0.2–1.6
PG	Damaged grains with circular, brown to black spot	%	0.8–7.6	0.2–3

Note:

Quality variables (QV) are grouped according to their effect on cooking or milling suitability of rice grains. BU = Brabender units: a procedure defined unit to measure the viscosity of a slurry (Karim et al., 2000).

Abbreviations: AC, amylose concentration; BDV, breakdown viscosity; BK, broken kernels; GT, gelatinization temperature; HRY, head rice yield; MWK, milky white grains; PC, protein concentration; PG, pecky grains; PV, peak viscosity; SC, setback viscosity.

The model parameters driving rice phenological development (i.e., the thermal time thresholds between consecutive phases and cardinal temperatures for development) and grain quality were taken from Cappelli, Pagani, et al. (2018) and Cappelli, Confalonieri, et al. (2018), where a multi-site, multi-year (2006–2014) evaluation of rice quality models was performed on Loto and Gladio cultivars in the same study area. For each rice quality variable considered, the best model considering the accuracy of fit (modeling efficiency, EF, Nash & Sutcliffe, 1970; minimum = −∞, optimum, and maximum = 1, unitless; relative root mean square error, RRMSE, Jørgensen et al., 1986; minimum and optimum = 0%; maximum = +∞) and the complexity (Akaike Information Criterion index, AIC, Akaike, 1974; the lower, the better) was selected. RRMSE measures the relative mean difference between simulated and observed data, EF indicates the rate of variation explained by the model, AIC considers the trade-off between the accuracy of fit (explanatory power) and the complexity (number of variables and parameters). The ten models are presented in Table 2, together with their inputs and calibrated values of parameters for Loto and Gladio cultivars. The adopted parameterization reflected the high susceptibility of Loto cultivar to the weather variability (i.e., heat stress, high windiness, and sudden air humidity changes at suboptimal temperatures) in the post-flowering/heading period, and reproduced the lower size, the shorter ripening, and the higher thermal requirements of Gladio development. The algorithmic description of the ten models is provided in Text S1, whereas the calibrated values of the WARM parameters targeting the simulation of crop phenology and the synthetic description of models' performances in reproducing phenological development and grain quality are reported in Table S1.

TABLE 2 Inputs and calibrated parameters of the quality models used in this study

QV	Unit	Inputs	Biological parameters	Unit	Value
					Loto	Gladio
Cooking quality
AC	% grain−1	Tavg, Rad, WS, Daysff, DVS	Tavg above which starch begins to be sensitive to sub-opt radiative and wind conditions	°C	22.9	23.0
			Daysff at which sensitive period to Temp. begins	d	10.0	10.0
PC	% grain−1	Tavg, Tnight, sunriseH, sunsetH, Daysff., DVS	Sensitivity to sub-opt T for protein synthesis (0 = max; 1 = none)	—	0.54	0.75
			Duration of sensitive period to T	d	16	16
			Average grain protein content		6.5	6.8
			Tmax for protein synthesis	°C	35	36
			Tmin for protein synthesis	°C	12	12.5
			Higher Tnight for opt. protein synthesis	°C	27.7	28.9
			Lower Tnight for opt. protein synthesis	°C	23.8	23.6
GT	°C	Tavg, RHmax, Daysff., DVS	Lower Tavg for max GT increases in response to T	°C	19	19
			Higher Tavg for max GT increases in response to T	°C	23	23
			RHmax threshold above which the GT response to T is linear	%	73	90
BDV	BU	Tavg, DVS	Average BDV	BU	613.2	597.5
			Lower Tavg for max BDV increases in response to T	°C	19	19
			Higher Tavg for max BDV increases in response to T	°C	23	23
PV	BU	Tmin, Tmax, Daysff, DVS	Daysff at which the sensitive period to temperature ends	d	20	20
			Higher opt. Tavg for granule-bound starch synthase I	°C	24	24
			Lower opt. Tavg for granule-bound starch synthase I	°C	33	33
SV	BU	Tavg, Daysff, DVS	Maximum rate of amylose accumulation	% d−1	0.55	1.16
			Sensitivity to sub-opt T for amylose synthesis (0 = max;1 = none)	—	0.55	0.5
			Daysff at which sensitive period to T begins	d	15	4
			Daysff at which sensitive period to T ends	d	34	17
			Tmax for amylose synthesis	°C	37.4	39.9
			Tmin for amylose synthesis	°C	16.3	16.5
			Higher T for opt. amylose synthesis	°C	26.2	28.7
			Lower T for opt. amylose synthesis	°C	21.6	25.0
Milling quality
HRY	%	Tnight, sunriseH, sunsetH, rainD, WS, DVS	Higher Tnight for opt. starch synthesis	°C	26.5	26.8
			Lower Tnight for opt. starch synthesis	°C	25.5	25.8
			Tmin for starch synthesis	°C	17	22.9
			Windy days threshold for damage	day	10	10
			WS speed threshold for damage	m s−1	2	2
			Cumulative rain threshold for damage	mm	7	14
			HRY reduction due to critical humidity, wind, and T	%	5.4	1.7
BK	%	Tmin, Tmax, ET0, RHmin, RHmax, rainD, WS, Daysfh, DVS	Average BK	%	8.9	6.0
			WS speed threshold for cracking	m s−1	1.7	1.7
			Humidity threshold for cracking	%	50	50
			Daysfh at which sensitive period to T ends	d	10	10
MWG	%	Tavg, Daysfh, DVS	Tavg above which damage occurs	°C	22	23
			Daysfh at which sensitivity to T progressively decreases	d	8; 25; 40	11; 29; 40
			Sensitivity to temperature stress (0 = none; 1 = max)	—	0.6; 0.39; 0	0.59; 0.30; 0
PG	%	Tmin, Tmax, RHmin, RHmax, Daysfh, DVS	Average PG	%	8.2	3.4
			Thermal excursion threshold for PG damage	°C	10	10
			Daysfh at which sensitive period to weather conditions ends	d	20	20

Note:

Abbreviations: AC, amylose concentration; BDV, breakdown viscosity; BK, broken kernels; BU, Brabender Units; GT, gelatinization temperature; HRY, head rice yield; MWK, milky white grains; PC, protein concentration; PG, pecky grains; PV, peak viscosity; SC, setback viscosity.

1 = emergence, 2 = flowering, 2.5 = full ripening, 3 = physiological maturity); Days_ff = days from flowering (days); Days_fh = days from heading (days); DVS = development stage code (unitless); ET₀ = potential evapotranspiration (mm); Rad = daily global solar radiation (MJ m⁻²); rain_D = daily rain (mm); RH_max maximum air relative humidity (%); RH_min = minimum air relative humidity (%); sunrise_H = sunrise hour; sunset_H = sunset hour; T = Air temperature; T_avg = average air daily temperature (°C); T_max = maximum air daily temperature (°C); T_min = minimum air daily temperature (°C); T_night = hourly night temperature (°C); WS = daily wind speed (m s⁻¹).

QV = Quality variable; BU = Brabender Units, arbitrary unit to measure the viscosity of a slurry (Karim et al., 2000).

Models assumptions

The assumptions in the modeling study presented here entail the choice of the time step, and the model structure and parameters. The temporal resolution of the rice quality models is daily or even higher, although the simulated biological processes can occur at sub-daily time resolution (e.g., 1 h or even less). We used semi-empirical models driven by meteorological inputs, using parameters to differentiate cultivar-specific sensitivity to environmental stresses during ripening (Table 2). The complex biological processes involved in rice quality formation and their mutual interactions as affected by pre- and post-harvest management operations are not explicitly considered. Grain quality is simulated using ten models which do not consider any feedback between variables and related processes. Grain length, size, and shape are important parameters determining grain aspect and milling quality which are not explicitly accounted for by our modeling approaches. Also, the effect of low radiation in exacerbating rice grains susceptibility to temperature damage is not implemented. Despite these assumptions, the ten rice models require inputs that are commonly available in agrometeorological databases and demonstrated to be capable to explain about 60%–90% of inter-annual rice quality variability when calibrated using reliable reference data. The parameter sets developed in this study are valid for sub-groups of varieties that share similar grain quality traits with Loto and Gladio cultivars. The application of rice quality models outside the conditions for which they are developed (e.g., areas with different pedo-climatic conditions) and/or to other cultivars should be preceded by a novel calibration, which requires multi-year experiments where plant phenology and grain quality are measured as reference datasets.

Input data for the spatially distributed simulations

A georeferenced database including information on rice spatial distribution (Figure 1a), sowing dates (Figure 1b), current climate (Figure 1c,d), and future scenarios was built at 2 × 2 km spatial resolution to reproduce the spatial heterogeneity of the study area. The input data were univocally assigned to each 2 × 2 simulation unit depending on spatial attributes (i.e., geographic coordinates of centroids) via the use of GIS-based software (Ginaldi et al., 2019). The meteorological data used as input to the rice quality models are presented in Table 3.

TABLE 3 Meteorological data as input to the rice quality models

[Table omitted. See PDF]

Note:

Res = Resolution; O = Original; S = This study; H = Hourly; D = Daily; ET0 = reference evapotranspiration. Nighttime air temperatures were defined as the hourly temperatures between sunset and sunrise hours.

References: 1 = www.ecmwf.int; 2 = Pagani et al. (2019); 3 = https://esgf-node.llnl.gov/search/cmip5/; 4 = Donatelli et al. (2009); 5 = Campbell (1985); 6 = Spitters et al. (1986); 7 = Hargreaves and Samani (1982).

Daily maximum and minimum air temperature (°C), precipitation (mm), wind speed (m s⁻¹), and reference evapotranspiration (mm) in the period 2003–2014 were retrieved from the database of the ERMES project (http://www.ermes-fp7space.eu/it/homepage/; Pagani et al., 2019). Daily global solar radiation (MJ m⁻²) was estimated as a function of the daily thermal excursion according to Hargreaves and Samani (1982), whereas daily maximum and minimum air relative humidity (%) and hourly air temperatures were estimated via the CLIMA software library (Donatelli et al., 2009). In particular, hourly temperature was estimated starting from daily data based on the equation of Campbell (1985). Nighttime air temperatures were defined as the hourly temperatures between sunset and sunrise hours, which were computed according to Spitters et al. (1986) (Text S2).

Future climatic series were generated by the stochastic weather generator CLIMAK (Danuso, 2002), following the procedure adopted by Fernandes et al. (2012), Cappelli et al. (2015), Paleari et al. (2015), Bregaglio et al. (2017), and Giuliani et al. (2019) (Text S3). Two contrasting Representative Concentration Pathways (RCPs) scenarios (AR5; Stocker et al., 2013), that is, RCP 2.6 (low impact, with radiative increase up to 2.6 W m⁻² and CO₂ concentration up to 420 ppm in 2100) and RCP 8.5 (high impact, with radiative increase up to 8.5 W m⁻² and CO₂ concentration up to 936 ppm in 2100) were considered. In agreement with Bregaglio et al. (2017) and Giuliani et al. (2019), four general circulation models (GCMs) were selected among those produced in the Coupled Model Intercomparison Project (CMIP5, https://pcmdi.llnl.gov/mips/cmip5/; Stocker et al., 2013) to take into account the uncertainty in future weather projections. They were the Norwegian Earth System Model (NOResm, Tjiputra et al., 2013), the Model for Interdisciplinary Research on Climate (MIROC-ESM, Watanabe et al., 2011), the Hadley Centre Global Environmental Model version 2 (HadGEM2-ES, Collins et al., 2011), and the GCM developed by the Goddard Institute for Space Studies (GISS-ES, Schmidt et al., 2006). For each combination of RCP (2) × GCM (4) × grid cell (1345), monthly temperature and precipitation anomalies (i.e., difference between the future and the current value of the weather variables) were computed by means of data available on the Earth System Grid Federation data portal (https://esgf-node.llnl.gov/search/cmip5/) and then used to generate 20-year synthetic weather series centered on 2030 and 2070, leading to a total of 21,520 climatic series. Projected thermal increases (Figure S1) and rainfall decreases (Figure S2) fluctuated around 1.6℃ and −16 mm in 2030 and 2.8℃ and −37 mm in 2070 (Figure 1), and followed a north-west to south-east gradient, in agreement with Cappelli et al. (2015), who assessed the future trends of maize and giant reed productivity in the same area. Information on rice distribution (i.e., area covered by rice in each 2 × 2 km simulation unit) and spatially distributed farmers' sowing dates was derived via multi-temporal analysis of MODIS composites images at 500 m spatial resolution (Boschetti et al., 2009).

Analysis of results

The outputs of spatially distributed simulations were averaged on each grid, considering the 12-years of the baseline and the 20 years corresponding to future scenarios for each combination of GCM × RCP × time frame × cultivar. Models results were used (i) to draw boxplots presenting the variability of the ten quality variables in the baseline and (ii) to evaluate the impact of climate change as percentage difference compared to the baseline. Rice quality variables were then aggregated considering milling (i.e., MS: head rice yield, milky white, pecky and broken kernels) and cooking suitability (i.e., CS, amylose, protein, setback, gelatinization temperature, breakdown, and peak viscosity), in order to give synthetic indicators of their expected future trends. For each combination of GCM × RCP × time frame × cultivar, percentage changes of milling and cooking suitability (S_celli, %, Equation 1) were computed for each grid cell i and used to draw boxplots and maps showing their spatio-temporal variability on the rice district.[Image Omitted. See PDF]where ${\overline{\mathit{QV}}}_{\mathit{jF}}$, ${\overline{\mathit{QV}}}_{\mathit{jB}}$ are the average values (in the unit of the variable) of the j-th simulated quality variable under the 20 years future (F) and 12 years baseline (B) scenarios, respectively; ${\overline{\mathit{QV}}}_{\mathit{jallcell}}$ is the average value of the j-th quality variable across all the 1345 grid cell of the rice district in the baseline scenario (unitless); ± indicates either positive (+) or negative (−) effect of the relative change (with respect to the baseline) of the j-th variable on rice quality, based on the preferences of the Italian consumers. Increases in breakdown, peak viscosity, gelatinization, broken, milky white, and pecky grains were considered negative (vice versa for decreases); increases in head rice yield, setback, protein, and amylose content were considered positive (vice versa for decreases).

The rationale for the aggregation of quality variables was driven by the preference of the Italian consumers for non-glutinous rice cultivars with intermediate amylose and medium to high setback, which tends to become hard, firm, and non-sticky upon cooking, that is, presenting the characteristic “al dente” texture, which is particularly suited for traditional rice dishes preparation as risotto and timbales. Increases in broken, milky white, and pecky grains (for MS) and in breakdown, peak viscosity, and gelatinization temperature (for CS) were then penalized, because high values of these variables are associated with reduced number of whole grains and therefore lower yield, edibility, and market value of rice (Siebenmorgen et al., 2013). The opposite criteria were followed for head rice yield, that is, the relative weight presence of whole grain after milling. Higher breakdown and peak viscosity make the boiled rice softer, more chewy and clingy (Liu et al., 2007), which are negative aspects for domestic consumers. Likewise, high gelatinization temperature causes longer cooking time (Fitzgerald et al., 2009), makes grains more prone to disintegrate when overcooked (Mackill et al., 1996) and negatively affects energy cost, hydration, and quality (e.g., color and yellowness) of parboiling process (Taghinezhad et al., 2016). Conversely, high protein values were considered favorable since they ameliorate the nutritional properties (Liu et al., 2007) while improving rice sensory and cooking quality: Indeed, protein content is negatively correlated with unpleasant characteristics such as grain adhesiveness (i.e., the tendency of the grains to form a compact mass after boiling; Mestres et al., 2011), peak viscosity and breakdown and positively with setback (Champagne et al., 2007).

The internal relationships between quality variables in the current and future climatic scenarios were evaluated with Pearson's correlation analysis (p < 0.01). A contour plot analysis was carried out to inspect the effects of temperature and cycle duration on the variability of milling and cooking suitability under climate change scenarios, given the primary role of these variables in influencing rice quality.

Boxplots were drawn with ggplot2 R package, version 0.9.3 (Wickham, 2009), whereas contour plots were drawn with the Interpolation of Irregularly and Regularly Spaced Data—akima - R package, version 0.6-2 (Akima et al., 2016). Correlation matrices were computed and graphically displayed using the corrplot R package, version 0.84 (Wei et al., 2021).

RESULTS Simulated rice quality variables in the baseline scenario

The variability of the ten quality variables simulated for Loto and Gladio in the baseline scenario is presented in Figure 2.

FIGURE 2. Boxplots of quality variables simulated in baseline climate (2003–2014) for Loto and Gladio cultivars across the Lombardo-Piemontese district. Each boxplot consists of 1345 data corresponding to the average simulated value for each grid cell in the baseline period. Boxes extend from 25th to 75th percentiles (interquartile distance, IQD), centerline is set at 50th percentile, whiskers extend up to the extreme data point which is 1.5 times higher/lower than the length of the box away from the box; outliers are represented as points. AC, amylose concentration; BDV, breakdown viscosity; BK, broken kernels; BU, Brabender Units; GT, gelatinization temperature; HRY, head rice yield; MWG, milky white grains; PC, protein concentration; PG, pecky grains; PV, peak viscosity; SC, setback viscosity

Considering CS, Gladio results showed higher and less variable average amylose content (26.28%, IQD of 0.65%), setback (765.78 BU, IQD of 11.41), and gelatinization temperature (78℃, IQD of 0.22℃) than Loto (AC = 15.68%, IQD of 1.37%, SV = 634 BU, IQD of 11.67 BU, GT = 70.72, IQD of 0.32℃). Breakdown (mean of 637.31 BU, IQD of 37.72 BU for Gladio; mean of 732.33 BU, IQD of 114.94 BU for Loto) and peak viscosity (mean of 1076.85 BU, IQD of 20.32 BU for Gladio; mean of 1198.26 BU, IQD of 45.34 BU for Loto) were lower and more uncertain. Despite the average protein content was very similar in the two cultivars (mean of 6.27% for Gladio and of 6.22% for Loto), its variability in Gladio (IQD of 0.36%) was three times higher than in Loto (IQD of 0.11%).

Regardless of the variety, simulated amylose (mean IQD of 1.01%) was more variable than protein (mean IQD of 0.24%) whereas setback (mean IQD of 11.54 BU) was the less uncertain variable among those connected to starch viscosity profile (14.94 < IQD < 45.34 BU).

The mean percentage of damaged grains was always higher in Loto (mean of 7.05%) compared to Gladio (mean of 3.35%), with highest values for broken kernels (mean of 10.73% for Loto and 5.69% for Gladio), followed by pecky (mean of 6.63% for Loto and 3.78% for Gladio) and milky white grains (mean of 3.79% for Loto and 0.57% for Gladio). The same consideration holds for the related variability (0.50 < IDQ < 1.80% for Loto and 0.06 < IDQ < 0.82% for Gladio), with pecky (0.82 < IQD < 1.80%) and milky white grains (0.06 < IQD < 0.52%) resulting as the most and less variable quality variables, respectively. Head rice yield assumed average values of 63% for Loto and 64% for Gladio, with similar mean IQD for both varieties (mean of 1.39% for Loto and 1.51% for Gladio). Correlations among the ten quality variables in the baseline scenario are presented in Figure 3 as a correlation matrix.

FIGURE 3. Linear correlations between simulated rice quality variables in the baseline. Each boxplot consists of 1345 values data achieved by corresponding to the averaging the average simulated data value for each grid cell across in the baseline period. Significant correlations (p [less than] 0.01) are coloured either in blue (positive) or red (negative) hues, while not significant correlations are represented in white. The size of each coloured circle is proportional to the value of the Pearson's correlation coefficients (r). AC, amylose concentration; BDV, breakdown viscosity; BK, broken kernels; GT, gelatinization temperature; HRY, head rice yield; MWK, milky white grains; PC, protein concentration; PG, pecky grains; PV, peak viscosity; SC, setback viscosity

The analysis of results revealed 42 significant values out of 45, with different strength: 28 Pearson r values were higher than 0.6 (strong correlation), nine values ranged between 0.2 and 0.6 (moderate correlation), and the last eight r values were lower than 0.2 (weak correlation). Strong correlations emerged on variables with effect on both milling and cooking properties, whereas weaker correlations were found with head rice yield and protein content with respect to the other quality variables.

Considering future simulations, in 2030, all the correlations were significant (Figure S3a), that is, also concerning the relationships of protein with milky white grains (r value of 0.71 instead of −0.02), peak viscosity (r value of 0.6 instead of 0.03), and pecky grains (r value of −0.72 instead of 0.01). Pearson r values between quality variables in 2070 were very similar than in 2030 (Figure S3b), with marginal deviations (<0.1 Pearson r) in 34 out of 45 correlations.

The percentage differences of rice quality variables in 2030 and 2070 with respect to the baseline are presented in Table S2, separately for the different combinations of GCM × RCP.GISS GCM and RCP 2.6 led to the lowest quality decline in most combinations of cultivar × quality variable × time frame, whereas the higher decreases were obtained in the combination HadGEM and RCP 8.5, especially in 2070. Loto showed higher values and variability compared to Gladio, and differences between the two cultivars followed the trend RCP 2.6 > 8.5, 2030 > 2070, GISS >HadGEM. Among the variables related to milling yield, milky white grains showed the highest increase (16%–27% in 2030, 14%–62% in 2070 for Loto, 4%–6% in 2030, 3%–6% in 2070 for Gladio), whereas head rice yield denoted the most marked decline (−3% to −7%, −2% to −20% for Loto in 2030 and 2070; −4% to −7%, −3% to −15% for Gladio in 2030 and 2070). The variations of broken and pecky grains in the future were below ±3%. Considering cooking suitability, major differences were simulated for setback, with divergent trend between cultivars, positive on Gladio (mean increase always in the range 0%–2%) and negative on Loto (mean of −3% in 2030 and of −6% in 2070). Other differences were observed on protein content (mean change of −1% in 2030 and −4% in 2070 for Loto; mean change of −6.5% in 2030 and −10% in 2070 for Gladio), breakdown (mean change of +3% in 2030 and +4.5% in 2070 for Loto; mean change of +7% in 2030 and +8% in 2070 for Gladio), and peak viscosity (mean change of +5% in 2030 and +7.5% in 2070 for Loto; mean change of +3.5% in 2030 and +5% in 2070 for Gladio). Expected future changes in amylose and gelatinization temperature were minimal.

Future trends in milling and cooking suitability

The percentage changes of milling and cooking suitability in future weather scenarios with respect to the baseline, considering the two rice cultivars (Loto and Gladio), the two RCPs (2.6 and 8.5), and the four GCMs (GISS, HadGEM, MIROC, and NOResm) are shown as boxplots in Figure 4.

FIGURE 4. Simulated percentage changes in milling and cooking suitability for Loto (light blue boxes) and Gladio (coral boxes) in future simulations with respect to the baseline. Data are distinctly presented for the two Representative Concentration Pathways (RCP, rows) 2.6 and 8.5, the four GCMs (x-axis; GIS = GISS, Nrs = NOResm, Mrc = MIROC, Hdl = HadGEM), and the two time frames (2030 and 2070, columns). The boxplots correspond to the variability in the 1345 grid cells corresponding to the rice area in the Lombardo-Piemontese district

For Gladio, simulated future milling suitability was stable across RCP, time frame, and GCM, and percentage changes with respect to the baseline ranged from 0% to +5% with peaks of ±10%. Loto simulations depicted different results, with a decline of milling suitability in 2030 (−2.5% to −10%), with the following trend: GISS < NOResm < MIROC < HadGEM and RCP 2.6 < RCP 8.5. This pattern was exacerbated in 2070, with average milling suitability declines comprised between −2.5% and −22.5%, with maximum value of −30% in the most extreme scenario (i.e., HadGEM 8.5). The variation of cooking suitability in 2030 was negative (0 and −5%), with similar values between cultivars, GCM, and time frames. In 2070, the cooking suitability decline followed the same pattern outlined by milling suitability in Loto, varying from 0 to −2.5% for both varieties in GISS 2.6, to −5% (Gladio), and to −12.5% (Loto) in HadGEM, RCP 8.5.

IQDs of milling suitability varied in the range 3.80%–6.32% (average 4.73%) and in the range 0.89%–1.57% for cooking suitability (average 1.17%). While the variability associated with Loto, RCP 8.5 and 2070 time frame fluctuated around a mean IQD value of 3.05% (0.89 < IQD < 6.32%), the uncertainty connected to Gladio, RCP 2.6 and 2030 time frame ranged between 0.92% and 5.28%, with a mean IQD of 2.85%.

The variability due to GCM progressively increased from GISS (0.89 < IQD < 5.46%, average 2.81%), to NOResm (0.92 < IQD < 5.53%, average 2.89%), to MIROC (0.92 < IQD < 5.34%, average 2.96%), and to HadGEM (0.92 < IQD < 6.32%, average 3.14%).

The combined effect of air temperature and length of simulated crop cycle on the variability of milling and cooking suitability was further investigated in both cultivars using contour plots; results are presented and discussed in Figure S4.

Spatio-temporal variability of rice milling and cooking suitability

The percentage differences of simulated milling and cooking suitability in the future with respect to the baseline are presented as maps in Figure 5, considering the two varieties Loto and Gladio and the two combinations GCM × RCP leading to extreme results (low impact scenario, GISS, RCP 2.6; high impact scenario, HadGEM, RCP 8.5).

FIGURE 5. Percentage differences of milling (a, b, c, d, Loto; e, f, g, h, Gladio), and cooking suitability (i, j, k, l, Loto; m, n, o, p, Gladio) in 2030 and 2070 with respect to the baseline. Results are presented for the two extreme scenarios GISS 2.6 (low impact) and HadGEM 8.5 (high impact) as percentage difference compared to baseline

Milling suitability simulations revealed a general decline over the rice district, with large heterogeneity according to the cultivar and time frame (Figure 5a–h), whereas cooking suitability remained almost stable, without large differences between cultivars, except for the high impact scenario in 2070 (Figure 5i–p).

In 2030, milling suitability for Loto decreased along a north-west to south-east gradient, according to higher thermal increases, with peaks of −15% for the high impact (Figure 5c) and −13% for the low impact scenario (Figure 5a). The only increases were concentrated in the coolest and rainiest areas of the Western part of the district—that is, with average daily temperature lower than 18.5℃ and cumulative precipitation of about 415 mm in the period March-September in the baseline. The simulations in 2070 depicted the same trend, with major declines in the areas characterized by highest thermal anomalies (−28% for high impact scenario; −11% for low impact scenario).

The spatial representation of milling suitability for Gladio (Figure 5e–h) did not show any clear pattern along the study area, with negative values in the western part and null/positive changes in the central and eastern part of the district.

Compared to milling suitability, the spatial representation of cooking suitability in 2030 highlighted less variability, with major declines in the range −0.1% to −6% regardless of the combination of cultivar ×scenario (Figure 5i,k,m,o). The higher thermal increases projected in 2070 led to a general worsening of cooking suitability with respect to the baseline (Figure 5j,l,n,p), with average value of −7%, and peaks of −12% for Loto in the high impact scenario.

DISCUSSION Consistency of simulated quality variables in the baseline scenario

The results of the simulations in the baseline scenario were consistent with the available information on the variety-specific quality variables in the study area (Table 1; Cappelli, Confalonieri, et al., 2018), confirming that the upscaling from field to large scale did not alter models validity. Gladio obtained high amylose, gelatinization temperature, and setback, whereas Loto showed higher breakdown, peak viscosity, and higher susceptibility to broken, milky white, and pecky grains formation. Average simulated values of protein content and head rice yield were very similar in the two cultivars. The correlation analysis revealed the internal relations between the ten simulated quality variables and provided a multifaceted evaluation of rice grain quality besides focusing on single (i.e., protein, Li et al., 2005; amylose, Chen et al., 2011; chalky grains, Masutomi et al., 2015) or partial groups (Lanning et al., 2012; Shen et al., 2007) of quality aspects. We conducted this analysis because the quality variables in our study were simulated without considering interactions and feedbacks neither between variables nor between processes (e.g., protein effect on reducing water availability for starch gelatinization and the impact of altered rates of starch biosynthesis/deposition on the incidence of grain damages are not implemented). Correlation analysis showed significant relationships between variables connected to milling and cooking suitability, and our results were consistent with previous studies, with few exceptions. In agreement with literature, amylose content was positively correlated with setback (Pang et al., 2016), gelatinization temperature (Varavinit et al., 2003), and head rice yield (Yun et al., 2015) and negatively with breakdown and peak viscosity (Pang et al., 2016). Consistent relationships emerged also for the indices related to starch viscosity properties, except for the relationship between gelatinization temperature and setback, which is contrasting in available studies (significant positive, Yun et al., 2015; non-significant negative, Pang et al., 2016). In our study, breakdown was positively associated with peak viscosity and negatively to setback and gelatinization temperature; peak viscosity was negatively correlated with setback and gelatinization temperature, and gelatinization temperature positively related to setback.

As expected, head rice yield was positively influenced by protein content (Balindong et al., 2018) and negatively by milky white and broken grains, which both make grains more liable to breakage during milling processes (Siebenmorgen et al., 2013). Unexpectedly, head rice yield was significantly positively correlated with pecky rice, which is often structurally damaged and more breakable than whole and healthy grains (Wang et al., 2002). This inconsistent correlation may be the result of the pecky rice data used to calibrate the model or of differences in the duration of the sensitive period to weather conditions. In our study, pecky rice was simulated as a function of daily thermal excursion, relative humidity, and consecutive drought days in 20 days after heading (BBCH 55), whereas the model of head rice was driven by hourly nighttime temperature during late maturity (BBCH 83–89), and decreased when cumulative rainfall and daily wind speed exceeded critical thresholds for starch synthesis after full flowering (BBCH 65–83). However, it must be pointed out that pecky rice is a common definition for many grain imperfections, which can be due to insect bites, diseases and virus infections, nematode attacks, suboptimal thermal/humidity conditions, and high windiness during grain filling (Groth & Lee, 2002; Groth et al., 1991; Wang et al., 2002). Protein content showed consistent significant positive correlation with amylose (Kakar et al., 2019), and weak or no correlations with the other quality variables. This suggests that the balance between C and N metabolisms (Lin et al., 2014) and protein composition (Balindong et al., 2018), rather than total protein content, may be associated with the occurrence of grain damages and/or modifications of starch viscosity during cooking. The weak correlation between head rice yield/protein content and the other quality variables may be also due to the different model input requirements, or to model structure and parameterization. Indeed, differently from the other eight models, the head rice yield and protein models use response functions based on hourly nighttime air temperature, considering a higher number of parameters (Table 2). The average simulated values for Loto and Galdio were very similar with differences lower than 1% for head rice and 0.2% for protein content. These similarities may also explain the significant positive correlation between head rice and protein content (average Pearson r = 0.48).

Simulated changes in rice quality variables in future climate

The future trends in rice quality in Northern Italy are projected to be deeply affected by climate change, with high variability and positive and negative variations according to the variable, climate scenario, and area. In general, Loto resulted more sensitive than Gladio to the warmer and drier conditions in the rice ripening period, due to its higher susceptibility to temperature stress in post-flowering (Table 1; Cappelli, Pagani, et al., 2018). The combination HadGEM GCM and RCP 8.5 was associated with the highest thermal increases, in turn leading to less favorable conditions for grain filling. In this climate scenario, temperatures frequently exceeded the optimum temperature for storage biosynthesis (22–29℃), packing of starch granules, and increment of grain fissuring (22–27℃), especially in 2070. Conversely, the combination GISS GCM RCP 2.6 led to minor changes with respect to baseline climate, and this reflected into a slight negative/positive conditions for starch and protein biosynthesis, thus limiting the occurrence of grain defects, especially in Gladio. The largest decreases in rice quality were projected for milky white grains, head rice yield, and protein content. These results are probably due to the projected higher nighttime/daily average temperatures and to the decrease in precipitations, leading to the non-homogeneous packing of starch granules in the caryopsis and favoring breakages in post-harvest processing. The higher temperature in post-flowering led to the shortening of grain filling and to suboptimal conditions for protein synthesis, thus reducing the N accumulation in grains and its conversion to protein (Cheng et al., 2010; DaMatta et al., 2010). While the projected incidence of milky white grains in Gladio slightly increased (average +4.5%), the Loto results depicted a huge variation (14% and 60%), consistently with climate change studies carried out in Japan using comparable scenarios from AR4 (Okada et al., 2011) to AR5 (Masutomi et al., 2019; Takimoto, Masutomi, & Tamura et al., 2019) IPPC SRES. The magnitude of protein and head rice yield decreases was coherent with experiments carried out under FACE conditions (Cheng et al., 2010; DaMatta et al., 2010; Taub et al., 2008; Yang et al., 2007) revealing an average decline up to 15% and 23.5%, respectively. Amylose content and gelatinization temperature remained stable in future simulations, due to the dominant role of genetics in controlling these traits (Cappelli et al., 2014). For pecky grains, the negative effects of thermal increases were counterbalanced by the decline in daily temperature and humidity excursion, which reduces the rates of moisture adsorption-desorption at the grain surface causing fissures in the kernels (Lan et al., 1999). The changes in the storage composition due to raising temperatures will lead to less favorable pasting properties of rice flours, by increasing peak viscosity, breakdown, and gelatinization temperature. These results suggest that, behind amylose modification, changes in the starch structure of the caryopsis could be connected to its viscoelastic and rheological properties. Indeed recent studies showed that the elevated thermal regimes during ripening decreased amylopectin branching (Lanning et al., 2012), also favoring the links between protein and amylose (Hamaker & Griffin, 1990) and between lipid and amylose (Liang et al., 2002), as well as starch pre-gelatinization due to the increased activity of α-amylase (Yamakawa et al., 2007). Consistently with experimental evidences, our simulations denoted greater gelatinization temperature, peak viscosity, breakdown, and lower setback in future climate, thus reducing texture and sensory attributes and increasing temperature requirements to gelatinize starch.

Future trends of rice milling and cooking suitability

The aggregation of the rice quality variables in milling and cooking suitability according to consumer preferences allowed synthesizing model outputs and simplifying the interpretation of our results. The largest impact of climate change on milling suitability was due to the high number of quality variables affected by temperature changes. The corresponding models indeed use daily response functions reducing the cultivar-specific susceptibility as the ripening progresses. Cooking suitability was less affected by climate change, showing very narrow relative percentage variations (mostly in the range 0%–5%) that were almost constant across the study area, scenario, and time frame considered (Figure 5). Besides the prevalent role of genetics versus environmental factors in controlling cooking quality indicators, the less responsiveness to climate anomalies was partly due to the fact that the related models do not consider the effect of precipitation changes on cooking quality (except for gelatinization temperature; Table 2) and implement static heat stress factors computed at the end of the rice season. Indeed, averaging temperature responses over the entire ripening instead of using daily or short period factors (e.g., 10 or 20 days after heading/flowering) as in the case of milling suitability, leads to smooth the effects of daily abrupt changes and/or extremes in temperature on quality variables determining cooking suitability trends.

The spatial representation of rice quality simulations revealed that the milling suitability progressively declined along a north-west to south-east gradient for Loto, according to the temperature increases and rainfall decreases. This generalized decline suggests that current thermal conditions are already close to the optimum for this cultivar, with few exceptions in the eastern and south-eastern areas. Conversely, Gladio results showed a different pattern, with positive changes in milling suitability in the central and eastern part of the district, and lowest values in the western areas. Future cooking suitability was much less affected in our simulations, with the exception of the combination HadGEM GCM and RCP 8.5 in 2070, thus confirming a predominant role of genetics versus environmental factors in regulating starch-related indicators (Siebenmorgen et al., 2013).

Our simulation study could be used to make economic projections on the future profitability of rice farmers, millers, and buyer businesses, because rice quality deeply influences the price of rice productions, with damaged grains redirected toward livestock feed or pet foods producers (Paranthaman et al., 2009). According to current prices of milled rice in the local market (approximately 650 € t⁻¹ for Gladio, 780 € t⁻¹ for Loto, www.enterisi.it, updated 18 February 2020), the projected decrease in milling suitability for Loto could translate into a reduction in the range 38 € t⁻¹/63 € t⁻¹ in 2030 and 32 € t⁻¹/154 € t⁻¹ in 2070. Considering the highest impact climatic scenario, changes in rice price could reach peaks of −117 € t⁻¹ in 2030 and −220 € t⁻¹ in 2070. Gladio results on the contrary may lead to average revenue increases of about 11 € t⁻¹ in all scenarios (with peaks of 60–65 € t⁻¹), except for the combination of HadGEM GCM and RCP 8.5 in 2070, in which the rice price is projected to slightly decrease by 6 euros € t⁻¹, with peaks of about −58 € t⁻¹. Considering cooking suitability, average millers income may be reduced of −15 € t⁻¹ in 2030 and −2 € t⁻¹ in 2070 regardless of the cultivar and climate scenario, with the exception of the combination Loto × HadGEM GCM × RCP 8.5 in 2070, in which projected price decrease was −72 € t⁻¹. Accordingly, farmers' income could also be reduced following the same trends, although changes will be smaller since the reference rough rice price is nearly half of the milled rice (i.e., 2.4 times lower).

Perspectives and limitations of the study

The results of this study are based on the European rice production context and cannot be generalized to other areas, because rice quality and related standards, end uses and consumer preferences largely vary across countries (Unnevehr et al., 1992). However, the consideration of globally recognized characteristics of rice milling and cooking quality makes our procedure adaptable to markets with very different standards. For instance, in Japan consumers prefer soft, sticky, short, rice with low amylose (i.e., 2%–19%), low protein content (i.e., <5%), and percentage of fissured and chalky grains <4% and <0.1%, respectively (Calingacion et al., 2014; Lee, 1998; Mutters & Thompson, 2009). Conversely, in China consumer preference focuses on soft short bold (northern and south-western provinces) to soft long slender (provinces in the central and southern China) rice with low amylose (i.e., 2%–19%), low to intermediate gelatinization temperature (i.e., <70–74℃) and high protein content (i.e., at list 7% for Japonica and 8% for Tropical Japonica types cultivars) (Calingacion et al., 2014; Juliano & Villareal, 1993). In any case, the selection of the rice quality variables and their effect on milling and cooking suitability, must be preceded by the calibration of the rice quality models according to the targeted environmental and management conditions. This requires the availability of field experiments where main phenological phases are observed (i.e., heading, flowering, and physiological maturity) and rice quality variables are collected. This currently represents the major limitation to a further extension of our procedure even in Northern Italy, where some pilot field trials have been conducted only on Selenio, a widespread round grain cultivar with high susceptibility to breakage and pecky grain formation (Miniotti et al., 2016).

We did not consider in our simulation study any water and nitrogen stress on rice growth and quality, according to standard rice-growing conditions in Northern Italy, where the crop is grown under continuous flooding with intensive use of nitrogen and chemical inputs. Where these assumptions do not hold, additional data will be needed for models execution, possibly considering experiments where rice quality is collected under different irrigation and fertilization strategies. This will also require the implementation of the rice quality models to consider these factors.

We did not consider in our study the possible implementation of adaptation strategies (Ishimaru et al., 2016), intended as both short- (e.g., shifting of sowing dates) and long term (e.g., adoption of high temperature tolerant cultivar). Such a work would have required the definition of the most promising traits for rice breeding, and the establishment of a tighter link between model parameters and related traits. Preliminary studies on model-based rice ideotyping in the same area (Paleari et al., 2017, 2019) indicate that variations of the parameters connected to temperature-related traits affecting chalkiness and grain breakage did not influence head rice yield and milky white grains when applied under current and climate change scenarios. This suggests that the impact of environmental factors could be larger than the one of the genotype. In any case, this analysis should be extended to other traits and qualitative aspects to give useful indications to breeders. Nevertheless, the significant relationships between quality variables simulated under current and climate change conditions could be used to support breeding programs. Indeed, the proxy tests that are usually performed to analyze rice grain quality preferences mainly rely on the estimation of three quality indicators: amylose content, gelatinization temperature, and gel consistency (Butardo et al., 2019). However, samples in the same amylose content class could have different sensory profiles which result in markedly different cooking and eating quality characteristics (Cuevas et al., 2018). Additional tests that characterize starch structure, paste viscosity, and grain texture are needed in order to gain a deeper understanding of rice cooking quality. These tests are time-consuming and expensive, and are thus not routinely used to develop breeding targets for specific markets (Butardo et al., 2019). In this context, rice quality models emerge as time- and cost-effective tools to perform integrated a priori evaluations of multifaceted aspects of rice grain quality in a changing climate, due to their ability in reproducing crop responses to boundary and unexplored conditions at an optimal spatial/temporal resolution. Site-specific model-based simulation experiments may support breeders in designing the best-suited selections programs to be tested in vivo with the aim of improving the milling and cooking quality in response to evolving consumer preference, market demand, and changing climate.

Our projections may be slightly underestimated because we did not consider the negative effects of the increased atmospheric CO₂ concentration on rice quality variables (Jing et al., 2016). These effects were not simulated because of the lack of information of the CO₂ responses on Italian cultivars, being those limited to some Japanese and Chinese cultivars. Recent studies demonstrated that rice grown under enriched CO₂ environments markedly reduced protein content and head rice yield, while significantly increasing milky white grains percentage, breakdown, and peak viscosities (Jing et al., 2016; Usui et al., 2014, 2016; Yang & Wang, 2019). The only positive changes will involve the decrease in the gelatinization temperature, which would result in decreased cooking time and improved color and yellowness of the parboiled rice (Taghinezhad et al., 2016). Another critical point in our approach is the formalization of the models for milling suitability, which currently do not consider different key factors in the simulation of grain imperfections. Despite grain length, size, and shape are important parameters determining grain aspect and milling quality, they are not explicitly considered in the rice quality models. In this context, recent studies reported a significant negative correlation between grain length and head rice yield, whereas grain width and length–breadth ratio are poorly associated with milling quality (Suman et al., 2020; Xie et al., 2013; Zhou et al., 2015). In our study, the varietal diversity in grain morphology and related influence on milling quality is reproduced via model parameterization, that is, changing values of those parameters describing variety-specific susceptibility to the formation of milky white and broken kernels (e.g., average BK, average PG, HRY reduction due to critical humidity, wind, and T; T_avg for MW damage and related sensitivity to temperature stress; Table 2). In this context, Gladio, which has an average grain length of 7.74 mm and is generally lower sensitive to grain imperfections, is characterized by lower values of all these parameters compared to Loto, which has an average grain length of 6.43 mm. Post-harvest operations, such as drying, storage, and milling, play a key role in determining rice milling yield and overall quality. Drying is a post-harvest operation aimed at reducing moisture of paddy rice to a suitable level, in order to minimize respiration rates and fungi/insect growth during storage (Tong et al., 2019). Drying and milling conditions cause abiotic and mechanical stresses in rice kernels, which may result in kernel fissuring and breakage and consequently lowering head rice yield. Nonetheless, we decided to focus our analysis on the pre-harvest period for four main reasons: (1) pre-harvest factors such as environmental conditions during growing season have crucial effects on rice grain quality, especially during ripening; (2) climate change impacts on crop development, growth, and quality occur during growing season; (3) available pre-harvest quality models can be easily coupled with crop models, allowing for a dynamic simulation of crop phenology and quality traits in current and future climate scenarios with relatively few inputs and model parameters, (4) models for post-harvest quality (e.g. Bunyawanichakul et al., 2007; Franco et al., 2020) require engineering and process background, are completely isolated from the outside environment (i.e., do not consider the effect of pedo-climatic conditions) and are too sophisticated to be effectively coupled with rice models (consist of chains of interdependent differential equations that require huge number of input factors—variables and parameters—and measured data for calibration). Lastly, Takimoto, Masutomi, Tamura, Nitta, et al. (2019) demonstrated that the simulation of chalky grains is more accurate when both air temperature and solar radiation are considered together and that low radiation levels can exacerbate the sensitivity of rice grains to temperature damage even in cultivars ranked as intermediate heat tolerant. This consideration is even more important in view of the models' application in a climate change impact assessment, since increasing air pollution and aerosol concentration are affecting the amount of solar radiation reaching the Earth's surface (Stanhill & Cohen, 2001; Thakur, 2018), causing the so called “global dimming” phenomenon.

CONCLUSIONS

Few studies on the prediction of the future trends of crop quality under climate change scenarios have been carried out so far. We focused here on rice quality, considering the multiple aspects determining rice milling and cooking suitability, considering two widely cultivated varieties in the main European rice-producing area and an ensemble of climate change projections. Our study indicates a general decline of milling suitability with a large heterogeneity between cultivars and spatial areas, whereas cooking suitability is predicted to remain almost stable in the future. Despite the explicit limitations and assumptions, our work provides plausible indications on rice quality trends and connected economic projections, which can be of interest for different stakeholders to preserve and enhance the sustainability of European rice systems. Therefore, the model-based approach herein presented should be considered as a reproducible procedure that enables quantifying shifts in the vulnerability of rice grain quality as a result of climate change and evolving local/regional consumer preferences, breeding targets and needs of crop/food scientists, food processors, crop growers, and market. In this light, future model developments should focus on a more detailed and explicit formalization of the genetic regulation of grain (i) length/size, (ii) composition, and (iii) susceptibility to heat stress damage during ripening (i.e., milky white, broken, and pecky grains formation), and their mutual dynamic interactions with pedo-climatic and management factors.

ACKNOWLEDGEMENTS

This study received financial support from (i) the EC Seven Framework Programme (FP7-SPACE-2013-1) under Grant Agreement No. 606983, 2013–2016. An Earth obseRvation Model based RicE information Service (ERMES, http://www.ermes-fp7space.eu/it/homepage/) and (ii) the Italian Ministry of Agricultural, Food and Forestry Policies (MiPAAF) under the AgriDigit-Agromodelli project (DM n. 36502 of 20/12/2018).