Dynamics of H3K4me3 Chromatin Marks Take the Lead over H3K27me3 for Gene Regulation during Flower Morphogenesis in Arabidopsis thaliana

Plant life-long organogenesis involves sequential, time and tissue specific expression of developmental genes. This requires activities of Polycomb Group (PcG) and trithorax Group complexes, respectively responsible for repressive Histone 3 trimethylation at lysine 27 (H3K27me3) and activation-related H3K4me3. However, the genome-wide dynamics in histone modifications that occur during developmental processes have remained elusive. Here, we report the distributions of H3K27me3 and H3K4me3 along with transcriptional changes, in a developmental series including Arabidopsis leaf and three stages of flower development. We found that chromatin mark levels are highly dynamic over the time series on nearly half of all Arabidopsis genes. Moreover, during early flower morphogenesis, changes in H3K4me3 prime over changes in H3K27me3 and quantitatively correlate with transcription changes, while H3K27me3 changes occur after prolonged expression changes. Notably, early activation of PcG target genes is dominated by increases in H3K4me3 while H3K27me3 remains present at the locus. Our results reveal H3K4me3 as greater predictor over H3K27me3 for transcription dynamics, unveil unexpected chromatin mechanisms at gene activation and underline the relevance of tissue-specific temporal epigenomics.


Introduction
Morphogenesis in plants relies on continuous and iterative production of new organs with specific functions. In annual plants like Arabidopsis thaliana, upon transition to the reproductive 2 of 32 phase, flower primordia are produced at the flanks of the inflorescence meristem. In these primordia, stem cell activity is re-established and the newly formed flower meristems produce founder cells for the different floral organs [Alvarez-Buylla, 2010. These cell-fate specifications coincide with changes in the expression pattern of developmental regulators . The LEAFY (LFY), APETALA 1 (AP1) and CAULIFLOWER (CAL) transcription factors (TF) specify floral meristem identity via activation of downstream regulators, among which the MADS domain TF encoding genes. Spatially restricted expression of these genes then induces the differentiation of stem cells into specific cell types for formation of the different floral organs, as described in the ABC(D)E model [Chahtane, 2017. Several reports on epigenetic mutants indicate that this regulation occurs through changes in the chromatin state of the corresponding loci, such as post-translational modifications of histones [Engelhorn, 2014 #169]. Histone 3 trimethylation at lysine 27 (H3K27me3) is a repressive modification targeting developmental regulators in all higher eukaryotes and is catalysed by Polycomb Group (PcG) proteins [Pirrotta, 1998Schuettengruber, 2007 #518; Morey, 2010 #395]. Genome-wide studies of H3K27me3 distribution in whole seedlings revealed that this mark targets more than one fourth of all Arabidopsis genes [Zhang, 2007Oh, 2008 #2408;Roudier, 2011 #2385], comprising tissue-or stage-specific regulators of the vegetative to reproductive transition, floral meristem identity genes and all above-mentioned floral homeotic genes. Activation of PcGrepressed genes is conducted by the conserved trithorax Group (trxG) proteins, involved in three main activation processes: (i) removal of repressive H3K27me3 marks, (ii) deposition of activationassociated H3K4me3 marks and (iii) chromatin remodelling to allow access of the transcriptional machinery [Schuettengruber, 2011 #519; #169; Kingston, 2014 #303].
Thus, PcG and trxG antagonistic factors ensure correct spatial and temporal expression of developmental regulators, such as those involved in flowering-time and flower development [Alvarez-Venegas, 2010 #2412;Kohler, 2010 #2450;Muller, 2011 #401;de la Paz Sanchez, 2015 #143; Pu, 2015 #458]. As a consequence, H3K27me3 and H3K4me3 chromatin marks should be highly dynamic throughout developmental stages and cell-types in plants. However, most genome-wide studies were so far performed for one developmental stage and often in a mixture of tissues during vegetative development, yielding only a static view of chromatin states. Some comparative studies used endosperm-specific nuclei [Weinhofer, 2010, meristem-enriched tissue [Lafos, 2011, roots [Roudier, 2011, or mature versus senescing leaves [Brusslan, 2012Brusslan, 2015]. All showed differential H3K27me3 marking between tissues, at thousands of genes. Additionally, several hundreds of genes differ in H3K27me3 intensity between undifferentiated and differentiated tissue, without loosing the mark completely [Lafos, 2011He, 2012 #632]. Several genome-wide analyses revealed an anti-correlation between changes in gene expression and changes in H3K27me3, and oppositely, a positive correlation for the H3K4me3 mark [Lafos, 2011Brusslan, 2012]. Studies at single gene level provided evidence that correlative changes in these histone marks are rather a cause than a consequence of the transcriptional state [Adrian, 2010Buzas, 2011.
Nevertheless, many questions remain open such as (i) the dynamics of histone mark between consecutive developmental stages; (ii) the correlation extent between changes in H3K27me3, H3K4me3 and gene expression; (iii) and the chronology of changes in H3K27me3, H3K4me3 and transcription. To obtain answers to these questions, we used the Arabidopsis "cauliflower" inducible system, which allows collecting meristem tissue from synchronized floral stages, in sufficient amounts for genome-wide analyses of chromatin and expression [Wellmer, 2006 #2396;Smaczniak, 2012 #1999; Pajoro, 2014. We resolved the dynamics of the H3K4me3 and H3K27me3 marks and their correlation with transcriptional changes in vegetative, inflorescence meristem, early developing flower and developed inflorescence tissues. This allowed establishing a comprehensive view of antagonistic histone modification dynamics during flower development in relation with gene activation or repression. We found that both chromatin marks are highly dynamic over development, with thousands of genes changing in their marking. H3K4me3 dynamics correlate quantitatively with expression changes, while pronounced H3K27me3 changes mostly occur after transcriptional changes, with the exception of genes that are activated and required very early during flower morphogenesis, like SEPALLATA1-3 (SEP1-3), LFY and an organ boundary specifying gene (CUP-SHAPED COTYLEDON 1 (CUC1)). However, also for these genes, early changes in H3K27me3 are mild and less pronounced than H3K4me3 changes. Thus we find that within the first two days of flower initiation, gene activation is dominated by deposition of H3K4me3, while strong reduction of H3K27me3 occurs at later stages. Together, our integrative analysis of histone marks and transcriptome provides a novel view on the mechanisms of gene regulation at the chromatin level.

RNA-seq reveals specific changes in gene expression during the developmental time series
To assess the correlation between changes in expression patterns and dynamics of histone marks during early flower development, we employed the ap1cal 35S::AP1-GR floral induction system introduced in Wellmer et al. (2006) [Wellmer, 2006 (Fig. 1A; Supplementary Fig. 1). This system allows harvesting synchronized flower buds along a developmental time course and insures (i) tissue homogeneity along the developmental time course as well as (ii) sampling in sufficient amounts for ChIP Smaczniak, 2012Pajoro, 2014].
The early time points we chose cover early flower development during which flower buds emerge from the inflorescence meristem and further initiate the flower patterning program. During this period, cells undergo fate transitions, being first re-allocated a flower stem cell fate, followed by differentiation into founder cells for specific floral organs. A suitable time point to study the early events of flower formation is one at which (i) changes in expression are already detectable for important regulators of flower development and (ii) floral organ identity genes did not reach their full expression level, to enable observation of transition states from repression to activation. In our experimental setup, two days after induction (t2) stands as the appropriate time point for investigating activation of floral organ identity genes and at which FM identity is still present (Fig. 1B, [Smaczniak, 2012). Moreover we previously showed using ChIP-PCR, that reduction of repressive H3K27me3 marks at the SEP3 locus is detectable at t2 [Smaczniak, 2012 and thus decided to perform ChIP-seq for H3K27me3 and H3K4me3 and RNA-seq on a tissue series consisting of wild-type leaf (L) and inflorescence (I), and ap1cal 35S::AP1-GR t0 and t2 samples.
Clustering analysis of expression data for selected floral regulators (flowering time, flower development) showed the presence of five major patterns of expression during the series (Fig. 1C-E). The assignment of the genes to the clusters revealed that the observed expression patterns are in accordance with known functions of floral regulator genes. It also confirms the t2 time-point as an intermediate stage during the activation of floral organ and meristem identity genes (Supplementary Text 1).
Genome-wide analysis of differentially expressed genes (DEG) showed that most genes changed from L to I, emphasizing the different nature of the two tissues. 251 genes are up-regulated from t0 to t2, among them many developmental regulators like SEP2, SEP3, FLOWERING LOCUS T (FT), UNUSUAL FLORAL ORGANS (UFO), AP3, AG, PI, AGAMOUS-LIKE 15 (AGL15), LFY and CUC1, while 448 are down-regulated (Table 1, Supplementary Table 1). K-means clustering and Gene Ontology (GO) analysis revealed that the gene expression is highly dynamic throughout the series and correlates with expected functionality of the tissues on a genome-wide scale (Supplementary Text 2. and 3.). Thus RNA-seq analysis demonstrates transcriptome dynamics throughout the tested tissue series and validates the suitability of the selected tissue types and time-points to study the influence of chromatin marks on genome expression. analysis (L, I, t0 and t2) (Supplementary Text 4.). The average expression values of H3K27me3 target genes are lower compared to all Arabidopsis genes, while H3K4me3 target genes display a higher average expression. This is true for all analysed tissue types and in accordance with H3K27me3 and H3K4me3 being associated with repressed and active expression states, respectively (Supplementary  Table 3).
Mark abundance distributions over loci for both marks match with previous descriptions [Zhang, 2007Oh, 2008Zhang, 2009Roudier, 2011: H3K27me3 is spread over the whole transcribed region, into the promoter and the 3' UTR, while H3K4me3 is more concentrated to the TSS, covering the proximal promoter and the 5' region of the gene ( Fig. 2A, Supplementary Fig. 7). These overall distributions are observed for all tissues including the meristematic stages t0 and t2, indicating that plant inflorescence meristems share similar H3K27me3 and H3K4me3 marking with differentiated cells.
When sorted by H3K4me3 abundance and assessed for transcription in the four tissues ( Fig.2A), three main groups of genes emerge: First, highly expressed H3K4me3-marked genes with a strong correlation between expression levels and both intensity and spread of the H3K4me3 signal. Second, lowly expressed genes that carry low H3K4me3 levels and high H3K27me3 levels. Many genes in this group display higher expression in the differentiated L and I samples than in the meristematic tissue. This is in line with the previous observation that expression patterns of H3K27me3 target genes are more tissue specific than other genes [Zhang, 2007. Third, genes with low marking in both marks. Interestingly, this is the group of genes with the lowest expression, probably representing genes carrying other repressive marks like DNA methylation or H3K9me2.

H3K27me3 and H3K4me3 levels are highly dynamic over the time series, in accordance with tissue functionality
Determination of differentially marked genes (DMGs) for the L-t0, L-I, t0-t2, t0-I and t2-I comparisons revealed that quantitative changes in H3K27me3 and H3K4me3 marks affect thousands of genes (Fig. 2C, Supplementary Table 4). In the case of very different tissue types (L-t0 or L-I), changes affect more than 60 % (H3K27me3) or 45 % (H3K4me3) of the corresponding target genes, respectively. Noticeably, the vast majority of genes showing an elevation of any of the two marks from L to t0 shows a reduction of the same mark from t0 or t2 to I, and vice versa (Supplementary Table 4), indicating that these changes are meristem specific.
For several hundred genes we detected changes in both antagonistic marks in opposite directions (H3K27me3 elevated and H3K4me3 reduced or vice-versa), while only few genes (none for the t0-t2 comparison) are changed in the same direction for both marks (Supplementary Table 4).
The examples of histone mark distribution over the SEP1 and SEP2 loci for t0, t2 and I (Fig. 2B) indicate that the changes, especially affecting H3K27me3, are subtle but span a wide range of the marked region.
Our quantitative analysis allowed to specifically identify genes with marks changing in the two first days of flower initiation, from t0 to t2 (early flower morphogenesis). Among these, 13 DMGs change in both marks towards activation (H3K27me3 reduced, H3K4me3 elevated) and 15 towards repression (H3K27me3 elevated, H3K4me3 reduced). These changes relate to coordinated activation or repression of key genes for flower development (Supplementary text 5).
A GO analysis further revealed that the functions of the H3K27me3 and H3K4me3 DMGs in the series correlate with the functionality of the tissues at the analysed time-points (Supplementary Text 6). To evaluate if changes in repressive or active histone marks correlate with changes in expression of the corresponding loci, we determined the fractions of differentially up-and down-regulated genes that display changes in histone mark abundance (Fig. 3A).

DMGs largely overlap with DEGs and histone mark changes are stronger among DEGs
For most comparisons, including all for H3K4me3, expression changes occur more frequently in the expected direction (up-regulation for elevated H3K4me3 and reduced H3K27me3 and vice-versa). This is especially true when comparing t0 to t2. In all cases, the number of DEG in the expected direction is significantly higher among DMGs than for all targets (based on a hypergeometric distribution with probability cut-off of 0.01). With the exception of the L to I comparison, DMGs in both marks display a higher percentage of DEGs compared to DMGs in one mark. These findings are significant for all cases in the expected direction with the exception of t0-t2 in the activation direction (based on a hyper-geometric distribution with probability cut-off of 0.01), indicating that combined changes in both marks are more prone to coincide with changes in transcriptional activity.
Since transcription is initiated by the assembly of the pre-initiation complex (PIC) near the transcription start site (TSS) of a gene, histone marks might primarily change at this region to allow activation or repression of transcription. We thus analysed the marks average distribution along the locus for genes displaying changes during induction of flower development (DMG from t0 to t2 or I) (Figs. 3B and S7). The t0-t2 and t0-I differences are clearly visible among the DMGs and occur rather uniformly over the covered regions than locally, except for a stronger change at the peak close to the TSS for H3K4me3 in the t0-t2 activation comparison. Changes in H3K4me3 are stronger when DEG subsets are considered, while for H3K27me3 this is only true for the t0-I activated genes and the t0-t2 repressed genes.
Notably, DMGs in one mark (regardless if H3K27me3 or H3K4me3) show a higher abundance of the other mark compared to all target genes, and especially DMGs changing in expression. This indicates that these DMGs genes are in a transition state. For H3K27me3 DMGs genes, a change in H3K4me3 abundance is visible in the opposite direction for all comparisons, and this change is stronger for genes changing in expression. This is in accordance with our observation that genes changing in both marks more frequently display changes in expression than genes changing in one mark.
Taken together, our analysis of DMGs and DEGs established a relation between H3K27me3/H3K4me3 dynamics and changes in transcriptional states during flower development. Furthermore, it shows that during early stages of flower morphogenesis, H3K4me3 dynamics affects more genes and is in general more strongly correlated with transcriptional changes than H3K27me3 dynamics.

Expression changes during early flower morphogenesis quantitatively correlate with changes in H3K4me3 while changes in H3K27me3 mostly occur after prolonged expression changes
Next, we aimed to test for a quantitative correlation between histone mark changes and expression changes. For this, we mainly focused our analysis on the genes changing during early flower morphogenesis (t0 to t2) as these constitute the most homogenous tissue samples, and because genes changing in this time frame are highly dynamic and comprise many important developmental regulators. A clear correlation was observed between H3K4me3 changes and expression changes: mark changes are of higher intensity, of wider spread, and greater towards the TSS for genes with stronger expression changes (Fig. 4A).
As H3K4me3 changes are on average slightly stronger at the TSS, we further considered all DEGs that change early but remain DEGs in the same direction for the t0 to I comparison, thus allowing to follow activation at early and later stages. For those genes, we analysed the average distribution of changes from t0 to t2 and I and t2 to I over the locus (Supplementary Fig. 8). A clear shift from the TSS towards further downstream region is visible from early to late changes for up-regulated genes, indicating that H3K4me3 changes occur mainly at the TSS during early activation but then spread more into the locus. For down-regulation, the opposite effect is observed but less pronounced.
For H3K27me3, no correlation between expression changes and mark changes is observed for all DEGs from t0 to t2. This could be due to only some of the DEGs being H3K27me3 targets, while H3K4me3 marks target a broader range of genes. Therefore, we considered only those among the early DEGs that are H3K27me3 targets. We again could observe only a slight preference for changes in the expected direction and no correlation with the intensity in expression change (Fig. 4B). Only when considering H3K27me3 target genes that had been up-or down-regulated for a longer time (i.e. DEGs which expression changed early but remained in the same direction for the t0 to I comparison), we could observe a clear prevalence of changes in the expected direction (Fig. 4C). Especially for activated genes of this category, the H3K27me3 signal changes over the whole locus including the promoter region for most genes and the change is more pronounced for genes with higher expression changes. This indicates that changes in H3K27me3 rather occur after prolonged changes in expression and is in accordance with a low number of H3K27me3 DMGs during early flower morphogenesis. However, there are some changes in H3K27me3 already observed for the t0 to t2 comparison (e.g. the MADS developmental regulators SEP1-3). We hypothesised that these genes might be very early activated/repressed genes. Indeed, we found that most DEGs with high average H3K27me3 change in the expected direction displayed already an increase/decrease in expression during the first day of flower development as determined by Wellmer et al. (2006) [Wellmer, 2006] (Supplementary Fig. 9).
Taken together, during flower morphogenesis, we observe that H3K4me3 dynamics are strongly predictive of expression changes, while H3K27me3 changes are rather predicted by already on-going transcriptional changes (Supplementary text 7).

Early activation of PcG target genes is dominated by H3K4me3 changes
Many regulators of flower development are known to be H3K27me3 targets and regulated by antagonising trxG and PcG factors [Zhang, 2007Gan, 2013 #198]. Therefore, we were interested to see whether changes in genes that will eventually significantly loose H3K27me3 and gain H3K4me3 are also dominated by H3K4me3 dynamics during early flower morphogenesis. We analysed DMGs that are targets of H3K27me3 at t0 in more details (Fig. 4D). Of these 8273 H3K27me3 target genes, 2776 become differentially methylated in one of the marks or both during flower morphogenesis i.e. from t0 to t2 and/or from t0 to I. The majority of them show only reduction in H3K27me3 (2221 genes), while 210 show only elevation in H3K4me3. The 345 genes that show a change in both marks are suitable to address the question of which mark is changed first. We thus further considered the subpopulation of these 345 genes, which show changes in marks during early activation, i.e. from t0 to t2 (48 genes). Among these, the proportions of changes in marks are inverted compared to all genes changing during flower morphogenesis: 6 genes change only in H3K27me3, 13 change in both marks (corresponding to the 13 genes discussed above and mentioned in Table 2) and the majority (29 genes) changes only in H3K4me3. This again indicates that during initiation of flower development (t0-t2) genes carrying H3K27me3 can display an increase in H3K4me3 without significantly changing in their H3K27me3 level within the observed tissue. Such effects have been shown also for stress-induced transcription at single gene level [Liu, 2014]. On the other hand, only few genes display a reduction in H3K27me3 without significant changes in H3K4me3.

Relative dynamics of histone marks, MADS TF binding and chromatin opening during early stages of gene activation
Since flower morphogenesis is mainly regulated by activation of ABC(D)E genes, and especially MADS domain TFs, we were interested in the correlation between MADS binding to chromatin [Pajoro, 2014] and histone mark changes (this study). We thus analysed the fractions of t0-t2 DMGs, which are targets of the MADS TFs SEP3 and AP1 in the floral induction system as determined by Pajoro et al. (2014) [Pajoro, 2014 at two, four and eight days after induction (Fig. 5, Supplementary Fig. 11A). At t2, slightly more than 20 % of genes with reduced H3K27me3 or elevated H3K4me3 are targeted by the MADS TFs, while almost 40 % of the genes changing in both marks are. At t4, the percentage of genes being MADS targets raises to around 60 % for those changing in a single mark and to more than 90 % for those changing in both marks (12 of 13 genes are MADS targets at t4). Towards t8 only minute changes occur in this distribution. These percentages of DMGs bound by the MADS TFs are higher than that observed for all H3K27me3 or H3K4me3 target genes: only 6.7 % of H3K27me3 (8.9 % for H3K4me3) target genes are MADS targets at t2, increasing to 23.4 % (33 % H3K4me) at t4 and 26.9 % (35.7 % H3K4me3) at t8.
As MADS TFs were proposed to assist in opening chromatin, we analysed if the genes changing in histone marks from t0 to t2 would also change in their chromatin conformation as determined in Pajoro et al. (2014) [Pajoro, 2014] by assessing DNase I hypersensitivity (DHS, Fig. 5, Supplementary Fig. 11B). Interestingly, about 40% or 60% of those DMGs are already located in DNase I hypersensitive (accessible) regions at t0. This is quite unchanged in t2 and the percentage of genes in open regions only increases slightly towards t4 and more towards t8. MADS binding seems to mainly affect other, still closed genes at t2, while a larger overlap between MADS binding and DHS is observed at t4 and even more at t8 (Fig. 5).
This dataset comparison indicates that the majority of genes changing in H3K27me3 and H3K4me3 histone marks during early steps of flower morphogenesis are targets of the floral organ identity regulators SEP3 and AP1, and that binding of these MADS TFs occurs at the same time as histone mark changes or slightly later. Furthermore, it would appear as if opening of previously closed chromatin occurs rather after histone mark changes and MADS binding.    The upper Venn-diagram shows all genes changing during flower morphogenesis (that is genes that change either from t0 to t2, t0 to I or both) with fraction of those genes changes in expression (dark red numbers in ()), in H3K27me3 (blue fraction), in H3K4me3 (red fraction) or both (orange).
The lower Venn-diagram shows all genes that change in both marks (345 genes) which already show activation related changes in marks and expression from t0 to t2.

Chromatin marks are highly dynamic and tissue fate-associated during flower morphogenesis
Several mutants in PcG and trxG components display defects in initiation of flowering and in floral architecture [Chahtane, 2017. Moreover, almost all floral organ identity genes and many flowering time regulators are targeted by H3K27me3, and detailed studies at single gene level (e.g. FT and FLOWERING LOCUS C) have shown that activation of these genes is accompanied by histone mark dynamics. These facts led to the assumption that PcG-associated H3K27me3 and trxG-associated H3K4me3 dynamics operate at the onset of flowering and flower development, affecting a wide range of target genes. However, the extent of the chromatin marks remodelling remained unknown. Through the floral induction system we were able to generate a genome-wide map of dynamics in histone mark abundance and of transcriptional activity, during early steps of flower development. Our study demonstrates that histone mark dynamics between vegetative and reproductive structure, and during flower initiation and morphogenesis truly exists and affects a total of 14826 genes (changed at least in one mark in one comparison, Tab S4) in Arabidopsis thaliana.
Several thousands of genes changed in histone marks from leaf towards inflorescence meristems and then changed in the opposite direction from inflorescence meristems towards fully expanded inflorescences. This indicates that histone marks antagonistically change in the different tissue types according to the underlying meristematic activities. These reverted patterns of histone mark changes throughout the time series thus likely reflect the chromatin features of differentiated versus meristemlike tissues and highlight the continuous and iterative mode of plant development.

Integrative analysis of chromatin and transcription datasets brings new insights on the mechanistic of gene activation, revealing unexpected order of events
Combination of the transcriptomics and chromatin dynamics datasets revealed that changes in histone marks are widely correlated with changes in expression in the anticipated direction. Especially, changes in both marks are more prone to coincide with expression changes.
The temporal resolution brought by our dataset revealed that activation of a locus does rather not follow the canonically assumed order of events, with repressive marks being removed and active marks being added during transcriptional activation. Some studies had challenged this view in the past, at single gene level in plants. They report that transcription can take place while repressive marks are still present, and reduction in these marks was proposed to rather be a consequence of expression of the target gene [Schubert, 2006;Adrian, 2010 #2508;Buzas, 2011 #2509]. Our genome-wide study shows that H3K4me3 is both quantitatively and temporally closely correlated with transcription changes, while changes in H3K27me3 most likely require prior transcriptional changes to occur (Fig. 6). This is resonant with the recent report in mouse that re-establishment of H3K4me3 occurs more rapidly than that of H3K27me3 during activation of the zygote genome [Liu, 2016. In mammalian cells, deposition of H3K27me3 by PcG proteins is hindered by the presence of H3K4me3 [Schmitges, 2011. If this would be the case in plants as well, deposition or removal of H3K4me3 would hinder or allow PcG proteins to access the respective histone and thus lead to decrease or increase in H3K27me3. In the case of removal of H3K27me3, that would require cell division to occur and thus could explain the temporal delay in H3K27me3 removal after gene activation. Moreover, we found that during flower formation, transcriptional activation events of PcG target genes are predominantly accompanied by changes in H3K4me3 while H3K27me3 marks remain at the loci, and only decline later during flower morphogenesis. This indicates that initiation of floral morphogenesis involves more activation than de-repression at the chromatin. The fact that decrease in H3K27me3 may be a consequence of gene transcription, could explain our observation of histone mark abundance changes over the complete gene body (even for a short t0-t2 transition) and the fact that genes changing in expression show a stronger change in histone mark abundance. On the other hand, early H3K4me3 changes are stronger at the TSS where the pre-initiation complex (PIC) is assembled during initiation of transcription.
Interestingly, genes changing in one histone mark display higher abundance of the respective other mark, even preceding the changes in one mark. This could indicate that these genes are kept in an intermediate state to allow fast regulation during development and raises the question whether active and repressive marks co-exist at a locus in a same cell on the same chromosome, at the time of flower bud initiation. Such so-called bivalent loci were described in the specific case of animal embryonic stem cells [Bernstein, 2006 #741] and more recently reported in plant studies using Arabidopsis seedlings [Luo, 2013;Sequeira-Mendes, 2014 #1955]. Thus far, we could not isolate bivalent states in inflorescence meristems and early developing flowers from the ap1cal 35S::AP1-GR tissue, by re-ChIP experiments. Whether they do not or very sporadically exist in a restricted amount of cells of the flower meristem, deserves to be further questioned in studies using technologies specifically designed for this purpose.
Finally, integrative analysis of histone marks, transcriptomics, profiling of MADS-domain proteins and chromatin accessibility, indicates that the majority of genes changing in H3K27me3 and H3K4me3 at early steps of flower morphogenesis are targets of the floral organ identity regulators SEP3 and AP1 and display local chromatin conformation change. Moreover, binding of the MADSdomain proteins occurs at the same time as histone mark changes or slightly later. This observation brings H3K4me3 and H3K27me3 dynamics at the forefront of molecular switches for gene expression, together with TFs that have been proposed to act as pioneer factors during flower development [Pajoro, 2014].
After 2-3 days of stratification at 4°C, plants were germinated on Murashige and Skoog medium supplemented with 0.3 % sucrose in growth cabinets (Percival) in a controlled environment in long days (LD) of 16 hours light at 18°C and 8 hours darkness at 16°C. At about ten days after germination, plants were transferred to a 1:1 mixture of soil/vermiculite.

Dexamethasone treatment and harvesting of plant material
The ap1cal 35S::AP1-GR system employs double mutants of the floral meristem identity genes Replicates for biological sampling were grown independently (at different times) and harvested from the same groups of plants for ChIP-seq (two replicates) and RNA-seq (three replicates).

Scanning Electron Microscopy
For scanning electron microscopy [Bowman, 1989, untreated inflorescence-like meristems of ~4 week old ap1cal or ap1cal 35S:AP1-GR plants were harvested, as were 5 d-post dexamethasonetreated ap1cal 35S:AP1-GR inflorescence-like meristems, and fixed using 1% glutaraldehyde in 0.025 M sodium phosphate buffer. Samples were then washed three times with 0.025 M sodium phosphate and further dehydrated in gradually concentrated alcohol solutions. The material was then critical point dried in liquid carbon dioxide, mounted on stubs and coated with platinum. Imaging was performed using a Hitachi 5500 scanning electron microscope.

Chromatin immunoprecipitation
Formaldehyde cross-linking and extraction of nuclei was performed as described in [Carles, 2009 #838] with the following modifications: infiltration time was prolonged (6 min and 24 min), filtration of the extract after grinding was performed through 75 µm and 50 µm nylon meshes and sonication was done for 13x 30 sec on, 1 min off, position H (this longer sonication was only performed for biological replicate 2). Immunoprecipitation reactions were set up as described in [Carles, 2009 #838] except that 15 µg of cauliflower and 10 µg of wild type chromatin was used and the chromatin was incubated with 5 µl of the antibody alone for 5 hours prior to addition of Protein A Agarose beads without salmon sperm DNA (Millipore).
Washing of the ChIP reaction was performed using binding/wash buffer as described in [Carles, 2009 #838] but for four times 10 min at 4°C and once 10 min at room temperature (RT).
Elution and reverse cross-linking was performed as described in ] except that the elution was done by vortex mixing for 1 min at RT and input samples (1/3 of the chromatin used for one ChIP sample) were brought to the same volume by adding elution buffer and Tris pH9 in the same ratio as to the ChIP samples.
Reverse cross-linked samples were purified using the MinElute Reaction Cleanup Kit (Quiagen) and DNA was eluted in 30 µl Elution buffer.
For biological replicate 2, three technical replicates of ChIP on the same nuclear extract were performed and pooled on one MinElute column.

RNA preparation, reverse transcription and gene expression analysis by semi-qPCR
Total RNA was isolated with the RNeasy Plant Mini Kit (Qiagen) and 8 µg of RNA were DNase treated with DNA-freeTM Kit (Ambion). For semi-quantitative PCR (semi-qPCR), 2,2 µg of DNase treated RNA were used for reverse transcription with SuperScript® III RNase H-reverse transcriptase (Life Technologies) and an oligo dT primer (18mer), according to the manufacturer's instructions. From 20 µl of the reverse transcription (RT) product, 1 µl was used for each PCR reaction. The annealing temperature was 54°C for all primer pairs and 25 cycles of PCR were performed for all genes. Primers for semi-qPCR were previously described [Carles, 2009 #838]. Elongation factor 1α (EF1α) was used as an internal control.

Library preparation and NGS
Libraries from ChIP samples of replicate 1 were prepared and sequenced on a GAIIx (Illumina) sequencer as described previously [Yant, 2010.
For biological replicate 2, ChIP-seq library preparation was performed using the Truseq ChIP sample preparation kit (Illumina, ref. IP-202-1012) according to the manufacturer's instructions except that PCR amplification was performed prior to size selection. Briefly, sonicated chromatin was repaired prior to 3' ends adenylation. Illumina's indexed adapters were ligated to the adenylated double-stranded DNA fragments. A 18-cycles PCR was performed on the ligated DNA using Illumina's PCR primers. The libraries were then size-selected between 300 bp and 700 bp on a 2% agarose gel. Each library was validated using high sensitivity chip (Agilent, ref. 5067-4626) on an Agilent Bioanalyzer and quantified using the Kapa library quantification kit (Clinisciences, ref. For RNA-seq, library preparation was performed using the Truseq stranded mRNA sample preparation kit (Illumina, ref.RS-122-2101 ) according to the manufacturer instructions. Briefly, polyadenylated RNAs were purified using oligo-d(T) magnetic beads, fragmented, and reverse transcribed using random hexamers, Super Script II (Life Technologies, ref. 18064-014) and Actinomycin D. During the second strand generation step, dUTP substituted dTTP. This prevented the second strand to be used as a matrix during the final PCR amplification. Double-stranded cDNAs were adenylated at their 3' termini before ligation was performed using Illumia's indexed adapters. Ligated cDNAs were amplified following 15 PCR cycles and PCR products were purified using AMPure XP Beads (Beckman Coulter Genomics, ref.A63881). Libraries were validated using a DNA1000 chip (Agilent, ref. 5067-1504)

ChIP-seq analysis
Reads were mapped to the Columbia reference genome (Tair10) using BWA , files were converted into bam format using samtools  and duplicated reads were removed using Picard-tools (picard.sourceforge.net). Pearson product-moment correlation coefficient (PCC) was calculated for sums of reads per million in 200 bp windows of the genome using a customized c++ script to assess similarity between samples. Simply defining genes as targets or non targets of a chromatin marks is not an adequate approach to study chromatin dynamics between tissues/time points because changes might occur only in a subset of cells and thus a reduction or elevation of a mark would be seen but not a total gain or loss of the respective mark. This would lead to detection of false negative genes in respect of changes. On the other hand, targets are defined by a threshold and thus genes which are very close to the threshold in both conditions, one above and one under it, might lead to the detection of false positives. Thus, we determined differentially enriched regions between two samples using the SICER-df function of SICER V1.1 [Zang, 2009 with window size 200 bp, gap size 400 bp and effective genome size 90 %. This function first determines significantly enriched regions in each sample (ChIP) compared to the respective control. As control, input samples of I were used for the wt samples and input samples of t2 were used for all cauliflower tissues. For all comparisons of samples (ChIP) and control (input), a FDR of 0.0001 was employed. Afterwards, the program determines all regions that are positive in one or both of the tested conditions (the two ChIPs to be compared) and then determines differentially enriched regions among those union islands. Due to the lower number of reads in the first biological replicate, a more lenient FDR was used for the detection of the small quantitative changes between the samples (0.05). For replicate 2, FDR 0.0001 was used. For the H3K4me3 samples, a strong variation in the background levels was observed between samples. To eliminate biases in the normalisation process caused by this fact, SICER-df was applied to files containing only genes on union islands filtered by filter_raw_tags_by_islands.py from SICER as recommended in [Zang, 2009. Genes were considered as significantly enriched (target genes) or significantly differentially enriched when overlapping with the respective regions with their gene body. Only genes found in both replicates were considered for further analysis. Mapping of reads to tair10 genome and determination of differentially expressed genes was performed using TopHat and Cufflinks with default parameters as described in [Trapnell, 2012. Three biological replicates were analysed together.

Further bioinformatics analyses
• Clustering Heat map visualisation and clustering was performed with the software Genesis 1.7.6 [Sturn, 2002]. Normalised expression values in each tissue/ time point (denoted as Fragments Per Kilobase Of Exon Per Million Fragments Mapped (FPKM)) were submitted to the analysis. To visualise patterns of expression dynamics over the time series rather than absolute values, z-scores were employed and calculated from FPKM values of a gene in a given sample as z=(FPKM-µ)/σ, where µ denotes the mean FPKM value for the gene over the whole tissue/time series and σ denotes the standard deviation of this mean. Hierarchical average linkage and k-means clustering were performed as described [Engelhorn, 2010. For the subset of floral regulators (Fig. 1), 5 main branches are visible in the hierarchical clustering tree. Therefore, k-means clustering was performed with k=5. For all DEGs, different numbers of clusters were tested and k=8 yielded the highest diversity of patterns. Distance for all clustering analyses was Euclidean.
• Gene Ontology analysis For a general overview of over-represented functions (GO slim) in gene sets, the Classification Super Viewer from the BAR website (bar.utoronto.ca) was employed [Provart, 2003. For detailed functional analysis over-represented GO terms were determined with the FatiGO tool [Al-Shahrour, 2007 #2524] of Babelomics 4 webservice (babelomics.bioinfo.cipf.es) [Medina, 2010.
• Visualisation of data in a genome browse For visualisation, bam files were converted to bigwig files and normalised by the number of reads in the respective library using awk, bedtools [Quinlan, 2010 and bedGraphToBigWig [Kent, 2010. Bigwig files were visualised in a locally running version of the Generic Genome Browser (GBrowse 2.54) [Stein, 2013.
• Calculation of average histone mark abundance patterns and heat maps over a locus For upstream and downstream regions of a locus, normalised histone mark abundance (expressed in reads per million) was calculated with base resolution and then averaged either in 200bp bins for heat maps or over all loci of interest with base pair resolution for average patterns. Transcribed regions were partitioned in 1500 bins (10 bins for heat maps) to account for variation in gene length and average normalised histone mark abundance was calculated for each bin.

AVAILABILITY OF SUPPORTING DATA
ChIP-seq and RNA-seq data were deposited at the NCBI Gene Expression Omnibus (GEO) under the SuperSeries GSE71583. Normalised H3K27me3 and H3K4m3 profiles expressed as reads per million can be viewed at https://gbrowse.cea.fr (for the moment protected by user:flowdev pw: Ai6aqu8r).

Conclusions
Epigenomic analyses with temporal and tissue-specific resolution has not been reported yet in planta, since most studies were performed at one developmental stage and often in a mixture of tissues, thus providing static and average views of chromatin landscapes. This lack was due to the challenging harvest of tissue-or stage-specific chromatin in sufficient amounts for comparison between samples. To overcome these limitations, we employed the inducible « cauliflower » system that allows production of flower meristems in large amounts and their synchronous development over time. Combination of the floral induction system with quantitative ChIP-seq analysis enabled us to detect changes in histone mark abundance with temporal resolution within the inflorescence, during early flower development. We found that quantitative changes in H3K27me3 and H3K4me3 marks take place within the first two days of flower initiation, affecting several hundreds of genes. During this period, we discovered that H3K4me3 is a stronger predictor of transcriptional changes than H3K27me3, for both activation and repression. Even for PcG target genes, early gene activation is dominated by increases in H3K4me3, and H3K27me3 remains present while H3K4me3 is elevated.
Although it was observed that many floral organ identity genes are up-regulated in the constitutive absence of H3K27me3 depositing enzymes, in the natural context H3K27me3 decline rather occurs after prolonged transcription. This indicates that presence of repressive H3K27me3 marks may not prevent early activation of genes, provided that active H3K4me3 marks are deposited.
In conclusion, we could establish the temporal resolution of chromatin events during the first steps of flower morphogenesis, and identified prevailing changes in H3K4me3 over H3K27me3.

Supplementary Materials:
The following are available online.

Acknowledgments:
We thank Jean-Luc Putaux and Christine Lancelon for access to critical point dryer equipment and Pierre-Henri Jouneau for assistance with scanning-electron microscopy imaging. We thank Béatrice Berthet for help with dissection of ap1cal tissue and Armelle Donglois for her investment in re-ChIP experiments. We are grateful to Laura Gregoire for critical discussion and input on data visualisation and to Leor Eshed Williams for critical reading of the manuscript. Author Contributions: CCC conceived and coordinated the study. JE, RB and CCC planned and performed all experiments except for sequencing library preparation, which was done by DP and HP. FO and MR performed base-calling and demultiplexing of sequencing reads. JE and RB analysed the ChIP-seq data, JE and CK analysed the RNA-seq data and were responsible for all other bioinformatics analyses. JE and CCC wrote the manuscript with the help of RB. All authors carefully read, edited and approved the final version of the manuscript.

Conflicts of Interest:
The authors declare no conflict of interest. The founding sponsors had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, and in the decision to publish the results.

Description of expression clusters among floral regulator genes
Cluster 1 corresponds to regulators of floral transition and inflorescence meristem identity. Their expression is highest in leaves, remains at almost the same level in t0 inflorescence meristem like tissue and then gradually decreases toward I. Cluster 2 correspond to genes regulating inflorescence meristem identity and specifying boundaries. Their expression is low in leaves, high at t0 and then gradually decreases. Cluster 3 comprises genes involved in floral transition, inflorescence meristem and boundaries, which also regulate floral organ development. Expression of these genes gradually rises from L to t2 and decreases in I. Clusters 4 and 5 carry late and early floral organ and meristem identity genes, respectively. Their expression rises from L to I with a main elevation of expression from t0 to t2 (Cluster 5) or from t2 to I (Cluster 4). For most of the floral regulator genes, changes in expression from t0 to t2 are similar to those observed by Wellmer et al. (2006) [Wellmer, 2006] in a micro-array experiment (Supplemental Figure 2).

K-means clustering of DEGs
K-means clustering of all DEGs (14459) with k=8 (Supplemental Figure 3) revealed that the most abundant pattern observed (cluster VII, more than one third of the DEGs) contains genes with a higher expression in leaves than in all other tissues, where their expression is uniformly lower. Three patterns include genes that are lower expressed in leaves than in all other tissues. This is to be expected as leaf tissue fundamentally differs from the three other tissue types (no floral organs, no meristematic tissue, fully differentiated). However, we also find a large number of genes that display elevated expression in inflorescence tissue compared to all other tissues (cluster VI, around one fifth of the DEGs), thus likely to be genes involved in later steps of flower development. The third most abundant group contains genes that are up-regulated in the t0 and t2 tissues and lower in both L and I, where L is lower than I (cluster V).

GO analysis of DEGs
To functionally characterize the DEGs, we performed a Gene Ontology (GO) analysis, starting with a broad GO slim analysis (Supplemental Figure 4 and S5) followed by a detailed, full GO analysis of categories overrepresented in the DEGs lists (Supplemental Table 2). This analysis revealed that the changes in expression pattern correlate with expected functionality of the tissues on a genome-wide scale. Genes activated during the series (from L to I) display over-representation for "developmental processes" in all comparisons except t0 to I. An over-representation for the "TFs" molecular function and the "transcription" biological process in t0 to t2 DEGs indicates that the 2-day induction of flower development affected gene regulatory mechanisms in a broad fashion. Terms related to floral organ development were expectedly enriched when considering genes upregulated from L to t0 or L to I. Importantly, those related to floral whorl, carpel or gynoecium development, already appeared during early induction of flower development (t0 to t2) further confirming the suitability of the chosen time points, while genes up-regulated later (t0 to I) are enriched for terms mainly referring to processes for later stages of development (e.g. pollen development). Plastid-localised proteins were over-represented in leaves and repressed elsewhere or present in "green" inflorescences compared to the white meristem structures present in t0 and t2, (Supplemental Figure 5). Many plastid-localised proteins encoded by genes down-regulated from L to other tissues refer to terms concerning light perception, photosynthesis, sugar metabolism and transport, reflecting a true difference in the photosynthetic capacity and yield. Similarly, terms referring to the response to stimuli are over-represented when comparing down-regulated DEGs in L vs t0 or vs I, but also t0 vs I, indicating that biotic and abiotic responses mainly occur in fully expanded tissues rather than in meristematic tissues.

Quality control for ChIP-seq experiments
We observed a high correlation between our two biological replicates, both at the read level (Pearson correlation coefficient > 0.75 for all time points, >0.9 for t0 and t2) and at the target gene level (overlap >77 % for all time points comprising the most distant L and I, >90% for most) (Supplemental Table 3). The core lists of H3K27me3-targets and H3K4me3-targets (6130 genes and 15088 genes respectively, Supplemental Figure 6) displayed overlaps for of >66 % with published ChIP-chip datasets and >88% with ChIP-seq datasets [Farrona, 2011Roudier, 2011Luo, 2013 #2511] (Supplemental Table 3), confirming the robustness of our dataset.

DMGs during early flower morphogenesis
Among the 13 genes that change in both marks towards activation, three are key regulators of flower development (SEP1, SEP2 and LFY), one is involved in flower organogenesis (SPATULA, SPT) and one participates in SAM development (AINTEGUMENTA-LIKE 6 (AIL6)) (

GO analysis of DMGs
To functionally characterize all DMGs, we determined GO terms that are significantly enriched among DMGs compared to all genes carrying the respective mark (Supplemental Table 5 and 6). The over-represented functions among H3K4me3-elevated genes vary strongly between the different comparisons: L-t0 and L-I are enriched in development-related functions including flower and meristem development, but also in DNA metabolism and cell cycle as expected for a comparison between a differentiated tissue and a meristematic (t0) or a meristem-containing tissue (I). In the two comparisons from t0 or t2 to I (more advanced in flower development), we found enriched terms related to response to light and other stimuli, photosynthesis and sugar metabolism, while during early flower morphogenesis (t0 to t2 comparison), DMGs changing towards activation are especially enriched in transcription and floral development (carpel and gynoecium development) specific terms. This result is in line with the fact that many flower regulators are TFs. Genes showing reduced H3K4me3 for L to t0 or L to I, are enriched in terms related to response to stimulus and depleted in terms related to development, DNA metabolism and transport. Genes with reduced H3K4me3 for t0 to I or t2 to I are enriched in terms related to development, including flower and floral organ development, suggesting that several genes involved in early flower development significantly loose activation marks at later stages. H3K27me3-changing targets mainly represent all types of H3K27me3-related functions (e.g. developmental processes, Supplemental Table 5), for H3K27me3 reduction, only the t0-to-t2 DMGs are strongly enriched in several specific terms for floral organ development. Genes displaying elevated H3K27me3 in L to t0 or L to I comparisons are enriched in metabolic processes, L-to-t0 genes are depleted in multicellular-organismal-development, while t2 to I is flower development.

Analysis of DMGs that are not DEGs
The observed delay in H3K27me3 changes for DEGs suggests that transcriptional changes might be rather a cause than a consequence of H3K27me3 dynamics. This raises the question why several genes show significant changes in H3K27me3 during early flower morphogenesis but not in transcription (Fig. 3A). It might be that changes in those genes were just below the threshold for differential expression. We therefore considered the absolute expression values and early expression changes for H3K27me3 DMGs (Supplemental Figure S10). We found that expression changes in the expected direction are indeed limited to a small proportion of the H3K27me3 DMGs. Interestingly, the majority of genes that do not change in expression are expressed at very low levels or not at all during the complete series. We hypothesise that these genes might gain or loose H3K27me3 maybe through the action of other repressive chromatin features such as DNA methylation or H3K9me3 marking. For H3K4me3, at least half of the DMGs show at least some change in the expected direction.
However, also few changes in the opposite direction are observed, maybe due to interplay with repressive TFs.  Supplementary Table 6. Functional enrichment analysis for H3K4me3 and H3K27me3 DMGs. Significantly enriched or depleted GO terms among DMGs were determined for the L-t0, L-I, t0-t2, t0-I and t2-I comparisons for changes in both directions for both marks (separately and common DMGs). All targets of the respective mark were employed as background distribution. FDR corrected p-values (adj_pvalue) are output of the FatiGO web tool, terms are considered significantly enriched (highlighted in green) or depleted (highlighted in red) compared to all targets of the mark for FDR corrected p-values lower than 0.05. Blank fields indicate that the term is not enriched in the respective sample. If a comparison is not listed for one direction, no significant terms were found for this specific set of DMGs. Elongation factor 1α