Integrated transcriptomic and metabolomic analysis provides insight into the regulation of leaf senescence in rice

Leaf senescence is one of the most precisely modulated developmental process and affects various agronomic traits of rice. Anti-senescence rice varieties are important for breeding application. However, little is known about the mechanisms underlying the metabolic regulatory process of leaf senescence in rice. In this study, we performed transcriptomic and metabolomic analyses of the flag leaves in Yuenong Simiao (YN) and YB, two indica rice cultivars that differ in terms of their leaf senescence. We found 8524 genes/204 metabolites were differentially expressed/accumulated in YN at 30 days after flowering (DAF) compared to 0 DAF, and 8799 genes/205 metabolites were differentially expressed in YB at 30 DAF compared to 0 DAF. Integrative analyses showed that a set of genes and metabolites involved in flavonoid pathway were significantly enriched. We identified that relative accumulation of PHENYLALANINE AMMONIA-LYASE (PAL), CINNAMATE 4-HYDROXYLASE (C4H), 4-COUMAROYL-COA LIGASE (4CL), CHALCONE SYNTHASE (CHS) and CHALCONE ISOMERASE (CHI) in YN30/0 was higher than that in YB30/0. Three flavonoid derivatives, including phloretin, luteolin and eriodictyol, showed lower abundances in YB than in YN at 30 DAF. We further revealed a MYB transcription factor, which is encoded by OsR498G0101613100 gene, could suppress the expression of CHI and CHS. Our results suggested a comprehensive analysis of leaf senescence in a view of transcriptome and metabolome and would contribute to exploring the molecular mechanism of leaf senescence in rice.

Leaf senescence, an orderly process of degradation and transformation, is the last stage of leaf development and is controlled by complex internal factors such as hormones, genetics, and external gene network regulatory factors, including light, water, temperature, mineral elements, and microorganisms [1][2][3] . Leaf senescence plays important roles in regulating nutrient distribution and improving the ability of plants to adapt to the environment. Through the process of senescence, leaves transfer nutrients to seeds or fruits at the end of development, providing necessary energy for offspring, seed germination and early growth of seedlings 4 . During senescence, leaf cells undergo dramatic changes in cell metabolism, structure, and gene expression. The most striking feature of these changes is that the decomposition of chlorophyll during chloroplast denaturation results in the yellowing of leaves, followed by hydrolysis of macromolecules such as lipids, proteins and nucleic acids, which in turn leads to separation of mitochondria and nuclei and ultimately cell death 5 . Together, the identification and characterization of hundreds of senescence-associated genes (SAGs) in plants and senescence-related mutants constitute one of the best approaches for understanding leaf senescence at the molecular level [5][6][7] . The leaf senescence database (LSD) also provides comprehensive information concerning SAGs and their corresponding mutants 8 .
Transcriptome sequencing is an important method for obtaining gene expression data in organisms 9 . Highthroughput profiling of transcripts is also an efficient method for investigating the process of plant leaf senescence, as has been demonstrated for Arabidopsis thaliana, rice, tobacco, soybean, wheat, cotton, and maize [10][11][12][13][14][15] . However, metabolites are the basis and direct manifestation of an organism's phenotype. Thus, metabolomics provides a global physiochemical view of cellular status during leaf senescence. By detecting the major carbohydrate and nitrogen metabolite markers of tobacco leaves at different developmental stages, researchers identified www.nature.com/scientificreports/ a sink-to-source transition at a particular stage in which an accumulation of carbohydrates and a depletion of both organic and inorganic nitrogen stores were observed 16 . By determining changes in chemical composition during leaf senescence in Arabidopsis thaliana, researchers also found that the C, Cr, Cu, Fe, K, Mo, N, P, S and Zn mobilized from senescing leaves 17 . In recent years, studies combining transcriptomics and metabolomics have been successfully used to research the biosynthesis pathway of medicinally important active ingredients in plants 18,19 , to study the fruit ripening and disease resistance of tomato 20,21 , to elucidate molecular mechanism underlying rice drought tolerance 22 , and so on. New genes related to senescence could still be discovered, and the metabolic pathways in which they are involved concerning the regulation of senescence should be explored by multi-omic approaches. Rice is one of the most important crop species in the world and is the staple food for more than half of the global population 23 . Rice yield and other agronomic traits are of concern to scientists. In the process of breeding, whether the process of functional leaves turning from green to yellow is stable and smooth in the late-mature stage of rice is also a concern of breeders. This stage is accompanied by the degradation of chlorophyll, pigment deposition and rapid changes in metabolite content. Leaves, especially flag leaves, which are the main functional organs of photosynthesis in rice, are the primary contributor to the accumulation of dry matter in grains. Therefore, premature leaf senescence is one of the main factors affecting rice yield stability. The leaf senescence database currently contains more than 180 SAGs experimentally identified in rice 8 ; these SAGs are involved in chloroplast degradation, hormones and transcription factors associated with photosynthesis, the energy metabolism pathway, nitrogen mobilization and so on. However, the biochemical basis and molecular regulatory mechanisms underlying these processes in rice are not clear.
Yuenong Simiao (YN) is a major rice cultivar in South China and exhibits a strong anti-senescence phenotype at the late-mature stage. To explore the underlying molecular mechanisms of the anti-senescence of YN, we performed transcriptomic and metabonomic analyses of the flag leaves of YN as well as YB, another rice variety that exhibits an early-senescence phenotype but is similar to YN in terms of flowering time, resistance, and yield. We identified 8524 differentially expressed genes (DEGs) in YN and 8799 DEGs in YB. By integrating transcriptome and metabolomic data, we found that a set of genes and metabolites involved in flavonoid pathway may exhibit different expression patterns among YN and YB. Interestingly, we found that one MYB TF, OsR498G0101613100, could inhibit the expression of CHALCONE ISOMERASE (CHI) and CHALCONE SYNTHASE (CHS), which were indispensable for the flavonoid pathway in rice. The findings can help to understand the regulatory mechanisms of flavonoid biosynthesis and may be useful for the further research on anti-senescence in rice.

Results
High-quality Yuenong Simiao rice exhibits a significant anti-senescence phenotype. Yuenong Simiao is a major rice variety with good quality in South China, with a total planting area of more than 200,000 ha. This variety has high yields, strong disease resistance and other excellent traits. In addition, the leaf colour changes slowly from green to yellow in the late stage of development, which is an important characteristic selected by breeders. To study the characteristics of leaf colour transformation, we chose YB as a control variety. The flowering period of YB was the same as that of YN, but the leaf colour changed from green to yellow and withered rapidly. Given that the leaf colour change from green to yellow is a notable indicator of senescence of plants, we performed detailed physiological, biochemical and molecular assays to determine whether YN has anti-senescence traits.
After flowering, YN plants exhibited a delayed green-to-yellow phenotype, becoming senescent approximately 10 days later than the YB plants did (Fig. 1a, b). In addition, compared with the YB plants, the YN plants underwent slower chlorophyll degradation. The most pronounced difference between the YN and YB plants was that, at 30 days after flowering (DAF), the chlorophyll content of YB decreased dramatically, while that of YN showed a slower decrease (Fig. 1c). Darkness is one of the most powerful external stimuli of leaf senescence and is used frequently as an effective method of simulating synchronous senescence 24 . Leaf segments of YN and YB were subjected to darkness for 5 days. We found that, compared with YB, YN showed significantly delayed dark-induced leaf senescence (Fig. 1d). Reactive oxygen species (ROS) constitute one of the factors that promote senescence 25 . We used 3,3′-diaminobenzidine (DAB) staining to examine H 2 O 2 levels in the two rice cultivars, and the results showed that the YN leaves displayed less H 2 O 2 accumulation than the YB leaves did (Fig. 1e). Similarly, chlorophyll degradation-related genes (Red Chlorophyll Catabolite Reductase 1 (RCCR1), NYC1 Like (NOL) and Non Yellow Colouring 3 (NYC3)) and three other senescence-associated genes (SAGs) (Oryza sativa NAC-like, activated by apetala 3/pistillata (OsNAP), Osh36 and OsI57) were expressed at higher levels in the fully expanded leaves of YB than in those of YN (according to the qPCR results) (Fig. 1f). These results clearly demonstrated that YN is an elite rice cultivar with a strong anti-senescence phenotype at the filling stage.
Transcriptional analysis reveals a set of differentially expressed genes involved in leaf senescence. To understand the transcriptional landscape in leaves during senescence, we performed an RNA-seq analysis of the flag leaves of YN and YB. We evaluated the flag leaves at 0 DAF and 30 DAF. Pearson's correlation coefficients between biological replicates indicated high reproducibility of the biological replicates and the array platform (Supplemental Table S1). Nearly 80% of the transcripts were covered more than 60% of the reads, indicating substantial coverage of the transcriptome (Supplemental Table S1). A snapshot of the transcriptome was generated by clustering all the samples using principal component analysis (PCA). Consistent with their distinct developmental programmes and cultivars, the samples of different cultivars at different developmental stages were clustered (Supplemental Figure S1). Compared to those in the 0 DAF samples, the transcripts in 30 Table S2). GO and KEGG analyses were performed for the DEGs, revealing information on the annotations of genes involved in biological processes, molecular functions and cellular components (GO analysis) and those involved in cellular processing, environmental information processing, genetic information processing, metabolism and organismal systems (KEGG analysis) ( Fig. 2a, b, Supplemental Figure S2, Supplemental Table S3). We analysed the top 15 metabolic pathways that were significantly enriched and found that genes differentially expressed in YN at 30 DAF compared to 0 DAF were significantly enriched in the following: photosynthesisantenna proteins, porphyrin and chlorophyll metabolism, arachidonic acid metabolism, glutathione metabolism, flavonoid biosynthesis, the mRNA surveillance pathway, benzoxazinoid biosynthesis, carbon metabolism, the pentose phosphate pathway, ribosome biogenesis in eukaryotes, isoflavonoid biosynthesis, aminoacyl-tRNA biosynthesis, fructose and mannose metabolism, fatty acid elongation, taurine and hypotaurine metabolism. Similarly, and found that genes differentially expressed in YB at 30 DAF compared to 0 DAF were significantly enriched in the following: the mRNA surveillance pathway, glutathione metabolism, photosynthesis-antenna proteins, isoflavonoid biosynthesis, fructose and mannose metabolism, arachidonic acid metabolism, the pentose phosphate pathway, flavonoid biosynthesis, carbon metabolism, porphyrin and chlorophyll metabolism, glycerophospholipid metabolism, MAPK signalling pathway-plant, alpha-linolenic acid metabolism, carotenoid biosynthesis, and glycerolipid metabolism (Fig. 2c, d, Supplemental Table S3). Interestingly, YN exhibited more upregulated genes at 30 DAF than at 0 DAF, and YB exhibited more downregulated genes. We identified 3016 genes whose expression was upregulated and 452 genes whose expression was downregulated only in YN0 versus YN30, and we identified 1217 genes whose expression was upregulated and 2526 genes whose expression was downregulated only in the YB0 versus YB30 comparison. A common set of 4851 (3103 upregulated and 1748 downregulated) genes was detected between the YN0 versus YN30 and YB0 versus YB30 comparisons (Fig. 2e). Taken together, these results suggest that at 30 DAF, when YN is still in a period of strong gene expression and metabolism, the gene expression in YB is attenuated (Fig. 2a, b), which is consistent with YB showing an increased-senescence phenotype. Transcription factors play important regulatory roles in plant growth and development. In total, 385 and 448 transcription factors (TFs) were differentially expressed in YN and YB, respectively ( Fig. 2f, g, Supplemental

Metabolome profiling analysis reveals a set of changes in metabolites involved leaf senescence.
To investigate the main metabolic changes reflecting the variation during leaf senescence and to evaluate their putative role in metabolic signalling, a widely targeted metabolome method was used to quantify primary and secondary metabolites identified at different leaf stages. Flag leaves of YN and YB sampled at 0 DAF and 30 DAF were subjected to UPLC-ESI-MS/MS analysis. In this work, 512 and 510 metabolites were identified and quantified in YB and YN, respectively (Supplemental Table S5). PCA showed that YB and YN were clearly separated in the PC1 × PC2 score plots (Supplemental Figure S3). Of the metabolites, 204 markedly changed in YN at 30 DAF compared to 0 DAF, in which the accumulation of 121 metabolites was upregulated and that of  Table S5). Cluster analysis was carried out for different metabolic components, including phenolic acids, amino acids and their derivatives, alkaloids, flavonoids, lipids, nucleotides and their derivatives, organic acids, lignans and coumarins, tannins, etc., which revealed different degrees of expression in the two rice cultivars at 30 DAF compared to 0 DAF (Fig. 3a). Among the content of these metabolites, the content of organic acids and amino acids decreased at the later stage of ageing, and flavonoids, alkaloids and lipids accumulated more in YN at 30 DAF than in YB at 30 DAF. We then classified the different KEGG pathways (Supplemental Table S6). We analysed the top 20 metabolic pathways with the largest number of differentially accumulated metabolites (Fig. 3b, c) and found that, in the two rice cultivars, metabolites involved in the following accumulated more at 30 DAF than at 0 DAF: metabolic pathways; biosynthesis of secondary metabolites; ABC transporters; aminoacyl-tRNA biosynthesis; arginine and proline metabolism; 2-oxocarboxylic acid metabolism; alanine, aspartate and glutamate metabolism; flavonoid biosynthesis; etc. In recent years, we have noticed that secondary metabolism, such as phenylpropanoid metabolism and flavonoid metabolism, has been reported to be involved in reducing oxidative damage 26 . Our results are consistent with those in the literature, further suggesting that the substances involved in these metabolic processes may be involved in the regulation of plant senescence.
Integrated analysis reveals substantial changes in the flavonoid biosynthesis pathway during leaf senescence. To explore the differences in regulatory mechanisms underlying leaf senescence between the two rice cultivars, differentially accumulated metabolites and DEGs in KEGG pathways were selected for integrated analysis. We selected the metabolic pathways of differentially accumulated metabolites and their corresponding genes (Table 1). It is worth noting that the flavonoid synthesis pathway in both rice cultivars was significantly enriched in differentially accumulated metabolites and DEGs. We rearranged 10 major compounds of flavonoid biosynthesis and 5 DEGs to corresponding positions of the flavonoid biosynthesis pathway (Fig. 4a,  b). In the flavonoid biosynthesis pathway, the abundance of 6 compounds increased at 30 DAF compared with 0 DAF between the two cultivars (Fig. 4b, d). Our metabolic analysis revealed that flavonoids and flavonoid derivatives, including phloretin, luteolin and eriodictyol, showed lower abundances in YB than in YN at 30 DAF (Fig. 4c, d). Secondary metabolites, including flavonoids, play an important role in maintaining redox homeostasis in cells by quenching reactive oxygen species 27,28 . The transcriptomic data showed that the relative expression of 5 DEGs (PAL, C4H, 4CL, CHS, CHI) in YN30/YN0 were generally higher than in YB30/YB0. To confirm this, we detected the expression of these five genes by qRT-PCR. The results were consistent with those obtained from RNA sequencing, the relative expression of C4H, 4CL, CHS and CHI were significantly increased, relative expression of PAL increased slightly, but not significantly, indicating that the sequencing data sets are reliable.  Figure S4). Taken together, these results suggest that the increased accumulation of flavonoids may be due to alterations in the expression of these DEGs in YN.
A MYB transcription factor negatively regulates the expression of OsCHS and OsCHI. Transcription factors play important roles in plant growth, development, reproduction and senescence. As shown above, a large number of MYB transcription factors were differentially expressed at the later filling stage between the two rice cultivars and thus may be involved in important biological functions. Previous studies have reported that MYB TFs are involved in regulating the expression of genes involved in the flavonoid synthesis pathway in plants [29][30][31] . According to our results, the expression of OsCHS and OsCHI, which encode chalcone synthase and chalcone isomerase, respectively, was significantly downregulated in YB at 30 DAF, which did not significantly differ from the results of YN (Supplemental Figure S4). We speculate that MYB TFs may be involved in the differential expression of CHS and CHI between the two rice varieties. To examine the relationship between MYBs  www.nature.com/scientificreports/ and flavonoid synthesis, we first analysed and screened potential MYBs that positively or negatively regulated the expression of OsCHS and OsCHI. For example, the transcripts of OsR498G0102785500 (LOC_Os01g74410, MYB48) and OsR498G1221897800(LOC_Os12g37970), whose abundance decreased more significantly in YB at 30 DAF, may be involved in the positive regulation of flavonoid synthesis (Fig. 5a, b). By contrast, the transcripts of OsR498G0102429100 (LOC_Os01g65370, MYB3), OsR498G0408899800 (LOC_Os04g43680, MYB4), OsR498G0714565800(LOC_Os07g37210) and OsR498G0101613100(LOC_Os01g45090), whose abundance increased more significantly in YB at 30 DAF, may be involved in the negative regulation of flavonoid synthesis (Fig. 5a, b). We subsequently used protoplast transient assays to determine the activity of different MYB TFs in the activation of OsCHS and OsCHI transcription, with blank GFP vectors used as negative controls. The expression of endogenous OsCHS and OsCHI was significantly activated when upon transfection with the 35S:OsR498G0101613100 plasmid compared with other plasmids (Fig. 5c, Supplemental Figure S5). Comparison of OsR498G0101613100, OsCHS and OsCHI expression between YN and YB showed that expression of OsCHS and OsCHI were upregulated at 30 DAF in the YB 30 versus YN 30 comparison, while expression of OsR498G0101613100 was downregulated at 30 DAF in this comparison (Fig. 5d). In addition, the promoter activity analysis also showed that the presence of OsR498G0101613100 could significantly suppress the LUC/ REN ratio (Fig. 5e), suggesting that OsR498G0101613100 could essentially inhibit the transcriptional activity of OsCHS and OsCHI. In addition, we also generated knockout lines of OsR498G0101613100 in rice by using CRISPR-Cas9-mediated gene editing (Fig. 5f). Functional analysis of two homologous lines showed antisenescence phenotype after darkness treatment (Fig. 5g). These results indicated that OsR498G0101613100 plays a negative role in regulating the expression of OsCHS and OsCHI, which may regulate the accumulation of metabolites in flavonoid biosynthesis pathway.

Discussion
High-quality rice with high yield, resistance, good quality and other good characteristics is the goal of rice breeding. Leaf colour transformation is one of important traits that been considered by breeders in the late stage of crop development, because it is actually a process of leaf senescence, in which the leaf colour changes from green to yellow, the photosynthetic capacity of leaves gradually decreases and the chlorophyll constantly degrades. It is estimated that the rice yield can be increased by 2% in theory if the life of functional leaves is prolonged one day at maturity stage 32 . Therefore, it is of great importance to cultivate anti-senescence rice varieties for stable and increased yield of rice. Yuenong Simiao is a conventional rice variety with high yield, disease resistance www.nature.com/scientificreports/ and senescence resistance, and therefore, is of practical significance for theoretical research and application of rice breeding to analyze its anti-senescence mechanism and explore the regulatory factors in the process of senescence. Leaf senescence is the last stage of plant leaf development, but it is not the least important stage. During this period, a large number of functional genes and metabolites change rapidly, which makes the plant leaves turn from green to yellow, accompanied by the gradual transfer of nutrients to the seeds. Our study combined metabolomics with the transcriptomics to analyse metabolism-related genes and metabolites involved in the process of rice leaf senescence. We found that, at the late stage of ageing, the contents of organic acids and amino acids in the plants decreased, which was consistent with findings associated with plant ageing and development 33 . At the late stage of development, nitrogen is used for the sink/source transition, and the supply of nitrogen from leaves is weakened; the levels of the majority of amino acids and organic acids during this time have been studied 16 . On the other hand, at 30 DAF, the anti-senescence variety YN accumulated more flavonoids, alkaloids and lipids than did YB, which may be due to the specificity of the YN variety. In the senescence process of Arabidopsis rosette leaves, contents decreased for chloroplast-localized lipid species; for example, the monogalactosyl diacylglycerols (MGDGs) and digalactosyl diacylglycerols (DGDGs) in the chloroplasts were degraded continuously during this process. Other lipid species, such as diacylglycerols, triacylglycerides (TAGs), ceramides, and glucosylceramides, and lysolipids such as lysophosphatidylethanolamines and lysophosphatidylcholines continuously increased during leaf senescence, which showed much higher accumulation in siliques than in senescent leaves 34 . In our study, the content of several free fatty acids increased at 30 DAF, including the following: lauric acid and 9,10-dihydroxy-12-octadecenoic acid; some sphingolipids, such as D-erythro-dihydrosphingosine and 4-hydroxysphinganine; and lysophosphatidylcholines, such as LysoPC 17:0 and LysoPC 18:0. However, lysophosphatidylcholines such as LysoPC 16:0, LysoPC 18:4, and LysoPC 18:3 increased the most at 30 DAF in YN (Supplemental Table S5). These differences in lipid metabolite content patterns may be related to the difference in senescence patterns between monocotyledons and dicotyledons.
Through integrated analysis of transcriptional and metabolic data, we focused on changes in the flavonoid biosynthesis pathway. Flavonoids have been reported to be involved in reducing oxidative damage in plants 35,36 . The levels of flavonol glycosides and anthocyanins can be increased by overexpressing flavonoid structural genes such as CHS and DIHYDROFLAVONOL-4-REDUCTASE (DFR), which resulted in a reduction in ROS accumulation and enhanced tolerance to salt stress in rice 37 . In this pathway, phloretin, luteolin and eriodictyol, which were previously reported to be involved in attenuating oxidative stress, were detected to accumulate the most in YN 30 DAF (Fig. 4c, d) [38][39][40][41] . These compounds are antioxidant regulators in human and mouse cells, but their activity in plants has not been reported. Nevertheless, we speculate that these compounds may be associated with the anti-senescence phenotype of YN. However, apigenin, naringenin chalcone, and epigallocatechin also accumulated at 30 DAF. Due to their small differences between the two varieties, we did not focus on these compounds.
MYB  [42][43][44] . GbMYBFL has been reported to enhance the accumulation of flavonoids 45 , DmMYB1 and CsMYB4a negatively regulate the synthesis of flavonoids 46 , GtMYBP3 and GtMYBP4 activate the expression of flavonol synthesis genes 47 , and so on. However, research on the relationship between MYBs and leaf senescence is limited. OsMYB102 inhibits ABA accumulation and downregulates ABA signalling responses to delay rice leaf senescence 48 . We found that the MYB TF OsR498G0101613100 may negatively regulate flavonoid biosynthesis to affect the antioxidant content. Furthermore, relationships between MYB TFs and the flavonoid pathway have been reported in citrus species. CsMYBF1 functions in controlling flavonol and hydroxycinnamic acid biosynthesis; interestingly, CsMYBF1 could activate the CHS gene promoter in Citrus but not in tomato 49 . We found that the MYB TF OsR498G0101613100 may inhibit the expression of CHS and CHI, although whether this regulation occurs directly requires further experimental verification. CHI was previously reported to be involved in the regulation of hull and internode colour 50 . However, we detected the expression of CHI in flag leaves. Further functional experiments are needed to confirm whether OsR498G0101613100 regulates CHI directly during leaf senescence.

Methods
Plant materials and growth conditions. Leaf senescence assays were conducted under field conditions. For RNA-seq and gene expression analysis, various flag leaves were collected starting from the heading stage. To determine the precise expression patterns in YN and YB, flag leaves were harvested every 10 days after flowering. To generate the knockout lines of OsR498G0101613100, the gRNA was designed to target the Exon 1 of this gene. The target sequence (+ PAM) is 5′-CAC GGG CCC TTG CGC ACC GCCGG-3′. The constructs were then introduced into pYLCRISPR/Cas9Pubi-H vector and then transformed into the japonica rice Zhonghua11 cultivar by the Agrobacterium-mediated method 51 . Chlorophyll and leaf SPAD measurements. Chlorophyll was extracted from 100 mg of leaf tissue (fresh weight), after which its content was determined by measuring the absorbance at 470 nm, 645 nm and 663 nm using a Varioskan Flash spectrophotometer (Thermo Scientific) as described previously 52 . The SPAD values of the flag leaves were determined by a SPAD-502 m. RNA-seq analysis. Total RNA was extracted from the flag leaves of rice plants at 0 DAF and 30 DAF.
Three replicate RNA-seq libraries were prepared from YB0, YB30, YN0 and YN30. A total of 12 libraries were sequenced separately using a BGISEQ-500 sequencer (BGI-Shenzhen, China). Raw sequencing reads were cleaned by removing adaptor sequences, reads containing poly-N sequences, and low-quality reads. Approximately 40 M clean reads were mapped to the R498 reference genome (http:// mbkba se. org/ R498/) using HISAT2/ Bowtie2 tools 53,54 . After the data were mapped, normalization was performed, and then FPKM (fragments per kilobase per million mapped reads) values were calculated using RESM software 55 . As previously described 56  Sample extraction and analysis of metabolites by ultra-performance liquid chromatography and tandem mass spectrometry. The flag leaves of rice plants at 0 DAF and 30 DAF were first freeze dried. Sample extract analysis, metabolite identification and quantification were performed at Wuhan MetWare Biotechnology Co., Ltd. (www. metwa re. cn), using a widely targeted metabolome method and in accordance with their standard procedures, as described in full previously 57,58 . Metabolite data analysis was conducted with Analyst 1.6.1 software (AB SCIEX, Ontario, Canada). The supervised multivariate method of partial least squares-discriminant analysis (PLS-DA) was used to maximize the metabolome differences between the pairs of samples. The relative importance of each metabolite to the PLS-DA model was checked using the variable importance in projection (VIP). Metabolites with a VIP ≥ 1 and a fold-change ≥ 1.5 or a fold-change ≤ 0.67 were considered differentially accumulated metabolites for group discrimination in the YB0 versus YB30 and YN0 versus YN30 sample comparisons.

Integrated metabolome and transcriptome analysis. Differentially accumulated metabolites and
DEGs in KEGG pathways were selected for integrated analysis. Differentially accumulated metabolites used for the integrated analysis were selected according to the information in Supplemental Table S6, and the DEGs used for the integrated analysis were selected from the top 15 metabolic pathways that were significantly enriched. We then analysed the overlapping metabolites and differentially expressed genes within these KEGG pathways.
Protoplast isolation and transfection. Ten-day-old rice plants were used for protoplast isolation according to a previously described procedure 59,60 . Briefly, the rice seedlings were cut into 0. Staining for hydrogen peroxide DAB-based analysis. Hydrogen peroxide was visually detected in leaves using DAB staining following the method described before 61 . The leaf samples were soaked in 1 mg/ml 3,3′-diaminobenzidine (DAB) for 8 h in the dark at 25 °C before being heated in a boiling 95% ethanol bath for 10 min for destaining. Finally, the samples were equilibrated for 4 h in ethanol at 25 °C.
Promoter activity analysis. The 1384 bp and 2000 bp DNA fragments upstream of the starting codons of CHS and CHI genes, respectively, were amplified by polymerase chain reaction. The resulting promoter frag-