Regulation mechanisms of flavonoids biosynthesis of Hancheng Dahongpao peels (Zanthoxylum bungeanum Maxim) at different development stages by integrated metabolomics and transcriptomics analysis

Flavonoids have strong free radical scavenging and antioxidant capacity. The high abundance of flavonoids in Chinese prickly ash peels have many benefits to human health. In this study, ‘Hancheng Dahongpao’, a main cultivar, was taken as materials to investigate the flavonoids biosynthesis mechanism of Zanthoxylum bungeanum Maxim at three key development stages by integration of metabolomics and transcriptomics analysis. A total of 19 differentially accumulated metabolites were identified, the key flavonoids compounds were kaempferol, quercetin and their glycoside derivatives, and two major anthocyanins (peonidin O-hexoside and peonidin 3-O-glucoside). 5 gene networks/modules including 15 important candidate genes were identified, which was highly correlated with flavonoids. Among these genes, ZM-163828 and ZM-184209 were strongly correlated with kaempferol and quercetin, and ZM-125833 and ZM-97481 were controlled the anthocyanins biosynthesis. Moreover, it was shown that MYB-ZM1, MYB-ZM3, MYB-ZM5, MYB-ZM6 and MYB-ZM7 coordinately controlled flavonoids accumulation through regulating the structural genes. Generally, this study systematically revealed the flavonoids metabolic pathways and candidate genes involved in flavonoids biosynthesis and laid a foundation for the potential targets for the breeding of new valuable Chinese prickly ash cultivars.


Background
Flavonoids are a large class of secondary metabolites with strong biological activities generated in the process of long-term ecological adaptation of plants [1,2]. Based on their structural differences, flavonoids are generally divided into nine main types, including chalcones, flavonoids, flavonols, dihydroflavonoids, flavanols, flavonoid glycosides, isoflavones, proanthocyanin and anthocyanins [3][4][5]. In plants, most flavonoids exist in the form of Open Access *Correspondence: 284835155@qq.com; zhengtyhy@163.com glycosides, and a few in the form of free state. The glycosylation sites are usually C-3 and C-7, and the glycosylation sites include glucose, raffinose, galactose, rhamnose, among which glucose is the most common. A variety of modifications such as glycosylation and methylation of flavonoids greatly enriched the derivatives varieties, and also improved the stability of these substances in plants [6,7]. Flavonoids are widely distributed in flowers, fruits and leaves of many plants, about 4000 flavonoids have been found [8]. Flavonoids, with many types and complex structures, exhibit tissue-specificity and are regulated by genes and biotic and abiotic factors in environment [9].
Flavonoids in plants have enhanced plant resistance and chemical defense functions, and play an active role in plant ecological adaptation. Flavonoids play an important role in plant growth, development, flowering, fruit and antibacterial disease prevention. For example, flavonoids in roots can regulate root growth and affect mineral absorption, thereby affecting plant growth [10][11][12]. Flavonoids are usually related to plant stress resistance, which is beneficial to the absorption of reactive oxygen species (ROS) and prevent UV damage to plants [13,14]. Flavonoids have many effects, such as anti-cardiovascular and cerebrovascular diseases, antioxidant, anti-aging, analgesic, anti-inflammatory, antitumor, antiviral and so on, which play an important role in human health. Numerous studies have confirmed that flavonoids have health functions, and further were processed into products with special functions such as health food and drugs [15].
Chinese prickly ash (Zanthoxylum bungeanum Maxim) is a kind of plant resource with high edible and medicinal value, and has great economic and ecological value. In recent years, it was widely planted in China, and its planting area expanded rapidly [28,29]. Its peel contains a variety of active ingredients with antibacterial, antitumor, antioxidant and other health effects, which are beneficial to human body [30,31]. As the main polyphenol secondary metabolites in Chinese prickly ash plants, flavonoids have great benefits for many physiological activities of the human body [32]. However, the investigation on gene regulatory networks of flavonoids in Chinese prickly ash peels was lacking. In this study, the fresh fruits of Chinese prickly ash peels (Zanthoxylum bungeanum Maxim. 'Hancheng Dahongpao') at the three key development stages were picked as experimental materials, and the genes involved in flavonoids metabolism and the differential metabolite in Chinese prickly ash peels were systematically studied to reveal the gene networks related to flavonoids metabolism.

Peel material
The 'Hancheng Dahongpao' fruits at different development stages: green stage (G), yellow stage (Y) and fullred stage (R) (Fig. 1) were selected as test materials. The 'Hancheng Dahongpao' was planted in the experimental garden of college of forestry, Northwest A & F University, Yangling City (34°15'N, 108°03'E), Shaanxi Province, China. All peels uniformity in color without signs of mechanical damage or disease (approximately 2 mm) were washed with distilled water, frozen in liquid nitrogen and stored at −80 °C for further use. Three individual biological replicates (9 samples in total) were used for metabolism and transcriptome analysis of peels at each stage.

Sample preparation and metabolite extraction
All peels were treated by vacuum freeze-drying, grinded to powder (30 Hz, 1.5 min) by a mixer mill (MM 400, Retsch). The powder (100 mg) was dissolved in 1.0 mL extract (70% methanol aqueous solution), and the dissolved samples were vortexed (10 s, 40 Hz) every 10 minutes three times and kept at 4 °C overnight. Following centrifugation at 10, 000 g for 10 min, the supernatant was filtered by microporous membrane (0.22 μm pore size) and stored in the sampling bottle for UPLC-MS/MS analysis.
The mass spectrometry data were processed by the software analyst 1.6.3 software (AB SCIEX, Ontario, Canada). The partial least squares-discriminant analysis results (OPLS-DA) could maximize the distinction between groups and help to find differential metabolites. Based on the OPLS-DA results, the differential metabolites were screened by combining fold change and the variable importance in projection values (VIP). VIP ≥ 1 and fold change ≥2 or fold change ≤0.5 were set as the selection standard differential metabolites.

RNA extraction, cDNA library construction and bioinformatics analyses
Total RNA extraction, detection and library construction from green, yellow and full-red peels were prepared by a previously described method. The cDNA libraries preparation was sequenced on the Illumina HiSeq 2000 platform, and 150 bp paired-end reads were generated. For gene-expression analysis, the clean reads were mapped to the reading of each gene by HTSeq v0.5.4p3 and the FPKM (fragments per kilobase of transcript per million mapped reads) method was used to estimate the reads count and expression levels of each gene, respectively.

Analysis of differentially expressed genes
Differentially expressed genes (DEGs) were identified based on the fold change of the FPKM values from the different samples and performed using DESeq2 v1.22.1. Following the difference analysis, the false discovery rate (FDR) control was applied to calculate the threshold of the P-value. The screening conditions for DEGs were |log2Fold Change| > = 1 and FDR < 0.05. The pathways of differentially expressed genes were annotated Based on the KEGG comprehensive database (Kyoto Encyclopedia of Genes and Genomes). WGCNA analysis (weighted gene coexpression network analysis) was carried out by WGCNA v1.69 to build a gene coexpression network based on genes coexpressed modules.

Experimental validation of candidate genes
Candidate genes expression analysis was performed with three independent biological replicates by using real-time quantitative PCR (RT-qPCR), which contained 10 differentially expressed structure genes and 3 MYB transcription factors genes. RT-qPCR analysis was performed on a Step One PULS real-time detection system detection system (ABI, Foster, CA, USA) using 2 × SG Fast qPCR Master Mix (High Rox, B639273, BBI, ABI). The candidate genes expression data was quantitatively calculated by 2 -ΔΔCT method with the reference gene UBQ. The primers and UBQ gene were shown in Table S1.

Metabolome profiling in 'Hancheng Dahongpao' peels at different development stages
Based on UPLC-MS/MS results, the metabolomics analysis revealed significant differences in in flavonoids metabolites in peels at different developmental stages. HCA (Hierarchical cluster analysis) on the accumulation patterns of flavonoids metabolites among different samples demonstrated the good repeatability in the sample group. Heatmap analysis of flavonoids metabolites showed that 9 samples were significantly separated into three clusters corresponding to three consecutive ripening stages (green stage, yellow stage and full-red stage) ( Fig. 2A). Principal component analysis (PCA) results exhibited that the G and Y stages samples were clustered together and significantly separated from the R stage samples, indicating that there were different accumulation and expression patterns of flavonoids in green, yellow and red fruits (Fig. S1). The HCA and PCA results presented that there were great differences in each treatment, and the difference between different repetitions of the same treatment was small, and the repeatability of the sample was good.

Differentially accumulated flavonoids metabolites in 'Hancheng Dahongpao' peels
To investigate the differentially accumulated metabolites (DAMs) of flavonoids in the 'Hancheng Dahongpao' fruits at different development stages (G-vs-Y, Y-vs-R and G-vs-R), the flavonoids compounds with VIP value ≥1 and fold change ≥2 or fold change ≤0.5 were selected. Based on pair-wise comparisons of samples (G-vs-Y, Y-vs-R and G-vs-R), 27, 36 and 53 DAMs in G-vs-Y (Table S2), G-vs-R (Table S3) and Y-vs-R (Table S4) were found, respectively. Due to the interaction of flavonoids DAMs in 'Hancheng Dahongpao' peels, different pathways were formed, and the flavonoids DMs were annotated and assigned to the KEGG pathways ( Fig. 2B-D). KEGG pathway enrichment analysis revealed that biosynthesis of secondary metabolites, phenylpropanoid biosynthesis, flavonoid biosynthesis, flavone and flavonol biosynthesis and anthocyanins biosynthesis were the main flavonoids enrichment pathways, which were enriched in these core-conserved metabolites. Therefore, we speculated that the DAMs in the above pathways might lead to flavonoids cause variations in in 'Hancheng Dahongpao' peels in the developmental process.

Transcriptome analysis of peels at different developmental stages
Transcriptome sequencing was performed using the 'Hancheng Dahongpao' peels in three developmental stages. After raw data filtering, sequencing error rate and GC content distribution inspection, 9 samples of 64.53 Gb clean data were obtained, and the data of each sample were higher than 6 Gb. The sequencing base Q30 was more than 92.00% (Table 2). 9 sequencing libraries were constructed and named G1, G2, G3, Y1, Y2, Y3, R1, R2 and R3. PCA results of the genes expression based on the number of fragments per kilobase of exon per million fragments mapped (FPKM) values showed that all the biological replicates clustered together, indicating the high reliability of our sequencing data. The contribution rate of PC1 and PC2 was 51.99 and 39.13%, respectively (Fig. S2). Based on correlation analysis between samples, it can be observed that there were reliable biological replicates within the group between samples (Fig. S3). PCA and correlation analysis results unraveled that the samples within the group had good biological repeatability, and there were significant differences between samples.
The transcriptome data of 'Hancheng Dahongpao' peels during the three development stages exhibit that 17,799 genes in common were expressed. There were 13,290, 13,623 and 15,140 DEGs in the comparison of G-vs-Y, Y-vs-R and G-vs-R respectively. Based on the log 2 |Fold Change| > 1 and FDR < 0.05, 8360, 4800 and 7092 genes were upregulated and 4930, 8823 and 8048 genes were downregulated in G vs Y, Y vs R, G vs R, respectively ( Fig. 3A and Table S5). In order to further understand the differential genes involved in flavonoids biosynthesis, DEGs were further assigned to KEGG pathway. In KEGG analysis, DEGs mapped to metabolic pathways occupied the largest proportion, with "biosynthesis of secondary metabolites" and "plant home signal transduction" ranking second and third (Fig. 3B-D), which demonstrated that the DEGs were mainly concentrated in the secondary metabolism process. Plant home signal transduction induced the endogenous hormones variation such as ABA and brassinosteroid (BRs), and the production and accumulation of small molecule substances such as jasmonic acid (JA) and salicylic acid (SA), and converted environmental signals into endogenous signals, and made the signal transduction and amplification. Moreover, phenylalanine and small molecules induced the expression and activity of flavonoids biosynthesis pathway enzyme genes, and further promoted the flavonoids synthesis and accumulation. In detail, the pathways of the flavonoid biosynthesis and phenylpropanoid biosynthesis involved in flavonoid accumulation were apparently enriched in the mature peels. Furthermore, a quite large number of starch and sucrose metabolism encoding genes were also enriched evidently. Soluble sugar was the basic synthetic material of plant secondary metabolism, and was the direct precursor and indirect precursor of flavonoids glycosides synthesis, regulating many plant life activities through signaling pathways. In addition, the ABC transporter encoding genes are also enriched in the G-vs-Y, Y-vs-R peels, indicating that the enhanced synthesized flavonoids and needs to be transferred to a suitable organelle for storage. However, note that antenna proteins involved in photosynthesis showed the highest rich factor during all the pathways. For each group, DEGs enrichment associated with flavonoid metabolism were mapped to flavonoid biosynthesis and phenylpropanoid biosynthesis pathway, which displayed significant upregulation in the comparison of G-vs-Y, Y-vs-R and G-vs-R.

Metabolic and gene coexpression networks in 'Hancheng
Dahongpao' peels at different developmental stages WGCNA analysis was carried out to identify coexpressed gene modules based on genes fpkm values and explore the correlation between genes in these modules and metabolites. A total of 17,799 genes were selected to use in WGCNA analysis, and the corresponding relationship between the scale independence and the mean connectivity under different thresholds were measured by the power values from 1 to 30. The scale independence and the mean connectivity of all gene adjacency functions in the corresponding gene modules showed that the best power value was 30. 7 distinct gene modules and clustering tree were built based on the correlation of gene expression (Fig. 4A). Each color in the module-level clustering tree graph represented a color corresponding to each gene on the clustering tree that belonged to the same module. The sample dendrogram and trait heatmap were constructed using the 19 differential metabolites contents at three development stage as the phenotypic data to illustrate the differential metabolites variation (Fig. 4B). The heatmap of the correlation between gene modules and differential metabolites revealed that five gene modules were highly correlated with flavonoids. A detailed description of the correlation between 19 DAMs and gene modules was present in Fig. 4C and Fig. 4D  was positively linked to gallocatechin (r 2 = 0.90) and catechin (r 2 = 0.94).

Differentially expressed genes linked to flavonoids biosynthesis in 'Hancheng Dahongpao' peels
To explore the difference of flavonoid biosynthesis in 'Hancheng Dahongpao' peels at different developmental stages, transcriptome and metabolomics data were integrated for joint analysis (Fig. 5).  UFGT (ZM-140426) in blue module were strongly upregulated, which resulted in the significant production of various anthocyanidin glycosides finally. Peonidin o-hexoside and peonidin 3-o-glucoside were extensively generated through the biosynthetic pathway, demonstrating that two compounds were the main anthocyanins in 'Hancheng Dahongpao' peels resulting in the peel color change from green to full-red. However, neither increase nor evident signals were found for some important intermediates (such as leucoanthocyanidins and anthocyanidins) in the anthocyanin pathway, indicating that these Fig. 4 A Network heatmap plot of genes subjected to co-expression module calculation. B Heat map of 19 differential metabolites contents at three development stage. C Module hierarchical clustering tree based on the correlation of gene expression. Each color in the graph represents a color corresponding to each gene on the cluster tree that belongs to the same module. D Correlation heat map between 19 differential metabolites and gene modules. Color from blue to red indicated r 2 values from-1 to 1 metabolites might be unstable in 'Hancheng Dahongpao' peels.
In G-vs-Y, CHI (ZM-141078) was the most significantly upregulated DEGs in the flavonoid pathway, increasing 84.729-fold, followed by C4H (ZM-136670 and ZM-137279), which showed 32.890-fold and 23.605fold upregulation in the yellow peels. In Y-vs-R, F3'H (ZM-158354) stood out from the DEGs with 7.655-fold change, followed by DFR (ZM-17311) and ANS (ZM-97481). In G-vs-R, the UFGT (ZM-140426) ranked first among the DEGs, increasing 13.183-fold, followed by ANS (ZM-97481) with 8.788-fold change and DFR (ZM-17311) with 6.083-fold change. These results indicated that the expression of the above structural genes was an important prerequisite for the flavonoids accumulation. In 'Hancheng Dahongpao' peels, flavonoids were usually converted to their glycoconjugates and then compartmentalized in vacuoles for storage, and these glycosylation reactions were catalyzed by glycosylation transferases. UDP-glycose/flavonoid glycosyltransferases (UFGT), a kind of glycosyltransferase, catalyzed the production of a wide range of substrates, such as anthocyanin and glycosides derivatives. The above results showed that the combined analysis of expression profile and metabolites constructed a clear flavonoid biosynthesis pathway during 'Hancheng Dahongpao' fruit development stage, and provided a strong research basis for further exploration of flavonoid metabolism in 'Hancheng Dahongpao' peels.

Analysis of transcription factors in 'Hancheng Dahongpao' peels
The flavonoids biosynthesis in 'Hancheng Dahongpao' peels was not only regulated by the structural genes discussed above, but also was controlled by transcription factors. A total of 418,382 and 604 differentially expressed transcription factors were identified in G-vs-Y, Y-vs-R and G-vs-R, respectively (Table S6). The differentially  expressed transcription factors were divided into various  families, such as WRKY, TCP, NAC, MYB, MADs-box,  HSF, HD-ZIP, GRAS, C2C2, bHLH and AP2/ERF. MYB transcription factors belonged to the R 2 R 3 MYB family, and were closely associated with flavonoids biosynthesis, which possessed obvious role in regulating flavonoids accumulation. 14 upregulated and 4 downregulated MYB transcription factors were identified in G-vs-Y, 13 upregulated and 6 downregulated MYB transcription factors were selected as significant transcription factors in Y-vs-R (Fig. 6). A sequence analysis of these MYB proteins presented that five transcription factors (MYB-ZM1, MYB-ZM3, MYB-ZM5, MYB-ZM6 and MYB-ZM7) were closely associated with flavonoids regulation (Fig. S4 and Fig. S5), in which the MYB-ZM5 (ZM-65290) exhibited a high sequence homology with the essential anthocyanin regulatory gene VvMYBA1, MYBA1 and MYBA2, involved in regulating anthocyanin synthesis. Multiple alignment analysis of the MYB proteins from various species showed that MYB-ZM1 (ZM-10723), MYB-ZM6 (ZM-166024), and MYB-ZM7 (ZM-71403) shared the same conserve R2 and R3 motif with the MYB transcription factors regulating the anthocyanin and flavonoids biosynthesis. MYB-ZM1and MYB-ZM7 showed a close relationship to ATMYB11, MYBF1, ATMYB12, and ATMYB111, which were involved in flavonoids biosynthesis. MYB-ZM6 exhibit a high sequence homology with VvMYBF1, which positively regulated flavonols biosynthesis. MYB-ZM3 (ZM-213334) was closely togethered with flavonoid MYB repressors, such as FaMYB1 and AtMYB4 [33]. Coincidentally, the MYB-ZM1, MYB-ZM3, MYB-ZM5, MYB-ZM6, and Fig. 6 A Differentially MYB expressed transcription factor encoding genes between G-vs-Y (above) and Y-vs-R (below). B All the mutual differentially expressed MYB genes between G, Y and R. C Log2 FC values of the mutual MYB genes in G-vs-Y and Y-vs-R MYB-ZM7 were selected to be the candidate transcription factors involved in flavonoids accumulation. Altogether, these results strongly indicated that MYB-ZM1, MYB-ZM3, MYB-ZM5, MYB-ZM6 and MYB-ZM7 participated in the transcriptional regulation of anthocyanins, flavonols and flavonoids in 'Hancheng Dahongpao' peels.

Correlation between DEGs and flavonoids DAMs
To further study the relationship between DEGs and DAMs in 'Hancheng Dahongpao' peels, the Pearson correlation coefficient between genes and metabolites were calculated. The results of correlation network showed a strong correlation (coefficients r > 0.8 or < −0.8) between flavonoids and their regulatory structural genes and transcription factors (Fig. S6). For example, MYB-ZM5 and the two structural genes (ANS and UFGT) were significantly positively correlated with peonidin o-hexoside and peonidin 3-o-glucoside. However, the kaempferol glucuronic acid was negatively regulated by MYB-ZM1, MYB-ZM3 and MYB-ZM6. The above results revealed that these structural genes and transcription factors played an important regulatory role in flavonoid synthesis in 'Hancheng Dahongpao' peels.

RT -qPCR validation of the transcriptomic data
To further verify the candidate genes involved in flavonoids biosynthesis, 10 structure differentially expressed genes and 5 MYB transcription factors related to flavonoids biosynthesis were selected to investigate and verify their expression patterns. Consistent with the transcriptome data, both the candidate genes, including C4H ( , were in concordance with the transcriptome analysis, which further confirmed that these genes may be associated with flavonoids accumulation (Fig. 7 and Table S7).

Discussion
Flavonoids have many physiological activities and are widely used in food. Chinese prickly ash has been widely concerned as a medicinal and food homologous plant, and its peels are rich in flavonoids and have great development and utilization value. The functional gene mining and molecular regulation mechanism of flavonoid biosynthesis have not been elucidated. It was expected to preliminarily reveal the related genes and metabolic pathways of flavonoids metabolism in Chinese prickly ash peel by multi-omics methods, so as to lay the foundation for the study of flavonoids in plants.
Flavonoids could be used as regulators in plant growth and development to protect plants against UV damage and prevent the invasion of pathogenic microorganisms. Isoflavones are also important phytoalexins [34]. Usually, flavonoids accumulated in Chinese prickly ash peels to protect chloroplasts from photoinhibitory and photooxidative damage by absorbing high-energy quanta [35]. In the whole development of Chinese prickly ash peel, the anthocyanins content was increased continuously, while most flavonoid differential metabolites were increased and decreased in green to yellow stage and yellow to red stage, respectively. Analysis of the differentially accumulated metabolites over the ripening stages revealed that flavonoids were the main metabolites modulating Chinese prickly ash peel reddening, and we identified several groups of flavonoids, including flavanone, flavone, flavonol, flavonoid and anthocyanins, that differentially accumulated in peels over the ripening stages. Kaempferol, quercetin, peonidin and their glycoside derivatives were the most important flavonoid compounds in Chinese prickly ash peels. In the present study, we found that peonidin o-hexoside and peonidin 3-o-glucoside were all up-regulated over the fruit ripening stages and therefore may be considered as the key anthocyanins conferring the red pigmentation of 'Hancheng Dahongpao' peels.
Flavonoid biosynthesis was mainly controlled by structural gene, which directly encoded various enzymes related to flavonoids secondary metabolism biosynthesis, and also was regulated by transcription factors, which were involved in regulating the expression of structural genes. Flavonoid biosynthesis were closely related to the expression of PAL, CHS, CHI, FLS, F3H, F3'H and other genes in pomegranate [36], strawberry [21], tea [37], apple [38], grape [39], asparagus [40], citrus [41] and other horticultural plants. In kiwifruit, F3H played a regulatory role in flavonoid biosynthesis by binding to transcription factors [42]. In the study of flavonoid metabolism in Ginkgo biloba, CHS, FLS, CHI, F3'H and DFR gene family, different genes of the same gene family expressed differently under dark and light conditions [43]. In this study, transcriptome data analysis showed that 26 differentially expressed genes were directly involved in flavonoid biosynthesis in 'Hancheng Dahongpao' peels at different development stages, The differentially expressed genes could be divided into three categories (Cluster I-III). The genes in cluster I were involved in phenylpropanoid biosynthesis, including three differentially expressed genes (C4H, CHS, CHI). Cluster II was linked to the synthesis of flavonoids and flavonols, and was composed of F3H, F3'H and FLS genes. The genes in cluster III were related to flavanols and anthocyanin synthesis, and there were four differentially expressed genes, including DFR, LAR, LDOX / ANS and UFGT. Through KEGG analysis, the core genes involved in flavonoids synthesis were C4H ( In grape, UFGT was reported to regulate anthocyanin biosynthesis during reddening, similarly as in our study [44]. The flavonoid metabolism was also controlled by transcription factors such as MYB-ZM1 (ZM-10723), MYB-ZM3(ZM-213334), MYB-ZM5 (ZM-65290), MYB-ZM6 (ZM-166024) and MYB-ZM7(ZM-71403). On the one hand, the unigenes MYB-ZM5 (ZM-65290) was clustered with VvMYBA1, MYBA1 and MYBA2, which regulated all anthocyanin biosynthesis in the various tissues of different plant species [45], and might be the candidate genes responsible for anthocyanin accumulation in Chinese prickly ash peel. On the other hand, MYB-ZM1, MYB-ZM3, MYB-ZM6 and MYB-ZM7 were closely gathered with flavonoid MYB transcription factors, such as FaMYB1 and AtMYB4, and regulated flavonoid accumulation via the transcriptional regulation of structural genes [33]. Combined metabolome and transcriptome analysis revealed the variations of flavonoids components and related regulatory genes in 'Hancheng Dahongpao' peels. A detailed pathway diagram of flavonoids biosynthesis of 'Hancheng Dahongpao' peels was constructed to better understand the relationship between genes and metabolites, and the flavonoids compounds formation. In addition, we would focus on how some changes in those gene networks regulated the biosynthesis of other nutrients and fruit economic traits in future studies.

Conclusion
In this study, the fresh fruits of Chinese prickly ash ('Hancheng Dahongpao') were used as the materials to reveal the flavonoids compounds formation mechanism of 'Hancheng Dahongpao' fruits based on the variation pattern of flavonoids metabolites at different developmental stages. 19 differential metabolites were identified in peels. 23 differentially expressed genes identified in the 5 genes modules were highly correlated with flavonoids at different developmental stages. Through the combined analysis of differential expressed genes and differential metabolites, 15 key candidate genes involved in flavonoid biosynthesis were finally screened. Altogether, this work lay an excellent foundation for the molecular mechanism of flavonoids biosynthesis and accumulation in 'Hancheng Dahongpao' peels at different development stages and provided a series of candidate genes with applications in the breeding of flavonoids-rich cultivars.