- Research article
- Open Access
Transcriptome profiling of ‘Kyoho’ grape at different stages of berry development following 5-azaC treatment
BMC Genomics volume 20, Article number: 825 (2019)
5-Azacytidine (5-azaC) promotes the development of ‘Kyoho’ grape berry but the associated changes in gene expression have not been reported. In this study, we performed transcriptome analysis of grape berry at five developmental stages after 5-azaC treatment to elucidate the gene expression networks controlling berry ripening.
The expression patterns of most genes across the time series were similar between the 5-azaC treatment and control groups. The number of differentially expressed genes (DEGs) at a given developmental stage ranged from 9 (A3_C3) to 690 (A5_C5). The results indicated that 5-azaC treatment had not very great influences on the expressions of most genes. Functional annotation of the DEGs revealed that they were mainly related to fruit softening, photosynthesis, protein phosphorylation, and heat stress. Eight modules showed high correlation with specific developmental stages and hub genes such as PEROXIDASE 4, CAFFEIC ACID 3-O-METHYLTRANSFERASE 1, and HISTONE-LYSINE N-METHYLTRANSFERASE EZA1 were identified by weighted gene correlation network analysis.
5-AzaC treatment alters the transcriptional profile of grape berry at different stages of development, which may involve changes in DNA methylation.
Grape (Vitis vinifera L.) is one of the most important perennial woody fruit crops in the world. The grape berry is consumed whole or in the form of raisins or wine and has high nutritional, medicinal, and economic value , making it one of the most popular fruits. Grape berry exhibits change in pigmentation, sugar and organic acid contents, and other quality components during development and ripening  and is a useful model for studying fruit development.
Transcriptome sequencing is the main technology for investigating genome-wide changes in gene expression patterns, and has been used to study metabolic pathways and gene expression during fruit development in many plants. Most of the research has focused on climacteric fruits such as bayberry , pear [4, 5], kiwifruit , peach , tomato , and apricot , although recent studies have also investigated non-climacteric fruits such as sweet orange  and strawberry . For example, cell wall biosynthesis, carbohydrate metabolism, the tricarboxylic acid cycle, and carotenoid biosynthesis were shown to be differentially regulated during fruit development and ripening of the sweet orange variety ‘Anliu’ and its red-fleshed mutant ‘Hong Anliu’ . Metabolic shifts occurred in the green-white-red stages of strawberry that were associated with differential gene expression, and it was found that oxidative phosphorylation plays an important role in the regulation of fruit maturation .
Whole-genome sequencing of the PN40024 genotype of grapevine, originally derived from Pinot Noir, was completed in 2007 and has provided a useful resource for functional genomic studies . A transcriptome analysis revealed that reduced biosynthesis, photosynthesis, and transport was the main reason for delayed senescence of the peel . Most genes showed comparable expression levels between ‘Kyoho’ berry and its early-ripening mutant ‘Fengzao’ , and an analysis of differentially expressed genes (DEGs) revealed that those related to oxidative stress genes likely promote the early ripening of ‘Fengzao’ grape berry. Genes involved in carbohydrate metabolism and regulation of flavonoid metabolism and those of the solute carrier family showed the most marked changes in expression in ‘Kyoho’ and transgenic berry peels , and it was later reported that V. vinifera VACUOLAR H+-PPASE 1 was activated by the MYB transcription factor MYBA1 and that hexokinase-mediated glucose signaling increased the expression of anthocyanin biosynthesis and transport-related genes to promote anthocyanin accumulation in grape peel. In addition, differences in the levels of microRNAs (miR169-NF-Y subunit, miR398-CSD, miR3626-RNA helicase, miR399-phosphate transporter, and miR477-GRAS transcription factor) and their targets have been observed in ‘Kyoho’ and ‘Fengzao’ during berry development and ripening .
DNA methylation is a mitotically reversible and meiotically heritable epigenetic modification  that is important in plant growth and development [18,19,20]. Recent studies have shown that DNA methylation is associated with fruit development and ripening [21,22,23,24,25,26]. Treatment with 5-azacytidine (5-azaC), a methyltransferase inhibitor, was shown to affect the development of tomato , strawberry , and Acca sellowiana  fruit by decreasing DNA methylation levels, resulting in an early ripening phenotype. Although 5-azaC treatment delayed fruit ripening in sweet orange , it also had a genome-wide demethylating effect . 5-AzaC promoted the early ripening of grape berry and reduced global methylation level at a concentration of 100 μΜ in our previous study . However, the mechanism by which 5-azaC alters gene expressions to accelerate berry ripening remains unknown.
To answer this question, in this study we carried out RNA-sequencing (RNA-seq) analysis of ‘Kyoho’ grape berry at five different stages of fruit development after 5-azaC treatment. The results provide novel insight into the molecular basis of grape berry ripening and a basis for future molecular studies.
Analysis of RNA-seq libraries
To identify the genes involved in grape berry development, we performed transcriptome sequencing of ‘Kyoho’ grape berry with or without 5-azaC treatment at different developmental stages. The RNA-seq data have been uploaded to the National Center for Biotechnology Information Sequence Read Archive under the accession number PRJNA542248. A total of 30 cDNA libraries were constructed comprising 1.37 billion raw reads; 1.33 billion clean reads (accounting for 96.74% of raw reads) were recorded after removing adapter sequences and reads of low quality and those with more than 5% N bases. The average number of clean reads per sample was about 45.76 million and the clean Q30 (sequencing error rate < 0.1%) base rate was > 93.6% for each sample. Ultimately, 1.21 billion high-quality reads (accounting for 91.32% of clean reads) were mapped to the grape reference genome; 29.56 million of these were mapped to multiple locations in the genome at a ratio of 2.23% (Additional file 1).
In the 5-azaC-treated and untreated control samples, more genes were expressed at the A3 (23883) stage than at the C3 (22710) stage, whereas fewer genes were expressed at the other four stages. We also analyzed the number of genes expressed at different levels (fragments per kilobase million [FPKM] ≥ 50, 50 > FPKM ≥10, 10 > FPKM ≥2, 2 > FPKM ≥0.1, FPKM < 0.1) and found that the number of genes with FPKM ≥10 was higher in berries at A2 and A3 stages than in berries at stages C2 and C3; the number of genes with different expression levels was greater at C2 than at A2 (Additional file 2).
Gene expression profile following 5-azaC treatment
Principal component analysis (PCA) revealed the heterogeneity of grape samples at different developmental stages based on gene expression in all samples. Dim1 and Dim2 had values of 23.3 and 18.8%, respectively, and accounted for 42.1% of the principal components (Fig. 1). PCA also revealed a consistency (i.e., no differences) between the three replicates at each developmental stage. Samples in the treatment and control groups at the same developmental stage clustered together, reflecting a lack of difference between them. On the other hand, samples at different developmental stages were dispersed irrespective of treatment condition, indicating that they differed significantly.
We carried out a Gene Ontology (GO) enrichment analysis in order identify the biological processes in berry development that were affected by 5-azaC treatment and identified 11 enriched GO terms including those related to Zinc ion binding, Pyrophosphatase activity, Nucleoside triphosphatase activity, Nuclease activity, Hydrolase activity, and Endonuclease activity (Fig. 2). The number and significance of genes related to Nuclease activity, Endonuclease activity, and Isomerase activity were similar for control and treatment groups. Pyrophosphatase activity and Hydrolase activity (acting on acid anhydrides) were significantly enriched after 5-azaC treatment at stages C1, C2, and C3 (with the highest fold enrichment at C3) as well as at stage A4. Meanwhile, Nucleoside triphosphatase activity was enriched at C2, C3, and A4 (Fig. 2).
Comparison of overall expression patterns by time course sequencing (TCseq) analysis
To determine the overall expression patterns of genes common to the treatment and control groups, categories with different expression patterns were identified by TCseq analysis. Genes with similar expression patterns clustered together, with the highest Calinski criterion value occurring in eight clusters, suggesting that this was the optimal number of clusters (Fig. 3). The gene numbers for clusters 1–8 ranged from 1414 (cluster 8) to 6213 (cluster 5). Genes in each cluster showed very similar expression patterns overall in the treatment and control groups, whereas those in different clusters showed distinct expression patterns (Fig. 4).
The genes in clusters 1–8 showing similar expression patterns between treatment and control groups (Additional file 3) were divided into the following four classes. The expression of genes in cluster 2/4/6 first increased and then decreased with berry development; gene expression in cluster 3/5 gradually decreased before reaching a plateau; genes in cluster 7 showed relatively stable expression levels at the early stage of berry development followed by gradual upregulation; and the level of genes in cluster 8 remained constant across developmental stages. Cluster 1 comprising 4261 genes was exceptional; gene expression at C1 decreased gradually with berry development, but the level at A1 was lower than at C1. GO enrichment analysis of genes in cluster 1 revealed significant grouping of 16 GO terms (Fig. 5) including Structural molecule activity, Structural constituent of ribosome, Pyrophosphatase activity, Nucleoside triphosphatase activity, Hydrolase activity, Protein heterodimerization activity, and Tubulin binding (Additional file 4).
Analysis of differentially expressed genes (DEGs)
We compared the transcriptional profiles of the treatment and control groups at the various stages of berry development and identified DEGs at each stage except for A1_C1. The number of DEGs between treatment and control groups at a given stage varied from 9 (A3_C3) to 690 (A5_C5). The expression of all 11 DEGs in A2_C2 was decreased (Table 1). The number of DEGs between successive developmental stages was 605 and 2188 for the control group and 104 and 2929 for the treatment group. There were fewer DEGs in A1_A2, A2_A3, and A3_A4 compared to C1_C2, C2_C3, and C3_C4, and the number of DEGs (up- and downregulated) was greater in A4_A5 than in C4_C5. The number of DEGs between successive developmental stages within treatment and control groups was greater than that between treatment and control groups at the same stage (Table 1).
There were very few DEGs at the early stage of development; notably, no DEGs were identified in A1_C1. All DEGs (11) identified in A2_C2 were significantly downregulated after 5-azaC treatment. Functional annotation showed that these genes were Non-specific lipid-transfer protein A (VIT_12s0028g01180, LTA), Endochitinase (VIT_00s1290g00010), Purple acid phosphatase 15 (VIT_05s0029g00200, PAP15), Probable xyloglucan endotransglucosylase/hydrolase protein 23 (VIT_11s0052g01230, XTH25), and Basic 7S globulin (VIT_14s0128g00200). One of the genes, XTH25, is involved in the xyloglucosyl transferase pathway (Table 2).
Nine DEGs in A3_C3 were annotated; two of these—encoding Probable pectinesterase/pectinesterase inhibitor 36 (VIT_15s0048g00500, PE) and Endonuclease 1 (VIT_00s0301g00100)—were downregulated after 5-azaC treatment whereas seven genes encoding UDP-glycosyltransferase 89B2 (VIT_17s0000g04750, UGT89B2), Probable galacturonosyltransferase-like 1 (VIT_18s0001g11860, GTL), Histidine kinase CKl1 (VIT_07s0005g01380), Beta-glucosidase 13 (VIT_13s0064g01760, BGLU13), and Probable sarcosine oxidase (VIT_04s0069g00860) were upregulated. BGLU13 and VIT_04s0069g00860 are involved in the β-glucosidase and sarcosine oxidase/l-pipecolate oxidase pathways, respectively (Table 2).
A total of 19 DEGs were identified in A4_C4: Probable sucrose-phosphate synthase 1 (VIT_04s0008g05730), polygalacturonase (PG) QRT3 (VIT_01s0011g01300, QRT3), and Probable flavin-containing monooxygenase 1 (VIT_18s0122g01430) were upregulated; VIT_18s0122g01430, is involved in the flavin monooxygenase and dimethylaniline monooxygenase (NO-forming) pathways. Downregulated genes were Non-specific lipid-transfer protein 2 (VIT_14s0006g02570), Acidic endochitinase (VIT_15s0046g01570), Pleiotropic drug resistance protein 2 (VIT_13s0074g00700), GDSL esterase/lipase 1 (VIT_09s0002g00550), Cationic peroxidase 1 (VIT_18s0001g06840), Glutathione S-transferase (VIT_07s0005g00030), and 23.6-kDa heat shock protein (VIT_16s0022g00510, HSP 23.6). The HSP23.6 gene belongs to the HSP20 family (Table 2). A5_C5 had the most DEGs. The expression levels of 28 DEGs encoding heat shock proteins and belonging to HSP20, HSP70, HSP90, and HSF_DNA-binding gene families were downregulated after 5-azaC treatment, as were all 28 DEGs related to photosynthesis and some methyltransferase genes including VIT_04s0023g02290, VIT_05s0049g01650, VIT_12s0028g02370, and VIT_08s0007g08470 (Additional file 5).
We performed a Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis of the DEGs in the treatment and control groups at the same developmental stage and found that only DEGs in A5_C5 were significantly enriched in KEGG pathways—namely, Protein processing in endoplasmic reticulum, Photosynthesis, Photosynthesis antenna proteins, Galactose metabolism, Flavone and flavanol biosynthesis, Diterpenoid biosynthesis, and ABC transporters; most genes were involved in Protein processing in endoplasmic reticulum (Additional file 6) and DEGs related to Photosynthesis were downregulated. The expression patterns and details of representative genes in key pathways are shown in Additional file 5.
Weighted gene correlation network analysis (WGCNA)
To gain insight into gene regulatory networks involved in the development of grape berry, we carried out a WGCNA of the transcriptome data of the five developmental stages of grape berry with or without 5-azaC treatment (19,387 genes, FPKM ≥0.5). In the initial module division, we set a soft threshold of 1 and used dynamic pruning to combine modules with a high similarity of characteristic genes (Fig. 6a). We obtained 20 gene modules with similar expression patterns; the total number of genes in each module ranged from 38 (palevioletred) to 6017 (darkolivegreen4). We analyzed the correlation between the characteristic genes of each module and berry development stage by calculating the Pearson correlation coefficient. The WGCNA results revealed that eight of the 20 modules were significantly correlated with a specific developmental stage (P < 0.05)—i.e., lightblue3, darkolivegreen1, powderblue, palevioletred, lightcyan, royalblue, blue3, and deeppink1 were significantly correlated with stages A1, A2, A4, A5, C2, C3, C4, and C5, respectively (Fig. 6b). Four of the modules were highly correlated with developmental stage (Pearson correlation coefficient > 0.9; P < 0.01), and the other modules were weakly correlated with five developmental stages in the treatment and control groups (Fig. 6b).
WGCNA can also be used to construct gene networks, where each node represents a gene and the connecting lines between genes represent co-expression correlations. Cytoscape v.3.6.1 was used to analyze the gene interaction network within the module and identify hub genes. The functional annotation revealed genes related to GDSL esterase/lipase, heat shock protein, caffeic acid 3-O-methyltransferase, and endoglucanase in the modules. We examined the gene regulatory networks in these modules and found that lightcyan was significantly correlated with stage C2 (r2 = 0.75), and included the VIT_07s0005g01490 (CLFSWN) gene encoding histone lysine N-methyltransferase EZA1 (Fig. 7a). Module royalblue, which was significantly correlated with C3 stage, included the VIT_19s0135g00030 (OMT1) gene encoding caffeic acid 3-O-methyltransferase (Fig. 7b). However, these two genes were not expressed after 5-azaC treatment. Module darkolivegreen1 was significantly correlated with A2 stage (r2 = 0.91) and included the hub gene VIT_06s0004g07770 encoding POD4 (Fig. 7c).
Validation by quantitative real-time (qRT-)PCR
We verified the accuracy and reproducibility of the results of the transcriptome analysis by qRT-PCR. We selected 15 genes—including nine that showed differential expression between the control and treatment groups at the same development stage, two (ATP and EFTU) from cluster 1 in the TCseq analysis, and three candidate hub genes (CLFSWN, POD4, and OMT1) identified in the WGCNA-for validation. ATP and EFTU showed similar expression levels in treated and control berries. The expression patterns of CLFSWN, POD4, and OMT1 were similar to those obtained in the transcriptome analysis. The expression levels of PAP15, XTH25, and LTA were lower than those in the control group at 35 days post-anthesis (dpa); at 45 dpa, the levels of UGT89B2, BGLU13, and GTL after treatment were higher whereas that of PE was lower than in the control group. Treatment with 5-azaC decreased HSP 23.6 expression at 55 dpa; QRT3 expression was lower at 25 and 35 dpa but higher at 45 and 55 dpa than in the control group. QRT3 expression during berry development was consistent with the transcriptome data, whereas the level of MENG in treated berries was lower than in the control group at 65 dpa. Overall, the relative expression of these genes during grape berry development as determined by qRT-PCR was consistent with and confirmed the accuracy and reliability of the transcriptome sequencing data (Fig. 8).
Tomatoes treated with 5-azaC ripen prematurely, and whole-genome bisulfite sequencing of the fruit at four stages of development has revealed many differentially methylated regions, suggesting that DNA methylation is an important regulator of fleshy fruit ripening . As a DNA methyltransferase inhibitor, 5-azaC has been used to reduce methylation levels in other plants [28, 30, 33]. We previously showed that 5-azaC decreased methylation in developing ‘Kyoho’ grape berry . 5-AzaC not only affects multiple physiological processes during plant development such as enhancing resistance and oxidation, inducing flowering, and influencing phenotype, but also affects gene expression during fruit maturation. RNA-seq is a useful tool for whole-genome expression profiling and identification of candidate genes related to this process. The TCseq results showed that the expression patterns of most genes were similar between the treatment and control groups (Fig. 4). No DEGs were identified in A1_C1 and only a few were detected in A2_C2, A3_C3, and A4_C4, confirming that 5-azaC treatment had little effect on the expression of most genes. However, there were 690 DEGs in A5_C5, which was more than the number in earlier developmental stages. Based on physiological indices and developmental state , 5-azaC was found to promote the early ripening of grape berry. Methylation normally leads to gene silencing but demethylation activates gene expression. 5-AzaC inhibited methylation but the DEGs expressed at A2 were downregulated, indicating that 5-azaC can also inhibit gene expression as previously reported . Thus, 5-azaC treatment has pleiotropic effects on gene expression and related regulatory networks.
Fruit softening involves the dissolution of the cell wall due to polymer depolymerization . The degradation of cell wall starch was shown to be the main cause of banana fruit softening/ripening [35, 36]. ‘Kyoho’ grape berry has a very hard texture at the very early stage of development, reaching peak hardness at 45 dpa before softening, which reflects the breakdown of the cell wall . Our transcriptome analysis revealed that genes associated with cell wall softening such as UGT89B2, GTL, and BGLU13 were significantly upregulated at 45 dpa after 5-azaC treatment (Table 1). Similarly, QRT3 expression was increased at 55 dpa when the berries were at the softening stage, suggesting that it regulates berry ripening. Tomato genome sequencing revealed the expression of over 50 structural genes encoding known or putative cell wall-modifying proteins during fruit development and ripening . Inhibiting PG activity was found to slow the softening of bayberry . Phospholipase (PL) activity degrades cross-linked pectin polymers, and pectin polysaccharide in the cell wall is further degraded by PG, resulting in fruit softening; silencing PL expression slows the degree of fruit softening in strawberry . As a PG, QRT3 was upregulated after 5-azaC treatment in early-ripening berries in this study. BG participates in the synthesis of abscisic acid, which is associated with grape ripening; VvBG1 overexpression in strawberry increases BG activity and promotes fruit ripening . Interestingly, we found that BGLU13 was upregulated by 5-azaC treatment.
HSPs are responsible for protein folding, assembly, translocation, and degradation during normal cell growth and development . HSP70, HSP90, and their co-chaperones are involved in signal transduction and protein targeting and degradation . Small HSP is implicated in seed, pollen, and fruit development [43, 44]. MiHSP17.6 expression, which regulates the development of mango fruit, declined during the early stage of fruit development, increased during the middle stage, and then decreased at the late stage; the highest expression level was at 60 dpa . Class II sHSP 17.4 mRNA was detected at all stages of the ripening process of tomato fruit and reached a maximal level at the later stage, whereas Class II HSP 17.6 and 17.7 mRNA had the highest expression at the turning and pink stages, respectively . FaHSP17.4 expression was increased during cell division in strawberry fruit development, but gradually decreased during maturation after cell division ceased . In this study, the HSP 23.6 gene belonging to the HSP20 family was downregulated at 55 dpa following 5-azaC treatment, with other HSPs (e.g., HSP20, HSP70, and HSP90) showing a similar trend at 65 dpa. These results indicate that 5-azaC inhibits the expression of HSP genes although the underlying mechanism remains to be determined.
Protein phosphorylation plays an important role in the regulation of growth and development of strawberry fruit. Phosphorylated proteins not only participate in transcriptional regulation and cell division but also modulate the response to plant hormones and sugar metabolism, which are related to fruit ripening and softening . In this study, GO enrichment analysis of genes at different developmental stages revealed factors that control berry development; their expression was stable and was not readily perturbed by exogenous substances. Moreover, they had related functions during berry development such as nuclease, isomerase, and endonuclease activities. However, the expression of some genes was altered at A4 by 5-azaC treatment including those encoding pyrophosphatase, nucleoside triphosphatase, and hydrolase, possibly as a result of changes in their methylation status. Meanwhile, the results of the TCseq analysis showed that most of the genes showed the same or similar expression patterns between treatment and control groups at different developmental stages, although gene expression patterns in cluster 1 were altered at 25 dpa. The GO enrichment analysis of cluster 1 genes also revealed that many were related to pyrophosphatase, nucleoside triphosphatase, and hydrolase activities.
PAP is produced by plant secretion under low phosphorus conditions and can hydrolyze phosphate groups on organophosphorus substrates and produce phosphorus for absorption and utilization by plants . The PAP-related gene MdPAP10 cloned from apple was found to be expressed in root, stem, leaf, flower and fruit, and its overexpression induced phosphorus-related gene expression . In this study, 5-azaC treatment reduced the level of PAP15 at 35 dpa, suggesting that this gene is related to grape berry development. However, the specific functions of PAP remain to be determined. In summary, our results suggest that grape berry ripening involves changes in protein phosphorylation state.
DNA methylation is generally considered to be a marker of transcriptional repression. However, decreased DNA methylation levels are also associated with the repression of genes, for instance during grape berry ripening; this is especially true in the case of genes related to photosynthesis, which are not required by ripened berries. Photosynthesis is important for the rapid growth of young fruit and starch production. With fruit ripening, the chloroplast differentiates into pigment cells and promotes coloration. Hypomethylation of DEGs is increased during chlorophyll biosynthesis. Accordingly, many photosynthesis-related genes were downregulated at A5_C5 (65 dpa). Whereas the berries of control plants were still green at this stage, those in the 5-azaC treatment group had changed color. Thus, 5-azaC may inhibit the expression of photosynthesis-related genes. However, additional studies are needed to clarify the molecular link between DNA methylation and gene expression changes during grape berry development and ripening.
Transcriptome analysis of grape berry was performed at five developmental stages to elucidate the gene expression networks controlling ripening of ‘Kyoho’ grape berry following 5-azaC treatment. Most genes were detected at similar levels between the treatment and control groups in the time series analysis, indicating that 5-azaC treatment does not significantly influence the expression of most genes. Functional annotation of DEGs revealed that they were mainly related to fruit softening, photosynthesis, protein phosphorylation, and heat stress. The results provide insight into the mechanisms that regulate grape berry development, which is useful for establishing grape varieties with specific favorable characteristics.
Plant material and 5-azaC treatment
Plant samples were collected from 6-year-old hedgerow-cultivated ‘Kyoho’ grape vines growing in Yanshi County, Luoyang, China (34.41°N, 112.46°E) in 2017, and were identified and stored in our laboratory. Field management was conducted according to local agricultural practices. A total volume of 20 ml of the methyltransferase inhibitor—100 μM 5-azaC with Silwet-77 as the surfactant (0.3% v/v) (both from Solarbio, Beijing, China)—was evenly sprayed on each berry cluster while control plants were sprayed with water as previously described . Two treatments were administered. Grape berries of uniform size and without mechanical damages were collected at five developmental stages-i.e., 25, 35, 45, 55, and 65 dpa. The berries were labeled as shown in Table 3 in order to facilitate subsequent analysis. The samples were immediately frozen in liquid nitrogen and stored at − 80 °C until further analysis.
Library preparation and transcriptome sequencing
cDNA library preparation and transcriptome sequencing was performed as previously described . Briefly, mRNA was enriched with oligo (dT) magnetic beads from total RNA. First-strand cDNA was synthesized using random hexamers with the fragmented mRNAs as templates; second strand buffer, dNTPs, RNaseH, and DNA polymerase I were then added to the reaction. The QIAQuick PCR kit and EB buffer were used for purification and elution, and end repair of the double-stranded cDNA was performed using A bases and an adapter. Fragments of the appropriate size were recovered by agarose gel electrophoresis and PCR amplification was performed to complete whole library preparation. RNA-seq libraries for the control and treatment groups were labeled as C and A, respectively. Three biological replicates of each library were sequenced on a HiSeq 2500 platform (Illumina, San Diego, CA, USA) with 150-bp paired-ends by Annoroad Gene Technology Co. (Beijing, China). RNA-seq data were uploaded to the Sequence Read Archive of the National Center for Biotechnology Information (accession number: PRJNA542248).
RNA-seq data analysis
RNA-seq data analysis was performed according to previously published protocols . To ensure data quality, clean reads were obtained by removing reads contaminated with adapters, those of low quality, and those in which N bases constituted > 5% of total bases. The Q30 of clean reads was calculated.
Bowtie2 v.2.2.3 was used to build a genome index and clean reads were aligned to the grape reference genome  using HISAT2 v.2.1.0  with the BWT algorithm . The expression levels of genes in each sample were calculated using HTSeq 0.6.0 and normalized, and are expressed as FPKM values to facilitate comparisons between different genes .
Differential expression analysis
Linear regression was used to estimate gene expression levels in each sample and the probability value (P value) of each gene that was differentially expressed between the two groups was calculated with the Wald test. The P value was corrected with the BH method by multiple hypothesis testing to obtain the q value (q). DESeq2 v.1.6.3  was used to analyze differential gene expression between the two groups. Genes with |log2[foldchange]| ≥ 1 and q < 0.05 in each pairwise comparison were identified as DEGs. The UniProt database was used to annotate and obtain detailed information on the DEGs.
TCseq analysis was performed using the R packages Vegan and Cairo to evaluate trends in gene expression at different developmental stages of grape berry using FPKM values, and genes with similar expression patterns were divided into sets ; the gene expression trends in the same set were similar under the same conditions. To investigate gene function, GO and KEGG enrichment analyses were performed using R package (clusterProfiler). GO terms and KEGG pathways with q < 0.05 were considered as significantly enriched.
To elucidate gene regulatory networks involved in grape berry development, WGCNA was performed for all genes (19,387) with FPKM ≥0.5 , with the soft threshold set to 1. WGCNA calculates the correlation between genes using a topological overlap method, which yields more biologically meaningful results. Cytoscape v.3.6.1 was used to analyze gene interaction networks in the module. Hub genes in a given module were screened with the MCODE module.
Validation of DEGs by qRT-PCR
To validate the accuracy of RNA-seq data, the expression levels of some DEGs, genes of cluster 1 identified in the TCseq analysis, and candidate hub genes were evaluated by qRT-PCR. Total RNA was isolated from ‘Kyoho’ grape berries in the treatment and the control groups using the RNAprep Pure Plant Kit (Polysaccharides & Polyphenolics-rich) (DP441; Tiangen, Beijing, China) according to the manufacturer’s instructions. HiScript II 1st Strand cDNA Synthesis Kit (R211–01; Vazyme, Nanjing, China) was used for reverse transcription. Using a suitable amount of reverse transcribed cDNA as template, qRT-PCR amplification was performed on a CFX96 Real Time PCR Detection System (Bio-Rad, Hercules, CA, USA) using TransStart Top Green qRCP SuperMix (AQ131; Trans, Beijing, China) according to the manufacturer’s instructions. The reaction conditions were 95 °C (30 s), followed by 40 cycles at 95 °C (5 s) and 56 °C (30 s). The experiment was repeated three times and fold change in gene expression level was calculated with the 2−ΔΔCt method . The VvUbiquitin1 gene was used as an internal reference. Primers used in this study are shown in Additional file 7. qRT-PCR data are presented as mean ± SD; the significance of differences between groups was evaluated with the independent samples t test.
Availability of data and materials
The RNA-Seq data supporting the results of this article have been uploaded to the Sequence Read Archive of NCBI (National Center for Biotechnology Information). It could be accessed via the NCBI SRA database with accession numbers of PRJNA542248 from 8th May 2020 onwards, as until then there is an embargo due to a complementary manuscript. Until then, these sequences are available from the corresponding author upon reasonable request.
Histone-lysine N-methyltransferase EZA1
Differentially expressed genes
Days post anthesis
Fragments per kilobase per million mapped reads
Probable galacturonosyltransferase-like 1
Heat shock protein
Kyoto Encyclopedia of Genes and Genomes
Non-specific lipid-transfer protein A
Caffeic acid 3-O-methyltransferase
Purple acid phosphatase 15
Real-Time Quantitative PCR
Weighted Gene Co-Expression Network Analysis
Probable xyloglucan endotransglucosylase/hydrolase protein 23
Pilati S, Bagagli G, Sonego P, Moretto M, Brazzale D, Castorina G, Simoni L,Tonelli C, Guella G, Engelen K, et al. Abscisic acid is a major regulator of grape berry ripening onset: new insights into ABA signaling network. Front Plant Sci. 2017;8:1093.
Robinson SP, Davies C. Molecular biology of grape berry ripening. Aust J Grape Wine R. 2000;6:175–88.
Feng C, Chen M, Xu C, Bai L, Yin X, Li X, Allan AC, Ferguson IB, Chen K. Transcriptomic analysis of Chinese bayberry (Myrica rubra) fruit development and ripening using RNA-Seq. BMC Genomics. 2012;13:19.
Hao PP, Wang GM, Cheng HY, Ke YQ, Qi KJ, Gu C, Zhang SL. Transcriptome analysis unravels an ethylene response factor involved in regulating fruit ripening in pear. Physiol Plant. 2018;163(1):124–35.
Xie M, Huang Y, Zhang Y, Wang X, Yang H, Yu O, Dai W, Fang C. Transcriptome profiling of fruit development and maturation in Chinese white pear (Pyrus bretschneideri Rehd). BMC Genomics. 2013;14:823.
Zhang A, Wang W, Tong Y, Li M, Grierson D, Ferguson I, Chen K, Yin X. Transcriptome analysis identifies a zinc finger protein regulating starch degradation in kiwifruit. Plant Physiol. 2018;178(2):850–63.
Pei M, Gu C, Zhang S. Genome-wide identification and expression analysis of genes associated with peach (Prunus persica) fruit ripening. Sci Hortic. 2019;246:317–27.
Sacco A, Raiola A, Calafiore R, Barone A, Rigano MM. New insights in the control of antioxidants accumulation in tomato by transcriptomic analyses of genotypes exhibiting contrasting levels of fruit metabolites. BMC Genomics. 2019;20(1):43.
Zhang Q, Feng C, Li W, Qu Z, Zeng M, Xi W. Transcriptional regulatory networks controlling taste and aroma quality of apricot (Prunus armeniaca L.) fruit during ripening. BMC Genomics. 2019;20(1):45.
Yu K, Xu Q, Da X, Guo F, Ding Y, Deng X. Transcriptome changes during fruit development and ripening of sweet orange (Citrus sinensis). BMC Genomics. 2012;13:10.
Wang Q, Zhao C, Zhang M, Li Y, Shen Y, Guo J. Transcriptome analysis around the onset of strawberry fruit ripening uncovers an important role of oxidative phosphorylation in ripening. Sci Rep. 2017;7:41477.
Jaillon O, Aury J, Noel B, Policriti A, Clepet C, Casagrande A, Choisne N, Aubourg S, Vitulo N, Jubin C, et al. The grapevine genome sequence suggests ancestral hexaploidization in major angiosperm phyla. Nature. 2007;449(7161):463–7.
Lijavetzky D, Carbonell-Bejerano P, Grimplet J, Bravo G, Flores P, Fenoll J, Hellín P, Oliveros JC, Martínez-Zapater JM. Berry flesh and skin ripening features in Vitis vinifera as sessed by transcriptional profiling. PLoS One. 2012;7(6):e39547.
Guo DL, Xi FF, Yu YH, Zhang XY, Zhang GH, Zhong GY. Comparative RNA-Seq profiling of berry development between table grape ‘Kyoho’ and its early-ripening mutant ‘Fengzao’. BMC Genomics. 2016;17(1):795–811.
Sun T, Xu L, Sun H, Yue Q, Zhai H, Yao Y. VvVHP1; 2 is transcriptionally activated by VvMYBA1 and promotes anthocyanin accumulation of grape berry skins via glucose signal. Front Plant Sci. 2017;8:1811.
Guo D, Li Q, Lv W, Zhang G, Yu Y. MicroRNA profiling analysis of developing berries for ‘Kyoho’ and its early-ripening mutant during berry ripening. BMC Plant Biol. 2018;18(1):285.
Kakutani T, Munakata K, Richards EJ, Hirochika H. Meiotically and mitotically stable inheritance of DNA hypomethylation induced by ddm1 mutation of Arabidopsis thaliana. Genetics.1999;151(2):831–8.
Lauria M, Rossi V. Epigenetic control of gene regulation in plants. Biochimica et Biophysica Acta Gene Regulatory Mechanisms. 2011;1809(8):369–78.
Law JA, Jacobsen SE. Establishing, maintaining and modifying DNA methylation patterns in plants and animals. Nat Rev Genet. 2010;11(3):204–20.
Zhu J. Active DNA demethylation mediated by DNA glycosylases. Annu Rev Genet. 2009;43(1):143–66.
Ma C, Jing C, Chang B, Yan J, Liang B, Liu L, Yang Y, Zhao Z. The effect of promoter methylation on MdMYB1 expression determines the level of anthocyanin accumulation in skins of two non-red apple cultivars. BMC Plant Biol. 2018;18(1):108.
Gu T, Ren S, Wang Y, Han Y, Li Y. Characterization of DNA methyltransferase and demethylase genes in Fragaria vesca. Mol Gen Genomics. 2016;291(3):1333–45.
Liu R, How-Kit A, Stammitti L, Teyssier E, Rolin D, Mortain-Bertrand A, Halle S, Liu M, Kong J, Wu C, et al. A DEMETER-like DNA demethylase governs tomato fruit ripening. Proc Natl Acad Sci U S A. 2015;112(34):10804–9.
Xu J, Xu H, Xu Q, Deng X. Characterization of DNA methylation variations during fruit development and ripening of sweet orange. Plant Mol Biol Rep. 2015;33(1):1–11.
Chen W, Kong J, Qin C, Yu S, Tan J, Chen Y, Wu C, Wang H, Shi Y, Li C, et al. Requirement of CHROMOMETHYLASE3 for somatic inheritance of the spontaneous tomato epimutation Colourless non-ripening. Sci Rep. 2015;5:9192.
Manning K, Tör M, Poole M, Hong YG, Thompson AJ, King GJ, Giovannoni JJ, Seymour GB. A naturally occurring epigenetic mutation in a gene encoding an SBP-box transcription factor inhibits tomato fruit ripening. Nat Genet. 2006;38(8):948–52.
Zhong S, Fei Z, Chen YR, Zheng Y, Huang M, Vrebalov J, McQuinn R, Gapper N, Liu B, Xiang J, et al. Single-base resolution methylomes of tomato fruit development reveal epigenome modifications associated with ripening. Nat Biotechnol. 2013;31(2):154–9.
Cheng J, Niu Q, Zhang B, Chen K, Yang R, Zhu J, Zhang Y, Lang Z. Downregulation of RdDM during strawberry fruit ripening. Genome Biol. 2018;19(1):212.
Fraga HPF, Vieira LN, Caprestano CA, Steinmacher DA, Micke GA, Spudeit DA, Pescador R, Guerra MP. 5-Azacytidine combined with 2,4-D improves somatic embryogenesis of Acca sellowiana (O. berg) Burret by means of changes in global DNA methylation levels. Plant Cell Rep. 2012;31(12):2165–76.
Huang H, Liu R, Niu Q, Tang K, Zhang B, Zhang H, Chen K, Zhu J, Lang Z. Global increase in DNA methylation during orange fruit development and ripening. Proc Natl Acad Sci U S A. 2019;116(4):1430–6.
Xu J, Wang X, Cao H, Xu H, Xu Q, Deng X. Dynamic changes in methylome and transcriptome patterns in response to methyltransferase inhibitor 5-azacytidine treatment in citrus. DNA Res. 2017;24(5):509–22.
Guo D, Li Q, Zhao H, Wang Z, Zhang G, Yu Y. The variation of berry development and DNA methylation after treatment with 5-azaC on ‘Kyoho’ grape. Sci Hortic. 2019;246:265–71.
Lang Z, Wang Y, Tang K, Tang D, Datsenka T, Cheng J, Zhang Y, Handa AK, Zhu JK. Critical roles of DNA demethylation in the activation of ripening-induced genes and inhibition of ripening-repressed genes in tomato fruit. Proc Natl Acad Sci U S A. 2017;114(22):e4511–9.
Brummell DA, Harpster MH, Civello PM, Palys JM, Bennett AB, Dunsmuir P. Modification of expansin protein abundance in tomato fruit alters softening and cell wall polymer metabolism during ripening. Plant Cell. 1999;11:2203–16.
Fan Z, Ba L, Shan W, Xiao Y, Lu W, Kuang J, Chen J. A banana R2R3-MYB transcription factor MaMYB3 is involved in fruit ripening through modulation of starch degradation by repressing starch degradation-related genes and MabHLH6. Plant J. 2018;96(6):1191–205.
Xiao Y, Kuang J, Qi X, Ye Y, Wu Z, Chen J, Lu W. A comprehensive investigation of starch degradation process and identification of a transcriptional activator MabHLH6 during banana fruit ripening. Plant Biotechnol J. 2018;16(1):151–64.
Consortium TTG. The tomato genome sequence provides insights into fleshy fruit evolution. Nature. 2012;485(7400):635–41.
Luo Z, Xu T, Xie J, Zhang L, Xi Y. Effect of hot air treatment on quality and ripening of Chinese bayberry fruit. J Sci Food Agric. 2009;89(3):443–8.
Jiménez-Bermúdez S, Redondo-Nevado J, Muñoz-Blanco J, Caballero JL, López-Aranda JM, Valpuesta V, Pliego-Alfaro F, Quesada MA, Mercado JA. Manipulation of strawberry fruit softening by antisense expression of a pectate lyase gene. Plant Physiol. 2002;128(2):751–9.
Jia H, Wang C, Zhang C, Haider MS, Zhao P, Liu Z, Shangguan L, Pervaiz T, Fang J. Functional analysis of VvBG1 during fruit development and ripening of grape. J Plant Growth Regul. 2016;35(4):987–99.
Park C, Seo Y. Heat shock proteins: a review of the molecular chaperones for plant immunity. Plant Pathol J. 2015;31(4):323–33.
Jacob P, Hirt H, Bendahmane A. The heat-shock protein/chaperone network and multiple stress resistance. Plant Biotechnol J. 2017;15(4):405–14.
Reddy PS, Kavi Kishor PB, Seiler C, Kuhlmann M, Eschen-Lippold L, Lee J, Reddy MK, Sreenivasulu N. Unraveling regulation of the small heat shock proteins by the heat shock factor HvHsfB2c in barley: its implications in drought stress response and seed development. PLoS One. 2014;9(3):e89125.
Sun W, Montagu MV, Verbruggen N. Small heat shock proteins and stress tolerance in plants. Biochim Biophys Acta. 2002;1577:1–9.
Can V, Luo C, Dong L, Liu Z, He X. Cloning and expression analysis of a small heat shock protein gene (MiHSP17.6) from Mangifera indica. Acta Horticulturae Sinica. 2014;41(12):2383–92.
Ding C, Wang CY, Gross KC, Smith DL. Reduction of chilling injury and transcript accumulation of heat shock proteins in tomato fruit by methyl jasmonate and methyl salicylate. Plant Sci. 2001;161(6):1153–9.
Li G, Zhao X, Zhou L, Houcheng Z. Cloning and expression analysis of small heat shock protein (FaHSP17.4) gene from strawberry fruit (Fragaria ananassa). Mol Plant Breeding. 2016;14(2):328–36.
Lü X, Li Y, Miao Y, Chen L, Shen Y, Qin L, Xing Y. Analysis of protein phosphorylation level at different developmental stages of strawberry fruit. Sci Agric Sin. 2016;49(10):1946–59.
Wasaki J, Yamamura T, Shinano T, Osaki M. Secreted acid phosphatase is expressed in cluster roots of lupin in response to phosphorus deficiency. Plant and Soil. 2003;248:129–36.
Li R, An J, You C, Wang X, Hao Y. Molecular cloning and functional characterization of purple acid phosphatase related gene MdPAP10 of apple. Sci Agric Sin. 2018;51(6):1182–91.
Zhong S, Joung J, Zheng Y, Chen Y, Liu B, Shao Y, Xiang JZ, Fei Z, Giovannoni JJ. High-throughput Illumina strand-specific RNA sequencing library. Cold Spring Harb Protoc. 2011;8:940–9.
Trapnell C, Williams BA, Pertea G, Mortazavi A, Kwan G, van Baren MJ, Salzberg SL, Wold BJ, Pachter L. Transcript assembly and quantification by RNA-Seq reveals unannotated transcripts and isoform switching during cell differentiation. Nat Biotechnol. 2010;28(5):511–5.
Langmead B, Salzberg SL. Fast gapped-read alignment with bowtie 2. Nat Methods. 2012;9(4):357–9.
Sirén J, Välimäki N, Mäkinen V. Indexing graphs for path queries with applications in genome research. IEEE/ACM Trans Comput Biol Bioinf. 2014;11(2):375–88.
Love MI, Huber W, Anders S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. 2014;15(12):550–70.
Wu M, Gu L. TCseq: time course sequencing data analysis; 2019.
Langfelder P, Horvath S. WGCNA: an R package for weighted correlation network analysis. BMC Bioinformatics. 2008;9:559.
Livak KJ, Schmittgen TD. Analysis of relative gene expression uata using real-time quantitative PCR and the 2−ΔΔCT method. Methods. 2001;25(4):402–8.
This work was supported by National Natural Science Foundation of China (NSFC:31672106), Zhongyuan Science Technology Innovation Leaders (194200510007) and National Key Research and Development Program of China (2018YFD1000105). The funding bodies had no role in the design of the study, collection, analysis, or interpretation of data, or in writing the manuscript.
Ethics approval and consent to participate
Consent for publication
The authors declare that they have no competing interests.
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Summary of the read numbers aligned onto the grape reference genome.
Summary of the gene expression numbers in the control and the treatment.
The cluster information of the genes from TC-seq analysis.
GO function enrichment analysis of gene expression in C1 cluster from TC-seq analysis.
The information of DEGs for A5-C5.
KEGG pathway analysis of DEGs for the treatment and the control at the same developmental stages.
The sequences of primer used in the qRT-PCR.
About this article
Cite this article
Guo, DL., Li, Q., Ji, XR. et al. Transcriptome profiling of ‘Kyoho’ grape at different stages of berry development following 5-azaC treatment. BMC Genomics 20, 825 (2019). https://doi.org/10.1186/s12864-019-6204-1