Comparative transcriptomic analysis of white and red Chinese bayberry (Myrica rubra) fruits reveals flavonoid biosynthesis regulation

Scientia Horticulturae 235 (2018) 9–20

Contents lists available at ScienceDirect

Scientia Horticulturae journal homepage: www.elsevier.com/locate/scihorti

Comparative transcriptomic analysis of white and red Chinese bayberry (Myrica rubra) fruits reveals flavonoid biosynthesis regulation ⁎

Liyu Shia, Xin Chenb, Wei Chenb, Yonghua Zhenga, , Zhenfeng Yangb, a b

T

⁎

College of Food Science and Technology, Nanjing Agricultural University, Nanjing, 210095, People’s Republic of China College of Biological and Environmental Sciences, Zhejiang Wanli University, Ningbo, 315100, People’s Republic of China

A R T I C L E I N F O

A B S T R A C T

Keywords: Chinese bayberry Transcriptome RNA-seq Flavonoid Proanthocyanidin Transcription factors

To fully elucidate molecular mechanisms of flavonoid biosynthesis in Chinese bayberry, transcriptomes of two genotypes with different colours, white cultivar (Shuijing, SJ) and red cultivar (Biqi, BQ), were compared. During fruit development, ‘BQ’ exhibited substantial increase of total anthocyanins content together with dramatic decrease of total soluble proanthocyanidins (PAs) content, while ‘SJ’ displayed decrease of total soluble PAs levels but failed to accumulate anthocyanin. Based on the sequencing results, 124,265 unigenes were generated with an average length of 708 bp. All genes involved in anthocyanin biosynthesis and glycosylation were identified and their expression patterns were in accordance with total anthocyanin accumulation in developing fruits of ‘BQ’ and ‘SJ’. Expression levels of leucoanthocyanidin reductase (LAR) and anthocyanidin reductase (ANR) were down-regulated in ‘SJ’ in agreement with the decrease in PAs, while an opposite trend was observed in ‘BQ’. The higher expression levels of LAR and ANR in ‘BQ’ may be an important reason for the higher levels of total soluble PAs as compared to ‘SJ’. Phylogenetic analysis showed that seven MYBs could be identified as putative homologues of PA-specific regulator and exhibited considerable genotypic and temporal specificity of expression. In addition, three WD40 genes in Chinese bayberry clustered with the WD40s related to flavonoid biosynthesis in other species. Finally, RNA-Seq data was validated by using qRT-PCR analysis with a high correlation, suggesting that the RNA-Seq data here are credible. These results provide new insights into the regulation of the complex branching pathway leading to various flavonoid compounds biosynthesis in bayberries.

1. Introduction

studied. Recently, increasing interest has been attracted to PAs because of their bioactive functions. Indeed, PAs contribute to the quality of many important plant products, such as some berries, tea, and wine and also benefit human health by protecting against free radical injury and cardiovascular disease (Bagchi et al., 2000; Cos et al., 2004; Aron and Kennedy, 2008). However, to our knowledge, there are no published studies on PA biosynthesis in Chinese bayberry fruit. The flavonoids metabolism pathway has been well characterized over the past several decades (Supplementary Fig. S1). The biosynthesis of anthocyanins, PAs and flavonols shares a series of common enzymatic steps, which are catalyzed by phenylalanine ammonia lyase (PAL), cinnamate 4-hydroxylase (C4H), 4-coumarate CoA ligase (4CL), chalcone synthase (CHS), chalcone isomerase (CHI), flavanone 3-hydroxylase (F3H) and flavanone 3′-hydroxylase (F3′H). Activities of

Chinese bayberry (Myrica rubra Sieb. & Zucc.) is a subtropical fruit native to southern China and other Asian countries with high commercial value for its flavor and appealing color. The fruit is rich in anthocyanins as well as other flavonoids compounds and recognized as a good source of natural antioxidants (Bao et al., 2005). Flavonoids are major plant secondary metabolites that contribute in many ways to the growth and survival of plants and also serve as antioxidant and anticancer agents in the human diet (Winkel-Shirley, 2001; Tohge et al., 2005). Anthocyanins, proanthocyanidins (PAs) and flavonols are the predominant flavonoids in fruit (Saito et al., 2013). Since anthocyanin accumulation is responsible for fruit coloration, the genetics and biochemistry of the anthocyanin biosynthetic pathway has been well

Abbreviations: PAs, proanthocyanidins; ANR, anthocyanidin reductase; MYB, Myeloblastosis proto-oncogene; bHLH, basic helix-loop-helix; PAL, phenylalanine ammonia lyase; C4H, cinnamate 4-hydroxylase; 4CL, 4-coumarate CoA ligase; CHS, chalcone synthase; CHI, chalcone isomerase; F3H, flavonoid 3-hydroxylase; F3′H, flavonoid 3′-hydroxylase; DFR, dihydroflavonol 4-reductase; ANS, anthocyanidin synthase; UFGT, UDP-glucose flavonoid 3-O-glucosyl transferase; OMT, O-methyltransferase; ACT, anthocyanin acyltransferase; FLS, flavonol synthase; LAR, leucoanthocyanidin reductase; UGT72L1, uridine diphosphate glycosyl transferase; FPKM, fragments per kilobase of exon per million fragments; RNA-seq, RNA sequencing; qRT-PCR, quantitative real-time PCR ⁎ Corresponding authors. E-mail addresses: [email protected] (Y. Zheng), [email protected] (Z. Yang). https://doi.org/10.1016/j.scienta.2018.02.076 Received 29 October 2017; Received in revised form 22 February 2018; Accepted 26 February 2018 0304-4238/ © 2018 Elsevier B.V. All rights reserved.

Scientia Horticulturae 235 (2018) 9–20

L. Shi et al.

Fig. 1. Images and flavonoids contents of Chinese bayberry fruits. (a) Images of 57 DAFB ‘Biqi’ (BQ1), 85 DAFB ‘Biqi’ (BQ2), 113 DAFB ‘Biqi’ (BQ3), 57 DAFB ‘Shuijin’ (SJ1), 85 DAFB ‘Shuijin’ (SJ2) and 113 DAFB ‘Shuijin’ (SJ3). The white bar represents 1 cm. (b) Contents of total flavonoids, soluble PAs and anthocyanins in BQ and SJ cultivars of three stages. Bars indicate standard errors; different letters indicate statistically significant differences at P < 0.05 according to Duncan’s multiple range test.

2007; Espley et al., 2007). Other MYB transcription factors, PtMYB134, AtTT2 and DkMYB2 specifically regulate PA biosynthesis in persimmon fruit, Arabidopsis seed and poplar leaves, respectively (Nesi et al., 2001; Akagi et al., 2010; Gesell et al., 2014). However, for Chinese bayberry, only two genes, MrMYB1 and MrMYB2, have been isolated and reported to control fruit anthocyanin formation (Niu et al., 2010), while the specific MYB transcription factors involved in PA or flavonol regulation in Chinese bayberry remain unclear. The interaction between bHLH and MYB proteins arises early during land plant evolution and the bHLH member is essential for the activity of the MYB partner (Feller et al., 2011). It has been showed that a given bHLH can interact with several different MYBs, providing multiple functional combinations. For example, the bHLH factor VvMYC1 from grapevine is shown to interact with different MYB proteins (VvMYB5a/ 5b, VvMYBA1/A2 and VvMYBPA1) for inducing promoters of flavonoid pathway genes involved in biosynthesis of anthocyanin and/or PA (Hichri et al., 2010). Likewise, the WD40 factor TTG1 can interact with TT8-TT2 and TT8-AtPAP1 complex to form different ternary transcription complexes that regulate PA and anthocyanin biosynthesis, respectively (Baudry et al., 2004; Shin et al., 2015). Thus, the WD40 cofactors are also adaptable and can control different branches of the flavonoid pathway. The high-throughput RNA sequencing (RNA-Seq) based on nextgeneration sequencing (NGS) and de novo assembly have recently emerged as a powerful approach for gene discovery and expression profiling (Kalra et al., 2013; Chen et al., 2015). In this technique, there

these enzymes lead to the synthesis of dihydroflavonols, which are subsequently either converted by flavonol synthase (FLS) to the flavonols or reduced by dihydroflavonol 4-reductase (DFR) to the leucocyanidins. Anthocyanin synthase (ANS) then catalyses the last common step shared by the anthocyanin and PA biosynthesis pathways. The final modification steps in the anthocyanin pathway include glycosylation by UDP-flavonoid glucosyltransferase (UFGT), methylation by O-methyltransferase (OMT) and acylation by anthocyanin acyltransferase (ACT). While, leucoanthocyanidin reductase (LAR) and anthocyanidin reductase (ANR) are specific to the PA branch pathway and convert leucocyanidin and cyanidin to catechin and epicatechin, respectively (Abrahams et al., 2003; Xie et al., 2003). Subsequently epicatechin can be converted to epicatechin 3′-O-glucoside by the action of UGT72L1, which is an epicatechin-specific glycosyltransferase identified in Medicago (Pang et al., 2013). It has been reported that flavonoid biosynthesis can also be transcriptionally regulated by three types of transcription factors, including R2R3-MYB factors, basic helix–loop–helix (bHLH) proteins, and WD40 proteins (Ramsay and Glover, 2005). These regulators physically interact to form the MBW complex that binds to promoters and activates transcriptions of their target genes (Baudry et al., 2004). It is the MYB transcription factors that seem to impart specificity, and numerous MYB factors controlling different branches of the flavonoid pathway have been identified in many plant species. For instance, the biosynthesis of anthocyanins in grapevine is regulated by VvMYBA1 and VvMYBA2, while VvMYBPA1 is shown to control PAs accumulation (Bogs et al.,

10

Scientia Horticulturae 235 (2018) 9–20

L. Shi et al.

low-quality reads with more than 50% bases with quality value ≤20. The retained high quality reads were further assembled using a de novo assembly program Trinity (Grabherr et al., 2011). For each library, Trinity combined reads with certain length of overlap to form longer sequences called contigs and the reads were mapped back to the contigs. With paired-end reads, it was possible to detect contigs from the same transcript and calculate the distances among these contigs. Finally, the contigs were further assembled to obtain sequences that cannot be extended on either end, and defined as unigenes. Functional annotation of the assembled unigenes was performed by searching against various public databases, including the NCBI nonredundant protein (Nr), the NCBI non-redundant nucleotide (Nt), the Swiss-Prot protein, the Protein family (Pfam), the euKaryotic Ortholog Groups (KOG), the Gene Ontology (GO), and the Kyoto Encyclopedia of Genes and Genomes (KEGG).

is no strict requirement for a reference genome (Lister et al., 2009); therefore, it is particularly attractive for non-model species without available genomic sequences (Wang et al., 2010; Hao et al., 2011; Tang et al., 2011). To date, genes involved in anthocyanin biosynthesis have been identified from the transcriptome of a single genotype, ‘Biqi’ (BQ) with RNA-Seq and the results showed that the expression of these genes were up-regulated during fruit development, concurrent with color change (Feng et al., 2012). Nevertheless, there are extensive gaps in our understanding of PA and flavonol biosynthesis in Chinese bayberry and the differences in flavonoids accumulation amongst different bayberry cultivars. It has been suggested that the molecular mechanisms of flavonoid biosynthesis and regulation might be diverse due to the different distribution of flavonoids among different species (Sun et al., 2015). Therefore, in current study, we compared the transcriptome profiling during the development and ripening of fruit from two Chinese bayberry varieties, white cultivar (Shuijing, SJ) and red cultivar (Biqi, BQ). The structural and regulatory genes involved in flavonoid biosynthesis and accumulation were identified successfully, and their expression patterns were analyzed as well. This work provides an important genetic resource for understanding the molecular mechanism of flavonoid accumulation in Chinese bayberry and the results may be helpful for further molecular breeding of the bayberry.

2.5. Differential gene expression analysis FPKM (fragments per kilobase per million fragments) was used to quantify gene expression, which was able to eliminate the effect of sequence length and sequencing levels on the calculation of gene expression (Mortazavi et al., 2008). Differential expression analysis between two samples (each three biological replicates) was performed by the DESeq R package (Anders and Huber, 2010). P-value corresponds to differential gene expression test. In this experiment, the value of adjusted P-value < 0.05 was set as the threshold to judge the significant DEGs. GO enrichment analysis of DEGs was performed using the GOseq R package, with adjusted P values being less than 0.05 (Young et al., 2010). The statistical enrichment of the DEGs in KEGG pathways was tested using the KOBAS software (Xie et al., 2011). Pathways with corrected P values < 0.05 were considered to be significantly enriched.

2. Materials and methods 2.1. Plant materials and RNA extraction Two types of Chinese bayberry (Myrica rubra Sieb. and Zucc.) varieties, ‘Shuijing’ (SJ) and ‘Biqi’ (BQ), were cultivated in the city of Shaoxing (30°01′N, 120°52′E) and Cixi (30°12′N, 121°24′E), respectively. Fruit were collected at 57, 85, and 113 days after flower bloom (DAFB) from three different trees, as shown in Fig. 1a. At each stage, three biological replicates were designed and the samples of each biological replicate were mixed samples of 45 fruit picked from the three trees (15 fruits per tree). All samples were transported to the laboratory within one hour and selected for uniformity without any damage, then immediately frozen in liquid nitrogen and kept at −80 °C until use. Total RNA was extracted using the Plant RNA kit (OMEGA, Norcross, GA) and purified with DNase I (Omega, Norcross, GA) to eliminate any DNA contamination. The RNA quality was assessed by denaturing agarose gel electrophoresis coupled with NanoDrop ND-1000 (Wilmington, DE, USA). RNA concentration was measured using Qubit RNA Assay Kit in Qubit 2.0 Flurometer (Life Technologies, CA, USA). Further, RNA integrity was determined using a Bioanalyzer 2100 (Agilent Technologies, CA, USA).

2.6. Validation of gene expression profile by qRT-PCR Total RNA was extracted using the same method mentioned above. Approximately 2 μg of total RNA was reverse transcribed into cDNA using the SuperRT First Strand cDNA Synthesis Kit (CWBIO, Beijing, China) following the manufacturer’s instructions. Sixteen unigenes were selected for qRT-PCR analysis to verify the expression patterns determined by RNA-seq study. All analysis was repeated three times using three biological replicates, and the specific primers for the selected genes are listed in Supplementary Table S1. qRT-PCR was performed using SYBR Green I Master Mix (Roche, Penzberg, GER) on a CFX96 Touch Real-Time PCR System (Bio-Rad, California, USA), according to the manufacturer’s protocol. The amplification program consisted of an initial denaturation step at 95 °C for 7 min, followed by 40 cycles of 15 s at 95 °C, and combined with each primer specific annealing temperature ranged from 50 to 60 °C for 30 s. Melting curve analysis was performed at the end of 40 cycles to confirm the specificity of amplification. PCR efficiency was evaluated in preliminary experiments to ensure % efficiency between 90 and 110. The relative expression level of each gene was normalized to MrACT (GenBank accession no. GQ340770) and calculated using the 2−ΔCt method (Livak and Schmittgen, 2001).

2.2. Flavonoid measurements Total anthocyanin contents were determined as previously described (Shi et al., 2014). Total flavonoid levels were measured using the aluminum ion colorimetric (Chang et al., 2002). Total soluble PAs were extracted and analysed as previous described (Porter et al., 1985). 2.3. RNA-seq library construction and sequencing RNA-seq libraries were constructed by using three biological replicates for ‘SJ’ and ‘BQ’ fruit of three different developmental stages, which were named SJ1, SJ2, SJ3 and BQ1, BQ2, BQ3 respectively, according to the manufacturer’s instructions (NEB, Ispawich, USA) (Wang et al., 2009). The libraries were sequenced using Illumina HiSeq™ 2000 and 100 bp single-end reads were generated.

2.7. Phylogenetic analysis and statistical analysis Phylogenetic trees were constructed using MEGA5 with neighbourjoining method and 1000 bootstrap replicates after alignment by ClustalW. The evolutionary distances were computed using the p-distance method. All values are shown as the mean ± standard errors. Statistical analysis was performed using SAS 9.2 software (SAS Institute, North Carolina, USA). Duncan’s multiple range test was used to compare the means at P < 0.05.

2.4. Reads assembly and functional annotation Before assembly, the raw sequences were filtered to remove adaptor sequences, reads in which unknown bases were more than 10%, and 11

Scientia Horticulturae 235 (2018) 9–20

L. Shi et al.

3. Results 3.1. Flavonoid levels of two cultivars among three developmental stages As the Chinese bayberry fruit developed, the levels of the total flavonoids decreased slight firstly and then increased dramatically in both two cultivars with significantly higher total flavonoids contents detected in ‘BQ’ when compared to ‘SJ’ (Fig. 1b). The contents of total soluble PAs in two cultivars decreased during fruit development. Higher levels of total soluble PAs were observed in ‘BQ’ fruit at 85 and 113 DAFB as compared to ‘SJ’ fruit. Total anthocyanin content in ‘BQ’ increased continuously during fruit development giving rise to a characteristic red pigmentation in this cultivar; however, no anthocyanins were detected in ‘SJ’ throughout fruit development (Fig. 1b). 3.2. De novo assembly and annotation Fig. 2. Histogram of GO classification of assembled unigenes.

Eighteen cDNA libraries were built for three developing fruits (57, 85 and 113 DAFB) from the two genotypes of ‘SJ’ and ‘BQ’, with three biological replicates. The major characteristics of these 18 libraries are outlined in Supplementary Table S2. The total number of raw reads in each library was approximately 41.2–60.7 million. All of the raw reads are available in the NCBI SRA database under SRA ID SRP107611 (Bioproject PRJNA386720). After filtering the raw data, the total number of clean tags per library ranged from 39.2 to 59.2 million with a Q20 percentage (percentage of bases with a Phred value > 20) over 93% and a GC percentage about 47% (Supplementary Table S2). These results suggest that the sequencing data were sufficient to ensure accurate de novo assembly. Then the unigene assembly was performed using the Trinity tool. As a result, 185,769 transcripts were obtained, with a mean length of 1170 bp and a N50 of 2475 bp. All the transcripts were then assembled to give 124,265 unigenes with an average size of 708 bp and a N50 of 1460 bp (Table 1). Approximately 124,265 assembled unigenes were functionally annotated using BLASTx against the NCBI non-redundant protein (Nr), Swiss-Prot protein, Protein family (Pfam), euKaryotic Ortholog Groups (KOG), Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) databases, and using BLASTn against the nucleotide (Nt) database, with an E-value cut-off of 1e−5. A total of 71,773 unigenes were matched to the known genes in the public databases, accounting for 57.75% of the total unigenes (Supplementary Table S3). Based on the NCBI Nr data base, 51,948 unigenes (41.8%) showed significant hits. Among them, 59.9% of the matched sequences showed homology (1e−45 < E-value < 1e−5), 19.2% were strong homology (1e−100 < E-value < 1e−45) and the remaining 20.9% (E-value < 1e−100) had very strong homology to the top matched sequences (Supplementary Fig. S2a). With respect to species, 7.5% of these unigenes had top matches with genes from Vitis vinifera, followed by Theobroma cacao (5.2%), Hordeum vulgare (4.7%), Oryza sativa (4.4%) and Prunus mume (4.2%) (Supplementary Fig. S2b). As shown in Fig. 2, 49,321 unigenes (39.69% of all the unigenes) were subjected to GO assignments for three main categories: cellular

component, molecular function, and biological process. In cellular component terms, genes related to ‘cell part’ (22,267, 45.1%) were abundantly expressed. Within the molecular function category, ‘binding’ (25,619, 51.9%) and ‘catalytic activity’ (24,420, 49.5%) were observed to be the two most abundant classes. The dominant subcategories within the biological process category were ‘metabolic process’ (28,866, 58.5%) and ‘cellular process’ (27,521, 55.8%). Based on sequence homology, 31,142 unigenes (25.06% of all the unigenes) were grouped into 30 KOG categories (Supplementary Table S3 and Fig. 3). Out of the 30 KOG categories, group O (Post-translational modification, protein turnover, and chaperones) was the largest category (4215, 13.53%). Totally, 28,390 unigenes (22.84% of all the unigenes) were mapped into 282 KEGG pathways (Supplementary Table S3 and Supplementary File 1). 3.3. Analysis of DEGs across developmental stages and genotypes To identify significant DEGs during fruit development, the transcript abundance of each gene in three developmental stages was compared pairwise and filtered with padj < 0.05 for both BQ and SJ. In total, 3259 (BQ2 vs BQ1), 6172 (BQ3 vs BQ1), and 3897 (BQ3 vs BQ2) genes were either up- or down-regulated (Supplementary File 2), reaching a total of 7617 DEGs in the BQ samples. A total of 5548 DEGs were detected among three SJ libraries, with 1353 DEGs between the SJ2 and SJ1 libraries, 1618 DEGs between the SJ3 and SJ2 libraries, and 4965 DEGs between the SJ3 and SJ1 libraries (Supplementary File 3). Also, a total of 10,437 DEGs were identified between the two genotypes in at least one developmental stage with padj < 0.05 as the threshold (Supplementary File 4). In order to further explore the functions of DEGs, the DEGs identified in the comparisons between different cultivars and different stages of ripening were subjected to KEGG pathway enrichment analysis, respectively (Supplementary File 5). The top eight KEGG pathways enriched for the DEGs identified in most comparison are shown in Supplementary Table S4. Biosynthesis of secondary metabolites such as ‘flavonoid biosynthesis’ and ‘phenylpropanoid biosynthesis’ represented the dominant pathways associated with most of gene upregulated during development in both cultivars. Notably, ‘flavonoid biosynthesis’ was significantly enriched in the up-regulated genes between BQ2 vs BQ1, BQ3 vs BQ1, BQ3 vs BQ2, SJ3 vs SJ1, and SJ3 vs SJ2 groups but with more gene members between stages in BQ, reflecting a distinctly more active flavonoid synthesis metabolism during fruit development of BQ as compared to SJ. Consistent with this, ‘flavonoid biosynthesis’ was significantly down-regulated in SJ2 vs BQ2 and SJ3 vs BQ3 comparisons. During fruit development, ‘photosynthesis-antenna proteins’ and ‘photosynthesis’ was included in the down-

Table 1 Summary of the assembly statistics. Category

Transcript

Unigene

Total number (≥200 bp) ≥500 bp ≥1000 bp Total length Min length Mean length Max length N50 N90

185,769 92,553 65,041 217,310,781 201 1170 17,341 2475 398

124,265 38,544 20,703 87,935,581 201 708 17,341 1460 254

12

Scientia Horticulturae 235 (2018) 9–20

L. Shi et al.

Fig. 3. Histogram of KOG classifications of assembled unigenes.

regulated pathways. Also, some metabolism such as ‘cyanoamino acid metabolism’ and ‘starch and sucrose metabolism’ were significantly enriched in the down-regulated genes between BQ2 vs BQ1, BQ3 vs BQ1 and BQ3 vs BQ2, but showed a different picture in SJ fruit. Moreover, BQ exhibited higher enrichment in ‘plant hormone signal transduction’ during development as compared to SJ. 3.4. Candidate genes involved in flavonoid biosynthesis In plants, flavonoids are a class of important secondary metabolites, which play roles in diverse physiological processes (Buer et al., 2010). The flavonoid biosynthetic pathway has been well characterized and a brief schematic is shown in Supplementary Fig. S1. We identified 38 candidate transcripts encoding almost all known enzymes involved in flavonoid biosynthesis (Fig. 4). The biosynthesis of anthocyanins, PAs and flavonols share the upstream phenylpropanoid pathway, comprising PAL, C4H, 4CL, CHS, CHI, F3H and F3’H. Three PAL genes were identified, of which one unigenes (c37792_g1) exhibited the highest expression level of the three PAL genes and showed increased expression in both ‘BQ’ and ‘SJ’ during fruit development. Up-regulation of unigene (c50014_g1) expression was evident during fruit development, while the other unigene (c33196_g1) showed a pattern of deceasing expression. Only one unigene (c47448_g2) encoding C4H was found, which was up-regulated throughout fruit development in both bayberry cultivars and showed highest the expression at 113 DAFB in ‘BQ’. Three 4CL genes (c42854_g1, c46876_g1 and c47756_g1) showed two distinct expression patterns between ‘BQ’ and ‘SJ’ cultivars; the expression levels were high and increased throughout fruit development in ‘BQ’, while they were low and no significant changes were found in ‘SJ’ during development. Two CHS genes (c38429_g1 and c39831_g1) displayed different expression patterns, with the expression level of the former being approximately 100 times higher than that of the latter during fruit development in both cultivars. CHI (c43930_g1), F3H (c47207_g1) and F3’H (c49836_g1) had similar expression patterns in

Fig. 4. Heat map diagram of the expression of genes related to flavonoid biosynthesis in two Chinese bayberry cultivars during development. Each stage of fruit development in ‘Biqi’ (BQ) and ‘Shuijin’ (SJ) are listed horizontally and the numbers 1–3 indicate 57, 85 and 113 DAFB, respectively. The unigenes are listed on the vertical line. The annotations are displayed on the right side. All of the FPKM (fragments per kilobase per million fragments) values of the unigenes are shown as logarithms.

both cultivars during ripening and showed highest the expression at 85 and 113 DAFB in ‘BQ’. In the anthocyanin branch, two unigenes (c48342_g1 and c52293_g3) encoding DFR were identified, both of which showed high expression during the later stages of fruit development (85 and 113 13

Scientia Horticulturae 235 (2018) 9–20

L. Shi et al.

DAFB) in both cultivars. A higher expression level of these two DFR genes was observed in ‘BQ’ fruit than ‘SJ’ during fruit development. One unigene (c51485_g1) encoding ANS was significantly up-regulated throughout fruit development in ‘BQ’. In contrast, its expression remained relatively constant and was much lower in ‘SJ’ during fruit development. The final modification steps in the anthocyanin pathway include glycosylation by UDP-flavonoid glucosyltransferase (UFGT), methylation by O-methyltransferase (OMT) and acylation by anthocyanin acyltransferase (ACT). Two genes (c43204_g1 and c50569_g3) encoding UFGT were identified and they showed low expression in ‘SJ’ at all stages. In ‘BQ’ fruit, the highest expression of c43204_g1 gene was detected at 85 and 113 DAFB, while c50569_g3 gene at 57 and 85 DAFB. The transcript abundance of OMT (c38496_g1) increased slightly during fruit development in two cultivars. Meanwhile, the highest expression of ACT (c12459_g1) was detected at 113 DAFB in ‘BQ’. For the PA pathway, two genes encoding leucoanthocyanidin reductase (LAR, c53873_g1) and anthocyanidin reductase (ANR, c35023_g1) were identified in the Chinese bayberry database. LAR revealed the highest expression level at 57 DAFB in ‘SJ’ and showed a gradual increase with ripening in ‘BQ’ fruit. One ANR gene was expressed at higher levels in ‘BQ’ than in ‘SJ’. UGT72L1, a glycosyltransferase involved in the production of epicatechin 3′-O-glucoside as a key step in PA biosynthesis or its regulation, was also investigated in present study. One unigene (c20738_g1), with 69% identity to the previously reported Medicago UGT72L1 gene (Pang et al., 2008), displayed low expression at 57 and 85 DAFB and highest expression at 113 DAFB in ‘BQ’, whereas it remained at high level and did not experienced significant changes in three stages in ‘SJ’. Two flavonol synthase genes (FLS, c33081_g1 and c47053_g2) were also identified: the expression of the former exhibited a greater increase in ‘BQ’ than in ‘SJ’ during fruit development while the latter lower in ‘BQ’ than in ‘SJ’ at all stages.

Table 2 Summary of transcription factor related unigenes identified in the Chinese bayberry fruit transcriptome. Family

Numbers of unigenes

Family

Numbers of unigenes

C2H2 WD40 ERF MYB-related bZIP C3H bHLH Orphans HB WRKY NAC MYB GRAS zn-clus C2C2-GATA FAR1 C2C2-Dof B3 LOB Trihelix G2-like SBP HSF TCP M-type Tify TUB OFP ARF zf-HD EIL

189 135 93 86 79 79 74 70 65 62 56 47 40 40 30 29 27 26 26 26 25 21 19 18 17 14 14 11 10 10 9

NF-YA NF-YC BES1 C2C2-YABBY DBB PLATZ CSD MIKC NF-YB AP2 BBR-BPC CPP E2F-DP GeBP SRS Alfin-like RAV RWP-RK S1Fa-like Whirly ARR-B C2C2-LSD DBP GRF HRT NF-X1 VOZ BSD C2C2-CO-like STAT ULT

8 8 7 7 7 7 6 6 6 5 5 5 5 5 5 4 3 3 3 3 2 2 2 2 2 2 2 1 1 1 1

values of these genes were further analysed. As shown in Fig. 5b, c28754_g2 exhibited the highest expression level of among the three Sg4 members throughout fruit development and showed highest the expression at 57 DAFB in ‘SJ’, while the rest two genes (c24596_g1 and c48297_g1) peaked at 57 DAFB and at 113 DAFB in ‘BQ’, respectively. Among the seven Sg5 members, c14334_g1, c90137_g1 and c40743_g1 exhibited a dramatic decrease in expression throughout fruit development and higher expression levels were detected in ‘SJ’ than ‘BQ’. In contrast, c26286_g1 and c16738_g1 were up-regulated at all stages and the highest expression levels were found at 113 and 85 DAFB in ‘BQ’, respectively. Additionally, c101161_g1 and c93152_g1 displayed highest expression at 57 DAFB in ‘BQ’, but with lower levels observed in ‘SJ’. Two Sg6 genes (c46923_g1 and c42295_g1) showed up-regulated expression pattern during fruit development in two cultivars with higher expression levels in ‘BQ’ than ‘SJ’. The expression of c40428_g1 decreased during fruit development in ‘BQ’, while an opposite trend in ‘SJ’. Based on the previous classification system of Arabidopsis thaliana (Heim et al., 2003; Carretero-Paulet et al., 2010), 129 AtbHLHs, 47 Chinese bayberry bHLHs and 8 flavonoids-related bHLHs of other plants were constructed in to a phylogenetic tree and divided into 10 groups (inside number) and 20 subfamilies (outside number; Fig. 6a). Notably, members of subfamily 2, 5, and 24 were identified as regulators in flavonoid metabolism (Carretero-Paulet et al., 2010). In the present study, there were four Chinese bayberry bHLHs in subfamily 2 (c50876_g4, c34055_g1, c44081_g1 and c45031_g1), three in subfamily 5 (c37607_g1, c43002_g1 and c49156_g1) and four in subfamily 24 (c52223_g3, c66477_g1, c5659_g1 and c42632_g1). Moreover, c37607_g1, c43002_g1 and c49156_g1 belonging to subfamily 5 were closely related to bHLHs for anthocyanin and (or) proanthocyanidin biosynthesis regulation in other plant species. The expression profiles of the 11 Chinese bayberry bHLHs from subfamily 2, 5, and 24 are presented in Fig. 6b. Remarkably, two subfamily 5 members, c37607_g1 and

3.5. Candidate transcription factors for flavonoid biosynthesis In order to identify transcription factor encoding genes in Chinese bayberry transcriptome, BLASTx analysis was performed against the PlnTFDB database. A total of 1466 unigenes showed similarity to transcription factors and could be assigned to 62 families (Table 2). Out of these, some transcription factors are known to be involved in the regulation of flavonoid biosynthesis, such as R2R3-MYB, basic helixloop-helix (bHLH) proteins and WD-repeat-containing proteins (Hichri et al., 2011). In this study, a total of 126 Arabidopsis MYBs, 32 Chinese bayberry MYBs and 13 flavonoids-related MYBs in other plant species were used to construct a phylogenetic tree divided into 25 MYB subgroups overall (Fig. 5a). MYBs that control anthocyanin, proanthocyanidin and flavonol accumulation were mainly identified in Sg4, Sg5, Sg6 and Sg7 reported by Hichri et al. (2011) previously. Three unigenes (c24596_g1, c28754_g2 and c48297_g1) were clustered with transcriptional repressors (FaMYB1 and VvMYBC2-L1) of Sg4 (Aharoni et al., 2001; Huang et al., 2014), while seven unigenes (c26286_g1, c14334_g1, c101161_g1, c90137_g1, c93152_g1, c40743_g1 and c16738_g1) were closely related to PtMYB134, AtTT2, DkMYB2 and VvMYBPA1, candidate genes for proanthocyanidin biosynthesis regulation (Nesi et al., 2001; Bogs et al., 2007; Akagi et al., 2010; Gesell et al., 2014) within Sg5. Similarly, c46923_g1 and c42295_g1 were in the same subgroup (Sg6) with MdMYB1, MdMYB10, VvMYBA1 and VvMYBA2, known to regulate anthocyanins (Kobayashi et al., 2002; Espley et al., 2007; Walker et al., 2007; Hu et al., 2016). And the two unigenes encode proteins with 99% amino acid sequence similarity to MrMYB1 and MrMYB2, respectively, both of which have been reported to control Chinese bayberry fruit anthocyanin formation (Niu et al., 2010). Moreover, members of Sg7 included c40428_g1 and three flavonolregulating MYBs AtPFG1, AtPFG2 and AtPFG3 (Stracke et al., 2010). To investigate the expression patterns of these thirteen Chinese bayberry MYBs during fruit development in ‘BQ’ and ‘SJ’, the FPKM 14

Scientia Horticulturae 235 (2018) 9–20

L. Shi et al.

Fig. 5. Evolutionary relationships of R2R3-MYBs and validation of Chinese bayberry MYBs. (a) Evolutionary relationships of 126 AtMYBs, 32 Chinese bayberry MYBs marked with the black dot, and 13 flavonoids-related MYBs of other plants highlighted with the blue dot. All ambiguous positions were removed for each sequence pair. (b) Heat map for the Chinese bayberry MYBs genes in subgroups 4, 5, 6 and 7. Each stage of fruit development in ‘Biqi’ (BQ) and ‘Shuijin’ (SJ) are listed horizontally and the numbers 1–3 indicate 57, 85 and 113 DAFB, respectively. The unigenes are listed on the vertical line. The scale represents the logarithms of the FPKM values of the unigenes.

15

Scientia Horticulturae 235 (2018) 9–20

L. Shi et al.

Fig. 6. Evolutionary relationships of bHLHs and validation of Chinese bayberry bHLHs. (a) Evolutionary relationships of 129 AtbHLHs, 46 Chinese bayberry bHLHs marked with the black dot, and 8 flavonoids-related bHLHs of other plants highlighted with the red dot. All ambiguous positions were removed for each sequence pair. (b) Heat map for the Chinese bayberry bHLHs genes in subfamilies 2, 5 and 24. Each stage of fruit development in ‘Biqi’ (BQ) and ‘Shuijin’ (SJ) are listed horizontally and the numbers 1–3 indicate 57, 85 and 113 DAFB, respectively. The unigenes are listed on the vertical line. The scale represents the logarithms of the FPKM values of the unigenes.

16

Scientia Horticulturae 235 (2018) 9–20

L. Shi et al.

compounds (Zhang et al., 2008). However, the genomic data available for Chinese bayberry was still limited, which hinders the study of the molecular mechanisms of flavonoids accumulation in Chinese bayberry. Here, two types of Chinese bayberry varieties with distinctively different colours, ‘BQ’ and ‘SJ’ were analyzed by RNA-seq to explore the underlying differences responsible for flavonoid content during fruit development and ripening. Based on the sequencing results, 124,265 unigenes were generated with a mean length of 708 bp, which is considerably longer than the previously reported 531 bp for Chinese bayberry (Feng et al., 2012). Because the expression patterns of genes determined by RNA-Seq and qRT-PCR analyses were largely consistent (Fig. 8b), we are confident that our RNA-seq data set will be broadly useful for studies of fruit transcriptome dynamics. Moreover, all transcripts involved in flavonoid biosynthesis, modification and regulation were identified in the transcriptome. Obvious cultivar specificity for ‘BQ’ and ‘SJ’ was observed in the anthocyanin biosynthesis and regulation. Specifically, ‘BQ’ displayed a significant up-regulated expression of many anthocyanin biosynthetic genes during fruit development, including PAL (c37792_g1), C4H (c47448_g2), 4CL (c42854_g1 and c46876_g1), CHS (c38429_1), CHI (c43930_g1), F3H (c47207_g1), F3’H (c49836_g1), DFR (c48342_g1 and c53393_g3), ANS (c51485_g1) and UFGT (c43204_g1), which was not observed in ‘SJ’ (Fig. 4). These transcript profiles coincided with the total anthocyanin content in these two cultivars, which accumulated in ‘BQ’ to a significantly higher extent compared to ‘SJ’ in all stages of fruit development (Fig. 1b). These findings confirm the high correlation of fruit anthocyanin content and the transcript levels of anthocyanin biosynthetic genes during fruit development in Chinese bayberry as previously reported by Feng et al. (2012) and Shi et al. (2014). The final modification steps in the anthocyanin pathway include methylation by O-methyltransferase (OMT) and acylation by anthocyanin acyltransferase (ACT) (Supplementary Fig. S1). The expression of OMT has been reported to be in accordance with anthocyanin content in grape (Azuma et al., 2009). Our study here reported that one OMT (c38496_g1) was up-regulated in ‘BJ’ and expressed in higher level in ‘BQ’ than in ‘SJ’ during fruit development, which is also probably consistence with the anthocyanin content (Fig. 4). Additionally, one ACT (c12459_g1) displayed the highest expression at 113 DAFB in ‘BQ’, the stage with the highest anthocyanin concentration which implied that c12459_g1 may be involved in the modification of anthocyanin in Chinese bayberry as well. PAs are synthesized via a branch of anthocyanin biosynthesis pathway under the catalyzation of two enzymes, leucoanthocyanidin reductase (LAR) and anthocyanidin reductase (ANR). In the present study, ‘BQ’ displayed substantial increase of total anthocyanin together with dramatic decrease of total soluble PAs during fruit development (Fig. 1b). Similar to ‘BQ’, the content of total soluble PAs in ‘SJ’ was also reduced throughout fruit development. However, both LAR (c53873_g1) and ANR (c35023_g1) showed two distinct expression patterns between the two genotypes. The expression levels of the two genes were down-regulated in ‘SJ’ totally in agreement with the decrease in PAs, while an opposite trend was observed in ‘BQ’ (Fig. 4). The reasons for the discrepancy between LAR and ANR expression and PAs content in ‘BQ’ are still needed further study. Despite this, the higher expression levels of LAR and ANR at 113 DAFB in ‘BQ’ may be an important reason for the higher levels of total soluble PAs detected at 113 DAFB in ‘BQ’ as compared to ‘SJ’. Moreover, one unigene (c20738_g1), the homolog of UGT72L1, which was active specifically toward the PA precursor (-)-epicatechin in Medicago (Pang et al., 2013), appeared to have higher expression levels in ‘BQ’ at 113 DAFB. The metabolic flux in the flavonoid biosynthetic pathway is controlled by substrate competition between the FLS and the DFR (Gou et al., 2011). In a recent study, the competing expression patterns of FLS and DFR were shown to be important determinants of flower color phenotype, via the regulation of the metabolic fate of the common substrates (Luo et al., 2016). Similar results were also observed in the

Fig. 7. Validation of Chinese bayberry WD40s. (a) Phylogenetic tree of flavonoids-related WD40s of Chinese bayberry and other plants. The blue dot represents WD40s in the Chinese bayberry transcriptome dataset. (b) Heat map for 4 flavonoids-related WD40s genes of Chinese bayberry. Each stage of fruit development in ‘Biqi’ (BQ) and ‘Shuijin’ (SJ) are listed horizontally and the numbers 1–3 indicate 57, 85 and 113 DAFB, respectively. The unigenes are listed on the vertical line. The scale represents the logarithms of the FPKM values of the unigenes.

c49156_g1, exhibited the highest expression level among the 11 Chinese bayberry bHLHs and remained relatively constant at all stages. Downregulation of c43002_g1 expression was evident during fruit development in two cultivars, with little difference between ‘BQ’ and ‘SJ’. The same pattern was observed for all genes in subfamily 24, whose expressions were low and decreased during fruit development with no difference between the two cultivars. All genes of subfamily 2 remained relative stable throughout fruit development, except for c34055_g1, which peaked at 57 DAFB in ‘BQ’ and at 85 DAFB in ‘SJ’. A total of 135 WD genes were isolated from the Chinese bayberry transcriptome dataset. Of these, only three genes (c49129_g1, c55619_g1 and c49942_g1) were found in the same group as WD40s involved in the regulation of the flavonoid pathway, such as PgWD40 (Ben-Simhon et al., 2011), FaTTG1 (Schaart et al., 2013), MdTTG1 (Brueggemann et al., 2010), and InWDR1 (Morita et al., 2006) (Fig. 7a). As shown in Fig. 7b, the expressions of c49129_g1 and c49942_g1 were both slightly up-regulated at all stages and higher expression levels were detected in ‘BQ’ than ‘SJ’. However, the expression level of c55619_g1 was low in two cultivars and only detected at 57 and 85 DAFB in ‘BQ’ and at 113 DAFB in ‘SJ’. 3.6. qRT-PCR confirmation of RNA-seq data To validate the accuracy of the RNA-Seq results, we selected 10 structural genes and 6 transcription factors involved in flavonoid biosynthesis and analyzed their expression levels at three stages in two cultivars by qRT-PCR (Fig. 8a). The expression levels of these 16 transcripts were consistent with those obtained by using RNA-Seq (Figs. 4 and 5b). The Pearson correlation coefficient was used to determine the correlation between different platforms. As a result, high correlation (the correlation coefficient of 0.863) was found between the RNA-seq data and the qRT-PCR results (Fig. 8b). Thus, the RNA-Seq data here are credible and can be used for subsequent experiments. 4. Discussion Chinese bayberry (Myrica rubra Sieb. and Zucc.), an important subtropical fruit crop, is a fruit enriched in flavonoid bioactive 17

Scientia Horticulturae 235 (2018) 9–20

L. Shi et al.

Fig. 8. qRT-PCR validation of 16 selected transcripts. (a) Expression patterns of 10 structural genes and 6 transcription factors involved in flavonoid biosynthesis. (b) Correlation of the expression levels of 16 selected genes measured by qRT-PCR and RNA-seq.

WD40 proteins that form a MYB-bHLH-WD40 (MBW) transcriptional activation complex (Dubos et al., 2010). In order to identify putative MYBs that regulate flavonoid biosynthesis in Chinese bayberry, 126 Arabidopsis MYBs, 32 Chinese bayberry MYBs and 13 identified flavonoids-related MYBs in other plant species were used to construct a phylogenetic tree (Fig. 5a). The Phylogenetic analyses showed that 13 Chinese bayberry MYBs were found to cluster with the identified flavonoids-related MYBs within Sg4, Sg5, Sg6 and Sg7, which have been reported to include the MYBs members controlling different branches of the flavonoid pathway (Hichri et al., 2011). Among the 13 Chinese

red (‘BQ’) and white (‘SJ’) colored Chinese bayberry fruits in our study. In ‘BQ’, expression of the DFR (c48342_g1) was higher than that of FLS (c47053_g2) during fruit development, which could lead to the increase in anthocyanin levels. However, the opposite trend was observed in ‘SJ’. Additionally, as compared to ‘BQ’ fruit, higher expression of FLS (c47053_g2) and lower transcripts of DFR (c48342_g1) were found in ‘SJ’ throughout fruit development. Altogether, these results might account for the no anthocyanin accumulation in this white cultivar. Biosynthesis of the three diverse groups of flavonoids, anthocyanins, PAs, and flavonols, is controlled by MYB combined with bHLH and 18

Scientia Horticulturae 235 (2018) 9–20

L. Shi et al.

patterns of genes involved in anthocyanin biosynthesis, and glycosylation coincided with the total anthocyanin content in bayberries. In addition, the transcript profiles of LAR and ANR were associated with differences in the PA content between the two cultivars. Moreover, 13 MYB genes, 11 bHLH genes and 3 WD40 genes could be identified as putative homologues of transcription factors controlling different branches of the flavonoid pathway in other species. Their expression profiles were analysis and their possible roles in flavonoid biosynthesis regulation were discussed. This study facilitates our understanding of the molecular mechanism of flavonoid accumulation in Chinese bayberry and provides an important datasets for further related studies.

bayberry MYBs, two unigenes, c46923_g1 and c42295_g1, encoding MrMYB1 and MrMYB2, respectively, were clustered with MdMYB1, MdMYB10, VvMYBA1 and VvMYBA2, known to regulate anthocyanin biosynthesis (Kobayashi et al., 2002; Espley et al., 2007; Walker et al., 2007; Hu et al., 2016). The two unigenes were up-regulated in two cultivars and higher expression levels were detected in ‘BQ’ than ‘SJ’ (Fig. 5b), as well as many anthocyanin pathway genes. Seven unigenes were closely related to PtMYB134, AtTT2, DkMYB2 and VvMYBPA1, candidate genes for PAs biosynthesis regulation (Nesi et al., 2001; Bogs et al., 2007; Akagi et al., 2009; Gesell et al., 2014), and collectively showed three different expression patterns (Fig. 5b), which implied exhibit considerable genotypic and temporal specificity of expression and each of them may play major roles in the differential expression of LAR and ANR in Chinese bayberry. Among the three Sg4 members (c24596_g1, c28754_g2 and c48297_g1), c28754_g2 was clustered with transcriptional repressors reported to be involved in flavonoid biosynthesis, FaMYB1 (Aharoni et al., 2001) and VvMYBC2-L1 (Huang et al., 2014), and showed the highest expression level throughout fruit development, indicating that c28754_g2 is most likely a negative MYB regulator of flavonoid accumulation. bHLH proteins form one of the largest transcription factor families in plants, and are structurally organized into basic helix-loop-helix DNA-binding conserved motifs (Jones, 2004). In Chinese bayberry transcriptome dataset, three unigenes (c37607_g1, c43002_g1 and c49156_g1) were closely related to the known plant bHLHs involved in anthocyanin and (or) PA biosynthesis regulation. Unigene c49156_g1 and Unigene c37607_g1 were found to encode proteins with 99% amino acid sequence similarity to MrbHLH1 and MrbHLH2, respectively. Previous study has reported that MrbHLH1, but not MrbHLH2, is the essential partner of MrMYB1 during anthocyanin biosynthesis regulation in bayberry (Liu et al., 2013). However, in our study, the two unigenes were most closely related to the grapevine VvMYC1 and VvMYCA1, both of which have been claimed to control anthocyanin and PA biosynthesis (Hichri et al., 2010; Matus et al., 2010). Besides, the two unigenes were highly expressed during fruit development in two cultivars, with little difference between ‘BQ’ and ‘SJ’, which was not reflected in the respective levels of anthocyanin (Fig. 6b). Thus, future studies will focus on the potential of c49156_g1 and c37607_g1 for engineering the PA pathway in Chinese bayberry. As a member of the MBW transcription complex, WD40 proteins are not thought to participate in the specific recognition of a target gene promoter, but rather seem to enhance gene activation in the MBW complex regulating flavonoid biosynthesis (Hichri et al., 2011). For instance, the WD40 factor TTG1 is required for normal expression of DFR and BAN and, therefore, for anthocyanin and PA synthesis in Arabidopsis (Walker et al., 1999). For Chinese bayberry, MrWD40-1 and MrWD40-2 have been isolated and, of them, only MrWD40-1 has been shown to interact with MrMYB1 and MrbHLH1 to enhance anthocyanin accumulation (Liu et al., 2013), but their function in PA biosynthesis is not clear yet. In this study, c49129_g1 and c49942_g1 encoding MrWD40-1 and MrWD40-2, respectively, as well as c55619_g1 were clustered with the WD40s related to flavonoid biosynthesis in other species (Fig. 7a). Unigene c49129_g1 was most closely related to pomegranate PgWD40, but was distant from c49942_g1 and c55619_g1. The expressions of c49129_g1 and c49942_g1 exhibited slightly upregulated trends during fruit development and higher expression levels were detected in ‘BQ’ than ‘SJ’. Further studies are still needed to determine whether these three genes are related to the PA biosynthesis regulation in Chinese bayberry.

Acknowledgements We thank the faculty and staff of Biomarker Technologies Corporation at Beijing for help in RNA library construction and sequencing. This study was supported by the National Natural Science Foundation of China (31671898) and the Innovation Research Foundation from Zhejiang Top Key Discipline of Bioengineering (CX2016002). Appendix A. Supplementary data Supplementary material related to this article can be found, in the online version, at doi:https://doi.org/10.1016/j.scienta.2018.02.076. References Abrahams, S., Lee, E., Walker, A.R., Tanner, G.J., Larkin, P.J., Ashton, A.R., 2003. The Arabidopsis TDS4 gene encodes leucoanthocyanidin dioxygenase (LDOX) and is essential for proanthocyanidin synthesis and vacuole development. Plant J. 35, 624–636. Aharoni, A., De Vos, C.H.R., Wein, M., Sun, Z.K., Greco, R., Kroon, A., Mol, J.N.M., O’connell, A.P., 2001. The strawberry FaMYB1 transcription factor suppresses anthocyanin and flavonol accumulation in transgenic tobacco. Plant J. 28, 319–332. Aron, P.M., Kennedy, J.A., 2008. Flavan-3-ols: nature, occurrence and biological activity. Mol. Nutr. Food Res. 52, 79–104. Akagi, T., Ikegami, A., Tsujimoto, T., Kobayashi, S., Sato, A., Kono, A., Yonemori, K., 2009. DkMyb4 is a Myb transcription factor involved in proanthocyanidin biosynthesis in persimmon fruit. Plant Physiol. 151, 2028–2045. Akagi, T., Ikegami, A., Yonemori, K., 2010. DkMyb2 wound-induced transcription factor of persimmon (Diospyros kaki Thunb.), contributes to proanthocyanidin regulation. Planta 232, 1045–1059. Anders, S., Huber, W., 2010. Differential expression analysis for sequence count data. Genome Biol. 11, R106. Azuma, A., Kobayashi, S., Goto-Yamamoto, N., Shiraishi, M., Mitani, N., Yakushiji, H., Koshita, Y., 2009. Color recovery in berries of grape (Vitis vinifera L.)’ Benitaka’, a bud sport of ‘Italia’, is caused by a novel allele at the VvmybA1 locus. Plant Sci. 176, 470–478. Bagchi, D., Bagchi, M., Stohs, S.J., Das, D.K., Ray, S.D., Kuszynski, C.A., Joshi, S.S., Pruess, H.G., 2000. Free radicals and grape seed proanthocyanidin extract: importance in human health and disease prevention. Toxicology 148, 187–197. Bao, J., Cai, Y., Sun, M., Wang, G., Corke, H., 2005. Anthocyanins, flavonols, and free radical scavenging activity of Chinese bayberry (Myrica rubra) extracts and their color properties and stability. J. Agric. Food Chem. 53, 2327–2332. Baudry, A., Heim, M.A., Dubreucq, B., Caboche, M., Weisshaar, B., Lepiniec, L., 2004. TT2, TT8, and TTG1 synergistically specify the expression of BANYULS and proanthocyanidin biosynthesis in Arabidopsis thaliana. Plant J. 39, 366–380. Ben-Simhon, Z., Judeinstein, S., Nadler-Hassar, T., Trainin, T., Bar-Ya’akov, I., BorochovNeori, H., Holland, D., 2011. A pomegranate (Punica granatum L.) WD40-repeat gene is a functional homologue of Arabidopsis TTG1 and is involved in the regulation of anthocyanin biosynthesis during pomegranate fruit development. Planta 234, 865–881. Bogs, J., Jaffe, F.W., Takos, A.M., Walker, A.R., Robinson, S.P., 2007. The grapevine transcription factor VvMYBPA1 regulates proanthocyanidin synthesis during fruit development. Plant Physiol. 143, 1347–1361. Brueggemann, J., Weisshaar, B., Sagasser, M., 2010. A WD40-repeat gene from Malus domestica is a functional homologue of Arabidopsis thaliana TRANSPARENT TESTA GLABRA1. Plant Cell Rep. 29, 285–294. Buer, C.S., Imin, N., Djordjevic, M.A., 2010. Flavonoids: new roles for old molecules. J. Integr. Plant Biol. 52, 98–111. Carretero-Paulet, L., Galstyan, A., Roig-Villanova, I., Martinez-Garcia, J.F., Bilbao-Castro, J.R., Robertson, D.L., 2010. Genome-wide classification and evolutionary analysis of the bHLH family of transcription factors in Arabidopsis, poplar, rice, moss, and algae. Plant Physiol. 153, 1398–1412. Chang, C.C., Yang, M.H., Wen, H.M., Chern, J.C., 2002. Estimation of total flavonoid content in propolis by two complementary colorimetric methods. J. Food Drug Anal

5. Conclusions RNA-seq-based comparative transcriptome strategy was used to obtain information regarding the molecular mechanisms underlying flavonoid accumulation in Chinese bayberry. All genes involved in flavonoid biosynthesis and modification were identified. The expression 19

Scientia Horticulturae 235 (2018) 9–20

L. Shi et al.

Nesi, N., Jond, C., Debeaujon, I., Caboche, M., Lepiniec, L., 2001. The Arabidopsis TT2 gene encodes an R2R3 MYB domain protein that acts as a key determinant for proanthocyanidin accumulation in developing seed. Plant Cell 13, 2099–2114. Niu, S.S., Xu, C.J., Zhang, W.S., Zhang, B., Li, X., Lin-Wang, K., Ferguson, I.B., Allan, A.C., Chen, K.S., 2010. Coordinated regulation of anthocyanin biosynthesis in Chinese bayberry (Myrica rubra) fruit by a R2R3 MYB transcription factor. Planta 231, 887–899. Pang, Y., Cheng, X., Huhman, D.V., Ma, J., Peel, G.J., Yonekura-Sakakibara, K., Saito, K., Shen, G., Sumner, L.W., Tang, Y., Wen, J., Yun, J., Dixon, R.A., 2013. Medicago glucosyltransferase UGT72L1: potential roles in proanthocyanidin biosynthesis. Planta 238, 139–154. Pang, Y., Peel, G.J., Sharma, S.B., Tang, Y., Dixon, R.A., 2008. A transcript profiling approach reveals an epicatechin-specific glucosyltransferase expressed in the seed coat of Medicago truncatula. Proc. Natl. Acad. Sci. U. S. A. 105, 14210–14215. Porter, L.J., Hrstich, L.N., Chan, B.G., 1985. The conversiton of procyanidins and prodelphinidins to cyanidin and delphindin. Phytochemistry 25, 223–230. Ramsay, N.A., Glover, B.J., 2005. MYB-bHLH-WD40 protein complex and the evolution of cellular diversity. Trends Plant Sci. 10, 63–70. Saito, K., Yonekura-Sakakibara, K., Nakabayashi, R., Higashi, Y., Yamazaki, M., Tohge, T., Fernie, A.R., 2013. The flavonoid biosynthetic pathway in Arabidopsis: structural and genetic diversity. Plant Physiol. Biochem. 72, 21–34. Schaart, J.G., Dubos, C., Romero De La Fuente, I., Van Houwelingen, A.M., De Vos, R.C., Jonker, H.H., Xu, W., Routaboul, J.M., Lepiniec, L., Bovy, A.G., 2013. Identification and characterization of MYB-bHLH-WD40 regulatory complexes controlling proanthocyanidin biosynthesis in strawberry (Fragaria ananassa) fruits. New Phytol. 197, 454–467. Shi, L., Cao, S., Shao, J., Chen, W., Zheng, Y., Jiang, Y., Yang, Z., 2014. Relationship between sucrose metabolism and anthocyanin biosynthesis during ripening in Chinese bayberry fruit. J. Agric. Food Chem. 62, 10522–10528. Shin, D.H., Cho, M., Choi, M.G., Das, P.K., Lee, S.K., Choi, S.B., Park, Y.I., 2015. Identification of genes that may regulate the expression of the transcription factor production of anthocyanin pigment 1 (PAP1)/MYB75 involved in Arabidopsis anthocyanin biosynthesis. Plant Cell Rep. 34, 805–815. Stracke, R., Jahns, O., Keck, M., Tohge, T., Niehaus, K., Fernie, A.R., Weisshaar, B., 2010. Analysis of PRODUCTION OF FLAVONOL GLYCOSIDES-dependent flavonol glycoside accumulation in Arabidopsis thaliana plants reveals MYB11-, MYB12-and MYB111independent flavonol glycoside accumulation. New Phytol. 188, 985–1000. Sun, H., Liu, Y., Gai, Y., Geng, J., Chen, L., Liu, H., Kang, L., Tian, Y., Li, Y., 2015. De novo sequencing and analysis of the cranberry fruit transcriptome to identify putative genes involved in flavonoid biosynthesis, transport and regulation. BMC Genomics 16, 652. Tang, Q., Ma, X., Mo, C., Wilson, I.W., Song, C., Zhao, H., Yang, Y., Fu, W., Qiu, D., 2011. An efficient approach to finding Siraitia grosvenorii triterpene biosynthetic genes by RNA-seq and digital gene expression analysis. BMC Genomics 12, 343. Tohge, T., Matsui, K., Ohme-Takagi, M., Yamazaki, M., Saito, K., 2005. Enhanced radical scavenging activity of genetically modified Arabidopsis seeds. Biotechnol. Lett. 27, 297–303. Walker, A.R., Davison, P.A., Bolognesi-Winfield, A.C., James, C.M., Srinivasan, N., Blundell, T.L., Esch, J.J., Marks, M.D., Gray, J.C., 1999. The TRANSPARENT TESTA GLABRA1 locus, which regulates trichome differentiation and anthocyanin biosynthesis in Arabidopsis, encodes a WD40 repeat protein. Plant Cell 11, 1337–1349. Walker, A.R., Lee, E., Bogs, J., Mcdavid, D.A.J., Thomas, M.R., Robinson, S.P., 2007. White grapes arose through the mutation of two similar and adjacent regulatory genes. Plant J. 49, 772–785. Wang, Z., Fang, B., Chen, J., Zhang, X., Luo, Z., Huang, L., Chen, X., Li, Y., 2010. De novo assembly and characterization of root transcriptome using Illumina paired-end sequencing and development of cSSR markers in sweetpotato (Ipomoea batatas). BMC Genomics 11, 726. Wang, Z., Gerstein, M., Snyder, M., 2009. RNA-Seq: a revolutionary tool for transcriptomics. Nat. Rev. Genet. 10, 57–63. Winkel-Shirley, B., 2001. Flavonoid biosynthesis. A colorful model for genetics, biochemistry, cell biology, and biotechnology. Plant Physiol. 126, 485–493. Xie, C., Mao, X., Huang, J., Ding, Y., Wu, J., Dong, S., Kong, L., Gao, G., Li, C.Y., Wei, L., 2011. KOBAS 2.0: a web server for annotation and identification of enriched pathways and diseases. Nucleic Acids Res. 39, W316–W322. Xie, D.Y., Sharma, S.B., Paiva, N.L., Ferreira, D., Dixon, R.A., 2003. Role of anthocyanidin reductase, encoded by BANYULS in plant flavonoid biosynthesis. Science 299, 396–399. Young, M.D., Wakefield, M.J., Smyth, G.K., Oshlack, A., 2010. Gene ontology analysis for RNA-seq: accounting for selection bias. Genome Biol. 11, R14. Zhang, W., Li, X., Zheng, J., Wang, G., Sun, C., Ferguson, I.B., Chen, K., 2008. Bioactive components and antioxidant capacity of Chinese bayberry (Myrica rubra Sieb. and Zucc.) fruit in relation to fruit maturity and postharvest storage. Eur. Food Res. Technol. 227, 1091–1097.

10, 178–182. Chen, J., Wu, X.T., Xu, Y.Q., Zhong, Y., Li, Y.X., Chen, J.K., Li, X., Nan, P., 2015. Global transcriptome analysis profiles metabolic pathways in traditional herb Astragalus membranaceus Bge. var. mongolicus (Bge.) Hsiao. BMC Genomics 16. Cos, P., De Bruyne, T., Hermans, N., Apers, S., Berghe, D.V., Vlietinck, A.J., 2004. Proanthocyanidins in health care: current and new trends. Curr. Med. Chem. 11, 1345–1359. Dubos, C., Stracke, R., Grotewold, E., Weisshaar, B., Martin, C., Lepiniec, L., 2010. MYB transcription factors in Arabidopsis. Trends Plant Sci. 15, 573–581. Espley, R.V., Hellens, R.P., Putterill, J., Stevenson, D.E., Kutty-Amma, S., Allan, A.C., 2007. Red colouration in apple fruit is due to the activity of the MYB transcription factor, MdMYB10. Plant J. 49, 414–427. Feller, A., Machemer, K., Braun, E.L., Grotewold, E., 2011. Evolutionary and comparative analysis of MYB and bHLH plant transcription factors. Plant J. 66, 94–116. Feng, C., Chen, M., Xu, C.J., Bai, L., Yin, X.R., Li, X., Allan, A.C., Ferguson, I.B., Chen, K.S., 2012. Transcriptomic analysis of Chinese bayberry (Myrica rubra) fruit development and ripening using RNA-Seq. BMC Genomics 13, 19. Gesell, A., Yoshida, K., Tran, L.T., Constabel, C.P., 2014. Characterization of an apple TT2-type R2R3 MYB transcription factor functionally similar to the poplar proanthocyanidin regulator PtMYB134. Planta 240, 497–511. Gou, J.Y., Felippes, F.F., Liu, C.J., Weigel, D., Wang, J.W., 2011. Negative regulation of anthocyanin biosynthesis in Arabidopsis by a miR156-targeted SPL transcription factor. Plant Cell 23, 1512–1522. Grabherr, M.G., Haas, B.J., Yassour, M., Levin, J.Z., Thompson, D.A., Amit, I., Adiconis, X., Fan, L., Raychowdhury, R., Zeng, Q., Chen, Z., Mauceli, E., Hacohen, N., Gnirke, A., Rhind, N., Di Palma, F., Birren, B.W., Nusbaum, C., Lindblad-Toh, K., Friedman, N., Regev, A., 2011. Full-length transcriptome assembly from RNA-Seq data without a reference genome. Nat. Biotechnol. 29 644-U130. Hao, D.C., Ge, G., Xiao, P., Zhang, Y., Yang, L., 2011. The first insight into the tissue specific Taxus transcriptome via illumina second generation sequencing. PLoS One 6, e21110. Heim, M.A., Jakoby, M., Werber, M., Martin, C., Weisshaar, B., Bailey, P.C., 2003. The basic Helix-Loop-Helix transcription factor family in plants: a genome-wide study of protein structure and functional diversity. Mol. Biol. Evol. 20, 735–747. Hichri, I., Heppel, S.C., Pillet, J., Leon, C., Czemmel, S., Delrot, S., Lauvergeat, V., Bogs, J., 2010. The basic Helix-Loop-Helix transcription factor MYC1 is involved in the regulation of the flavonoid biosynthesis pathway in grapevine. Mol. Plant 3, 509–523. Hichri, I., Barrieu, F., Bogs, J., Kappel, C., Delrot, S., Lauvergeat, V., 2011. Recent advances in the transcriptional regulation of the flavonoid biosynthetic pathway. J. Exp. Bot. 62, 2465–2483. Hu, D.G., Sun, C.H., Ma, Q.J., You, C.X., Cheng, L., Hao, Y.J., 2016. MdMYB1 regulates anthocyanin and malate accumulation by directly facilitating their transport into vacuoles in apples. Plant Physiol. 170, 1315–1330. Huang, Y.F., Vialet, S., Guiraud, J.L., Torregrosa, L., Bertrand, Y., Cheynier, V., This, P., Terrier, N., 2014. A negative MYB regulator of proanthocyanidin accumulation, identified through expression quantitative locus mapping in the grape berry. New Phytol. 201, 795–809. Jones, S., 2004. An overview of the basic helix-loop-helix proteins. Genome Biol. 5, 226. Kalra, S., Puniya, B.L., Kulshreshtha, D., Kumar, S., Kaur, J., Ramachandran, S., Singh, K., 2013. De novo transcriptome sequencing reveals important molecular networks and metabolic pathways of the plant, Chlorophytum borivilianum. PLoS One 8, e83336. Kobayashi, S., Ishimaru, M., Hiraoka, K., Honda, C., 2002. Myb-related genes of the Kyoho grape (Vitis labruscana) regulate anthocyanin biosynthesis. Planta 215, 924–933. Lister, R., Gregory, B.D., Ecker, J.R., 2009. Next is now: new technologies for sequencing of genomes, transcriptomes, and beyond. Curr. Opin. Plant Biol. 12, 107–118. Liu, X., Feng, C., Zhang, M., Yin, X., Xu, C., Chen, K., 2013. The MrWD40-1 gene of Chinese bayberry (Myrica rubra) interacts with MYB and bHLH to enhance anthocyanin accumulation. Plant Mol. Biol. Rep. 31, 1474–1484. Livak, K.J., Schmittgen, T.D., 2001. Analysis of relative gene expression data using realtime quantitative PCR and the 2(T)(-Delta Delta C) method. Methods 25, 402–408. Luo, P., Ning, G., Wang, Z., Shen, Y., Jin, H., Li, P., Huang, S., Zhao, J., Bao, M., 2016. Disequilibrium of flavonol synthase and dihydroflavonol-4-reductase expression associated tightly to white vs. red color flower formation in plants. Front. Plant Sci. 6, 1257. Matus, J.T., Poupin, M.J., Canon, P., Bordeu, E., Alcalde, J.A., Arce-Johnson, P., 2010. Isolation of WDR and bHLH genes related to flavonoid synthesis in grapevine (Vitis vinifera L.). Plant Mol. Biol. 72, 607–620. Morita, Y., Saitoh, M., Hoshino, A., Nitasaka, E., Iida, S., 2006. Isolation of cDNAs for R2R3-MYB, bHLH and WDR transcriptional regulators and identification of c and ca mutations conferring white flowers in the Japanese morning glory. Plant Cell Physiol. 47, 457–470. Mortazavi, A., Williams, B.A., Mccue, K., Schaeffer, L., Wold, B., 2008. Mapping and quantifying mammalian transcriptomes by RNA-Seq. Nat. Methods 5, 621–628.

20