DNA Methylation Heterogeneity Induced by Collaborations Between Enhancers

Abstract

During mammalian embryo development, reprogramming of DNA methylation plays important roles in the erasure of parental epigenetic memory and the establishment of naive pluripotent cells. Multiple enzymes that regulate the processes of methylation and demethylation work together to shape the pattern of genome-scale DNA methylation and guide the process of cell differentiation. Recent availability of methylome information from single-cell whole genome bisulfite sequencing (scBS-seq) provides an opportunity to study DNA methylation dynamics in the whole genome in individual cells, which reveal the heterogeneous methylation distributions of enhancers in embryo stem cells. In this study, we developed a computational model of enhancer methylation inheritance to study the dynamics of genome-scale DNA methylation reprogramming during exit from pluripotency. The model enables us to track genome-scale DNA methylation reprogramming at single-cell level during the embryo development process and reproduce the DNA methylation heterogeneity reported by scBS-seq. Model simulations show that DNA methylation heterogeneity is an intrinsic property driven by cell division along the development process, and the collaboration between neighboring enhancers is required for heterogeneous methylation. Our study suggests that the mechanism of genome-scale oscillation might not be necessary for the DNA methylation heterogeneity during exit from pluripotency.

1. Introduction

In mammalian development, reprogramming of DNA methylation (5-methylcytosine) patterns plays a crucial role in defining cell fate. Upon fertilization, DNA methylation marks represent an epigenetic barrier that restricts mammalian development and, hence, need to be restored and subsequently rebuilt with the commitment to particular cell fates (Seisenberger et al., 2013). The segregation of cell lineages gives rise to different somatic tissues associated with tissue-specific DNA methylation patterns (Styblo et al., 2000; Hon et al., 2013; Greenberg and Bourchis, 2019). The genome-wide DNA methylation reprogramming events coincide with the changes in concentrations of DNA methyltransferases (DNMTs) and the enzymes that initiate the removal of DNA methylation (ten-eleven-translocation family proteins, TETs) (Seisenberger et al., 2013; Greenberg and Bourchis, 2019). The recent maturation of single-cell sequencing technologies has enabled us to observe a variety of sequencing information at individual cell level, such as the genome, transcriptome, and epigenome (Stuart and Satija, 2019). It becomes a challenge issue in computational biology to develop single-cell based computational model that can help us to better understand the process of DNA methylation pattern formation, as well as cell fate decision during early embryo development.

The availability of methylome information from single-cell whole genome bisulfite sequencing (scBS-seq) provides an opportunity to study DNA methylation patterns in the whole genome in individual cells (Smallwood et al., 2014; Farlik et al., 2015; Rulands et al., 2018). A recent study applied scBS-seq to embryo stem cells (ESCs) cultured under naive [two chemical inhibitors (“2i”) of MEK1/2 and GSK3 $α$ /β] and primed (“serum”) conditions to explore DNA methylation dynamics in cells undergoing a biological transition (Rulands et al., 2018), primed ESCs had increased variance at several genomic annotations associated with active enhancer elements, including H3K4me1 and H3K27ac sites and low methylated regions. Analysis of scBS-seq data showed that individual primed ESCs have average DNA methylation levels varying between 17% and 86% at enhancers, while naive ESCs showed minimal cell-to-cell variability, and DNA methylation heterogeneity was resolved upon differentiation to embryoid bodies (Rulands et al., 2018).

In Rulands et al. (2018), it was proposed that the DNA methylation heterogeneity is associated with coherent, genome-scale oscillations in DNA methylation, and the amplitude is dependent on the CpG density. Moreover, a mathematical model of delay differential equation with autocatalytic de novo methylation was proposed to show that global oscillations may emerge from the biochemistry of methylation turnover due to Hopf bifurcation with increasing values of the time delay, and a Kuramoto model was applied to describe the global heterogeneous coupling of CpGs through DNMT3a/b binding. Nevertheless, genome-scale oscillations in DNA methylation is a very strong assumption, which may imply global oscillations in transcriptions of most genes, and is not supported by the experimental data of DNA methylation dynamics during transition from naive to primed pluripotency in vitro (Fig. 1F). Hence, we asked how the transitions of DNA methylation heterogeneity should be explained through a simple mechanism?

FIG. 1.

DNA methylation transition in embryo cells. (A) DNA methylation variance in naive and primed ESCs compared to the enhancer elements defined using published H3K4me1 ChIP-seq data (Creyghton et al., 2010) and the average methylation levels in individual cells. (B) Distributions of enhancer methylation levels in two cells from the naive condition. (C) Distributions of enhancer methylation levels in three cells from the primed condition. (D) DNA methylation variance of epiblast cells at E4.5, E5.5, and E6.5 of mice embryo and the average methylation levels in individual cells. (E) Heterogeneity (H) indexes of cells under different conditions. (F) Dynamics of average DNA methylation from 0 to 8 hours in the “2i release” experiments. Here the methylation levels were calculated from the average over the enhancers at ch1. All data were obtained from Rulands et al. (2018). ChIP-seq, chromatin immunoprecipitation sequencing; ESCs, embryo stem cells; H, heterogeneity index.

DNA methylation and chromatin dynamics have been modeled quantitatively in various genomic contexts of biological significance (Sneppen and Dodd, 2011; Haerter et al., 2013; Berry et al., 2017; Huang and Lei, 2017; Song et al., 2017). The collaboration between neighboring CpGs was highlighted in recent studies (Haerter et al., 2013; Song et al., 2017), which play essential roles in the formation of global patterns and the genome-scale transitions of DNA methylation. The collaboration may directly come from the binding of DNMT3a/b to neighboring CpGs (Rulands et al., 2018) or indirectly through the interaction with methylations in the histones H3K9 (Lehnertz et al., 2003) and H3K36 (Weinberg et al., 2019). In addition to the de novo methylation and demethylation, dilution and maintenance of methylated marks during DNA replication may also contribute to the cell-to-cell variance of DNA methylation.

In this study, motivated by the scBS-seq data of heterogeneous methylation at genomic annotations associated with active enhancer elements, we developed a model of DNA methylation, considering the stochastic dynamic of methylation levels of enhancers over cell divisions and the collaboration between neighboring enhancers, to investigate the transition of DNA methylation from naive to primed ESCs. The model focuses on the random inheritance of DNA methylations during cell cycling and can automatically reproduce the DNA methylation heterogeneity on enhancers during embryonic development. Our results suggest that the mechanism of genome-scale oscillation proposed by Rulands et al. (2018) may not be required for the observed heterogeneity during exit from pluripotency.

2. Results

2.1. Transition of DNA methylation patterns from naive to primed ESCs

We analyzed the scBS-seq data separately for ESCs cultured under naive (“2i”) and primed (“serum”) conditions (Rulands et al., 2018). Similar to the analysis in Rulands et al. (2018), taking published H3K4me1 chromatin immunoprecipitation sequencing (ChIP-seq) data from primed ESCs (Creyghton et al., 2010) as a definition of enhancer elements, the methylation levels of enhancers in primed ESCs increase comparing with naive ESCs (Fig. 1A). In this study, the methylation level of an enhancer is defined as the average level of all CpG sites contained in the enhancer. For each CpG site, we assigned a value 0 for unmethylated, 0.5 for half-methylated, and 1 for full methylated; hence, the methylation level of an enhancer takes value from the interval [0, 1] (or from 0% to 100% methylated). Moreover, we calculated the distribution of methylation levels of all enhancers in individual cells. Primed ESCs showed higher cell-to-cell variability at the distribution patterns than the naive ESCs (Fig. 1B, C). We also analyzed the parallel scM&T sequencing of in vivo epiblast cells at E4.5, E5.5, and E6.5 (Rulands et al., 2018), which showed an increase in the methylation level in enhancers from E4.5 to E5.5 (Fig. 1D). We note that there are a few cells that show low methylation levels at E5.5; however, all cells have high methylation at E6.5, which suggests a transition dynamics of methylation levels (Fig. 1D).

To quantify the heterogeneity of DNA methylations among different cells, we proposed a definition of heterogeneity index based on the methylation levels of enhancers in each cell. Assuming that there are n cells, and $p_{i} (1 \leq i \leq n)$ the distribution of all enhancer methylation levels of the i'th cell, we defined the heterogeneity index (H) as the average of Kullback–Leibler divergence between any two cells. Mathematically, the heterogeneity index is formulated as $H = \frac{1}{n (n - 1)} \sum_{i, j = 1}^{n} K L (p_{i} | | p_{j}),$ (1)

where $K L (p_{i} | | p_{j})$ means the Kullback–Leibler divergence between the distributions of enhancer methylation levels for the two cells i and j, $K L (p_{i} | | p_{j}) = \int_{0}^{1} p_{i} (x) log \frac{p_{i} (x)}{p_{j} (x)} d x .$

We calculated the heterogeneity index based on the above data from naive and primed ESCs and the cells at E4.5, E5.5, and E6.5 mice embryo. There are no significant changes in the heterogeneity of naive ESCs in comparing with the primed ESCs (Fig. 1E). For the mice embryo cells, the heterogeneity index increases from E4.5 to E5.5 and decreases from E6.6 to E6.5 (Fig. 1E).

To validate the genome-scale DNA methylation oscillations in Rulands et al. (2018), we analyzed the same data from an in vitro “2i release” model in which cells were transferred from naive 2i to primed serum culture conditions. The average methylation levels of enhancers from the first chromosome (ch1) were calculated to obtain the dynamics from 0 to 8 hours after 2i release. The results do not show significant oscillations in the methylation level (Fig. 1F). These findings suggest that the assumption of genome-scale DNA methylation oscillations might not be necessary to explain the transition of methylation and heterogeneity from naive to primed ESCs.

2.2. Stochastic dynamics of enhancer methylation levels

The dynamics of DNA methylation/demethylation have been modeled quantitatively with exquisite details at biochemistry of single CpG sites (Sontag et al., 2006; Haerter et al., 2013; Song et al., 2017; Lei et al., 2018). The methylation state of a single CpG site is often random due to the stochastic biochemical reactions. Nevertheless, the average methylation level of CpG sites associated with a genomic segment is more predictable. During embryo development, the most significant changes in DNA methylation occur during DNA replication when the 5-methylcytosine marks are diluted to two daughter strains and are restored through enzymes DNMT1 and nuclear protein 95 (NP95 or UHRF1). Correlating global DNA methylation with replication timing repli-seq data has shown that late-replicating regions did not have lower DNA methylation than early-replicating regions (Rulands et al., 2018). Thus, while we omit the detailed dynamics between DNA replications, we can represent the methylation level of an enhancer by the average methylation level at late-replication stage of each cell cycle.

2.2.1. Formulation

To consider the dynamics of enhancer methylation levels, assuming that there are N enhancers in a chromatin and letting $β_{i}^{t}$ the methylation level of the i'th enhancer at cycle $t (0 \leq β_{i}^{t} \leq 1)$ , we only need to formulate the dynamics of the states ${\vec{β}}^{t} = (β_{1}^{t}, β_{2}^{t}, \dots, β_{N}^{t})$

over cell cycles t. During cell cycling, the methylation states update as a consequence of the regulations through enzymes DNMT1, DNMT3a/b, and TETs, which lead to the following iteration

Hence, while we omit the biochemistry details, the stochastic dynamics of enhancer methylation levels can be formulated as an iteration ${\vec{β}}^{t + 1} = T ({\vec{β}}^{t})$ (3)

for each cell cycle. In this study, ${\vec{β}}^{t}$ represents the methylation state of a cell before cell division, and ${\vec{β}}^{t + 1}$ is the state of one daughter cell after cell division.

We note that the iteration Equation (3) is usually a random map. Given the state of mother cells, the methylation state of the daughter cell is a random valuable whose probability density is dependent on the state of the mother cell. Hence, to formulate the iteration map, we need to write down the conditional probability density function $P r o b ({\vec{β}}^{t + 1} = \vec{x} | {\vec{β}}^{t}),$ (4)

the probability of ${\vec{β}}^{t + 1}$ given the state of the mother cell ${\vec{β}}^{t}$ . While the probability of each enhancer, given the state of the mother cell, is independent to each other, we have $P r o b ({\vec{β}}^{t + 1} = \vec{x} | {\vec{β}}^{t}) = \prod_{i = 1}^{N} P r o b (β_{i}^{t + 1} = x_{i} | {\vec{β}}^{t}) .$ (5)

The probability density $P r o b (β_{i}^{t + 1} = x_{i} | {\vec{β}}^{t})$ is usually not known and may depend on the biochemistry details of methylation/demethylation. Nevertheless, it is possible to write down the phenomena formulation if we overlook the detail process. If there are m_i CpGs in the i'th enhancer, each CpG has a probability $p_{i}^{t}$ to be methylated (here the superscript t specified the dependence with the time t) and $1 - p_{i}^{t}$ to be unmethylated after cell division (here we omitted the state of half-methylation), the probability to have k_i methylated CpGs is given by a binomial distribution $C_{m_{i}}^{k_{i}} {(p_{i}^{t})}^{k_{i}} {(1 - p_{i}^{t})}^{m_{i} - k_{i}}$ with the parameters $p_{i}^{t}$ . To extend the probability to a more general situation, the binomial distribution can be generalized to a beta-binomial distribution through two shape parameters $a_{i}^{t}$ , $b_{i}^{t}$ . Moreover, while we are only interested at the probability of the methylation level $x_{i} = k_{i} ∕ m_{i}$ , we replaced the beta-binomial distribution with the beta distribution and hence $P r o b (β_{i}^{t + 1} = x_{i} | {\vec{β}}^{t}) = \frac{x^{a_{i}^{t} - 1} {(1 - x_{i})}^{b_{i}^{t} - 1}}{B (a_{i}^{t}, b_{i}^{t})}, B (a, b) = \frac{Γ (a) Γ (b)}{Γ (a + b)},$ (6)

where $Γ (\cdot)$ is the gamma function, and $a_{i}^{t}$ and $b_{i}^{t}$ are shape parameters depending on ${\vec{β}}^{t}$ . Equations (5) and (6) together define the conditional probability $P r o f ({\vec{β}}^{t + 1} = \vec{x} | \vec{β})$ . The beta distribution is one of the few common “named” distributions that give probability 1 to a finite interval. As the shape parameters a and b vary, the beta distribution can take different shapes, including strictly decreasing ( $a \leq 1, b > 1$ ), strictly increasing ( $a > 1, b \leq 1$ ), U-shaped ( $a < 1, b < 1$ ), or unimodal ( $a > 1, b > 1$ ). Thus, the dependences of parameters $a_{i}^{t}, b_{i}^{t}$ on the state ${\vec{β}}^{t}$ are important to define the shape of the distribution function.

Now, to define the dependences $a_{i}^{t} ({\vec{β}}^{t})$ and $b_{i}^{t} ({\vec{β}}^{t})$ , assuming that there are functions $ϕ_{i}^{t} ({\vec{β}}^{t})$ and $η_{i}^{t}$ so that the average of $β_{i}^{t + 1}$ , given ${\vec{β}}^{t}$ , is ${⟨ β_{i}^{t + 1} ⟩|}_{{\vec{β}}^{t}} = ϕ_{i}^{t} ({\vec{β}}^{t}),$

and the variance ${V a r (β_{i}^{t + 1})|}_{{\vec{β}}^{t}} = \frac{1}{1 + η_{i}^{t}} (1 - ϕ_{i}^{t} ({\vec{β}}^{t})) ϕ_{i}^{t} ({\vec{β}}^{t}),$

than¹ $a_{i} = η_{i}^{t} ϕ_{i}^{t} ({\vec{β}}^{t}), b_{i}^{t} = η_{i}^{t} (1 - ϕ_{i}^{t} ({\vec{β}}^{t})) .$ (7)

Hence, we only need to identify the functions $ϕ_{i} ({\vec{β}}^{t})$ and $η_{i}^{t}$ , and the parameters $a_{i}^{t}$ and $b_{i}^{t}$ can be defined accordingly. Specifically, we take $η_{i}^{t}$ as the number of CpGs in the i'th enhancer, that is, $η_{i}^{t} = m_{i} .$ (8)

This simple assumption means the inverse proportion of the variance of enhancer methylation levels with the number of CpGs.

To define the function $ϕ_{i}^{t}$ , we assumed that the methylation level of daughter cell depends on that of the mother cell through three components: basal methylation level, autocatalytic effect, and collaboration between enhancers. Hence, the function $ϕ_{i}^{t}$ was formulated as

where $H_{[0, 1]} (z) = \{\begin{matrix} 0, & z < 0, \\ z, & 0 \leq z \leq 1 \\ 1, & z > 1 . \end{matrix}$ (10)

Here $μ_{0}, μ_{1}, n, v, α, L$ are parameters, with $μ_{0}$ the basal level, $μ_{1}$ the coefficient for autocatalysis, n the Hill coefficient, v the parameter for the autocatalytic efficiency, and $α$ the coefficient for long distance collaboration, and L defines the range of collaboration between enhancers. The autocatalytic efficiency υ usually depends on the activity of de novo methylation/demethylation regulated by DNMT3a/b and TETs and, hence, can change with time t during embryo development. Here, we always have $0 \leq ϕ_{i}^{t} \leq 1$ due to the function $H_{[0, 1]} (z)$ .

The third term in Equation (9) shows the collaboration between neighboring enhancers. Here, we assumed the coherent collective behaviors of DNA methylation/demethylation when the average coupling through enzyme binding is sufficiently strong so that the nearby enhancers tend to the same trends of either methylation or demethylation. The similar mechanism was introduced previously to reproduce the long-distance correlation of DNA methylation between CpGs (Song et al., 2017). Here the collaboration effect is limited by cooperative range L and the coefficient $α$ .

2.2.2. Numerical scheme and parameters

To model the methylation dynamics following developmental process with the above iteration Equation (2), we first initialized the methylation level of N enhancers ${\vec{β}}^{0} = (β_{1}^{0}, β_{2}^{0}, \dots, β_{N}^{0})$ (here $t = 0$ ). Here, each enhancer i associates with an integer m_i for the number of CpGs. Next, at each step, we calculated the shape parameters $a_{i}^{t}$ and $b_{i}^{t}$ using Equations (7)–(9). Finally, for each enhancer i generated a random number according to the beta distribution Equation (6) and set $t = t + 1$ . We repeated the above scheme to generate the dynamics of methylations in each enhancer following multiple cell cycles.

In model simulations, we took $μ_{0} = 0.01, μ_{1} = 0.9, n = 3, α = 0.5, L = 50$ and varied the autocatalysis ability v to represent different culture conditions (v = 0.1, 0.04, 0.01 for situations of high, mediate, and low level DNA methylation, respectively). Moreover, to mimic enhancers in a chromosome, we performed simulations with $N = 1000$ enhancers; the CpG numbers for each enhancer were taken following an exponential distribution $P r o b (m_{i} = m) = \frac{e^{(- m ∕ 3)}}{1.6}$ (11)

in accordance to the statistics from mouse genome (Fig. 2A). In simulations, we can refer the distribution of methylation levels from a naive cell as the initial condition (Fig. 2B).

FIG. 2.

Number of CpG sites and the initial methylation level. (A) Distribution of number of CpG sites in enhancers and the fitting curve [Eq. (11)]. (B) Distributions of methylation at a naive ESC.

2.3. Transition of DNA methylation heterogeneity during the development process

To verify the proposed model, we varied the parameter v (v = 0.1, 0.04, 0.01) to mimic different conditions. For each value v, we initialized a cell with an initial state of a naive cell and ran the model simulation for 15 cell cycles to mimic the transition from naive to primed condition. The simulated distribution of enhancer methylation levels at each cell cycle was calculated and is shown by Figure 3A–C. The enhancer methylation distributions depend on the parameter v: when $v = 0.1$ , the enhancers were homogeneous with low level methylation; when $v = 0.04$ , the cells showed obvious methylation heterogeneity, the methylation levels transfer from low to high along with cell cycling; and when $v = 0.01$ , most enhancers change to high level methylations in a few cycles. When $v = 0.04$ , we have low, mediate, or high methylation levels in the enhancers at different cycles (Fig. 3D). These results are in agreement with experimental observations of ESCs cultured under serum condition, and hence, the model can be used to mimic the DNA methylation heterogeneity of ESCs during exit from pluripotency.

FIG. 3.

DNA methylation heterogeneity from model simulation. (A–C) Distributions (violin plots) of methylation levels in enhancers obtained from model simulation with $v = 0.1$ (A), $v = 0.04$ (B), and $v = 0.01$ (C). (D) Histogram methylation levels in enhancers for cycles 0, 5, and 14. (E) Transition of DNA methylation from low to high level. Each dot represents the average methylation level in a cell. (F) Evolution of the heterogeneity index. The heterogeneity index was calculated from the cells shown by (E). Here $v = 0.04$ in (E, F).

To further examine the transition dynamics of methylation heterogeneity in ESCs, we set $v = 0.04$ and initialized a population of cells according to the methylation distribution at Figure 2B. Next, we performed the simulation scheme to mimic a development process of 48 hours. In simulations, each cell divides with a probability of 0.3 per hour and collected a subpopulation of cells to calculate the average methylation level and the heterogeneity index every 3 hours. Simulations showed that the average methylation level increased from 0 to 48 hours, but showed obvious diversity from 12 to 30 hours (Fig. 3E). The heterogeneity index that showed nonmonotonous with the development process first increase from 0 to 27 hours to reach a high level of 0.12 and then decrease to a low level heterogeneity at 48 hours (Fig. 3F). In embryo development, DNA methylation heterogeneity in ESCs increased from naive to primed, and the heterogeneity was resolved upon differentiation to embryoid bodies (Rulands et al., 2018). These results showed similar dynamics in both experiments and our model simulations.

2.4. Collaboration and methylation heterogeneity

The proposed model includes collaboration between neighboring enhancers so that there are coherent collective behaviors during enzyme binding. To investigate the effects of neighboring collaboration, we varied the parameters $α$ and L and examined the changes in methylation heterogeneity. Here, the parameter $α$ measures the collective strength, and L gives the regions of neighboring collaboration. Simulations showed that the heterogeneity increased with $α$ for different values of L and increased with L for a large value $α$ (Fig. 4A).

FIG. 4.

Collaboration and methylation heterogeneity. (A) Heterogeneity with different values of the parameter α and L. (B–E) Transition dynamics of the average methylation level with varied parameter values of α and L.

To further examine how collaborations may affect the DNA methylation dynamics, we varied the parameters $α$ and L and calculated the evolution dynamics over a period of 48 hours. When $α = 0.5$ and $L = 10$ , the cells showed as highly heterogeneous during the intermediate transition region (Fig. 4B). When either L or $α$ decreases, the cells showed less heterogeneity (Fig. 4C, D). In particular, when $α = 0$ , which represents the situation without collaboration, all cells showed similar DNA methylation dynamics during simulation from 0 to 48 hours, with the average methylation increasing from a low level to an intermediate level of 50% (Fig. 4E). These results suggest that the collaboration between neighboring enhancers is required to produce DNA methylation heterogeneity.

3. Discussion

Reprogramming of DNA methylation plays important roles in mammalian early embryo development. ESCs show heterogeneous methylation distributions under primed conditions. To understand the mechanism of methylation heterogeneity, previous studies suggest a mechanism of genome-scale oscillations in DNA methylation. Nevertheless, experiment data did not support the assumption of genome-scale oscillations. In this study, we proposed a computational model for the stochastic transitions of enhancer methylations during cell cycling. The model combines random distribution of methylation marks in DNA replication, autocatalysis of DNA methylation due to the binding of DNMT3a/b, and the collaboration between neighboring enhancers during the reconstruction of methylation marks. The proposed model can nicely explain the transition of methylation level and heterogeneous methylation distributions. Model simulations showed that the proper values of the autocatalysis are important for the heterogeneity between different cells, and increasing the collaboration between neighboring enhancers can promote the heterogeneity. Our model suggests that methylation heterogeneity is a natural consequence of stochastic transition of DNA methylation between cell cycles and the collaboration between CpGs; however, the assumption of genome-scale oscillations might not be necessary for the observed heterogeneous methylation distributions.

The proposed model mainly considers the dynamics of enhancer methylation levels in a cell and omits the biochemical reactions involved in methylation or demethylation, which may be regulated by various enzymes. In the model, we assumed a beta distribution that connects the methylation level in daughter cells with those of the mother cell. The beta distribution can take different forms based on the shape parameters that are defined by the state of mother cells. Thus, the proposed model framework mainly focuses at the general effect of methylation state transition between cell cycles, while omitting the details of biochemical reactions. In contrast, despite the complex biochemical reactions, they mainly affect the methylation levels through the based methylation/demethylation processes and the collaboration between CpGs and, hence, may end up to the function $ϕ$ in the model. Hence, the proposed model provides a general framework to sum up different level biochemical reaction model for DNA methylation.

Footnotes

Author Disclosure Statement

The authors declare they have no conflicting financial interests.

Funding Information

This work is supported by the National Natural Science Foundation of China (11831015 and 11872084).

References

Berry

, Dean

, and Howard

2017. Slow chromatin dynamics allow polycomb target genes to filter fluctuations in transcription factor activity. Cell Syst. 4, 445–457.e8.

Creyghton

M.P.

, Cheng

A.W.

, Welstead

G.G.

, et al. 2010. Histone H3K27ac separates active from poised enhancers and predicts developmental state. Proc. Natl. Acad. Sci. U S A. 107, 21931–21936.

Farlik

, Sheffield

N.C.

, Nuzzo

, et al. 2015. Single-cell DNA methylome sequencing and bioinformatic inference of epigenomic cell-state dynamics. Cell Rep. 10, 1386–1397.

Greenberg

M.V.C.

, and Bourchis

2019. The diverse roles of DNA methylation in mammalian development and disease. Nat. Rev. Mol. Cell Biol. 328, 916–18.

Haerter

J.O.

, Lövkvist

, Dodd

I.B.

, et al. 2013. Collaboration between CpG sites is needed for stable somatic inheritance of DNA methylation states. Nucleic Acids Res. 42, 2235–2244.

Hon

G.C.

, Rajagopal

, Shen

, et al. 2013. Epigenetic memory at embryonic enhancers identified in DNA methylation maps from adult mouse tissues. Nat. Genet. 45, 1198–1206.

Huang

, and Lei

2017. Dynamics of gene expression with positive feedback to histone modifications at bivalent domains. Int. J. Mod. Phys. B, 4, 1850075.

Lehnertz

, Ueda

, Derijck

A.A.H

.A., et al. 2003. Suv39h-mediated histone H3 lysine 9 methylation directs DNA methylation to major satellite repeats at pericentric heterochromatin. Curr. Biol. 13, 1192–1200.

Lei

, Nie

, and Chen

2018. A single-cell epigenetic model for paternal psychological stress-induced transgenerational reprogramming in offspring. Biol. Reprod. 98, 846–855.

10.

Rulands

, Lee

H.J.

, Clark

S.J.

, et al. 2018. Genome-scale oscillations in DNA methylation during exit from pluripotency. Cell Syst. 7, 63–76.

11.

Seisenberger

, Peat

J.R.

, Hore

T.A.

, et al. 2013. Reprogramming DNA methylation in the mammalian life cycle: Building and breaking epigenetic barriers. Phil. Trans. R. Soc. B Biol. Sci. 368, 20110330.

12.

Smallwood

S.A.

, Lee

H.J.

, Angermueller, et al. 2014. Single-cell genome-wide bisulfite sequencing for assessing epigenetic heterogeneity. Nat. Methods, 11, 817–820.

13.

Sneppen

, and Dodd

I.B.

2011. A simple histone code opens many paths to epigenetics. PLoS Comput. Biol. 8, e1002643.

14.

Song

, Ren

, and Lei

2017. Collaborations between CpG sites in DNA methylation. Int. J. Mod. Phys. B, 31, 1750243.

15.

Sontag

L.B.

, Lorincz

M.C.

, and Georg Luebeck

2006. Dynamics, stability and inheritance of somatic DNA methylation imprints. J. Theor. Biol. 242, 890–899.

16.

Stuart

, and Satija

2019. Integrative single-cell analysis. Nat. Rev. Genet. 20, 257–272.

17.

Styblo

, Del Razo

L.M.

, Vega

, et al. 2000. Comparative toxicity of trivalent and pentavalent inorganic and methylated arsenicals in rat and human cells. Arch. Toxicol. 74, 289–299.

18.

Weinberg

D.N.

, Papillon-Cavanagh

, Chen

, et al. 2019. The histone mark H3K36me2 recruits DNMT3A and shapes the intergenic DNA methylation landscape. Nature, 69, 915–916.