To probe the nucleation region dynamics at FLC, we first used our existing model developed for Polycomb silencing generally including mammalian targets (Berry et al., 2017a). Briefly, the existing model was constructed around a series of stepwise transitions between low (H3K27me0) and high (H3K27me3) methylated states at each histone. Transitions to more highly methylated states were catalysed through read/write feedbacks by existing nearby H3K27me2/me3. High levels of H3K27me2/me3 silenced transcription, while transcription antagonised H3K27me2/me3, causing methylation loss by, for example, nucleosome eviction (Schlissel and Rine, 2019). DNA replication caused loss of pairs of histone modifications on a single nucleosome, with a probability of one half. More details can be found in Materials and methods. We found previously that this model could robustly generate bistable expression states, where either low expression/high H3K27me2/me3 or high expression/low H3K27me2/me3 states are stable (Berry et al., 2017a). Note that there are some differences between the mammalian and plant systems, particularly in the rapid addition of H3K27me1 by specialised methyltransferases in plants (Jiang and Berger, 2017). However, we found that incorporating a separate, rapid H3K27me0 to me1 transition in the model that was not catalysed by PRC2 did not substantially affect our results, as our existing model already led to rapid H3K27me1 recovery. We therefore initially retained our original model.

We first simulated our original model, but including a specific, small number of nucleosomes representing the nucleation region, which previous work had identified as consisting of about three nucleosomes (Yang et al., 2017). To confirm the nucleosome positioning within this region, we performed MN (Micrococcal Nuclease) together with ChIP experiments, MNase-ChIP (Figure 1B, Figure 1—figure supplement 1). These experiments revealed the presence of two well-positioned nucleosomes with high MNase-ChIP peaks, and then a further less-well positioned nucleosome ~500 bp downstream of the transcription start site, with a broader MNase-ChIP peak, making up the three nucleosome nucleation region. Most other model parameter values were taken from our previous modelling work (Berry et al., 2017a), supplemented with previous experimental measurements specific to the FLC locus. These latter measurements included the FLC transcription rates in the fully active (Ietswaart et al., 2017) and fully silenced (Yang et al., 2014) states. We also used a cell cycle duration of about one day in the warm, with the cell cycle seven times longer in cold conditions. These and other parameter choices are fully specified and justified in the Materials and methods.

We initially examined the persistence of the metastable memory. Prior experimental data (Yang et al., 2017) for the spreading mutant, lhp1, suggested a metastable memory timescale of tens of days. We collected additional experimental data in the lhp1 mutant, measuring the fraction of root cells with reactivated FLC-Venus levels in the warm after 10 weeks cold treatment. Our results, with higher time resolution in the warm (7 days, 10 days, 14 days) than in Yang et al., 2017 confirmed the long reactivation timescale (Figure 1C), with only ~30% reactivation after 14 days in the warm. We next simulated the metastable memory dynamics in our above model, conservatively assuming that only one sixth of the nucleation region histones in the H3K27me2/me3 configuration are sufficient for full silencing. We assumed here and later that reactivation of only one of the two (independent) copies of FLC in a simulated cell would be sufficient for overall cellular reactivation, with reactivation occurring if the fraction of H3K27me2/me3 dropped below the fully silenced threshold. We started from an initial simulation configuration with three fully methylated nucleosomes in the nucleation region, and all other nucleosomes unmethylated, at the start of an initial post-cold time lag of one day before cell division resumes, and incorporating a further lag of 1 day for reactivation to generate visible protein (Yang et al., 2017) (Materials and methods). To match the observed lack of spreading in the lhp1 mutant, we set the PRC2-mediated methylation rate outside the nucleation region to zero to prevent non-nucleation region methylation. With these rules, we found that the simulated metastable memory did not last nearly as long as in our experiments (Figure 1C), with ~50% reactivation occurring in much less than 5 days. If one third of the nucleation region histones in the H3K27me2/me3 configuration are needed for full silencing, reactivation is even more rapid (Figure 1—figure supplement 2A). Even if we allow our experimental quantification of the FLC-Venus signal to be much less strict (with fewer off cells, Figure 1—figure supplement 3 and Materials and methods), the resulting experimental metastable memory duration is still too long to agree with our simulation results. This simulation analysis improved on previous considerations by simulating the full model with three nucleation region nucleosomes rather an idealised ‘best case’ scenario for memory maintenance (Yang et al., 2017). Indeed, the new simulations showed that the three nucleation-region nucleosomes were even further away from providing adequate memory maintenance than previously reported (Yang et al., 2017). We found that 8–10 nucleation region nucleosomes were required to match the experimental metastable memory duration (Figure 1C), a number which is much too large. Even in this case, the shape of the simulated reactivation profile did not match that found experimentally (Figure 1C). In essence, because of the small number of memory-carrying nucleosomes, the probability that the memory of silencing is lost at DNA replication due to unequal partitioning between the two daughter strands is too large, leading to rapid transcriptional reactivation. The existing model cannot therefore explain how the three-nucleosome nucleation region can metastably store the memory of silencing for tens of cell cycles (Berry et al., 2017a; Yang et al., 2017). Taken together, these results show that the existing model must be amended to handle metastable memory maintenance at FLC.

In order to explain the persistence of the FLC metastable memory observed above, there are two possibilities: either the known memory elements are more stable than previously thought, or there are additional memory elements acting to boost the stability of the memory. For the first possibility, the partitioning of the H3K27me3 memory elements at DNA replication would need to be more ordered, and not entirely random as previously assumed. This alternative is at least in principle feasible, due to the existence of specific factors at DNA replication forks, promoting either leading or lagging strand deposition (Gan et al., 2018; Petryk et al., 2018; Yu et al., 2018), and by evidence for regulated biased segregation in Drosophila germline stem cells (Wooten et al., 2019). For the second possibility, there would need to be additional memory storage elements to H3K27me2/me3. An optimal example of the first possibility would be a strictly alternating distribution of successive nucleosomes from the parental to the two daughter strands, so that at FLC each daughter strand would receive at least one of the three parental nucleosomes, thus more stably preserving the memory. The ordered distribution hypothesis does, however, conflict with previous data on nucleosome dynamics at DNA replication (Annunziato, 2005). These data were not collected specifically for Polycomb targets, so it remains formally a possibility that at appropriate Polycomb-target genes the nucleosome distribution after DNA replication could be more ordered. However, as yet we know of no evidence specifically pointing in this direction. For completeness, we analyse this alternative model in the Materials and methods (Figure 1—figure supplement 2B). Even in the best-case alternating distribution case, we find that more than three nucleation region nucleosomes are needed (Figure 1—figure supplement 2B), more than observed experimentally, and again the shape of the simulated reactivation profile does not match that found experimentally. For these reasons, we do not favour models based solely on ordered nucleosome partitioning.

The second possibility is that there are extra, non-histone modification memory storage elements in the FLC nucleation region. These elements could stabilise the metastable memory due to their greater number, which would reduce the probability that they will be so unevenly partitioned at DNA replication that memory on one daughter strand would be lost. For these extra elements to contribute, they cannot be other histone modifications, as at DNA replication if a nucleosome is distributed to the other daughter strand then all histone marks are lost.

A further conceptual problem concerns how histone modifications can first be deposited during nucleation in the cold, as the read/write mechanism requires H3K27me2/me3 to already be present before feedbacks can add more. A low level of H3K27me3 can be generated by noisy, background processes, but these could be insufficient to effectively nucleate a H3K27me3 peak in the cold. Accordingly, we hypothesise that the new elements are also needed to initially create a H3K27me3 nucleation peak in the cold and must therefore be deposited in the nucleation region before H3K27me3.

To satisfy these stringent constraints, we propose a protein assembly forming at the nucleation region, as we detail below. Clear candidates that could make up such an assembly are PRC2-interacting proteins, potentially the VIN3, VERNALIZATION 5 (VRN5) and VEL1 proteins that immunoprecipitate with PRC2 (Wood et al., 2006; De Lucia et al., 2008). These proteins contain PHD fingers proposed to facilitate recruitment at FLC (Kim and Sung, 2013) (see below) and a protein interaction (VEL) domain at their C-terminus, which facilitates homo and heterodimerization in yeast and plant cells, which was linked to histone modification levels (Greb et al., 2007). VIN3 expression is upregulated by cold (Sung and Amasino, 2004), and both VIN3 and VRN5 proteins associate with the nucleation region at FLC (Sung and Amasino, 2004; Yang et al., 2017). Longer cold leads to higher VIN3 levels and hence a potentially higher local concentration, leading to establishment of the H3K27me3 nucleation peak at that locus. In the following, we develop a model based on VIN3/VRN5, although we emphasise that other proteins may also be involved.

Polycomb proteins are recruited to their targets (here, the nucleation region) through a variety of mechanisms: sequence-specificity, transcription factors, lack of CG methylation, histone modifications and DNA structure (Skourti-Stathaki et al., 2019; Yu et al., 2019). We propose that a combination of these factors brings the putative assembly subunits to the nucleation region. Here, for example, the upstream factor VIVIPAROUS1/ABI3-LIKE 1 (VAL1) can bind to the nucleation region due to sequence specificity and indirectly interact with VIN3 (Qüesta et al., 2016). When a sufficiently large number of VIN3 proteins (here, four) are present by chance simultaneously, a nucleation threshold is overcome, and a metastable, dynamically exchanging assembly forms. The assembly can then attract more proteins through strong positive feedback, counterbalancing monomer loss in a dynamic steady state. Hence, the assembly is a digital phenomenon: either present or not present at any given locus, consistent with the properties of metastable silencing in individual cells (Angel et al., 2015; Yang et al., 2017). Once formed, the assembly stimulates catalysis of H3K27 methylation through stimulation and/or recruitment of PRC2 factors (Figure 2A), thereby forming and then maintaining the H3K27me3 nucleation peak.

Figure 2 with 2 supplements see all

Download asset Open asset

The act of nucleation thus involves formation of the assembly and subsequent addition of H3K27me3. If the nucleation-region histone modifications are lost at DNA replication and replaced by newly incorporated, unmarked nucleosomes, then surviving members of the assembly can recreate the H3K27me3 peak. We assume here that components of the assembly persist through DNA replication at the locus and are also randomly distributed to the daughter strands at DNA replication similar to, but separate from, the random nucleosome redistribution. Less random segregation might also be possible (as we consider explicitly for the nucleosomes). However, an equal distribution to the two daughter strands is less likely as this could generate memory that is too persistent. In this way, assembly components also constitute epigenetic memory storage elements, conceptually similar to histone modifications. Passage of the replication fork will certainly disrupt the assembly, but if members of the assembly can be transiently retained around the replication fork, rather than diffusing away, then seeding of two daughter assemblies on the two daughter strands could be efficient. For this reseeding process to be effective, we note that the assembly size before replication must be more than twice the minimal size needed for a stable assembly, due to an approximate two-fold dilution at DNA replication. In our fitted model below this constraint is satisfied, as the assembly size (around 17) is much bigger than twice the minimum assembly size (fixed at 4). Only in rare cases where enough assembly members are lost at DNA replication, so that the assembly copy number drops close to or below the nucleation threshold, will nucleation sometimes be lost. Note that members of the assembly feedback on themselves and onto nucleation region H3K27me3 (on all nucleation region histones equally), but we assume that H3K27me3 does not feedback to recruit extra members of the assembly. The rationale for this choice will be explained later. We also allow for H3K27me3 to feedback to add more H3K27me3 on nearest-neighbour histones in the nucleation region in an allosteric read/write mechanism, but at lower rates. This could occur through PRC2 being constitutively recruited close to the nucleation region independently of H3K27me3 consistent with the methyltransferase CLF being detectable by ChIP in that region before cold (Yang et al., 2017) and then being allosterically activated by nearby H3K27me3 (Materials and methods).

We next used existing data on the rise of VIN3 protein levels during cold to help parameterise our model of nucleation region dynamics, which showed that VIN3 levels rise linearly with time in the cold (see Figure 2B and Hepworth et al., 2018). We also allowed for a reduction in the push to transcribe in the cold to represent the non-epigenetic antisense ('COOLAIR') pathway which downregulates FLC expression (Materials and methods). This effect encourages H3K27me3 nucleation by reducing the antagonising effects of transcription (Berry et al., 2017a). We also assumed a short lifetime for VIN3 in the nucleation region, with a timescale of minutes, qualitatively motivated by the rapid recovery of PRC2 subunits previously observed by FRAP (Youmans et al., 2018). Hence, the nucleation region assembly, once formed, is a highly dynamic object. Early in the cold, VIN3 levels are low, which leads to a low probability that sufficient VIN3 proteins come together to initiate the assembly; naturally explaining the low nucleation probability at that time. Interestingly, we found that later in the cold, we could only fit the further rise of the H3K27me3 nucleation peak by assuming that eventually VIN3 is no longer limiting for nucleation and that the probability of nucleation saturates at a high enough VIN3 concentration. Without this assumption, it was not possible to fit the sigmoidal form of the population-level H3K27me3 nucleation peak rise. The model therefore predicts that other components are also needed for nucleation that eventually become limiting after a long enough duration of cold. The nucleation region model (see Figure 2A, Figure 2—figure supplement 1 and Materials and methods for a full description) could then explain the population-level rise in H3K27me3 levels in the nucleation region during the cold (Figure 2C): as time passed in the cold, more and more loci digitally nucleate the assembly helped by rising VIN3 levels, leading to a digital H3K27me3 nucleation peak, and therefore a smoothly rising H3K27me3 peak at a cell population level. In these simulations, we used a reduced value of the PRC2-mediated methylation rate $k_{me}^{low}$ outside the nucleation region to inhibit non-nucleation region methylation in the cold, as observed experimentally.

We could also now explain the persistence of the metastable memory in the warm following cold in the lhp1 mutant. During this phase, no further de novo nucleation was permitted in the model, consistent with VIN3 being depleted in experiments within hours of reintroduction to warm temperatures (see below) (Hepworth et al., 2018). We also set the PRC2-mediated methylation rate outside the nucleation region to zero to prevent non-nucleation region methylation and initialised the simulation with three fully H3K27me3 occupied nucleation region nucleosomes, as we did previously for Figure 1C, as well as the maximum number of assembly members (Materials and methods). We also assumed that one third of the nucleation region histones in the H3K27me2/me3 configuration are sufficient for full silencing. With the extra non-histone modification memory elements, we found that the nucleation region memory was preserved for long enough to match our experimental data (Figure 2D), with the best fit requiring 17 extra memory elements in the assembly. The fit of the simulated reactivation profile to the data was also much better than that found previously in Figure 1C, with slower initial reactivation which then accelerates with time in the warm. This shape is due to the need to first lose the assembly before H3K27me3 can then be lost; such a two-step process leads to a sigmoidal reactivation profile. In addition, the model was also able to explain the low levels of H3K27me3 experimentally observed in the nucleation region prior to cold (and arrival of VIN3): the initially low levels of H3K27me3 were unable to generate enough allosteric activation through constitutive PRC2 activity to allow these levels to increase. Note that the fit in Figure 2D while good is not perfect, indicating that our model may still be too simple. One possible extension is a combined model with both ordered nucleosome distribution and extra memory elements. This generalised model can also explain our data, using a strictly alternating nucleosome distribution at DNA replication and a best fit reduced number of eight memory elements (Figure 1—figure supplement 2C). However, we emphasise that the extra memory elements are still required.

One potential difficulty of the extra memory element model for FLC is that the VIN3 protein disappears rapidly after the end of cold treatment, yet the metastable memory is persistent for many days in experiments in the subsequent warm, requiring that the assembly memory component also persists in the warm. Hence, other factors must be involved, potentially the constitutively expressed VIN3 homologues, VRN5 and VEL1 or components of a PRC1 complex. For simplicity in the model, we model a single protein unit capable of positive feedback once a critical local concentration has been reached at the nucleation region. However, consistent with nucleation only occurring during the cold, spontaneous binding to the nucleation region is only allowed in the model during the cold, but with positive feedback possible at all times both in the cold and warm. With this reasoning, the hybrid assembly/H3K27me3 model can explain all our data on nucleation and persistence of metastable memory.

So far, we have dissected the dynamics of the nucleation region, in the cold and then in the warm in the lhp1 mutant. However, in the wild type, the H3K27me3 domain expands from the nucleation region to cover the entire locus after transfer to warm conditions, in a process we term spreading from the nucleation region to the body region. Up to now, we have limited this spread, as appropriate for the wild-type in the cold and lhp1 mutant in the warm. We now relax this restriction for the wild-type, using a value $k_{me}^{high}$ of the PRC2-mediated methylation rate outside the nucleation region in the warm, and examine in more detail how this spreading process occurs.

Interestingly, at all times when present, both in the cold and in the subsequent warm, VIN3, VRN5 and the PRC2 methyltransferases CLF and SWINGER (SWN) are predominantly located in the nucleation region (Yang et al., 2017). There is some limited spreading of these factors but only to significantly lower levels. Although alternatives such as nucleosome mobility have previously been proposed (Chory et al., 2019), we consider here the possibility that spreading is mediated by a looping interaction between the nucleation region-located methyltransferases and nucleosomes in the rest of the FLC gene body. This interaction would be allosterically activated by H3K27me2/me3 on the target nucleosome. During this interaction, the key proteins maintain their presence in the nucleation region consistent with previous protein ChIP data. Such a mechanism is also consistent with PRC2 being present at unexpectedly low levels, even in cases with active cell cycling, with only one methyltransferase for every ~20 H3K27me3 marks in Drosophila (Bonnet et al., 2019). Localising the methyltransferases in discrete regions, and then allowing for a longer-range looping interaction to a larger domain, can clearly make efficient use of limited number of methyltransferases to reliably coat extended domains with H3K27me3 (Oksuz et al., 2018).

We therefore assumed the constitutively present methyltransferases CLF and SWN close to the nucleation region were able to loop to contact with other body-region nucleosomes at the FLC locus and, if H3K27me2/me3 was present on the body region nucleosome, then more H3K27me could be allosterically added to that nucleosome or its nearest neighbours. Although, in this interpretation, the addition of H3K27me is meditated by looping from the nucleation region, mathematically this process turns out to be identical to the local spreading process implemented in Berry et al., 2017a. This is because H3K27me2/me3 must be present on the target body nucleosome, before more H3K27me can be added to that nucleosome or its neighbours. Hence, effectively it is this neighbouring H3K27me2/me3 that catalyses the addition, with the looping interaction simply enabling the reaction to proceed. Note that we assume here, for simplicity, that such looping could occur anywhere within the locus with equal probability (see below), but not outside, and that this probability was much higher in the warm than in the cold in the wild-type. Such looping could be mediated, at least in part, by the formation of discrete LHP1 foci, as mutations in lhp1 (in the RNA binding domain) that compromise the foci also compromise efficient spreading (Berry et al., 2017b). Within the LHP1 foci, we suppose that looping interactions are relatively unconstrained, facilitating the spreading process, but that these interactions are much reduced outside. We emphasise that the looping interaction and addition of H3K27me in the body is not directly enhanced by the assembly: the assembly’s only effect on histones is to enhance H3K27me3 in the nucleation region. Hence, with the exception of the assembly dynamics itself, and its role in enhancing nucleation region H3K27me2/me3, the underlying model is similar to that of Berry et al., 2017b.

With these considerations, we then simulated the full nucleation and spreading dynamics in the ColFRI wild-type (Figure 3A). Beginning in the warm, the initially low level of H3K27me3 ensured limited allosteric activation and therefore continued low levels of H3K27me3 across the whole locus. In the cold, nucleation occurred, as discussed above. With a low level of looping in the cold, low levels of spreading across the locus could be generated, as observed experimentally. On return to simulated warm, looping was then much more extensive, and it was then possible to generate H3K27me3 spreading dynamics in the warm similar to previously published ChIP data. Note that we do not explicitly fit to warm ChIP data, because the spreading dynamics at a whole plant level (as assayed by ChIP) is complex, as only cells with an active cell cycle can spread (Yang et al., 2017). At subsequent DNA replication events, any dilution of H3K27me3 in the gene body was filled in by the looping interaction/allosteric activation before the next round of replication (Figure 3A,B). In parameterising these processes, we considered prior data from tobacco cells where H3K27me3 dilution at DNA replication was reversed after only ~2 hr (Jiang and Berger, 2017), much faster than in mammalian systems (Alabert et al., 2015). For our Arabidopsis FLC simulations, our best fit used a timescale intermediate between these cases. Slower rates resembling mammalian systems were unable to generate effective spreading, with limited spreading diluted away by DNA replication, contrary to experiments.

Figure 3

Download asset Open asset

Our previous ChIP results in the warm demonstrate that even in a fully spread state in the standard laboratory wild-type ColFRI, the population-average H3K27me3 levels in the nucleation region drop to those of the body (Yang et al., 2017; Qüesta et al., 2020). We previously referred to this new state, with lowered nucleation region H3K27me3, as the ‘perpetuated’ state (Qüesta et al., 2020), as in ColFRI, this state is stable.

These results are straightforward to understand within our model. As mentioned previously, assembly positive feedback assists in the addition of H3K27me but not vice versa. Hence, eventually, the assembly will suffer stochastic subunit loss at DNA replication to below the nucleation threshold. Due to the lack of feedback from H3K27me, the remaining subunits will then disperse (Figure 3A,B,top). This then leads to a reduction in the H3K27me2/me3 feedbacks and levels in the nucleation region (as observed experimentally Yang et al., 2017; Qüesta et al., 2020) potentially due to the lack of assembly proteins to boost PRC2 activity (Figure 3A,B, middle). Note that this reduction is only experimentally observed to occur in the nucleation region not elsewhere: hence, the model was constructed such that the assembly only positively feedbacks onto H3K27me3 in the nucleation region and not elsewhere, as mentioned previously. The methyltransferases CLF and SWN are, however, still present in the nucleation region, due to constitutive binding, as revealed by ChIP data (Yang et al., 2017). Therefore, the addition of H3K27me across the locus through the looping interaction is still present. Due to the higher levels of H3K27me2/me3 across the locus, we found that CLF/SWN activity from the nucleation region can still perpetuate the silent state using long-distance looping/allosteric activation, even without assembly proteins replenishing the nucleation region H3K27me2/me3. Such feedback is not present at such high levels during warm periods prior to cold exposure due to the lack of sufficient H3K27me2/me3 to permit significant allosteric activation, as discussed previously. This inherent bistability ensures that there is no premature rise in H3K27me3 levels before cold.

There are many parallels between the FLC nucleation region and Polycomb Response Elements (PREs), as both are required in their respective contexts for epigenetic state switching. We therefore asked if further insight could be gained by treating the FLC nucleation region as an effective PRE. If so, our model should capture the results of PRE (nucleation region) excision experiments, where epigenetic memory is lost after 1–6 further cell divisions after excision (Coleman and Struhl, 2017; Laprell et al., 2017). Indeed, we recover very similar results in our model (Figure 4A). The nucleation region, with its constitutively bound CLF/SWN, and its looping interaction, are essential for maintenance of the silenced state in the spread and perpetuated states. As soon as the PRE nucleation region (with any assembly elements) is removed, the feedbacks are greatly diminished (Figure 4A, red, yellow curves) and the silenced state decays over the subsequent few cell cycles. In this case, the decay of the silencing is somewhat slower than would be predicted just from nucleosome dilution at DNA replication (Figure 4A, pink curve) as we allow here for the effect of limited PRC2 binding throughout the gene body which is able, at a low level, to add further H3K27me through standard read/write feedbacks. This scenario has also been reported experimentally (Coleman and Struhl, 2017). However, this feedback mechanism can easily be adjusted to be too weak to perpetuate the state fully, in agreement with PRE excision experiments.

Figure 4 with 1 supplement see all

Download asset Open asset

To explore the model further we also simulated the Arabidopsis variety Lov-1 from northern Sweden, that reactivates FLC expression most markedly if it is not exposed to sufficiently long cold (Qüesta et al., 2020), in contrast to the standard laboratory wild-type ColFRI. In previous work, we traced this reactivation to an instability in the perpetuated state. One possible explanation for this reactivation would be that FLC in Lov-1 has an intrinsically higher transcriptional activity in the warm, with the higher antagonistic activity eventually overwhelming the H3K27me3-mediated silencing. In support of this idea, Lov-1 differs in cis from ColFRI at FLC through Single Nucleotide Polymorphisms (SNPs) specifically in the nucleation region, SNPs that could alter trans factor recruitment strengths. To further test this hypothesis, we therefore simulated stronger transcriptional activation (α) in warm conditions (Figure 4B,C) for Lov-1, as well as slightly altered nucleation region parameters in the cold (Materials and methods), again consistent with nucleation region SNPs. The stronger transcriptional activation, together with the shorter cell cycle in the warm versus cold, caused a loss of silencing at the population level in the warm. The eventual H3K27me3 loss occurred through H3K27me3 levels dropping too low at a given locus after DNA replication to be recoverable in the subsequent cell cycle, due to the stronger antagonistic effect of transcription (Figure 4C). Consistently, in simulated non-replicating Lov-1 cells (Figure 4—figure supplement 1), the silencing was stable. However, this reactivation could only occur after loss of the assembly, as otherwise the assembly memory elements could simply re-nucleate H3K27me3. Accordingly, we found there was a delay before transcriptional activation could occur, as proposed in Qüesta et al., 2020. In this way, through subtly impaired bistability due to increased transcriptional activation, the model was able to give a simple explanation for the eventual reactivation behaviour observed in replicating cells in Lov-1, but not in ColFRI.

We simulate the status of H3K27 methylation and of the protein assembly at the FLC locus by expanding a previous mathematical model of a generic Polycomb target gene (Berry et al., 2017a; Figure 2—figure supplement 1). The simulations are conducted similarly, using the direct Gillespie algorithm (Gillespie, 1977). The details of the simulations are elaborated in this section.

In the extended model of this work, in addition to the status of the histone modifications, we also simulated how many proteins are bound in an assembly in the nucleation region of the gene. We therefore simulated both histone modifications and proteins (Table 1). The histone modifications are further divided into two regions, whether they are located in the nucleation region (NR) or body region (B).

In the simulations, there are six possible events: (1) the addition of a methyl group to a histone H3K27, (2) the removal of a methyl group from a histone H3K27, (3) the binding of a protein to the assembly (4), the unbinding of a protein from the assembly, (5) transcription, and (6) DNA replication.

When a methyl group is added to a histone H3K27 in the simulation, an unmethylated histone transitions to monomethylated, monomethylated to dimethylated, etc. The methylation propensity $r_i^{\text{me}}$ (Table 3) is divided into two parts. $b_i^{\text{me}}$ describes methylation addition from methyltransferases bound to the nucleation region that can loop to interact with distant H3K27me2/me3 and then add more H3K27me to the same or nearest neighbour nucleosomes. This term also includes addition from methyltransferases that can bind to H3K27me2/me3 and add more H3K27me to the same or nearest neighbour nucleosomes. In the model of Berry et al., 2017a, only the second of these possibilities was considered. Here, we generalise to add the first possibility, noting that the mathematical implementation of the model remains unchanged by this reinterpretation. We also allow for 'noisy' addition of methylation without the need for pre-existing H3K27me2/me3, but at a much lower level, as in Berry et al., 2017a. The second component of the methylation propensity, with parameter k_p-me, is due to feedback from the protein assembly to nucleation region H3K27 methylation, which boosts H3K27 methylation addition in the nucleation region only, provided the assembly is present.

The removal of methyl groups, demethylation, is described by $r_i^{\text{dem}}$. When demethylation is carried out a trimethylated histone is converted to dimethylated, dimethylated to monomethylated etc.

The binding of a protein to the assembly is described by the propensity $r^{bind}$ (Table 3). The assembly proteins could consist in part of VRN5 and, in the cold, VIN3. VIN3 levels are known to rise approximately linearly with increasing duration of cold, and are hypothesised here to be required for assembly formation, as VIN3 mutants are defective in nucleation. However, once the assembly has formed, we expect that VRN5 and potentially other proteins can also participate in assembly dynamics, as discussed in the main text. In the warm, VIN3 levels drop very rapidly and so VRN5 and others are likely principal constituents of the assembly in the warm following cold. To simplify the model, we consider only a single composite protein with the properties of a combined VIN3/VRN5. In the cold only, the binding propensity builds up with the number of weeks of cold (representing the VIN3 contribution, ${\gamma }_{bind}\left(w\right)$ in Table 2) and saturates with a timescale $w_0$ weeks in the cold. Once formed, with the number of protein subunits above a critical threshold (equal to v in the model), positive feedback arises and contributes a much larger constant binding term (potentially representing the VRN5 contribution, $k_{bind}G$ in Table 3, with $k_{bind}$ a constant in all conditions). In this way nucleation of the assembly is a slow process in the cold (as observed experimentally), but once formed, strong positive feedback gives effective maintenance. When the binding of a protein to the assembly occurs in the simulations, the status of one site in the protein assembly occupancy changes from unbound to bound. The size of the assembly is the number of sites with a bound protein.

The unbinding of a protein from the assembly is described by the propensity $r_k^{unbind}$, and is the same in all conditions. For unbinding of protein k, the appropriate assembly occupancy status changes from bound to unbound. If sufficient proteins are lost from the assembly through $r_k^{unbind}$ or DNA replication (see below), the protein number may drop below the critical threshold v, thereby abolishing the strong feedback. In the warm, the remaining bound proteins will then rapidly disperse due to the lack of binding and the assembly is lost. In the cold, an assembly can still potentially reform due to the background cold-induced binding potentially via VIN3.

Transcription is described by the stepwise linear function f (Table 3), with maximum value f_lim. Previous work has shown that full FLC repression can be achieved based only on nucleation region (NR) nucleosomes (Yang et al., 2017). However, in the warm, once the assembly has dispersed, silencing must also be capable of being maintained by methylation on non-NR nucleosomes. This must be the case, as silencing would be deficient if only NR region nucleosomes enforced silencing in the case where the small number of NR-nucleosomes were not sufficiently inherited at DNA replication. Hence, we used a combination of both NR and non-NR region nucleosomes to control silencing. Whichever of the two regions had the higher H3K27me2/me3 level, that was the level that dictated the level of FLC silencing achieved, as described in Table 3. When transcription is carried out in the simulations, the histones are exposed to demethylation and histone exchange. The histones are demethylated (as described above) with a probability of $p_{dem}$ per histone and exchanged with a probability of $p_{ex}$ per histone. When a histone is exchanged, both histones of a nucleosome are changed to unmethylated (histone 𝑖 and 𝑖+1 if 𝑖 is odd, histone 𝑖 and 𝑖−1 if 𝑖 is even).

When replication is carried out, with a probability of 0.5 per nucleosome, the nucleosome is exchanged to a nucleosome with both histones unmethylated. In addition, with probability 0.5 for each site with a bound protein in the assembly, the assembly occupancy status is changed from bound to unbound. Thus, each individual nucleosome, and separately each individual protein in the assembly, is inherited through DNA replication with 50% probability. DNA replication is carried out when the next event in the Gillespie simulation is drawn to occur at an elapsed time since the last DNA replication equal to, or greater than, the cell cycle duration. This duration differs between cold and warm conditions (see Tables 4 and 5).

In Figure 1—figure supplement 2B, replication is simulated with an ordered, equal distribution of the nucleosomes to the daughter strands. After replication every second nucleosome is replaced with a nucleosome with both histones unmethylated. Whether odd or even numbered nucleosomes are replaced is chosen randomly with a probability of 0.5. In Figure 1—figure supplement 2C, ordered nucleosome distribution is combined with the assembly model described above.

The propensities are described in detail in Tables 2 and 3. Further justifications of the equations are described in Berry et al., 2017a. Note that with the exception of the assembly dynamics, and how the assembly affects H3K27me addition in the nucleation region, together with the control of silencing by both NR and non-NR region nucleosomes, the underlying model is identical to that in Berry et al., 2017a, although with altered parameter values (see next section).

Some parameters of the model vary between cold and pre-/post-cold conditions; these variations are specified in Table 4. The cell cycle length is different between cold and warm, consistent with the observed decrease in growth when plants are in vernalization conditions (Zhao et al., 2020).

The strength of trans-activation ($\alpha$) is also changed between the pre-cold/post-cold and cold, where, post-cold, $\alpha$ is changed after a time lag of 1 day. The change in $\alpha$ in the cold is due to an observed decrease in transcription in the cold that is not caused by an increase in H3K27 methylation (Hepworth et al., 2018). This effect is believed to be mediated by an antisense transcription mechanism which we do not model in detail here. Instead we simply incorporate its outcome, a reduction in the push to transcribe in the cold.

The methyltransferases CLF and/or SWN are always present in the nucleation region in all conditions (Yang et al., 2017). Therefore, the methylation rate $k_{me}$ is set to a fixed value both in the NR and outside at all times, with one exception. The exception is a reduced rate $k_{me}^{low}$ in the cold outside of the NR. This is required to supress premature spreading of H3K27me2/me3 away from the NR in the cold. The origin for this suppression is not yet mechanistically clear. However, in our interpretation, the NR can loop over to non-NR nucleosomes and catalyse the spread of H3K27me2/me3 from a non-NR nucleosome with H3K27me2/me3 to the same nucleosome or its nearest neighbours. Either the looping interaction itself, or the strength of the catalytic reactions, may be reduced by cold temperatures.

The model parameters that are fixed throughout for all simulations are listed in Table 5, with references for where they are justified. Most parameter values are found in Berry et al., 2017a; however, a few of them are changed to adapt to the features observed at FLC. The parameters $f_{min}$ and $f_{max}$ are redefined (Table 5), and therefore new values of $p_{dem}$ and $p_{ex}$ used to ensure maintenance of low H3K27me3 levels pre-cold, but high H3K27me3 levels post-cold (Figure 2—figure supplement 2).

Free parameters are listed in Table 6. If $N_{\text{NR}}$, $N_{\text{p}}$, $N_{gen}^x$ and $N_{gen}^y$ are not specified, the values in Table 6 are used. The values of $k_{\text{me}}$ were chosen in order to ensure no spreading to the body region in the cold, but to efficiently spread post-cold. The value of N_NR was chosen to fit the width of the nucleation region (Figure 1B; Figure 1—figure supplement 1), while N_P and the other parameters were chosen to fit the rise of the nucleation peak in the cold, as well as the stability of the nucleation peak/spread state after the cold in the lhp1 mutant and wild-type, respectively. We also varied our fitted parameters by 10% and found that our conclusions were all qualitatively unchanged by these variations.

The definition and simulations used to find the bistability measure $Bi$ (Figure 2—figure supplement 2) are the same as in Berry et al., 2017a. $Bi$ is calculated from the combined output from the simulations of a system of 60 histones initialised in each of the uniform H3K27me0 or H3K27me3 states and simulated for 50 cell cycles. $P_{OFF}$ is the proportion of the time where the system is in a high H3K27me2/me3 state and $P_{ON}$ is the proportion of time it is in a high H3K27me0/me1 state (definitions in Table 7). The value of $Bi$ is averaged over 4000 simulations for each parameter set. A value of $Bi$ close to 1 is regarded as a bistable simulation. The simulations are without the protein assembly to ensure the histone feedbacks alone are able to maintain the spread state after the assembly has disappeared.

The system is initialised with no proteins and all histones in me0. Before cold, the system was equilibrated for 10 days. This duration was sufficiently long to allow the system after 10 days to generate methylation and protein levels as a function of time from the previous DNA replication event, that were invariant from one cell cycle to the next. To match our experiments, the value of $N_{gen}^x$ was chosen to match a total of 8 weeks cold exposure, with $N_{gen}^y$ representing 20 further days in the warm.

When the lhp1 mutant is simulated starting at the end of the cold, the nucleation region is initialised with all histones with me3, and the body region with all histones me0 (Figure 1C, Figure 2D, Figure 1—figure supplements 2 and 3). When the protein assembly is included (Figure 2D, Figure 1—figure supplement 2C), the simulations are additionally initialised with the maximum possible number of proteins bound in the assembly. In these cases, the noisy methylation transition rates are also reduced (see Table 9) to prevent noisy addition of methylation from causing silencing, as would otherwise be possible due to the small size of the nucleation region. For these reactivation simulations, we simulated two independent copies of FLC for a cell, and assumed that reactivation of only one of the two copies would be sufficient for overall cellular reactivation, with reactivation (FLC-ON) occurring if the fraction of H3K27me2/me3 dropped below the fully-silenced threshold.

In Figure 4, PRE excision and the Lov-1 accession are simulated. For the PRE excision, the simulations are mostly the same as earlier; however, there are changes after cold treatment: for all remaining time after 12 days (including the 1 day without replication, see below), the nucleation region is changed to unmethylated, the number of proteins bound to the assembly is set to zero, $k_{me}$ in the nucleation region is set to zero and the $k_{me}$ in the body region is decreased to $k_{me}^{ex}$, as specified in Figure 4A. This is summarised in Table 8. For the simulations of Lov-1, the parameters that are different from ColFRI are specified in Table 9.

When we show the simulated H3K27me3 levels and fraction of proteins in the assembly, we average the time course over $N_{sim}$ simulations (Table 7). In the Gillespie algorithm, the time steps are not constant. We therefore convert the time steps and the corresponding H3K27me3/protein assembly levels to constant time steps. We do so by choosing a short time step $t_{out}$ and assign it to the H3K27me3 levels that was recorded at the end of that timestep. To avoid having all the simulations synchronised, so that, for example, DNA replication does not occur at the same time for all simulations, each individual simulation is started at a random time in the cell cycle. When the simulation shifts from warm to cold conditions, we convert where the system is in the cell cycle to where it would be in the cold cell cycle, that is if there is 2 hr left before replication in the warm cell cycle, it is then 14 hr left before replication in the cold cell cycle. An analogous rule is implemented upon return to warm conditions, but only after a time lag of 1 day without replication is first simulated, to implement a delay while the plant adjusts to warm temperatures. Then, the cell cycle is converted to the warm cell cycle, that is if there is 14 hr left before replication in the cold cell cycle, it is then 2 hr left before replication in the warm cell cycle. During the 1-day lag, all parameters adopt their warm values (Table 4) except for the transcription activity parameter α which retains its cold value, to allow sufficient post-cold H3K27me2/me3 spreading.

For the lhp1 mutant simulations, the one-day time lag without replication is also implemented, as above, with the reduced value of α, after which each simulated locus enters the standard warm cell cycle at a random phase. An additional 1-day lag for protein to become visible is also included.

The source code of the model simulations is available at https://github.com/ceclov/hybrid-PRC2-model copy archived at swh:1:rev:953304772e7cec326736b0c72693ba698bb2540a (Lövkvist, 2016). 

The ChIP data in Figures 2C and 3A are from Yang et al., 2014. The data in Figure 4B are from Qüesta et al., 2020. To measure H3K27me3 levels for the nucleation region and the body region separately, primers for those specific regions were used (Table 10) and the average value separately calculated over each region. The sequences and positions of these primers are specified in the corresponding Supplementary Information of the above references. For the nucleation region and body regions separately, the lowest value before or during the cold is then subtracted from each time point. The resulting values for the nucleation and body regions are then normalised to the value of the nucleation region at the longest cold treatment (8 weeks in ColFRI, 12 weeks in Lov-1). This procedure accounts for loci in our experiments that have nucleated/spread before cold treatment and allows us to focus on loci that are un-nucleated at the start of cold, as appropriate for our simulations.

Note that we only use ChIP data during but not after the cold to constrain the models. This is because the increase in H3K27me3 levels in the warm seen in the body region requires an active cell cycle (Yang et al., 2017). Since only a subset of the cells in the plant are actively cycling and dividing, fitting the rise in body H3K27me3 requires knowledge not only of the intrinsic spreading process at a single locus, but also of the relative fraction of cells that are dividing and therefore capable of undergoing spreading. Due to this additional complexity, which cannot yet be fully controlled, we only fit to ChIP data in the cold, where there is less division and so these complications are reduced. Note that this complication does not apply to the root imaging data, since all the cells being imaged are actively cycling and dividing.

Experimental data in Figure 1C, Figure 2D, Figure 1—figure supplements 2 and 3 were processed as described in Yang et al., 2017 and as described below. VIN3 expression levels in Figure 2B are from ColFRI with 8°C cold treatment (Hepworth et al., 2018).

In this section, we describe experimental details of the imaging and quantification of FLC expression in lhp1 (Figure 1C, Figure 2D, Figure 1—figure supplements 2 and 3), as well as for the ChIP and MNase experiments (Figure 1B and Figure 1—figure supplement 1).

When lhp1 is used it refers to lhp1-3 which is described in Larsson et al., 1998. To image FLC levels in lhp1, a transgenic line with FLC-Venus in the FRI^sf2 background (described in Yang et al., 2017) was used. For the MN-ChIP experiment, ColFRI was used.

For both experiments, the seeds were vapour-phase sterilised as described in Clough and Bent, 1998 and sown on Murashige and Skoog (MS) medium plates without glucose (see Table 11 for individual experiments). The plates were kept at 4°C in the dark for 2 days for stratification. Before vernalization, or for non-vernalized conditions (NV), the plants were grown in long-day conditions (16 hr light, 8 hr darkness with constant 20°C). When plants were vernalized, they were first pre-grown in NV conditions (number of days specified in Table 11) and then transferred to vernalization conditions (8 hr light, 16 hr darkness with constant 5°C). For the lhp1-imaging, the plants were further transferred back to long-day conditions.

The roots were imaged on microscope slides (VWR, 631–0117) with No. 1,5 coverslips (VWR, 631–0138). The slips were sealed using a frame of adhesive tape (0.12 mm thickness, Sigma-Aldrich, GBL620001-1EA). The coverslip was pressed on the adhesive tape to prevent drying during imaging. The primary roots were chosen and incubated on the slide in 30 µL 1xPBS containing 2 µg/mL propidium iodide (PI, Sigma-Aldrich, P4864) and incubated for 1 min before imaging. The imaging was performed on a Leica SP8 X confocal microscope using a 20X/0.75NA objective lens with water immersion, with illumination at 514 nm (Argon ion laser, 5% intensity). For the Venus channel, emissions were detected between wavelengths 518 nm and 550 nm by a cooled detector (Leica HyD SMD) in photon-counting mode. Propidium Iodide (PI) was detected simultaneously with FLC-Venus, by collecting emissions between wavelengths 610 nm and 680 nm in standard imaging mode. For each root, a z-step size of 0.95 μm was used, starting from the upper surface and then 22–40 z-slices depending on the thickness of the root. An image format of 1024 x 512 pixels with a 1.75X zoom factor was used.

The images were analysed to generate mean FLC-Venus intensities per cell using a processing pipeline, previously used in Yang et al., 2017 and available at https://github.com/JIC-Image-Analysis/root-3d-segmentation (Hartley and Olsson, 2016). The mean FLC-Venus intensity for a cell was determined by dividing the total FLC intensity by the total cell volume (in voxels). For the last time point, 14 days after a 10-week cold treatment, the pipeline was adapted to allow precise parameter tuning of the segmentation to correct for difficulties in imaging. This new pipeline is available at https://github.com/JIC-Image-Analysis/lhp1-analysis (Hartley, 2020). In this pipeline, segmentations, PI cell wall information and FLC-Venus fluorescence are combined to generate mean intensity measurements. Segmentations are loaded and reconstructed cells touching the border of the image are removed. The remaining cells are then rank ordered by size, and the top 300 cells (determined by examination of numbers of reconstructed cells in earlier time points) used to measure FLC-Venus intensity. Mean FLC-Venus intensity was then calculated for each of these remaining reconstructed cells by dividing total FLC intensity by the cell volume in voxels.

For each timepoint except the last at 14 days post cold, the background levels were estimated as at or below the 98-th percentile of mean cellular FLC-Venus intensity in roots at that timepoint where all cells exhibited silenced FLC-Venus. Cells above this threshold were defined as FLC-ON, and the proportion of FLC-ON cells was calculated for each time point (Figure 1C, Figure 2D, Figure 1—figure supplements 2 and 3). For the last time point, 14 days after 10 weeks of cold, all roots had some FLC-ON cells and therefore the 85-th percentile was used to conservatively define the FLC-ON cells. In Figure 1—figure supplement 3, we compare the histone-feedback only model to this data, concluding that even when using the 40-th percentile as a cut-off, it is still not possible to satisfactorily fit the model. Note that to ensure that the argon laser of Leica SP8X was stable after turning on, the five first roots imaged at each time-point were excluded. The statistics of the cells included in Figure 1—figure supplement 3 are summarised in Table 12.

FLC expression is low after 6 weeks of cold treatment in lhp1 and FRI^sf2 (Yang et al., 2017). However, in Berry et al., 2015, it was shown that even though there is FLC-Venus signal post 6 weeks vernalization, there is very little FLC-Venus protein present after 6 weeks cold in FRI^sf2 in an immunoblot. Here, we therefore treated the plants with 10 weeks of cold, to fully silence FLC. Some Venus intensity is nevertheless still observed 1, 2, and 4 days following return to warm after 10 weeks of cold, particularly at the 1-day time point. Furthermore, the Venus signal observed at these early time points, 1, 2, and 4 days after cold, is not located in the nucleus (Figure 1—figure supplement 4, top). This finding indicates that the early accumulation is a cold-induced phenomenon that might be Venus-specific. For these reasons, we therefore exclude these time points from our analysis, as they do not reflect transcriptional reactivation of FLC in the roots. In contrast, the Venus signal at later timepoints (7, 10, and 14 days after cold) is nuclear-localised (Figure 1—figure supplement 4, bottom), indicating that data from these time points do represent genuine FLC reactivation.

The images in Figure 1—figure supplement 4 were processed with a Gaussian blur (Fiji, 0.2 μm filter size) to the FLC-Venus channel. A sum projection over 9 z-slices (8.55 μm) was performed. For the cell wall stain, the central z-plane was extracted and overlaid on the FLC-Venus sum projection.

MN-ChIP experiments were performed as described in Yang et al., 2017 with modifications and omitting a strong sonication step. After washes with Honda buffer, the nuclei were washed once with Resuspension buffer (10 mM Tris-HCl pH=8; 5 mM NaCl; 2.5 mM CaCl2; 2 mM MgCl2; Protease inhibitor cocktail (Sigma, #5056489001)). Next, the pellet was weighted and resuspended in Resuspension buffer in a 1:5 ratio (mg pellet: µL buffer). 150 µL nuclei were added to commercial 10x Takara MNase buffer (Takara Bio, #2910A) and Micrococcal nuclease (Takara Bio, #2910A) to obtain 200 U/ml enzyme concentration and incubated at 37°C for 45 min. Suitable enzyme concentrations were determined experimentally to obtain chromatin fragmented to mononucleosomal resolution. The reaction was stopped by the addition of 25 mM EGTA (final concentration). Afterwards, the samples were sonicated twice using Bioraptor Pico (Diagenode, #B01060010) in mild conditions: 10 s ON/50 s OFF (low settings) and used for immunoprecipitation as described in Yang et al., 2017. The antibodies were: anti-H3 (Abcam, #ab1791) and anti-H3K27me3 (Millipore, #07–44).

The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.

We thank Rea Antoniou-Kourounioti, Geng-Jen Jang, Silvia Costa and Govind Menon for critical reading of the manuscript, and all members of the Howard and Dean groups for discussions.

1. Received: March 5, 2021
2. Accepted: August 3, 2021
3. Version of Record published: September 2, 2021 (version 1)
4. Version of Record updated: September 8, 2021 (version 2)

© 2021, Lövkvist et al.

This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.