In order to validate our protocol, we analysed H3K27me3 enrichment in the gene body of a DAM gene cluster on peach (Prunus persica L Batch) non-dormant buds (Fig. 2). H3K27me3 enrichment is known to be associated with a repressed transcriptional state. It has been shown that DORMANCY-ASSOCIATED MADS-BOX (DAM)-related genes display contrasting H3K27me3 profiles in dormant and non-dormant buds, with H3K27me3 signal only observed in non-dormant buds (de la Fuente et al. 2015). We observe a higher enrichment for H3K27me3 at the peach PpeDAM5 and PpeDAM6 loci compared with the PpeDAM3 gene in non-dormant buds, or compared to the gene PpeEF1, known to be enriched in H3K4me3 and depleted in H3K27me3 (Fig. 2; Fig. 3a). The replication of previously published results at the DAM genes confirms that our improved ChIP-seq protocol is working on tree buds.Fig. 2

H3K27me3 ChIP-seq profile at the DAM gene cluster in peach. IGV screenshot of ChIP-seq data for H3K27me3 and its corresponding input performed on non-dormant peach buds (Prunus persica) at PpeDAM3, PpeDAM5 and PpeDAM6 genes. Genes are represented by black lines, with arrows indicating gene directionality and rectangles representing exons

To demonstrate the versatility of this protocol, we performed ChIP-seq for two histone marks (H3K27me3 and H3K4me3) on buds of three tree species: peach, apple (Malus x domestica Borkh.) and sweet cherry (Prunus avium L.) (Fig. 3). We analysed the signal at the genes ELONGATION FACTOR 1 (EF1), known to be enriched in H3K4me3, and AGAMOUS (AG), known to be enriched in H3K27me3 (Saito et al. 2015). While H3K27me3 is associated with a repressive transcriptional state, H3K4me3 is on the contrary associated with transcriptional activation. We observed a strong H3K4me3 signal at EF1 and enrichment for H3K27me3 at AG locus for the three species (Fig. 3). This result confirms that this ChIP-seq protocol works on buds for several tree species. To test the adaptability of this protocol for other biological system, we also performed ChIP-qPCR in Arabidopsis thaliana and in Saccharomyces cerevisiae. We reproduced previously published results for the histone variant H2A.Z at HSP70 locus in Arabidopsis thaliana (Cortijo et al. 2017; Kumar and Wigge 2010), and for the binding of the TF Hsf1 at the SSA4 promoter in yeast (Erkina and Erkine 2006), Suppl. Fig. 1). This protocol is thus versatile and can be used on buds for several tree species as well as in other biological systems such as Arabidopsis thaliana and yeast.

To test if chromatin state and gene expression can be directly compared using this protocol, we carried out RNA-seq and ChIP-seq on the same starting material of sweet cherry floral buds. To start with, we analysed expression together with enrichment for H3K27me and H3K4me3 for PavEF1 and PavAG genes. We observe that PavEF1, marked by H3K4me3, is highly expressed, while PavAG, marked by H3K27me3, is very lowly expressed (Fig. 3c, Suppl. Fig. 4).

We then compared the abundance of H3K4me3 histone mark to the expression patterns of PavDAM5 and PavDAM6 during the dormancy period at three different dates for sweet cherry floral buds (October, December and January; Fig. 4). Firstly, we defined if the flower buds harvested are in endodormancy or ecodormancy at each time-point (Fig. 4a). The time of dormancy release, after which the buds are in ecodormancy, is defined when the percentage of bud break reaches 50% at BBCH stage 53 (Meier 2001). We observe that samples harvested in October and December are in endodormancy and the ones harvested in January are in ecodormancy (Fig. 4a). PavDAM5 and PavDAM6 are two key genes involved in sweet cherry dormancy corresponding to the peach PpeDAM5 and PpeDAM6, respectively. We find that PavDAM6 is highly expressed in October at the beginning of endodormancy and that its expression decreases in December and January (Fig. 4b), and that PavDAM5 is highly expressed in deep dormancy (December) and less expressed in ecodormancy (January, Fig. 4b). These results are in agreement with previous observations showing that DAM genes are upregulated in dormant peach buds (de la Fuente et al. 2015). We observe H3K4me3 enrichment at the beginning of these genes, at the level of the first exon (Fig. 4c, Suppl. Fig. 5) as expected for this chromatin mark. A low or no enrichment was found for the two controls (H3 and INPUT, Fig. 4c) meaning that the enrichment seen for H3K4me3 at these genes is relevant. We observe that H3K4me3 enrichment at PavDAM5 is higher in December compared to October and January (Fig. 4c). This is in agreement with the higher PavDAM5 expression in December. We observe that H3K4me3 enrichment at PavDAM6 is also higher in December compared to October and January, while its peak of expression is in October (Fig. 4c).Fig. 4

Expression and H3K4me3 profiles of PavDAM5 and PavDAM6 genes during dormancy in sweet cherry flower buds. a Evaluation of bud break percentage under forcing conditions (n = 5). b Transcriptional dynamics of PavDAM5 and PavDAM6 genes at three different dates (October 21, 2014, December 5, 2014 and January 27, 2015; n = 3). Expressions are represented in TPM (transcripts per kilobase million). c IGV screenshot of the first biological replicate of ChIP-seq data for H3K4me3 and H3 at three different dates (October 21, 2014, December 5, 2014 and January 27, 2015) and their corresponding inputs at PavDAM5 and PavDAM6 loci. Genes are represented by black rectangles, with arrows indicating gene directionality and bigger rectangles representing exons

To further investigate the link between gene expression and H3K4me3 enrichments during bud dormancy, we analysed genome-wide changes between time-points in expression level as well as H3K4me3 signal. For this, we employed a gene centric approach by measuring the strength of the H3K4me3 signal at each gene and identified genes showing significant changes between at least two of the three time-points (see “Material and methods” for more detail). We identified 671 genes that show significant changes in H3K4me3 between at least two of the sampling dates. We then performed a hierarchical clustering of these genes based on a Z score [(signal for a time-point - average over all time-points)/STD across all time-points], which normalises for differences in H3K4me3 signal between genes, thus allowing direct gene-by-gene comparisons of their patterns across the studied time-points. We observe a reduction in H3K4me3 signal over time, and in particular between endodormancy and ecodormancy for most genes (Fig. 5a, see Suppl. Figure 7b for examples of H3K4me3 signal at a few genes). We then compared changes in H3K4me3 and expression over time, for the genes in the main clusters depicted in Fig. 5a (Fig. 5b, c). Genes in the purple (233 genes) and blue (313 genes) clusters are associated with simultaneous reductions in H3K4me3 signal and expression over time (Fig. 5b, c). Genes in the green (53 genes), the gold (40 genes) and the orange (13 genes) clusters show increases in H3K4me3 signal and in expression over time (Fig. 5b, c). Genes in the red cluster (19 genes) are characterised by a transient reduction in H3K4me3 signal and expression in December (Fig. 5b, c).Fig. 5

Genome wide comparison of H3K4me3 and expression in sweet cherry flower buds in endodormancy and ecodormancy. a Hierarchical clustering of genes showing a change in H3K4me3 level between the time points, based on the Z score of the normalised H3K4me3 level around the TSS of each gene [Z score = (signal for a time point − average over all time points)/STD across all time points]. The result is represented as a heatmap where red indicates a high Z score and blue a low Z score. The genes were separated into six clusters, indicated by the side coloured bar. b Average signal in each time point for genes in clusters showing changes in H3K4me3 over time (identified in Fig. 5a) for expression (red) and H3K4me3 (blue). Mean normalised expression and H3K4me3 levels was used (level in one time-point/average across the entire time course for a given gene). The standard deviation is also shown as a transparent ribbon. c Expression level (TPM) and normalised H3K4me3 signal over time for a few individual genes in each of the 6 clusters identified on Fig. 5a. Each point represents one biological replicate. Normalised H3K4me3 signal corresponds to the number of H3K4me3 ChIP reads in a 2000 bp window around the TSS of each gene normalised by the number of H3 ChIP reads in the same window.

We explored signalling pathways that were represented in the different clusters (Fig. 5c, Suppl. Fig. 7b, Suppl. Table 2). Among the genes classed in the green, gold and orange clusters, that increase expression over time, we identified the GLUTATHION S-TRANSFERASE19 (PavGSTU19), LATE EMBRYOGENESIS ABUNDANT14 (PavLEA14) and GALACTINOL SYNTHASE 4 (PavGolS4) genes (Fig. 5c), potentially involved in the response to drought and oxidative stresses (Hara et al. 2010; Nishizawa et al. 2008; Singh et al. 2005), together with genes associated with growth and cellular activity such as PavSWEET2, GIBBERELLIN 20 OXIDASE 1 (PavGA20ox1) and AMINO ACID PERMEASE3 (PavAAP3). On the other hand, genes from clusters purple and blue that are downregulated during ecodormancy include XYLOGLUCAN ENDOTRANSGLUCOSYLASE/HYDROLASE6 (PavXTH6), DOF AFFECTING GERMINATION2 (PavDAG2), GA-STIMULATED IN ARABIDOPSIS6 (PavGASA6), TREHALOSE-PHOSPHATASE/SYNTHASE7 (PavTPS7), PavHVA22C and PHYTOCHROME-ASSOCIATED PROTEIN1 (PavPAP1) (Fig. 5c). Homologues of XTH6, DAG2, GASA6 and HVA22C were found to be activated under cold, drought and other stress-associated stimulus while TPS7 and PAP1 are more likely associated with metabolism and auxin pathway. We also found RESPIRATORY BURST OXIDASE HOMOLOGUE C (PavRBOHC) in the cluster red, transiently downregulated in December, of which the homologue in Arabidopsis thaliana is involved in reactive oxygen species production.

In this study, we described a combined ChIP/RNA-seq protocol for low abundance and complex plant tissues such as tree buds. This method allows a robust comparison of epigenetic regulation and gene expression as we use the same starting material. More notably, this protocol could permit to perform ChIP/RNA-seq for kinetic experiment with short intervals (every minute or less) and to collect samples in the field.

Several studies have led to the identification of molecular mechanisms involved in dormancy, including a cluster of DAM genes (Bielenberg et al. 2008). DAM-related genes are upregulated in dormant buds in peach (de la Fuente et al. 2015; Leida et al. 2012; Yamane et al. 2011), leafy spurge (Horvath et al. 2010), pear (Saito et al. 2015), apple (Mimida et al. 2015) and apricot (Sasaki et al. 2011). Conversely, DAM-related genes are downregulated in non-dormant buds. In particular, it was shown that H3K27me3 abundance is increased at DAM5 and DAM6 loci in non-dormant buds compared to dormant buds (de la Fuente et al. 2015). Using the proposed ChIP-seq protocol, we found similar results for H3K27me3 abundance in the DAM genes in peach in non-dormant buds (Fig. 2), thus validating our improved ChIP-seq method. In particular, we observe a higher enrichment for H3K27me3 at the DAM5 and DAM6 loci compared with DAM3 in non-dormant buds in peach, suggesting that not all DAM genes are regulated the same way during dormancy. Additionally, we demonstrated the correlation between the presence of histone marks and gene expression in sweet cherry (Prunus avium L.) for two control genes PavAG and PavEF1 known to be under control of H3K27me3 and H3K4me3, respectively (Fig. 3c). As the ChIP-seq and RNA-seq were performed on the same biological material, the level of gene expression and the presence/absence of particular histone marks can be directly compared with confidence. In the last two decades, ChIP has become the principal tool for investigating chromatin-related events at the molecular level such as transcriptional regulation. Our protocol will allow analysis of chromatin and expression dynamics in response to abiotic and biotic stresses, and this for trees in controlled conditions as well as growing in fields. Improvements to the ChIP-seq approach are still needed and will include an expansion of available ChIP-grade antibodies and a reduction of the hands-on time required for the entire procedure. A remaining challenge is to further decrease the amount of starting material without compromising the signal-to-noise ratio.

Previous studies highlighted the importance of DAM genes as key components of dormancy in perennials. We have shown that, in sweet cherry, PavDAM6 is highly expressed in October, at the beginning of dormancy, and then downregulated over time (Fig. 4). Conversely, we found that PavDAM5 is highly expressed at the end of endodormancy (December), and then downregulated during ecodormancy (Fig. 4). Both PavDAM5 and PavDAM6 genes were downregulated in ecodormancy, which suggests an important role in the maintenance of endodormancy as observed in Chinese cherry (Zhu et al. 2015), peach (Bielenberg et al. 2008; de la Fuente et al. 2015; Yamane et al. 2011) and Japanese apricot (Saito et al. 2015; Sasaki et al. 2011). However their timing of expression during endodormancy was different, suggesting that they have non-redundant roles and that their expression might be regulated by different factors during dormancy.

We conducted a ChIP-seq in sweet cherry for H3K4me3, which is associated with gene activation, in order to link the abundance of histone marks to expression patterns at three different dates along dormancy (October, December and January; Fig. 4). We found an H3K4me3 enrichment around the translation start site of both PavDAM6 and PavDAM5 during dormancy (Fig. 4), as reported in peach for PpeDAM6 (Leida et al. 2012) (Leida et al. 2012) and in leafy spurge for DAM1 (Horvath et al. 2010). Generally, H3K4me3 enrichment is more abundant at the PavDAM5 locus than at the PavDAM6 locus, and this is associated with the difference in their expression levels (Fig. 4). We find a positive relation between changes over time for H3K4me3 and expression level at PavDAM5 (Fig. 4). Similarly, we find a general trend that genes with an increase over time in H3K4me3 level also show an increase in expression, while genes with a reduction in H3K4me3 over time also show a reduction in expression (Fig. 5). Despite these general correlations, changes in H3K4me3 signal that occur at a significant magnitude have been observed for only 671 genes. These results suggest that while chromatin marks might be involved in transcriptional regulation of some genes in sweet cherry buds during dormancy, other genes might exhibit changes in expression that are not associated with chromatin changes. In other genes, like PavDAM6, any relation between chromatin marks and expression seems to be more complex and possibly integrates other chromatin marks than H3K4me3. Performing a similar analysis for other chromatin marks, such as H3K27me3 that is associated with transcriptional repression, might reveal complex links between chromatin and expression dynamics during bud dormancy in sweet cherry.

Our analysis has been done with only 3 time-points spanning a period of 4 months. A finer time resolution would be needed to better identify relations between chromatin expression dynamics during dormancy, and in particular to define if chromatin changes happen before or after expression changes. Together with our results, the observation from Leida et al. (2012) and De la Fuente et al. (2015) of the presence of repressive histone mark such as H3K37me3 in peach DAM loci during dormancy support the hypothesis of a balance of histone mark enrichment controlling the dormancy process in sweet cherry and probably more largely in perennials.

In addition to DAM genes, we identified genes that showed differential H3K4me3 enrichment between the different bud dormancy stages. We found genes involved in the response to drought, cold and oxidative stresses that were either highly expressed at the beginning of endodormancy, such as PavGSTU19 and PavLEA14, or expressed during ecodormancy (PavXTH6, PavDAG2, PavGASA6 and PavHVA22C). This is consistent with previous studies showing the key role of redox regulation (Ophir et al. 2009; Considine and Foyer 2014; Beauvieux et al. 2018) and stress-associated stimulus (Maurya and Bhalerao 2017) in bud dormancy. Interestingly, we highlighted PavGASA6 and PavGA20ox1, two genes associated with the gibberellin pathway. In particular, PavGA20ox1, that encodes an enzyme required for the biosynthesis of active GA (Plackett et al. 2012), is markedly upregulated after endodormancy release, associated with H3K4me3 enrichment at the locus, therefore, suggesting that epigenetic and transcriptomic states facilitate an increase in active GA levels during ecodormancy, as previously shown in hybrid aspen (Rinne et al. 2011) and grapevine (Zheng et al. 2018).

Finally, our results highlighted the CHROMATIN REMODELING5 (PavCHR5) gene, upregulated during ecodormancy. In Arabidopsis, CHR5 is required to reduce nucleosome occupancy near the transcriptional start site of key seed maturation genes (Shen et al. 2015). We can therefore hypothesise that chromatin remodelling, other than post-transcriptional histone modification, occurs during dormancy progression. Our results confirm that the analysis of differential H3K4me3 states between endodormant and ecodormant buds might lead to the identification of key signalling pathways involved in the control of dormancy progression.

Sweet cherry trees (cultivar ‘Burlat’), apple trees (cultivar ‘Choupette’) and peach trees (unknown cultivar) were grown in an orchard located at the Fruit Experimental Unit of INRA in Toulenne (France, 44° 34′ N 0° 16′ W) under standard agricultural practices. Sweet cherry flower buds used for the RNA-seq and ChIP-seq experiment were collected on October 21, 2014 and December 5, 2014 for dormant buds and January 27, 2015 for non-dormant buds. Apple buds were collected on January 25, 2016 and peach buds on February 5, 2016. Buds were harvested from the same branches, flash frozen in liquid nitrogen and stored at − 80 °C prior to performing ChIP-seq and RNA-seq.

Three branches bearing floral buds were randomly chosen from the sweet cherry cultivar ‘Burlat’ trees at different dates. Branches were incubated in water pots placed in forcing conditions in a growth chamber (25 °C, 16 h light/8 h dark, 60–70% humidity). The water was replaced every 3–4 days. After 10 days under forcing conditions, the total number of flower buds that reached the BBCH stage 53 (Meier 2001) was recorded. We estimate that endodormancy is released when the percentage of buds at BBCH stage 53 is above 50% after 10 days under forcing conditions.

Timing 2–3 h. Work on ice, except when specified otherwise. Add 25 ml of ice-cold Extraction buffer 1 [0.4 M sucrose, 10 mM HEPES pH 7.5, 10 mM MgCl₂, 5 mM β-mercaptoethanol, 1 mM PMSF, 1% PVP-40 (polyvinylpyrrolidone), 1 tablet of complete protease inhibitor EDTA free for 50 ml of buffer from Sigma cat no. 11836170001] to the powder. For each buffer, the protease inhibitor (PMSF), the tablet of complete protease inhibitor and β-mercaptoethanol should be added directly before using the buffer. PVP-40 is a chelator used to remove phenolic derivatives as well as polysaccharides and improve the quality of the chromatin extraction. It is commonly used in RNA and chromatin extractions in buds sampled from fruit trees (Gambino et al. 2008; Leida et al. 2012; Ionescu et al. 2017). When extracting chromatin in other biological systems such as Arabidopsis thaliana or Saccharomyces cerevisiae, the PVP-40 in Extraction buffer 1 may optionally be removed.

Immediately add 675 μl of 37% formaldehyde solution (1% final concentration) and invert the tube several times to resuspend the powder. Cross-link the samples by incubating at room temperature for 10 min and then quench the formaldehyde by adding 1.926 ml 2 M of fresh glycine solution (0.15 M final concentration). Invert the tube several times and incubate at room temperature for 5 min. Filter the homogenate through Miracloth (Millipore cat no. 475855) in a funnel and collect in a clean 50-ml Falcon tube placed on ice. Repeat the filtration step once more. Centrifuge the filtrate at 3200×g for 20 min at 4 °C. Discard the supernatant by inverting the tube, being careful not to disturb the pellet. Gently resuspend the pellet in 1 ml of Extraction buffer 2 [0.24 M sucrose, 10 mM HEPES pH 7.5, 10 mM MgCl₂, 1% Triton X-100, 5 mM β-mercaptoethanol, 0.1 mM PMSF, 1 tablet protease inhibitor EDTA free for 50 ml of solution], without creating bubbles, and transfer the solution to a clean 1.5-ml tube. Centrifuge at 13,500×g for 10 min at 4 °C. Carefully remove the supernatant by pipetting. If the pellet is still green, repeat the resuspension in 1 ml of Extraction buffer 2, centrifuge at 13,500×g for 10 min at 4 °C, and remove the supernatant. In a new 1.5-ml tube, add 300 μl of Extraction buffer 3 (1.7 M sucrose, 10 mM HEPES pH 7.5, 2 mM MgCl₂, 0.15% Triton X-100, 5 mM β- mercaptoethanol, 0.1 mM PMSF, 1 mini-tablet protease inhibitor EDTA free for 50 ml of solution). Slowly resuspend the pellet in 300 μl of Extraction buffer 3 to prevent the formation of bubbles. Take the resuspended pellet and carefully layer it on top of the 300 μl Extraction buffer 3. Centrifuge at 21,200×g for 1 h at 4 °C. During this process, nuclei are pelleted through a sucrose cushion to remove cellular contaminants.

From this step, chromatin fragmentation can be performed in two different ways: (A) sonication, to shear the chromatin into 100–500 bp fragments, or (B) MNase (Micrococcal nuclease) digestion, to enrich for mono-nucleosomes (~ 150–200 bp). Results shown in this paper come from ChIP-seq performed on sonicated chromatin.

We recommend for the RNA extraction the use of RNeasy® Plant Mini kit from Qiagen (cat no. 74904) for less than 50 samples with the following minor modifications:

Start from 50 to 70 mg of buds powder with scales. Only remove the tubes from the liquid nitrogen when the RNA Extraction buffer is prepared.

Alternatively, RNA can be extracted using the MagMAX™-96 Total RNA Isolation Kit from Thermo Fisher (cat no. AM1830) for more than 50 samples following manufacturer’s instructions. Store RNA at − 80 °C.

We recommend carrying out RNA-seq library preparation using the Truseq Stranded mRNA Library Prep Kit from Illumina (96 samples, 96 indexes, Illumina cat no. RS-122-2103). Check the quality of libraries using 4200 TapeStation or Bioanalyzer instruments (Agilent) following manufacturer’s instructions. See Suppl. Fig. 3 for profiles of libraries with and without adapter contaminations. If the libraries are contaminated with adapters, continue with a size selection step with SPRI beads to remove them, otherwise, proceed directly with ‘the quantification and pool of libraries’ section.