Improved inference of population histories by integrating genomic and epigenomic data

Abstract
eLife assessment
Introduction
Results
Discussion
Materials and methods
Appendix 1
Appendix 2
Data availability
References
Article and author information
Metrics

Abstract

With the availability of high-quality full genome polymorphism (SNPs) data, it becomes feasible to study the past demographic and selective history of populations in exquisite detail. However, such inferences still suffer from a lack of statistical resolution for recent, for example bottlenecks, events, and/or for populations with small nucleotide diversity. Additional heritable (epi)genetic markers, such as indels, transposable elements, microsatellites, or cytosine methylation, may provide further, yet untapped, information on the recent past population history. We extend the Sequential Markovian Coalescent (SMC) framework to jointly use SNPs and other hyper-mutable markers. We are able to (1) improve the accuracy of demographic inference in recent times, (2) uncover past demographic events hidden to SNP-based inference methods, and (3) infer the hyper-mutable marker mutation rates under a finite site model. As a proof of principle, we focus on demographic inference in Arabidopsis thaliana using DNA methylation diversity data from 10 European natural accessions. We demonstrate that segregating single methylated polymorphisms (SMPs) satisfy the modeling assumptions of the SMC framework, while differentially methylated regions (DMRs) are not suitable as their length exceeds that of the genomic distance between two recombination events. Combining SNPs and SMPs while accounting for site- and region-level epimutation processes, we provide new estimates of the glacial age bottleneck and post-glacial population expansion of the European A. thaliana population. Our SMC framework readily accounts for a wide range of heritable genomic markers, thus paving the way for next-generation inference of evolutionary history by combining information from several genetic and epigenetic markers.

eLife assessment

This important study extends existing sequentially Markovian coalescent approaches to include the combined use of SNPs and hypervariable loci such as epimutations. This is an intriguing addition to infer population size history in the recent past, and the authors provide solid validation of their methods via simulation and analysis of empirical data in Arabidopsis thaliana. Given the increasing availability of such data, this work is a timely contribution and represents a foundation for further developments to explore when and where these methods will be best used.

https://doi.org/10.7554/eLife.89470.4.sa0

Introduction

A central goal in population genetics is to reconstruct the evolutionary history of populations from patterns of genetic variation observed in the present. Relevant aspects of these histories include past demographic changes as well as signatures of selection. Inference methods based on deep learning (DL, Korfmann et al., 2023), approximate Bayesian computation (ABC, Boitard et al., 2016), or sequential Markovian coalescent (SMC, Li and Durbin, 2011; Schiffels and Durbin, 2014) aim to infer this information directly from full genome sequencing data, which is becoming rapidly available for many (non-model) species due to decreasing costs. The SMC, in particular, offers an elegant theoretical framework that builds on the classical Wright-Fisher and the backward-in-time Kingman coalescent stochastic models (e.g. Kingman, 1982; Charlesworth and Charlesworth, 2010; Wakeley, 2008). Both models conceptualize Mendelian inheritance as generating the genealogy of a population (or a sample), that is, the unique history of a fragment of DNA passing from parents to offspring. When this genealogy includes the effect of recombination, it is called the ancestral recombination graph (ARG, Hudson, 1983; Wiuf and Hein, 1999).

Under the Kingmann coalescent model, the true genealogy of a population (or sample) is defined by its topology and branch length, and contains the information on past demographic changes and life history traits (Nordborg, 2000; Sellinger et al., 2020; Strütt et al., 2023; Tellier et al., 2011) as well as selective events (Charlesworth and Charlesworth, 2010; Wakeley, 2008). The genealogical and the mutational processes of any heritable marker can, therefore, be disentangled, and the frequency of any given marker state is given by the shape of the genealogy in time (see Figure 1A). A central assumption about heritable genomic markers is that they are generated by two homogeneous Poisson mutation processes along the genome as well as through time. This entails that mutations in different genealogies are independent due to the effect of recombination (Wiuf and Hein, 1999; McVean and Cardin, 2005), and that there are no time periods with a large excess, or a severe lack, of mutations along a genealogy (mutations are independently distributed in time within a DNA fragment). In other words, the frequencies of polymorphisms at DNA markers observed across a sample of sequences are constrained by, as well as inform on, the underlying genealogy at this locus (Figure 1A). To clarify these assumptions, we present a schematic representation of marker 1 (yellow in Figure 1), which fulfils both homogeneous Poisson processes in time and along the genome. We also present cases applicable to a second genomic marker 2 that violates the model assumptions, namely by not being heritable (Figure 1B) or not following a non-homogeneous Poisson process in the genome (Figure 1C) or in time (Figure 1D).

Figure 1 with 1 supplement see all

Download asset Open asset

Schematic distribution of two markers along the genealogy and four genomes.

(A) Schematic distribution of marker 1 (yellow star) and marker 2 (green star) along the genealogies in a sample of four genomes, both following a homogeneous Poisson process. (B) The green marker 2 is not heritable so its distribution is independent of the genealogy. (C) The green marker 2 is spatially structured along the genome, violating the distribution of the Poisson process along the genome and conflicting with the genealogy. (D) The green marker 2 does not follow the Poisson process through time, for example burst of mutations at a specific time point represented by given branches of the genealogies in green. The yellow marker 1 has an identical Poisson process along the genome and the genealogy in all four panels, and for readability, marker 2 exhibits light and dark green states.

Despite the power of the SMC, well-known model violations such as variation in recombination and mutation rates along the genome (Barroso et al., 2019; Barroso and Dutheil, 2023) or pervasive selection (Schraiber and Akey, 2015; Johri et al., 2021; Johri et al., 2020) can compromise the accuracy of demographic and selective inference (Gattepaille et al., 2013; Sellinger et al., 2021). There are two other important issues that have received less attention in the literature. The first issue occurs when the population recombination rate (ρ) is higher than the population mutation rate (θ). In such cases, inferences can be biased if not erroneous (Terhorst et al., 2017; Sellinger et al., 2021; Sellinger et al., 2020), because several recombination events cannot be inferred due to the lack of single-nucleotide polymorphisms (SNPs for point mutations). This problem affects many species, though interestingly not humans which have a ratio $ρ / θ \approx 1$ . A second issue occurs when the mutational process along the genealogy is too slow to be informative about sudden and strong variation in population size (i.e. population bottlenecks), such as during colonization events of novel habitats. The typical low mutation rate of 10⁻⁹ up to 10⁻⁸ (per base, per generation) found in most species, therefore, places strong limitations on SMC analysis of recent bottleneck events (up to ca. 10⁴ generations ago) when inference is based solely on SNP data. Indeed, bottlenecks are often either not found, or when inferred, their timing and magnitude are not well estimated (inferred smoother than in reality, Johri et al., 2021; Sellinger et al., 2021), even when a large number of samples is used. A typical example is the large uncertainty of the timing and magnitude of the population size bottleneck during the last glacial maximum (LGM) and post-LGM expansion in Arabidopsis thaliana European populations based on several studies using different accessions and SMC inference methods (Alonso-Blanco et al., 2016; Durvasula et al., 2017).

Nonetheless, current SMC, DL, or ABC inference methods making use of full genome sequence data rely almost exclusively on SNPs for inference (Schiffels and Durbin, 2014; Terhorst et al., 2017; Sellinger et al., 2020; Boitard et al., 2016; Korfmann et al., 2024). There are both practical and theoretical reasons for using SNPs: They are easily detectable from short-read re-sequencing data and their mutational process is well approximated by the infinite site model (Charlesworth and Charlesworth, 2010; Wakeley, 2008), simplifying the inference of the underlying genealogy. However, other heritable genomic markers exist whose mutation rates can be several orders of magnitude higher than that of SNPs and could thus be more informative about recent demographic events. These include microsatellites, insertions, deletions, and transposable elements (TEs). Although those heritable markers are not necessarily neutral (such as TEs, which are likely to be under weak purifying selection), they contain information on the evolutionary history of the population. Current technological limitations still impede the easy detection and estimation of allele frequencies for many of these markers (Yang et al., 2018; Ou et al., 2019; Wang, 2018). For example, identifying insertion/excision variation of transposable elements or copy number variation of microsatellites requires a high-quality reference genome and ideally long-read sequencing approaches (Ou et al., 2019). In addition to these genomic markers, DNA cytosine methylation is emerging as a potentially useful epigenetic marker for phylogenetic inference in plants (Yao et al., 2021; Yao et al., 2023). Stochastic gains and losses of DNA methylation at CG dinucleotides, in particular, arise at a rate of ca. 10⁻⁴ up to 10⁻³ per site per generation (that is 4–5 orders of magnitude faster than DNA point mutations, van der Graaf et al., 2015), and can be inherited across generations (Pisupati et al., 2023; Weigel and Colot, 2012). These so-called spontaneous epimutations are likely neutral at the genome-wide scale (Vidalis et al., 2016; Johannes and Schmitz, 2019b, but see Muyle et al., 2021; Pisupati et al., 2023), and can be easily detected from bisulphite converted short read sequencing data (Lister et al., 2008; Schmitz et al., 2013). Recent work suggests that CG methylation data can be used as a molecular clock for timing divergence between pairs of lineages over timescales ranging from years to decades (Yao et al., 2023).

However, the theoretical integration of the above-mentioned (epi)genomic markers into a population genomics and SMC inference framework is not trivial. Because of the high mutation rate, the mutational process at these (hyper-mutable) markers is reversible and more consistent with a finite site, rather than infinite site, model, which can result in extensive homoplasy (as known for microsatellite markers, Estoup et al., 2002). Indeed, classic expectations of population genetics diversity statistics, mostly built for SNPs, need to be revised for these hyper-mutable markers (Charlesworth and Jain, 2014; Wang and Fan, 2014). Here, we develop the theoretical and methodological inference framework named SMCtheo to include additional (potentially hyper-mutable) markers in the SMC. We showcase our model using extensive simulations as well as application to published DNA cytosine methylation data (in genic regions) from local populations of A. thaliana (Schmitz et al., 2013; Vidalis et al., 2016). We demonstrate that integration of hyper-mutable genomic markers into SMC models significantly improves the inference accuracy of past variation of population size or can even uncover demographic events not uncovered using SNPs alone. Our proof-of-principle approach opens up novel avenues for studying population genetic processes over time scales that have been largely inaccessible using traditional SNP-based approaches. This may prove particularly useful when exploring recent demographic changes of endangered species as a way to assess their potential for extinction in the context of biodiversity loss and global change.

Results

Theoretical results with two markers underlying the SMC computations

We study polymorphic sites across genomes of several sampled individuals which exhibit several possible markers (DNA nucleotides, methylation, TEs, indels, microsatellites,…). We define any marker by (1) its maximum number of possible states ( $n b_{s}$ ), for example nucleotide sites have four states (A, T, C and G) while a methylation site has two states (methylated or unmethylated), and (2) its mutation rate μ, that is the rate at which the state of a marker changes into another state per position and per generation Anzai et al., 2003; for simplicity, we assume equal mutation rates between all bases, known as the Jukes-Cantor model. More specifically, we are interested in two rates: the DNA mutation rate for changes in DNA nucleotides and the epimutation rate for changes in methylation state. Furthermore, we assume that at each position on the genome, only one type of marker can occur and be observed. We obtain as a first theoretical result the probability for a given site in the genome to be identical ( $P (i d)$ ) or segregating ( $P (s e g)$ ) (i.e. polymorphic) in a sample of size two ( $n = 2$ , two sampled chromosomes are compared):

\begin{array}{ll} P (i d, n = 2) = \frac{1}{n b_{s}} + \frac{(n b_{s} - 1)}{n b_{s}} e^{- 2 μ t_{M} \frac{(n b_{s})}{(n b_{s} - 1)}} \\ P (s e g, n = 2) = \frac{(n b_{s} - 1)}{n b_{s}} - \frac{(n b_{s} - 1)}{n b_{s}} e^{- 2 μ t_{M} \frac{(n b_{s})}{(n b_{s} - 1)}} \end{array}

This probability is a function of the time to the most recent common ancestor (TMRCA in text and $t_{M}$ in Equation 1, details in Appendix 1 and 2 B). The probability for a mutation to occur for a given marker increases with an increased TMRCA (Charlesworth and Charlesworth, 2010; Wakeley, 2008), but under high mutation rates (and high effective population size) the marker may not be polymorphic in the sample as mutations may be reversed (so-called homoplasy, Estoup et al., 2002; Charlesworth and Jain, 2014). In Figure 1—figure supplement 1, we illustrate these properties by computing the probability in Equation 1 for different mutation rates. The inference of recent demographic events and bottlenecks relies on the presence of polymorphic sites to detect recent coalescent events (TMRCA), and should be improved by using markers with high (or fast) mutation rate (e.g. hyper mutable).

In the following, we simulate data under different demographic scenarios using the sequence simulator program msprime (Baumdicker et al., 2022; Kelleher et al., 2016), which generates the ARG of $n$ sampled diploid individuals (set to $n = 5$ throughout this study, leading to 10 haploid genomes). This ARG contains the genealogy of a given sample at each position of the simulated chromosomes. We then process the ARG to create DNA sequences according to the model parameters and the type of marker considered. We first assume a set of genomic markers obtained for a sample size $n$ and mutating according to a homogeneous Poisson process along the genome and in time (along the genealogy) as in Figure 1A. To simulate the sequence data, we define the number of marker types (any number between 1 and the sequence length) and the proportion of sites of each marker type in the sequence. Each marker is characterized by both parameters $n b_{s}$ and μ. For simplicity, we simulate sequences with two markers but note that the method can be easily extended to additional markers. Marker 1 represents 98% of the sequence and has a per-site mutation rate $μ_{1} = 10^{- 8}$ mimicking nucleotide SNP markers under an infinite site model (thus considered as bi-allelic at a given DNA site, Yang, 1996). By contrast, marker 2 composes the complementary 2% of the sequence length, with a per-site mutation rate of $μ_{2} = 10^{- 4}$ per generation between two possible states. Marker 2 is thus hyper-mutable compared to marker 1 and mimics methylation/epimutation sites. Note that mutation events at Markers 1 and 2 are simulated under a finite site model.

We use different SMC-based methods throughout this study. These methods include: (1) MSMC2 used as a reference method (Malaspinas et al., 2016), (2) SMCtheo is an extension of the PSMC’ (Li and Durbin, 2011; Schiffels and Durbin, 2014) accounting for any number of heritable theoretical markers, and (3) eSMC2 which is equivalent to SMCtheo but accounting only for SNPs markers (Sellinger et al., 2021) (to avoid any bias in implementation differences between SMCtheo and MSMC2). All methods are hidden Markov models (HMM) derived from the pairwise sequentially Markovian coalescent (PSMC’) (Schiffels and Durbin, 2014) and assume neutral evolution and a panmictic population. The hidden states of these methods are the coalescence time of a sample of size two at a position on the sequence. From the distribution of the hidden states along the genome, all methods can infer population size variation through time as well as the recombination rate (Schiffels and Durbin, 2014; Malaspinas et al., 2016; Sellinger et al., 2021).

The inclusions of hyper-mutable genomic markers improves demographic inference

We assume that the mutation rate of marker 1 is $μ_{1} = 10^{- 8}$ per generation per bp. We use this information to estimate the mutation rate of marker 2, which we vary from $μ_{2} = 10^{- 8}$ to $μ_{2} = 10^{- 2}$ per generation per bp. The estimation results based on simulated data under a constant population size of $N = 10, 000$ are displayed in Table 1. We find that our approach is capable of inferring $μ_{2}$ with high accuracy for rates up to $μ_{2} = 10^{- 4}$ . However, when the mutation rate $μ_{2}$ is 10^-2, our approach underestimates it by a factor three, suggesting the existence of an accuracy limit. To demonstrate that information can be gained by integrating marker 2 (with $μ_{2} = 10^{- 4}$ ), we compared the ability of several inference methods to recover a recent bottleneck (Figure 2A). All methods correctly infer the amplitude of population size variation. When accounting only for marker 1 (with $μ_{1} = 10^{- 8}$ ), MSMC2 and eSMC2 fail to infer accurately the sudden variation of population size. However, with the inclusion of hyper-mutable marker 2, our SMCtheo approach correctly infers the rapid change of population size of the bottleneck (Figure 2A, green). It is encouraging that an accurate estimation of the demography is obtained, even when the mutation rate of marker 2 is unknown (Figure 2A, blue).

Figure 2 with 2 supplements see all

Download asset Open asset

Performance of SMC approaches using different markers.

Estimated demographic history of a bottleneck (black line) by SMC approaches using two genomic markers. In orange and red are the estimates by MSMC2 and eSMC2 based on only marker 1. Estimates from SMCtheo integrating both markers are in green (with known $μ_{2}$ ) and in blue with unknown $μ_{2}$ . The demographic scenarios are (A) 10-fold recent bottleneck with an ancestral population size $N = 10, 000$ , (B) 10-fold recent bottleneck with an ancestral population size $N = 1, 000$ , (C) 10-fold bottleneck with an ancestral population size $N = 10, 000$ , and (D) a very severe (1000 fold) and very recent bottleneck with incomplete size recovery. In A, B, and D, we assume $r / μ_{1} = 1$ (with $r = μ_{1} = 10^{- 8}$ , $μ_{2} = 10^{- 4}$ per generation per bp) and in C, $r / μ_{1} = 10$ (with $r = 10^{- 7}$ , $μ_{1} = 10^{- 8}$ , and $μ_{2} = 10^{- 4}$ per generation per bp). In all cases (**A, B, C and D**) 10 sequences (5 diploid indivudals) of 100 Mb were used as input.

Table 1

Average estimated mutation rate of the second theoretical genomic marker.

Average estimated values of the mutation rate of marker 2 ( $μ_{2}$ ), knowing that of marker 1. We use 10 sequences (5 diploid individuals) of 100 Mb ( $r = μ_{1} = 10^{- 8}$ per generation per bp) under a constant population size fixed at $N = 10, 000$ . The coefficient of variation over 10 repetitions is indicated in parentheses.

True $μ_{2}$ value	Estimated value of $μ_{2}$
10⁻⁸	9.9×10⁻⁹ (0.02)
10⁻⁶	1.0×10⁻⁶ (0.008)
10⁻⁴	1.4×10⁻⁴ (0.01)
10⁻²	3.05×10⁻³ (0.41)

Furthermore, some species or populations might feature small effective population sizes (ca. $N = 1, 000$ ), potentially resulting in reduced genomic diversity. In such cases, the inclusion of hyper-mutable markers should also improve demographic inference. We present the results of such a scenario in Figure 2A, where the population size was divided by a factor 10 compared to the previous scenario in Figure 2A. We find that in the absence of the hyper-mutable marker 2, no approach can correctly infer the variation of population size. From the shape of the inferred demography, methods using only marker 1 do not suggest the existence of a bottleneck followed by recovery (the ‘U-shaped’ demographic scenario is not apparent with the orange and red lines, Figure 2B). Yet, when integrating both markers, the population size can be recovered, even if the mutation rate of marker 2 is not a priori known. In both Figure 2A and B, we assume that marker 2 occurs at a frequency of 2% in the genome. This percentage may be unrealistically high depending on the marker and the species. To test the impact of reducing marker 2 frequency, we repeat the simulations shown in Figure 2A but set its frequency to as low as 0.1% (a 20-fold reduction). We find that the inclusion of the hyper-mutable marker 2 continues to improve inference accuracy in very recent times, albeit less pronounced than in Figure 2A (see Figure 2—figure supplement 1). This suggests that a very small proportion of hyper-mutable genomic sites is sufficient to significantly improve the accuracy of inferences.

All full genome inference methods, especially SMC approaches, display lower accuracy when the population recombination rate ( $ρ = 4 N r$ ) is larger than the population mutation rate of marker 1 ( $θ_{1} = 4 N μ_{1}$ ). We simulate sequence data under a bottleneck scenario slightly more ancient than in Figure 2A and assume that $ρ / θ_{1} = r / μ_{1} = 10$ and $ρ / θ_{2} = r / μ_{2} = 10^{- 3}$ . Our results show that by integrating the genomic marker 2 which mutation rate is larger than the recombination rate, estimates of the recombination rate as well as past population size variation are substantially improved (Table 2, Figure 2C). Indeed, analyzing only marker 1, eSMC2 and MSMC2 identify the bottleneck (albeit smoothed) and only slightly overestimate recent population size (Figure 2D). By integrating the hyper-mutable marker 2, our SMCtheo approach correctly infers the strength and time of the bottleneck when $μ_{1}$ and $μ_{2}$ are known (Figure 2D, green line), while the timing of the bottleneck is slightly shifted in the past when $μ_{2}$ is unknown and estimated by our method (Figure 2D, blue line). When $μ_{2}$ is unknown, SMCtheo additionally infers a spurious sudden variation of population size between 10,000 and 100,000 generations ago. Using only marker 1, the estimates of the recombination rate are inaccurate (Table 2). To complete the visual representation and provide a quantitative assessment of inference accuracy, we compute the root mean square error (RMSE) values for demographic inference (Supplementary file 1a). We further improve the accuracy of estimation by optimizing the likelihood (LH) to estimate the recombination rate and demography compared to the classically used Baum-Welch (BW) algorithm (Table 2 and Figure 2—figure supplement 1, Figure 2—figure supplement 2). Our results demonstrate that SNPs are limiting and insufficient for accurate inferences in recent times and that the inclusion of an additional marker with a mutation rate higher than the recombination rate generates significant improvements in demographic inference. However, by directly optimizing the likelihood, the true recombination rate can be recovered well even with marker 1 only.

Table 2

Estimates of recombination rates with one or both markers.

For SMCtheo, BW stands for the use of the Baum-Welch algorithm to infer parameters, and LH for the use of the likelihood. We use 10 sequences of 100 Mb with $r = 10^{- 7}$ , $μ_{1} = 10^{- 8}$ and $μ_{2} = 10^{- 4}$ per generation per bp in a population with a past bottleneck event. The coefficient of variation over 10 repetitions is indicated in brackets.

Method	True recombination rate	Average estimated recombination rate
MSMC2 (BW)	10⁻⁷	0.23×10⁻⁷ (0.017)
1 Marker: BW	10⁻⁷	0.25×10⁻⁷ (0.012)
2 Marker: BW	10⁻⁷	0.90×10⁻⁷ (0.004)
1 Marker: LH	10⁻⁷	0.84×10⁻⁷ (0.036)
2 Marker: LH	10⁻⁷	0.94×10⁻⁷ (0.01)

Integrating DNA methylation improves the accuracy of inference

Definition of the theoretical model for DNA methylation

Following the previously encouraging results of demographic inference with SNPs and a hyper-mutable marker under the specific assumptions of Figure 1A, we develop a specific SMCm method to jointly analyse SNPs and CG methylation as an epigenetic hyper-mutable marker. Since our SMCm stems from the eSMC (Sellinger et al., 2020; Strütt et al., 2023), it corrects for the effect of self-fertilization when applying to A. thaliana. We focus here on methylation located in CG contexts within genic regions as these have been found to evolve neutrally (Vidalis et al., 2016; Yao et al., 2021; Yao et al., 2023). The methylation of individual CG dinucleotides produces a biallelic heritable marker with a finite number of (epi)mutable sites (Figure 3). In a sample of several sequences from a population, variation in the methylation status of individual CGs is known as single methylation polymorphism (SMP, Figure 3A), which could be used for demographic and divergence inference (van der Graaf et al., 2015; Vidalis et al., 2016). However, CG methylation sites can also be organized in spatial clusters (of similar state) due to region-level epimutation (Figure 3B, Weigel and Colot, 2012; Denkena et al., 2021; Muyle et al., 2021). Region-level epimutations can have different epimutation rates than individual CG sites. Population-level variation in the methylation status of these clusters is known as differentially methylated regions (DMRs). Furthermore, when integrating SMP and DMR epimutational processes (i.e. what we here call region-level epimutation), the methylation status of CG sites is therefore affected by the superposition of both processes. Therefore, the simulation and modeling of epimutation processes of SMPs are more complex than in our previous model as we need to account for the effect of region methylation as well as for methylation and demethylation epimutation rates to be different and asymmetrical (van der Graaf et al., 2015; Denkena et al., 2021).

Figure 3

Download asset Open asset

Schematic representation of site and region epimutations.

Schematic representation of a sequence undergoing epimutation at (A) the cytosine site level and (B) at the region level. A methylated cytosine in CG context is indicated in black, and an unmethylated cytosine in white.

To make our simulations realistic, we use the A. thaliana genome sequence as a starting point, and focus on CG dinucleotides within genic regions. To that end, we select random 1 kb regions within genes and choose only those CG sites that are clearly methylated or unmethylated in A. thaliana natural populations based on whole genome bisulphite sequencing (WGBS) measurements from the 1001 G project (SI text). Our simulator for CG methylation is built in a similar way as the one described above, but the epimutation rates are allowed to be asymmetric with the per-site methylation rate ( $μ_{S M}$ ) and demethylation ( $μ_{S U}$ ). Region-level epimutations are also implemented, setting the region length to either 1 kb (Muyle et al., 2021) or 150 bp (Denkena et al., 2021). The region-level methylation and demethylation rates are defined as $μ_{R M}$ and $μ_{R U}$ , respectively. We assume that site-level and region-level epimutation processes are independent. Making this assumption explicit later allows us to test if it is violated in comparisons with actual data. Our simulator also assumes that DNA mutations and epimutations are independent of one another. That is, for simplicity, we ignore the fact that methylated cytosines are more likely to transition to thymines as a result of spontaneous deamination (Johannes, 2019a). We also ignore the possibility that new DNA mutations could act as CG methylation quantitative trait loci and affect CG methylation patterns in both cis and trans. Such events are extremely rare, so the above assumptions should hold reasonably well over short evolutionary time scales. As the goal is to apply our approach to A. thaliana, we simulate sequence data for a sample size $n = 10$ (but considering A. thaliana haploid) from a population displaying 90% selfing (Sellinger et al., 2020) under a recent severe population bottleneck demographic scenario. We simulate data assuming previous estimates of the rates of recombination (Salomé et al., 2012), DNA mutation (Ossowski et al., 2010), and site- and region-level methylation (van der Graaf et al., 2015; Denkena et al., 2021).

As guidance for future analyses of demographic inference using SNPs and DNA methylation data, the theoretical and empirical analysis of A. thaliana methylomes consist of the following five steps: (1) assessing the relevance of region-level methylation (DMRs) for inference, (2) inference of site and region epimutation rates, (3) comparing statistics for the SNPs, SMPs, and DMRs distributions, (4) demographic inference using SNPs with SMPs or DMRS, and (5) demographic inference using SNPs with SMPs and DMRs.

Step 1: assessing the relevance of region-level methylation (DMRs) for inference

We determine our ability to detect the existence of spatial correlations between epimutations. That is, we asked if site-specific epimutations can lead to region-level methylation status changes across a range of epimutation rates (assuming two sequences of 100 Mb, $r = μ_{1} = 10^{- 8}$ per generation per bp and a constant population size $N = 10, 000$ , results in Supplementary file 1b). If site-specific epimutations are independently distributed, the probability of a given site being in a given (methylated or unmethylated) state should be independent of the state of nearby sites (knowing the epimutation rate per site). Conversely, if there is a region effect on epimutation (DMRs), two consecutive sites along the genome would exhibit a positive correlation in their methylated states. We therefore calculate from the per-site (de)methylation rates $μ_{S M}$ and $μ_{S D}$ the probability that two successive cytosine positions are identical in their methylation assuming they are independent. This probability can be compared to the one observed from methylation data (here simulated) so that we obtain a statistical test for the existence of a positive correlation in the methylation status of nearby sites, interpreted as a regional-level epimutation process (p-value = 0.05) according to Figure 1A. A small p-value of the test (<0.05) suggests the existence of a region effect for methylation/demethylation affecting neighbouring cytosines, contrary to a high p-value indicating no spatial structure of methylation distribution. We find that when region epimutation rates are higher than (or similar to) site-level epimutation rates, namely $μ_{R M} ⪆ μ_{S M}$ and ( $μ_{R U} ⪆ μ_{S U}$ ), the existence of regions of consecutive cytosines is detected with high accuracy. However, when site-level epimutation rates are higher ( $μ_{S U} > μ_{R U}$ and $μ_{S M} > μ_{R M}$ ) than region-level epimutation rates, region-level changes cannot be readily detected (Supplementary file 1b). When methylated regions are detected, we can further determine their length using a specifically developed HMM using all pairs of genomes (similarly to Shahryary et al., 2020; Denkena et al., 2021; Taudt et al., 2018). While the length of the methylated region is pre-determined in our simulations (1 kb or 150 bp), site-level epimutation occur which can change the distribution of methylation states in that region and across individuals, thus DMR regions can vary in length along the genome and between pairs of chromosomes.

Step 2: inference of site- and region-level epimutation rates

As the epimutation rates of most plant species remain unknown, we assess the accuracy of SMCm to infer epimutation rates at the site- and region-level directly from simulated data. We first assume that either only site- or only region epimutations can occur and infer their respective rates (see Supplementary file 1c and d). Our SMCm approach can accurately recover these rates except when these are higher than 10^-4. Next, we assess the accuracy of our approach to simultaneously infer site- and region-level epimutation rates, assuming that region and site epimutation rates are equal (Supplementary file 1e and Figure 4—figure supplement 1). Similar to our previous observation, we find that when the epimutation rates are very high (e.g. close to 10^-2), accuracy is lost compared to slower epimutation rates. Nonetheless, our average estimated rates are off from the true value by less than a factor 10. Hence, under our model assumptions, we can recover the correct order of magnitude for site- and region-level methylation and demethylation rates.

Step 3: distribution of statistics for SNPs, SMPs, and DMRs

To gain insights on the distribution of epimutations under the described assumptions, we look at key statistics from our simulations: the distribution of distance between two recombination events versus the distribution of the length of estimated DMR regions (Figure 4A), and the LD decay for SMPs (in genic regions) and SNPs (in all contexts) (Figure 4C and D). In our simulations, DMR regions have a maximum fixed size, but their length depends on the interaction between the region- and site-level epimutation rates. As mentioned in step 1, the methylated/demethylated regions are detected using the binomial test and their length estimated by the HMM. Therefore, while variation exists for the length of these regions (Figure 4A), regions are on average shorter than the span of genealogies along the genome, which are defined by the frequency of recombination events along the genome ( $r = 3.5 \times 10^{- 8}$ as in A. thaliana). There is virtually no linkage disequilibrium (LD) between epimutations due to the high epimutation rate (Figure 4C), while the LD between SNPs can range over few kbp (Figure 4D, as observed in A. thaliana [Cao et al., 2011; Schmitz et al., 2013]). Note, however, that the region methylation process in itself does not generate LD because this measure can only be computed if SMPs are present in frequency higher than $2 / n$ in the sample, that is there is no LD measure defined for monomorphic methylated/unmethylated regions. In other words, our simulator generates SNPs, SMPs, and DMRs, which fulfil the three key assumptions from Figure 1A. We note that by using a constant population size $N = 10, 000$ , the LD decay for SNPs is higher than in the A. thaliana data which exhibit an effective population size of ca. $N = 250, 000$ (Cao et al., 2011) and past changes in size.

Figure 4 with 1 supplement see all

Download asset Open asset

Key statistics for epimutations and mutations.

(A) Histogram of the length between two recombination events (genomic span of a genealogy) and DMRs size in bp of the simulated data. (B) Histogram of genealogy span and DMRs size in bp from the *A. thaliana* data (10 German accessions). (C) Linkage disequilibrium decay of epimutations in our samples of *A. thaliana* (red) and simulated data (blue). (D) Linkage disequilibrium decay of mutations in our *A. thaliana* samples (red) and simulated data (blue). The simulations reproduce the outcome of a recent bottleneck with sample size $n = 5$ diploid of 100 Mb, the rates per generation per bp are $r = 3.5 \times 10^{- 8}$ , $μ_{1} = 7 \times 10^{- 9}$ , $μ_{S M} = 3.5 \times 10^{- 4}$ , $μ_{S U} = 1.5 \times 10^{- 3}$ , and per 1 kb region $μ_{R M} = 2 \times 10^{- 4}$ and $μ_{R U} = 1 \times 10^{- 3}$ .

Step 4: demographic inference based on SNPs with SMPs or DMRs

We test the usefulness of either SMPs or DMRs for demographic inference. Simulations under the demographic model from steps 1–3 assume DNA mutations (SNPs) and only site epimutations (SMPs), that is no region-level methylation ( $μ_{R M} = μ_{R U} = 0$ ). We perform inference of past demographic history under different amounts of potentially methylated sites with and without a priori knowledge of the methylation/demethylation rates (Figure 5A and B). When the site epimutation rates are a priori known, the sharp decrease of population size can be accurately detected. When epimutation rates are unknown, the shape of the past demographic history is also well inferred except for a scaling issue (a shift along the x- and y-axes similar to that in Figure 5D). When we vary the amount of potentially methylated sites (2%, 10%, and 20%) our inference results remain largely unchanged. This suggests that having methylation measurements for as low as 2% of all CG sites being epimutable in the genome is entirely sufficient to improved SNP-based demographic inference (eSMC2 in Figure 5A). The RMSE values for demographic inference are computed for all cases in Figure 5 to provide an additional quantitative understanding of our results (Supplementary file 1f).

Figure 5 with 5 supplements see all

Download asset Open asset

Performance of SMC approaches using site epimutations (SMPs) and mutations (SNPs) under a bottleneck scenario.

Estimated demographic history by eSMC2 (blue) and SMCm assuming the epimutation rate is known (**B and D**) or not (**A and C**) where the percentage of CG sites with methylated information varies between 20% (red), 10% (orange) and 2% (green) using 10 sequences of 100 Mb in **A and B** (with 10 repetitions) and 10 sequences of 10 Mb in **C and D** (three repetitions displayed) under a recent severe bottleneck (black). The parameters are: $r = 3.5 \times 10^{- 8}$ per generation per bp, mutation rate $μ_{1} = 7 \times 10^{- 9}$ , methylation rate to $μ_{S M} = 3.5 \times 10^{- 4}$ and demethylation rate to $μ_{S U} = 1.5 \times 10^{- 3}$ per generation per bp.

The amount of sequence data used in Figure 5A, B is fairly large compared to real datasets (10 haploid genomes of length 100 Mb). We therefore ran the SMCm and eSMC2 on sequence data simulated under the same scenario but with a reduced sequence length of 10 Mb ( $n = 5$ diploid, Figure 5C and D, only 3 repetitions are presented for visibility). In this case, we found that inference is significantly affected when using only SNPs (eSMC2 in blue), as we are unable to correctly recover the demographic scenario. However, incorporating SMPs with known site-level epimutations into the model leads to substantial inference improvements (Figure 5C, D, Supplementary file 1f). We additionally quantify the accuracy gain in ARG inference by inferring the expected coalescent time (TMRCA) at each position in the genome by the three approaches (eSMC2, SMCm with unknown epimutation rates and SMCm with known epimutation rates) under the same scenario from Figure 5. The RMSE values of the TMRCA inference are presented in Supplementary file 1g. We confirm our intuition that integrating epimutations slightly improves the accuracy of TMRCA when the epimutation rates are known but does not when the rates are unknown.

To quantify the effect of DMRs on inference, we simulate data under the same demographic scenario but assume only region-level epimutations (DMRs, $μ_{S M} = μ_{S U} = 0$ ). The results for DMR region sizes 1 kb and 150 bp are displayed in Figure 5—figure supplement 1 and Figure 5—figure supplement 1, respectively. As in Figure 5, we observed a gain of accuracy in inference when region-level epimutation rates are known, while the length of the region (1 kb or 150 bp) does not seem to affect the result. However, no significant gain of information is observed when integrating DMR data with unknown epimutation rates (Figure 5—figure supplement 1, Figure 5—figure supplement 2). In summary, CG methylation SMPs and to a lesser extend DMRs, can be used jointly with SNPs to improve demographic inference (Supplementary file 1h presents the corresponding RMSE values for demographic inference shown in Figure 5—figure supplement 1, Figure 5—figure supplement 2), especially in recent times (Supplementary file 1f and h).

Step 5: demographic inference based on SNPs with SMPs and DMRs

Since site- and region-level methylation processes can occur in real data, we run SMCm on simulated data under the same demographic scenario, but now using both site (SMPs) and region (DMRs) epimutations and accounting for both mutation processes (with rates similar to the one found in A. thaliana). Inference results are displayed in Figure 5—figure supplement 3 (RMSE values in Supplementary file 1i). When the epimutations rates are unknown, we observe a gain of accuracy when integrating epimutations, especially in recent times. However, when epimutation rates are a priori known, we observe a loss of accuracy when accounting for epimutations. This loss of accuracy is due to the mislabeling of the methylation region status (in step 1) when site and region-level epimutations occur jointly at similar rates (as there will be methylated sites in unmethylated regions and unmethylated sites in methylated regions).

Finally, we assess the inference accuracy when using SNPs and SMPs but ignoring in SMCm the region methylation effect (DMRs), even though this latter process takes place (Figure 5—figure supplement 4, RMSE values in Supplementary file 1j). The inference accuracy decreases compared to the previous results (Figure 5—figure supplements 1–4), and while the sudden variation of population is somehow recovered, the estimates of the time and magnitude of size change are not well recovered in recent times. Hence, those results demonstrate the importance of accounting for site and region-level epimutations processes in steps 1–5. We demonstrate that our SMCm can exhibit, to some extent, an improved statistical power for demographic inference using SNPs and SMPs while accounting for site and region-level methylation processes under the assumptions of Figure 1A. We show that (1) using SMPs we can unveil past demographic events hidden by limitations in SNPs, (2) the correct demography can be uncovered irrespective of knowing a priori the epimutation rates, (3) ignoring site or region-level processes can decrease the accuracy of inference, and (4) knowing the epimutation rates may improve the estimate of demography compared to simultaneously estimating them with SMCm.

Joint use of SNPs and SMPs improves the inference of recent demographic history in A. thaliana

Step 1: assessing the strength of region-level methylation process in A. thaliana

We apply our inference model to genome and methylome data from 10 A. thaliana plants from a German local population (Cao et al., 2011). We start by assessing the strength of a region effect on the distribution of methylated CG sites along the genome. As expected from Denkena et al., 2021, for all 10 individual full methylomes, we reject the hypothesis of a binomial distribution of methylated and unmethylated sites along the genomes, suggesting the existence of region effect methylation (yielding DMRs), meaning that CG are more likely to be methylated if in a highly methylated region, and conversely for unmethylated CG. This is consistent with the autocorrelations in mCG found in Denkena et al., 2021; Briffa et al., 2023; Lyons et al., 2023. As a first measure of methylated region length, we test the independence between two annotated CG methylations given a minimum genomic distance between them (within one genome). We observe an average p-value smaller than 0.05 for distances up to 2,000 bp but then the p-value rapidly increases (>0.4) (Figure 5—figure supplement 5). As a second measure, our HMM (based on pairs of genomes) yields a DMR average length of 222 bp (distribution in Figure 4B).

We conclude that the minimum distance for epimutations to be independent along a genome is over 2 kb and spans larger distances than the typically proposed DMR size (ca. 150 bp in Denkena et al., 2021 and 222 bp in our analysis) and can therefore cover the size of a gene (see Muyle et al., 2021; Briffa et al., 2023). The simulations and data from A. thaliana indicates that the epimutation processes that produce DMRs at the population level in plants cannot simply result from the cumulative action of single-site epimutations. This insight is consistent with recent analyses of epimutational processes in gene bodies, which seems to indicate that the autocorrelation in CG methylation is a function of cooperative methylation maintenance and the distribution of histone modifications (Briffa et al., 2023; Lyons et al., 2023).

Step 2: site- and region-level epimutation rates

We use the rates empirically estimated in A thaliana and taken in the above simulations ( $μ_{S M} = 3.5 \times 10^{- 4}$ and $μ_{S U} = 1.5 \times 10^{- 3}$ per bp per generation and $μ_{R M} = 2 \times 10^{- 4}$ and $μ_{R U} = 1 \times 10^{- 3}$ per region per generation, van der Graaf et al., 2015; Denkena et al., 2021).

Step 3: distribution statistics for SNPs, SMPs and DMRs in A. thaliana

Since our SMC model assumes that DNA, SMP and DMR polymorphisms are determined by the underlying population/sample genealogy, DMR which span long genomic regions may spread across multiple genealogies and thus violates our modeling assumptions. We thus further investigate the potential discrepancies between the data and our model (Figure 4). We infer the DMR sizes from all 10 A. thaliana accessions using our ad hoc HMM, and measure the bp distance between a change in the expected hidden state (i.e. coalescent time) along the genome, which we interpret as recombination events (called the genomic span of a genealogy). The resulting distributions are found in Figure 4B. We observe that both distributions have a similar shape, but DMRs are, on average, twice as large as the inferred genomic genealogy span: average length of 222 bp (DMR) vs 137 bp (genealogy) and median length of 134 bp (DMR) vs 62 bp (genealogy). This means that, on average, DMRs are larger than the average distance between two recombination events, thus violating the homogeneous distribution of epimutations along the genome (Figure 1C).

To further unveil potential non-homogeneity of the distribution of epimutations, we assess the decay of LD of mutations (SNPs) and epimutations (SMPs) (Figure 4C and D), confirming the results in Schmitz et al., 2013. We find the LD between SMPs in the data to be high (and higher than LD between SNPs) for distances smaller than 100 bp (red line in Figure 4C, D). The LD decay of SMPs is much faster than for SNPs (no linkage disequilibrium between epimutations for distances >100 bp), likely stemming from (1) epimutation rates being much higher than the DNA mutation rate, and (2) the high per-site recombination rate in A. thaliana. Moreover, the LD between SMPs at distances smaller than 100 bp in A. thaliana being much higher compared to our simulations (Figure 4C), we suggest that additional local mechanisms of epimutation processes may not be accounted for in our model of the region-level methylation process.

Step 4: demographic inference for A. thaliana based only on SNPs and SMPs

Finally, we apply the SMCm approach to data from the German accessions of A. thaliana. When using SNP data only, the demographic results are similar to those previously found (Sellinger et al., 2020; Strütt et al., 2023; Figure 6 purple lines), with no strong evidence for an expansion post-Last Glacial Maximum (LGM) (Cao et al., 2011). We then sub-sample and analyze segregating SMPs, which exhibit both methylated and unmethylated states in our sample (as in van der Graaf et al., 2015). Here we ignore DMRs and account only for SMPs. When we use as input the methylation and demethylation rates that have been inferred experimentally (van der Graaf et al., 2015), a mild bottleneck post-LGM is followed by recent expansion (Figure 6 blue lines). By contrast, letting our SMCm estimate the epimutation rates, we find in recent times a somehow similar but stronger demographic change post-LGM. We find a strong bottleneck event occurring between ca. 5,000 and 10,000 generations ago, followed by an expansion until today (Figure 6 green lines). The inferred site epimutation rates are 10,000 faster than the DNA mutation rate (Supplementary file 1k), which is close to the expected order of magnitude from experimental measures with and without DMR effects (van der Graaf et al., 2015; Denkena et al., 2021). Both estimates thus yield a post-LGM bottleneck followed by a recent population expansion. These results indicate that the inclusion of DNA methylation data can aid in the accurate reconstruction of the evolutionary history of populations, particularly in the recent past, where SNPs reach their resolution limit. This is made possible by the fact that the DNA methylation status at CG dinucleotide undergoes stochastic changes at rates that are several orders of magnitude higher than the DNA mutation rate and can be inherited across generations similar to DNA mutations.

Figure 6 with 6 supplements see all

Download asset Open asset

Integrating epimutations and mutations on German accessions of A. *thaliana*.

Estimated demographic history of the German population by eSMC2 (only SNPs, purple) and SMCm when keeping polymorphic methylation sites (SMPs) only: green with epimutation rates estimated by SMCm, blue with epimutation rates fixed to empirical values. The region epimutation effect is ignored. The parameters are $r = 3.6 \times 10^{- 8}$ , $μ_{1} = 6.95 \times 10^{- 9}$ , and when assumed known, the site methylation rate is $μ_{S M} = 3.5 \times 10^{- 4}$ and demethylation rate is $μ_{S U} = 1.5 \times 10^{- 3}$ .

Step 5: demographic inference correcting for DMRs in A. thaliana

To assess the robustness of our inference results, we run SMCm using all cytosines (CG) sites with an annotated methylation status (segregating or not) while accounting or not for DMRs (Figure 6—figure supplement 1). We fix epimutation rates to the empirically estimated values and confirm the estimates from Figure 6. When the region-level methylation process is ignored, the inferred demography (blue lines in Figure 6—figure supplement 1) is similar to the estimates from SMPs with fixed rates in Figure 6 (blue lines). When the region-level methylation process is taken into account (orange lines in Figure 6—figure supplement 1), the inferred demography is similar to that of Figure 6 (green lines). In the case where we infer the epimutation rates (sites and region) the demographic history inference is not improved compared to that estimated using SNPs only (Figure 6—figure supplement 1, green and red lines) while the inferred epimutation rates are lower than expected (Supplementary file 1k and l), but the ratio of site to region epimutation rates is consistent with empirical estimates (Denkena et al., 2021).

Discussion

Current approaches analyzing whole genome sequences rely on statistics derived from the distribution of ancestral recombination graphs (Gattepaille et al., 2016; Sellinger et al., 2021; Korfmann et al., 2024; Strütt et al., 2023; Brandt,, 2022; Wohns et al., 2022; Speidel et al., 2019; Kelleher et al., 2019). In this study we present a new SMC method that combines SNP data with other types of genomic (TEs, microsatallites) and epigenomic (DNA methylation) markers. We focus mainly on the inclusion of genomic markers whose mutation rates exceed the DNA point mutation rate, as such (hyper-mutable) markers can provide increased temporal resolution in the recent evolutionary past of populations, and aid in the identification of demographic changes (e.g. population bottlenecks). We demonstrate that by integrating multiple heritable genomic markers, the population size variation in very recent time can be more accurately recovered (outperforming any other methods given the amount of data used in this study [Terhorst et al., 2017; Speidel et al., 2019]). Our results indicate that correctly integrating multiple genomic marker can improve TRMCA inference, which is becoming a field of high interest (Korfmann et al., 2024; Hubisz et al., 2020; Mahmoudi et al., 2022). Our simulations demonstrate that if the SNP mutation rate is known, the mutation rate of other markers can be recovered (under the condition that the marker follow all hypotheses described in Figure 1). Moreover, our method accounts for the finite site problem that arises at reversible (hyper-mutable) markers and/or where effective population size is high (Tellier et al., 2011; Upadhya and Steinrücken, 2022). Overall, the simulator and SMC methods presented here therefore pave the way for a rigorous statistical framework to test if a common ARG can explain the observed diversity patterns under the model hypotheses laid out in Figure 1. We find that comparisons of LD for different markers along the genome is a useful way to assess violations of our model assumptions.

As proof of principle, we apply our approach on data originating from whole genome and methylome data of A. thaliana natural accessions (focusing on CG context in genic regions, as in Vidalis et al., 2016; Yao et al., 2021; Yao et al., 2023). Indeed, A. thaliana presents the largest genetic and epigenetic data-set of high quality. Additionally, the methylation states in CG context has been proven mainly heritable and is well documented (Denkena et al., 2021; Hazarika et al., 2022; van der Graaf et al., 2015). We first investigate the distribution of epimuations along the genomes. Our model-based approach provides strong evidence that DMRs cannot simply emerge from spontaneous site-level epimutations that arise according to a Poisson processes along genome. Instead, stochastic changes in region-level methylation states must be the outcome of spontaneous methylation and demethylation events that operate at both the site- and region-level (as corroborated by Pisupati et al., 2023; Briffa et al., 2023; Lyons et al., 2023). Our epimutation model cannot fully describe the observed diversity of epimutations along the genome (Denkena et al., 2021), meaning that the epimutation processes may indeed be more complex than expected (Denkena et al., 2021; Hazarika et al., 2022; Briffa et al., 2023; Lyons et al., 2023). We observe non-independence between annotated methylation sites spanning genomic regions larger than the span of the underlying genealogy (determined by recombination events) which no model can currently describe. Additionally, we find high LD between SMPs over short distances which does not appear in our simulations (simulation performed under the current measures of epimutation rates). Thus, methylation probably violate the assumptions of a Poisson process distribution along the genome and in time (i.e. Figure 1), in line with recent functional studies (Pisupati et al., 2023; Hazarika et al., 2022; Lyons et al., 2023). We thus further caution against conclusions on the role of natural (purifying) selection (Muyle et al., 2021) or its absence (Vidalis et al., 2016) based on population epigenomic data due to the violation of the above-mentioned assumptions. Additionally, we suspect those model violations to explain the discrepancy between the epimutation rates we inferred and the ones measured experimentally (van der Graaf et al., 2015; Denkena et al., 2021). To solve this discrepancy, one would need to develop a theoretical epimutation model capable of describing the observed diversity at the evolutionary time scale and then use this model to reanalyse the sequence data from the biological experiment to re-estimate the epimutation rates. We thus suggest a possible way forward for modeling epimutations through an Ising model (Zhang et al., 2018) to account for the heterogeneous methylation process. However, our preliminary work and the simulation results in Briffa et al., 2023, indicate that such model generates non-homogeneous mutation process in space (i.e. along the genome) and time, violating our current SMC assumptions (Figure 1C, D). Hence, there is a need to develop a more realistic methylation model for epimutations. A model accounting for heterogeneous rates would probably need to rely on a more sophisticated HMM (e.g. continuous time Markov chains Ki and Terhorst, 2023 for SMC approaches) than what is presented here or to use other full genome inference methods (see Korfmann et al., 2024) which are not constrained by the SMC assumptions (Figure 1) but depends on simulations.

Interestingly, the distance of LD decay for SMPs matches quite well the estimated distance between recombination events (Figure 4). In addition to our theoretical results in Table 2, this observation reinforces the usefulness of using SMPs (or any hyper-mutable marker) to improve estimates of the recombination rate along the genome in species where the per site DNA mutation rate (μ) is smaller than the per site recombination rate ( $r$ ) as in A. thaliana.

Nonetheless, we find that a restricted focus on segregating SMPs in genic regions could meet our model assumptions reasonably well, and thus provides a promising way forward. Using these segregating SMPs, we recover a past demographic bottleneck followed by an expansion which could fit the post- Last Glacial Maximum (LGM) colonization of Europe (although caution must be taken concerning the reliability of those results as pointed above), a hypothesized scenario (François et al., 2008) which could not be clearly identified using SNPs only from European (relic and non-relic) accessions (Cao et al., 2011). Currently strong evidence from inference methods are lacking (Cao et al., 2011, Figure 4 in Durvasula et al., 2017). Indeed, beyond the limits of using SNPs only, current results are limited by theoretical frameworks unable to simultaneously account (and disentangle) for extensive background selection (reinforced by very high selfing), population structure and variation in molecular rates (e.g. mutation rates, Monroe et al., 2022), which are all known to be present in A. thaliana. Those various forces are known to bias inference results when non-accounted for (Charlesworth and Jensen, 2023; Rodríguez et al., 2018), and may explain the variance in our demographic estimates. We also note that using CG methylated sites in genic regions may be problematic as the typical genealogies at these loci could be shorter than the genome average due to the presence of background selection, thus making the inference of such short TMRCA more difficult (even with SMPs) than in non-coding regions (which do not harbour desirable CG methylation sites, van der Graaf et al., 2015; Vidalis et al., 2016; Yao et al., 2021).

We suggest that simultaneously accounting for multiple heritable markers can help disentangle different evolutionary forces, such as between selection and variation in mutation rate: selection has a local effect on the population genealogy, while the mutation rate variation would only locally affect that given marker but not the genealogy (Charlesworth and Jensen, 2023). The absence of conflicting demography inferred from SNPs and from methylation confirms at the time scale of thousands of generations, CG methylation sites are mainly heritable and can be modeled using population genetics theory (Charlesworth and Jain, 2014; Vidalis et al., 2016 but see Pisupati et al., 2023) and used to estimate divergence between lineages (Yao et al., 2023; Yao et al., 2021). In other words, fast ecological local adaptation (Schmid et al., 2018) and response to stresses (Srikant and Drost, 2020) may likely not be prominent forces endlessly reshaping CG methylation patterns (non-heritability in Figure 1B).

Overall, our results demonstrate that our approach can be used in different cases. If the epimutations/genomic markers evolutionary mechanisms are not well understood (Pisupati et al., 2023; Briffa et al., 2023; Lyons et al., 2023), our approach provides inference tools to study the markers’ rates and distribution process along the genome, without requiring additional experimental data. If the evolution of epimutations/genomic markers are well understood (including a measure of the mutation rates) and can be modeled to describe the observed intra-population diversity, these can be integrated to improve the SMC performance. Hence, when applying our approach to genome-wide genetic and epigenetic data, it is advisable to use accurately annotated markers with, if possible, information regarding their inheritance and mutational properties. Regarding methylation specifically, while the set of gene body methylated genes previously used (Vidalis et al., 2016; Yao et al., 2023) are likely the optimal choice (Yao et al., 2021), these are too few and too scattered across the genome to maximize the statistical power of SMC methods. We, therefore, use methylation sites at all genic regions. Yet, despite the wealth of functional studies and data on methylation in A. thaliana, the distribution of epimutations is not fully understood (Hazarika et al., 2022; Pisupati et al., 2023), but independent rates for sites and region-level have been estimated (van der Graaf et al., 2015; Denkena et al., 2021; Yao et al., 2023). We note here the promising methylation modeling framework by Briffa et al., 2023; Lyons et al., 2023, albeit it does not yet consider evolutionary processes at the population level. Our results shed light on the inference accuracy in presence of site and region-level epimutations when occurring at similar rates (Figure 5—figure supplement 3). When accounting for region-level epimutations, our algorithm requires first inferring via an HMM the methylation status of a region in order to later on compute the epimutation probabilities (i.e. the emission matrix of the SMC HMM). Hence, in the presence of site and region-level epimutations occurring at similar rates, recovering the region methylation status becomes harder as methylated sites are observed in the unmethylated regions (and unmethylated sites observed in the methylated regions). The mislabelling of the region methylation status lead to accuracy loss due to the use of the wrong emission probability at the later steps of the SMC inference (Forward-Backward algorithm). When epimutation rates are freely inferred, their values are based on the estimated methylation region status. Therefore, even if the inferred rates are incorrect, these are sufficiently consistent with the inferred region methylation status to contain information and slightly improve inference accuracy. Additionally, extra care must be taken when dealing with epigenomic data in other species as the SMP calling might not be as simple as for Arabidopsis thaliana due to potential difference of methylation between different tissues or pool of cells. Similarly, we ignore here the potential dependence between SNPs and SMPs, as more empirical evidence (and modeling) is required to quantify the potential interaction between both mutational processes.

On a brighter note, with the release of new sequencing technology (Lang et al., 2020), long and accurate reads are becoming accessible, leading to the availability of highquality reference genomes for model and non-model species alike (Nurk et al., 2020). Additionally, the quality of re-sequencing (population sample) genome data and their annotations is enhanced so that additional markers such as transposable elements, insertion, deletion or microsatellites can be called with increasing confidence. These accurate genomes will provide access to new classes of genomic markers that span the entire mutational spectrum. We therefore suspect that, in the near future, there will be an improvement in our understanding of the heritability of many markers besides SNPs. Adding other genomic markers besides SNPs will improve full genome approaches, which are currently limited by the observed nucleotide diversity (Kelleher et al., 2019; Speidel et al., 2019; Schweiger and Durbin, 2023). Additionally, the potential complexity resulting by integrating multiple independent markers could be tackled by the use of continuous time Markov chains for the emission matrix. We predict that our results pave the way to improve the inference of (1) biological traits or recombination rate through time (Deng et al., 2021; Strütt et al., 2023), (2) multiple merger events (Korfmann et al., 2024), and (3) recombination and mutation rate maps (Barroso et al., 2019; Barroso and Dutheil, 2023). Our method also should help to dissect the effect of evolutionary forces on genomic diversity (Johri et al., 2022; Johri et al., 2021), and to improve the simultaneous detection, quantification and dating of selection events (Albers and McVean, 2020; Bisschop et al., 2021; Johri et al., 2020).

Hence, there is no doubt that extending our work, by simultaneously integrating diverse types of genomic markers into other theoretical framework (e.g. ABC approaches), likely represents the future of population genomics, especially to study species for which many thousands of samples cannot be obtained. We believe our approach helps to develop more general classes of models capable of leveraging information from any type and amount of diversity observed in sequencing data, and thus to challenge our current understanding of genome evolution.

Materials and methods

Simulating two genomic markers

Request a detailed protocol

The sequence is written as a sequence of markers with a given state. Each site is annotated as M $X$ S $Y$ , where $X$ indicates the marker type and $Y$ is the current state of that marker: for example, M1S1 indicates a marker of type 1 in state 1 at this position. To simulate sequence of theoretical marker we start by simulating an ARG which is then split in a series of genealogies (i.e. a sequence of coalescent trees) along the chromosome and create an ancestral sequence (based on equilibrium probability of marker states). Mutation events (nucleotides or epimutations for methylable cytosine) are added along the sequence, that is along the series of genealogies. The ancestral sequence is thus modified by mutation event assuming a finite site model (Yang, 1996) conditioned to the branch length and topology of the genealogies. Each leaf of the genealogy is one of the $n$ samples. Our model has thus two important features: (1) markers are independent from one another, and (2) a given marker has a polymorphism distribution between samples (frequencies of alleles) determined by one given genealogy. The simulator can be found in the latest version of eSMC2 R package (https://github.com/TPPSellinger/eSMC2; copy archived at Sellinger, 2024a).

Simulating methylome data

Request a detailed protocol

We now focus on methylation data located at cytosine in CG context within genic regions. Only, CG sites in those regions are considered ‘methylable’, and CG sites outside those defined genic regions do not have a methylation status and are considered ‘unmethylable’. We vary the percentage of CG site with methylation state annotated from 2 to 20% of the sequence length. The simulator can in principle simulate epimutations in different methylation context and different rates (Lister et al., 2008; Cokus et al., 2008; Zilberman et al., 2007; Zhang et al., 2006). We simulate epimutations as described above but with asymmetric rates: the methylation rate per site is $μ_{S M} = 3.5 \times 10^{- 4}$ , and the demethylation rate per site is $μ_{S M} = 1.5 \times 10^{- 3}$ (van der Graaf et al., 2015; Denkena et al., 2021). For simplicity and computational tractability, we assume that when an epimutation occurs, it occurs on both DNA strands which then present the same information. In other words, for a haploid individual, a cytosine site can only be methylated or unmethylated (as in Taudt et al., 2018). For region level epimutations, the region length is either 1kbp (Muyle et al., 2021) or 150 bp (Denkena et al., 2021). The region level methylation and demethylation rates are set to $μ_{R M} = 2 \times 10^{- 4}$ and $μ_{R U} = 10^{- 3}$ respectively (similar to rates measured in A. thaliana, Denkena et al., 2021). In addition to this, unlike for theoretical marker described above, mutations, site and region epimutations can occur at the same position of the sequence.

To simulate methylation data, we start with an ancestral sequence of random nucleotide and then randomly select regions in which CG sites have their methylation state annotated (representing the genic regions). Cytosine in CG context in those regions are either methylated or unmethylated (noted as M or U). Cytosine in other context or regions are considered as non-methylable (and noted as C). The ancestral methylation state is then randomly attributed according to the equilibrium probabilities. Our simulator then introduces DNA mutations, site- and region-epimutations in a similar way as described above.

SMC methods

Request a detailed protocol

All three methods (eSMC2, SMCtheo, and SMCm) are based on the same mathematical foundations and implemented in a similar way within the eSMC2 R package (https://github.com/TPPSellinger/eSMC2; Strütt et al., 2023; Korfmann et al., 2024; Sellinger et al., 2021; copy archived at Sellinger, 2024a). This allows to specifically quantify the accuracy gained by accounting for multiple genomic markers.

SMC optimization function

All current SMC approach rely on the Baum-Welch (BW) algorithm for parameter estimation in order to reduce computational load (as described in Terhorst et al., 2017). Yet, the Baum-Welch algorithm is an Expectation-Maximization algorithm, and can hence fall in local extrema when optimizing the likelihood. We alternatively extend SMCtheo to estimate parameters by directly optimizing the likelihood (LH) at the greater cost of computation time (even when using the speeding techniques described in Sand et al., 2013). We run this approach on a sub-sample of size six haploid genomes to limit the required computational time.

eSMC2 and MSMC2

Request a detailed protocol

SMC methods based on the PSMC’ (Schiffels and Durbin, 2014), such as eSMC2 and MSMC2, focus on the coalescent events between two individuals (i.e. two haploid genomes or one diploid genome). The algorithm moves along the sequence and estimates the coalescence time at each position by assessing whether the two sequences are similar or different at each position. If the two sequences are different, this indicates a mutation took place in the genealogy of the sample. The intuition being that the absence of mutations (i.e. the two sequences are identical) is likely due to a recent common ancestor between the sequences, and the presence of several mutations likely reflects that the most recent common ancestor of the two sequences is distant in the past. In the event of recombination, there is a break in the current genealogy and the coalescence time consequently takes a new value according to the model parameters (Marjoram and Wall, 2006; Schiffels and Durbin, 2014). A detailed description of the algorithm can be found in Malaspinas et al., 2016; Sellinger et al., 2020.

SMCtheo based on several genomic markers

Request a detailed protocol

Our SMCtheo approach is equivalent to PSMC’ but takes as input a sequence of several genomic markers. The algorithm goes along a pair of haploid genomes and checks at each position which marker is observed and whether both states of the marker are identical or not. The approach is identical to the one described above, except that the probability of both sequences being identical at one site depends on the mutation rate of the marker at this site (Equation 1). While the mutation rates for many heritable genomic markers are unknown, there is an increasing amount of measures of the DNA (SNP) mutation rate for many species. Our SMCtheo approach is able to leverage the information from the distribution of one theoretical marker (e.g. mutations for SNPs) to infer the mutation rate of the other marker 2 (assuming both mutation rates to be symmetrical). If more than 1% of sites are polymorphic in a sequence, we use the finite site assumption. If not, then from the diversity observed, the different mutation rates can be recovered by simply comparing Waterson’s theta ( $θ_{W}$ ) between the reference marker (i.e. with known rate) and the marker with the unknown rates. For example, if the diversity ( $θ_{W}$ ) at marker 2 is smaller by a factor ten than the reference marker 1 (and no marker violates the infinite site hypothesis), the mutation rate of marker 2 is inferred to be ten times smaller (corrected by the number of possible states). However, if the marker 2 violates the infinite site hypothesis, a Baum-Welch algorithm is run to infer the most likely mutation rates under the SMC to overcome this issue (the Baum-Welch algorithm description can be found in Sellinger et al., 2020).

SMCm

Request a detailed protocol

When integrating epimutations, the number of possible observations increases compared to eSMC2. As in eSMC2, if the two nucleotides (DNA mutation) at one position are identical at a non-methylable site, we indicate this as 0. If the two nucleotides differ, it is indicated as 1 (i.e. a DNA mutation occurred). When assuming site-level epimutation only, three possible observations are possible at a given methylable posisiton: (1) if the two cytosines from the two chromosomes are unmethylated, it is indicated as a 2, (2) if the two cytosines are methylated, it is indicated as a 3, and (3) if at a position a cytosine is methylated and the other one unmethylated, it is indicated as a 4. Depending on the mutation, methylation and, demethylation rates, different frequencies of these states are possible in the sample of sequences, which provide information on the emission rate in the SMC method. When both site- and region-level methylation processes occur, the methylation state is conditioned by the region-level methylation state (increasing the number of possible observations to 9).

To choose the appropriate settings for SMCm (i.e. if there are region-level epimutations), we test if the methylation state is distributed independently from one another along one genome. In the absence of a region methylation effect, the probability at each site (position) to be methylated or unmethylated should be independent of the previous position (or any other position). Conversely, if there is a region effect on epimutation, two consecutive sites along one genome would exhibit a positive correlation in their methylated states (and across pairs of sequences). We, therefore, calculate the probability that two successive positions with an annotated methylation state would be identical under a binomial distribution of methylation along a given genome. We then compare theoretical expectations to the observed data and build the statistical test based on a binomial distribution of probabilities. If the existence of region-level epimutation is detected, the regions level methylation states are recovered through a HMM similar to Shahryary et al., 2020; Denkena et al., 2021; Taudt et al., 2018. Note that this HMM model does not include information from epimutation rates known from empirical studies. The complete description of the mathematical models and probabilities are in the Appendix 1 and 2.

We postulate that the epimutation rates remain unknown in most species, while the DNA mutation rate may be known (or approximated based on a closely related species). Hence, we develop an approach based on the SMC capable of leveraging information from the distribution of DNA mutations to infer the epimutation rates (similar to what is described above). Our approach first tests if epimutations violate or not the infinite site assumptions. If less than 1% of sites with their methylation state annotated are polymorphic in a sequence, we use the infinite site assumption: the site and region level epimutation rates can be recovered straightforwardly from the observed diversity ( $θ_{W}$ , see above). Otherwise, a Baum-Welch algorithm is run to infer the most likely epimutation rates (site rate for SMP and region rates for DMRs; van der Graaf et al., 2015; Vidalis et al., 2016; Taudt et al., 2018).

Calculation of the root mean square error (RMSE)

Request a detailed protocol

To quantify the accuracy of each demographic inference, we evaluate the root mean square error (RMSE). To do so, we choose a hundred points uniformly spread across the time window (in $\log_{10}$ scale), and compare the actual population size and the one estimated by a given method at each of these points. We thus have the following formula:

R M S E = \sqrt{\frac{\sum_{i = 1}^{10^{2}} (y_{i} - y_{i}^{*})^{2}}{10^{2}}},

where $y_{i}$ is the true population size at the time point i, and $y_{i}^{*}$ is the estimated population size at the time point $i$ .

Inference of the time to the most recent common ancestor (TMRCA)

Request a detailed protocol

To infer the TMRCA at each position of the genome, we use an approach similar to the PSMC’ described in Schiffels and Durbin, 2014. We first run a forward and backward algorithm on our sequence data (see appendix of Sellinger et al., 2020; Terhorst et al., 2017 for computation details). From the output results, we calculate the probability to be in each hidden state at each position of the genome (note that the output product of the forward-backward algorithm is rescaled so that the sum of probability is one), which we use to compute the expected coalescent time at each position on the genome using the following formula:

T M R C A_{i} = \sum_{j = 1}^{n} f o_{i, j} \times b a_{i, j} \times T c_{j},

with $i$ is the position on the genome, $j$ is the hidden state index, $n$ is the number of hidden state, $f o$ is the output from the forward algorithm, $b a$ is the output from the backward algorithm, $\sum_{j = 1}^{n} f o_{i, j} \times b a_{i, j} = 1$ , and $T c$ is a vector containing all the hidden states (i.e. coalescent times).

Sequence data of Arabidopsis thaliana

Request a detailed protocol

We download genome and methylome data of A. thaliana from the 1001 genome project Cao et al., 2011. We select 10 individuals from the German accessions, respectively, corresponding to the accession numbers: 9783, 9794, 9808, 9809, 9810, 9811, 9812, 9816, 9813, 9814. We only keep methylome data in CG context and in genic regions (Vidalis et al., 2016; Denkena et al., 2021). The genic regions are based on the current reference genome TAIR 10.1. The SNPs and epimutations are called according to previously published pipeline (Taudt et al., 2018; Denkena et al., 2021). As in previous studies Sellinger et al., 2020; Fulgione et al., 2018; Durvasula et al., 2017, we assume A. thaliana data to be haploid due to high homozygosity (caused by high selfing rate). The resulting files are available on GitHub at https://github.com/TPPSellinger. To perform analysis, we chose $μ = 6.95 \times 10^{- 9}$ per generation per bp as the DNA mutation rate Ossowski et al., 2010 and $r = 3.6 \times 10^{- 8}$ as the recombination rate Salomé et al., 2012 per generation per bp. To have the most realistic model, we assume that the methylome of A. thaliana undergoes both region (RMM) and site (SMM) level epimutations (Denkena et al., 2021). When fixed, we respectively set the site methylation and demethylation rate to $μ_{S M} = 3.48 \times 10^{- 4}$ and $μ_{S U} = 1.47 \times 10^{- 3}$ per generation per bp according to van der Graaf et al., 2015. We additionally set the region level methylation and demethylation rate to $μ_{R M} = 1.6 \times 10^{- 4}$ and $μ_{R U} = 9.5 \times 10^{- 4}$ per generation per bp according to Denkena et al., 2021. Because we do not account for the effect of variable mutation or recombination rate along the genome, we cut the five chromosomes of A. thaliana into eight smaller scaffolds (Barroso and Dutheil, 2023; Barroso et al., 2019). By doing this, we remove centromeric regions and limit the effect of the variation of mutation and recombination rate along the genome. The selected regions and the SNP density (from the German accessions) are represented in Figure 6—figure supplements 2–6.

Appendix 1

Sequentially Markovian Coalescent for several markers

Our approach is a re-implementation of the PSMC’,MSMC2, eSMC and eSMC2 but accounting for different genomic markers. Hence, the Hidden Markov model is exactly the same as previously described, but with a different emission matrix. For each site, we first check what marker is present. We then set the correct substitution rate and number of hidden states. We note as $i d$ (identical) the event where both markers are in the same state, and $s e g$ if both are in different states (polymorphic). Extending the work in Yang, 1996, we have the following formula:

\begin{aligned} P (i d | γ) = \frac{1}{n b_{s}} + \frac{(n b_{s} - 1)}{n b_{s}} * e^{- 2 μ t_{γ} \frac{(n b_{s})}{(n b_{s} - 1)}} \\ P (s e g | γ) = \frac{(n b_{s} - 1)}{n b_{s}} - \frac{(n b_{s} - 1)}{n b_{s}} * e^{- 2 μ t_{γ} \frac{(n b_{s})}{(n b_{s} - 1)}} \end{aligned}

Where $μ$ is the substitution rate of the marker per N generation, $t_{γ}$ is the average coalescent time in state γ and $n b_{s}$ is the number of possible stat the marker can take.

Appendix 2

Sequentially Markovian Coalescent with DNA methylation (SMCm)

SMCm is similar to PSMC’,MSMC2, eSMC and eSMC2 but additionally accounting for epimutations.

Accounting only for site DNA epimutations

Here we assume epimutations occur similarly as mutations (under a finite site model). Since the model accounts for sequence and DNA methylation polymorphisms, there are at each position 5 different possible observations when comparing two sequences. The first observation is 0, corresponding to a non-methylable site where the two nucleotides are identical. 1, if the two nucleotides are different. 2 if it is a methylable site and both are unmethylated. 3, if the site is methylable and both are methylated. Finally, 4 is it is a methylable site and one cytosine is methylated and the other unmethylated. Therefore, after approximating the formula assuming DNA methylation state is not affected by DNA mutations, we find:

\begin{aligned} P (0 | γ) = e^{- 2 μ t_{γ}} \\ P (1 | γ) = 1 - e^{- 2 μ t γ} \\ P (2 | γ) = p_{u} \times p_{m 1} \times p_{m 1} + (1 - p_{u}) \times (1 - p_{m 2}) \times (1 - p_{m 2}) \\ P (3 | γ) = p_{u} \times (1 - p_{m 1}) \times (1 - p_{m 1}) + (1 - p_{u}) \times p_{m 2} \times p_{m 2} \\ P (4 | γ) = p_{u} \times 2 \times p_{m 1} \times (1 - p_{m 1}) + (1 - p_{u}) \times 2 \times p_{m 2} \times (1 - p_{m 2}) \\ p_{u} = \frac{μ_{u}}{μ_{u} + μ_{m}} \\ θ_{m} = (μ_{u} + μ_{m}) \times t_{γ} \\ p_{m 1} = p_{u} + (1 - p_{u}) \times e^{- θ_{m}} \\ p_{m 2} = (1 - p_{u}) + p_{u} \times e^{- θ_{m}} \end{aligned}

Where $μ$ is the mutation rate per nucleotide per N generation, $μ_{m}$ the methylation rate per generation, $μ_{u}$ the demethylation rate per generation and $t_{γ}$ the average coalescent time in state γ. Additionally, we define $p_{u}$ as the probability to be unmethylated at equilibrium, as well as $p_{m 1}$ and $p_{m 2}$ the respective probability to stay unmethylated or methylated after a time $t_{γ}$ .

Accounting only for region epimutations

Here, we assume region epimutation occurs similarly to mutations (under a finite site model). However, unlike previously, the epimutations affect multiple sites; hence, only the first position of a methylated region is considered, and the following positions will be considered as missing data (because it is one block of information). Therefore, there are at each position 6 different possible observations when comparing two sequences. The first observation is 0, corresponding to a non-methylable site where the two nucleotides are identical. 1, if the two nucleotides are different. 2 if it is a region with methylation state annotated and both regions are unmethylated. 3, if it is a region with methylation state annotated and both regions are methylated. 4 if it is a region with methylation state annotated and both regions are in different methylation states. 5 is missing data. Therefore, after approximating the formula, assuming the methylation state is not affected by mutations, we have the following formula:

\begin{aligned} P (0 | γ) = e^{- 2 μ t_{γ}} \\ P (1 | γ) = 1 - e^{- 2 μ t γ} \\ P (2 | γ) = p_{u} \times p_{m 1} \times p_{m 1} + (1 - p_{u}) \times (1 - p_{m 2}) \times (1 - p_{m 2}) \\ P (3 | γ) = p_{u} \times (1 - p_{m 1}) \times (1 - p_{m 1}) + (1 - p_{u}) \times p_{m 2} \times p_{m 2} \\ P (4 | γ) = p_{u} \times 2 \times p_{m 1} \times (1 - p_{m 1}) + (1 - p_{u}) \times 2 \times p_{m 2} \times (1 - p_{m 2}) \\ P (5 | γ) = 1 \\ p_{u} = \frac{μ_{u}}{μ_{u} + μ_{m}} \\ θ_{m} = (μ_{u} + μ_{m}) \times t_{γ} \\ p_{m 1} = p_{u} + (1 - p_{u}) \times e^{- θ_{m}} \\ p_{m 2} = (1 - p_{u}) + p_{u} \times e^{- θ_{m}} \end{aligned}

Where $μ$ is the mutation rate per nucleotide per N generation, $μ_{m}$ the region methylation rate per generation, $μ_{u}$ the region demethylation rate per generation and $t_{γ}$ the average coalescent time in state $γ$ . Additionally, we define $p_{u}$ as the probability for the region to be unmethylated at equilibrium, as well as $p_{m 1}$ and $p_{m 2}$ the respective probability for the region to stay unmethylated or methylated after a time $t_{γ}$ . To recover the region epimutations, we use a hidden Markov model (HMM). The HMM takes as input genome and methylome data described above and compares two sequences to recover epiregion. The hidden Markov model has 9 hidden states: 1 (regions with no methylation information), 2 (no methylation information in individual 1& mainly methylated region in individual 2), 3 (no methylation information in individual 1& mainly unmethylated region in individual 2), 4 (no methylation information in individual 2& mainly methylated region in individual 1), 5 (mainly methylated region in individual 1& mainly methylated region in individual 2), 6 (mainly methylated region in individual 1& mainly unmethylated region in individual 2), 7 (mainly unmethylated region in individual 1& no methylation information in individual 2), 8 (mainly unmethylated region in individual 1& mainly methylated region in individual 2), 9 (mainly unmethylated region in individual 1& mainly unmethylated region in individual 2). Additional parameters are necessary to define the transition rate. The user needs to define the minimum number of annotated methylable sites to form a region (by default 4) and minimum size of a region in bp (100 by default). Our approach then defines the transition rate as the transition rate maximizing the number of regions respecting the defined criteria. The emission matrix is defined by the user, but by default, the probability of observing methylated and unmethylated in a methylated region is set to 0.8 and 0.2 respectively. In an unmethylated region, the probabilities are respectively set to 0.2 and 0.8. We use those values throughout the study.

Accounting for site and region epimutations

Here, we assume site and region epimutation occur similarly to mutations (under a finite site model). Like previously, the epimutations affect multiple sites; hence, only the first position of a methylated region is considered, and the following positions will be considered as missing data (because it is one block of information). However, the observation will depend on the region. Therefore, there are at each position 10 different possible observations when comparing two sequences. The first observation is 0, corresponding to a non-methylable site where the two nucleotides are identical. 1, if the two nucleotides are different. 2 if it is a region with methylation state annotated and both regions are unmethylated, and all sites are unmethylated. 3 if it is a region with methylation state annotated, both regions are unmethylated, and all sites are methylated. 4 if it is a region with methylation state annotated, and both regions are unmethylated, and at least one site is segregating. 5 if it is a region with methylation state annotated and both regions are methylated, and all sites are unmethylated. 6 if it is a region with a methylation state annotated, both regions are methylated, and all sites are methylated. 7 if it is a region with methylation state annotated, both regions are methylated, and at least one site is segregating. 8 if it is a region with a methylation state annotated and both regions are in different methylation states. 9 is missing data. Therefore, we have the following formula:

\begin{aligned} P (0 | γ) = e^{- 2 μ t_{γ}} \\ P (1 | γ) = 1 - e^{- 2 μ t γ} \\ P (2 | γ) = (P_{O} (P_{O O} P_{O O})) + (P_{P} (P_{P O} P_{P O})) + (P_{M} (P_{M O} P_{M O})) + (P_{U} (P_{U O} P_{U O})) \\ P (3 | γ) = (P_{O} (P_{O P} P_{O P})) + (P_{P} (P_{P P} P_{P P})) + (P_{M} (P_{M P} P_{M P})) + (P_{U} (P_{U P} P_{U P})) \\ P (4 | γ) = (P_{O} (2 P_{O P} P_{O O})) + (P_{P} (2 P_{P P} P_{P O})) + (P_{M} (2 P_{M P} P_{M O})) + (P_{U} (2 P_{U P} P_{U O})) \\ P (5 | γ) = (P_{O} (P_{O U} P_{O U})) + (P_{P} (P_{P U} P_{P U})) + (P_{M} (P_{M U} P_{M U})) + (P_{U} (P_{U U} P_{U U})) \\ P (6 | γ) = (P_{O} (P_{O M} P_{O M})) + (P_{P} (P_{P M} P_{P M})) + (P_{M} (P_{M M} P_{M M})) + (P_{U} (P_{U M} P_{U M})) \\ P (7 | γ) = (P_{O} (2 P_{O U} P_{O M})) + (P_{P} (2 P_{P U} P_{P M})) + (P_{M} (2 P_{M U} P_{M M})) + (P_{U} (2 P_{U U} P_{U M})) \\ P (8 | γ) = (p_{u} 2 p_{m 1} (1 - p_{m 1})) + (1 - p_{u}) 2 p_{m 2} (1 - p_{m 2}) \\ P (9 | γ) = 1 \end{aligned}

Where we have:

\begin{aligned} p_{u} = \frac{μ_{u}}{μ_{u} + μ_{m}}; p_{r e g_d} = \frac{μ_{r e g_d}}{μ_{r e g_d} + μ_{r e g_m}} \\ θ_{m} = (μ_{u} + μ_{m}) t_{γ}; θ_{r e g_m} = (μ_{r e g_d} + μ_{r e g_m}) t_{γ} \\ p_{m 1} = (p_{u} + ((1 - p_{u}) e^{(- θ_{m})})); p_{m 2} = ((1 - p_{u}) + (p_{u} e^{(- θ_{m})})) \\ p_{r e g_m 1} = (p_{r e g_d} + ((1 - p_{r e g_d}) e^{(- θ_{r e g_m})})); p_{r e g_m 2} = ((1 - p_{r e g_d}) + (p_{r e g_d} e^{(- θ_{r e g_m})})) \\ t_{m} = \frac{1}{(μ_{r e g_d} + μ_{r e g_m})} \\ p_{e q 1} = (p + ((1 - p) e^{(- (μ_{u} + μ_{m}) * (t_{m}))})); p_{e q 2} = ((1 - p) + (p e^{(- (μ_{u} + μ_{m}) * (t_{m}))})) \\ P_{n o_r e g_e v e n t} = (e x p (- (μ_{r e g_d} + μ_{r e g_m}) * t_{γ})) \\ p_{m 1_s m a l l} = (p_{u} + ((1 - p_{u}) e^{(- (μ_{u} + μ_{m}) * ((\frac{1}{(μ_{r e g_d} + μ_{r e g_m})}) - (\frac{(t_{γ} e^{(- (μ_{r e g_d} + μ_{r e g_m}) t_{γ})})}{(1 - e x p (- (μ_{r e g_d} + μ_{r e g_m}) t_{γ}))})))})) \\ p_{m 2_s m a l l} = ((1 - p_{u}) + (p_{u} e^{(- (μ_{u} + μ_{m}) * ((\frac{1}{(μ_{r e g_d} + μ_{r e g_m})}) - (\frac{(t_{γ} e^{(- (μ_{r e g_d} + μ_{r e g_m}) t_{γ})})}{(1 - e x p (- (μ_{r e g_d} + μ_{r e g_m}) t_{γ}))})))})) \\ P_{O} = (p_{r e g_d}) (1 - p_{e q 1}) \\ P_{P} = (p_{r e g_d}) (p_{e q 1}) \\ P_{U} = (1 - p_{r e g_d}) (1 - p_{e q 2}) \\ P_{M} = (1 - p_{r e g_d}) (p_{e q 2}) \\ P_{O O} = p_{r e g_m 1} (((1 - P_{n o_r e g_e v e n t}) (1 - p_{m 1_s m a l l})) + (P_{n o_r e g_e v e n t} (p_{m 2}))) \\ P_{O P} = p_{r e g_m 1} (((1 - P_{n o_r e g_e v e n t}) (p_{M} 1_{s} m a l l)) + (P_{n o_r e g_e v e n t} (1 - p_{m 2}))) \\ P_{O M} = (1 - p_{r e g_m 1}) p_{m 2_s m a l l} \\ P_{O U} = (1 - p_{r e g_m 1}) (1 - p_{m 2_s m a l l}) \\ P_{P O} = p_{r e g_m 1} (((1 - P_{n o_r e g_e v e n t}) (1 - p_{m 1_s m a l l})) + (P_{n o_r e g_e v e n t} (1 - p_{m 1}))) \\ P_{P P} = p_{r e g_m 1} (((1 - P_{n o_r e g_e v e n t}) (p_{m 1_s m a l l})) + (P_{n o_r e g_e v e n t} (p_{m 1}))) \\ P_{P M} = (1 - p_{r e g_m 1}) p_{m 2_s m a l l} \\ P_{P U} = (1 - p_{r e g_m 1}) (1 - p_{m 2_s m a l l}) \\ P_{M O} = (1 - p_{r e g_m 2}) (1 - p_{m 1_s m a l l}) \\ P_{M P} = (1 - p_{r e g_m 2}) (p_{m 1_s m a l l}) \\ P_{M M} = p_{r e g_m 2} (((1 - P_{n o_r e g_e v e n t}) (p_{m 2_s m a l l})) + (P_{n o_r e g_e v e n t} (p_{m 2}))) \\ P_{M U} = p_{r e g_m 2} (((1 - P_{n o_r e g_e v e n t}) (1 - p_{m 2_s m a l l})) + (P_{n o_r e g_e v e n t} (1 - p_{m 2}))) \\ P_{U O} = (1 - p_{r e g_m 2}) (1 - p_{m 1_s m a l l}) \\ P_{U P} = (1 - p_{r e g_m 2}) (p_{m 1_s m a l l}) \\ P_{U M} = p_{r e g_m 2} (((1 - P_{n o_r e g_e v e n t}) (p_{m 2_s m a l l})) + (P_{n o_r e g_e v e n t} (1 - p_{m 1}))) \\ P_{U U} = p_{r e g_m 2} (((1 - P_{n o_r e g_e v e n t}) (1 - p_{m 2_s m a l l})) + (P_{n o_r e g_e v e n t} (p_{m 1}))) \end{aligned}

Where $μ$ is the mutation rate per nucleotide per N generation, $μ_{m}$ the site methylation rate per generation, $μ_{u}$ the site demethylation rate, $μ_{r e g_m}$ the region methylation rate per generation, $μ_{r e g_d}$ the region demethylation rate per generation and $t_{γ}$ the average coalescent time in state $γ$ . Please refer to the R function build_emi_m from the Package eSMC2 for a clearer description of the probabilities. Note that in the text, these probabilities are written as $μ_{1}$ for mutation rate per nucleotide, $μ_{S M}$ the site methylation rate, $μ_{S U}$ the site demethylation rate, $μ_{R M}$ the region methylation rate, and $μ_{R U}$ the region demethylation rate.

Data availability

The eSMC2 R package can be found at: https://github.com/TPPSellinger/eSMC2 (copy archived at Sellinger, 2024a). The input files created from Arabidopsis thaliana sequence data are available on GitHub at https://github.com/TPPSellinger/Arabidopsis_thaliana_methylation (copy archived at Sellinger, 2024b).

The following previously published data sets were used

1. Cao J
2. Schneeberger K
3. Ossowski S
4. Günther T
5. Bender S
6. Fitz J
7. Koenig D
8. Lanz C
9. Stegle O
10. Lippert C
11. Wang X
12. Ott F
13. Müller J
14. Borgwardt K
15. Schmid KJ
16. Weigel D
17. Alonso-Blanco C
(2011) Genomes
ID ProjectMPICao2010. Arabaidopsis thaliana methylomes.

https://1001genomes.org/projects/MPICao2010/index.html

References

1. Albers PK
2. McVean G
(2020) Dating genomic variants and shared ancestry in population-scale sequencing data
PLOS Biology 18:e3000586.

https://doi.org/10.1371/journal.pbio.3000586
- PubMed
- Google Scholar
1. Alonso-Blanco C
2. Andrade J
3. Becker C
4. Bemm F
5. Bergelson J
6. Borgwardt KM
7. Cao J
8. Chae E
9. Dezwaan TM
10. Ding W
11. Ecker JR
12. Exposito-Alonso M
13. Farlow A
14. Fitz J
15. Gan X
16. Grimm DG
17. Hancock AM
18. Henz SR
19. Holm S
20. Horton M
21. Jarsulic M
22. Kerstetter RA
23. Korte A
24. Korte P
25. Lanz C
26. Lee CR
27. Meng D
28. Michael TP
29. Mott R
30. Muliyati NW
31. Nägele T
32. Nagler M
33. Nizhynska V
34. Nordborg M
35. Novikova Py
36. Picó FX
37. Platzer A
38. Rabanal FA
39. Rodriguez A
40. Rowan BA
41. Salomé PA
42. Schmid KJ
43. Schmitz RJ
44. Seren Ü
45. Sperone FG
46. Sudkamp M
47. Svardal H
48. Tanzer MM
49. Todd D
50. Volchenboum SL
51. Wang C
52. Wang G
53. Wang X
54. Weckwerth W
55. Weigel D
56. Zhou X
(2016) 1,135 genomes reveal the global pattern of polymorphism in Arabidopsis thaliana
Cell 166:481–491.

https://doi.org/10.1016/j.cell.2016.05.063
- PubMed
- Google Scholar
1. Anzai T
2. Shiina T
3. Kimura N
4. Yanagiya K
5. Kohara S
6. Shigenari A
7. Yamagata T
8. Kulski JK
9. Naruse TK
10. Fujimori Y
11. Fukuzumi Y
12. Yamazaki M
13. Tashiro H
14. Iwamoto C
15. Umehara Y
16. Imanishi T
17. Meyer A
18. Ikeo K
19. Gojobori T
20. Bahram S
21. Inoko H
(2003) Comparative sequencing of human and chimpanzee MHC class I regions unveils insertions/deletions as the major path to genomic divergence
PNAS 100:7708–7713.

https://doi.org/10.1073/pnas.1230533100
- PubMed
- Google Scholar
(2019) Inference of recombination maps from a single pair of genomes and its application to ancient samples
PLOS Genetics 15:e1008449.

https://doi.org/10.1371/journal.pgen.1008449
- PubMed
- Google Scholar
1. Barroso GV
2. Dutheil JY
(2023) The landscape of nucleotide diversity in Drosophila melanogaster is shaped by mutation rate variation
Peer Community Journal 3:267.

https://doi.org/10.24072/pcjournal.267
- Google Scholar
1. Baumdicker F
2. Bisschop G
3. Goldstein D
4. Gower G
5. Ragsdale AP
6. Tsambos G
7. Zhu S
8. Eldon B
9. Ellerman EC
10. Galloway JG
11. Gladstein AL
12. Gorjanc G
13. Guo B
14. Jeffery B
15. Kretzschumar WW
16. Lohse K
17. Matschiner M
18. Nelson D
19. Pope NS
20. Quinto-Cortés CD
21. Rodrigues MF
22. Saunack K
23. Sellinger T
24. Thornton K
25. van Kemenade H
26. Wohns AW
27. Wong Y
28. Gravel S
29. Kern AD
30. Koskela J
31. Ralph PL
32. Kelleher J
(2022) Efficient ancestry and mutation simulation with msprime 1.0
Genetics 220:iyab229.

https://doi.org/10.1093/genetics/iyab229
- PubMed
- Google Scholar
(2021) Sweeps in time: leveraging the joint distribution of branch lengths
Genetics 219:iyab119.

https://doi.org/10.1093/genetics/iyab119
- PubMed
- Google Scholar
1. Boitard S
2. Rodríguez W
3. Jay F
4. Mona S
5. Austerlitz F
(2016) Inferring Population size history from large samples of genome-wide molecular data - an approximate bayesian computation approach
PLOS Genetics 12:e1005877.

https://doi.org/10.1371/journal.pgen.1005877
- PubMed
- Google Scholar
1. Brandt, DY
(2022) Evaluation of methods for the inference of ancestral recombination graphs
Genetics 221:iyac044.

https://doi.org/10.1093/genetics/iyac044
- Google Scholar
1. Briffa A
2. Hollwey E
3. Shahzad Z
4. Moore JD
5. Lyons DB
6. Howard M
7. Zilberman D
(2023) Millennia-long epigenetic fluctuations generate intragenic DNA methylation variance in Arabidopsis populations
Cell Systems 14:953–967.

https://doi.org/10.1016/j.cels.2023.10.007
- PubMed
- Google Scholar
1. Cao J
2. Schneeberger K
3. Ossowski S
4. Günther T
5. Bender S
6. Fitz J
7. Koenig D
8. Lanz C
9. Stegle O
10. Lippert C
11. Wang X
12. Ott F
13. Müller J
14. Alonso-Blanco C
15. Borgwardt K
16. Schmid KJ
17. Weigel D
(2011) Whole-genome sequencing of multiple Arabidopsis thaliana populations
Nature Genetics 43:956–963.

https://doi.org/10.1038/ng.911
- PubMed
- Google Scholar
Book
1. Charlesworth B
2. Charlesworth D
(2010)
Elements of Evolutionary Genetics

Roberts and Company.
- Google Scholar
1. Charlesworth B
2. Jain K
(2014) Purifying selection, drift, and reversible mutation with arbitrarily high mutation rates
Genetics 198:1587–1602.

https://doi.org/10.1534/genetics.114.167973
- PubMed
- Google Scholar
1. Charlesworth B
2. Jensen JD
(2023) Population genetic considerations regarding evidence for biased mutation rates in Arabidopsis thaliana
Molecular Biology and Evolution 40:msac275.

https://doi.org/10.1093/molbev/msac275
- PubMed
- Google Scholar
1. Cokus SJ
2. Feng S
3. Zhang X
4. Chen Z
5. Merriman B
6. Haudenschild CD
7. Pradhan S
8. Nelson SF
9. Pellegrini M
10. Jacobsen SE
(2008) Shotgun bisulphite sequencing of the Arabidopsis genome reveals DNA methylation patterning
Nature 452:215–219.

https://doi.org/10.1038/nature06745
- PubMed
- Google Scholar
1. Deng Y
2. Song YS
3. Nielsen R
(2021) The distribution of waiting distances in ancestral recombination graphs
Theoretical Population Biology 141:34–43.

https://doi.org/10.1016/j.tpb.2021.06.003
- PubMed
- Google Scholar
(2021) Region-level epimutation rates in Arabidopsis thaliana
Heredity 127:190–202.

https://doi.org/10.1038/s41437-021-00441-w
- PubMed
- Google Scholar
(2017) African genomes illuminate the early history and transition to selfing in Arabidopsis thaliana
PNAS 114:5213–5218.

https://doi.org/10.1073/pnas.1616736114
- PubMed
- Google Scholar
(2002) Homoplasy and mutation model at microsatellite loci and their consequences for population genetics analysis
Molecular Ecology 11:1591–1604.

https://doi.org/10.1046/j.1365-294X.2002.01576.x
- Google Scholar
(2008) Demographic history of european populations of Arabidopsis thaliana
PLOS Genetics 4:e1000075.

https://doi.org/10.1371/journal.pgen.1000075
- PubMed
- Google Scholar
(2018) Madeiran Arabidopsis thaliana reveals ancient long-range colonization and clarifies demography in eurasia
Molecular Biology and Evolution 35:564–574.

https://doi.org/10.1093/molbev/msx300
- PubMed
- Google Scholar
(2013) Inferring population size changes with sequence and SNP data: lessons from human bottlenecks
Heredity 110:409–419.

https://doi.org/10.1038/hdy.2012.120
- PubMed
- Google Scholar
(2016) Inferring past effective population size from distributions of coalescent times
Genetics 204:1191–1206.

https://doi.org/10.1534/genetics.115.185058
- PubMed
- Google Scholar
1. Hazarika RR
2. Serra M
3. Zhang Z
4. Zhang Y
5. Schmitz RJ
6. Johannes F
(2022) Molecular properties of epimutation hotspots
Nature Plants 8:146–156.

https://doi.org/10.1038/s41477-021-01086-7
- PubMed
- Google Scholar
(2020) Mapping gene flow between ancient hominins through demography-aware inference of the ancestral recombination graph
PLOS GENETICS 16:e1008895.

https://doi.org/10.1371/journal.pgen.1008895
- PubMed
- Google Scholar
1. Hudson RR
(1983) Properties of a neutral allele model with intragenic recombination
Theoretical Population Biology 23:183–201.

https://doi.org/10.1016/0040-5809(83)90013-8
- Google Scholar
1. Johannes F
(2019a) DNA methylation makes mutational history
Nature Plants 5:772–773.

https://doi.org/10.1038/s41477-019-0491-z
- Google Scholar
1. Johannes F
2. Schmitz RJ
(2019b) Spontaneous epimutations in plants
The New Phytologist 221:1253–1259.

https://doi.org/10.1111/nph.15434
- PubMed
- Google Scholar
(2020) Toward an evolutionarily appropriate null model: Jointly inferring demography and purifying selection
GENETICS 215:173–192.

https://doi.org/10.1534/genetics.119.303002
- PubMed
- Google Scholar
1. Johri P
2. Riall K
3. Becher H
4. Excoffier L
5. Charlesworth B
6. Jensen JD
(2021) The impact of purifying and background selection on the inference of population history: Problems and prospects
MOLECULAR BIOLOGY AND EVOLUTION 38:2986–3003.

https://doi.org/10.1093/molbev/msab050
- PubMed
- Google Scholar
1. Johri P
2. Aquadro CF
3. Beaumont M
4. Charlesworth B
5. Excoffier L
6. Eyre-Walker A
7. Keightley PD
8. Lynch M
9. McVean G
10. Payseur BA
11. Pfeifer SP
12. Stephan W
13. Jensen JD
(2022) Recommendations for improving statistical inference in population genomics
PLOS Biology 20:e3001669.

https://doi.org/10.1371/journal.pbio.3001669
- PubMed
- Google Scholar
(2016) Efficient coalescent simulation and genealogical analysis for large sample sizes
PLOS Computational Biology 12:e1004842.

https://doi.org/10.1371/journal.pcbi.1004842
- PubMed
- Google Scholar
1. Kelleher J
2. Wong Y
3. Wohns AW
4. Fadil C
5. Albers PK
6. McVean G
(2019) Inferring whole-genome histories in large population datasets
Nature Genetics 51:1330–1338.

https://doi.org/10.1038/s41588-019-0483-y
- PubMed
- Google Scholar
1. Ki C
2. Terhorst J
(2023) Exact decoding of a sequentially markov coalescent model in genetics
Journal of the American Statistical Association 1:1–14.

https://doi.org/10.1080/01621459.2023.2252570
- Google Scholar
1. Kingman JFC
(1982) The coalescent
Stochastic Processes and Their Applications 13:235–248.

https://doi.org/10.1016/0304-4149(82)90011-4
- Google Scholar
(2023) Deep learning in population genetics
Genome Biology and Evolution 15:evad008.

https://doi.org/10.1093/gbe/evad008
- PubMed
- Google Scholar
(2024) Simultaneous inference of past demography and selection from the ancestral recombination graph under the beta coalescent
Peer Community Journal 4:397.

https://doi.org/10.24072/pcjournal.397
- Google Scholar
1. Lang D
2. Zhang S
3. Ren P
4. Liang F
5. Sun Z
6. Meng G
7. Tan Y
8. Li X
9. Lai Q
10. Han L
11. Wang D
12. Hu F
13. Wang W
14. Liu S
(2020) Comparison of the two up-to-date sequencing technologies for genome assembly: HiFi reads of pacific biosciences sequel II system and ultralong reads of oxford nanopore
GigaScience 9:giaa123.

https://doi.org/10.1093/gigascience/giaa123
- PubMed
- Google Scholar
1. Li H
2. Durbin R
(2011) Inference of human population history from individual whole-genome sequences
Nature 475:493–496.

https://doi.org/10.1038/nature10231
- PubMed
- Google Scholar
(2008) Highly integrated single-base resolution maps of the epigenome in Arabidopsis
Cell 133:523–536.

https://doi.org/10.1016/j.cell.2008.03.029
- PubMed
- Google Scholar
1. Lyons DB
2. Briffa A
3. He S
4. Choi J
5. Hollwey E
6. Colicchio J
7. Anderson I
8. Feng X
9. Howard M
10. Zilberman D
(2023) Extensive de novo activity stabilizes epigenetic inheritance of CG methylation in Arabidopsis transposons
Cell Reports 42:112132.

https://doi.org/10.1016/j.celrep.2023.112132
- PubMed
- Google Scholar
1. Mahmoudi A
2. Koskela J
3. Kelleher J
4. Chan YB
5. Balding D
(2022) Bayesian inference of ancestral recombination graphs
PLOS Computational Biology 18:e1009960.

https://doi.org/10.1371/journal.pcbi.1009960
- PubMed
- Google Scholar
1. Malaspinas AS
2. Westaway MC
3. Muller C
4. Sousa VC
5. Lao O
6. Alves I
7. Bergström A
8. Athanasiadis G
9. Cheng JY
10. Crawford JE
11. Heupink TH
12. Macholdt E
13. Peischl S
14. Rasmussen S
15. Schiffels S
16. Subramanian S
17. Wright JL
18. Albrechtsen A
19. Barbieri C
20. Dupanloup I
21. Eriksson A
22. Margaryan A
23. Moltke I
24. Pugach I
25. Korneliussen TS
26. Levkivskyi IP
27. Moreno-Mayar JV
28. Ni S
29. Racimo F
30. Sikora M
31. Xue Y
32. Aghakhanian FA
33. Brucato N
34. Brunak S
35. Campos PF
36. Clark W
37. Ellingvåg S
38. Fourmile G
39. Gerbault P
40. Injie D
41. Koki G
42. Leavesley M
43. Logan B
44. Lynch A
45. Matisoo-Smith EA
46. McAllister PJ
47. Mentzer AJ
48. Metspalu M
49. Migliano AB
50. Murgha L
51. Phipps ME
52. Pomat W
53. Reynolds D
54. Ricaut FX
55. Siba P
56. Thomas MG
57. Wales T
58. Wall CM
59. Oppenheimer SJ
60. Tyler-Smith C
61. Durbin R
62. Dortch J
63. Manica A
64. Schierup MH
65. Foley RA
66. Lahr MM
67. Bowern C
68. Wall JD
69. Mailund T
70. Stoneking M
71. Nielsen R
72. Sandhu MS
73. Excoffier L
74. Lambert DM
75. Willerslev E
(2016) A genomic history of Aboriginal Australia
Nature 538:207–214.

https://doi.org/10.1038/nature18299
- PubMed
- Google Scholar
1. Marjoram P
2. Wall JD
(2006) Fast “coalescent” simulation
BMC Genetics 7:16.

https://doi.org/10.1186/1471-2156-7-16
- PubMed
- Google Scholar
1. McVean GAT
2. Cardin NJ
(2005) Approximating the coalescent with recombination
Philosophical Transactions of the Royal Society of London. Series B, Biological Sciences 360:1387–1393.

https://doi.org/10.1098/rstb.2005.1673
- PubMed
- Google Scholar
1. Monroe JG
2. Srikant T
3. Carbonell-Bejerano P
4. Becker C
5. Lensink M
6. Exposito-Alonso M
7. Klein M
8. Hildebrandt J
9. Neumann M
10. Kliebenstein D
11. Weng ML
12. Imbert E
13. Ågren J
14. Rutter MT
15. Fenster CB
16. Weigel D
(2022) Mutation bias reflects natural selection in Arabidopsis thaliana
Nature 602:101–105.

https://doi.org/10.1038/s41586-021-04269-6
- PubMed
- Google Scholar
(2021) Gene body methylation is under selection in Arabidopsis thaliana
Genetics 218:iyab061.

https://doi.org/10.1093/genetics/iyab061
- PubMed
- Google Scholar
1. Nordborg M
(2000) Linkage disequilibrium, gene trees and selfing: an ancestral recombination graph with partial self-fertilization
Genetics 154:923–929.

https://doi.org/10.1093/genetics/154.2.923
- PubMed
- Google Scholar
1. Nurk S
2. Walenz BP
3. Rhie A
4. Vollger MR
5. Logsdon GA
6. Grothe R
7. Miga KH
8. Eichler EE
9. Phillippy AM
10. Koren S
(2020) HiCanu: accurate assembly of segmental duplications, satellites, and allelic variants from high-fidelity long reads
Genome Research 30:1291–1305.

https://doi.org/10.1101/gr.263566.120
- PubMed
- Google Scholar
(2010) The rate and molecular spectrum of spontaneous mutations in Arabidopsis thaliana
Science 327:92–94.

https://doi.org/10.1126/science.1180677
- PubMed
- Google Scholar
1. Ou S
2. Su W
3. Liao Y
4. Chougule K
5. Agda JRA
6. Hellinga AJ
7. Lugo CSB
8. Elliott TA
9. Ware D
10. Peterson T
11. Jiang N
12. Hirsch CN
13. Hufford MB
(2019) Benchmarking transposable element annotation methods for creation of a streamlined, comprehensive pipeline
Genome Biology 20:275.

https://doi.org/10.1186/s13059-019-1905-y
- PubMed
- Google Scholar
(2023) On the causes of gene-body methylation variation in Arabidopsis thaliana
PLOS Genetics 19:e1010728.

https://doi.org/10.1371/journal.pgen.1010728
- PubMed
- Google Scholar
1. Rodríguez W
2. Mazet O
3. Grusea S
4. Arredondo A
5. Corujo JM
6. Boitard S
7. Chikhi L
(2018) The IICR and the non-stationary structured coalescent: towards demographic inference with arbitrary changes in population structure
Heredity 121:663–678.

https://doi.org/10.1038/s41437-018-0148-0
- PubMed
- Google Scholar
1. Salomé PA
2. Bomblies K
3. Fitz J
4. Laitinen RAE
5. Warthmann N
6. Yant L
7. Weigel D
(2012) The recombination landscape in Arabidopsis thaliana F2 populations
Heredity 108:447–455.

https://doi.org/10.1038/hdy.2011.95
- PubMed
- Google Scholar
(2013) zipHMMlib: a highly optimised HMM library exploiting repetitions in the input to speed up the forward algorithm
BMC Bioinformatics 14:339.

https://doi.org/10.1186/1471-2105-14-339
- PubMed
- Google Scholar
1. Schiffels S
2. Durbin R
(2014) Inferring human population size and separation history from multiple genome sequences
Nature Genetics 46:919–925.

https://doi.org/10.1038/ng.3015
- PubMed
- Google Scholar
(2018) Contribution of epigenetic variation to adaptation in Arabidopsis
Nature Communications 9:4446.

https://doi.org/10.1038/s41467-018-06932-5
- PubMed
- Google Scholar
1. Schmitz RJ
2. Schultz MD
3. Urich MA
4. Nery JR
5. Pelizzola M
6. Libiger O
7. Alix A
8. McCosh RB
9. Chen H
10. Schork NJ
11. Ecker JR
(2013) Patterns of population epigenomic diversity
Nature 495:193–198.

https://doi.org/10.1038/nature11968
- PubMed
- Google Scholar
1. Schraiber JG
2. Akey JM
(2015) Methods and models for unravelling human evolutionary history
Nature Reviews. Genetics 16:727–740.

https://doi.org/10.1038/nrg4005
- PubMed
- Google Scholar
1. Schweiger R
2. Durbin R
(2023) Ultrafast genome-wide inference of pairwise coalescence times
Genome Research 33:1023–1031.

https://doi.org/10.1101/gr.277665.123
- PubMed
- Google Scholar
(2020) Inference of past demography, dormancy and self-fertilization rates from whole genome sequence data
PLOS Genetics 16:e1008698.

https://doi.org/10.1371/journal.pgen.1008698
- PubMed
- Google Scholar
(2021) Limits and convergence properties of the sequentially markovian coalescent
Molecular Ecology Resources 21:2231–2248.

https://doi.org/10.1111/1755-0998.13416
- PubMed
- Google Scholar
Software
1. Sellinger TPP
(2024a) eSMC2, version swh:1:rev:66e60b339556a1b9faa63126e0237288770612e9
Software Heritage.

https://archive.softwareheritage.org/swh:1:dir:08f86e3a805a5d8e24c407c8f6fa3043e437453b;origin=https://github.com/TPPSellinger/eSMC2;visit=swh:1:snp:a18ca35636ebb0c123f3516b38023effdf6de73e;anchor=swh:1:rev:66e60b339556a1b9faa63126e0237288770612e9
Software
1. Sellinger TPP
(2024b) Arabidopsis_thaliana_methylation, version swh:1:rev:324a6df6b612c48a3cd5eb919f2fefb212dd5a06
Software Heritage.

https://archive.softwareheritage.org/swh:1:dir:22ffaf5c4a1fc683e669cbf18aac5bf59d3fed7a;origin=https://github.com/TPPSellinger/Arabidopsis_thaliana_methylation;visit=swh:1:snp:148971ffe4fd957e87d2f408bb6c789d5fc4f356;anchor=swh:1:rev:324a6df6b612c48a3cd5eb919f2fefb212dd5a06
(2020) AlphaBeta: computational inference of epimutation rates and spectra from high-throughput DNA methylation data in plants
Genome Biology 21:260.

https://doi.org/10.1186/s13059-020-02161-6
- PubMed
- Google Scholar
1. Speidel L
2. Forest M
3. Shi S
4. Myers SR
(2019) A method for genome-wide genealogy estimation for thousands of samples
Nature Genetics 51:1321–1329.

https://doi.org/10.1038/s41588-019-0484-x
- PubMed
- Google Scholar
1. Srikant T
2. Drost HG
(2020) How stress facilitates phenotypic innovation through epigenetic diversity
Frontiers in Plant Science 11:606800.

https://doi.org/10.3389/fpls.2020.606800
- PubMed
- Google Scholar
(2023) Joint inference of evolutionary transitions to self-fertilization and demographic history using whole-genome sequences
eLife 12:e82384.

https://doi.org/10.7554/eLife.82384
- PubMed
- Google Scholar
(2018) METHimpute: imputation-guided construction of complete methylomes from WGBS data
BMC Genomics 19:444.

https://doi.org/10.1186/s12864-018-4641-x
- PubMed
- Google Scholar
(2011) Inference of seed bank parameters in two wild tomato species using ecological and genetic data
PNAS 108:17052–17057.

https://doi.org/10.1073/pnas.1111266108
- PubMed
- Google Scholar
(2017) Robust and scalable inference of population history from hundreds of unphased whole genomes
Nature Genetics 49:303–309.

https://doi.org/10.1038/ng.3748
- PubMed
- Google Scholar
1. Upadhya G
2. Steinrücken M
(2022) Robust inference of population size histories from genomic sequencing data
PLOS Computational Biology 18:e1010419.

https://doi.org/10.1371/journal.pcbi.1010419
- PubMed
- Google Scholar
(2015) Rate, spectrum, and evolutionary dynamics of spontaneous epimutations
PNAS 112:6676–6681.

https://doi.org/10.1073/pnas.1424254112
- PubMed
- Google Scholar
(2016) Methylome evolution in plants
Genome Biology 17:264.

https://doi.org/10.1186/s13059-016-1127-5
- PubMed
- Google Scholar
Book
1. Wakeley J
(2008)
Coalescent Theory: An Introduction. Roberts and Company

Greenwood Village, USA, CO.
- Google Scholar
1. Wang J
2. Fan C
(2014) A neutrality test for detecting selection on DNA methylation using single methylation polymorphism frequency spectrum
Genome Biology and Evolution 7:154–171.

https://doi.org/10.1093/gbe/evu271
- PubMed
- Google Scholar
1. Wang C
(2018) MSIpred: a python package for tumor microsatellite instability classification from tumor mutation annotation data using a support vector machine
Scientific Reports 8:17546.

https://doi.org/10.1038/s41598-018-35682-z
- PubMed
- Google Scholar
1. Weigel D
2. Colot V
(2012) Epialleles in plant evolution
Genome Biology 13:1–6.

https://doi.org/10.1186/gb-2012-13-10-249
- PubMed
- Google Scholar
1. Wiuf C
2. Hein J
(1999) Recombination as a point process along sequences
Theoretical Population Biology 55:248–259.

https://doi.org/10.1006/tpbi.1998.1403
- Google Scholar
1. Wohns AW
2. Wong Y
3. Jeffery B
4. Akbari A
5. Mallick S
6. Pinhasi R
7. Patterson N
8. Reich D
9. Kelleher J
10. McVean G
(2022) A unified genealogy of modern and ancient genomes
Science 375:eabi8264.

https://doi.org/10.1126/science.abi8264
- PubMed
- Google Scholar
1. Yang Z
(1996) Statistical properties of a DNA sample under the finite-sites model
Genetics 144:1941–1950.

https://doi.org/10.1093/genetics/144.4.1941
- PubMed
- Google Scholar
(2018) Indel detection from DNA and RNA sequencing data with transIndel
BMC Genomics 19:270.

https://doi.org/10.1186/s12864-018-4671-4
- PubMed
- Google Scholar
(2021) Epimutations define a fast-ticking molecular clock in plants
Trends in Genetics 37:699–710.

https://doi.org/10.1016/j.tig.2021.04.010
- PubMed
- Google Scholar
1. Yao N
2. Zhang Z
3. Yu L
4. Hazarika R
5. Yu C
6. Jang H
7. Smith LM
8. Ton J
9. Liu L
10. Stachowicz JJ
11. Reusch TBH
12. Schmitz RJ
13. Johannes F
(2023) An evolutionary epigenetic clock in plants
Science 381:1440–1445.

https://doi.org/10.1126/science.adh9443
- PubMed
- Google Scholar
1. Zhang X
2. Yazaki J
3. Sundaresan A
4. Cokus S
5. Chan SW-L
6. Chen H
7. Henderson IR
8. Shinn P
9. Pellegrini M
10. Jacobsen SE
11. Ecker JR
(2006) Genome-wide high-resolution mapping and functional analysis of DNA Methylation in Arabidopsis
Cell 126:1189–1201.

https://doi.org/10.1016/j.cell.2006.08.003
- Google Scholar
1. Zhang Y
2. Wang S
3. Wang X
(2018) Data-driven-based approach to identifying differentially methylated regions using modified 1D ising model
BioMed Research International 2018:1070645.

https://doi.org/10.1155/2018/1070645
- PubMed
- Google Scholar
(2007) Genome-wide analysis of Arabidopsis thaliana DNA methylation uncovers an interdependence between methylation and transcription
Nature Genetics 39:61–69.

https://doi.org/10.1038/ng1929
- Google Scholar

Article and author information

Author details

Thibaut Sellinger
1. Professorship for Population Genetics, Department of Life Science Systems, Technical University of Munich, Munich, Germany
2. Department of Environment and Biodiversity, Paris Lodron University of Salzburg, Salzburg, Austria
Contribution
Conceptualization, Software, Formal analysis, Investigation, Visualization, Methodology, Writing - original draft

Competing interests
No competing interests declared
Frank Johannes

Professorship for Plant Epigenomics, Department of Molecular Life Sciences, Technical University of Munich, Freising, Germany

Contribution
Conceptualization, Writing - review and editing

Competing interests
No competing interests declared

"This ORCID iD identifies the author of this article:" 0000-0002-7962-2907
Aurélien Tellier

Professorship for Population Genetics, Department of Life Science Systems, Technical University of Munich, Munich, Germany

Contribution
Conceptualization, Supervision, Funding acquisition, Writing - original draft, Project administration, Writing - review and editing

For correspondence
aurelien.tellier@tum.de

Competing interests
No competing interests declared

"This ORCID iD identifies the author of this article:" 0000-0002-8895-0785

Funding

Deutsche Forschungsgemeinschaft (317616126 (TE809/7-1))

Aurélien Tellier

Austrian Science Fund (project no. TAI 151-B)

Thibaut Sellinger

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

Acknowledgements

We thank Zhilin Zhang and Rashmi Hazarika for giving and processing the data of Arabidopsis thaliana. TS is supported by the Deutsche Forschungsgemeinschaft, project number 317616126 (TE809/7-1) to AT, and the Austrian Science Fund (project no. TAI 151-B) to Anja Hörger.

Version history

Preprint posted: May 16, 2023
Sent for peer review: May 19, 2023
Reviewed Preprint version 1: August 24, 2023
Reviewed Preprint version 2: February 27, 2024
Reviewed Preprint version 3: August 5, 2024
Version of Record published: September 12, 2024

Cite all versions

You can cite all versions using the DOI https://doi.org/10.7554/eLife.89470. This DOI represents all versions, and will always resolve to the latest one.

Copyright

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.