Integrative analysis of transcriptome and metabolome reveals flavonoid biosynthesis regulation in Rhododendron pulchrum petals

Background Color is the major ornamental feature of the Rhododendron genus, and it is related to the contents of flavonoid in petals. However, the regulatory mechanism of flavonoid biosynthesis in Rhododendron pulchrum remains unknown. The transcriptome and metabolome analysis of Rhododendron pulchrum with white, pink and purple color in this study aimed to reveal the mechanism of flavonoid biosynthesis and to provide insight for improving the petal color. Results Flavonoids and flavonols are the major components of flavonoid metabolites in R.pulchrum, such as laricitrin, apigenin, tricin, luteolin, isoorientin, isoscutellarein, diosmetin and their glycosides derivatives. With transcriptome and metabolome analysis, we found CHS, FLS, F3’H, F3′5’H, DFR, ANS, GT, FNS, IFR and FAOMT genes showed significantly differential expression in cultivar ‘Zihe'. FNS and IFR were discovered to be associated with coloration in R.pulchrum for the first time. The FNS gene existed in the form of FNSI. The IFR gene and its related metabolites of medicarpin derivatives were highly expressed in purple petal. In cultivar ‘Fenhe', up-regulation of F3’H and F3′5’H and down-regulation of 4CL, DFR, ANS, and GT were associated with pink coloration. With the transcription factor analysis, a subfamily of DREBs was found to be specifically enriched in pink petals. This suggested that the DREB family play an important role in pink coloration. In cultivars ‘Baihe', flavonoid biosynthesis was inhibited by low expression of CHS, while pigment accumulation was inhibited by low expression of F3′5'H, DFR, and GT, which led to a white coloration. Conclusions By analyzing the transcriptome and metabolome of R.pulchrum, principal differential expression genes and metabolites of flavonoid biosynthesis pathway were identified. Many novel metabolites, genes, and transcription factors associated with coloration have been discovered. To reveal the mechanism of the coloration of different petals, a model of the flavonoid biosynthesis pathway of R.pulchrum was constructed. These results provide in depth information regarding the coloration of the petals and the flavonoid metabolism of R.pulcherum. The study of transcriptome and metabolome profiling gains insight for further genetic improvement in Rhododendron. Supplementary Information The online version contains supplementary material available at 10.1186/s12870-022-03762-y.

The genus Rhododendron contains brightly colored flowers and is of great use in horticulture [21]. Flavanols and anthocyanins are the main pigments in Rhododendron [22], their quantities determine the petal color ranging from light pink to violet [23][24][25][26][27][28]. At present, several studies have been done on genes and transcription factors involved in the flavonoid biosynthesis pathway in Rhododendron [29][30][31]. Rhododendron pulchrum Sweet is a species of Rhododendron belonging to the Ericaceae family and is widely cultivated in temperate Europe, Asia and North America [32]. It is highly resistant and widely cultivated in the Yangtze Delta region [33]. Based on the study of Qian et al., the main anthocyanin in purple petals is Peonidin, while Pelargonidin is the main component in pink petals, and no anthocyanin is detected in white petals of R.pulchrum [34]. There were 149 flavonoids and their derivatives identified in R.pulchrum cultivars with different colors, among them flavone and C-glycosylated flavone were the most important flavonoid metabolites [35]. In one transcriptome analysis of R.pulchrum petals, it was found that most differentially expressed genes (DEGs) were located between flower buds and early flowering, while transcripts for the genes MYC2, TIR1, CYCD3, COL-1, and EIN3 peaked at flower buds [36]. However, the regulatory mechanism of flavonoid biosynthesis in R.pulchrum remains unclear.
Recent studies have used integrated multi-omics analyses to identify and analyze traits related genes and pathways in plants [5,6,[37][38][39][40][41]. A combined analysis of transcriptomes and metabolomes can directly and accurately reflect organismal change. We performed the analysis of metabolome and transcriptome to understand the flavonoid biosynthesis pathway in petals of R.pulchrum for the first time. The findings of this study explain the coloring mechanism and provide a theoretical basis for color breeding of R.pulchrum.

Seven thousand seven hundred thirteen differentially expressed genes were found among three cultivars of R. pulchrum
To further investigate the molecular basis for flavonoid biosynthesis, transcriptome analysis was carried out to identify differentially expressed genes (DEGs) in petals of R.pulchrum. 61.84 Gb of clean data was produced from the petals of R.pulchrum (Table S3). The transcriptome data were then mapped to the genome of Rhododendron simsii and the mapping rate ranged from 84.21% to 85.27%. A total of 7713 DEGs were identified using the cut-off values of |log 2 Fold change|≥ 1 and p Value < 0.05. In the comparison group between cultivars 'Baihe' (BMJ) and 'Fenhe' (FMJ), 381 DEGs were upregulated while 1320 were down-regulated. Between the cultivars 'Zihe'(ZMJ) and 'Baihe'(BMJ) group, 1390 were up-regulated and 1278 were down-regulated. There were 2117 up-regulated genes and 1227 down-regulated genes in cultivar 'Fenhe' (FMJ) and 'Zihe' (ZMJ) group (Fig. 4).

Comparative analysis of genes enriched by GO and KEGG in three comparison groups
To understand their biological functions and gene interactions, 5156 of 7713 DEGs were annotated to GO database. 1208 and 1739 DEGs were categorized into 50 functional groups when cultivar 'Baihe' was compared with 'Fenhe' and 'Zihe' , respectively. The number become to 2209 DEGs and 52 functional groups when comparing cultivar 'Fenhe' and 'Zihe' . Three major GO categories were 'Biological process' , 'Cellular component' , and 'Molecular function' . In detail, 'Cellular process' , 'Singleorganism process' , and 'Metabolic process' were the most enriched terms in 'Biological process' , 'Cell' , 'Cell part' , and 'Organelle' were the most enriched terms in 'Cellular component' , while 'Binding' , 'Catalytic activity' , and 'Transporter activity' were the most enriched terms in 'Molecular function' (Fig. 5). 22 anthocyanin biosynthetic processes and 29 flavonoid biosynthetic processes related genes were identified in the 'Biological process' , while two anthocyanin biosynthetic processes and one flavonoid biosynthetic process related genes were identified in the 'Molecular function' .
Moreover, to identify the metabolic pathways involved in the flavonoids' biosynthesis process, 1433 DEGs were annotated with the KEGG database. 326, 472, and 635 DEGs were enriched in 19 pathways belonging to 6 categories ( Fig. 6). According to KEGG analysis, the top pathways were the metabolic pathways for carbohydrates, lipids, and other secondary metabolites. The classification showed that a large number of genes were enriched in the 'Phenylpropanoid biosynthesis'(ko00940) and 'Flavone and  (Table S4). Our study found DEGs enriched in the 'Brassinosteroid biosynthesis pathway'(ko00905) in three comparison groups, and 22 DEGs were enriched in 'Plant hormone signal transduction pathway'(ko04075) between cultivar 'Baihe' and 'Fenhe' group. Also, 9 DEGs were enriched in 'Photosynthesis'(ko00195) between cultivars 'Zihe' and 'Baihe' group (Table S4). It was reported that plant hormones and light signals were closely related to the synthesis of anthocyanins [44,45]. Our results indicated that plant hormones play an important role in white to pink coloration, while light signals play an active role in white to purple coloration in R.pulchrum.

Modeling of a new color-related flavonoid biosynthesis pathway in R.pulchrum
Based on the expression patterns of key genes and metabolites, a model of the flavonoid biosynthesis pathway in R.pulchrum was constructed. In the modeling of cultivar 'Baihe' , we speculated the inhibition of CHS gene led to a decrease in chalcone and naringin at the early stage (Fig. 8a). FLS, F3′5'H, and F3'H competed for the same substrate of dihydrokaempferol. Then, up-regulated of FLS and F3'H genes accompanied by down-regulated of F3′5'H gene drive the synthesis of kaempferol and quercetin derivatives in cultivar 'Baihe' (Fig. 8). Inhibition of F3′5'H, DFR and GT genes led to an inactivation of pigment accumulation, which resulted in the appearance of white petals (Fig. 9).
In the modeling of cultivar 'Zihe' , CHS gene showed a higher expression level in 'Zihe' than in other two cultivars. Previous studies have reported that the high expression of CHS gene is responsible for purple anthocyanins in Freesia hybrids [47] and wheat [48].  (Fig. 9).
In the modeling of 'Fenhe' , coumaric acid expression was higher, but the low expression of the 4CL gene restrained the accumulation of 4-Coumaroyl-CoA at the early stage of flavonoid synthesis pathway (Fig. 8b). The expression of FLS gene is inhibited, meanwhile the F3'H and F3′5'H genes were highly expressed, which led to overaccumulation of dihydromyricetin and dihydroquercetin (Fig. 8a). Because of the inhibition of the DFR, ANS, and GT genes, the delphinidin and cyanidin synthesis pathways are shut down leading to the formation of light purplish pink petals in cultivar 'Fenhe' . (Fig. 9).
To validate the key genes related to flavonoid biosynthesis, nine genes were randomly selected for further verification with RT-PCR (Table S6). It was found that six genes (4CL, CHS, ANS, F3 ′5' H, GT, Viodelphin-GT) were up-regulated in cultivar 'Zihe' , three genes (F3H, FLS, LAR) were up-regulated in cultivars 'Baihe' and 'Fenhe' . These results were consistent with what we discovered from RNA-seq (Fig. S4).   [49,50]. CHS, F3′5'H, ANS, and Viodelphin-GT were highly expressed, and LAR had a low expression in purple petals, and the metabolites associated with these genes were also enriched in the petals (Table S5). Our integrative analysis suggested that differential expression of CHS, FLS, F3' H, F3 ′5 'H, DFR, GT, FAOMT in the flavonoid biosynthesis caused different colors between white and purple. In addition, F3' H and F3 ′5 'H upregulated, 4CL, DFR, ANS and GT down-regulated are significantly associated with pink coloration.

Color-specific modules and hub genes were identified by using WGCNA
To further validated the flavonoid biosynthesis related key genes, a WGCNA (Weighted Correlation Network Analysis) was performed [51]. After filtering out the genes with a low expression (FPKM < 0.05), 4066 genes were obtained for the WGCNA. 12 distinct modules were identified (labeled with different colors) shown in the dendrogram in Fig. 10a. Heat map of the 12 modules are shown in Fig. 10b. The lightcyan1, mediumpurple3 and grey60 modules with 903, 41 and 658 genes, were highly associated with the cultivars 'Zihe' (r = 0.97, 0.

Discussion
Multi-omics analysis can be used to study the molecular mechanisms of plant coloration [36]. Multiple metabolites are mainly accumulated in flavonoid biosynthesis, flavonoid and flavonol biosynthesis, and anthocyanin biosynthesis pathways in R.pulchrum. We provided a model illustrating the color related key genes and metabolites and their regulatory mechanism in R. pulcherum.

Floral color and key genes in flavonoid biosynthesis of R. pulchrum
At the early stage of flavonoid biosynthesis, 4CL and CHS genes showed differential expression levels among three cultivars, but not happened with CHI gene. 4-Coumarate CoA Ligase (4CL) participates in the metabolism of a variety of secondary compounds in plants, such as lignin, coumarin, flavonoids, anthocyanins, and phenols [52][53][54]. We had identified four genes that encode 4-Coumarate CoA Ligase, three of them (Rhsim10G0185800, Rhsim10G0199800, and Rhsim10G0199000) are highly expressed in cultivar 'Zihe' , and the other one (Rhsim10G0200300) is highly expressed in cultivar 'Fenhe' . The high expression of these genes contributed to the accumulation of more precursor compounds in cultivar 'Zihe' . Chalcone synthase (CHS), as a key initiating enzyme in the flavonoid synthesis pathway, exists in a multi-gene family in most species [55]. CHS expression was significantly lower in cultivar 'Baihe' compared with 'Fenhe' and 'Zihe' , which resulted in the inhibition of flavonoid biosynthesis in 'Baihe' at the early stage. CHI gene did not show differential expressions in different color petals. As Takos et al. discovered there was no difference in CHI expression between the peels of red and non-red apple varieties, but when CHI was completely absent, anthocyanin synthesis was limited [56]. Therefore, they believed CHI was not a key enzyme for the catalytic synthesis of chalcone naringin, but rather acted as rate-limiting enzyme in the flavonoid pathway.
At the late stage of flavonoid biosynthesis, FLS, F3H, F3 'H, DFR, LAR, F3 ′5' H, ANS, GT, and violdelphin-GT genes showed differential expression in three petals. Flavonol synthase (FLS), a key enzyme in the biosynthesis of flavonols, catalyzes the production of kaempferol and myricetin from dihydrokaempferol and dihydromyricetin respectively. A deletion of FLS or an alteration of gene expression can affect the coloration of plant organs [57,58]. As a result of high FLS expression in cultivar 'Baihe' , high levels of kaempferol and quercetin derivatives were found. FLS play a more active role in catalyzing the dihydrokaempferol than F3'H [59], resulting in the accumulation of kaempferol and quercetin derivatives in 'Baihe' . It was reported that FLS and DFR compete for substrates to synthesize flavonols, and an upregulation of FLS and a downregulation of DFR may result in a reduction of anthocyanin accumulation [60], this is in line with our findings. There are MYB transcription factor binding sites and light signal response sites in most FLS gene promoters, indicating that FLS gene transcription is regulated by MYB transcription factors and light signals [61,62]. In this study, the DEGs are enriched in the 'light signaling pathway' in 'Fenhe' and 'Zihe' group (Table S4). However, whether the FLS gene is regulated by light signal and MYB transcription factors in R.pulchrum needs further study. The flavonoid 3 'hydroxylase (F3'H) is a critical enzyme in producing Purplish Red [63]. As well as catalyzing kaempferol to form quercetin, it also catalyzes dihydrokaempferol to form dihydroquercetin. There are two genes that encode F3 'H (Rhsim01G0176800 and Rhsim08G01436 00). Rhsim01G0176800 was highly expressed in cultivar 'Baihe' , whereas Rhsim08G0143600 was highly expressed in cultivar 'Fenhe' . We speculate that these two genes may have distinct regulatory pathways, Rhsim01G0176800 catalyzes kaempferol to produce quercetin, while Rhsim08G0143600 catalyzes dihydrokaempferol to produce dihydroquercetin, further validations are required. F3 ′5' H is a key gene associated with the synthesis of blue color. In our study, F3 ′5' H (Rhsim13G0208100, Rhsim13G0208200) are highly expressed in cultivar 'Zihe' , leading to the synthesis of dihydromyricetin, a substrate of delphinidin and malvidin glycosides. Both LAR and ANS play an important role in the synthesis of proanthocyanins and anthocyanins, which are downstream genes in the flavonoid biosynthesis pathway. It is known that leucocyanidin is the substrate of LAR and ANS genes, which could catalyze leucocyanidin into catechins and cyanins, respectively [45]. LAR gene expression is affected by temperature, light, and hormones [64], and ANS gene expression is closely linked to plant organ coloration [65,66]. In our finding, LAR was highly expressed in cultivars 'Fenhe' and 'Baihe' , resulting in the synthesis of more Catechin derivatives using leucocyanidin. ANR and GT were competing for substrate cyanins. The enzyme ANR catalyzes the production of epicatechin from cyanins, which is also affected by light, temperature, and hormones [64]. By catalyzing the glycosylation of anthocyanin, flavonoid 3-O-glucosyltransferase (GT) is an important enzyme in the last stage of anthocyanin synthesis [67]. Interestingly, there was no difference in ANR expression among the three cultivars, but cultivar 'Zihe' exhibited a higher concentration of epicatechin derivatives, suggesting that the use of cyanidin as a substrate to synthesize epicatechin derivatives resulted in accumulation of less red pigmentation.

Floral color and MYBs, AP2/ERFs transcription factors in R. pulchrum
MYB is one of the largest transcription factor families in plants. It mainly regulates the flavonoid biosynthesis via two pathways. Firstly, MYBs activate the expression of flavonoid synthesis genes by interacting with bHLHs and WD40s to form the MBW complex [68]. FaMYB5 may act as negative regulators of PA/ anthocyanin biosynthesis in strawberry [69]. MYB5 (Rhsim01G0174700) was highly expressed in cultivars 'Baihe' and 'Fenhe' , suggesting that it may negatively regulate anthocyanin synthesis. Secondly, MYBs can regulate the expression of structural genes independently to influence flavonoid synthesis. As reported by Luo et al., FTMYB15 was transferred from Tartary Buckwheat to Arabidopsis thaliana, the pigment accumulation was found in the leaves and seed coat of the transgenic plant [70]. MYB15 (Rhsim04G0061700) was discovered with downregulation in cultivar 'Zihe' . It is suggested that it may inhibit purple coloration in R. pulcherum. According to previous studies, MYB6 regulates flavonol biosynthesis through activating the promoters of F3H and FLS [71]. Based on our findings, MYB6 (Rhsim07G0002700) was highly expressed in 'Fenhe' , and rutin content in 'Fenhe' was significantly higher than in 'Zihe' . It is indicated MYB6 was also involved in the pink coloration of petals through the regulation of flavonol synthesis. The MYB101 participated in pollen tube reception and leaf development [72][73][74]. According to our results, MYB101 (Rhsim06G0163100) was highly expressed in cultivar 'Fenhe' , suggesting that MYB101 may be responsible for regulating flavonoids in R.pulchrum with pink petals.
UV-B can be used as a signal to regulate plant development, such as regulating photomorphogenesis and promoting the accumulation of flavonoids and anthocyanins [75][76][77]. MYB73 is an important transcription factor associated with plant hormones [78]. UVR8 interacts with MYB73/MYB77 (MYB DOMAIN PROTEIN 73/77) in a UV-B-dependent manner, it can affect plant growth and development through light and auxin signal [79]. When MYB73 (Rhsim05G0088000) was highly expressed in cultivar 'Zihe' , MYB73 may contribute to the accumulation of flavonoids and anthocyanins in R.pulchrum by regulating the photomorphogenesis process. MYB61 (Rhsim11G0106400) was highly expressed in dark flowers, suggesting that it may regulate the coloration in R.pulchrum.
The AP2/ERF family is one of the largest transcription factors families in plants and regulate primary and secondary metabolism [80][81][82]. It has been shown that AP2/ERF can regulate the biosynthesis of anthocyanins in plants. MdMYB1 is a positive regulator of anthocyanin biosynthesis, and MdERF38 interacts with MdMYB1 to promote the binding of MdMYB1 to its target gene in the drought-induced conditions [83]. DREB subfamily, which belongs to AP2/ERF family, was only found to be highly expressed in cultivar 'Fenhe' . But further research is needed to confirm whether DREB can regulate flavonoid biosynthesis in pink petals by affecting the interactions of TFs under certain induction conditions. Rhsim07G0115900 and Rhsim03G0161100, which encode ERF5, were discovered in this study. Rhsim07G0115900 was lowly expressed in cultivar 'Zihe' when compared with 'Fenhe' and 'Baihe' , while Rhsim03G0161100 was highly expressed in cultivar 'Baihe' when compared with 'Fenhe' . It is suggested that ERF5 might be involved in regulation of flavonoid biosynthesis in R. pulcherum.

The FNS and IFR genes were identified in R.pulchrum
Flavone synthases had two different subtypes, Flavone synthase I (FNSI) and Flavone synthase II (FNSII). FNSI is a nonheme iron-dependent a-ketoglutarate dioxygenase, which is found in the cytosol [84]. As a member of the P450 family, FNSII is found in the inner membrane of cells. Both enzymes can convert naringin to apigenin, and they are not simultaneously in the same species. Researcher has reported that Medicago sativa and soybean contain only FNSII [71,85]. Metabolome analysis revealed that the apigenin derivatives had a lower expression in cultivar 'Zihe' than in 'Baihe' and 'Fenhe' . It was found that Rhsim03G0216000, which encodes FNS, was lowly expressed in cultivar 'Zihe' . Sequence alignment results revealed that Rhsim03G0216000 had a great similarity (74%) to FNSI(AT5G24531.1). Based on these findings, it is likely that FNS is in the form of FNSI in R.pulchrum.
Isoflavone is a type of phytoestrogens, which have only been found in some species of Papilionoideae [86,87]. The isoflavone reductase (IFR) catalyzes the conversion of isoflavones to Medicarpin, which is the first step in isoflavone decomposition [88,89]. The IFR gene (Rhsim07G0136900) was highly expressed in cultivar 'Zihe' , meanwhile metabolite of isoflavones (Medi-carpin3-O-(6'-malonylglucoside)) highly enriched in cultivar 'Zihe' . Therefore, it is necessary to investigate whether the isoflavone derivatives cause the purple coloring of R.pulchrum. Medicalarpin is a natural phytoalexin of the pterocarpan subfamily. It has been found in many traditional Chinese medicines, such as Pueraria Lobata [90]and Hedysari Radix [91]. The studies have reported that it has various biologically benefits, including stimulation of bone regeneration, induction of apoptosis, and induction of lipolysis [92]. Interesting, many medicarpin derivatives have been found in cultivar 'Zihe' for the first time, suggesting that this cultivar is valuable for medicine research. It has been observed that a few purple spots are sometimes on the white petals of R.pulchrum, we discovered delphinidin and malvidin were accumulated in the petals. The blockage of glycosylation related genes in the flavonoid synthesis pathway may explain why white flowers accumulated some delphinidin and malvidin.

Conclusion
The transcriptome and metabolome analysis were used to identify genes, transcription factors, and metabolites associated with the coloration of R.pulchrum. For known color-related transcription factors, the DREB family were found to play a critical role in the formation of the pink color. Additionally, our results have shown that FNS and IFR were highly expressed in the flavonoid biosynthesis pathway, and that FNS exists in the form of FNSI in R.pulchrum. The IFR gene and its related metabolites of medicarpin derivatives were highly expressed in cultivar 'Zihe' , indicating that IFR involved in the purple coloration of R.pulchrum. CHS, FLS, F3' H, F3 ′5 'H, DFR, GT and FAOMT genes were correlated with white and purple coloration, whereas F3' H and F3 ′5 'H up-regulated, 4CL, DFR, ANS and GT down-regulated were associated with pink coloration. We provided a model illustrating the color related key genes and metabolic and their regulatory mechanism of R. pulcherum. The study provided some basis for breeding new cultivars with different colors of petals.

Plant materials and treatment
The white-flowered R.pulchrum Sweet ('Baihe'), pinkflowered R.pulchrum Sweet ('Fenhe') and purpleflowered R.pulchrum Sweet ('Zihe') were grown at the Rhododendron resource center in the Shanghai Botanical Garden (31°15′ N, 121°45′ E), Shanghai, China (Fig. 12). After collecting the petals from 'Baihe' , 'Fenhe' and'Zihe' , all materials were frozen in liquid nitrogen and stored at − 80 °C for RNA and metabolite extraction. In the transcriptome analysis of 'Baihe' , 'Fenhe' and 'Zihe' , for each color, we gathered three samples as biological replicates, then each sample randomly selected three blooms petals from the same plant and pooled into one for RNA-Seq. For metabolism analysis, we had made five biological replicates from five individual plant for each color.

Sample preparation and metabolite extraction
Each color of petal was collected from five individual plants. 1 mL of methanol: water (7:3, v:v) was added to 80 mg samples. 2-chloro-l-phenylalanine (20μL, 0.3 mg / mL) was added as an internal standard. All samples were precooled at -20 °C for 2 min, pulverized with grinder at 60 HZ for 2 min, then ultrasonicated for 30 min and placed at -20 °C for 20 min. The samples were centrifuged for 10 min (13,000 rpm, 4 °C). The supernatants (150μL) from each tube were collected with crystal syringes, filtered through 0.22 μm microfilters, and injected into LC vials. Until the liquid chromatography-mass spectrometry (LC-MS) analysis was performed, the vials were stored at -80 °C. A pooled QC sample was prepared by mixing aliquots of all the samples. For the metabolic profiling, an ACQUITY UHPLC system (Waters Corporation, Milford, USA) in conjunction with an AB SCIEX Triple TOF 5600 System (AB SCIEX, Framingham, MA) was used. ACQUITY UPLC BEH C18 columns (1.7 μm, 2.1 × 100 mm) were used in both positive and negative mode. The binary gradient elution system consisted of (A) water (containing 0.1% formic acid, v/v) and (B) acetonitrile (containing 0.1% formic acid, v/v). The gradient program consisted of 0 min, 95% A; 2 min, 80% A; 4 min, 75% A; 9 min, 40% A; 17 min, 0% A; 19 min, 0% A; 19.1 min, 95% A; and 20.1 min, 95% A. The flow rate was 0.4 mL/min and the injection volume were 5μL. The temperature of the column was maintained at 45 °C. During the LC-MS analysis, samples were maintained in the instrument chamber at 4 °C.

RNA-seq and annotation
For each sample, 1 g fully opened stage petals were used for extracting the RNA using the mirVana miRNA Isolation Kit (Ambion). RNA purity and quantification were assessed using a NanoDrop 2000 spectrophotometer (Thermo Scientific, USA). RNA integrity was evaluated using an Agilent 2100 Bioanalyzer (Agilent Technologies, Santa Clara, CA, USA). Following the manufacturer's instructions, libraries were constructed using TruSeq Stranded mRNA LTSample Prep Kit (Illumina, San Diego, CA, USA). Afterwards, these libraries were sequenced on the Illumina sequencing platform (Illumina HiSeq X Ten) and 150 bp paired-end reads were obtained. Raw reads were first processed using Trimmomatic [93], and low-quality reads and adaptors were filtered out. HISAT2 [94] was used to map the remaining reads to the reference genome. The differential expression analysis was carried out using the DESeq R package [95]. According to the DESeq analysis, genes with adjusted P-values of < 0.05 were significantly differentially expressed. A hierarchical cluster analysis of DEGs was performed to analyze gene expression patterns. GO enrichment and KEGG [96] pathway enrichment analysis of DEGs were performed with R based on the hypergeometric tests. The annotation of gene functions was carried out using the following databases: Kyoto Encyclopedia of Genes and Genomes (KEGG), eukaryotic Clusters of Orthologous Groups (KOG), NCBI non-redundant (NR), Swiss PROT sequence protein, Gene Ontology (GO) and the homologous protein family (Pfam). Three biological replicates of each sample were performed.

RT-qPCR verification
A total of nine structural genes related to flavonoid biosynthesis were selected for validation by quantitative real-time PCR (RT-qPCR). Total RNA was extracted with Total RNA Kit (Ambion).The RNA was treated with DNase I(Invitrogen) to remove trace genomic DNA, then reverse-transcribed in a GeneAmp ® PCR System 9700 (Applied Biosystems, USA).The specific primers were designed in the laboratory and synthesized by TsingKe Biotech (Supplementary Table S6). RT-qPCR was performed using a LightCycler ® 480 II Real-time PCR with Instrument (Roche, Swiss).The PCR amplification was 94℃ for 30 s, followed by 45 cycles of 94℃ for 5 s, 60℃ for 30 s. The expression levels of mRNAs were normalized to GAPDH and were calculated using the 2 − ΔΔCt method [97]. The primers of RT-qPCR are showed in Table S6.

Metabolites-transcripts correlation and weighted gene coexpression network analysis
The Pearson correlation algorithm was used to calculate the correlation between genes and metabolites. DEGs and DEMs have been mapped simultaneously to the KEGG Pathway database to obtain information on their pathways. The correlation analysis was performed by R 3.5.1 using Pearson correlation analysis. Coexpression networks was performed with WGCNA analysis [51].