Metabolomics Study of Flavonoids and Anthocyanin-Related Gene Analysis in Kiwifruit (Actinidia chinensis) and Kiwiberry (Actinidia arguta)

This study investigated the flavonoid compounds in Actinidia chinensis and Actinidia arguta fruits. A total of 125 flavonoids, including 9 anthocyanins, 12 catechins, 17 flavanones, 48 flavones (including 14 flavone C-glycosides), 29 flavonols, 6 isoflavones, and 4 proanthocyanidins, were identified in “Hongyang” kiwifruit (red flesh), “Jintao” kiwifruit, “Mini Amethyst” kiwiberry (purple flesh), and “Kuilv” kiwiberry. Thirty-nine metabolites showed significantly different contents between “Hongyang” and “Jintao,” and 38 of them showed higher content in “Hongyang,” whereas 39 metabolites showed significantly different contents between “Mini Amethyst” and “Kuilv,” and 31 of them showed higher content in “Mini Amethyst.” This result indicates the superior nutritional value of the pigmented kiwi cultivars in terms of flavonoids. Multivariate statistical analysis indicates that the variation in flavonoid profiles contributes to the pigmentation phenotypes of “Hongyang” and “Mini Amethyst.” Further comparative transcriptomic analysis revealed that structural genes in the anthocyanin synthesis pathway (AcF3H, AcF3′H, AcDFR, AcUFGT) and transcription factors (AcMYB10, AcbHLH5) may be involved in the pigmentation of the red-fleshed A. chinensis, whereas AaF3H, AaF3GT, and AaMYB110 may play important roles in the pigmentation of the purple-fleshed A. arguta. This study provides broader insight into the variation in flavonoid profiles among kiwifruit/berry, evaluates the flavonoid nutrition of the four cultivars, and provides additional evidence for the correlation between the genes and metabolites involved in flavonoid synthesis.


Introduction
Kiwifruit (Actinidia chinensis) is a high-value fruit crop that is recognized for its nutrient content and attractive appearance (Li et al. 2007). Kiwiberry (Actinidia arguta), also called "hardy kiwifruit," a semiwild kiwifruit relative, has attracted increasing attention from breeders and customers due to its smaller fruit, glabrous peel, and distinct flavor Latocha and Jankowski 2011). Most kiwifruit and kiwiberry varieties have green or pale-yellow flesh and are rich in vitamin C, dietary fiber, and flavonoids that can contribute to human diet (Richardson et al. 2018). Meanwhile, few colored varieties, such as the red-pigmented A. chinensis "Hongyang" and purple-pigmented A. arguta var. arguta, have exhibited great commercial potentials in recent years due to their gorgeous color and antioxidant function provided by accumulated anthocyanins (Jaeger and Harker 2005). These qualities all share at least a partial basis in the flavonoid pathway, which gives rise to fruit colors, flavors, and healthful compounds (Winkel-Shirley 2001). Thus, previous works have highlighted anthocyanin metabolites and pigmentation mechanisms of pigmented kiwifruit and kiwiberry (Montefiori et al. 2005;Montefiori et al. 2009). However, large-scale evaluation of anthocyanin-related flavonoid compounds has seldom been reported despite colorless/pale-color or yellow flavonoids predominating among the total flavonoids of kiwi species, influencing fruit pigmentation and supporting health functions as well.
Flavonoids represent one of the largest classes of plant secondary metabolites and are involved in performing various physiological functions, including UV protection, insect attraction, and pathogen defense as well as leading to variation in plant color (Winkel-Shirley 2001). Flavonoids not only contribute to the flavor of fruit crops but also act as a powerful antioxidant in the human diet (Crozier et al. 2009). Flavonoids are synthesized via the flavonoid branch of the phenylpropanoid metabolic pathway. Based on the substituents, flavonoids can be subclassified into anthocyanidins, flavanones, flavones, isoflavones, flavanols, and flavonols (Tohge et al. 2017), which are derived through the flavonoid biosynthetic pathway using a set of common enzymatic reactions. Anthocyanin is a kind of end-product of the flavonoid pathway and is synthesized by a series of enzymes, including chalcone synthase (CHS), chalcone isomerase (CHI), flavanone 3-hydroxylase (F3H), flavonoid 3′-hydroxylase (F3′H), flavonoid 3′5′-hydroxylase (F3′5′H), dihydroflavonol 4reductase (DFR), leucoanthocyanidin dioxygenase (LDOX), and UDP-glucose flavonoid 3-O-glucosyltransferase (UFGT). Since anthocyanin is the most important water-soluble pigment in plants (He and Giusti 2010), members of the anthocyanin chemical family, including pelargonidin, cyanidin, delphinidin, peonidin, petunidin, and malvidin, have been well characterized in many important fruit crops (Jaakola 2013). However, according to anthocyanins biosynthesized through the specific branched pathway of the flavonoid pathway (Zhang et al. 2014), the metabolism of anthocyanin was also dynamically correlated with several other flavonoids located upstream or on other branches of the metabolic pathways (Nabavi et al. 2018). Thus, only highlighting anthocyanidins or a few anthocyanin-related flavonoids, such as dihydromyricetin, bracteatin, luteolin, and (−)-epigallocatechin (Li et al. 2018b), is limiting to reveal the dynamic metabolic pattern underlying the pigmentation of kiwifruit and kiwiberry and is also not sufficient to pave the way for further gene function study and genetic improvement of these taxa.
The genetic mechanisms underlying anthocyanin accumulation have been repeatedly studied in many horticultural crops via transcriptome analysis (El-Sharkawy et al. 2015;Cao et al. 2018;Lin et al. 2018;Bai et al. 2017). Although transcriptomic studies were conducted on kiwifruit and kiwiberry, and some candidate genes were screened out (Li et al. 2018a;Li et al. 2018b;Li et al. 2015b;Wang et al. 2019;Peng et al. 2019), these results were mainly based on limited numbers of anthocyanidin-related metabolites. Furthermore, pigmented kiwifruit and kiwiberry also showed distinct pigmentation patterns. For example, the pigmentation of red kiwifruit always centered around the fruit core and showed a bright-red color, whereas the pigmentation of purple kiwiberry showed a dark-red or purple color with a global pattern. However, due to the different genetic backgrounds of these two species, the pigmentation mechanisms of these two species have not been discussed together despite the previous studies having selected analogous candidate genes (Li et al. 2018a;Li et al. 2018b;Li et al. 2015b). Currently, a combination of metabolic and transcriptomic profiles has provided important information for understanding gene-to-metabolite relations, which have been successfully applied in many plants, such as potato (Cho et al. 2016;Stushnoff et al. 2010), tea (Li et al. 2015a), and fig . Hence, more comprehensive flavonoid identification together with transcriptomic analysis would not only deepen our understanding of kiwifruit and kiwiberry pigmentation mechanisms but also lay the foundation for the evolutionary study of anthocyanidin biosynthesis in the Actinidia genus.
In the present study, a large-scale flavonoid characterization was carried out on red-fleshed A. chinensis "Hongyang" (AcH), A. chinensis "Jintao" (AcG), purple A. arguta "Mini Amethyst" (AaP), and A. arguta "Kuilv" (AaG) (Fig. 1a). Then, multiple statistical analyses were conducted to address the differences in flavonoid accumulation between pigmented and nonpigmented individuals and between kiwifruit and kiwiberry. Finally, the transcriptomic differences of pigmented and nonpigmented kiwifruit and kiwiberry were investigated respectively according to the differential flavonoid metabolites.

The Flavonoid Metabolic Profiling of Four Kiwi Cultivars
To profile the flavonoid metabolites and evaluate the nutritional values of the four kiwi cultivars with different flesh colors, a newly developed LC-MS-based, widely targeted metabolomics method was employed to identify and relatively quantify the flavonoid metabolites.
Intriguingly, all four cultivars shared four anthocyanidins (cyanidin 3-O-glucoside, malvidin 3-O-galactoside, malvidin 3-O-glucoside, and cyanidin 3-O-rutinoside) ( Fig. 1b; Table S1), indicating that the unpigmented kiwifruit/berry can also synthesize anthocyanins and suggesting that their accumulation may be limited, which is in accordance with a previous literature (Li et al. 2018b). Three anthocyanins (delphinidin 3-O-glucoside, cyanidin O-syringic acid, and pelargonidin 3-Obeta-D-glucoside) were only detected in red-fleshed AcH and purple-fleshed AaP ( Fig. 1b; Table S1), suggesting a key role in the fruit pigmentation of these kiwi species. This result is in agreement with previous reports that cyanidin derivatives have been found in all Actinidia species that accumulate anthocyanins (Montefiori et al. 2009). However, delphinidin 3-Oglucoside was detected both in AcH and AaP ( Fig. 1b; Table S1), which differs from the results of a previous study that reported that delphinidin derivatives only existed in A. melanandra and A. arguta var. purpurea (Montefiori et al. 2009). In addition, pelargonidin 3-O-beta-D-glucoside was detected in AcH and AaP (Table S1), which have not been reported in the genus Actinidia. In addition, cyanidin Oacetylhexoside was only detected in AaP (Table S1), suggesting its contribution to the red color of AaP. Collectively, the above results not only indicate that kiwifruit and kiwiberry are abundant in flavonoid components but also suggest a wider range of nutritional flavonoid-derived than anthocyanin-derived nutritional compounds.
To further assess the differences in metabolic profiles among the four cultivars, hierarchical cluster analysis (HCA) was performed on the above profiles. As shown in Fig. 2a, the metabolite profiles fall into four main clusters: clusters I, II, III, and IV represent the highly accumulated flavonoids in AaG, AaP, AcG, and AcH, respectively. In addition, principal component analysis (PCA) was performed on the relative content of flavonoids in the four cultivars to address the internal structure of flavonoid variation. As shown in Fig. 2b, quality control (QC) samples ("mix01," "mix02," and "mix03") grouped closely together, indicating the integrity and accuracy of the following analysis. The two main PCs (PC1 and PC2) accounted for 37.13% and 24.75% of the total system variability, respectively, and the four cultivars grouped into three main groups. AcH (red fleshed) and AaP (purple fleshed) were distinguished from the other three, whereas AcG and AaG clustered fairly close together. These results suggest that the nature and the variation in flavonoids in both species may contribute to the pigmentation of AcH and AaP. However, the variation in flavonoid profiles did not clearly separate between AcG and AaG but did clearly separate between AcH and AaP, suggesting that distinct genetic mechanisms underlie the metabolite profiles of AcH and AaP.
To further assess the differential flavonoid metabolites that induced variation in metabolite profiles, partial least squares discriminant analysis (PLS-DA) was performed on the flavonoid profiles between AcH and AcG, AaP and AaG, AcH and AaP, and AcG and AaG. As shown in Fig. 3a, high predictability (Q2) and strong goodness of fit (R2X, R2Y) of the PLS-DA models were observed in the comparisons between AcH and AcG (Q2 = 0.994, R2X = 0.845, R2Y = 0.999), between AaG and AaP (Q2 = 0.995, R2X = 0.9, R2Y = 0.999), between AcH and AaP (Q2 = 0.998, R2X = 0.944, R2Y = 1), and between AcG and AaG (Q2 = 0.986, R2X = 0.846, R2Y = 997). Screening of differential metabolites was performed based on fold change (FC) analyses, and variables determined to be important in the projection scores (VIP) were identified (VIP ≥ 1.0 and fold change ≥ 2 or ≤ 0.5 were set as thresholds). As shown in Fig. 3b, between AcH and AcG, a total of 39 differential metabolites (including 6 anthocyanins, 2 catechins, 10 flavanones, 8 flavones, 2 flavone Cglycosides, 9 flavonols, and 2 isoflavones) showed significantly different content, 38 of which showed higher content in AcH than that in AcG, and only syringetin showed higher content in AcG, suggesting a superior nutritional value of AcH in terms of flavonoids (Table 1). According to the previous reports, anthocyanins have antioxidant, angiocardiopathy, and cancer prevention capacities and are recognized as natural health ingredients with beneficial effects on humans (Pojer et al. 2013). In our study, anthocyanin content varied greatly between AcH and AcG. For example, delphinidin 3-O-glucoside, cyanidin O-syringic acid, and pelargonidin 3-O-beta-Dglucoside were only detected in AcH, and the contents of cyanidin 3-O-glucoside, cyanidin 3-O-rutinoside, and pelargonidin 3-O-beta-D-glucoside were 1753-, 9-, and 371fold higher in AcH than those in AcG (Table 1). These results agree with the finding that cyanidins are the predominant pigment in kiwifruit with red colors (Montefiori et al. 2009;Liu et al. 2018).
For kiwiberry, 39 differential metabolites were identified between AaP and AaG. Among them, 31 showed higher content in AaP (6 anthocyanins, 5 catechins, 1 flavanone, 6 flavones, 4 flavone C-glycosides, 8 flavonols, and 1 isoflavone), whereas 8 showed higher content in AaG (1 flavanone, 6 flavones, and 1 flavone C-glycosides) ( Table 2), suggesting a more complicated regulation of metabolic flux in kiwiberry compared with that in kiwifruit. Among the six differentially accumulated anthocyanins, the content of cyanidin 3-Oglucoside in AaP was 1358-fold higher that that present in AaG, and the other five anthocyanins (delphinidin 3-O-  The content of each metabolite comes from the average content of three biological replicates. b Score plots for principle components 1 and 2. "Mix" corresponds to the quality control sample. The same color represents three biological replications  (Table 2). Collectively, the content level of flavonoids was higher in pigmented kiwifruit/berry than that in unpigmented kiwifruit/berry either in kind or in quantity, which suggests that these cultivars may differ in the expression of flavonoid biosynthetic or regulatory gene expression. Between AaP and AcH, 43 differential metabolites were identified. In addition, 21 of them (2 catechins, 3 flavanones, 7 flavones, 3 flavone C-glycosides, 6 flavonols, and 1 isoflavone) showed higher content in AaP, whereas 22 (1 anthocyanin, 6 flavanones, 8 flavones, 1 flavone C-glycoside, 2 flavonols, and 3 isoflavones) showed higher content in AcH (Table S1). Between AcG and AaG, 37 metabolites (1 anthocyanin, 3 catechins, 5 flavanones, 12 flavones, 2 flavones C-glycosides, 11 flavonols, and 3 isoflavones) showed significantly differential content, among which 13 showed higher and 24 showed lower content in AcG (Table S1). These results indicate a great difference in flavonoid composition in pigmented kiwifruit/kiwiberry and unpigmented kiwifruit/kiwiberry, suggesting the distinct flavonoid synthesis and regulation mechanism of A. chinensis and A. arguta.
Among these DEGs involved in anthocyanin biosynthesis, F3H, F3′H, DFR, and UFGT showed significantly upregulated expression in AcH. Combined with the metabolome results, although cyanidins are the major anthocyanin category in redfleshed kiwifruit, only the contents of dihydrokaempferol and the final cyanins were significantly higher in AcH than those in AcG (Table S1). Since F3H catalyzes the conversion of naringenin to dihydrokaempferol and UFGT glycosylates anthocyanidins, these results suggest a regulatory scenario: the upregulated expression of F3H promotes the accumulation of dihydrokaempferol, providing sufficient substrates for the biosynthesis of downstream metabolites, and UFGT increases the hydrophilicity and stability of cyanidins, mediating the flux of flavonoid intermediates towards anthocyanin biosynthesis. Thus, we speculate that F3H and UFGT perform key functions in the pigmentation of red-fleshed kiwifruit. In addition, F3′5′H showed low expression levels in AcH and no expression in AcG ( Fig. 4; Table S2). To our knowledge of anthocyanin biosynthesis, dihydrokaempferol would be hydroxylated by F3′ H to produce dihydroquercetin leading to the accumulation of cyanidins, which are responsible for red color (Montefiori et al. 2009); alternatively, dihydrokaempferol could be hydroxylated by F3′5′H to produce dihydromyricetin leading to accumulation of delphinidins, which are responsible for the purple color (Dixon et al. 2013). Since delphinidin 3-O-glucoside was also detected in AcH, and F3′H showed a higher expression level than F3′5′H in AcH, we speculate that F3′H also plays an important role in cyanidin accumulation in AcH by competing with F3′5′H. In addition, DFR and FLS showed higher expression in AcH, and flavonols and anthocyanins also showed higher content in AcH (Table 1). Coupled with the fact that FLS and DFR compete for common substrates to direct biosynthesis to flavonols or anthocyanins, we speculate that the anthocyanin accumulation of AcH is also influenced by DFR.
In addition, it has been proven that transcription factors (TF), such as R2R3-MYB, basic helix-loop-helix (bHLH), and WD40, are involved in flavonoid biosynthesis through regulating the transcription of the above structural genes (Zhao et al. 2013;Lin-Wang et al. 2010). In our transcriptome dataset, 82 MYB TFs were identified in the DEGs (Table S2), and nine of the TFs belonged to the R2R3-MYB subfamily. Among the nine R2R3-MYB, five (Acc29865, Acc00493, Acc07258, Acc10925, and Acc02217) showed upregulated expression in AcH, and Acc00493 (AcMYB10) showed the greatest fold change (~32-fold) between AcH and AcG, suggesting that this gene regulates the activity of the above structural genes involved in anthocyanin biosynthesis (F3H, F3′H, DFR, and UFGTs). In addition, 47 bHLH TFs showed markedly different expression between AcH and AcG; 11 of these TFs were MYC-like bHLH TFs, and eight of them showed higher expression levels in AcH (Table S2). Intriguingly, one bHLH5 (Acc19563) was reported to be the partner protein of AcMYB10 and AcMYB110 (Brendolise et al. 2017), which play an important role in the regulation of anthocyanin synthesis. In addition, thirteen WD40 TFs were also differentially expressed between AcH and AcG (Table S2).
To further confirm the results of the transcriptomic analysis, the expression of the 16 structure genes, 4 MYBs, and 7 bHLHs was measured by qRT-PCR in the two kiwifruit cultivars. The qRT-PCR results showed that most of the differences in expression were consistent with the expression levels obtained by RNA-seq (Fig. 5). The qRT-PCR data for 5GT (Acc20256) and FLS (Acc24372), however, did not differ significantly between AcH and AcG. As expected, the expression levels of F3H (Acc01070), F3′H (Acc25698), DFR (Acc01005), and 3GT (Acc15351, Acc21776, Acc26546, and Acc26544) werẽ 5.2-,~28.8-,~3.8-fold, and > 15-fold higher in AcH than those in AcG (Fig. 5), respectively, and the expression levels of AcMYB10 (Acc00493) and bHLH5 (Acc19563) were~7.9-and~15.8-fold higher in AcH than those in AaG, respectively. The above results indicated the reliability of the transcriptomic analysis and suggested the involvement of these candidate genes in anthocyanin biosynthesis.

Transcriptome Analysis and Candidate Genes Involved in Anthocyanin Accumulation of Purple-Fleshed Kiwiberry
A total of 44.11 Gb clean data was obtained from the RNAseq libraries prepared from AaP and AaG, with 89.70-90.81% of Q30 rate (Table 4). Since the mapping rate of A. arguta reads was not satisfactorily high (~20% to the two A. chinensis kiwifruit genomes (Huang et al. 2013;Pilkington et al. 2018)), we de novo assembled A. arguta transcripts. A total of 129,537 unigenes were assembled with an N50 of 861. Then, all the assembled unigenes were annotated by public databases (Nr, Nt, Swiss-Prot, GO, and KEGG), and a total of 35,407 unigenes were annotated by using BLAST (E-value ≤ 1e-5) and HMMER (E-value ≤ 1e-10). A total of 7686 DEGs were detected, with 3721 upregulated and 3965 downregulated in AaP (Table S4). These DEGs were then subjected to GO analysis, and all the DEGs were classified into 54 functional groups, including 20 groups in "Biological processes," 17 groups in "Cellular components," and 17 in "Molecular functions." Among them, "Metabolic process," "Cell," and "Binding" were the most abundant terms under the above three categories (Fig. S1c). KEGG analysis revealed that "Metabolism" was also the most enriched category, followed by "Genetic information processing" (Fig. S1d). There were 29 DEGs enriched in "Phenylpropanoid biosynthesis (ko00940)" and 14 DEGs enriched in "Flavonoid metabolism (ko00941)." Almost all the structural genes in anthocyanin biosynthesis showed higher expression in AaP, such as one PAL, one C4H, two 4CLs, two CHSs, two CHIs, one F3H, three F3′H, one DFR, one LDOX (ANS), two 3GTs, two 5GTs,  one FLS, one COMT, one ANR, one LAR, and one HCT gene ( Fig. 4; Table S3).
In the flavonoid biosynthesis network, naringenin acts as a junction point: it can be either catalyzed by F3H to produce dihydrokaempferol leading to anthocyanin biosynthesis or catalyzed by flavone synthase (FNS) to produce flavone derivatives ( Fig. 4) (Zhang et al. 2014). It is worth noting that F3H (c96449.g0) showed~600-fold higher expression in AaP than that in AaG; the metabolites downstream of dihydrokaempferol, especially anthocyanidin derivatives, showed higher accumulation in AaP, whereas one flavanone and seven flavone derivatives showed higher accumulation in AaG ( Table 2). The above results indicate that more naringenins were allocated to the synthesis of flavone and flavanone in AaG, whereas more naringenins were allocated to synthesis of dihydrokaempferol and anthocyanins in AaP, suggesting that F3H (c96449.g0) can switch the metabolic flux towards anthocyanin biosynthesis in AaP. On the other hand, in our metabolite profile, delphinidin 3-O-glucoside showed the highest content in AaP among all four cultivars. Delphinidin has been reported to be responsible for purple color in the plant kingdom (Dixon et al. 2013). In our transcriptomic dataset, the expression of F3′5′H (c109457.g0) in AaP was higher than that in AaG, indicating its contribution to the unique purple color of AaP. In addition, the unigene (c77921.g0) encoding a protein that is highly homologous with the reported AcF3GT1 and AcF3GT1 has been proven to be a key enzyme that regulates the accumulation of anthocyanin in red-fleshed kiwifruit cultivars (Montefiori et al. 2011). This gene (c77921.g0) showed2 70-fold higher expression in AaP than that in AaG. Thus, we designate c77921.g0 as AaF3GT and speculate that it plays an important role in the anthocyanin biosynthesis of A. arguta. Furthermore, 54 MYB TFs were differentially expressed between AaP and AaG, and 32 of them were upregulated in AaP (Table S3), of which 10 belong to R2R3 MYB-like genes. Thirty-one bHLHs showed differential expression between AaP and AaG, and four of them (c115933.g0, c103501.g2, c110375.g0, and c115920.g0) were MYC-like bHLH TFs and showed higher expression in AaP (Table S3). Only one WD40 was upregulated in AaP (Table S3). It is worth noting that one of these, R2R3-MYB, AaMYB (c94087.g0), showed an~300-fold higher expression in AaP than that in AaG. AaMYB (c94087.g0) encodes a protein that is highly homologous to AcMYB110, which has been reported to regulate anthocyanin production in the petals of A. chinensis and A. eriantha (Fraser et al. 2013). Thus, we designated this gene as AaMYB110, and we speculate that it plays an analogous role in the anthocyanin accumulation of A. arguta. To confirm the results of the transcriptomic analysis, 25 candidate genes were selected from DEGs and were subjected to qRT-PCR. The expression trend was representative of all the selected genes (except DFR), and these findings were consistent with the transcriptomic analysis results (Fig. 6). As expected, AaF3H, AaF3GT, and AaMYB110 showed~148-, 308-, and~1418-fold higher expression in AaP than that in AaG, respectively.

Conclusion
We identified 125 flavonoid metabolites in A. chinensis "Hongyang" (AcH) and "Jintao" (AcG), purple A. arguta accession "Mini Amethyst" (AaP), and A. arguta "Kuilv" (AaG); this study represents the broadest characterization of the flavonoids in kiwi species to date. Combined with flavonoid metabolic profiles, transcriptome analysis screened out candidate genes, AcF3H, AcF3′H, AcDFR, AcF3GT, AcMYB10, and bHLH5, which may be responsible for the anthocyanin accumulation in "Hongyang," and AaF3H, AaF3GT, and AaMYB110, which may be responsible for anthocyanin accumulation in A. arguta "Mini Amethyst." These findings improved our understanding of the flavonoid metabolism and pigmentation mechanism in kiwi species and laid the foundation for related gene function studies in the future.

Plant Materials
The red-fleshed A. chinensis "Hongyang" (AcH) and "Jintao" (AcG), purple A. arguta accession "Mini Amethyst" (AaP), and A. arguta "Kuilv" (AaG) were collected from the germplasm field of Wuhan Botanical Garden. Three vines of each accession were treated as three biological replicates, and 15-20 fruits were sampled from each vine. All the fruits were sampled at a ripe-for-eating stage. The endocarp of A. chinensis (AcH and AcG) and the flesh of A. arguta (AaP and AaG) were sampled. The fleshes of sampled fruits were first cut into pieces, immediately frozen in liquid nitrogen, and then stored at − 80°C until RNA and metabolite extraction.

Metabolite Extraction and Identification
The freeze-dried samples were crushed to a uniform powder using a mixer mill (1.5 min, 30 Hz; MM 400, Retsch). A measure (100 mg) of powder from each sample was suspended in 1.0 mL 70% aqueous methanol, kept at 4°C for 12 h in the dark, and centrifuged at 10,000 rpm for 10 min at 4°C. Finally, the supernatant solution was absorbed (CNWBOND Carbon-GCB SPE Cartridge, 250 mg, 3 mL; ANPEL, Shanghai, China) and filtered (SCAA-104, 0.22-μm pore size; ANPEL) before LC-MS analysis. LC-MS analysis was performed on an LC-ESI-MS/MS system (HPLC, Shim-pack UFLC SHIMADZU CBM30A system; MS, Applied Biosystems 4500 Q TRAP). Five microliters of the extracted solution was injected into an HPLC system (Shimpack UFLC SHIMADZU CBM30A) equipped with a C18 column (Waters ACQUITY UPLC HSS T3, 1.8 μm, 2.1 mm × 100 mm). The binary solvent system was ultrapure water containing 0.04% acetic acid as the mobile phase A and acetonitrile containing 0.04% acetic acid as the mobile phase B. The A:B (v/v) gradient was 95:5 at 0 min, 5:95 at 11.0 min, 5:95 at 12.0 min, 95:5 at 12.1 min, and 95:5 at 15.0 min. The flow rate was kept at 0.40 mL/min, and the column temperature was maintained at 40°C.
The HPLC effluent was connected to an electrospray ionization (ESI)-triple quadrupole-linear ion trap-MS/MS system (Applied Biosystems 4500 Q TRAP). The ESI source operation parameters were as follows: ion source, turbo spray; source temperature, 550°C; ion spray voltage (IS), 5500 V; ion source gas I (GSI), gas II (GSII), and curtain gas (CUR) were set at 55, 60, and 25.0 psi, respectively; and collision gas (CAD) was high. Instrument tuning and mass calibration were performed with 10 and 100 μmol/L polypropylene glycol solutions in QQQ and LIT modes, respectively. The monitoring mode was set to multiple reaction monitoring (MRM). Analyst 1.6.1 was used to process the mass spectrometry data. Peak area integration was performed on all mass spectral peaks. The area of each chromatographic peak (Area) represented the relative amount of the corresponding substance. Integral calibration was performed for the peaks in the mass spectrum of the same metabolite in different samples.

Multivariate Statistical Analysis and Targeted Metabolite Identification
Qualitative analysis of primary and secondary MS data was carried out by comparing the accurate precursor ions (Q1), product ions (Q3), retention time (RT), and fragmentation patterns with those obtained by the injection of standards under the same conditions if the standards were available (Sigma-Aldrich, USA; http://www.sigmaaldrich.com/united-states.html), or by a selfcompiled database MWDB (MetWare biological science and Technology Co., Ltd. Wuhan, China) and publicly available metabolite databases if the standards were unavailable. PCA and orthogonal partial least-squares discriminant analysis OPLS-DA with data from 6 samples (two cultivars × three biological replicates) were performed to evaluate the differences in metabolite composition. The reliability correlation (p(corr)) values of all metabolites from the S-plot of the OPLS-DA were obtained using the first component. The metabolites satisfying the following criteria were selected as potential markers: fold change ≥ 2 or ≤ 0.5 and VIP ≥ 1 (variable importance in project (VIP) of the OPLS-DA model). The metabolite content data were normalized using the range method, and the clustering analysis of metabolites between different samples was performed with R software. Molecular formulas and assignments of the compounds to flavonoid pathways were established through KEGG (Kyoto Encyclopedia of Genes and Genomes) (http:// www.genome.jp/kegg).

RNA Extraction, Library Construction, and Illumina Sequencing
A 100 mg frozen sample was ground to a fine powder in liquid nitrogen. Then, the total RNA was extracted using TRIzol reagents (Invitrogen, USA). Each RNA sample was subjected to DNase digestion (TaKaRa, Dalian, China) to remove the remaining DNA. RNA degradation and contamination were monitored on 1% agarose gels. The purity was checked using a NanoPhotometer® spectrophotometer (IMPLEN, CA, USA). The concentration was measured using a Qubit® RNA Assay Kit with a Qubit®2.0 Fluorometer (Life Technologies, CA, USA), and its integrity was assessed using the RNA Nano 6000 Assay kit with the Agilent Bioanalyzer 2100 system (Agilent Technologies, CA, USA).
Three micrograms of total RNA from each sample was prepared for construction of a sequencing library. Sequencing libraries were generated using NEBNext®Ultra™ RNA Library Prep kit for Illumina® (NEB, USA), and index codes were added to attribute sequences of each sample. All procedures for cDNA library construction were performed following the standard Illumina sample preparation protocol. The library quality was assessed on the Agilent Bioanalyzer 2100 system. The library preparations were sequenced on an Illumina HiSeq 2000 platform, and the paired-end reads were generated.

De Novo Transcriptome Assembly and Functional Annotation
Adapter-containing, ploy-N, or low-quality reads of raw data were removed by in-house Perl scripts to generate clean data. Q20, Q30, and GC content and sequence duplication level were calculated based on the clean data. Then, transcriptome assembly was accomplished using Trinity. Gene function was annotated based on NR (NCBI nonredundant protein sequences), Pfam (protein family), KOG/COG/eggNOG (Clusters of Orthologous Groups of proteins), Swiss-Prot (a manually annotated and reviewed protein sequence database), and KEGG and GO (Gene Ontology) databases. The MEME 5.0.2 online program (http://meme-suite.org/) was used for motif identification of the annotated MYB protein sequences.

Screening and Enrichment of Differentially Expressed Genes
The gene expression levels of each sample were estimated by RSEM. Differential expression analysis of each of the two sample groups was performed using DESeq. Genes with an adjusted P value of 0.05 and abs (log2 (fold change)) of 1 found by DESeq were designated as differentially expressed genes (DEGs). Gene Ontology (GO) enrichment analysis of the DEGs was implemented by the topGO R-packages based on the Kolmogorov-Smirnov test. KOBAS software was used to test the statistical enrichment of differentially expressed genes in the KEGG pathways.
Validation of RNA-seq Results by Quantitative Real-time PCR Genes related to anthocyanin biosynthesis were selected from the DEG results. The primers used in RT-PCR are shown in Table S4. RNA was extracted from samples transformed using HiPure Plant RNA kits (Magen, R4151-03). The cDNA libraries were constructed by using HiScript II QRT SuperMix for qPCR (Vazyme, R223-01) and were then used as the input for quantitative PCR (qPCR) experiments. The RT-PCR was performed on a Roche (Indianapolis, IN, USA) LightCycler®480, following the manufacturer's instructions for LightCycler®480 SYBR Green I Master mix (Vazyme). The reactions were carried out with the following cycling program: 95°C for 30 s, followed by 45 cycles of 95°C for 10 s, 60°C for 10 s, and 72°C for 20 s. Melt-curve analyses were carried out as follows: 15 s at 95°C, 1 min at 60°C, 30 s at 95°C, and 15 s at 60°C. Each sample was amplified in triplicate, and Roche LightCycler 480 software version 1.5.1.62 was used to perform the data analysis. The 2 −ΔΔCt method was employed with Achn107181 (kiwifruit actin gene) as an endogenous control (Petriccione et al. 2015).
need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.