RNA-Seq Time Series of Vitis vinifera Bud Development Reveals Correlation of Expression Patterns with the Local Temperature Profile

Plants display sophisticated mechanisms to tolerate challenging environmental conditions and need to manage their ontogenesis in parallel. Here, we set out to generate an RNA-Seq time series dataset throughout grapevine (Vitis vinifera) early bud development. The expression of the developmental regulator VviAP1 served as an indicator of the progression of development. We investigated the impact of changing temperatures on gene expression levels during the time series and detected a correlation between increased temperatures and a high expression level of genes encoding heat-shock proteins. The dataset also allowed the exemplary investigation of expression patterns of genes from three transcription factor (TF) gene families, namely MADS-box, WRKY, and R2R3-MYB genes. Inspection of the expression profiles from all three TF gene families indicated that a switch in the developmental program takes place in July which coincides with increased expression of the bud dormancy marker gene VviDRM1.


Introduction
Plants are sessile organisms that cannot escape from herbivores or changes in environmental conditions. As a consequence, various stress response mechanisms [1,2] and complex specialized metabolite pathways evolved to counteract adverse situations and conditions [3,4]. These response mechanisms also need to protect plant embryogenesis, as well as vegetative developmental processes like an outgrowth of side shoots from resting or latent buds. All these developmental processes, including the establishment of axillary buds, have to proceed undisturbed despite potentially challenging and/or unfavorable environmental conditions. Similar to other woody perennial plants like, e.g., apple or poplar, grapevine (V. vinifera) bud development spans over two years between bud initiation and growth of new side shoots. Newly formed buds enter a dormancy phase in the wintertime between the two growing seasons before buds sprout in the second season [5,6]. In the spring of the first season (April/May on the northern hemisphere), new axillary buds are formed on young grapevine shoots. These new buds initially contain meristems that develop into embryonic shoots with their shoot apical meristems (SAM) and containing primordia for leaves, tendrils and inflorescences. This implies that different types of meristems, including lateral and inflorescence meristems, co-exist in the buds. The floral transition takes place at about June, five to seven weeks after the burst of "old" buds (i.e., the buds that are one year ahead in development). Inflorescence primordia differentiate from uncommitted primordia formed within the new buds. Due to further differentiation of inflorescence meristems into inflorescence branch meristems (about July), the compound buds finally contain the embryonic version of next year's shoots, each with tissues for first leaves, inflorescences and tendrils [6,7]. The buds enter endodormancy which passes over into ecodormancy depending on the environmental conditions of fall and winter [8][9][10]. In the early spring of the second season, ecodormancy is released and inflorescence branch meristems produce single flower meristems in swelling buds (April) and flower organ development begins [6]. It is important to note that the precise timing of floral transition and development including the establishment of phenological differences [11] strongly depend on environmental conditions and genotype [12].
Heat-shock proteins (HSPs) are a group of proteins, which were initially detected due to their accumulation in response to quickly increased temperature. HSPs are central parts of a molecular mechanism to endure higher temperatures. First reports of HSPs in plants reach back to the 1980s when they were described based on cell culture experiments with tobacco and soybean [13]. HSPs are assumed to support several physiological functions under normal growth conditions. This includes folding, unfolding, localization, accumulation and degradation of other proteins [14,15]. Additionally, irreversible aggregation of other proteins is prevented and refolding is facilitated under heat stress [16]. Several categories of HSPs based on sequence homology and typical molecular weight have been defined [14], thus leading to multiple polyphyletic groups of HSPs.
WRKY transcription factors (TFs) are a family of TFs, which play an important role in the regulation of responses to environmental stress conditions [17][18][19]. R2R3-MYB TFs are often responsible for controlling the formation of specialized metabolites in response to environmental triggers, but also regulate several plant-specific processes including root hair and trichome differentiation [20][21][22]. MADS-box TFs are typically involved in the regulation of developmental processes like the determination of plant organ identity [23][24][25]. One especially important developmental regulator is APETALA1 (AP1), also a MADS-box factor, which connects signals received from the environment with initiation and/or progress of developmental processes [26,27]. VviAP1 and VviAIL2, a V. vinifera homolog of the MADS-box gene AINTEGUMENTA-like (AtAIL1, At1g72570), have been postulated to be involved in the photoperiodic control of seasonal growth [28]. In addition, marker genes for the dormant state of buds have been described. One such marker gene is DRM1, a gene that has been found initially in Pisum sativum to encode a dormancy-associated protein [29]. Subsequently, DRM1 homologs have been identified in many species in the context of bud dormancy, including V. vinifera [9,30].
In the model plant Arabidopsis thaliana, phylotranscriptomic evidence for a molecular embryonic hourglass was published [31][32][33]. We attempted to create an RNA-Seq dataset to examine the early bud development of V. vinifera for a similar general pattern. While an hourglass pattern was not detected in the data (Marcel Quint, Martin-Luther-University Halle-Wittenberg, Germany, personal communication), we harnessed the time series of V. vinifera RNA-Seq samples to investigate changes in gene expression during early bud development throughout the first season at a fine scale and observed a strong influence of high temperatures on the expression of HSP genes of field-grown plants. In addition, a switch in the expression patterns of various TF genes was observed that happens in parallel to or shortly after the switch from uncommitted primordia to inflorescence primordia. This switch in expression pattern coincides with the onset of expression of the dormancy marker gene VviDRM1.

RNA-Seq Time Series of Early Bud Development and Transcript Accumulation Patterns of Selected Marker Genes
Young buds of V. vinifera 'Calardis Musqué' were harvested over a period of 156 days of the first season of development ( Figure 1, File S1). RNA from sampled buds was subjected to RNA-Seq analyses in triplicates. Values for gene expression, inferred from values for transcript abundance in the triplicates, were calculated for all genes (Files S2-S4; see Methods for details).

RNA-Seq Time Series of Early Bud Development and Transcript Accumulation Patterns of Selected Marker Genes
Young buds of V. vinifera 'Calardis Musqué' were harvested over a period of 156 days of the first season of development ( Figure 1, File S1). RNA from sampled buds was subjected to RNA-Seq analyses in triplicates. Values for gene expression, inferred from values for transcript abundance in the triplicates, were calculated for all genes (Files S2-S4; see Methods for details). Figure 1. Documentation of the grapevine material used for sampling. The upper row of pictures displays the bud stages during the first year of development, which was the target of this study. The lower row of pictures displays the growth status of the vines from which the young buds were taken. Pictures were taken in parallel to harvest. In the center (grey bar), the sampling timeline is depicted. Above and below the pictures, the main developmental phases are mentioned, with the BBCH stages (see Methods) at the bottom for the second year of the two-year reproductive cycle of grapevine.
The expression pattern of the MADS-Box genes VviAP1 (CRIBI2.1 ID VIT_201s0011g00100), VviAIL2 (VIT_209s0002g01370), and VviSOC1a (VIT_215s0048g01250) are displayed in Figure 2a. VviAP1 transcript levels were zero or very low until the end of June and rise until September. Due to the time distance between the sampling points, we interpret the data as essentially one peak in September. VviAP1 and VviAIL2 display quite similar expression patterns. The increase in VviAP1 transcript levels at the end of June coincides with the time when floral transition, the differentiation of uncommitted primordia into inflorescence primordia, took place. There is no overall correlation between the transcript levels of VviAP1 and VviAIL2 and day length. However, the rise of transcript levels coincides with the beginning of the reduction of day length after midsummer. We also checked the expression of the dormancy marker gene VviDRM1 (VIT_210s0003g00090) and found high transcript levels of this gene in the buds with a clear increase starting in July (Figure 2b). The lower row of pictures displays the growth status of the vines from which the young buds were taken. Pictures were taken in parallel to harvest. In the center (grey bar), the sampling timeline is depicted. Above and below the pictures, the main developmental phases are mentioned, with the BBCH stages (see Methods) at the bottom for the second year of the two-year reproductive cycle of grapevine.
The expression pattern of the MADS-Box genes VviAP1 (CRIBI2.1 ID VIT_201s0011g00100), VviAIL2 (VIT_209s0002g01370), and VviSOC1a (VIT_215s0048g01250) are displayed in Figure 2a. VviAP1 transcript levels were zero or very low until the end of June and rise until September. Due to the time distance between the sampling points, we interpret the data as essentially one peak in September. VviAP1 and VviAIL2 display quite similar expression patterns. The increase in VviAP1 transcript levels at the end of June coincides with the time when floral transition, the differentiation of uncommitted primordia into inflorescence primordia, took place. There is no overall correlation between the transcript levels of VviAP1 and VviAIL2 and day length. However, the rise of transcript levels coincides with the beginning of the reduction of day length after midsummer. We also checked the expression of the dormancy marker gene VviDRM1 (VIT_210s0003g00090) and found high transcript levels of this gene in the buds with a clear increase starting in July (Figure 2b).

Average Gene Expression Values of HSP Genes Reflect the Local Temperature Profile
We made use of the weather data recorded at the vineyard from which the samples for RNA-Seq were derived. Since the recorded temperature profile during 2016 displayed significant oscillation, we tested the hypothesis that a heat-shock response may take place in the buds. A total of 131 putative HSP genes were identified based on annotation data (File S5), but only 81 of these show substantial transcript abundance (see Methods). The averaged gene expression pattern of this set of 81 HSP genes (median values) shows a good overall correlation (r = 0.51, p = 0.0001; Spearman rank correlation coefficient) with the temperature profile. We observed clear HSP gene transcript level peaks at time points with high temperatures (Figure 3; see File S5 for the correlation values of individual genes, a plot with the arithmetic mean of expression values and the expression pattern of a single HSP gene). This is especially noticeable on 24 June when the highest average temperature was recorded. A strong correlation of the day length/photoperiod with the expression pattern of these genes was not observed.

Average Gene Expression Values of HSP Genes Reflect the Local Temperature Profile
We made use of the weather data recorded at the vineyard from which the samples for RNA-Seq were derived. Since the recorded temperature profile during 2016 displayed significant oscillation, we tested the hypothesis that a heat-shock response may take place in the buds. A total of 131 putative HSP genes were identified based on annotation data (File S5), but only 81 of these show substantial transcript abundance (see Methods). The averaged gene expression pattern of this set of 81 HSP genes (median values) shows a good overall correlation (r = 0.51, p = 0.0001; Spearman rank correlation coefficient) with the temperature profile. We observed clear HSP gene transcript level peaks at time points with high temperatures (Figure 3; see File S5 for the correlation values of individual genes, a plot with the arithmetic mean of expression values and the expression pattern of a single HSP gene). This is especially noticeable on 24 June when the highest average temperature was recorded. A strong correlation of the day length/photoperiod with the expression pattern of these genes was not observed.

Investigation of Transcription Factor Gene Families: WKRY, MADS-Box, and R2R3-MYBs
We harnessed the presented RNA-Seq time series for the analyses of expression patterns of three TF gene families. Heatmaps display transcript levels of genes encoding MADS-box (File S6), R2R3-MYB (File S7), and WRKY (File S8) TFs. Only genes that display detectable transcript accumulation values were considered (see Methods for the threshold). As can be seen from all three heatmaps, the gene activity patterns of quite some of the transcription factor genes change quite dramatically with the onset of VviAP1 transcript accumulation between 28 June and 26 July.
While VviSVP1 (VIT_200s0313g00070), an ortholog of the A. thaliana MADS-box gene SHORT VEGETATIVE PHASE/AGL22 (At2g22540), shows transcript levels with almost constant values, the gene VviFLC2 (VIT_214s0068g01800) displays a time course quite similar to those of VviAP1 and VviAIL2. VviFLC2 is, like its paralog VviFLC1 (VIT_201s0010g03890), closely related to the A. thaliana MADS-box gene FLOWERING LOCUS C/AGL25 (At5g10140) which encodes a central repressor of the floral transition. In contrast, VviTM8a (VIT_217s0000g01230), which was named according to a gene initially detected in Solanum lycopersicum (TOMATO MADS 8) that became "founder" of a specific sub-clade of evolutionary related MADS-box genes, shows a transcript accumulation peak at the end of June. Finally, VviSOC1a (VIT_215s0048g01250), a homolog of the A. thaliana MADS-box gene SUPPRESSOR OF OVEREXPRESSION OF CO 1/AGL20 (At2g45660) displays transcript accumulation values that decline after middle of August (Figure 2a).
With regard to the R2R3-MYB genes, VviMYBC2-L1 (VIT_201s0011g04760), VviMYB4A (VIT_203s0038g02310), and VviMYBPAR (VIT_211s0016g01300) are prominent examples of genes that support the gene expression pattern switch during July. In addition to the July switch, several R2R3-MYB genes display high transcript accumulation values specifically in November. Clear examples are VviMYB15 (VIT_205s0049g01020), VviMYB14 (VIT_205s0049g01020), and VviMYB30A (VIT_217s0000g06190). Based on its high transcript levels, VviMYBPA1 (VIT_215s0046g00170), a homolog of the A. thaliana R2R3-MYB genes TRANSPARENT TESTA 2/AtMYB123 (At5g35550) and AtMYB5 (At3g13540), appears as an important regulator. It is worth noting that another homolog of AtMYB123 and AtMYB5, namely VviMYBPAR, also shows high transcript levels from end of August to September. Another R2R3-MYB gene that stands out due to high transcript levels is VviMYB174 (VIT_218s0001g09850), a homolog of AtMYB73 (At4g37260) and AtMYB77 (At3g50060). As expected for organs/tissues not accumulating anthocyanins, VviMYBA1 (VIT_202s0033g00410), an ortholog of

Investigation of Transcription Factor Gene Families: WKRY, MADS-Box, and R2R3-MYBs
We harnessed the presented RNA-Seq time series for the analyses of expression patterns of three TF gene families. Heatmaps display transcript levels of genes encoding MADS-box (File S6), R2R3-MYB (File S7), and WRKY (File S8) TFs. Only genes that display detectable transcript accumulation values were considered (see Methods for the threshold). As can be seen from all three heatmaps, the gene activity patterns of quite some of the transcription factor genes change quite dramatically with the onset of VviAP1 transcript accumulation between 28 June and 26 July.
While VviSVP1 (VIT_200s0313g00070), an ortholog of the A. thaliana MADS-box gene SHORT VEGETATIVE PHASE/AGL22 (At2g22540), shows transcript levels with almost constant values, the gene VviFLC2 (VIT_214s0068g01800) displays a time course quite similar to those of VviAP1 and VviAIL2. VviFLC2 is, like its paralog VviFLC1 (VIT_201s0010g03890), closely related to the A. thaliana MADS-box gene FLOWERING LOCUS C/AGL25 (At5g10140) which encodes a central repressor of the floral transition. In contrast, VviTM8a (VIT_217s0000g01230), which was named according to a gene initially detected in Solanum lycopersicum (TOMATO MADS 8) that became "founder" of a specific sub-clade of evolutionary related MADS-box genes, shows a transcript accumulation peak at the end of June. Finally, VviSOC1a (VIT_215s0048g01250), a homolog of the A. thaliana MADS-box gene SUPPRESSOR OF OVEREXPRESSION OF CO 1/AGL20 (At2g45660) displays transcript accumulation values that decline after middle of August (Figure 2a).
With regard to the R2R3-MYB genes, VviMYBC2-L1 (VIT_201s0011g04760), VviMYB4A (VIT_203s0038g02310), and VviMYBPAR (VIT_211s0016g01300) are prominent examples of genes that support the gene expression pattern switch during July. In addition to the July switch, several R2R3-MYB genes display high transcript accumulation values specifically in November. Clear examples are VviMYB15 (VIT_205s0049g01020), VviMYB14 (VIT_205s0049g01020), and VviMYB30A (VIT_217s0000g06190). Based on its high transcript levels, VviMYBPA1 (VIT_215s0046g00170), a homolog of the A. thaliana R2R3-MYB genes TRANSPARENT TESTA 2/AtMYB123 (At5g35550) and AtMYB5 (At3g13540), appears as an important regulator. It is worth noting that another homolog of AtMYB123 and AtMYB5, namely VviMYBPAR, also shows high transcript levels from end of August to September. Another R2R3-MYB gene that stands out due to high transcript levels is VviMYB174 (VIT_218s0001g09850), a homolog of AtMYB73 (At4g37260) and AtMYB77 (At3g50060). As expected for organs/tissues not accumulating anthocyanins, VviMYBA1 (VIT_202s0033g00410), an ortholog of the A. thaliana R2R3-MYB gene PRODUCTION OF ANTHOCYANIN PIGMENT 1/AtMYB75 (At1g56650), is not significantly expressed in young buds.
Several genes encoding TFs of the WRKY type show a substantial increase in transcript levels at the gene expression pattern switch during July (e.g., VviWRKY20 (VIT_207s0005g02570) or VviWRKY31 (VIT_210s0003g02810)). Examples for the opposite change, i.e., reduction of transcript levels during July, are VvWRKY23 (VIT_207s0031g01840) and VviWRKY41 (VIT_213s0067g03140). High transcript levels that are strongly reduced towards winter were detected for VviWRKY25 (VIT_208s0058g01390), that is, together with VviWRKY41, homologous to AtWRKY54 (At2g40750), AtWRKY70 (At3g56400) and AtWRKY46 (At2g46400). Like some of the R2R3-MYB genes, also several VviWRKY genes display high transcript accumulation values specifically in November and/or an increase towards winter. These include VviWRKY16 (VIT_206s0004g07500), VviWRKY45 (VIT_214s0108g01280) and VviWRKY33 (VIT_211s0037g00150) that are all homologous to the A. thaliana WRKY genes of group I-C including AtWRKY33 (At2g38470), AtWRKY58 (At3g01080), and AtWRKY32 (At4g30935).

Identification of qRT-PCR Reference Genes
The quite long RNA-seq time series from the tissue of compound buds, covering changing day length and oscillating weather conditions in the field, allowed the identification of candidate reference genes for quantitative Real-Time PCR (qRT-PCR) experiments. The 20 best candidates for reference genes in qRT-PCR experiments were predicted based on an overall high expression and a low variation in steady-state transcript abundance (see Methods, File S9). A manual inspection of the functional annotation of these genes supported the quality of this dataset since it covers well-known qRT-PCR reference genes encoding polyubiquitin, glyceraldehyde-3-phosphate dehydrogenase, elongation factor Tu, and actin as top candidates. The most promising candidate, i.e., position one of the list of candidate genes, is VviUBQ10 (VIT_219s0177g00040) which is homologous to the five A. thaliana genes encoding polyubiquitin (At4g05320 and others). The second position is occupied by VIT_219s0015g01090, a homolog of A. thaliana HEAT SHOCK PROTEIN 81.4 (At5g56000, see Discussion).

Discussion
This time series of 18 RNA-Seq data points throughout the first year of development of V. vinifera compound buds allows the investigation of gene expression patterns throughout early bud development at relatively high resolution. A similar time-series analysis was performed previously with Affymetrix arrays to determine gene expression patterns for the cultivar 'Tempranillo' [8]. While the data from 'Tempranillo' that was grown near Madrid cover the time from May (first year) to April (second year) with eight time points, the focus of the time series presented here was the early bud development in June until November of the first year. Nevertheless, the previously reported expression pattern changes in July and also towards winter (November time point) [8] were essentially matched by our dataset. The main difference of our data is, in addition to the digital quantification and wider dynamic range intrinsic to RNA-Seq when compared to array hybridisation, that there are 11 samples from the month June.
Our chronologically dense collection of samples allows the detection of small developmental differences between time points, but it is affected by a high variation through individual differences between plants, sampled buds, and varying environmental conditions like temperature and other weather conditions. Harvesting the buds in the field is technically challenging since the buds are small and need to be cut out of the axil between shoot and leaf or from hard wood. The sampling has to be performed quickly and there was no time to remove adhering tissue from shoots before freezing the samples in liquid nitrogen. The results presented were derived from data of the year 2016 and from time points for which three independent biological replicates were available. However, there are more data on time points for which individual replicates were lost mainly during the RNA extraction procedures due to the problematic technical properties of the respective samples (File S1). While these time points should not be used for statistical analyses, they still provide additional support for patterns observed or may increase the power of co-expression analyses in the future.
Predicted reference genes for qRT-PCR experiments of bud samples contain commonly used reference genes like GAPDH, actin, polyubiquitin, and elongation factor Tu. Additionally, novel candidate genes were identified which displayed a constant expression level. This aligns well with a previous study, which reported that novel reference genes identified by genome-wide in silico analysis outperformed typical reference genes in common wheat [34].
We selected three well-known and-in the model plant A. thaliana-well-characterised TF gene families for more detailed analyses of gene expression patterns. MADS-box genes were selected because of their relation to development, and also because quite some of these genes were analysed in the 'Tempranillo' study of Diaz-Riquelme et al. [8]. R2R3-MYB genes and WRKY genes were selected because of their link to stress responses as well as control of the accumulation of specialized metabolites.
The gene VviSOC1a, a potential integrator of multiple flowering signals that cumulate in the establishment of inflorescence meristems [35], is expressed in June and expression declines after middle of August when VviAP1 and also VviLFY/VviFL (VIT_217s0000g00150, see File S3; ortholog of AtLEAFY [36], At5g61850) show increasing expression levels. This increased expression of potential inflorescence meristem identity genes (AP1, LEAFY) coincides with the proliferation of inflorescence primordia, giving rise to inflorescence branch primordia in the developing compound buds [5,6,37]. Based also on the increasing expression of the dormancy marker gene VviDRM1, it can be postulated that other parts of the compound bud are already in July on their way to endodormancy. Additionally, the gene VviTM8a, a homolog of TOMATO MADS 8 that plays a role in tomato flower development [38], displays an interesting expression pattern that hints at inflorescence developmental processes taking place during June and the beginning of July.
With respect to the expression of the two R2R3-MYB factors known to control proanthocyanidin (PA) biosynthesis that displayed conspicuous expression patterns, namely VviMYBPA1 and VviMYBPAR [39,40], three prominent potential target genes VviLAR1 (VIT_201s0011g02960) VviANS (or VviLDOX, VIT_202s0025g04720) and VviANR (VIT_200s0361g00040) that all encode enzymes of PA biosynthesis show expression patterns expected for targets (File S3), also with a reduction of expression levels towards winter. The three genes encode the enzymes leucoanthocyanidin reductase (LAR), anthocyanidin synthase (ANS, also referred to as LDOX) and anthocyanidin reductase (ANR) which are required for biosynthesis of catechin and epicatechin that are precursors of PAs [40]. This indicates that the compound buds accumulate PAs during summer and fall in preparation for winter. Similarly, it is conceivable that the activation of VviMYB14 and VviMYB15 results in the synthesis of stilbenes [41,42], which fits the activation of three of the V. vinifera stilbene synthase genes [43] that are clustered on chromosome 16. While VviSTS35 (VIT_216s0100g01070) and VviSTS41 (VIT_216s0100g01130) might be targets of VviMYB15 based on strong co-expression in November, VviSTS36 (VIT_216s0100g01100) fits better as a potential target of VviMYB14. The analysis of VviSTS genes [43] also covered three VviCHS genes (VviCHS3: VIT_205s0136g00260, VviCHS1: VIT_214s0068g00920 and VviCHS2: VIT_214s0068g00930). The three VviCHS genes are co-expressed with a pattern similar to that of VviLAR1, VviANS and VviANR which fits the substrate requirement for catechin/epicatechin biosynthesis and to VviMYBPA1 and/or VviMYBPAR as potential regulators. Future studies should measure the level of phenylpropanoid compounds in grapevine buds to investigate the correlation of metabolite accumulation with gene expression of biosynthesis genes and regulators. The homologs of the highly expressed gene VviMYB174, AtMYB73 and AtMYB77, have been implicated in root-related auxin responses [44]. A similar function of VviMYB174 would fit its high expression throughout the compound but development from June to November.
Of the VviWRKY genes that have been implicated in the control of VviSTS's, namely VviWRKY03, VviWRKY24, VviWRKY43 and VviWRKY53 [42], only VviWRKY03 (VIT_201s0010g03930) displayed expression with a pattern that fits that of VviMYB14 and VviSTS36, allowing to hypothesize that STS36 expression is under combinatorial control of MYB14/WRKY03 in compound buds. VviWRKY25, and to some extent VviWRKY41 as well, that are both homologous to the AtWRKY genes implicated in brassinosteroid-regulated plant growth [45], show expression patterns that fit growth actions until September that are then abandoned towards November and winter.
To the best of our knowledge, this is the first report of a correlation between HSP gene expression patterns with the environmental temperature in a comprehensive time series of field (vineyard) samples. However, the repeated formation of HSPs was previously described in the seeds, seed pods, and flowers of Medicago sativa [46]. The occurrence of HSPs at standard (not stressed) growth conditions indicated a potential role of HSPs in development [46]. The annotation "heat-shock" was initially introduced based on upregulation of genes in heat stress experiments [13,47]. HSPs were also detected at substantial levels in field-grown Gossypium hirsutum under increased temperature and drought stress [48]. Reports from Oryza sativa support a stress signal integration function of HSPs [49]. Therefore, we also checked for other stress factors like documented pathogen attack, crop protection treatments, or drought stress, but the temperature was the only factor with a substantial correlation. The detection of a gene (VIT_219s0015g01090, see File S5) that encodes an HSP81.4 variant among the potential reference genes for qRT-PCR is explained by a high expression that is relatively constant but follows the temperature profile closely in a very limited fluctuation range. Considering the large number of physiological functions of HSPs, more or less constantly expressed chaperons, that are related in sequence to proteins which are really responsive (at the level of gene expression) to quickly increasing temperatures, might also be annotated as HSPs. This assumption aligns well with previous findings that HSPs can have functions in the integration of stress signals [50]. From checking single genes among the 81 well-expressed HSP genes (File S5), it is evident that the gene expression correlation of bona fide HSP genes with the temperature profile is even stronger than described for the gene set. Moreover, it is possible that additional factors like an underlying developmental pattern or UV-B exposure have an additional influence on the observed heat-shock gene expression profile.
We investigated grapevine early bud development and differentiation with respect to changes in the transcriptome under field conditions. We describe expression level coincidence of selected developmental marker genes as well as several well-known TF genes with developmental switches of V. vinifera bud development. More importantly, we found a correlation of expression levels of genes encoding HSPs with the local daily average temperature that is very obvious for genes like VIT_202s0025g00280 which encodes an HSP90 variant. Follow-up studies may focus on the protective function of HSP genes that display temperature responsiveness under field conditions, especially regarding the negative effect of increasing temperatures due to climate change that will probably cause serious damages to viticulture.  (Figure 1), bud samples were taken in triplicates. From three different vines, four buds each were harvested in a batch and immediately frozen in liquid nitrogen. The buds were taken from the fourth to eight node of the shoots emerging from the middle section of the cane (Figure 4). Vines that appeared to be equal in their overall developmental stage were chosen. Those that showed symptoms of nutrient deficiency or diseases were excluded as a sample source.

Biological Material
The growth stage of 'Calardis Musqué' on every sampling date was described using the BBCH scale [51]. The BBCH stage systematically classifies the developmental state of crop plants, including grapevine. If more than 60% of the organs of the vines (e.g., buds, leaves, inflorescences, berries) were classified to a certain BBCH stage, this stage was assigned to the respective vines and recorded together with the harvested samples. In parallel to harvest and throughout the season, pictures were taken from individual buds and organs relevant for BBCH stage determination.
Plants 2020, 9, x FOR PEER REVIEW 9 of 14 grapevine. If more than 60% of the organs of the vines (e.g., buds, leaves, inflorescences, berries) were classified to a certain BBCH stage, this stage was assigned to the respective vines and recorded together with the harvested samples. In parallel to harvest and throughout the season, pictures were taken from individual buds and organs relevant for BBCH stage determination.

RNA Extraction, Library Preparation, and Sequencing
Total RNA was extracted, from four buds each, in triplicate per time point. Up to 100 mg of liquid nitrogen ground tissue was applied to the Spectrum™ Plant Total RNA kit (Sigma-Aldrich, Taufkirchen, Germany) according to the manufacturer's instructions for protocol B. After on-column DNase treatment with the DNase I Digest Set (Sigma-Aldrich, Taufkirchen, Germany), the RNA was quantified. An amount of 500 ng total RNA per sample was used to prepare sequencing libraries according to the Illumina TruSeq RNA Sample Preparation v2 Guide. Purification of the polyAcontaining mRNA was performed using two rounds of oligo(dT) oligonucleotides attached to magnetic beads. During the second elution of the polyA+ RNA, the RNA was fragmented and primed for cDNA synthesis. After cDNA synthesis, the DNA fragments were end-repaired and A-tailing was performed. Multiple indexing adapters, specific for each library and sample, were ligated to the ends of the cDNA fragments and the adapter-ligated fragments were enriched by 12 cycles of PCR. After qualification and quantification, the resulting sequencing libraries were equimolarly pooled and sequenced generating 100 nt single-end reads on eight lanes of an Illumina HiSeq1500 flowcell at the Sequencing Core Facility of the Center for Biotechnology (CeBiTec) at Bielefeld University.

Bioinformatic Analysis of RNA-Seq Data
All RNA-Seq read datasets generated were submitted to the European Nucleotide Archive (ENA; for accession numbers see File S1, the ENA study ID is PRJEB35820). Time points with only two successful biological replicates were included in the submission to the ENA database (File S1), but excluded from the investigations presented here. Python scripts developed for customized analyses, including the generation of the plots in Figures 2a,b and 3, are available at Github: https://github.com/bpucker/vivi-bud-dev. RNA-Seq reads were mapped to the CRIBI2.1 reference genome sequence of PN40024 [52] via STAR v.2.51b [53] with previously optimized parameters including a minimal alignment length cutoff of 90% and a minimal similarity cutoff of 95% of the read length [54]. FeatureCounts v1.5.0-p3 [55] was deployed for quantification of steady-state transcript levels at the gene level based on these mappings and the CRIBI2.1 annotation [56]. Previously developed Python scripts [54] were applied to merge the resulting count tables and to calculate counts per million (CPMs) and reads per kb per million mapped reads (RPKMs). We

RNA Extraction, Library Preparation, and Sequencing
Total RNA was extracted, from four buds each, in triplicate per time point. Up to 100 mg of liquid nitrogen ground tissue was applied to the Spectrum™ Plant Total RNA kit (Sigma-Aldrich, Taufkirchen, Germany) according to the manufacturer's instructions for protocol B. After on-column DNase treatment with the DNase I Digest Set (Sigma-Aldrich, Taufkirchen, Germany), the RNA was quantified. An amount of 500 ng total RNA per sample was used to prepare sequencing libraries according to the Illumina TruSeq RNA Sample Preparation v2 Guide. Purification of the polyA-containing mRNA was performed using two rounds of oligo(dT) oligonucleotides attached to magnetic beads. During the second elution of the polyA+ RNA, the RNA was fragmented and primed for cDNA synthesis. After cDNA synthesis, the DNA fragments were end-repaired and A-tailing was performed. Multiple indexing adapters, specific for each library and sample, were ligated to the ends of the cDNA fragments and the adapter-ligated fragments were enriched by 12 cycles of PCR. After qualification and quantification, the resulting sequencing libraries were equimolarly pooled and sequenced generating 100 nt single-end reads on eight lanes of an Illumina HiSeq1500 flowcell at the Sequencing Core Facility of the Center for Biotechnology (CeBiTec) at Bielefeld University.

Bioinformatic Analysis of RNA-Seq Data
All RNA-Seq read datasets generated were submitted to the European Nucleotide Archive (ENA; for accession numbers see File S1, the ENA study ID is PRJEB35820). Time points with only two successful biological replicates were included in the submission to the ENA database (File S1), but excluded from the investigations presented here. Python scripts developed for customized analyses, including the generation of the plots in Figure 2a,b and Figure 3, are available at Github: https://github.com/bpucker/vivi-bud-dev. RNA-Seq reads were mapped to the CRIBI2.1 reference genome sequence of PN40024 [52] via STAR v.2.51b [53] with previously optimized parameters including a minimal alignment length cutoff of 90% and a minimal similarity cutoff of 95% of the read length [54]. FeatureCounts v1.5.0-p3 [55] was deployed for quantification of steady-state transcript levels at the gene level based on these mappings and the CRIBI2.1 annotation [56]. Previously developed Python scripts [54] were applied to merge the resulting count tables and to calculate counts per million (CPMs) and reads per kb per million mapped reads (RPKMs). We attempted to include VviFT (GSVIVT00012870001 in the Vv8x genome sequence) in the analyses of selected target genes, but the corresponding sequence region is not included in the genome sequence version (file Vv12x_CRIBI.fa) on which the CRIBI2.1 annotation is based. To include the three VviCHS genes [43], structural gene annotation was optimised for VIT_214s0068g00920 and VIT_214s0068g00930.
Average day temperature values were retrieved from the weather station in the vineyard for the time from 1 November 2015 to 31 December 2016 (File S10). Air temperature values at 2 m above ground were used. Day lengths for Siebeldingen (Germany) during the sample period were calculated from publicly available times for sunrise and sunset for the geographic location of Siebeldingen.
HSP genes in CRIBI2.1 were identified based on the annotation text of homologs in A. thaliana by filtering for the strings "heat" and "shock" occurring together in the functional annotation text of the genes. The assignment between A. thaliana (Araport11 annotation, [57]) and V. vinifera genes is based on reciprocal best BLAST hits or at least best BLAST hits of the encoded peptide sequences as described before [54]. Lowly expressed genes were excluded from downstream analyses by applying a minimal CPM cutoff of 10 (applied to the CPM sum per gene over all samples). The Python package matplotlib v2.1.0 [58] was used for visualization of the data. The Spearman rank correlation coeffient of the HSP average transcript abundance per sample (median in Figure 3, arithmetic mean in File S5) and the average day temperature was calculated using the Python package SciPy [59]. As control, random sets of the same number of genes with the same lower CMP cutoff value were analysed in parallel.
Members of the transcription factor families WRKY [18], MADS-box [25], and MYB [21,41] were identified based on the published gene family analyses. The V. vinifera WRKY gene family has also been characterised by Guo et al. [60] which, unfortunately, resulted in conflicting gene designations. For consistency with Vannozzi et al. [42] we only used the VviWRKY gene designations of Wang et al. [18]. To allow the expression analysis of all previously described MADS-box genes, the CRIBI v2.1 annotation was extended with corresponding gene models using "VIT_230_" as prefix for the additional locus IDs. The Python packages matplotlib v2.1.0 [58] and seaborn v0.8.1 (https://github.com/mwaskom/seaborn) were used for visualization of RPKM values of selected TF genes in heatmaps (Files S6-S8).
Candidates for reference genes suitable for qRT-PCR experiments in the future were identified based on our comprehensive set of RNA-Seq samples. First, genes with a substantial expression level defined as the sum of all samples greater or equal to 500 [CPM] were selected. Second, these candidate set was filtered for a low variation defined as small standard deviation values across all samples normalized by the median of all values.

Funding:
We acknowledge the financial support of the German Research Foundation (DFG, http://www.dfg.de) to BW (WE1576/16-2) and to RT (TO152/5-2) in the context of SPP-1530. In addition, this paper is also based on work from COST Action CA 17111 INTEGRAPE, supported by COST (European Cooperation in Science and Technology).
the German Research Foundation and the Open Access Publication Fund of Bielefeld University for the article processing charge.

Conflicts of Interest:
The authors declare no conflict of interest.