Skip to main content

Transcriptomic analysis of temporal shifts in berry development between two grapevine cultivars of the Pinot family reveals potential genes controlling ripening time

Abstract

Background

Grapevine cultivars of the Pinot family represent clonally propagated mutants with major phenotypic and physiological differences, such as different colour or shifted ripening time, as well as changes in important viticultural traits. Specifically, the cultivars ‘Pinot Noir’ (PN) and ‘Pinot Noir Precoce’ (PNP, early ripening) flower at the same time, but vary in the beginning of berry ripening (veraison) and, consequently, harvest time. In addition to genotype, seasonal climatic conditions (i.e. high temperatures) also affect ripening times. To reveal possible regulatory genes that affect the timing of veraison onset, we investigated differences in gene expression profiles between PN and PNP throughout berry development with a closely meshed time series and over two separate years.

Results

The difference in the duration of berry formation between PN and PNP was quantified to be approximately two weeks under the growth conditions applied, using plant material with a proven PN and PNP clonal relationship. Clusters of co-expressed genes and differentially expressed genes (DEGs) were detected which reflect the shift in the timing of veraison onset. Functional annotation of these DEGs fit to observed phenotypic and physiological changes during berry development. In total, we observed 3,342 DEGs in 2014 and 2,745 DEGs in 2017 between PN and PNP, with 1,923 DEGs across both years. Among these, 388 DEGs were identified as veraison-specific and 12 were considered as berry ripening time regulatory candidates. The expression profiles revealed two candidate genes for ripening time control which we designated VviRTIC1 and VviRTIC2 (VIT_210s0071g01145 and VIT_200s0366g00020, respectively). These genes likely contribute the phenotypic differences observed between PN and PNP.

Conclusions

Many of the 1,923 DEGs show highly similar expression profiles in both cultivars if the patterns are aligned according to developmental stage. In our work, putative genes differentially expressed between PNP and PN which could control ripening time as well as veraison-specific genes were identified. We point out connections of these genes to molecular events during berry development and discuss potential candidate genes which may control ripening time. Two of these candidates were observed to be differentially expressed in the early berry development phase. Several down-regulated genes during berry ripening are annotated as auxin response factors / ARFs. Conceivably, general changes in auxin signaling may cause the earlier ripening phenotype of PNP.

Peer Review reports

Background

Vitis vinifera subsp. vinifera (grapevine) belongs to the family Vitaceae. With 6,000 to 11,000 cultivars, it is one of the most important perennial crops worldwide [1]. Grapevine fruit development can be divided into two physiological phases, berry formation and berry ripening. Veraison refers to the transition from berry formation to berry ripening, and each of the two phases is represented by a sigmoidal growth curve of development [2]. The progress through development is described by stages referred to as “BBCH stages” (acronym derived from the names of the coordinating institutions involved in stage definition) that have been defined for several crops including grapevine [3, 4]. The first physiological phase is described as berry formation (berry initiation and growth with cell divisions) and lasts from the end of flowering (BBCH71) until ~60 days later when the majority of berries are touching each other (BBCH79). The developmental stage of veraison (BBCH81) is the end of berry formation and the start of berry ripening [2]. Phenotypically, veraison is the developmental switch when the berries start to soften, accompanied by the onset of accumulation of phenylpropanoids. In red grapevine cultivars, veraison is also indicated by a colour change of the berries that is caused by the beginning of accumulation of anthocyanins, a major class of phenylpropanoids. Members of the well-studied protein superfamily of R2R3-MYB transcription factors (TFs) are considered to be mainly accountable for controlling anthocyanin accumulation [5,6,7]. Berry ripening starts at veraison and continues until harvest (BBCH89), this phase includes cell enlargement, sugar accumulation and acidity decline.

Timing of veraison has also been studied at the level of genetic loci and genomic regions that control this trait. Since it is a quantitative trait influenced by several to many genetic loci, quantitative trait locus (QTL) analyses have been performed. These studies detected a major QTL for timing of the onset of veraison on chromosome 16, combined with a minor QTL on chromosome 18 [8]. By integrating a number of QTL studies, several meta-QTLs connected to genetic control of veraison time were detected, with the most relevant located on chromosome 14, 16, and 18 [9].

Anthropogenic climate change is resulting in successively earlier ripening of grapes with a significant impact on berry quality and consequently the expected flavours of a desired wine style [10]. In addition, the time of veraison and harvest of a given cultivar may differ greatly, driven by regional and/or year-specific differences in weather conditions. Obviously, this calls for a better molecular understanding of the control of ripening time in grapevine.

Comparison of different grapevine cultivars grown at the same environmental conditions often uncovers differences in the duration of berry formation, timing of veraison, duration of berry ripening, and ripening time in general. However, the underlying genetic factors are mostly unknown. Previous studies have elucidated how ripening time is affected by internal and external factors. For example, the effect of phytohormones on berry ripening has been widely studied [1]. In general, fruit growth is discussed to be controlled by several phytohormones which play essential roles to trigger or delay ripening processes [11]. Grapevine is a non-climacteric fruit and effects of abscisic acid (ABA) have been investigated in many studies as ABA is considered to trigger ripening [12,13,14]. Furthermore, it was shown that ABA is involved in controlling leaf senescence [15], responses to drought [16] and pathogen defense [17]. In grapevine, although not as central as in climacteric fruits like tomato (Solanum lycopersicum), the phytohormone ethylene is involved in the control of berry ripening [1, 13, 18, 19], while auxin has been shown to induce a delay of ripening [20, 21].

Fruit development of both, dry and fleshy fruits, has been studied very intensively for the obvious reason that fruits are central to human nutrition [22, 23]. The main model system for studies on fleshy fruits is tomato, because of established genetics and molecular biology, access to mutants, and well advanced transgenic approaches to gene function identification [24, 25]. Berry development of grapevines has also been studied intensively [1, 26] and often at the level of the transcriptome. In quite some of the studies, predominantly late berry development stages were sampled to bring the development stage of veraison into the focus [9, 27,28,29]. In addition, whole berry development was studied with coarse time point distribution [30,31,32,33,34].

To monitor gene expression changes at a high resolution throughout grapevine berry development, starting from flowering until berries are matured, we sampled a comprehensive time series from two Pinot cultivars across two years. The samples were collected from the grapevine cultivar 'Pinot Noir' (PN) and the comparably earlier ripening cultivar 'Pinot Noir Precoce' (PNP) that is expected to be closely related to PN. The cultivar PNP is listed in the Vitis International Variety Catalogue (VIVC; [35]) and described to flower at the same time as PN but to reach veraison significantly earlier than PN [36]. Quantitative data for transcript levels, interpreted as values for gene expression, were generated by RNA-Seq. We studied the general course of gene expression patterns throughout berry development in both years and cultivars, and identified a number of differentially expressed genes (DEGs) between PN and PNP prior to veraison. These DEGs can be considered as important candidates for either delaying or pushing forward berry development. Our main aim was the identification of genes controlling the speed of development, to offer an entry point into characterization of the relevant molecular functions in grapevine, and to facilitate future breeding strategies that address traits relevant to, and affected by, climate change.

Results

Phenotypical comparison between two Pinot cultivars

To study ripening shifts, we used samples of two closely related grapevine cultivars. The cultivar PNP is an earlier ripening clonal variant of its ancestor PN. Clonal relation of PN and PNP was confirmed by a set of 24 SSR markers that all displayed the identical allele status for both cultivars (Additional file 1: Table S1). To confirm and validate the phenotypic differences between PN and PNP, detailed BBCH developmental stages were determined and documented (Fig. 1). PN and PNP display similar phenotypic properties during development and flower (BBCH65) at the same time. However, veraison (BBCH81) is shifted to ~2 weeks earlier for PNP, and similar shifts were observed in four different years (Table 1). In addition, Fig. 1A shows an overview across the time points at which samples were taken. The phenotypic differences between PN and PNP are illustrated in images of developing berries taken between onset of berry formation and veraison (Fig. 1B and Additional file 1: Table S2 and Table S3). Veraison (BBCH81) is visible as the onset of anthocyanin accumulation and is detected ~2 weeks later in PN compared to PNP.

Fig. 1
figure 1

Phenotypical observations and sampling scheme. A Sampling time points and days after onset of flowering (DAF) are indicated in red. The developmental stage observed is shown in the BBCH stages [3, 4]. DAF zero (0) is set at BBCH61 (onset of flowering, 10% of flowerhoods fallen [3]). Berry development is depicted schematically and categorized into the phases flowering (yellow), berry formation (green), and berry ripening (purple) for both cultivars. The junction between green and purple indicates veraison (BBCH81). To orient for time of year, numbered days after January 1st are shown. B Images of grape bunches and developing berries taken in 2014 are shown to document the differences between PN and PNP. Images were taken 35, 41, 49 and 56 DAF. Scale bar: 50 mm

Table 1 Observed flowering- and berry development shifts between the cultivars PNP and PN in 2014, 2015, 2016 and 2017 at the Geilweilerhof vineyards, Siebeldingen, Germany (in days after January 1st)

Global view of gene expression patterns

We harvested triplicate samples in 2014 and 2017 from flowering until after veraison (for time points see Fig. 1 and Additional file 1: Tables S2 and S3), individual harvests are referred to below as subsamples, and analyzed them by RNA-Seq. After preprocessing the raw data (see Materials and methods), reads derived from each subsample were mapped to the reference sequence from PN40024 and analyzed with respect to the CRIBI V2.1 annotation dataset. For 2014, approximately 19.7 million reads per subsample were obtained from each of the 72 libraries. An overall alignment rate of 79% to the grape reference genome sequence was reached. For 2017, approximately 43.5 million reads per subsample were obtained from each of the 78 libraries. From these, an overall alignment rate of 92% to the reference was calculated. Expression values were initially detected as Transcripts Per Kilobase Million (TPM) and averaged over the three subsamples for each sample. Considering both years separately, a total of 28,692 genes were detected as expressed in both cultivars and in both years. In contrast, 2,152 CRIBI V2.1 genes were found to be not expressed.

The correlation between gene expression data, determined as TPM values per sample, of the datasets from both years over all genes was r = 0.5095 (Pearson correlation coefficient) for PN and r = 0.6557 for PNP, respectively. For PN and PNP, 10,205 and 16,226 genes, respectively, expression values were significantly correlated (p-value < 0.05) between the years 2014 and 2017. A list of the correlation strength of the eight time points with the same BBCH stage is provided in Additional file 1: Table S4.

To visualize global trends and similarity of the gene expression values obtained from all subsamples, a Principal Component Analysis (PCA) of both years was performed with keeping the subsamples separate. The first component PC1 explains 57% of the variance, while the second component PC2 explains 17% (Fig. 2). Almost all data points of the subsamples (triplicates within a sample) from both years cluster near to each other. The data follow a track of time in a nearly consecutive and continuous way. Main actors, which influence most of the variance in the data, were genes related to cell wall modification, secondary metabolism, wounding-response and hormone signaling. The top 500 genes responsible for most of the variance in PC1 and PC2 are listed in Additional file 1: Table S5 and Table S6, together with functional information for each listed gene from MapMan/Mercator and RefSeq.

Fig. 2
figure 2

Principle Component Analysis (PCA) of gene expression values from all subsamples. Each data point represents a single subsample of the triplicates for each time point of both years (2014 and 2017 as indicated by [DAF]_14 and [DAF]_17 with the colour code) and for both cultivars (PN as triangles, PNP as circles)

Cluster analysis for identification of co-expressed genes

Gene expression time series profiles combined data from the subsamples/triplicates for each time point across two years for both cultivars. Expression profiles were compared using the clustering tool CLUST. The goal was the characterization of similarity and/or differences in gene expression among years and cultivars throughout berry development. Over all four datasets, 13 PN/PNP clusters of genes with similar gene expression patterns (C1-C13) were obtained (Additional file 2: Figure S1A). In these clusters, 3,316 (12.2% ) of the 27,139 genes expressed during berry development were found co-expressed among both years and cultivars. It should be noted that CLUST uses criteria to define expressed and not expressed genes that differ from the ones applied above (see Materials and methods). The observed expression profiles differ clearly between the clusters, which was in part the result of the restricted number of clusters that CLUST extracts. Manual inspection of the clusters revealed little deviation of individual gene expression profiles within each individual cluster for a given year or cultivar. All cluster gene memberships, also those for the additional cultivar-specific cluster analyses (see below), are available in Additional file 1: Table S7-S9.

The PN/PNP clusters C2, C5, C6 and C12 (C12_PN/PNP selected as example, see Fig. 3) reveal a small but detectable difference in the gene expression profile between both sampled years, but are almost identical for both cultivars. The PN/PNP clusters C1, C7 and C11 (C1_PN/PNP selected as example, see Fig. 3) show similar expression profiles over the two years, but stand out by shifted expression peaks that distinguish PN and PNP.

Fig. 3
figure 3

Two selected gene expression profile clusters with either a cultivar-specific difference (C1_PN/PNP) or a weather/field condition-specific difference (C12_PN/PNP) after clustering all data (both cultivars and both years). Strength of gene expression (quantile normalization) was plotted over the time course of berry development. Sampling time points are detailed in Fig. 1 and were restricted to those eight equivalent time points at which the cultivars display the same BBCH stage (Additional file 1: Table S4). For all PN/PNP clusters see Additional file 2: Figure S1A. IDs of genes that make up the clusters are listed in Additional file 1: Table S7

To characterize the clusters with respect to potential functions of the co-expressed genes included in a given cluster, GO term enrichment for biological processes was calculated. The full list of enriched GO terms for all clusters is listed in Additional file 1: Table S10-S12. Two examples for GO terms appearing with highly significant incidence were ‘response to oxidative stress’ in cluster C11_PN/PNP (term GO:0051276) and ‘regulation of defense response’ in cluster C5_PN/PNP (term GO:0031347).

Two additional cluster analyses were performed, one for the PN data from both years (Additional file 2: Figure S1B) and one for the PNP data from both years (Additional file 2: Figure S1C). These analyses revealed a high abundance of genes from the expansin gene family with similar expression profiles in clusters C0_PNP and C6_PN. The cluster C16_PN was found to have a highly significant enrichment for 'vegetative to reproductive phase transition' (GO:0010228). Cluster-gene memberships for the cultivar-specific clustering are available in Additional file 1: Table S8-9, and the corresponding GO term enrichment is summarized in Additional file 1: Table S11-12.

Analyses of differentially expressed genes

The gene expression time series throughout berry development was analyzed for differentially expressed genes (DEGs) between the two cultivars PN and PNP with DESeq2. Genes with significantly differential expression were selected by using the filters adjusted p-value (PADJ) < 0.05 and log2fold change (LFC) > 2. The results are summarized in Table 2 and are detailed at the gene level per time point compared in Additional file 1: Table S13.

Table 2 Filtering steps applied for selecting DEGs, and the number of DEGs that were carried on after each selection step. For details see Materials and methods

In total, 8,206 and 4,419 DEGs were identified with PADJ greater than 0.05 for 2014 and 2017, respectively. Almost twice as many DEGs were initially detected in 2014 samples compared to 2017. By applying the filter for an at least 2-fold difference in expression level (LFC > 2), the number of significant DEGs decreased, mainly for the PN/PNP time series from 2014. Only few DEGs between PN and PNP were observed during flowering (BBCH61 to 69) at the beginning of both time series (see Fig. 4). Within the berry formation phase (BBCH71 to 79), the number of DEGs detected increased towards veraison (BBCH81) as both genotypes increasingly varied in physiological stages. The highest number of DEGs was observed in parallel to the time-shifted veraison of PNP relative to PN. This time point was also the most phenotypically different between the two cultivars (see Figs. 1 and 4). A set of veraison-specific genes was defined by selecting the DEGs from time points DAF35 and DAF41 from 2014 that were also observed to coincide with this phenotype at DAF42 and DAF49 in the 2017 gene expression data. These criteria identified 388 veraison-specific DEGs. This set of 388 genes was compared to results from similar studies and found to be in agreement (e.g. 81% [27] and 52% [28]; IDs of the 388 genes, the genes that match results from the other studies and their functional annotation, are included in Additional file 1: Table S14). During the subsequent phase of berry ripening (BBCH81 to 89, after around DAF56 in 2014 and around DAF60 in 2017), the number of DEGs detected decreased.

Fig. 4
figure 4

Visualization of the number of DEGs detected between PN and PNP in a logarithmic scale (log10). Results for 2014 are shown in purple, those for 2017 in orange. The time series from the two years were aligned at veraison of PNP; the timeline is given as days after onset of flowering (DAF). DEGs are counted for a pair of PN/PNP samples for each time point individually, the number above each column mentions the number of DEGs detected at the respective time point. Groups of newly appearing DEGs relative to an earlier time point are indicated by a new colour shade in the column (bar) for each time point. For members of a given group of DEGs, the attributed colour shade is kept for the subsequent time points (columns/samples). The pictures between the two column series display BBCH65 (full bloom, 50% of flowerhoods fallen [3]) and BBCH81 (veraison) of PNP and PN

We developed a visualization for the numbers of DEGs detected and the changes with respect to which genes are newly appearing as differentially expressed at a given time point (sample pair PN/PNP) in the time series (Fig. 4). Groups of newly appearing DEGs are containing only few genes early in berry formation, while numbers increase at veraison of PNP. After veraison of PN, the number of DEGs decreased. If DEGs appearing in several time points are counted only once, 3,342 and 2,745 unique DEGs (different genes) are detected from 2014 and 2017, respectively (compare Table 2).

To further increase the reliability, reproducibility, and relevance of the selected DEGs, the intersection between the DEGs identified in the two years studied was computed. In total, 1,923 unique DEGs were obtained (Table 2). To reveal DEGs potentially involved in the control of timing of ripening, i.e. genes that might be involved in the trait that mainly distinguishes PN and PNP, only intersecting DEGs which appeared at time points before veraison in PNP were picked. This resulted in a list of 12 DEGs that may control ripening time. It should be noted that these putative regulatory DEGs are supposed to be relevant before the set of veraison-specific genes implements the phenotypic changes at veraison. The full list of DEGs, their identity and annotation information as well as their fit to the selection criteria on the way from all (raw) DEGs to potentially regulatory DEGs are detailed in Additional file 1: Table S13. IDs of the 12 putative ripening time control genes, the genes that match results from related studies (7 DEGs [27], 4 DEGs [34] and 3 DEGs [28]) and their functional annotation are included in Additional file 1: Table S14, the most relevant data are summarized in Table 3.

Table 3 Collection of features of the set of 12 genes classified as potential regulators based on their differential expression before veraison of PNP. The detailed expression patterns are shown in Figure 5 and Additional file 2: Figure S3. Details on annotation are listed in Additional file 1: Table S13

Functional classification of DEGs

To complement the gene lists with functional information from grapevine that might potentially be informative for berry development, the 1,923 intersecting DEGs were analyzed with respect to enrichment of genes that have been assigned to biological pathways already established for grapevine (see Materials and methods). For 46 of the 247 defined grapevine pathways, significant enrichment (permuted p-value <0.1) was detected. The most reliable predictions (permuted p-value <0.001) for pathways that might be relevant were photosynthesis antenna proteins (vv10196; 9 DEGs); nitrogen metabolism (vv10910, 19 DEGs); phenylpropanoid biosynthesis (vv10940, 66 DEGs), tyrosine metabolism (vv10350, 33 DEGs); transport electron carriers (vv50105, 18 DEGs); phenylalanine metabolism (vv10360, 33 DEGs); brassinosteroid biosynthesis (vv10905, 8 DEGs) and flavonoid biosynthesis (vv10941, 30 DEGs). The enrichment results are provided in Additional file 1: Table S15. The same analysis was also carried out for the 12 putative ripening time control (Additional file 1: Table S16) and the 388 veraison-specific DEGs (Additional file 1: Table S17).

A check of the 1,923 intersecting DEGs revealed that 141 TF genes are included. Of these, 48 DEGs were clearly up- and 93 down-regulated at their first appearance in the time series. The full list of TF encoding genes that were higher expressed in PNP (up-regulated), or lower expressed in PNP (down-regulated), compared to PN, is shown in Additional file 1: Table S18.

For a more detailed view on the expression patterns of selected TF encoding genes, we generated for the TF gene family with the highest abundance among the 141 TF genes, namely the R2R3-MYB-type TFs with 22 cases in the MapMan functional assignment, an expression heatmap (Additional file 2: Figure S2). As a result, VviMYB24 (VIT_214s0066g01090), which is related to At3g27810/AtMYB21, At5g40350/AtMYB24 and At3g01530/AtMYB57 according to TAIR/PhyloGenes, was identified as an early appearing DEG that showed its highest expression level at flowering (BBCH61). Prominent R2R3-MYB genes known to be relevant for anthocyanin accumulation like VviMYBA1 (VIT_202s0033g00410), VviMYBA2 (VIT_202s0033g00390), VviMYBA3 (VIT_202s0033g00450) and VviMYBA8 (VIT_202s0033g00380) were detected as expressed starting from veraison (BBCH81) in both cultivars and with a time shift towards earlier expression in PNP. An additional R2R3-MYB gene with a similar expression pattern is VviMYB15 (VIT_205s0049g01020). Other R2R3-MYB genes are expressed early during berry formation, these include VviMYBF1 (VIT_207s0005g01210, related to At2g47460/AtMYB12/AtPFG1) as well as VviMYBPA5 (VIT_209s0002g01400) and VviMYBPA7 (VIT_204s0008g01800, both related to At5g35550/AtMYB123/AtTT2). According to their related expression patterns visualized in the heatmaps (Additional file 2: Figure S2), the R2R3-MYB genes fall into three groups that roughly fit to the three phases marked in Fig. 1B, namely flowering, berry formation, and berry ripening (see Discussion).

Putative candidates for ripening time control genes

As mentioned above, DEGs detected in both years at time points before veraison of PNP were selected and considered as putative genes that control ripening time (Table 3, Additional file 1: Table S13). The VitisNet enrichment analyses performed for these 12 candidates resulted in 2 pathways that showed significant (permuted p-value < 0.05) enrichment with two genes in the pathway: auxin signaling (vv30003 with VviEXPA5 (VIT_206s0004g00070) and VviEXPA14 (VIT_213s0067g02930)) and cell wall (vv40006 with VviPL1 (VIT_205s0051g00590) and VviGRIP28 (VIT_216s0022g00960)); see Additional file 1: Table S16).

A detailed check of the data presented in Fig. 4, together with results for these putatively ripening time control DEGs, resulted in the identification of two DEGs that stand out from the whole list of DEGs. Both genes almost completely lack expression in the early ripening cultivar PNP while there is clear expression in PN. Therefore, these two genes were detected as DEGs throughout the whole time series in both years. The first of the two, designated VviRTIC1 for "Ripening Time Control" (VIT_210s0071g01145, encoding a protein similar to "protein of unknown function DUF789"), is expressed during flowering (BBCH61 - 65) and is more or less continuously down-regulated over time in PN (Fig. 5A). The second of the two, designated VviRTIC2 (VIT_200s0366g00020, encoding a protein similar to “cysteine-rich receptor-like protein kinase”), displays expression in PN during berry formation as well as during berry ripening with a peak before veraison in 2017 (Fig. 5B). The expression patterns of the 10 remaining DEGs of the putative ripening time control gene set are shown in Additional file 2: Figure S3.

Fig. 5
figure 5

Expression patterns of VviRTIC1 (VIT_210s0071g01145) in (A) and VviRTIC2 (VIT_200s0366g00020) in (B) from RNA-Seq data of PN (blue) and PNP (red). Error bars display the standard deviation of triplicates. Left, expression profile from 2014. Right, expression profile from 2017. The y-axis represents the read counts from the output of DESeq2. The x-axis represents the development stages in days after onset of flowering (DAF)

Data on the set of 12 putative ripening time controlling genes are collected in Table 3, with a focus on relative up- or downregulation of expression before veraison of PNP. Also, a short description of the differences of the expression patterns before, at and after veraison of PN and PNP is included.

The expression patterns derived from RNA-Seq confirm each other and also allow digital quantification of low transcript accumulation levels. Nevertheless, three genes were chosen for confirmation via qRT-PCR, namely VviRTIC1, VviRTIC2 and VviERF027 (VIT 216s0100g00400). VviERF027 was included to cover a gene that displays, in different samples of the time series, differential expression as well as equally high expression in PN and PNP. The results obtained by qRT-PCR are fully congruent with the data from RNA-Seq (see Additional file 2: Figure S4). Thus, differential gene activity in a developmental pattern and in a genotype-specific way has been detected by RNA-Seq as well as qRT-PCR.

Discussion

One of the first detectable mentions of the cultivar 'Pinot Precoce' in connection with the synonym '(German) Früh Burgunder Traube, translated: Early Burgundy Grape' (PNP) is in the French book "Ampelographie retrospective" [37]. PNP is considered to be closely related to 'Pinot noir' (PN) grapes and, here, we have confirmed the clonal relationship of PN and PNP by 24 well distributed genomic SSR markers. Although this does not prove that PN is the ancestor, it is very likely that PNP was derived from the cultivar PN by somatic mutation as suggested earlier [36]. We used these two isogenic cultivars, that are distinguished by a clear duration of berry formation phenotype, to analyze changes in gene expression throughout berry development. The aim was to identify candidate genes that control the speed of berry development and veraison timing. Samples from inflorescences as well as from forming and ripening berries were collected from the onset of flowering until after PN and PNP veraison in 2014 and 2017. These samples were subjected to RNA-Seq analyses in two well resolved time series.

Phenotypic differences between the cultivars PN and PNP

The data from 4 years of careful assessment of the BBCH developmental stages of PN and PNP at the same location validate earlier observations from viticulture [36] that lead to establishment of PNP as a distinct grapevine cultivar in north European wine growing countries. Berry formation lasts about two weeks less in PNP, is clearly accelerated compared to PN and results in PNP entering veraison approximately two weeks earlier than PN (Fig. 1A, Table 1). It is reasonable to assume that this acceleration affects berry formation throughout, from immediately after fruit set until veraison. Functionally, this hypothesis implies that the genes that are responsible for the control of timing of berry development, and for the establishment of the phenotypic difference between PN and PNP, should be acting already very early in berry development, starting at least shortly after flowering and at or even before BBCH61 to BBCH79. At the end of ripening (harvest), the berries of PNP reach high sugar content earlier within the season when compared to PN.

General validation of the RNA-Seq dataset

To estimate overall data quality, the expression profiles obtained from PN and PNP were correlated for the two sampled years, 2014 and 2017. Pearson correlation was moderate, but this is expected considering the conditions of the free field environment. Exposure of the vines to external factors like biotic or abiotic stressors, including weather conditions that differ significantly between the years, also affect the transcriptome which reduces the level of correlation [38]. In a PCA, almost all data points lie on the same intended track, and biological replicates (subsamples) from both years are located close to each other. The main actors, which predominantly influence the variance in the dataset, are genes related to cell wall modification, secondary metabolism, wounding and hormone signaling. These gene categories fit logical expectations since berry development is known to (i) be controlled by hormones, (ii) require new cell walls, and (iii) be accompanied by accumulation of secondary/specialized metabolites [24, 26]. These initial results validated the quality of the dataset and indicated clearly that sampling of biologically closely related material for the subsamples/triplicates was successful.

Co-expression analysis shows similar gene expression clusters between cultivars and years

To further validate the data with respect to comparability as well as reproducibility between the two years, related gene expression profiles were identified among all genes by clustering the data from the four different time series. Generally, clusters of the same genes with similar expression patterns over time were observed for both cultivars and both years. Also, the cluster analyses for gene expression patterns among the years 2014 & 2017 in only PN and in only PNP, confirmed comparability of the gene expression patterns obtained in these two years. Detailed inspection revealed clusters representing expression profiles (and clusters of genes) with and without an environmental influence. Especially the cultivar-specific clusters C1_PN/PNP, C7_PN/PNP and C11_PN/PNP stand out. Comparison of the expression profiles for PN to those of PNP in these clusters identified a similar pattern that is moved to a different time in PNP. These findings are in agreement with the shifted ripening time phenology of the two cultivars discussed above. In contrast, the clusters C2_PN/PNP, C5_PN/PNP, C6_PN/PNP and C12_PN/PNP display more pattern similarity among the two years than among the two cultivars. Thus, the genes in these four clusters may display dependence on environmental factors in their expression patterns, potentially due to differences in the weather conditions between the two years studied. Strong environmental effects on gene expression patterns have also been described for grapevine berry development at 11 different environments (vineyards) from northern Italy [39]. The remaining other PN/PNP clusters C0, C3, C4, C8, C9 and C10 display highly similar expression profiles over all four conditions. The genes included in these clusters are probably less affected by environmental factors and/or the genotypic difference between PN and PNP. We conclude that our RNA-Seq results and expression level comparisons between two years are based on valid data.

When the genomic location of the DEGs is analyzed, a genome region on chromosome 16 comes into focus. In this region, 54 of the DEGs from the set of 1,923 intersecting DEGs (Table 2) are located. Of these, 28 encode stilbene synthases [40] that are all up-regulated after veraison of PNP (BBCH83). Stilbenes are a group of phenylpropanoid compounds (that includes resveratrol) which are detected in many plants, which often accumulate in response to biotic and abiotic stresses, and which are formed as a basic structuce by the key enzyme stilbene synthase. The genome region fits to a major QTL (Ver1) for "timing of the onset of veraison" on chromosome 16 [8]. It remains to be determined if this aggregation of DEGs is by chance. Potentially, the observation is biased by co-regulation of a large number of closely linked stilbene synthase genes.

Differentially expressed genes throughout berry development and identification of veraison-specific genes

Differential gene expression analysis and subsequent filtering revealed 1,923 DEGs between PNP and PN. DEG detection was based on a comparison of samples taken from the two cultivars at very similar DAF. As expected for the characterized phenotype, PNP reaches veraison when PN is still in the phase of berry formation. Consequently, the strong increase in the number of detected DEGs shortly before and at veraison of PNP results from the different developmental stage of PNP compared to the lagging PN. Subsequently, when also PN enters veraison, the number of DEGs declines (note that Fig. 4 uses a logarithmic scale). A list of 388 genes that show up in both years with a veraison-specific expression pattern was extracted and compared to published results. Interestingly, about 81.5% of the 388 PN/PNP veraison-specific genes were also described in the 4,351 differentially expressed genes between the table grape cultivar ‘8612-CTR’ (wild type) and its early ripening bud mutation ‘8612-ERN’ [27]. Also, analyses of berries from the cultivars 'Cabernet Sauvignon' and 'Pinot Noir' by RNA-Seq identified a gene set of 5,404 genes marking the onset of berry ripening [28]. This set covers 51.5% of the 388 PN/PNP veraison-specific genes (Additional file 1: Table S14). Several "switch genes" which are supposed to encode key regulators of the developmental transition at veraison [34, 41] are also included in the 388 veraison-specific gene set (Additional file 1: Table S14). We conclude that the PN/PNP veraison-specific set of 388 genes represents a core set of genes that are relevant for executing the switch from berry formation to berry ripening. A relatively small gene set was detected that still displays high overlap to those found by other studies that addressed similar biological questions, indicating that the specific experimental setup and implemented filters used here are appropriate to remove unrelated genes. Here, the comparison of "wildtype to mutant" RNA-Seq results in isogenic background between PN and PNP, reduced environmentally controlled transcriptome differences by sampling in the same vineyard/location, and dense time course sampling together with high RNA-Seq read coverage allowed good resolution power.

In order to check for potentially co-expressed genes within the veraison-specific gene set, the memberships of these genes in the PN/PNP cluster analysis were investigated. A total of 48 veraison-specific DEGs were detected in cluster C6_PN/PNP (contains 914 genes). These 48 genes include several prominent ripening-related genes like VviGRIP61 (VIT_201s0011g05110), VviMYBA8, VviMRIP1 (VIT_205s0049g00760, [42]), VviGRIP4 (VIT_205s0049g00520) and VviGRIP28. The ~20 VviGRIP genes were previously detected by differential cDNA screening as ripening-induced genes in grape [43]. Another relevant cluster is C5_PN/PNP (contains 263 genes) which includes 37 of the 388 veraison-specific DEGs. Among these are VviMYBA1, VviMYBA2, VviMYB15 and VviGRIP22 (VIT_206s0004g02560). The two clusters C5_PN/PNP and C6_PN/PNP show quite similar patterns (Additional file 2: Figure S1A). It was, at first, not obvious which difference has forced CLUST to put a given gene in either C5_PN/PNP or C6_PN/PNP. However, a comparison of the expression patterns of VviMYBA2 (in C5_PN/PNP) and VviMYBA8 (in C6_PN/PNP; see Additional file 2: Figure S2 for a heatmap) shows that there are borderline cases regarding assignment to either C5_PN/PNP or C6_PN/PNP.

In total, 22 genes encoding R2R3-MYB TFs were found among the 1,923 intersecting DEGs. Based on the timing of expression in PN and PNP, the 22 R2R3-MYB genes can be classified into three groups (Additional file 2: Figure S2). The first group is represented by VviMYB24 which is expressed during early flowering (BBCH61) but switched off already at the transition from flowering to berry formation (BBCH71). VviMYB24 is potentially orthologous to a group of three A. thaliana R2R3-MYB genes (AtMYB21/24/57) that are expressed in flowers and which function redundantly to regulate stamen development in the context of jasmonate action [44]. It is tempting to speculate that VviMYB24 has a similar function in grape.

The second group covers about 15 R2R3-MYB genes that are expressed during berry formation and pre-veraison (BBCH71 to 77). This group includes VviMYBF1 which regulates flavonol biosynthesis [45], and VviMYBPA5 as well as VviMYBPA7 which belong to the clade of AtTT2-related genes that control proanthocyanidin (PA, flavan-3-ol) biosynthesis [7, 46,47,48]. The other R2R3-MYB genes in this group are less well characterized although there are functions described for some of them, e.g. VviMYBC2-L3 (VIT_214s0006g01620) as repressor of specific branches of the phenylpropanoid pathway [49].

The third group of R2R3-MYB genes is active starting at veraison (after BBCH81) and covers about six genes. Among them are the anthocyanin accumulation controlling genes, VviMYBA1, VviMYBA2, VviMYBA3 and VviMYBA8 for which there is good evidence that they trigger anthocyanin biosynthesis [50]. Since PN and PNP are red berry cultivars, activity of the TF genes that direct anthocyanin accumulation is expected. In addition, this group includes VviMYB14 and VviMYB15 that are supposed to regulate the stilbene biosynthetic pathway [51]. With regard to the heatmaps (Additional file 2: Figure S2) and the analyses of the DEGs in this study in general, it should be noted that while the resolution within the developmental program and time is quite good, our data do not resolve the exact location of gene expression. Therefore, it remains to be determined if the expression detected is derived from berry skin, flesh, the seed or other tissues/cells.

Putative ripening time control genes acting early in berry development

To focus on genes that are contributing to the acceleration of berry formation in PNP, and/or to the control of timing of veraison, we selected DEGs detected at time points prior to veraison of PNP (Table 2, Fig. 4). This resulted in a set of only 12 genes that are potentially involved in the regulation of ripening time (Table 3). According to our hypothesis that the genes relevant for acceleration of berry formation in PNP, which cause the earlier onset of ripening in PNP, should be acting from at least shortly after flowering, we designated this set of genes as "ripening time controlling". However, genes that encode components of the respective regulatory networks and target genes of regulators including secondarily affected DEGs are surely included as well [52]. The 12 putative ripening time control genes, the DEGs detected before veraison of PNP, encode proteins related to auxin action, pectin processing enzymes related to cell wall modification, TFs from the HD-Zip as well as NF-Y/LEC families, a cysteine-rich receptor-like protein kinase, an oleosin, and proteins with domains of unknown function.

The expression patterns of this set of differentially expressed genes, and also the complete dataset of DEGs detected (Table 2), was screened for genes that were higher expressed in PNP than in PN. However, although there are many genes that display earlier upregulation of expression in PNP than in PN in the context of earlier veraison, none of the 12 "early differential" putative ripening time control genes was significantly higher expressed in PNP than in PN already at or early after flowering (Table 3). Also, among the other DEGs, no gene with such a differential expression pattern was detected. Such early gene activation in PNP compared to PN could hint at a dominant regulator that promotes faster ripening, but the data are more in favor of loss of an inhibitor of fast ripening.

The two genes assigned to auxin signaling by VitisNet (vv30003) encode expansins (VviEXPA5 and VviEXPA14, [53]). Expansins are known to be involved in fruit ripening through cell wall expansion and cell enlargement [54]. Auxin can delay the onset of veraison and ripening processes in grapevine [19,20,21]. Since reduced expression of genes from the auxin signaling pathway may indicate reduced auxin action due to lower auxin levels, the accelerated entry of PNP into veraison might be initiated by reduced auxin levels. Additionally, the genes VviPL1 (pectate lyase 1 [55]), VviPME10 (pectin methylesterase 10, VIT_206s0009g02560) and VviGRIP28 (encoding a pectin methylesterase inhibitor precursor-like protein) are also related to cell wall processes, indicating that cell wall modification is an important target process also prior to veraison [54]. The gene VviGRIP28 was also detected within a veraison-specific meta-QTL designated ver/ph16.1 [9]. It remains to be determined if this correlation has a functional basis.

The two genes in the set of 12 that encode TFs are VviHDZ28 (VIT_216s0098g01170, [56]) and VviLEC1 (VIT_200s0956g00020, [57]). The V. vinifera gene VIT_216s0098g01170 that has been designated VviHDZ28 has also been considered as a homolog of AtHB12 (At3g61890), but it seems that VviHDZ07 (VIT_202s0025g02590) and VviHDZ27 (VIT_215s0048g02870) are more similar to AtHB12. In these cases, which lack clearly assignable homologs, transfer of functional information reaches its limits and might be restricted to concluding that VviHDZ28 is important for organ development in Vitis. The gene VviLEC1 is one of three genes in V. vinifera which are homologs of AtLEC1 (At1g21970, NF-YB9) and AtL1L (LEC1-like, At1g21970, NF-YB6). LEC1 and L1L are central regulators of embryo and endosperm development. They control, among other processes, embryo morphogenesis and accumulation of storage reserve [58]. It is tempting to speculate that the reason for the detection of VviLEC1 among the 12 putative ripening time control genes is that also seed development needs to be accelerated in PNP compared to PN. This would explain earlier and higher expression of VviLEC1 in PNP compared to PN as observed (Additional file 2: Figure S3G). Consequently, VviOLE5 (VIT_216s0013g00880, encoding an oleosin involved in oil body formation [59]) would fit into the picture as relevant for lipid storage during seed development. According to the proposed enzyme function as an alcohol acyltransferase by the protein encoded by VIT_209s0018g01490 involved in volatile ester formation [60], this gene could play a similar role. For the gene VIT_205s0077g01980 no functional annotation is available (uncharacterized protein), although homologs exist throughout the Magnoliophyta.

Candidates for causal genes explaining the difference between PN and PNP

Among the 12 putative ripening time control genes, of which 10 are discussed above, two are especially interesting. Detailed analyses of the full set of DEGs, visualized in Fig. 4, resulted in the identification of VviRTIC1 and VviRTIC2, that could possibly be centrally involved in the accelerated berry development and earlier beginning of ripening in PNP compared to PN. The special feature of the expression patterns of the two genes (Fig. 5, Table 3) is that both are differentially expressed already at the first time point analyzed which was selected to hit the BBCH stage 61 (flowering before full bloom, DAF zero (0)). Also, both genes are only barely expressed in PNP in both years studied, while expression in PN is high at almost all time points. VviRTIC1 is annotated to encode a protein containing a domain of unknown function (DUF789), while VviRTIC2 is annotated to encode a "cysteine-rich receptor-like protein kinase" (acronym CRK). The best BLASTp hit to A. thaliana protein sequences indicates that it is related to At4g23180/AtCRK10, but a closer inspection shows that similarity to At4g05200/AtCRK25, At4g23160/AtCRK8 and At4g23140/AtCRK6 is almost as high. This ambiguity, and also the fact that the V. vinifera genome contains several genes related to VviRTIC2 (e.g. VIT_210s0071g01200, VIT_202s0087g01020 or VIT_203s0017g01550 as listed by PhyloGenes), complicates transfer of functional information. CRKs are a subgroup of plant receptor-like kinases [61] and are encoded by a family of 44 genes in A. thaliana. In a systematic analysis of the functions of A. thaliana CRKs, evidence was collected for involvement in the control of plant development, biotic and abiotic stress responses, photosynthesis as well as stomatal regulation [62]. This systematic phenotypic screen of a large set of T-DNA insertion mutants revealed distinct phenotypes for various of the CRK genes, but assignment of a molecular function to individual CRKs beyond recognition of unknown ligands and signal transmission by phosphorylation remains a large challenge.

As pointed out above, it is very possible that the genes we have identified are part of a genetic pathway that controls timing of berry development in V. vinifera, but that we have hit genes in a downstream part of this pathway. The relevance of the two candidate genes in the causal genetic difference between PN and PNP remains to be determined. Phase-separated genome sequences of the cultivars will be required to resolve the genome structure of both alleles of VviRTIC1 and 2 the genes in PN and PNP for an informative comparison. In future studies, we will address this question, for example by long read DNA sequencing.

Conclusions

This study detected 1,923 DEGs between the Pinot cultivars PN and PNP. The two clonal cultivars display a phenotypic difference in berry development timing where PNP reaches from full bloom to veraison faster than PN. We defined 388 DEGs as veraison-specific and 12 DEGs as putatively controlling ripening time. The relatively small number of veraison-specific genes displays a very high overlap with results published for similar studies (see Additional file 1: Table S14) and could be used for studying a phytohormone network that is acting similarly in PN and PNP, but accelerated in PNP. Additionally, the ripening time control genes identified here might offer access to a set of genes putatively important for triggering or delaying the start of berry ripening in grapevine. Further investigations are needed to elucidate structural differences in the genomes, the function of the observed DEGs, and their role in shifting the onset of ripening in grapevine.

Material and methods

Plant material and analysis of clonal relation

The grapevine (Vitis vinifera subsp. vinifera L.) cultivar PNP (Pinot Precoce Noir, VIVC No. 9280) is early ripening and has been described to be related to the cultivar PN (Pinot Noir, VIVC No. 9279) [35] that ripens later than PNP. To prove the clonal relation, DNA from both cultivars was genotyped utilizing 24 polymorphic SSR markers (VVS2, VVMD7, VVMD5, VVMD32, VVMD28, VVMD27, VVMD25, VVMD24, VVMD21, VVIV67, VVIV37, VVIQ52, VVIP60, VVIP31, VVIN73, VVIN16, VVIH54, VVIB01, VrZAG83, VrZAG79, VRZAG67, VrZAG62, VMC4F3.1, VMC1B11) as described [63]. The two cultivars used have been identified as accession DEU098_VIVC9280_Pinot_Precoce_Noir_DEU098-2008-076 and DEU098_VIVC9279_Pinot_Noir_DEU098-2008-075, respectively. The tissue used for harvest is indicated below and in Fig. 1. Both cultivars do not belong to an endangered species and were obtained and are grown in accordance with German legislation.

Phenotypical characterization and sampling of plant material

Plant material was harvested from PN and PNP grapevines trained in trellis. The plants are growing at the vineyards of JKI Geilweilerhof located at Siebeldingen, Germany (N 49°21.747, E 8°04.678). The grapevine plants were planted with an interrow distance of 2.0 m and spacing of 1.0 m in north-south direction. Inflorescences, developing and ripening berry samples of PNP and PN for RNA extraction were collected in two years with three independent biological replicates (subsamples) each. Sampling took place at systematic time points (12 time points in 2014, 13 time points in 2017), and at approx. 8 a.m. each day. In 2014, harvesting took place regularly every 7 days with only two exceptions (one day deviation, DAF13 and DAF27). In 2017, harvesting was adapted to BBCH stages (Fig. 1A). The timeline in both years is described as days after onset of flowering (DAF), with onset of flowering defined as the day at which 10% of the individual flowers have lost their caps (BBCH61 [3]). For each subsample within the triplicates, material from two neighboring grapevines was selected. Grapevine plants were weekly phenotyped according to BBCH stage [3, 4]. Phenotyping was performed repeatedly to ensure sampling from vines of the same development stage (e.g. percentage of open flowers during flowering, or berry development stage) to reach uniform subsamples. The phenotypical observations were summarized in Additional file 1: Table S2 and S3. From these, the durations of flowering, berry formation and berry ripening as well as the resulting shifts between the cultivars were calculated (Table 1). Furthermore, images from berry developmental stages of both cultivars were taken in 2014 for 35, 41, 49 and 56 DAF. The sampled material was directly frozen in liquid nitrogen and stored at -70°C until RNA extraction.

RNA extraction and cDNA library construction

Biological replicates, i.e. the subsamples, were ground separately under liquid nitrogen. Total RNA was extracted using an RNA Isolation Kit (Sigma-Aldrich Spectrum™ Plant Total RNA) according to suppliers' instructions. Quality control, determination of RIN numbers [64] and estimation of the concentrations of all RNA samples was done on a Bioanalyzer 2100 (Agilent) using RNA Nano 6000 Chips. For RNA-Seq, 500 ng total RNA per subsample were used to prepare sequencing libraries according to the Illumina TruSeq RNA Sample Preparation v2 Guide. For subsamples from 2014 and 2017, 72 and 78 libraries were constructed and sequenced, respectively. Enrichment of poly-A containing mRNA was performed twice, using poly-T oligos attached to magnetic beads included in the Illumina kit. During the second elution of the poly-A+ RNA, the RNA was fragmented and primed for cDNA synthesis. After cDNA synthesis, the fragments were end-repaired and A-tailing was performed. Multiple indexing adapters were ligated to the ends of the cDNA fragments and the adapter ligated fragments were enriched by 10 cycles of PCR. After quality check using Bioanalyzer 2100 HS-Chips (Agilent) and exact quantification by Quant-iT PicoGreen dsDNA assay on a FLUOstar Optima Plate Reader (BMG LABTECH), the libraries were pooled equimolarly.

RNA-Seq

Single end (SE) sequencing of the pooled barcoded libraries from 2014 was performed on an Illumina HiSeq1500 in High Output mode generating 100 nt reads. For samples from 2017, sequencing was done using an Illumina NextSeq500 generating 83 nt SE reads; two runs were performed with the same pool of barcoded libraries from 2017.

Processing of RNA-Seq read data

Raw reads were trimmed with Trimmomatic (version 0.36) [65]. For raw reads from the year 2014, the following settings were used: LEADING:10 TRAILING:10 SLIDINGWINDOW:4:15 MINLEN:50. In addition, a collection of all available Illumina adapter sequences was supplied to remove matches within the parameter 2:30:10. For raw reads from the year 2017, trimming settings were set to LEADING:6 TRAILING:6 SLIDING WINDOW:4:15 MINLEN:36. All trimmed reads were quality-checked via FastQC (version 0.11.8) [66]. Thus, possible adapter sequences and low-quality bases were removed. All trimmed reads passing QC were mapped to the reference genome sequence PN40024 (version 12Xv2) [67] using the graph-based alignment tool HISAT2 (version 2.1.0) [68, 69] with no additional soft clipping. Afterwards, all tagged genes (structural gene annotation: CRIBI v2.1) were counted as raw read counts with FeatureCounts (Bioconductor package Rsubread version 3.8 [70]). To estimate transcript abundance as a measure for gene expression, counts for Transcripts Per Kilobase Million (TPM, [71]) were determined.

Basic gene expression analyses

TPM counts from the various samples were used for manual gene expression inspection, for determination of the number of expressed and not expressed genes, and to calculate the correlation between gene expression values from both years. Genes with a TPM value > 0 added up over all samples from one year were classified as expressed, conversely genes with a TPM value = 0 added up over all samples as not expressed. A custom python script was applied utilizing the function pearsonr from SciPy python package (v. 1.2.3) [72], which calculates the Pearson correlation coefficient and the p-value for all year-to-year comparisons. Expression data pairs for TPM counts per gene from both sampled years, averaged over the three subsamples of each sample, were used. To test correlations and relationships between expression values from the two years, where samples were harvested with slightly different sampling patterns (see Fig. 1A), eight equivalent time points with the same BBCH stages between the cultivars of each year were selected (see Additional file 1: Table S4).

Principal component analysis

To explore data similarity, a Principal Component Analysis (PCA) was calculated over all gene expression values from both years and cultivars for all subsamples. All data points were normalized using variance stabilizing transformation function 'vst' from the R package DESeq2 (v. 1.12.4) [73]. Subsequently, the principal components were generated using 'prcomp' from the R package 'stats' (v. 3.5.2) [74]. The resulting PCA object, displaying the main components PC1 and PC2, was plotted and exported. Additionally, genes with the highest variance contribution to PC1 and PC2 were extracted separately.

Functional annotation of genes

Transfer of annotation information from other plant species, mainly A. thaliana, was calculated using MapMan's sequence annotation tool Mercator (v. 3.6) [75, 76]. Additionally, all open reading frame (CDS from V. vinifera/grapevine genes) sequences were aligned to the non-redundant protein sequence data base RefSeq [77] with the basic local alignment tool for proteins BLASTp [78] (e-value ≤ 0.001). Short descriptions of gene functions were extracted and added to the gene lists in Additional file 1: Table S5, S6, S13, S14.

GO term enrichment for biological processes was calculated via the R package 'topGO' (v. 2.38.1) [79]. Subsequently, statistical reliability was calculated using Fishers exact test. All Gene IDs and their corresponding GO terms were extracted from the CRIBI database (http://genomes.cribi.unipd.it/DATA/V2/annotation/bl2go.annot_with_GO_description.txt). All results of the GO term enrichment are deposited in Additional file 1: Table S10-12.

Cluster analysis

To reveal co-expressed genes over all four datasets, the tool CLUST (v. 1.10.8) was used with default parameters [80]. As input, raw read counts from eight time points were used. These time points were selected to cover the same BBCH stages of PN and PNP from the years 2014 and 2017 (Additional file 1: Table S4). First, all data were pre-processed as described in the CLUST manual. Values from corresponding subsamples (triplicates) were combined and averaged. To filter out uninformative (very low) gene expression values, an additional filter was applied: genes not reaching a sample expression value > 1 in at least three conditions and in at least one cultivar from one year were discarded (-fil-v 1 -fil-c 3 -fil-d 1). Afterwards, the data were quantile normalized according to the RNA-Seq defaults of CLUST. Genes showing a flat expression profile were filtered out by applying the default settings [80].

Differential gene expression analyses

For analyses of differentially expressed genes, DESeq2 (v. 1.12.4; R Bioconductor) was employed. To test if gene expression differs significantly between two samples, the likelihood ratio test nbinomLRT, included in the DESeq2 package, was used. As input, raw read counts from all time points were used. Normalization factors and dispersion estimates were used as described [73]. The output table contained all differentially expressed genes (DEGs) and the corresponding values for baseMean, log2FoldChange (LFC), lfcSE (LFC standard error), stat (difference in deviation between the reduced model and the full model), p-value and PADJ (adjusted p-value). To focus on significantly differentially expressed genes from the DESeq2 analyses, cut-off filters PADJ ≤ 0.05 and LFC > 2 were applied.

Confirmation of differential gene expression by qRT-PCR

To verify the RNA-Seq results, four time points from 2017 were selected (DAF zero (0), DAF28, DAF57 and DAF77) for qRT-PCR. Synthesis of cDNA from the RNA subsamples was carried out with First Strand cDNA Synthesis Kit (ProtoScript® II; NEB) according to the manufacturer's instructions. The qRT-PCR assay was performed using Luna Universal qPCR Master Mix (NEB) with a total volume of 20 µl. Sequences of the primers used are listed in Additional file 1: Table S23. Reaction products/amplicons were detected based on SYBR green via a CFX96 Real Time PCR Detection System (Bio-Rad). For each time point, three biological and three technical replicates (i.e., each subsample in triplicate) were measured. Cycling conditions after initial denaturation 2 min at 95°C: denaturation 5 sec at 95°C, annealing/extension 30 sec at 60°C, cycled 35 times. For QC, each reaction was controlled by product melt analysis (65°C - 95°C). As negative controls, no template control (NTC) and no reverse transcriptase control (-RT) were measured as well using three technical replicates. The polyubiquitin gene VviUbiquitin1 (VIT_216s0098g01190) was used for normalization [46]. Measurements were analyzed using the CFX Maestro V.4.1.2433.1219 (Bio-Rad) by normalization via the relative quantitative ΔΔCt method.

Selection of gene sets potentially relevant for ripening time and comparison with literature data

In order to identify gene sets from the DEGs relevant for control and implementation of ripening, an intersection between the DEGs detected at all time points between both years was built. To determine a subset of putatively ripening time control genes, the intersection between both years covering the development stages BBCH61 (onset of flowering) to BBCH79 (one developmental BBCH stage before veraison) was used (time points 2014: DAF 0-35; 2017: DAF 0-42). Furthermore, a set of veraison-specific genes was defined from the DEGs detected at the intersection of development stages BBCH79 (one developmental BBCH stage before veraison) to BBCH81 (onset of ripening / veraison; time points 2014: DAF35-41; 2017: DAF42-49). To test for biological relevance of the subsets, all DEGs were screened to their occurrence in similar relevant studies [9, 27, 28, 30, 34, 41, 81, 82].

Visualization of gene numbers newly appearing as differentially expressed

To visualize appearance of DEGs over time, a stacked bar plot script was set up using the R package ‘plotly’ (v. 4.9.2.1) [83]. Each bar represents the amount of DEGs of a given time point or condition. In order to track groups of DEGs newly appearing at a given time point throughout the following time points, the colour shade representing the group of DEGs remains the same.

Pathway enrichment analysis

To search for possible targets in known pathways of grapevine, a pathway enrichment analysis using the tool VitisPathways [84] was performed. To achieve a reliable enrichment, 1000 permutations, a Fisher’s exact test of p < 0.05 and permuted p-value < 0.1 were set. Thus, all significant enriched pathway genes and their relations can be displayed in VitisNet [85], a specific molecular network for grapevine.

Heatmaps

As an extension to assignment of genes to biosynthesis pathways, the genes were also filtered for annotation as coding for transcription factors (TFs). This filter was based on the annotation information transferred from Mercator and RefSeq (see above). To look at the entire family of R2R3-MYB TF genes, the list of MYB genes identified via MapMan was extended by additional grapevine R2R3-MYB gene family members that have been characterized [6, 7]. The R2R3-MYB genes detected among the intersecting DEGs were displayed in heatmaps addressing the four individual time series (2 cultivars, 2 years) using the R package 'pheatmap' (v. 2.1.3) [86]. Predictions for phylogenetic relationships were deduced from PhyloGenes v. 2.2 [87].

Availability of data and materials

The FASTQ files containing all RNA-Seq reads (PN2014, PNP2014, PN2017 and PNP2017) have been deposited at the European Nucleotide Archive (ENA) according to the INTEGRAPE guidelines under the accession numbers PRJEB39262, PRJEB39261, PRJEB39264 and PRJEB39263, respectively (see Additional file 1: Tables S19 to S22). All scripts developed for this study are available on GitHub [https://github.com/bpucker; https://github.com/jenthein].

Abbreviations

ABA:

Abscisic acid

CRK:

Cysteine-rich Receptor-like protein Kinase

DAF:

Days after onset of flowering

DEGs:

Differentially expressed genes

PNP:

Pinot Noir Precoce

PN:

Pinot Noir

LFC:

Log2fold change

PADJ:

Adjusted p-value

PCA:

Principle component analysis

QTL:

Quantitative Trait Locus

TF:

Transcription factor

TPM:

Transcripts per kilobase million

References

  1. Kuhn N, Guan L, Dai ZW, Wu BH, Lauvergeat V, Gomès E, Li SH, Godoy F, Arce-Johnson P, Delrot S. Berry ripening: recently heard through the grapevine. J Exp Bot. 2014;65(16):4543–59.

    Article  CAS  PubMed  Google Scholar 

  2. Coombe BG. Research on development and ripening of the grape berry. Am J Enol Viticult. 1992;43(1):101–10.

    Article  Google Scholar 

  3. Lorenz DH, Eichhorn KW, Bleiholder H, Klose R, Meier U, Weber E. Phenological growth stages of the grapevine (Vitis vinifera L. ssp. vinifera) - codes and descrptions according to the extended BBCH scale. Aust J Grape Wine Res. 1995;1(2):100–3.

    Article  Google Scholar 

  4. Bleiholder H, Weber E, Feller C, Hess M, Wicke H, Meier U, van den Boom T, Lancashire PD, Buhr L, Hack H et al. Growth stages of mono- and dicotyledonous plants - BBCH Monograph: Federal Biological Research Centre for Agriculture and Forestry; 2001.

  5. Dubos C, Stracke R, Grotewold E, Weisshaar B, Martin C, Lepiniec L. MYB transcription factors in Arabidopsis. Trends Plant Sci. 2010;15(10):573–81.

    Article  CAS  PubMed  Google Scholar 

  6. Matus JT, Aquea F, Arce-Johnson P. Analysis of the grape MYB R2R3 subfamily reveals expanded wine quality-related clades and conserved gene structure organization across Vitis and Arabidopsis genomes. BMC Plant Biol. 2008;8:83.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  7. Wong DCJ, Schlechter R, Vannozzi A, Holl J, Hmmam I, Bogs J, Tornielli GB, Castellarin SD, Matus JT. A systems-oriented analysis of the grapevine R2R3-MYB transcription factor family uncovers new insights into the regulation of stilbene accumulation. DNA Res. 2016;23(5):451–66.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  8. Zyprian E, Ochßner I, Schwander F, Šimon S, Hausmann L, Bonow-Rex M, Moreno-Sanz P, Grando MS, Wiedemann-Merdinoglu S, Merdinoglu D, et al. Quantitative trait loci affecting pathogen resistance and ripening of grapevines. Mol Genet Genomics. 2016;291(4):1573–94.

    Article  CAS  PubMed  Google Scholar 

  9. Delfino P, Zenoni S, Imanifard Z, Tornielli GB, Bellin D. Selection of candidate genes controlling veraison time in grapevine through integration of meta-QTL and transcriptomic data. BMC Genomics. 2019;20(1):739.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  10. Delrot S, Grimplet J, Carbonell-Bejerano P, Schwandner A, Bert P-F, Bavaresco L, Costa LD, Di Gaspero G, Duchêne E, Hausmann L, et al. Genetic and genomic approaches for adaptation of grapevine to climate change. In: Kole C, editor., et al., Genomic designing of climate-smart fruit crops. Cham: Springer; 2020. p. 157–270.

    Chapter  Google Scholar 

  11. Fenn MA, Giovannoni JJ. Phytohormones in fruit development and maturation. Plant J. 2020.

  12. Osorio S, Scossa F, Fernie AR. Molecular regulation of fruit ripening. Front Plant Sci. 2013;4:198.

    Article  PubMed  PubMed Central  Google Scholar 

  13. Sun L, Zhang M, Ren J, Qi J, Zhang G, Leng P. Reciprocity between abscisic acid and ethylene at the onset of berry ripening and after harvest. BMC Plant Biol. 2010;10:257.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  14. 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.

    Article  PubMed  PubMed Central  Google Scholar 

  15. Zhu Z, Li G, Yan C, Liu L, Zhang Q, Han Z, Li B. DRL1, encoding A NAC transcription factor, is involved in leaf senescence in grapevine. Int J Mol Sci. 2019;20(11):2678.

    Article  CAS  PubMed Central  Google Scholar 

  16. Rossdeutsch L, Edwards E, Cookson SJ, Barrieu F, Gambetta GA, Delrot S, Ollat N. ABA-mediated responses to water deficit separate grapevine genotypes by their genetic background. BMC Plant Biol. 2016;16(91):s12870-12016-10778–2874.

    Google Scholar 

  17. Coelho J, Almeida-Trapp M, Pimentel D, Soares F, Reis P, Rego C, Mithöfer A, Fortes AM. The study of hormonal metabolism of Trincadeira and Syrah cultivars indicates new roles of salicylic acid, jasmonates, ABA and IAA during grape ripening and upon infection with Botrytis cinerea. Plant Sci. 2019;283:266–77.

    Article  CAS  PubMed  Google Scholar 

  18. Chervin C, El-Kereamy A, Roustan J-P, Latché A, Lamon J, Bouzayen M. Ethylene seems required for the berry development and ripening in grape, a non-climacteric friut. Plant Sci. 2004;167(6):1301–5.

    Article  CAS  Google Scholar 

  19. Böttcher C, Burbidge CA, Boss PK, Davies C. Interactions between ethylene and auxin are crucial to the control of grape (Vitis vinifera L.) berry ripening. BMC Plant Biol. 2013;13:222.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  20. Böttcher C, Harvey K, Forde CG, Boss PK, Davies C. Auxin treatment of pre-veraison grape (Vitis vinifera L.) berries both delays ripening and increases the synchronicity of sugar accumulation. Aust J Grape Wine Res. 2011;17(1):1–8.

    Article  CAS  Google Scholar 

  21. Ziliotto F, Corso M, Rizzini FM, Rasori A, Botton A, Bonghi C. Grape berry ripening delay induced by a pre-véraison NAA treatment is paralleled by a shift in the expression pattern of auxin- and ethylene-related genes. BMC Plant Biol. 2012;12:185.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  22. Seymour GB, Ostergaard L, Chapman NH, Knapp S, Martin C. Fruit development and ripening. Annu Rev Plant Biol. 2013;64:219–41.

    Article  CAS  PubMed  Google Scholar 

  23. Giovannoni J, Nguyen C, Ampofo B, Zhong S, Fei Z. The epigenome and transcriptional dynamics of fruit ripening. Annu Rev Plant Biol. 2017;68:61–84.

    Article  CAS  PubMed  Google Scholar 

  24. Karlova R, Chapman N, David K, Angenent GC, Seymour GB, de Maagd RA. Transcriptional control of fleshy fruit development and ripening. J Exp Bot. 2014;65(16):4527–41.

    Article  CAS  PubMed  Google Scholar 

  25. Wang R, Tavano E, Lammers M, Martinelli AP, Angenent GC, de Maagd RA. Re-evaluation of transcription factor function in tomato fruit development and ripening with CRISPR/Cas9-mutagenesis. Sci Rep. 2019;9(1):1696.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  26. Serrano A, Espinoza C, Armijo G, Inostroza-Blancheteau C, Poblete E, Meyer-Regueiro C, Arce A, Parada F, Santibanez C, Arce-Johnson P. Omics approaches for understanding grapevine berry development: regulatory networks associated with endogenous processes and environmental responses. Front Plant Sci. 2017;8:1486.

    Article  PubMed  PubMed Central  Google Scholar 

  27. Ma Q, Yang J. Transcriptome profiling and identification of the functional genes involved in berry development and ripening in Vitis vinifera. Gene. 2019;680:84–96.

    Article  CAS  PubMed  Google Scholar 

  28. Fasoli M, Richter CL, Zenoni S, Bertini E, Vitulo N, Dal Santo S, Dokoozlian N, Pezzotti M, Tornielli GB. Timing and order of the molecular events marking the onset of berry ripening in grapevine. Plant Physiol. 2018;178(3):1187–206.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  29. Sun R, He F, Lan Y, Xing R, Liu R, Pan Q, Wang J, Duan C. Transcriptome comparison of Cabernet Sauvignon grape berries from two regions with distinct climate. J Plant Physiol. 2015;178:43–54.

    Article  CAS  PubMed  Google Scholar 

  30. Zenoni S, Ferrarini A, Giacomelli E, Xumerle L, Fasoli M, Malerba G, Bellin D, Pezzotti M, Delledonne M. Characterization of transcriptional complexity during berry development in Vitis vinifera using RNA-Seq. Plant Physiol. 2010;152(4):1787–95.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  31. Xu L, Yue Q, Xiang G, Bian F, Yao Y. Melatonin promotes ripening of grape berry via increasing the levels of ABA, H2O2, and particularly ethylene. Hortic Res. 2018;5:41.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  32. Griesser M, Martinez SC, Eitle MW, Warth B, Andre CM, Schuhmacher R, Forneck A. The ripening disorder berry shrivel affects anthocyanin biosynthesis and sugar metabolism in Zweigelt grape berries. Planta. 2018;247(2):471–81.

    Article  CAS  PubMed  Google Scholar 

  33. Zarrouk O, Brunetti C, Egipto R, Pinheiro C, Genebra T, Gori A, Lopes CM, Tattini M, Chaves MM. Grape ripening is regulated by deficit irrigation/elevated temperatures according to cluster position in the canopy. Front Plant Sci. 2016;7:1640.

    Article  PubMed  PubMed Central  Google Scholar 

  34. Massonnet M, Fasoli M, Tornielli GB, Altieri M, Sandri M, Zuccolotto P, Paci P, Gardiman M, Zenoni S, Pezzotti M. Ripening transcriptomic program in red and white grapevine varieties correlates with berry skin anthocyanin accumulation. Plant Physiol. 2017;174(4):2376–96.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  35. Maul E, Sudharma KN, Ganesch A, Brühl U, Hundemer M, Kecke S, Mahler-Ries A, Marx G, Schreiber T, Walk M et al. 30 years VIVC - Vitis International Variety Catalogue (www.vivc.de). In: XI International Conference on Grapevine Breeding and Genetics; Yanqing, Beijing, China. 2014.

  36. Robinson J, Harding J, Vouillamoz J. Wine grapes: a complete guide to 1,368 vine varieties, including their origins and flavours. London: Penguin Books; 2012.

    Google Scholar 

  37. Roy-Chevrier J. Ampelographie retrospective. Paris;1900.

  38. Pucker B, Schwandner A, Becker S, Hausmann L, Viehover P, Topfer R, Weisshaar B, Holtgrawe D. RNA-Seq time series of vitis vinifera bud development reveals correlation of expression patterns with the local temperature profile. Plants. 2020;9(11):1548.

    Article  CAS  PubMed Central  Google Scholar 

  39. Dal Santo S, Tornielli GB, Zenoni S, Fasoli M, Farina L, Anesi A, Guzzo F, Delledonne M, Pezzotti M. The plasticity of the grapevine berry transcriptome. Genome Biol. 2013;14(6):R54.

    Article  CAS  Google Scholar 

  40. Vannozzi A, Dry IB, Fasoli M, Zenoni S, Lucchin M. Genome-wide analysis of the grapevine stilbene synthase multigenic family: genomic organization and expression profiles upon biotic and abiotic stresses. BMC Plant Biol. 2012;12:130.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  41. Palumbo MC, Zenoni S, Fasoli M, Massonnet M, Farina L, Castiglione F, Pezzotti M, Paci P. Integrated network analysis identifies fight-club nodes as a class of hubs encompassing key putative switch genes that induce major transcriptome reprogramming during grapevine development. Plant Cell. 2014;26(12):4617–35.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  42. Burger AL, Zwiegelaar JP, Botha FC. Characterisation of the gene encoding the Merlot ripening-induced protein 1 (mrip1): evidence that this putative protein is a distinct member of the plant proline-rich protein family. Plant Sci. 2004;167(5):1075–89.

    Article  CAS  Google Scholar 

  43. Davies C, Robinson SP. Differential screening indicates a dramatic change in mRNA profiles during grape berry ripening. Cloning and characterization of cDNAs encoding putative cell wall and stress response proteins. Plant Physiol. 2000;122(3):803–12.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  44. Qi T, Huang H, Song S, Xie D. Regulation of jasmonate-mediated stamen development and seed production by a bHLH-MYB complex in arabidopsis. Plant Cell. 2015;27(6):1620–33.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  45. Czemmel S, Stracke R, Weisshaar B, Cordon N, Harris NN, Walker AR, Robinson SP, Bogs J. The grapevine R2R3-MYB transcription factor VvMYBF1 regulates flavonol synthesis in developing grape berries. Plant Physiol. 2009;151(3):1513–30.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  46. Bogs J, Downey MO, Harvey JS, Ashton AR, Tanner GJ, Robinson SP. Proanthocyanidin synthesis and expression of genes encoding leucoanthocyanidin reductase and anthocyanidin reductase in developing grape berries and grapevine leaves. Plant Physiol. 2005;139(2):652–63.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  47. Bogs J, Jaffe FW, Takos AM, Walker AR, Robinson SP. The grapevine transcription factor VvMYBPA1 regulates proanthocyanidin synthesis during fruit development. Plant Physiol. 2007;143(3):1347–61.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  48. Koyama K, Numata M, Nakajima I, Goto-Yamamoto N, Matsumura H, Tanaka N. Functional characterization of a new grapevine MYB transcription factor and regulation of proanthocyanidin biosynthesis in grapes. J Exp Bot. 2014;65(15):4433–49.

    Article  CAS  PubMed  Google Scholar 

  49. Cavallini E, Matus JT, Finezzo L, Zenoni S, Loyola R, Guzzo F, Schlechter R, Ageorges A, Arce-Johnson P, Tornielli GB. The phenylpropanoid pathway is controlled at different branches by a set of R2R3-MYB C2 repressors in grapevine. Plant Physiol. 2015;167(4):1448–70.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  50. Matus JT, Cavallini E, Loyola R, Holl J, Finezzo L, Dal Santo S, Vialet S, Commisso M, Roman F, Schubert A, et al. A group of grapevine MYBA transcription factors located in chromosome 14 control anthocyanin synthesis in vegetative organs with different specificities compared with the berry color locus. Plant J. 2017;91(2):220–36.

    Article  CAS  PubMed  Google Scholar 

  51. Vannozzi A, Wong DCJ, Holl J, Hmmam I, Matus JT, Bogs J, Ziegler T, Dry I, Barcaccia G, Lucchin M. Combinatorial regulation of stilbene synthase genes by WRKY and MYB transcription factors in grapevine (Vitis vinifera L.). Plant Cell Physiol. 2018;59(5):1043–59.

    Article  CAS  PubMed  Google Scholar 

  52. Forlani S, Masiero S, Mizzotti C. Fruit ripening: the role of hormones, cell wall modifications, and their relationship with pathogens. J Exp Bot. 2019;70(11):2993–3006.

    Article  CAS  PubMed  Google Scholar 

  53. Dal Santo S, Vannozzi A, Tornielli GB, Fasoli M, Venturini L, Pezzotti M, Zenoni S. Genome-wide analysis of the expansin gene superfamily reveals grapevine-specific structural and functional characteristics. PLoS One. 2013;8(4):e62206.

    Article  CAS  Google Scholar 

  54. Du M, Spalding EP, Gray WM. Rapid auxin-mediated cell expansion. Annu Rev Plant Biol. 2020;71:379–402.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  55. Guillaumie S, Fouquet R, Kappel C, Camps C, Terrier N, Moncomble D, Dunlevy JD, Davies C, Boss PK, Delrot S. Transcriptional analysis of late ripening stages of grapevine berry. BMC Plant Biol. 2011;11:165.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  56. Li Z, Zhang C, Guo Y, Niu W, Wang Y, Xu Y. Evolution and expression analysis reveal the potential role of the HD-Zip gene family in regulation of embryo abortion in grapes (Vitis vinifera L.). BMC Genomics. 2017;18(1):744.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  57. Cui M, Wang W, Guo F, Fan X, Guan L, Zheng T, Zhu X, Jia H, Fang J, Wang C, et al. Characterization and temporal–spatial expression analysis of LEC1 gene in the development of seedless berries in grape induced by gibberellin. Plant Growth Regul. 2020;90(3):585–96.

    Article  CAS  Google Scholar 

  58. Jo L, Pelletier JM, Harada JJ. Central role of the LEAFY COTYLEDON1 transcription factor in seed development. J Integr Plant Biol. 2019;61(5):564–80.

    Article  CAS  PubMed  Google Scholar 

  59. Liu Q, Sun Y, Su W, Yang J, Liu X, Wang Y, Wang F, Li H, Li X. Species-specific size expansion and molecular evolution of the oleosins in angiosperms. Gene. 2012;509(2):247–57.

    Article  CAS  PubMed  Google Scholar 

  60. Qian X, Liu Y, Zhang G, Yan A, Wang H, Wang X, Pan Q, Xu H, Sun L, Zhu B. Alcohol acyltransferase gene and ester precursors differentiate composition of volatile esters in three interspecific hybrids of Vitis labruscaxV. Vinifera during berry development period. Food Chem. 2019;295:234–46.

    Article  CAS  PubMed  Google Scholar 

  61. Liang X, Zhou JM. Receptor-like cytoplasmic kinases: central players in plant receptor kinase-mediated signaling. Annu Rev Plant Biol. 2018;69:267–99.

    Article  CAS  PubMed  Google Scholar 

  62. Bourdais G, Burdiak P, Gauthier A, Nitsch L, Salojarvi J, Rayapuram C, Idanheimo N, Hunter K, Kimura S, Merilo E, et al. Large-scale phenomics identifies primary and fine-tuning roles for CRKs in responses related to oxidative stress. PLoS Genetics. 2015;11(7):e1005373.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  63. Dallakyan M, Zyprian E, Maul E, Yesayan A, Hovhannisyan N. Genetic characterization of rape varieties in Armenia. Vitis. 2015;54:23–6.

    Google Scholar 

  64. Schroeder A, Mueller O, Stocker S, Salowsky R, Leiber M, Gassmann M, Lightfoot S, Menzel W, Granzow M, Ragg T. The RIN: an RNA integrity number for assigning integrity values to RNA measurements. BMC Mol Biol. 2006;7:3.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  65. Bolger AM, Lohse M, Usadel B. Trimmomatic: a flexible trimmer for Illumina sequence data. Bioinformatics. 2014;30(15):2114–20.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  66. FastQC, a quality control tool for high throughput sequence data. [http://www.bioinformatics.babraham.ac.uk/projects/fastqc/]. Accessed 16 Nov 2018.

  67. Jaillon O, Aury JM, 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.

    Article  CAS  PubMed  Google Scholar 

  68. Kim D, Langmead B, Salzberg SL. HISAT: a fast spliced aligner with low memory requirements. Nat Methods. 2015;12(4):357–60.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  69. Kim D, Paggi JM, Park C, Bennett C, Salzberg SL. Graph-based genome alignment and genotyping with HISAT2 and HISAT-genotype. Nat Biotechnol. 2019;37(8):907–15.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  70. Liao Y, Smyth GK, Shi W. featureCounts: an efficient general purpose program for assigning sequence reads to genomic features. Bioinformatics. 2014;30(7):923–30.

    Article  CAS  PubMed  Google Scholar 

  71. Zhao S, Ye Z, Stanton R. Misuse of RPKM or TPM normalization when comparing across samples and sequencing protocols. RNA. 2020;26(8):903–9.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  72. Virtanen P, Gommers R, Oliphant TE, Haberland M, Reddy T, Cournapeau D, Burovski E, Peterson P, Weckesser W, Bright J, et al. SciPy 1.0: fundamental algorithms for scientific computing in Python. Nat Methods. 2020;17(3):261–72.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  73. 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.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  74. R: A language and environment for statistical computing [https://www.r-project.org]. Accessed  29 Nov 2018.

  75. Thimm O, Bläsing O, Gibon Y, Nagel A, Meyer S, Krüger P, Selbig J, Müller LA, Rhee SY, Stitt M. MAPMAN: a user-driven tool to display genomics data sets onto diagrams of metabolic pathways and other biological processes. Plant J. 2004;37(6):914–39.

    Article  CAS  PubMed  Google Scholar 

  76. Lohse M, Nagel A, Herter T, May P, Schroda M, Zrenner R, Tohge T, Fernie AR, Stitt M, Usadel B. Mercator: a fast and simple web server for genome scale functional annotation of plant sequence data. Plant Cell Environ. 2014;37(5):1250–8.

    Article  CAS  PubMed  Google Scholar 

  77. Pruitt KD, Tatusova T, Maglott DR. NCBI reference sequences (RefSeq): a curated non-redundant sequence database of genomes, transcripts and proteins. Nucleic Acids Res. 2007;35(Database issue):D61-65.

    Article  CAS  PubMed  Google Scholar 

  78. Altschul SF, Gish W, Miller W, Myers EW, Lipman DJ. Basic local alignment search tool. J Mol Biol. 1990;215(3):403–10.

    Article  CAS  PubMed  Google Scholar 

  79. Alexa A, Rahnenfuhrer J, Lengauer T. Improved scoring of functional groups from gene expression data by decorrelating GO graph structure. Bioinformatics. 2006;22(13):1600–7.

    Article  CAS  PubMed  Google Scholar 

  80. Abu-Jamous B, Kelly S. Clust: automatic extraction of optimal co-expressed gene clusters from gene expression data. Genome Biol. 2018;19:172.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  81. Zenoni S, D’Incà E, Tornielli GB. Genetic dissection of grape berry ripening control: defining a role for NAC transcription factors. Acta Hortic. 2019;1248:387–402.

    Article  Google Scholar 

  82. Wei L, Cao Y, Cheng J, Xiang J, Shen B, Wu J. Comparative transcriptome analyses of a table grape ‘Summer Black’ and its early-ripening mutant ‘Tiangong Moyu’ identify candidate genes potentially involved in berry development and ripening. J Plant Interact. 2020;15(1):213–22.

    Article  CAS  Google Scholar 

  83. Sievert C. Interactive web-based data visualization with R, plotly, and shiny. Boca Raton: CRC Press; 2020.

    Book  Google Scholar 

  84. Osier MV. VitisPathway: gene pathway analysis for V. vinifera . Vitis. 2016;55:129–33.

    Google Scholar 

  85. Grimplet J, Van Hemert J, Carbonell-Bejerano P, Díaz-Riquelme J, Dickerson J, Fennell A, Pezzotti M, Martínez-Zapater JM. Comparative analysis of grapevine whole-genome gene predictions, functional annotation, categorization and integration of the predicted gene sequences. BMC Research Notes. 2012;5:213.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  86. pheatmap: Pretty Heatmaps. R package version 1.0.12 [https://CRAN.r-project.org/package=pheatmap]. Accessed 16 Dec 2019.

  87. Zhang P, Berardini TZ, Ebert D, Li Q, Mi H, Muruganujan A, Prithvi T, Reiser L, Sawant S, Thomas PD, et al. PhyloGenes: an online phylogenetics and functional genomics resource for plant gene function inference. Plant Direct. 2020;4(12):e00293.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

Download references

Acknowledgements

We would like to thank Boas Pucker for providing several python scripts for handling large datasets as well as helpful comments on the manuscript, Katharina Frey-Sielemann for help visualizing the DEG distributions and Wiebke Halpape for providing the GO annotation enquiry. We also thank Kyle J. Lauersen for higly valuable language editing, and Willy Keller for the excellent technical support. We are very grateful to the Bioinformatics Resource Facility support team of the CeBiTec for providing computing infrastructure and excellent technical support. This publication is based in part upon contributions from COST Action CA 17111 INTEGRAPE, supported by COST (European Cooperation in Science and Technology); the authors thank BMBF/PtJ for funding in the context of IPAS (project acronym novisys). The authors from UniBi also wish to thank the other members of the Chair of Genetics and Genomics of Plants for their support. We would also like to thank the Open Access Publication Funds of Bielefeld University for funding the publication charge.

Funding

Open Access funding enabled and organized by Projekt DEAL. This project was funded by BMBF/PtJ (FKZ 031A3496).

Author information

Authors and Affiliations

Authors

Contributions

Sampling and phenotyping were done at the JKI (Julius Kühn-Institute, Institute for Grapevine Breeding Geilweilerhof, Siebeldingen, Germany). Sequencing and data analysis were performed at Bielefeld University, Faculty of Biology & Center for Biotechnology (CeBiTec). Organization, supervision at Bielefeld University and parts of sampling were done by DH. KH, AK and FS coordinated weekly sampling, phenotyping, and image capture for documentation. Logistic work was done by KH. The design of the experiments was set up by FS, DH and KH. RNA isolation and cDNA synthesis were carried out by DH and PV. PV accomplished library preparation and sequencing. Quantitative RT-PCR was performed by JT and DH. RT and KH supervised the work at JKI Geilweilerhof. LH performed the SSR analysis. BW supervised the work at Bielefeld University. FS, KH, RT and BW acquired project funding and wrote the project proposal. All bioinformatic data analyses, creation of figures, tables and writing of the manuscript were performed by JT with the help of DH. JT and BW drafted the manuscript. All authors have read and approved the final manuscript.

Corresponding author

Correspondence to Bernd Weisshaar.

Ethics declarations

Ethics approval and consent to participate

Not applicable

Consent for publication

Not applicable

Competing interests

The authors declare that they have no competing interests.

Additional information

Publisher's Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Supplementary Information

Additional file 1: Table S1.

Clonal relation of PN and PNP confirmed by 24 SSR markers. Table S2. Sampling and phenotypical observations (BBCH) for 2014. Table S3. Sampling and phenotypical observations (BBCH) for 2017. Table S4. Pearsonr correlations of gene expression values over time between both years. Table S5. Top 500 genes influencing principal component PC1. Table S6. Top 500 genes influencing principal component PC2. Table S7. Cluster memberships of clusters obtained from the cluster analysis PN/PNP. Table S8. Cluster memberships of clusters obtained from the cluster analysis PN. Table S9. Cluster memberships of clusters obtained from the cluster analysis PNP. Table S10. GO term enrichment for the clusters from the cluster analysis of PN/PNP. Table S11. GO term enrichment for the clusters from the cluster analysis of PN. Table S12. GO term enrichment for the clusters from the cluster analysis of PNP. Table S13. Differentially expressed genes between PN and PNP in both years and filtering steps. Table S14. Overlap of 12 putatively ripening time control DEGs (time points: DAF 0-35 from 2014 intersecting with DAF 0-42 from 2017) and 388 veraison-specific DEGs (time points: DAF35+41 from 2014 intersecting with DAF42+49 from 2017) with previous transcriptomic grapevine studies. Table S15. VitisNet enrichment for intersecting DEGs. Table S16. VitisNet enrichment for regulatory DEGs. Table S17. VitisNet enrichment for veraison-specific DEGs. Table S18. Up- & Down-regulated TFs among the DEGs. Table S19. ENA sample identifier and metadata linked to the study for PNP RNA-Seq reads from 2014 (PRJEB39262). Table S20. ENA sample identifier and metadata linked to the study for PN RNA-Seq reads from 2014 (PRJEB39261). Table S21. ENA sample identifier and metadata linked to the study for PNP RNA-Seq reads from 2017 (PRJEB39264). Table S22. ENA sample identifier and metadata linked to the study for PN RNA-Seq reads from 2017 (PRJEB39263). Table S23. Sequences of primers used for qRT-PCR.

Additional file 2: Figure S1.

Cluster analyses of PN/PNP in both years (A), PN in both years (B) and PNP in both years (C). Figure S2. Heatmap of all MYB-TFs gene expression among the DEGs for PN and PNP in both years. Figure S3. Expression patterns of genes from the ripening-regulatory gene set. Figure S4. Confirmation of expression data from RNA-Seq by qRT-PCR

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Theine, J., Holtgräwe, D., Herzog, K. et al. Transcriptomic analysis of temporal shifts in berry development between two grapevine cultivars of the Pinot family reveals potential genes controlling ripening time. BMC Plant Biol 21, 327 (2021). https://doi.org/10.1186/s12870-021-03110-6

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doi.org/10.1186/s12870-021-03110-6

Keywords