Transcriptome analysis of genes involved in anthocyanins biosynthesis and transport in berries of black and white spine grapes (Vitis davidii)

Background The color of berry skin is an important economic trait for grape and is essentially determined by the components and content of anthocyanins. The fruit color of Chinese wild grapes is generally black, and the profile of anthocyanins in Chinese wild grapes is significantly different from that of Vitis vinifera. However, V. davidii is the only species that possesses white berry varieties among Chinese wild grape species. Thus, we performed a transcriptomic analysis to compare the difference of transcriptional level in black and white V. davidii, in order to find some key genes that are related to anthocyanins accumulation in V. davidii. Results The results of anthocyanins detection revealed that 3,5-O-diglucoside anthocyanins is the predominant anthocyanins in V. davidii. It showed obvious differences from V. vinifera in the profile of the composition of anthocyanins. The transcriptome sequencing by Illumina mRNA-Seq technology generated an average of 57 million 100-base pair clean reads from each sample. Differential gene expression analysis revealed thousands of differential expression genes (DEGs) in the pairwise comparison of different fruit developmental stages between and within black and white V. davidii. After the analysis of functional category enrichment and differential expression patterns of DEGs, 46 genes were selected as the candidate genes. Some genes have been reported as being related to anthocyanins accumulation, and some genes were newly found in our study as probably being related to anthocyanins accumulation. We inferred that 3AT (VIT_03s0017g00870) played an important role in anthocyanin acylation, GST4 (VIT_04s0079g00690) and AM2 (VIT_16s0050g00910) played important roles in anthocyanins transport in V. davidii. The expression of some selected DEGs was further confirmed by quantitative real-time PCR (qRT-PCR). Conclusions The present study investigated the transcriptomic profiles of berry skin from black and white spine grapes at three fruit developmental stages by Illumina mRNA-Seq technology. It revealed the variety specificity of anthocyanins accumulation in V. davidi at the transcriptional level. The data reported here will provide a valuable resource for understanding anthocyanins accumulation in grapes, especially in V. davidii. Electronic supplementary material The online version of this article (doi:10.1186/s41065-016-0021-1) contains supplementary material, which is available to authorized users.


Background
Grapevine (Vitis L.) is an important economic crop used as table fruit, dried raisins, and for wines or juice. The color of berry skin is an important economic trait for grapes. And it is closely related to the components and content of anthocyanins [1,2]. In recent years, the profile of anthocyanins in grape berry skin has been widely studied. The anthocyanins are glycosides and acylglycosides of anthocyanidins. The main anthocyanidins found in grapes are pelargonidin, cyanidin, delphinidin, peonidin, petunidin as well as malvidin which is usually the predominant anthocyanidin in most red grapes [3,4]. However, the kind of anthocyanins exhibits obvious differences in different grape species [5]. In the V. vinifera cultivars, only 3-O-monoglucoside anthocyanins are detected, and just a few cultivars produce pelargonidinbased anthocyanins [6,7]. In non-V. vinifera species, 3,5-O-disglucoside anthocyanins widely exist, even the pelargonidin-based anthocyanins are also detected [7,8]. In addition, most V. vinifera cultivars possess acylated forms of the anthocyanins, but some other grape species do not produce acylated forms of anthocyanins, such as V. rotundifolia [9].
Anthocyanins in grape berry skin are synthesized via the flavonoid pathways, which has been extensively studied. In fact, the synthesis of anthocyanins shares the same upstream pathways with proanthocyanidins and flavonol derivatives [10]. Phenylalanine is the precursor of flavonoid, which is used as substrate, phenylalanine ammonia lyase (PAL), cinnamate-4-hydroxylase (C4H) and 4-coumaroyl-CoA synthase (4CL) catalyze a series of reactions to produce 4-coumaroyl-CoA. The catalysis of chalcone synthase (CHS) is the first committed step in the flavonoid biosynthetic pathway, which can catalyze the synthesis of chalcones [10]. Subsequently, after the action of chalcone isomerase (CHI), the basic three rings of the general C6-C3-C6 flavonoid skeleton is produced [11]. The B ring of the naringenin flavanone can be further hydroxylated by flavonoid -3′-hydroxylase(F3′H) or flavonoid-3′5′-hydroxylase (F3′5′H) to form eriodictyol or pentahydroxyflavanone [12,13]. The naringenin, eriodictyol and pentahydroxyflavanone can be modified by the catalysis of flavanone-3β-hydroxylase (F3H) to form the corresponding dihydrokaempferol, dihydroquercetin and dihydromyricetin, respectively. Besides, the dihydrokaempferol can also be catalyzed by F3′H or F3′5′H to produce other two dihydroflavonols, dihydroquercetin or dihydromyricetin [1]. Then, dihydroflavonol-4-reductase (DFR) catalyzes these dihydroflavonols to form their corresponding leucoanthocyanidins [1]. In the past, leucoanthocyanidin dioxygenase (LDOX) used to be considered as the first key enzyme that could catalyze the formation of anthocyanidins [14,15] and lead the flavonoid flux into the anthocyanin branch. However, more studies showed that LDOX also play an important role in the biosynthesis of proanthocyanidins [16][17][18][19]. And the anthocyanidins which are the production of LDOX can also be catalyzed by anthocyanidin reductase (ANR) to produce the substrates for the proanthocyanidins synthesis [20,21]. The glycosylation, methylation and acylation of anthocyanidins are very important for their stabilization. Glycosylation of the anthocyanidins is the key step to produce anthocyanins which is catalyzed by the UDP-glucose: anthocyanidin: flavonoid glucosyltransferase (UFGT) [1,[22][23][24][25]. In V. vinifera, the anthocyanidins can only be catalyzed by UFGT to glycosylate at C3 position [24]. Thus, it is also called 3GT. The S-adenosyl-L-methionine (SAM) or Omethyltransferase (OMT) can catalyze the methylation of the hydroxyl groups at the C3 positions or both at the C3 and C5 positions on the B rings of the anthocyanins [26,27]. Acylation can greatly enhance the structural diversity and stability of anthocyanins [28]. It is catalyzed by the action of anthocyanin acyltransferases (AAT) [3,29]. There are mainly two types of ACTs that are classified based on the acyl group donors: the BAHD family using acyl-CoA and the serine carboxypeptidase-like(SCPL) group using acyl-activated sugars [30,31]. Anthocyanins are synthesized in the cytoplasm but accumulate in the vacuoles. Recent researches showed that glutathione Stransferase (GST), multidrug resistance-associated protein (MRP) and multidrug and toxic compound extrusion (MATE) were closely related to the transport of anthocyanins [32].
As part of the flavonoid pathway, the synthesis of anthocyanins are regulated by a complex regulation at the transcriptional level [4,10,15,[33][34][35]. Generally, the flavonoid pathway of anthocyanins biosynthesis is under the control of Myb transcriptional factors, basic helixloop-helix proteins(bHLH) and WD40-like proteins, which also play crucial roles in the regulation of flavonols and proanthocyanidins [10,[17][18][19][36][37][38]. The research showed that the bHLH gene, VvMYC1, was characterized as a component of the transcriptional complex regulating anthocyanins biosynthesis in grapevine [39]. In grapes, a series of R2R3-Myb transcriptional factors are related to the synthesis of anthocyanins [40]. The first Myb transcriptional factor Myba in grape was identified and isolated from V. labrusca hybrids [41], and the results suggested that VlmybA1-1, VlmybA1-2 and VlmybA2 transcriptional factors are involved in the regulation of anthocyanin biosynthesis in the grape by regulating the expression of the UFGT gene [41][42][43]. In red V. vinifera grapes, the functional VvmybA1 can regulate the expression of the UFGT gene to promote the synthesis of anthocyanin. However, in white V. vinifera grapes, a retrotransposon (Gret1) is inserted in the 5′-flanking region of VvmybA1 gene to form a non-functional VvmybA1a gene, resulting in the transcription factor losing its function [44]. Subsequently, two other VvmybA regulator genes, VvmybA2 and VvmybA3 have been cloned and identified [45]. The VvmybA1 and VvmybA2 are very similar, both of them can regulate the accumulation of anthocyanins in the grape berries, and the white-fruited grapes are caused by the mutation of these two similar and adjacent regulatory genes [46,47]. More recent researches suggest that variation in anthocyanins content in grapes is involved with the VvmybA gene cluster [48][49][50].
China is the concentrated distribution area of East Asian Vitis species, which has more than 35 Vitis species [51][52][53]. In recent years, the research of anthocyanins in Chinese wild grapes has been carried out, and the profile of anthocyanins in Chinese wild grapes is significantly different from V. vinifera [5,7,[54][55][56]. Generally, only 3-O-monoglucoside anthocyanins are detected in V. vinifera cultivars, and almost all of the cultivars of V. vinifera are devoid of pelargonidin-based anthocyanins. However, in Chinese wild grapes (such as V. amurensis, V. davidii and V. quinquangularis), the 3,5-O-disglucoside anthocyanins widely exist, and pelargonidin-based anthocyanins are detected in some wild species [5,7,55,56]. The berry color of most Chinese wild grapes is black, except for the white-fruited varieties of V. davidii. The subsequent research found that VvmybA1a gene was detected in white-fruited varieties of V. davidii [57]. The white fruit may be caused by this mutation of mybA1 gene. In order to find the causes that lead to the different profiles of anthocyanins between Chinese wild grapes and V. vinifera, we used the black and white varieties of V. davidii as the materials, and the method of transcriptome analysis to detect the differences at the transcriptional level that might be involved in anthocyanins accumulation of these two varieties.

Anthocyanin composition and content
A total of 24 kinds of anthocyanins were detected at three fruit developmental stages of black and white spine grapes (Table 1). Five categories of anthocyanins were detected: 11 kinds of anthocyanidin diglucosides and 13 kinds of anthocyanidin monoglucosides. And 14 kinds of coumaroylated anthocyanins were detected. The highest content of anthocyanins was detected in the third stage black spine grape (B3, DAF120), it was 1382.127 mg.kg −1 . In black spine grape, the major anthocyanin was Malvidin-3,5-O-diglucoside, which accounted for 87% of the total anthocyanins. Besides, we also detected trace amount of anthocyanins in white V. davidii.

mRNA sequencing
A total of 12 cDNA libraries were constructed from the total RNA of black and white spine grape berry skin in three fruit developmental stages (Fig. 1), two biological replicates were made at each stage. These cDNA libraries were subjected to pair-end reading with the Illumina HiSeq 2000 platform, generating from 56 million to 64 million pair-end raw reads of 100 bp in length, respectively. And the raw data have been submitted to SRA database of NCBI (Accession ID: SRP070860; Link: https:// trace.ncbi.nlm.nih.gov/Traces/sra/?study=SRP070860). After removing the low-quality reads and trimming the adapter sequences, we obtained from 27 million to 31 million clean reads ( Table 2). The subsequent analyses were based on the clean data. Though only two biological replicates were used for RNA-seq, the correlation analysis showed that the R 2 of the two replicates of all the samples were greater than 0.95. It fully illustrated the consistency of two biological replicates and the reliability of the RNAseq results.

Sequence alignment and mapping to the reference genome
We used the French-Italian Public Consortium for Grapevine Genome Characterization, publically accessible version of the complete V. vinifera genome at 12× coverage (ftp://ftp.ensemblgenomes.org/pub/release-23/ plants/fasta/vitis_vinifera/dna/) [58], as the reference genome. After quality control, the reads from black spine grape (64-74%) and white spine grape (63-76%) successfully aligned to the reference genome. Most of the reads from each fruit developmental stage for black and white spine grape aligned to a single position. These uniquely mapped reads account on average for approximately 67 and 71% of the total number of sequenced reads for the black and white spine grape, respectively. The number of reads from each fruit developmental stage for both black and white spine grape mapped to '+' and '-' strand were mostly equal. The number of nonsplice reads was approximately twice to splice reads (Table 3).
Approximately 80% of the reads from each fruit developmental stage for both black and white spine grape mapped to exons, 20% mapped to intergenic regions, and 1% mapped to introns (Fig. 2). The reads mapped to introns resulted from the residual of pre-mRNA and intron retention in the process of alternative splicing. However, the incomplete genome annotation led to the reads mapped to intergenic regions.

Differential expression analysis
We used the RPKM value to show the gene expression level (Table 4). About 45-50% of the genes fall into the RPKM range of 0-1, 9-10% of the genes fall into the RPKM range of 1-3, 18-20% of the genes fall into the RPKM range of 3-15, 14-17% of the gens fall into the RPKM range of 15-60, and approximately 7% of the genes that RPKM value is greater than 60.
There were thousands of differential expression genes in three fruit developmental stages. But, we focused on the genes which were related to the change of berry skin color in fruit development. Thus, we compared the expression of genes in different stages, such as B2vsB1, B3vsB2, B3vsB1, B1vsW1, B2vsW2 and B3vsW3. Then, we used the differential expression genes which were gained from the comparison to make venn diagrams (Fig. 3). In the venn B1(DAF40) B2(DAF80) B3(DAF120) W1(DAF40) W2(DAF80) W3(DAF120)  diagram of B2vsB1, B3vsB2 and B3vsB1, 1440 differential expression genes were found. And in the venn diagram of B1vsW1, B2vsW2 and B3vsW3, it was 1006. Accordingly, the scope of selecting the candidate differential expression genes is effectively narrowed.

Selection of candidate genes Functional enrichment analysis
Gene Ontology (GO, http://www.geneontology.org/) is an international standardized gene function classification system which can be used to classify the function of the predicted genes. GO include three parts: molecular  function, biological process and cellular component. We used the differential expression genes between B2 and W2 to make GO enrichment analysis. Thirty significant terms were found, which belonged to biological process and molecular function (Fig. 4). The color of the berry skin is determined by the pigment, and the anthocyanins are the main pigment in grapes [1,2]. The anthocyanins were synthesized in the endoplasmic reticulum, and transported into vacuole by a serious of anthocyanins transporters [32]. Thus, in these 30 terms, transmembrane transport (GO:0055085) and pigment catabolic process(GO:0046149), which possibly contained the candidate genes, came into our attention. The transmembrane transport term contained 107 up-regulated expression genes and 110 down-regulated expression genes. The pigment catabolic process contained 6 upregulated expression genes and 14 down-regulated expression genes (Additional file 1: Table S1).

Expression pattern analysis
We also used differential expression gene clustering methodology to find the genes that are related to the biosynthesis of anthocyanins in spine grape. The genes that showed same expression pattern were clustered together, and they may obtain the same function or belong to the same biological pathway. We used all of the differential expression genes to draw a heatmap (Fig. 6a). The genes that had the same or similar expression patterns were effectively clustered together. However, so many clusters were produced, and it was difficult to find the clusters that contained the genes closely related to the color of the berry skin. But, the previous studies have suggested that UFGT(VIT_16s0039g02230) [1,24,25] and MybA1(VIT_02s0033g00410) [41] were the key genes to regulate the synthesis of the anthocyanins.

B1
B2 B3 W1 W2 W3  Thus, we selected the genes whose expression patterns were similar to these two genes. Then, we obtained 100 genes and used them to make a heatmap (Fig. 6b). The genes were divided into two groups that contained 68 and 32 genes respectively. And the expression pattern of the group that contained 68 genes was very similar to UFGT and MybA1. Then, we can select the genes that were related to the biosynthesis and transport of anthocyanins.

The candidate genes
After these analyses, we finally obtained 41 differential expression genes, which were considered to be related to the biosynthesis and transport of anthocyanins in the berry skin of spine grape ( Table 6). The expression patterns of these genes were very similar, mainly expressed in B2 and B3 (Fig. 7). In these candidate genes, just one gene (VIT_14s0128g00600) was not annotated. Yet its expression pattern was very similar to the other candidate genes, and we inferred that it was probably related to the anthocyanins accumulation.

Selected candidate UGT genes in white spine grape
Trace amount of anthocyanins were detected in white spine grape (Table 1). Thus, we inferred that some other The venn diagrams of differential expression genes. Venn diagrams indicate the overlap of differential expression genes between B2vsB1, B3vsB1 and B3vsB2 (a), and between B1vsW1, B2vsW2 and B3vsW3 (b) Fig. 4 The most enriched GO terms (B2vsW2). Thirty most enriched GO terms were displayed, included 19 biological processes (green) and 11 molecular function (orange) UGTs are probably related to anthocyanins accumulation in white spine grape. We compared the third stage of white spine grape with the first stage, 62 significant differential expression genes were annotated as UGTs (Additional file 2: Table S2), containing 20 up-regulated and 42 down-regulated genes. From the 20 up-regulated genes, we selected the most 5 significant differential expression genes as the candidate genes ( Table 7). All of these five candidate genes were above 4.5-fold higher in W3 than in W1. Especially, VIT_00s0324g00050 showed highly significant differential expression (P = 3.79E-129).

qRT-PCR analysis
To confirm the results obtained from RNA-seq, relative expression profiles of 30 genes were analyzed by real-time RT-PCR at the three fruit developmental stages in black and white spine grapes. These 30 genes related to the biosynthesis and transport of anthocyanins were chosen for qRT-PCR analysis. The black spine grape showed much higher gene expression of the selected genes than that of the white spine grape. For most of the genes, the qRT-PCR results were consistent with those obtained from the expression profile determined from RNA-Seq data (Fig. 8).

Discussion
Anthocyanin composition and content in black and white V. davidii In V. vinifera, only 3-O-monoglucoside anthocyanins were detected [6]. However, in our research, 3,5-O-diglucoside anthocyanins was the predominant anthocyanins in V. davidii. The studies have shown that, the 3,5-O-diglucoside anthocyanins were usually detected in other Chinese wild grapes [5,7,55,56]. Therefore, the V. davidii was closer to the other Chinese wild grapes from the perspective of the composition of anthocyanins. We inferred that the high 3-5-O-disglucoside anthocyanins concentration in black spine grape was closely related to the expression of VIT_09s0002g06590 (5GT).
Interestingly, we also detected trace amount of anthocyanins in white V. davidii. It was generally accepted that anthocyanins were present only in red grapes, and it was used as a standard to define the difference between red grape and white grape [46,62]. But, a very recent research suggested that some white grapes contained measurable traces of anthocyanins [63]. Previous researches showed that UFGT gene was not expressed in white grape berries, and the expression of UFGT was regulated by two very similar regulatory genes, VvMYBA1 and VvMYBA2, which were not transcribed in white grape berries [46,62]. The present study revealed that UFGT (VIT_16s0039g02230) did not express in white V. davidii. However, why did we detect anthocyanins in white grape berries? A previous study showed that there were as many as 181 putative UDP-glycosyltranferases (UGTs) found in the genome of V. vinifera [64]. Arapitsas et al. [65] detected measurable trace amounts of anthocyanins in some white grape varieties, and they inferred that the other UGTs expressed in berry Statistics of KEGG pathway enrichment. The vertical axis refers to pathways, horizontal axis refers to rich factor. The sizes of the dots represent the numbers of differential expression genes in this pathway. The colors of the dots correspond to the scopes of Q-value skin accepts anthocyanidin as substrate and is therefore involved in the synthesis of trace amount anthocyanins detected in the berry skin of white grape. In our study, 5 candidate UGT genes were above 4.5-fold higher in W3 than in W1, especially VIT_00s0324g00050 showed highly significant differential expression (P = 3.79E-129). These genes are probably related to anthocyanins accumulation in white V. davidii.

The candidate genes related to anthocyanins biosynthesis
The synthesis of anthocyanins via the phenylpropanoid and flavonoid pathway, and it is catalyzed by a serious of enzymes ( Fig. 9) [1,10,37]. However, most of the genes coding for these enzymes are multi-copied in the grape genome [65]. Phenylalanine ammonia lyase (PAL) catalyzes the first step of the phenylpropanoid pathway. In our study, two differential expression genes (VIT_06s 0004g02620 and VIT_13s0019g04460) were annotated as PAL. VIT_06s0004g02620 was mainly expressed in B2 and B3, but just had a trace amount expression in W3. VIT_13s0019g04460 was also mainly expressed in B2 and B3, and had a high expression in W3. Thus, we inferred that VIT_06s0004g02620 played a major role in the pathway. Under the action of C4H and 4CL, 4-Coumaroyl-CoA was produced. We did not find any obvious differential expression genes annotated as C4H.  However, VIT_16s0039g02040 and VIT_16s0050g00390 were annotated as 4CL. Chalcone synthase (CHS) is the first committed enzyme of the flavonoid pathway. Previous studies suggested that 3 CHS genes (CHS1, CHS2 and CHS3) played a role in this pathway. Of the CHSs, the mRNAs of CHS1and CHS2 were detected in both leaves and berry skins of white and red grape cultivars, whereas the mRNA of CHS3 was mainly accumulated in the berry skin of red cultivars during coloration [65][66][67].
In our study, only CHS2 (VIT_14s0068g00920) and CHS3 (VIT_05s0136g00260) showed significant expression differences between black and white spine grapes. Especially, CHS3 (VIT_05s0136g00260) was detected as having a very high expressional level in B2 and B3, and significantly higher than that in white spine grape. It was consistent with the previous studies [65][66][67]. Two CHI genes (VIT_13s0067g03820 and VIT_13s0067g02870) showed obvious differential expression levels between black and white spine grapes. These two CHI genes (CHI1 and CHI2) have been reported by previous researches [65,67,68]. The researches suggested that two F3H genes (F3H1 and F3H2) played a role in the flavonoid pathway, and the mRNA of F3H2 was detected at a high level in grape berry skin during coloration [65,67]. However, in our study, only the expression of F3H1 (VIT_04s0023g03370) in black spine grape showed significantly higher than white spine grape. Maybe F3Hs play different roles in the flavonoid pathway between V. vinifera and V. davidii. F3′H and F3′5′H belong to the cytochrome P450 super family, and catalyze hydroxylation at the 3′ and 3′,5′ positions of the B-ring of the flavonoid to produce the precursors for cyanidin-based anthocyanins and delphinidin-based anthocyanins, respectively. Thus, components of anthocyanins in grape berry skins are closely related to the expression of F3′H and F3′5′H [62,69]. The cDNAs of F3′H and F3′5′H were first isolated from petunia [70,71]. The previous studies showed that grapevines contain two copies of F3′ H and sixteen copies of F3′5′H, both of the F3′Hs are located on chr17, fifteen of F3′5′Hs located on chr6, and one of F3′5′Hs located on chr8 [72]. The previous study showed that, the mRNA levels of F3′H and F3′5′H were high in grape berry skins at the harvest stage [73]. In our study, two (VIT_11s0016g01020 and VIT_17s00 00g07200) and four (VIT_06s0009g02830, VIT_06s00 09g03010, VIT_06s0009g02810 and VIT_06s0009g02920) differential expression genes were annotated as F3′H and F3′5′H, respectively. Of the two F3′Hs, one F3′H gene (VIT_17s0000g07200) was located on chr17. However, the other F3′H gene (VIT_11s0016g01020) was located on chr11. In addition, VIT_11s0016g01020 showed more significantly differential expression level than VIT_17s00 00g07200 between black and white grapes. We suggested that VIT_11s0016g01020 (F3′H) played a major role in anthocyanins accumulation in black spine grape. All of the four F3′5′H genes showed high expression level in B2 and B3, especially VIT_06s0009g02830 and VIT_06s0009 g03010. But F3′5′H genes were almost not expressed at any of the three stages in white spine grape. Besides, we could see that the expression level of F3′5′H was obviously higher than F3′H in black spine grape. This led to more pentahydroxy-flavone and leucodelphinidin production, which were substrates of malvidin-based anthocyanins. Thus, we detected malvidin-based anthocyanins as the predominant kind of anthocyanins in black spine grape. The dihydroflavonols were catalyzed by DFR and LDOX to produce anthocyanidin. We found the differential expression genes (VIT_18s0001g12800 and VIT_02s0025g04720) Fig. 6 The heatmap of the differential expression genes. a The heatmap of all the differential expression genes. b The heatmap of the selected 100 differential expression genes. Using log 10 (RPKM + 1) value to cluster, the red color referred to high expressed genes, the blue color refer to low expressed genes. From blue to red, refer to the log 10 (RPKM + 1) value gradually rised up. The vertical axis refer to the genes showed similar expressed pattern were clustered in a group

Regulation of anthocyanins synthesis
It has been suggested that there was a mybA gene cluster on chromosome 2 of the grape, which contained MYBA1, MYBA2 and MYBA3 [48]. The expression of MYBA1 and MYBA2 could promote the synthesis of anthocyanins through regulating the expression of UFGT [41,46,47]. In our study, MYBA1 (VIT_02s0033g00410) showed high expression levels in B2 and B3, but it did not express in white spine grape. MYBA2 (VIT_02s00 33g00390) and MYBA3 (VIT_02s0033g00450) showed high expression levels in B2, B3 and W3. Jiao et al. [57] detected Gret1 in white-fruited varieties of V. davidii. Thus, it was reasonable that MYBA1 did not express in white spine grape. MYBA2 could also promote the expression of UFGT [41]. Though MYBA2 showed high expression level in W3, UFGT (VIT_16s0039g02230) almost did not express and just trace amount anthocyanins were detected in white spine grape. We could infer from this observation that MYBA1 may play a major role in Table 6 Expression profiles of anthocyanin biosynthesis and transport candidate genes in spine grape (Continued)   Fig. 7 The heatmap of the 41 candidate genes. The denotations are the same as in Fig. 5. The corresponding gene IDs were noted on the right side of the picture regulating the synthesis of anthocyanins in spine grapes. Besides, another R2R3-MYB gene VIT_01s0011g04760 was annotated as MYB4 expressed both in black and white spine grapes. But it showed significantly differential expression between them. The expression level of VIT_01s0011g04760 in black spine grape was obviously higher than in white spine grape. Thus, it can be regarded as a new candidate gene and needs further studies to confirm whether this gene regulates the biosynthesis of anthocyanins in grapes. The flavonoid biosynthesis pathway is under the control of Myb transcriptional factors, basic helixloop-helix proteins (bHLH) and WD40-like proteins [7, 10, 17-19, 36, 38]. The first bHLH was submitted as VvMYCA1 (accession number EF193002; gene ID VIT_15s0046g02560) [76,77], which was considered to regulate the expression of UFGT. Subsequently, the second bHLH was called VvMYC1 (accession number EU447172; gene ID VIT_07s0104g00090) [39], which was characterized as a component of the transcriptional complex regulating anthocyanin biosynthesis in grapevine. However, both of these two genes did not show obvious differences at the transcriptional level between black and white spine grapes in our study. But, another gene (VIT_14s0060g01010) annotated as bHLH exhibited differential expression between black and white spine grapes. It was shown that WDR1 (Genbank accession number DQ517914) contributed positively to the accumulation of anthocyanins [77]. However, we did not find any WD40 gene having obvious differential expression. In recent years, researches showed that NAC TFs were involved in the regulation of anthocyanins accumulation. A NAC TF has been proposed to be involved in the regulation of anthocyanins accumulation during the response of blood orange to cold exposure [78]. In addition, PpNAC1 can activate the transcription of PpMYB10.1, resulting in anthocyanins pigmentation in bloodfleshed peach [79]. In grape, there was no report that the NAC TFs were related to the anthocyanins accumulation. In our study, a candidate gene (VIT_14s0108g01070) was annotated as NAC, which may be related to anthocyanins accumulation in spine grapes.

Modification of anthocyanins
Anthocyanidins need to be modified by glycosylation, methylation and acylation to form stabilized anthocyanins. In V. vinifera, the glycosylation was catalyzed by UDP-glucose: anthocyanidin: flavonoid glucosyltransferase (UFGT) at C3 position, and the gene encoding UFGT had been cloned [24,80]. The action of UFGT was crucial for anthocyanin accumulation in grape berry skin [1,22,23,25,50,67,81,82]. In our study, a differential expression gene (VIT_16s0039g02230) was annotated as UFGT, which had a high expression level in B2 and B3, and did not express in white spine grape. However, in non-V. vinifera species, the 3-5-O-disglucoside anthocyanins are widely present [7,8]. There must be another UDP-glucose: anthocyanidin: flavonoid glucosyltransferase to catalyze glycosylation at C5 position. In recent years, a few studies reported the isolation of 5GT genes in grapes [83,84]. Jánváry et al. [83] cloned functional Cha5GT and nonfunctional Dia5GT from the heterozygous hybrid cultivar 'Regent' , a cross of V. vinifera cv. 'Diana' and the interspecific hybrid cv. 'Chambourcin'. The functional analysis of Cha5GT and Dia5GT suggested that two mutations in the 5GT gene eliminated its enzymatic activity. Because of the absence of active 5GT, dis-glucosidic anthocyanins could not be produced in V. vinifera red grapes. He et al. [84] cloned the full-length cDNA of UDP-glucose: anthocyanin 5-Oglucosyltransferase (Va5GT) from V. amurensis Rupr. cv. 'Zuoshanyi'. The results suggested that Va5GT was a key enzyme in the biosynthesis of dis-glucosidic anthocyanins in V. amurensis grape berries. In our study, a differential expression gene (VIT_09s0002g06590) showed a high expressional level just in B2 and B3, which was annotated as UGT. It has high homology with Va5GT. Thus, we inferred that the expression of VIT_09s0002g06590 (5GT) led to a high 3-5-O-disglucoside anthocyanins concentration in black spine grape.
In plants, methylated anthocyanidins accounted for a large proportion of the total reported anthocyanidins [3]. Especially, three methylated anthocyanidins, peonidin, petunidin and malvidin, were commonly present in grape berry skin. The methylation was catalyzed by Sadenosyl-L-methionine (SAM) or O-methyltransferase  (OMT) at the C3 positions or both at the C3 and C5 positions on the B rings of the anthocyanins in grape [26,27]. And the cDNAs of several OMT had been cloned in grapes [85,86]. Subsequently, a QTL for anthocyanins methylation variation was identified that was colocalized with a cluster of three putative OMT genes: VIT_01 s0010g03470 (OMT3), VIT_01s0010g03490 (OMT2) and VIT_01s0010g03510 (OMT1) [87,88]. Fournier-Level et al. [87] reported that OMT2 gene presented two SNPs associated with methylation level. It probably led to a structural change of the OMT2 protein, thus significantly affected the enzyme specific catalytic efficiency for the 3′-O-methylation of delphinidin-3-glucoside. In this study, three differential expression genes (VIT_07s0031g00350, VIT_01s0010g03510 and VIT_01s0010g03490) were annotated as OMTs. Of the three genes, VIT_01s0010g03510 (OMT1) and VIT_01s0010g03490 (OMT2) expressions were consistent with the previous researches [87,88], which did not express in white spine grape. It is interesting that OMT1 showed higher activity than OMT2 in black spine grape. In addition, a new candidate gene (VIT_07s0 031g00350) was also annotated as OMT, and showed high (See figure on previous page.) Fig. 8 Quantitative real-time PCR validation of RNA-Seq data. Relative expression profiles of 30 genes showed the expression fold changes (FC) in comparison between the three fruit developmental stages in the black and white spine grapes. Histograms represent expression levels as assessed by RNA-Seq, data are reported as means ± SE of two biological replicates (left axis), the columns with diagonal lines or blank represent the RPKM of black or white spine grape, respectively. The line charts represent expression fold changes as assessed by qRT-PCR, data are reported as means ± SE of three replicates (right axis), the lines with black triangle or '×' mark represent black or white spine grape, respectively Acylation was a common modification of anthocyanins in grapes. It not only increased diversity of anthocyanins, but also improved the color stabilization and intensity for the anthocyanins [28]. Two enzyme families (BAHD-ATs and SCPL-ATs) have been reported as related to the acylation of anthocyanins [30,31]. A recent study identified a number of QTLs associated with variation in acylated anthocyanin levels in F1 progeny from a 'Syrah'x'-Pinot Noir' cross. The strongest candidate genes within these QTLs included those belonging to the BAHD and SCPL acyltransferase family [88]. But, no QTL was found to cause the presence/absence of acylation in berries. In our study, two differential expression genes (VIT_1 4s0068g01440 and VIT_03s0017g00870) were annotated as acyltransferase (AT). Especially VIT_03s0017g00870 did not express in white spine grape and at the first stage of black spine grape. Yet it showed high expressional level in B2 and B3. Simultaneously, we detected acylated anthocyanins in black spine grape. Thus, this candidate gene can be very possibly related with the acylation of anthocyanins. It is known that V. vinifera cv. Pinot Noir does not synthesise acylated anthocyanins [89]. After SNP analysis of the candidate genes, at the position 36 of the putative coding region of VIT_03s0017g00870, A is replaced by G, compared with the PN40024 grapevine genome [58], but it does not alter the predicted amino acid sequence. Rinaldo et al. [90] identified a gene, anthocyanin 3-O-glucoside-6-O-acyltransferase (Vv3AT), which encoded a BAHD acyltransferase protein. This protein can promote the synthesis of acylated anthocyanins, which was transcriptionally regulated by VvMYBA. This is consistent with our study.

Anthocyanins transport
Anthocyanins are synthesized in the endoplasmic reticulum and then transported into the vacuoles. Three kinds of anthocyanin transporters: glutathione Stransferase (GST), multidrug resistance-associated protein (MRP) and multidrug and toxic compound extrusion (MATE), were reported as related to the anthocyanins transport [32,91]. In grapes, a few candidate anthocyanin transporters have been reported.
Conn et al. [92] cloned five GSTs (VvGST1, VvG ST2, VvGST3, VvGST4 and VvGST5) from V. vinifera. cv. Gamay Fréaux and the study showed that VvGST1 and VvGST4 coded for the enzymes which had the function of anthocyanins transport. However, in our study, GST1 (VIT_19s0093g00320) was not detected, GST4 (VIT_04s0 079g00690) and GST5 (VIT_19s0015g02730) showed significantly differential expression between black and white spine grape. Especially, GST4 just express in B2 and B3 (Fig. 10). GST2 (VIT_07s0005g00030) showed the highest expressional level in the first stage, and it possessed higher transcriptional level in white spine grape than black one. GST3 (VIT_12s0028g00920) also showed the highest expressional level at the first stage, but did not show obvious differential expression between black and white spine  Fig. 10 Expressional levels of six GSTs in black and white spine grapes. The vertical axes refer to expression levels as assessed by RNA-Seq, horizontal axes refer to three berry developmental stages. The black columns represent black spine grape, the white columns represent white spine grape grapes. Thus, we inferred that GST4 (VIT_04s0079g00690) and GST5 (VIT_19s0015g02730) play an important role in anthocyanins transport in spine grapes. Besides, another two candidate genes (VIT_19s0015g02690 and VIT_19s 0015g02880) were annotated as GSTs, which showed higher expressional level in black spine grape than in white one (Fig. 10). Hence, these two candidate genes are also probably related to anthocyanins transport in spine grapes.
Gomez et al. [93] identified three multidrug and toxic compound extrusion (MATE) genes as candidate antho-MATEs(AMs). In their study, just AM1 and AM3 were cloned from V. vinifera L. cv. Syrah, AM2 was not successfully cloned. They inferred that AM2 was probably not expressed at detectable levels in mature berry or was a pseudogene. Accordingly, their results revealed that AM1 and AM3 just could transport acylated anthocyanins in the presence of Mg ATP [93]. Yet in our study, AM1 (VIT_16s0050g00900) was hardly expressed in spine grapes, the expressional level of AM3 (VIT_16s005 0g00930) in white spine grape was higher than in black spine grape. However, AM2 (VIT_16s0050g00910) was mainly expressed in B2 and B3 (Fig. 11). Thus, we infer that AM2 plays an important role in spine grape anthocyanins transportation.
The plant ATP binding cassette (ABC) transporters, in particular from the ABCC subfamily (formerly named multidrug resistance proteins [MRPs]) had been reported as related to the accumulation of flavonoid in vacuoles [94]. In grapes, an ABC transporter, ABCC1 was identified, which localizes to the tonoplast and is involved in the transport of glucosylated anthocyanidins [95]. In our study, four candidate genes were annotated as ABC transporters: VIT_16s0050g02480 (ABCC1), VIT_10s0003g 04390 (ABCC2), VIT_09s0002g02430 (ABCC8) and VI T_03s0017g01290 (ABCG11). Francisco et al. [95] were unsuccessful to clone the ABCC2 from the exocarp of grape, and inferred it might not be expressed in detectable amounts during the ripening stage. However, our study showed that both ABCC1 and ABCC2 were expressed in black and white spine grape berry skins. All of the four candidate ABCs were higher in expression in black spine grape than in white spine grape, but the differences were not very significant (Fig. 12). Thus, more studies are needed to confirm the key ABC transporters in grapes.

Conclusions
V. davidii is the only Chinese wild grape species which possesses white berry varieties. High levels of 3,5-Odiglucoside anthocyanins were detected in the black berry skin of V. davidii. The present study investigated the transcriptome profiles of the berry skin from black and white spine grapes at three fruit developmental stages by Illumina mRNA-Seq technology. The examination of absolute expression count for every gene has enabled us to carry out a global investigation of gene expression at these three key time-points in black and white spine grapes. The transcriptome analysis presented thousands of DEGs. We used gene clustering and the enrichment of GO and KEGG to describe the transcriptional patterns of genes involved in anthocyanins accumulation. We found 41 differentially expressed genes probably related to anthocyanins accumulation in V. davidii, including the genes that encode enzymes, transcription factors and transporters involved in anthocyanins biosynthesis, regulation and transport. Some genes were consistent with the previous studies in other grape species, some were newly found in our study. Trace amounts of anthocyanins were detected in berry skins of white V. davidii. Five candidate UGT genes were probably related to the biosynthesis of anthocyanins in white V. davidii. In conclusion, the present study provides new insights into the understanding of anthocyanins accumulation in grapes.

Sample collection
Berries of black and white spine grapes (V. davidii) were used in this study. Both of the black and white spine grapes were collected from the wild in Hunan, China. The plants were grown in China National Germplasm Resources Repository of Grape, Zhengzhou Fruit Research Institute, Chinese Academy of Agricultural Sciences. Skins of grape berries were collected at three different fruit developmental stages: before veraison  Fig. 11 Expressional levels of AM2 and AM3 in black and white spine grapes. The denotations are the same as in Fig. 10 (Veraison refers to the stage that the berries begin to color and be soft) (40 days after flowering [DAF]), at veraison (80 DAF) and fruit ripe period (120 DAF) (Fig. 1). Two biological replicates were made at each stage. At each stage, we harvested 6-8 clusters at each sampling date, and 5 berries were collected from each cluster. Thus, 30-40 berries were used to collect the skin at each sampling date. After collection, samples were flashfrozen in liquid nitrogen and stored at −80°C until further processing.

The extraction and determination of anthocyanins
Extraction of total anthocyanins was as described by He et al. [7] with some modifications. The peels were put into mortars and ground in liquid nitrogen. Aliquots of 0.5 g ground powder were added into 10 ml centrifuge tubes with 8 ml 2% formic acid-methanol solution. After ultrasonic oscillation for 10 min, the extracts were put on a shaker in the dark in the table concentrator at 25°C , 200 rpm for 30 min, followed by centrifugation at 4°C , 12000 rpm for 10 min. The supernatants were transferred into 50 ml centrifuge tubes. The residues were reextracted 3 times. The organic fractions were pooled, evaporated by a vacuum rotary evaporator (BUCH, USA) at 40°C. The residual parts were poured into activated solid phase extraction cartridges. The solid phase extraction cartridges were washed with 5 ml water for 2 times. After removing the leacheate, the solid phase extraction cartridges were eluted with 10 ml methanol for 2 times. The filtrates were collected and evaporated to dryness, then re-dissolved in 5 mL 0.5% hydrochloric acid-methanol solution. Finally, the solutions were filtered through a 0.22 μm Millipore filter for analysis.
The anthocyanins content was determined by an ACQUITY Ultra Performance Liquid Chromatography (UPLC) system (Waters, Milford, MA, USA) linked to both a PDAeλ photo diode array detector (Waters, Milford, MA, USA) and a Micromass Quattro micro TM API benchtop triple quadrupole mass spectrometer (Waters MS Technologies, Manchester, UK), with a electrospray ionization (ESI) source operating in multiple reaction monitoring (MRM) mode. Sample solutions were injected into a ACQUITY UPLC®HSS T3 column (2.1 × 150 mm i.d,with 1.8 μm particle size, Waters, Milford, MA, USA), which was maintained at 30°C. The mass spectrometric acquisition parameters were as follows: ESI source temperature 150°C, desolvation gas temperature 400°C, desolvation gas flow rate 800 L/h, cone gas flow rate 50 L/h, collision gas (high purity argon gas) flow rate 0.14 mL/min. The mobile phase had acetonitrile as solvent A, and 0.5% hydrochloric acid solution as solvent B. The gradient profile began with 5-10% A for 1 min, 10-25% A for 16 min, 25-40% A for 18 min, 40-100% A for 19 min, and then returned to initial conditions for 20 min and for 5 min. The flow rate was 1.0 ml min −1 and the column temperature was set at 40°C. The injection volume was 2.0 μl. The detection wavelength was 520 nm.

Total RNA extraction and qualification
Total RNA was extracted using TIANGEN RNAprep Pure Plant Kit (Tiangen Biotech Beijing, China). RNA degradation and contamination was monitored on 1% agarose gels. RNA purity was checked using the Nano-Photometer® spectrophotometer (IMPLEN, CA, USA). RNA concentration was measured using Qubit® RNA Assay Kit in Qubit® 2.0 Flurometer (Life Technologies,  Fig. 12 The expressional level of the ABCs in black and white spine grapes. The denotations are the same as in Fig. 10 CA, USA). RNA integrity was assessed using the RNA Nano 6000 Assay Kit of the Bioanalyzer 2100 system (Agilent Technologies, CA, USA).
cDNA library construction and transcriptome sequencing A total amount of 3 μg RNA per sample was used as the input material. Sequencing libraries were generated using NEBNext® Ultra™ Directional RNA Library Prep Kit for Illumina® (NEB, USA) following manufacturer's recommendations and index codes were added to attribute sequences to each sample. Briefly, mRNAs were purified from the total RNAs using poly-T oligoattached magnetic beads. Fragmentation was carried out using divalent cations under elevated temperature in NEBNext First Strand Synthesis Reaction Buffer (5X). First strand cDNA was synthesized using random hexamer primer and M-MuLV Reverse Transcriptase (RNa-seH − ). Second strand cDNA synthesis was subsequently performed using DNA Polymerase I and RNase H. In the reaction buffer, dNTPs with dTTP were replaced by dUTP. The remaining overhangs were converted into blunt ends via exonuclease/polymerase activities. After adenylation of 3′ ends of DNA fragments, NEBNext Adaptor with hairpin loop structure were ligated for hybridization. In order to select cDNA fragments of preferentially 150-200 bp in length, the library fragments were purified with AMPure XP system (Beckman Coulter, Beverly, USA). Then 3 μl USER Enzyme (NEB, USA) was used with size-selected, adaptor-ligated cDNAs at 37°C for 15 min followed by 95°C for 5 min before PCR. Then PCR was performed with Phusion High-Fidelity DNA polymerase, universal PCR primers and Index (X) Primer. The PCR products were purified (AMPure XP system) and library quality was assessed on the Agilent Bioanalyzer 2100 system.
The clustering of the index-coded samples was performed on a cBot Cluster Generation System using Tru-Seq PE Cluster Kit v3-cBot-HS (Illumia) according to the manufacturer's instructions. After cluster generation, the library preparations were sequenced by an Illumina Hiseq 2000 platform and 100-base paired-end reads were generated.

Sequencing data analysis
Raw data (raw reads) of fastq format were first processed through in-house PERL scripts. In this step, clean data (clean reads) were obtained by removing reads containing adapter, reads containing ploy-N and low quality reads from raw data. At the same time, Q20, Q30 and GC content of the clean data were calculated. All the downstream analyses were based on the clean data with high quality. Reference genome and gene model annotation files were downloaded from genome website directly (ftp://ftp.ensemblgenomes.org/pub/release-23/plants/ fasta/vitis_vinifera/ dna/). Index of the reference genome was built using Bowtie v2.0.6 and base paired-end clean reads were aligned to the reference genome using TopHat v2.0.9. We selected TopHat as the mapping tool since TopHat can generate a database of splice junctions based on the gene model annotation file and thus a better mapping result than other non-splice mapping tools.
HTSeq v0.6.1 was used to count the reads numbers mapped to each gene. And then RPKM of each gene was calculated based on the length of the gene and reads count mapped to this gene. RPKM (Reads Per Kilobase of exon model per Million mapped reads) was used for the effect of sequencing depth and gene length for the reads count at the same time, and is currently the most commonly used method for estimating gene expression levels [96].

Differential expression analysis
Differential expression analysis of two conditions/groups (two biological replicates per condition) was performed using the DESeq R package (1.10.1). DESeq provided statistical routines for determining differential expression in digital gene expression data using a model based on the negative binomial distribution. The resulting P-values were adjusted using the Benjamini and Hochberg's approach for controlling the false discovery rate. Genes with an adjusted P-value < 0.05 selected by DESeq were assigned as differentially expressed.

Selection of candidate genes
We used Gene Ontology (GO) and KEGG enrichment analysis of differentially expressed genes to select candidate genes. GO enrichment analysis of differentially expressed genes was implemented by the GOseq R package, in which gene length bias was corrected. GO terms with corrected P-value less than 0.05 were considered significantly enriched by differential expressed genes. KEGG is a database resource for understanding highlevel functions and utilities of the biological system, such as the cell, the organism and the ecosystem, from molecular-level information, especially large-scale molecular datasets generated by genome sequencing and other high-through put experimental technologies (http://www.genome.jp/kegg/). We used KOBAS software to test the statistical enrichment of differential expression genes in KEGG pathways. We also used the differentially expressed genes clustering methodology to find the candidate genes. The clustering analysis was made based on the expression pattern of the differentially expressed genes. We selected the genes that the expression patterns were very similar to the genes closely related to the anthocyanins biosynthesis as the candidate genes.