Combined Analysis of the Fruit Metabolome and Transcriptome Reveals Candidate Genes Involved in Flavonoid Biosynthesis in Actinidia arguta

1Kiwifruit (Actinidiaceae, genus Actinidia), a kind of perennial and deciduous plant, is one of the fruit trees that have been domesticated and cultivated successfully in the last century, and in recent years has been grown commercially worldwide. As a fruit tree originating from China, the genus Actinidia comprises a total of 76 species and about 125 known taxa worldwide, but among these, only two species—Actinidia chinensis and A. deliciosa—have been cultivated commercially.

2Flesh color is an important quality characteristic and is also an important index for evaluation of fruit quality in Actinida. The two species mentioned above, A. chinensis and A. deliciosa, have dominated the international market because of their large fruit, high vitamin C content, unique flavor, and long storage period properties. In addition, these two species are often used as experimental materials for the study of flesh coloring, which mainly concentrate on green- and yellow-fleshed coloring, corresponding to chlorophyll metabolism and carotenoid metabolism, respectively. In recent years, because of their rich functional components (phenols and flavonoids), red-fleshed kiwifruit has attracted more and more attention from consumers and researchers.

3Flavonoids, secondary metabolites found in plants, contribute to plant environmental adaptation, fruit development, and even human health. These compounds accumulate in various tissues at different developmental stages, and they are also the direct factors that cause red color in Actinida flesh. Flavonoid biosynthesis controlled by a diverse array of exogenous and endogenous factors is the most studied secondary metabolic pathway in plants, including Actinidia. Most structural gene-coding enzymes have been studied in Actinida, as well as transcription factors. However, most of the reports on red-fleshed Actinidia coloring are mainly concentrated on Actinidia chinensis, and rarely examine Actinidia arguta, which is a kind of all-red-fleshed mini-type Actinidia.

4The post-genomic era is densely populated with databases reporting the results of “-omics” studies, which provide insights into genome-wide approaches on the basis of comprehensive analysis of transcriptome, proteome, and metabolome, in addition to the traditional molecular biological, biochemical and reverse-genetic approaches. Accordingly, metabolomics is the attempt to identify and quantify all endogenous small molecule metabolites in an organism or biofluid sample. Furthermore, metabolomics analysis can be regarded as a technical means of association analysis together with other data to analyze the gene function involved in the metabolic pathway of interest and can also provide supporting information for gene mining. In recent years, combined high-throughput methods have been used to study color development. Combined metabolome and transcriptome analyses in fruits have clarified the relationship between the contents of various secondary metabolites and the corresponding differentially expressed genes, broadening the global view of plant color regulation.

5To investigate the metabolite changing rule associated with color change during Actinidia development, two different cultivars—“Hongbaoshixing” (“HB”) and “Yongfengyihao” (“YF”)—were selected to perform metabolome and transcriptome analysis. All metabolites that could be detected by instrument were identified and then used for differential analysis. Finally, 41 metabolites related to flavonoid biosynthesis were obtained, 7 of which presented significant difference in the comparison between HB100 and HB70. Through combined metobolome and transcriptome analysis, 2 pathways—flavonoid biosynthesis and galactose metabolism—were significantly different in terms of change during fruit development. In flavonoid biosyntheis, 8 differentially expressed structural genes were obtained. In addition, 7 transcription factor genes related to flavonoid biosynthesis were obtained by transcriptome analysis. RT-qPCR results and cluster analysis revealed the genes AaF3H, AaLDOX, AaUFGT, AaMYB, AabHLH and AaHB2 could be determined as candidate genes on which to conduct further study. Through integrated analysis of identified metabolites with significant difference and related genes involved in flavonoid biosynthesis, the association between metabolites and genes was explained so as to provide valuable information for revealing the red coloring of A. arguta.

6In order to explore the change of metabolites during A. arguta fruit development (Figure 1), a metabolome program was carried out. Only the m/z and chromatographic retention time were collected for each detected feature. Therefore, the quality control for features obtained by XCMS software was conducted by the overall MS signal intensity controlled by TIC (total ion chromatogram) (Figure S1), evaluation of width for m/z peak, and assessment of width for retention time peak for each feature (Figure S2).

7The original wiff format data was converted to the mzXML format using MSConvert software, followed by alignment and extraction of the peak and calculation of the peak area. The final statistics showed that 18,598 and 10,239 metabolites were obtained by the POS (positive) and NEG (negative) models, respectively, 9016 and 4699 of which were annotated (Table 1).

8The metabolites identified above were assigned to the KEGG and PLANTCYC databases. All 8070 and 4321 metabolites were classified into 18 KEGG second-grade pathways, 3580 and 2461 of which were classified into “metabolism” in the POS and NEG models, respectively. For the “metabolism” term, the top priority was “Global and overview maps”, followed by “biosynthesis of other secondary metabolites”, ”metabolism of terpenoids and polyketides”, “amino acid metabolism”, “metabolism of cofactors and vitamins” and “carbohydrate metabolism” (Figure 2). 480 and 340 metabolites were assigned to “biosynthesis of other secondary metabolites” in the POS and NEG models, respectively. Of these, 73 and 67 metabolites were involved in “flavonoid biosynthesis” and “anthocyanin biosynthesis” in the POS and NEG models, respectively (Figure 3).

9All 1267 and 360 metabolites obtained from level-two identification were assigned to the HMDB database, 550 and 250 of which were matched and classified into 11 HMDB super classes in POS and NEG, respectively. Of these, 305 metabolites were included in the “organic acids and derivatives” term, which presented the majority of all classes in the POS model. Meanwhile, “organic acids and derivatives”, which contained 59 metabolites, was the first class, followed by “organooxygen compounds” and “lipids and lipid-like molecules” in the NEG model (Figure 4).

10In order to obtain quantitative information on the metabolites, 14,132 and 8400 high-quality metabolites obtained from the 18,598 and 10,239 of all metabolites were used for differential analysis (Table 2). According to the quantitative results of the identified metabolites, the differential metabolites between different group comparisons were analyzed based on fold-change and p-value. In the comparison between HB100 and HB70, a total of 2421 and 1838 metabolites presented as being up-regulated and down-regulated, respectively (Table 3). In addition, a total of 41 metabolites were identified as being flavonoids among the high-quality metabolites, 7 of which were with significant difference, and were identified as bracteatin, luteolin, dihydromyricetin, cyanidin, pelargonidin, delphinidin and (−)-epigallocatechin (Table S1 and Table 4).

11A total of 2614 transcription factor genes were identified through transcriptome data analysis and could be grouped into 80 transcription factor families (Table S2), the majority of which were “AP2-EREBP”, “MYB” and “C2H2”, followed by “C3H“, ”Orphans“, “NAC”, “bHLH”, “WRKY” and “HB”. Of these 2614 transcription factor genes, 7 were annotated as being associated with flavonoid biosynthesis. The “c88340_g1”, “c105731_g1 (AaMYBC1)” and “c18002_g1 (named AaMYB)”, “c102250_g1” and “c86342_g1 (named AabHLH)”, and “c127253_g1 (named AaHB1)” and “c129627_g1 (named AaHB2)” belonged to “MYB”, ”bHLH“ and “HB” families, respectively (Figure 5).

12In order to investigate the association between metabolites and genes involved in the same biological process (KEGG Pathway), the comprehensive analysis of metabolome and transcriptome was performed using Pearson’s Correlation Coefficient. The results showed that 646, 156 and 484 differentially expressed genes participated in 13, 18 and 12 pathways among three groups HB100 vs. HB70, HB100 vs. YF10 and YF100 vs. YF70, respectively (Table S3). The 2 pathways, “flavonoid biosynthesis” and “galactose metabolism”, were a common presence among the three group comparisons. In flavonoid biosynthesis, a total of 8 candidate DEGs were found (Figure 6).

13A total of 14 structural genes comprising the 8 DEGs found above and 6 other structural genes involved in flavonoid biosynthesis were determined to conduct RT-qPCR, as well as 5 regulatory genes AaMYBC1, AaMYB, AabHLH, AaHB1 and AaHB2 (Figure 7). The expression patterns of these three structural genes AaF3H, AaLDOX and AaUFGT were similar to the two regulatory genes AaMYB and AaHB2. The expression levels of AaF3H, AaLDOX and AaUFGT in “HB100” was significantly higher than those in “HB70” and “YF100” (Figure 7A). A similar rule of expression was detected for AaMYB and AaHB2. The expression levels of AaMYB and AaHB2 in “HB100” were higher than those found in “HB70” and “YF100” (Figure 7B). Cluster analysis results showed that the genes AaF3H, AaLDOX, AaUFGT, AaMYB, AabHLH and AaHB2 were clustered into the same class (Figure 8), suggesting that the expression patterns of these genes were similar to each other. Because MYB is the biggest and most important transcription factor family involved in flavonoids, AaMYB was used for BLAST in order to search for other MYB sequences in other species. A phylogenetic tree was constructed with the BLAST results. The results showed that AaMYB and CsMYB5a (MYB5a of Camellia sinensis) were clustered into one class (Figure 9).

14In order to better understand the relationship between metabolites and genes in flavonoid biosynthesis, all results of metabolites and genes were combined to establish a network, aiming to show the relationship between gene expression and metabolite accumulation more intuitively (Figure 10). Among the 14 structural genes, AaF3H, AaLDOX and AaUFGT were highly expressed in HB100 vs. HB70 in comparison with the other 11 structural genes. The 7 metabolites with significant differences in HB100 vs. HB70 of flavonoid biosynthesis were: bracteatin, luteolin, dihydromyricetin, cyanidin, pelargonidin, delphinidin and (−)-epigallocatechin. It is generally known that transcription factor plays a role by binding the promoter of structural genes. Therefore, we speculate the hypothesis that the three transcription factors AaMYB, AabHLH, and AaHB2 activate expression of AaF3H, AaLDOX and AaUFGT by binding the promoter of these three structural genes, inducing the accumulation of the 7 metabolites described above in flavonoid biosynthesis. Specifically, AaMYB, AabHLH, AaHB2, AaF3H, AaLDOX and AaUFGT were the candidate genes obtained from this study and could be used for further study, revealing the red mechanism of A. arguta.

15Plant metabolomics, a new field in the post-genome era, presents the rule of change of metabolites in various tissues. Metabolites are the final products of cell biological regulation process, and their level can be regarded as the response of plant development to genetic and environmental changes. Therefore, metabolomic analysis can be used to investigate the relationship between biological processes and phenotypes; furthermore, some intuitive changes can also be observed at the metabolic level. Up until now, several reports have examined the metabolic responses of different species to flavonoid biosynthesis, such as Fagopyrum esculentum, Camellia sinensis, and Ficus carica L.. This research suggests that metabolome analysis plays a crucial role in explaining the molecular responses to flavonoid biosynthesis. In this study, through fruit metabolome analysis on two different A. arguta cultivars “HB” and “YF” at two developmental stages, a total of 28,837 metabolites were obtained, of which 13,715 had been annotated (Table 1). In addition, 22,532 high-quality metabolites were identified and selected for differential analysis (Table 2). The focus of our research was on flesh coloring, so the HB100 vs. HB70 comparison was selected as the object for further analysis. In the comparison HB100 vs. HB70, 41 metabolites were identified as being flavonoids, 7 of which, with significant difference, were identified as bracteatin, luteolin, dihydromyricetin, cyanidin, pelargonidin, delphinidin and (−)-epigallocatechin (Table S1). Such a result provides us with the insight that the difference between these 7 metabolites leads to the difference between red and green coloring in A. arguta flesh.

16Transcription factors are absolutely necessary for regulation of gene expression. Normally, transcription factor proteins play a role through the combination of their own DNA-binding domain and the cis-acting element of their target genes. In plants, transcription factors could participate in various biological process, including developmental regulation, defense elicitation, and stress responses. Numerous studies have shown that flavonoid biosynthesis is regulated by the MBW complex, comprising MYB, bHLH, WD40. Therefore, finding transcription factors related to flavonoid biosynthesis is essential for investigating fruit coloring. In this study, 5 transcription factor genes were obtained by transcriptome analysis and used for expression analysis. The results showed that AaMYB, AabHLH and AaHB2 were highly expressed at HB100, indicating that these three transcription factors might play a key role in fruit coloring in A. arguta. Phylogenetic analysis revealed that AaMYB was highly homologous to CsMYB5a, rather than the other MYB in A. chinensis, in which AcMYB75, AcMYB110 and AcMYBF110 were the key MYB transcription factors controlling fruit coloring, indicating that the MYB transcription factor that plays a key role in regulating flavonoid biosynthesis might be different in different Actinidia species.

17Genes, including structural and regulatory genes, involved in flavonoid biosynthesis and regulation have been found, studied and reported in many plants, including Actinidia. However, most genes have been obtained and identified using traditional study technology. Since transcriptome analysis has come to be regarded as a crucial way to study the expression level, structure and function of genes in order to revealing phenotypic traits, combined analysis of transcriptome and metabolome has increasingly become a popular and practical tool for the mining of new genes involved in various metabolic pathways. In this study, metabolome data, coupled with transcriptome profiling, was carried out in a combined analysis for the discovery of genes involved in flavonoid biosynthesis, thus searching for useful information to illustrate phenomenon of red coloring in A. arguta fruit. 19 genes, comprising 14 structural genes and 5 transcription factor genes, were obtained and used to analyze expression level; additionally, cluster analysis was conducted and the construction of phylogenetic tree was also performed. The results showed that structural genes including AaF3H, AaLDOX and AaUFGT and transcription factor genes including AaMYB, AabHLH and AaHB2 were highly expressed at HB100, when the flesh color significantly presented red (Figure 7). In addition, the cluster analysis results suggested that AaF3H, AaLDOX, AaUFGT, AaMYB, AabHLH and AaHB2 were clustered into one class, indicating that their expression patterns were similar to each other (Figure 8). Based on these results, a regulatory network of flavonoid biosynthesis was established to show the role of genes involved in pathways more intuitively (Figure 10). Thus, such a model of action was derived: the three transcription factors AaMYB, AabHLH and AaHB2 interact with promotors of the three structural genes aF3H, AaLDOX, and AaUFGT to control the expression, inducing the accumulation of metabolites and the appearance of the red phenotype of A. arguta fruit.

18This method—combined metabolome and transcriptome—is an effective analytical method for explaining the relationship between key genes and metabolites involved in biosynthesis pathways. Using this method, we determined the candidate genes and metabolites involved in the flavonoid biosynthesis pathway, providing valuable information and a useful reference for explaining the phenomenon of red coloring of A. arguta fruit. Nevertheless, the specific mechanism still needs further research, and remains to be explored.

19Two different types of kiwifruit (Actinidia arguta), “Hongbaoshixing” (“HB”, a kind of all-red-fleshed A. arguta cultivar) and “Yongfengyihao” (“YF”, a kind of all-green-fleshed A. arguta cultivar), were selected as experimental materials. Two different sampling stages of fruit development, defined in days after full bloom (DAFB), were 70 DAFB (70d, green fruit) and 100 DAFB (100d, the color-break stage, at which point they start to turn red). ‘YF’ fruit selected as a control case were green throughout the whole development process and were also sampled at 70d to 100d (Figure 1). Both of them are tetraploids with similar genetic backgrounds. The fruit samples at the two developmental stages described above were harvested from the National Kiwifruit Germplasm Garden of Zhengzhou Fruit Research Institute, CAAS, Henan province, China. Three biological replicates were collected per sample, each with 30 fruits randomly collected from 6 kiwifruit trees, every two of which were set as a biological replication; thus, all data were obtained based on three independent biological replicates. All flesh samples of the fresh fruits were dissected using a blade, frozen immediately in liquid nitrogen, and then stored at −80 °C until further use.

20The collected samples were thawed on ice, and metabolites were extracted with 50% methanol buffer (50% solution of methanol in distilled water). Briefly, 20 μL of sample was extracted with 120 μL of precooled 50% methanol, vortexed for 1 min, and incubated at room temperature for 10 min; the extraction mixture was then stored overnight at −20 °C. After centrifugation at 4000× g for 20 min, the supernatants were transferred into new 96-well plates. Three independent repetitions were executed for the extraction and subsequent analysis process.

21All chromatographic separations were performed using an ultra-performance liquid chromatography (UPLC) system (SCIEX, Cheshire, UK). An ACQUITY UPLC BEH Amide column (100 mm × 2.1 mm, 1.7 µm, Waters, UK) was used for the reversed phase separation. The column oven was maintained at 35 °C. The flow rate was 0.4 mL/min and the mobile phase consisted of solvent A (25 mM ammonium acetate + 25 mM NH4H2O) and solvent B (IPA:CAN = 9:1 + 0.1% formic acid). Gradient elution conditions were set as follows: 0~0.5 min, 95% B; 0.5~9.5 min, 95% to 65% B; 9.5~10.5 min, 65~40% B; 10.5~12 min, 40% B; 12~12.2 min, 40~95% B; 12.2~15 min, 95% B. The injection volume for each sample was 4 µL.

22A high-resolution tandem mass spectrometer TripleTOF5600plus (SCIEX, Cheshire, UK) was used to detect metabolites eluted form the column. The Q-TOF was operated in both positive and negative ion modes. XCMS software 3.2.0 (UC, Berkeley, CA, USA) was used to control the chromatograph and mass spectrometer.

23MSConvert was used to transform LC-MS raw data into the mzXML format, which was then processed by the XCMS, CAMERA and metaX toolbox, implemented in the R software. The combined retention time (RT) and m/z data were used to identify each ion.

24In order to explain the physical and chemical properties and biological functions of metabolites, the online PLANTCYC (http://www.plantcyc.org/), KEGG (http://www.kegg.jp/) and in-house (http://spldatabase.saskatoonlibrary.ca/), HMDB (http://www.hmdb.ca/) databases were used to perform level-one and level-two identification and annotation. Screening and quantitative analysis for differential metabolites were conducted using metaX software (http://metax.genomics.cn/).

25RNA isolation, purification and monitoring, and cDNA library construction and sequencing were performed as previously described. Briefly, RNA purity, concentration and integrity were checked, measured and assessed using the NanoPhotometer® spectrophotometer (IMPLEN, Westlake Village, CA, USA), Qubit® RNA Assay Kit in Qubit® 2.0 Flurometer (Life Technologies, Carlsbad, CA, USA) and RNA Nano 6000 Assay Kit of the Agilent Bioanalyzer 2100 system (Agilent Technologies, Santa Clara, CA, USA), respectively. Sequencing libraries were generated using NEBNext® Ultra™ RNA Library Prep Kit for Illumina® (NEB, Ipswich, MA, USA) following manufacturer’s recommendations and were sequenced on an Illumina Hiseq platform. According to the manufacturer’s instructions, NEBNext® Ultra™ RNA Library Prep Kit for Illumina® (NEB, Ipswich, MA, USA) was used to generate sequencing libraries, which were then sequenced on an Illumina Hiseq platform.

26Transcription factors were predicted using iTAK software. The identification and classification of transcription factors were conducted by previously described methods. On the basic of model of negative binomial distribution, DESeq R package (1. 10. 1) was used to perform differential expression analysis by providing statistical routines used to determine differential expression by means of the digital data of gene expression. The false discovery rate was deliberately controlled by the adjusted P values, which were adjusted by Benjamini and Hochberg’s approach. The differentially expressed genes were defined as whose adjusted p-value < 0.05 found by DESeq.

27Total RNA was extracted using a modified CTAB (cetyltrimethyl ammonium bromide) method. The testing of RNA quality and determination of RNA concentration were performed by 1.0% agarose gel electrophoresis and micro ultraviolet spectrophotometry (Thermo NanoDrop 2000, Thermo Fisher Scientific, Waltham, MA, USA), respectively. Approximately 1µg of total RNA was determined for cDNA synthesis using RevertAid™ First Strand cDNA synthesis kit (Thermo Fisher Scientific, Waltham, MA, USA). The LightCycler® 480 real-time PCR system with a 96-well plate was used to conduct an amplified reaction consisting of 95 °C for 5 min, followed by 45 cycles of 10 s at 95 °C, 20 s at 60 °C, and 20 s at 72 °C in a volume of 10 μL. At the end of each experiment, a melt-curve analysis was carried out using the default parameters (5 s at 95 °C and 1 min at 65 °C). The Actinidia β-actin was used for normalization. All analyses were repeated three times using biological replicates. Primer sequences are listed in Table S4.

28The relative expressions were calculated using the 2−ΔΔCt method, and GraphPad Prism 5 (GraphPad Software Inc., San Diego, CA, USA) was used for chart preparation. The R-3.4.2 and MEGA6 were used to conduct the heatmap and cluster analysis. IBM SPSS Statistics 20 was used to test significant differences.