Abstract
Tropospheric ozone (O3) significantly reduces rice yield. Koshihikari, a leading Japanese rice cultivar, has been recognized as an O3 tolerant cultivar; however, the mechanisms underlying its high O3 tolerance remain unknown. Therefore, we aimed to elucidate the genetic and physiological mechanisms underlying high O3 tolerance in Koshihikari. A series of chamber experiments were conducted to examine photosynthesis, growth, yield-related traits, and gene expression profiles under chronic O3 conditions in Koshihikari and Takanari, with contrasting O3 tolerance. Koshihikari showed no reduction in unhulled-grain weights due to higher total dry weight under chronic O3 conditions, whereas Takanari showed significantly lower grain weights. The high O3tolerance in Koshihikari was attributed to its highly stable photosynthetic performance and increased biomass allocation to its leaf blades. RNA-seq and gene co-expression network analyses revealed that the genes involved in photosynthesis and carbohydrate metabolism are associated with contrasting O3 tolerance.OsRbcS3 , encoding the RuBisCO small subunit, showed contrasting expression profiles between the two cultivars. In Koshihikari,OsRbcS3 was identified as a hub gene candidate in the gene co-expression network, which was highly correlated with photosynthetic performance. These results suggest that OsRbcS3 plays a key role in the genetic mechanisms underlying the high O3tolerance of Koshihikari.
Introduction
Tropospheric ozone (O3) is a major air pollutant that negatively impacts crop growth and yield (Fiscus, Booker & Burkey 2005; Ainsworth, Yendrek, Sitch, Collins & Emberson 2012). While tropospheric O3 concentrations ([O3]) exhibit spatial and temporal variability, most agricultural regions are exposed to elevated [O3] during the crop growing season, leading to a reduction in global crop production (Ainsworth 2008, 2017). In China, where the O3 is the major air pollutant, the relative yield losses in wheat (Triticum aestivum L.), rice (Oryza sativa L.), and maize (Zea mays L.) have been estimated to reach 33%, 23%, and 9%, respectively, during last decades (Feng et al. 2022). Tropospheric [O3] has more than doubled since the Industrial Revolution (Monks et al. 2015), and a further increase is predicted in Asia, the world’s largest producer of rice (Ainsworth, 2017). On this background, improving O3 tolerance in rice is a critical challenge for a secure and sustainable food production (Frei 2015). In field environments, chronic O3 stress (typically <150 nmol mol-1) can mainly occur during the crop growing season, which causes visible injury to rice in the short term and yield reduction in the long term (Frei 2015). Therefore, it is important to elucidate the mechanisms underlying tolerance to chronic O3 stress to develop a practical breeding strategy for rice.
Natural genetic variations in O3 tolerance have been observed in model plants (Brosché et al. 2010) and several crop plants such as maize (Choquette et al. 2019), rice (Sawada & Kohno 2009), soybean (Betzelberger et al. 2010), and wheat (Zhuet al. 2011). In rice, certain cultivars showed no yield reduction under chronic O3 conditions (Sawada & Kohno 2009). Genetic analyses such as quantitative trait loci analysis and genome-wide association studies using mapping populations have identified the chromosomal regions associated with the natural genetic variations in leaf bronzing (Ueda et al. 2015) and grain yield (Tsukahara, Sawada, Matsumura, Kohno & Tamaoki 2013) of rice under elevated O3 conditions. Tsukahara et al. (2015) indicated the involvement of the ABERRANT PANICLE ORGANIZATION 1in genotypic variations associated with O3-induced yield reduction in rice. Importantly, the pyramiding of two quantitative trait loci, OzT8 and OzT9 , resulted in improved O3 tolerance in terms of biomass and grain yield in rice (Wang et al. 2014). Although these results indicate the importance of utilizing genetic variations in O3tolerance, the associated genetic and physiological factors have been identified in only a small number of cultivars. For example, Koshihikari, a leading Japanese rice cultivar, has been recognized as an O3 tolerant cultivar; however, the mechanisms underlying its high O3tolerance remain unknown (Cho et al. 2013a). An understanding of the mechanisms underlying natural variations in such untapped genetic resources is essential to breed rice for higher O3tolerance.
Under elevated O3 conditions, the downregulation of gas diffusion and biochemical processes results in reduced photosynthetic performance and biomass production in crops (Frei 2015; Ainsworth 2017). O3 induces rapid stomatal closure within minutes to hours and prevents additional O3 uptake by plants (Kollist et al. 2007; Vahisalu et al. 2010). O3-induced stomatal closure reduces the CO2 assimilation rate (A ), because the conductance of gas diffusion via the stomata (g s) is a major limiting factor for leaf photosynthesis (Farquhar & Sharkey 1982; Ainsworth 2017). The absorbed O3 is immediately decomposed into reactive oxygen species (ROS) such as hydrogen peroxide (H2O2), superoxide anion (O2-), hydroxyl radical (· OH), and nitric oxide (NO) in the apoplast (Hoigne and Bader 1975). ROS damage membrane lipids, proteins, and nucleic acids and induce programmed cell death, causing visible injury to the plant body (Fuhrer and Booker 2003). In addition, ROS destroy photosynthesis-related pigments and proteins, which lead to the downregulation of photosynthesis. Thus, differences in the response of gas diffusion and biochemical processes could be responsible for the natural variations among rice cultivars in the O3tolerance of photosynthesis and biomass production (Wilkinson, Mills, Illidge & Davies 2012).
Genome-wide transcriptional profiles revealed by microarray and RNA-seq analyses help elucidate the genetic factors associated with plant responses to the environment. Previous studies using transcriptome analysis identified O3-responsive genes associated with various biological processes in plants (Miyazaki et al. 2004; Liet al. 2006; Ludwikow & Sadowski 2008; Yendrek, Koester & Ainsworth 2015). However, there is limited knowledge on gene expression profiles associated with natural genetic variations in O3 tolerance (Frei, Tanaka, Chen & Wissuwa 2010; Choet al. 2013b). Little is known about how the expression levels of such genes are organized and how they function in coordination to regulate O3 tolerance. In recent years, gene co-expression network analysis has become available to understand the key factors involved in the gene-expression regulation (Langfelder & Horvath 2008). Constructing gene co-expression networks with transcriptomic data obtained from O3 tolerant cultivars would reveal the genetic basis of mechanisms of O3tolerance and the associated natural variations in rice.
The present study aimed to elucidate the genetic and physiological mechanisms underlying high O3 tolerance in the rice cultivar, Koshihikari. We focused on two Japanese rice cultivars, Koshihikari and Takanari, which have been reported to have contrasting O3 tolerances (Sawada & Kohno 2009; Sawada, Komatsu, Nanjo, Khan & Kohno 2012; Cho et al. 2013b). We conducted chamber experiments over two consecutive years to examine the effects of chronic O3 conditions on photosynthesis, growth, and yield-related traits in the two cultivars. In addition, we combined RNA-seq and gene co-expression network analyses to elucidate gene expression profiles and their regulatory networks associated with the contrasting O3 tolerance between the two rice cultivars.
Materials and Methods
Plant materials and cultivation
The experiments were performed in 2019 and 2020 using two rice cultivars: a leading Japanese rice cultivar of the subspeciesjaponica , Koshihikari (O3 tolerant), and a high-yielding rice cultivar of the subspecies indica , Takanari (O3 sensitive). Seeds of the two cultivars were sown in a seedling box on May 9 in both 2019 and 2020 and grown in a sunlit growth chamber (Koito Electric Industries, Shizuoka, Japan) at the National Institute for Environmental Studies, Tsukuba-city, Ibaraki, Japan (36°02’56.2 N, 140°07’02.8 E). In the growth chamber, air temperature and relative humidity were maintained at 25℃ and 70%, respectively. The seedlings were transplanted 11 d after sowing to 1/5000 Wagner pots, with one seedling per pot. Pots were filled with a mixture of commercial nursery soil (Kumiai Ibaraki Baido [0.5 g N, 0.9 g P, 0.5 g K kg-1], Ibaraki Kumiai Baido Co. Ltd., Ibaraki, Japan) and vermiculite (Fukushima no vermiculite, Fukushima Vermi Co. Ltd., Fukushima, Japan). Pots were watered and fertilized using 0.4 g of granular chemical fertilizer (Kumiai Kasei 7 Go [8% N, 8% P, and 5% K] Katakura Chikkarin Co. Ltd., Tokyo, Japan).
For chronic O3 exposure, the seedlings at four days after transplanting were transferred to a growth chamber, where [O3] was maintained at 80/0 nmol mol-1 for 7:00/19:00 and air temperature was controlled at 28℃/22℃ for 7:00/19:00, with a relative humidity of 70% under natural sunlight. The same was replicated in another chamber but under control conditions where only filtered air was circulated with 6 nmol mol-1 O3 throughout the day on an average. We cultivated 6 plants in 2019 and 16 plants in 2020 for Koshihikari and 6 plants in 2019 and 8 plants in 2020 for Takanari per treatment. Plant appearance was captured 53 days after the beginning of the treatment (DAT) in 2019.
Analysis of growth and yield-related traits
The aboveground parts of the plants were harvested after reaching full maturity in 2019 and 2020. Unhulled grain weight and total aboveground dry weight were measured after drying at 70°C for 72 h. Harvest index was calculated as the ratio of unhulled grain weight to total aboveground dry weight. Using unhulled grains, the 1000-grain weight and total number of grains were determined, and the number of grains per panicle was calculated. In 2020, plant height and tiller number were measured at 66–67 DAT and 89–90 DAT, respectively. The aboveground parts of the plants were harvested at 90 DAT to measure the total dry weight, leaf blade area and weight, and leaf-sheath weight of the whole plant. The number of replications for each cultivar and treatment was 4–11.
Gas exchange and chlorophyll fluorescence measurements
Concurrent measurements of gas exchange and chlorophyll fluorescence were conducted using a portable gas exchange system, LI-6800 (LI-COR , Lincoln, NE, USA). Measurements were made with rice plants at four growth stages in 2019, 49–50, 68–69, 88–89, and 108–109 DAT, and at two growth stages in 2020, 66–67 and 89–90 DAT. The CO2 assimilation rate (A ), stomatal conductance (g s), and CO2concentration [CO2] at intercellular airspaces (C i) per unit leaf area were measured in the uppermost fully expanded leaf from 8:30 to 14:00 under a [CO2] of 400 µmol m–1, photosynthetic photon flux densivty (PPFD) of 1500 µmol m–2 s–1, and air temperature of 30°C. In 2020, we estimated A for a single leaf blade (A leaf) by multiplying the leaf blade area andA per unit leaf area. The quantum yield of PSⅡ (Φ PSⅡ) was estimated using the following equation, as described by Baker (2008):
Φ PSⅡ = (F m’ –F ’)/F m
The electron transport rate for photosystem II (ETRⅡ) was estimated using the following equation:
ETRⅡ = PPFD × 0.5 × 0.875 × Φ PSⅡ
In addition, A was measured under a series of [CO2] at 100, 200, 300, 400, 500, 600, 750, 900, 1200, and 1500 µmol mol–1 forA -C i analysis at 54–56 DAT in 2019. The maximum rate of ribulose-1,5-bisphosphate (RuBP) carboxylation (V cmax) was estimated according to the biochemical model of C3 photosynthesis (Farquhar, Caemmerer and Berry, 1980) as described by Sakoda et al. (2016). The number of replications for each cultivar and treatment was 4–6.
Analysis of the physiological and morphological traits related to leaf photosynthesis
The Soil Plant Analysis Development (SPAD) value was measured on the top 1st, 2nd, 3rd, and 4th leaves of 4–5 plants at 49–50 DAT in 2019 using an SPAD meter (SPAD-502, Konica Minolta , Tokyo, Japan). In addition, the contents of nitrogen, ribulose-1,5-bisphosphate carboxylase/oxygenase (RuBisCO), and chlorophyll a +b per unit leaf area were measured on the same leaf where gas exchange was measured in 2019. For nitrogen quantification, the sampled leaf tissue (> 2 cm2) was dried at 70°C for 72 h, weighed, and ground into a powder. The nitrogen content per unit leaf area was determined as described by Kimura, Hashimoto-Sugimoto, Iba, Terashima & Yamori (2020). For RuBisCO and chlorophyll quantification, the sampled leaf tissue (> 2 cm2) was homogenized using a pre-cold mortar and pestle in an extraction buffer described by Sakoda et al. (2016). RuBisCO and chlorophylla +b contents were quantified as described by Sakoda et al. (2016).
The RuBisCO activation state was measured on the same leaf used for gas exchange measurements at 68–69, 88–89, and 108–109 DAT in 2019. After gas exchange measurements, the leaf section clamped in the chamber (< 2 cm2) was immediately sampled and frozen in liquid nitrogen. The rate of the carboxylation reaction catalyzed by RuBisCO before and after the activation of RuBisCO by MgCl2 and NaHCO3 was measured by coupling 3-phosphoglyceric acid formation with NADH oxidation, as described by Qu et al. (2021).
Stomatal density, defined as the stomatal number per unit leaf area, was measured on the same leaf as the gas exchange measurements and was conducted in 2019, as described by Sakoda et al. (2019). A replica of the abaxial side of the leaf blade was prepared using the nail polish impression method (Ceulemans et al., 1995). The replica was observed at 200× magnification using an optical microscope, and a digital image was obtained (VHX-6000, Keyence, Osaka, Japan). The stomatal number and guard cell length was manually counted using the image analysis software ImageJ (Schneider, Rasband & Eliceiri 2012).
RNA-seq analysis
Total RNA was extracted from leaf tissues at 66–67 and 89–90 DAT in 2020 using the RNeasy Plant Mini Kit (QIAGEN, Hilden, Germany). The poly-A mRNA was purified, and then, the library was constructed using the Illumina Truseq Stranded mRNA Sample Preparation Kit (Illumina, Inc., San Diego, USA). A total of 37 libraries were multiplexed and sequenced using the Illumina HiSeq X platform (Table S5). These sequence data have been submitted to the Sequence Read Archive (SRA) of the DNA Data Bank of Japan (DDBJ) under accession number DRA014708. Raw sequence data were processed for trimming and mapping onto the reference sequence, as described by Hashida et al. (2022), with minor modifications. The count data of each transcript were normalized by the iDEGES normalization method of a Shiny-based application, TCC-GUI (Su, Sun, Shimizu & Kadota 2019), with default settings.
We performed hierarchical clustering and principal component analysis (PCA) with the normalized expression level of the top-1000 highly expressed genes. The differentially expressed genes (DEGs) under chronic O3 conditions were detected at 66–67 and 89–90 DAT for each cultivar (Fig. S1) with |log2(fold-change)| >1 and false discovery rate = 0.05, using the TCC-GUI. To visualize the intersection of up- or downregulated DEGs among four groups of two cultivars in two stages, we used the package UpSetR (Conway, Lex & Gehlenborg 2017) in R software version 3. 6. 1 (R Foundation for Statistical Computing, Vienna, Austria). In addition, gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses were performed on the web server application g:Profiler (Raudvereet al. 2019), with a false discovery rate of 0.05.
Weighted gene co-expression network analysis (WGCNA)
To construct the gene co-expression network, we compared the transcriptome data merging 66–67 and 89–90 DAT between the control and chronic O3 conditions in Koshihikari and Takanari, respectively (Fig. S1). The co-expression network of DEGs was constructed using the R package WGCNA (Langfelder & Horvath 2008), with the following parameters: softthreshold = 21, minModuleSize = 20, mergeCutHeight = 0.1. WGCNA assigned the DEGs to 26 modules in Koshihikari and 33 modules in Takanari. Thereafter, Pearson’s correlation co-efficient (r ) was calculated between the module eigengene (ME), the first component in PCA of the gene expression profiles for each module, and the physiological traits that were investigated concurrently with the sampling for RNA extraction.
The constructed networks were visualized using the open-source software platform, Cytoscape ver. 3.8.2 (Shannon et al. 2003), with the edge weight > 0.1 in Koshihikari and 0.2 in Takanari. To explore the hub genes that play a central role in the gene co-expression network regulation in each module, we used NetworkAnalyzer and MCDS in Cytoscape. The hub genes in each module were detected as the genes that showed the top-5% highest degree and betweenness centrality calculated by NetworkAnalyzer and were assigned to the dominator or connector by MCDS.
Statistical analysis
The significance of differences in the SPAD values of the leaves at each position was compared between the plants under control and chronic O3 conditions for each cultivar using one-way ANOVA atp <0.05. The significance of differences in all physiological traits was also compared between the plants under control and chronic O3 conditions of each cultivar using one-way ANOVA at p <0.05. Pot was regarded as the unit for ANOVA. Statistical analyses were conducted using R. The correlations among all physiological traits were analyzed using a correlation analysis tool in the Scikit-learn and Seaborn libraries of Python.
Results
Effects of chronic O3 conditions on yield-related traits
We examined the effects of chronic O3 on yield-related traits in Koshihikari and Takanari (Fig. 1). Chronic O3conditions significantly decreased the total dry weight, unhulled grain weight, harvest index, 1000-grains weight, and total grain number in Takanari during both experimental years, whereas its effect on panicle number was inconsistent between the two years. In contrast, chronic O3 conditions significantly increased the total dry weight of Koshihikari in both years, whereas an increase in grain yield, harvest index, and total grain number was observed in 2020 only. The effects of chronic O3 on harvest index, 1000-grains weight, and panicle number were inconsistent between the two years in Koshihikari.
Effects of chronic O3 conditions on the photosynthesis- and growth-related traits
The effects of chronic O3 exposure on gas exchange and chlorophyll fluorescence parameters were evaluated at the four growth stages in 2019 (Fig. 2). Chronic O3 conditions significantly decreased A in both cultivars at all stages, except for the 1st measurement in Koshihikari. In addition, chronic O3 conditions decreasedg s and ETRⅡ in both cultivars throughout all the stages, although the decrease was not significant at some stages. The decrease in photosynthetic parameters by chronic O3conditions was greater in Takanari than in Koshihikari at later growth stages. There was no significant effect of chronic O3conditions on C i in Takanari, whereasC i was significantly higher under chronic O3 conditions than under the control conditions in Koshihikari. In addition, A -C i analysis showed lower A under chronic O3 conditions than under control conditions at any C i value in Koshihikari. A was lower under chronic O3conditions at C i below 600 µmol mol-1, whereas it was similar between the control and chronic O3 conditions at C i above 700 µmol mol-1 in Takanari. In both cultivars, chronic O3 conditions significantly decreasedV cmax.
The appearance of the plants under control and chronic O3 conditions in 2019 is shown in Fig. 2G and H. The chronic O3 conditions seemed to have a positive effect on growth and induced minor visible injury in Koshihikari. In contrast, this condition induced severe growth inhibition and visible injury in Takanari plants. The magnitude of visible injury was evaluated at several leaf positions based on the SPAD values (Fig. 2I and J). The SPAD values at the 2nd, 3rd, and 4th leaves were significantly lower under chronic O3 conditions than under control conditions in both cultivars. The decrease in SPAD value was similar among the leaves at different positions in Koshihikari, whereas in Takanari, it was larger in the older leaves than that in the younger leaves.
In 2019, the nitrogen content significantly decreased under chronic O3 conditions at all stages in Koshihikari and at 88–89 and 108–109 DAT in Takanari (Fig. 3). In addition, chronic O3 conditions significantly decreased the chlorophyll and RuBisCO contents at most of the measurement stages but had no effect on the RuBisCO activation state in either cultivar. Chronic O3 conditions had no clear or consistent effect on stomatal density and guard cell length throughout all measurement stages.
In 2020, we evaluated A , a single leaf blade area, andA leaf at 66–67 and 89–90 DAT under control and chronic O3 conditions (Fig. 4). Chronic O3 conditions significantly decreased A at 89–90 DAT in Koshihikari and at both stages in Takanari. The single leaf-blade area was significantly higher in Koshihikari under chronic O3 conditions than under control conditions, whereas it was lower at both stages in Takanari under chronic O3conditions. Chronic O3 conditions had no significant effect on A leaf in Koshihikari at either stage, whereas they significantly decreased A leaf in Takanari at both stages. At 89–90 DAT, plant height, tiller number, total dry weight, and whole leaf blade area and weight were significantly lower in Takanari plants under chronic O3conditions than under control conditions. Plant height, whole leaf blade area and weight, and the ratio of whole leaf blade weight to leaf sheath weight (blade/sheath) were significantly higher in Koshihikari under chronic O3 conditions than that under control conditions.
Correlations among photosynthesis, growth, and yield-related traits
We investigated correlations among photosynthesis, growth, and yield-related traits in Koshihikari and Takanari under control and chronic O3 conditions (Fig. 5). For most trait combinations, Takanari showed higher r values than Koshihikari, suggesting a greater sensitivity to chronic O3conditions in Takanari. Grain weight showed a significant correlation with the single leaf-blade area, plant height, and C i in Koshihikari, whereas it showed a significant correlation with all traits except for A and ETRⅡ in Takanari. In addition, grain yield showed the highest r value with a single leaf blade area (r = 0.790) in Koshihikari and tiller number (r = 0.909) in Takanari. A significant correlation was observed between grain yield and single leaf-blade area (r = 0.590) in Takanari, whereas in Koshihikari, no clear relationship between grain yield and tiller number was noted.
Transcriptional profile and gene co-expression network under chronic O3 conditions
To elucidate the gene co-expression network associated with O3 response in Koshihikari and Takanari, we conducted RNA-seq and network analyses in 2020. Hierarchical clustering analysis and PCA showed the major effects of genotype and chronic O3 conditions on the variations in the transcriptional profiles (Figs. S2 and S3). In Koshihikari, chronic O3conditions upregulated 963 genes and downregulated 1542 genes at 66–67 DAT and upregulated 2653 genes and downregulated 2680 genes at 89–90 DAT (Fig. 6). In Takanari, chronic O3 conditions upregulated 5231 genes and downregulated 3948 genes at 66–67 DAT and upregulated 5548 genes and downregulated 4258 genes at 89–90 DAT, indicating that chronic O3 conditions had a larger impact on the transcriptional profiles in Takanari than in Koshihikari. Under chronic O3 conditions, 431 upregulated and 321 downregulated DEGs were detected in both cultivars and growth stages, which is important for the genotype-unspecific O3response of rice (Fig. 6, Table S6). In contrast, 126 upregulated and 525 downregulated DEGs were detected in Koshihikari but not in Takanari at both growth stages, indicating the possible association of these DEGs with a Koshihikari-specific O3 response.
We listed the up- and downregulated DEGs with the top-25 lowest adjustedp -values at each stage in Koshihikari (Table 1) and Takanari (Table 2). Among the 50 genes, seven upregulated and four downregulated DEGs were commonly detected at both stages in Koshihikari, whereas six upregulated and eight downregulated DEGs were commonly detected in Takanari. In addition, the upregulated gene, Os03g0103300 (hybrid proline- or glycine-rich protein, control of low-temperature germinability, and pre-harvest sprouting resistance), showed the lowestp values at both stages in Koshihikari.
The enriched GO terms and KEGG pathways under chronic O3conditions were investigated in Koshihikari and Takanari at the two growth stages (Table S7). In the present study, we focused on the enriched GO terms for biological processes and showed those with the top-10 lowest p -value for each cultivar and stage (Fig. 6). In Koshihikari, the top-10 enriched GO terms with upregulation were inconsistent between the two stages. Among the top-10 enriched GO terms with downregulation, “Energy reserve metabolic process” (GO:0006112) was consistent between the two stages. In Takanari, eight out of the top-10 upregulated GO terms, including transport and localization, and eight out of the top-10 downregulated GO terms, including photosynthesis, were commonly enriched at two stages. The GO terms associated with photosynthesis (GO:0015979, GO:0019684, GO:0009765, and GO:0009768) were highly enriched with upregulation in Koshihikari at 66 DAT and with downregulation in Takanari at both stages. In addition, the GO terms associated with carbohydrate metabolic process (GO:0005975 and GO:0044262) were highly enriched with downregulation in Koshihikari at both stages, whereas they were not affected by chronic O3 conditions in Takanari.
On comparing the transcriptome data merging 66–67 and 89–90 DAT between the control and chronic O3 conditions in Koshihikari and Takanari, respectively, we detected 5081 genes in Koshihikari and 13352 genes in Takanari as DEGs induced by chronic O3 conditions (Fig. S1). We performed WGCNA with these DEGs in Koshihikari and Takanari, respectively. The DEGs were assigned to 26 modules, including 27–663 genes, in Koshihikari (Fig. 7), and 33 modules, including 38–2827 genes, in Takanari (Fig. 8). We then analyzed the r value between the ME of each module and the single leaf-blade area in Koshihikari or tiller number in Takanari because these traits showed the highest r value with unhulled-grain weight (Fig. 5). The MEs of 17 modules showed a significant correlation with the single leaf-blade area in Koshihikari and those of 33 modules showed a significant correlation with the tiller number in Takanari. Topological analysis of these modules detected 34 hub gene candidates in Koshihikari and 83 hub gene candidates in Takanari (Table S9). The hub gene candidates with the highest degree in each module were listed for Koshihikari (Table 3) and Takanari (Table 4).
Among the constructed modules, the blue, brown, turquoise, and yellow modules play an important role in the gene regulatory mechanisms underlying O3 tolerance because these were the largest of all modules and those MEs showed a significant correlation between single leaf-blade area in Koshihikari and tiller number in Takanari. In Koshihikari, the following genes showed the highest degree and betweenness centrality in each of the four modules: Os02g0127900(“Conserved hypothetical protein”) in brown module,Os09g0572900 (Dynamin-related protein 1E, Negative regulation of programmed cell death, Control of cytochrome c release) in blue module,Os04g0482300 (“Kelch related domain containing protein”) in turquoise module, and Os12g0576600 (“Metallophosphoesterase domain containing protein”) in yellow module (Fig. 7B-E, Table 3). In addition, three hub gene candidates, Os06g0160700 (“Starch synthase”) in the turquoise module and Os07g0106000(“Metallophosphoesterase domain containing protein”) andOs11g0658900 (“Similar to Lipase family protein”) in the brown module, were included in the up- or down-regulated DEGs with the top-100 lowest adjusted p values. In Takanari, the turquoise module had the largest size and the highest r value (0.940) between its ME and tiller number of all modules. The following genes showed the highest degree in each of the blue, brown, turquoise, and yellow modules:Os03g0803000 (synaptobrevin domain containing protein) in the brown module, Os02g0819600 (serine/threonine protein kinase-related domain containing protein) in the blue module,Os02g0156600 (hypothetical conserved gene) in the turquoise module, and Os07g0685500 (alpha/beta hydrolase family protein) in the yellow module. The three hub gene candidates, Os02g0156600and Os04g0623800 (similar to aminomethyltransferase, mitochondrial precursor (EC 2.1.2.10) (glycine cleavage system T protein; GCVT) in the turquoise module and Os12g0182300 (protein kinase, core domain containing protein) in the blue module, were included in the up- or down-regulated DEGs with the top-100 lowestp values.
The enriched GO terms and KEGG pathways in each module were also investigated for both the cultivars (Table S10). Focusing on the GO terms associated with carbohydrate metabolic process and photosynthesis, which were differently enriched between the two cultivars, those terms were highly enriched in the turquoise and green modules in Koshihikari. In the turquoise module, two genes functioning in carbohydrate metabolic processes, Os01g0897600 (glycoside hydrolase, family 1 protein) and Os07g0662900 (similar to 4-alpha-glucanotransferase (EC 2.4.1.25)), were identified as hub genes (Fig. 7, Table S9). In the green module of Koshihikari, three genes, Os12g0291100 (similar to Petunia ribulose 1, 5-bisphosphate carboxylase small subunit mRNA (clone pSSU 51), partial cds. (fragment)].), Os01g0662700(similar to naphthoate synthase (EC 4.1.3.36)), and Os03g0563300(magnesium-chelatase subunit ChlI, Mg-chelatase I subunit, chlorophyll biosynthesis, chloroplast development), were identified as hub gene candidates (Fig. 9). The expression levels of these three genes were significantly upregulated at 66–67 DAT in Koshihikari but downregulated at both stages in Takanari under chronic O3 conditions.
Discussion
Physiological mechanisms underlying the high O3 tolerance in Koshihikari
Koshihikari and Takanari showed contrasting responses to chronic O3 conditions in the leaf blade area at the single-leaf and whole-plant levels (Fig. 4). The greater leaf area due to increased biomass allocation to leaf blades might enhance light interception efficiency at the whole-plant level, contributing to greater CO2 assimilation and biomass production in Koshihikari under chronic O3 conditions than under control conditions (Fig. 1–5). In contrast, the decreased single leaf blade area and tiller number resulted in a severe reduction in the whole leaf blade area and biomass production in Takanari under chronic O3 conditions (Fig. 1 and 4). It is commonly known that chronic O3 conditions affect biomass allocation in various plant species (Cooley & William 1987). Elevated O3 conditions resulted in greater biomass production in Japanese oak species owing to increased biomass allocation to leaves (Kitao et al., 2015). Kobayashi et al. (1995) also reported higher biomass allocation to leaf blades in Koshihikari but lower grain yields under elevated O3 conditions. In Koshihikari, the O3 effects on growth and yield varied under different nitrogen fertilization regimes (Tatsumi, Abiko, Kinose, Inagaki & Izuta 2019). These results suggest that increased biomass allocation to leaf blades can be an important factor for the high O3tolerance in Koshihikari, depending on environmental conditions such as soil nutrient levels.
While chronic O3 conditions significantly reducedA in Koshihikari and Takanari, a greater reduction in Awas observed in Takanari at the later growth stage (Fig. 2), which might contribute to a more severe reduction in biomass. O3-induced A reduction is attributable to the decreased activity of gas diffusion from the atmosphere to the chloroplast stroma and biochemical reactions related to CO2 fixation in chloroplasts (Ainsworth et al.2012). In the present study, chronic O3 conditions decreased g s in both cultivars; this decrease was much larger in Takanari than in Koshihikari. g scan be regulated by the size, depth, and opening of a single stoma and its density (Franks & Beerling 2009). There was no consistent reduction in stomatal density and size in either cultivar during the growth period, indicating that the g s reduction was due to O3-induced stomatal closure (Frei 2015). Chronic O3 conditions also significantly decreased the biochemical activity evaluated by ETRⅡ and V cmax, which was accompanied by reduced nitrogen and RuBisCO contents in the leaves of both cultivars (Fig. 2 and 3). In contrast, the single leaf-blade area substantially increased, andA leaf was not affected in Koshihikari under chronic O3 conditions (Fig. 4). These results suggest that increased leaf blade area with a decrease in leaf thickness resulted in lower N and RuBisCO contents per unit leaf area, which would cause an apparent reduction in A in Koshihikari. Therefore, it could be concluded that photosynthetic performance was maintained in Koshihikari under chronic O3 conditions, contributing to the superior O3 tolerance.
There was no clear difference in g s between Koshihikari and Takanari under chronic O3 conditions (Fig. 3 and 4), whereas visible injury was more severe in Takanari, especially in older leaves (Fig. 2). This suggests superior ROS detoxification ability after generation of ROS in the leaves of Koshihikari compared to those of Takanari. Iseki, Homma, Endo & Shiraiwa (2013) reported that the japonica rice cultivars such as Koshihikari show greater tolerance to oxidative stress than theindica rice cultivars such as Takanari. Moreover, Sawada and Kohno (2009) indicated that in rice, the japonica cultivars are superior to indica cultivars in terms of O3tolerance. Possibly, the high ROS detoxification ability may enable to maintain stomatal opening and biochemical activity for photosynthesis and biomass production under chronic O3 conditions in Koshihikari.
Genetic mechanisms underlying variations in O3 tolerance between Koshihikari and Takanari
Many studies have been conducted to elucidate O3tolerance mechanisms using mutant populations of model plants, deciphering the involvement of various genes (Kanna et al. 2003; Tamaoki et al. 2003; Saji et al. 2008, 2017; Nagatoshiet al. 2016). In addition, previous studies using transcriptome analysis have identified O3-responsive genes involved in transcription, signal transduction, transport, metabolism, redox control, senescence, and defense response (Ludwikow & Sadowski 2008). We confirmed that the genes associated with these biological processes were significantly affected by chronic O3 conditions in Koshihikari and Takanari (Fig. 6, Tables S6 and S7). Using microarray analysis, Frei et al. (2010) identified the gene encoding a putative ascorbate oxidase, which is associated with natural genetic variations in rice O3 tolerance. In the present study, RNA-seq analysis detected 124 upregulated and 525 downregulated DEGs in Koshihikari but not in Takanari under chronic O3 conditions, which possibly contributed to the high O3 tolerance in Koshihikari (Fig. 6). In addition, gene co-expression network analysis revealed the gene modules whose expression profiles corresponded to O3 response in each cultivar (Fig. 7 and 8). The hub gene candidates detected in the modules could have key functions in the chronic O3 responses of Koshihikari and Takanari (Tables 3, 4, S9). The hub gene candidates did not overlap between the two cultivars, suggesting that the gene regulatory mechanisms underlying O3 response can differ between rice cultivars.
In Koshihikari, Os03g0815800 (OsERF58 ; similar to ethylene-responsive transcription factor 5 (ethylene-responsive element binding factor 5; EREBP-5 ) was detected as the hub gene candidate in the brown module, whose ME showed the highest r value with a single leaf-blade area (Fig. 7B, Table S9). Using genome-wide association studies, Ueda et al. (2015) reported that polymorphisms in EREBP were associated with natural genetic variations in O3 tolerance in rice. It is possible thatOsERF58 plays an important role in the gene regulatory mechanisms underlying the O3 response of the leaf blade area in Koshihikari.
In Koshihikari, two hub gene candidates, Os07g0106000 in the brown module and Os12g0576600 in the yellow module, encoded metallophosphoesterase (MPE) domain-containing proteins (Tables 3 and S9). The MPE domain is included in diverse enzymes such as mammalian phosphoprotein phosphatases, acid sphingomyelinases, and purple acid phosphatases (Matange, Podobnik & Visweswariah 2015). In Arabidopsis, an activation-tagging mutant of purple acid phosphatase increased foliar ascorbate content and improved acute O3 tolerance (Zhang, Gruszewski, Chevone & Nessler 2008). In contrast,Os12g0576600 , annotated as purple acid phosphatase, andOs07g0106000 were significantly downregulated by chronic O3 conditions in Koshihikari. It is suggested that the change in the expression level of these genes might be associated with the O3 response of the leaf blade area in Koshihikari; however, these genes did not directly contribute to ROS detoxification.
The GO terms associated with carbohydrate metabolism were highly enriched and downregulated under chronic O3 conditions in Koshihikari at both growth stages but not in Takanari (Fig. 6, S7). The relationship between carbohydrate metabolism and grain quality under chronic O3 conditions has been previously discussed in rice (Frei 2015; Sawada et al. 2016). However, carbohydrate metabolism has not been recognized as an important factor underlying the O3 response in terms of growth, yield, and the associated genetic variation in crops. Genes associated with carbohydrate metabolism were highly enriched in the turquoise module of Koshihikari; three genes, Os01g0897600 , Os06g0160700 , andOs07g0662900 , were identified as hub genes (Table S9). Importantly, Os06g0160700 , soluble starch synthase I, was included among the DEGs with the top-100 lowest adjusted p -values at both growth stages in Koshihikari. It is suggested that the co-expression network of genes involved in carbohydrate metabolism is associated with superior O3 tolerance in Koshihikari compared to that in Takanari.
Typically, the genes associated with photosynthesis are downregulated under chronic O3 conditions (Ludwikow & Sadowski 2008), as observed in Takanari (Fig. 6). In contrast, these genes were significantly upregulated or unaffected in Koshihikari under chronic O3 conditions (Fig. 6), consistent with the lack of decrease in A leaf. This result is supported by Sawada et al. (2012), who reported the upregulation of photosynthesis-related proteins in Koshihikari under chronic O3 conditions. The photosynthesis-related genes were highly enriched in the green module in Koshihikari and the turquoise-module in Takanari, whose ME was significantly correlated withA leaf (r = 0.58, Koshihikari and 0.85, Takanari). Under chronic O3 conditions, the hub gene candidates Os12g0291100 , Os01g0662700 , andOs03g0563300 in the green module of Koshihikari were significantly upregulated or unaffected in Koshihikari, whereas their expression levels were downregulated at both stages in Takanari. This contrasting response is associated with the genotypic variation in the O3 response in terms of photosynthesis.
Among all the analyzed genes, the expression level ofOs12g0291100 (OsRbcS3 ), which encodes the RuBisCO small subunit (RbcS), was relatively high in both cultivars. The five genes (OsRbcS1-5 ) were identified to encode the RbcS in rice (Suzuki, Miyamoto, Yoshizawa, Mae & Makino 2009). OsRbcS3 is a representative RbcS gene in rice (Suzuki et al. 2009). In addition, the expression level of OsRbcS3 in photosynthetic tissues is typically high and has a significant effect on RuBisCO content and photosynthetic performance in rice (Morita, Hatanaka, Misoo & Fukayama 2014; Kanno, Suzuki & Makino 2017). Downregulation ofOsRbcS s has been reported in rice exposed to acute O3 conditions (Cho et al. 2008). In the present study, the expression level of OsRbcS3 was nearly 3-fold higher under chronic O3 conditions than under control conditions in Koshihikari at 66–67 DAT, whereas it decreased by up to 57.7% in Takanari (Fig. 9). It suggests that OsRbcS3 plays a central role in the regulatory network of photosynthesis-related genes under chronic O3 conditions, contributing to superior O3 tolerance in Koshihikari compared to that in Takanari.
In conclusion, Koshihikari showed no reduction in growth- and yield-related traits under chronic O3 conditions, whereas Takanari showed a substantial reduction. The greater O3 tolerance in Koshihikari than in Takanari could be attributed to (1) increased leaf blade area due to increased biomass allocation to leaf blades and (2) higher stability of photosynthetic performance in Koshihikari. RNA-seq and gene co-expression network analyses revealed significant differences in the gene expression profiles and regulatory networks between Koshihikari and Takanari. The different expression profiles of genes involved in carbohydrate metabolism and photosynthesis could be associated with the genetic variation in O3 tolerance between the two rice cultivars. OsRbcS3 , which encodes the RuBisCO small subunit, showed contrasting expression patterns between the two cultivars under chronic O3 conditions. In addition, this gene was identified as a hub gene candidate in the gene co-expression network, which was highly correlated with photosynthetic performance, suggesting that OsRbcS3 plays a key role in the genetic mechanisms underlying the O3-response in Koshihikari. The genes which were suggested to contribute to the high O3tolerance of Koshihikari in the present study may be a promising target in rice breeding towards improved O3 tolerance.