Improved inference of population histories by integrating genomic and epigenomic data

  1. Thibaut Sellinger
  2. Frank Johannes
  3. Aurélien Tellier

  • Curated by eLife

    eLife logo

    eLife assessment

    This important study extends existing sequentially Markovian coalescent approaches to include the use of hypervariable loci such as epimutations. This is an intriguing addition 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 -- and thus the potential use of this approach -- there would be additional value in more extensive consideration of when and where these methods are best used.

Reviewed by

Read the full article See related articles

Listed in

Log in to save this article

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, e.g. 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 A. thaliana using DNA methylation diversity data from 10 European natural accessions. We demonstrate that segregating Single Methylated Polymorphisms (SMPs) satisfy the modelling 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.

Article activity feed

  1. Version published to 10.7554/elife.89470.2 on eLife
  2. Version published to 10.7554/elife.89470 on eLife
  3. eLife

    Author Response

    The following is the authors’ response to the original reviews.

    First, we would like to thank you and all the reviewers for acknowledging the meaningful contribution of our manuscript to the field. Your useful comments helped us improve the manuscript's quality. We understood the key issues of the manuscript were the quantification of inference accuracy and applicability to methylome data. We here therefore present a revised version of the manuscript addressing all major comments.

    For each demographic inference we have added the root mean square error as demanded by the reviewers. These results confirm the previous interpretation of the graphs especially in recent times. We also added TMRCA inference analysis as requested by one reviewer as a proof of principle that integrating multiple markers can improve ARG inference.

    The discussion was rewritten to further discuss the challenges of application to empirical methylation data. We clarify that in the case epimutations are well understood and modelled, they can be integrated into a SMC framework to improve the approaches accuracy. When epimutations are not well understood, our approach can help understand the epimutations process through generations at the evolutionary time scale along the genome. Hence, in both cases our approach can be used to unveil marker evolution processes through generations, and/or deepen our understanding of the population past history. We hope our discussion underlies better how our approach is designed and can be used.

    eLife assessment

    This important study advances existing approaches for demographic inference by incorporating rapidly mutating markers such as switches in methylation state. The authors provide a solid comparison of their approach to existing methods, although the work would benefit from some additional consideration of the challenges in the empirical use of methylation data. The work will be of broad interest to population geneticists, both in terms of the novel approach and the statistical inference proposed.

    Public Reviews:

    Reviewer #1 (Public Review):

    The authors developed an extension to the pairwise sequentially Markov coalecent model that allows to simultaneously analyse multiple types of polymorphism data. In this paper, they focus on SNPs and DNA methylation data. Since methylation markers mutate at a much faster rate than SNPs, this potentially gives the method better power to infer size history in the recent past. Additionally, they explored a model where there are both local and regional epimutational processes.

    Integrating additional types of heritable markers into SMC is a nice idea which I like in principle. However, a major caveat to this approach seems to be a strong dependence on knowing the epimutation rate. In Fig. 6 it is seen that, when the epimutation rate is known, inferences do indeed look better; but this is not necessarily true when the rate is not known. A roughly similar pattern emerges in Supp. Figs. 4-7; in general, results when the rates have to be estimated don't seem that much better than when focusing on SNPs alone. This carries over to the real data analysis too: the interpretation in Fig. 7 appears to hinge on whether the rates are known or estimated, and the estimated rates differ by a large amount from earlier published ones.

    Overall, this is an interesting research direction, and I think the method may hold more promise as we get more and better epigenetic data, and in particular better knowledge of the epigenetic mutational process. At the same time, I would be careful about placing too much emphasis on new findings that emerge solely by switching to SNP+SMP analysis.

    Answer: We thank the reviewer 1 for his positive comments and acknowledging the future promises of our method as better and more reliable data will be available in different species. We appreciate the reviewer noticing the complete set of work undertaken here to integrate local and regional effects of methylation into a model containing as much knowledge of the epigenetics mutational processes as possible. Note that in Figure 2 of the manuscript we observed a gain of accuracy even when the rates are unknown. Our results thus suggests that the accuracy gain of additional marker with unknown rates is also possible, although it is most likely be scenario and rate dependent.

    At last, as noticed and highlighted by the very recent work of the Johannes lab (Yao et al. Science 2023) using phylogenetic methods, knowing the epimutation rate is essential at short time scale to avoid confounding effects of homoplasy. In our estimation of the coalescent trees, the same applies, though our model considers finite site markers. We now provide additional evidence for the potential gain of power to infer the TMRCA (Supplementary Table S7) when knowing or not the epimutation rates and revised the discussion to clarify the potential shortcomings/caveats for the analysis of real data.

    Reviewer #2 (Public Review):

    A limitation in using SNPs to understand recent histories of genomes is their low mutation frequency. Tellier et al. explore the possibility of adding hypermutable markers to SNP based methods for better resolution over short time frames. In particular, they hypothesize that epimutations (CG methylation and demethylation) could provide a useful marker for this purpose. Individual CGs in Arabidopsis tends to be either close to 100% methylated or close to 0%, and are inherited stably enough across generations that they can be treated as genetic markers. Small regions containing multiple CGs can also be treated as genetic markers based on their cumulative methylation level. In this manuscript, Tellier et al develop computational methods to use CG methylation as a hypermutable genetic marker and test them on theoretical and real data sets. They do this both for individual CGs and small regions. My review is limited to the simple question of whether using CG methylation for this purpose makes sense at a conceptual level, not at the level of evaluating specific details of the methods. I have a small concern in that it is not clear that CG methylation measurements are nearly as binary in other plants and other eukaryotes as they are in Arabidopsis. However, I see no reason why the concept of this work is not conceptually sound. Especially in the future as new sequencing technologies provide both base calling and methylating calling capabilities, using CG methylation in addition to SNPs could become a useful and feasible tool for population genetics in situations where SNPs are insufficient.

    Answer: We thank the reviewer 2 for his positive comments. Indeed, surveys of CG methylation in other plant species show that its distribution is clearly bimodal (i.e. binary). This is not the case for non-CG methylation, such as CHG and CHH (where H=C,T,A). However, these later types of methylation contexts are also not heritable across generations and can therefore not be used as heritable molecular markers.

    Reviewer #3 (Public Review):

    I very much like this approach and the idea of incorporating hypervariable markers. The method is intriguing, and the ability to e.g. estimate recombination rates, the size of DMRs, etc. is a really nice plus. I am not able to comment on the details of the statistical inference, but from what I can evaluate it seems sound and reasonable. This is an exciting new avenue for thinking about inference from genomic data. I have a few concerns about the presentation and then also questions about the use of empirical methylation data sets.

    I think a more detailed description of demographic accuracy is warranted. For example, in L245 MSMC2 identifies the bottleneck (albeit smoothed) and only slightly overestimates recent size. In the same analysis the authors' approach with unknown mu infers a nonexistent population increase by an order of magnitude that is not mentioned.

    Answer: We thank the reviewer 3 for his positive comments and refer to our answer to reviewer 1 above. We added RMSE (Root Mean Square Error) analyses to quantify the inference accuracy. We apologize for not mentioning this last point. Thank you for pointing this out and we have now fixed it (line 245-253).

    Similarly, it seems problematic that (L556) the approach requiring estimation of site and region parameters (as would presumably be needed in most empirical systems like endangered nonmodel species mentioned in the introduction) does no better than using only SNPs. Overall, I think a more objective and perhaps quantitative comparison of approaches is warranted.

    Answer : See answer to reviewer 1 above, and more elaborate answers below. We provide now new RMSE analyses to quantify the accuracy of our demographic inference (Supplementary Tables 1,6,7,8,9,10). We also discuss the validity and usefulness of our approach when the epimutation rates are unknown. In short, the discussion was rewritten to further discuss the challenges of application to empirical methylation data. We clarify that in the case epimutations are well known and modelled (as much is known in A. thaliana for example), they can be integrated into a SMC framework to improve the accuracy of the method approach. When epimutations are not well understood and rates unknown, our approach can help understand the epimutational process through generations at the evolutionary time scale. Hence, whether makers are understood or not, our approach can be used to study the marker evolutionary processes through generations and/or to deepen our understanding of the population past history. We hope our discussion underlies better how our approach is designed and can be used.

    The authors simulate methylated markers at 2% (and in some places up to 20%). In many plant genomes a large proportion of cytosines are methylated (e.g. 70% in maize: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC8496265/). I don't know what % of these may be polymorphic, but this leads to an order of magnitude more methylated cytosines than there are SNPs. Couldn't this mean that any appreciable error in estimating methylation threatens to be of a similar order of magnitude to the SNP data? I would welcome the authors' thoughts here.

    Answer : The reviewer is correct and this is an interesting question. First, studies show that heritable epimutations in plants are restricted to CG dinucleotides that are located well outside of the target regions of de novo methylation pathways in plants. Most of these CGs tend of fall within so-called gene body methylated regions. While it is true that plant species can differ substantially in their proportion of methylation at the genome-wide scale, the number of gene body methylated genes (i.e. genic CG methylation) is relatively similar, and at least well within the same order of magnitude (Takuno et al. Nature Plants 2016, review in Muyle et al. Genome Biol Evol 2022). Moreover, spontaneous CG epimutations in gene body methylated regions has been shown to be neutral (van Der Graaf et al. 2015, Vidali et al. 2016, Yao et al. 2023), which is an ideal property for phylogentic and demographic inference.

    Second, CG methylation calls are sometimes affected by coverage or uncertainty. Stringent filtering for reliable SMP calls typically reduces the total proportion of CG sites that can be used as input for demographic inference. Here we only kept CG sites where the methylation information could be fully trusted after SMP calling (i.e. >99.9% posteriori certainty). Overall, this explains why the percentage of sites with methylation information is so small, and why we have decided to work on simulation with 2% of reliable methylated markers.

    Nevertheless, for the sake of generality, it may be that in some species such as maize a higher percentage of polymorphic methylated sites can be used, and the number of SMPs could be higher than that of SNPs when the effective population size is very small (due to past demographic history and/or life history traits). In this case, any error in the epimutation rate and variance due to the finite site model estimation (and homoplasy) are not corrected by the lack of SNPs and can lead to mis-inference.

    A few points of discussion about the biology of methylation might be worth including. For example, methylation can differ among cell types or cells within a tissue, yet sequencing approaches evaluate a pool of cells. This results in a reasonable fraction of sites having methylation rates not clearly 0 or 1. How does this variation affect the method? Similarly, while the authors cite literature about the stable inheritance of methylation, a sentence or so more about the time scale over which this occurs would be helpful.

    Answer: We thank reviewer 3 for asking those very interesting questions, which we further developed below and mention in the discussion (lines 716-722).

    For Arabidopsis thaliana:

    Following up on our previous comment above, the majority of the CG sites that serve as input to our approach are located in body methylated genes. Previous work has shown that CG methylation in these regions shows essentially no tissue and cellular heterogeneity (e.g. Horvath et al. 2019). This means that bulk methylation measurements only show limited susceptibility to measurement error. That said, to guard against any spurious SMPs call that could arise from residual measurement variation, we applied stringent filtering of CG methylation. We have kept sites where the methylation percentage is close to either 0% or 100% (the rest being removed from the analysis). We have used similar filtering strategies in previous studies of epimutational processes in mutation accumulation lines and long-lived perennials (work of the Johannes lab). In these later studies we found that the SMP calls sufficiently accurate for inferences of phylogenetic parameters in experimental settings (Sharyhary et al. Genome Biology 2021, Yao et al. Science, 2023).

    For other species:

    It is true that currently, evaluating the methylation state of a site from a pool of cells may be problematic for some species for two main reasons: 1) it will add noise to the signal and SMP calling could be erroneous, and 2) the methylation state used in analysis might originate from different tissues at different location of the genome/methylome. Overall, this will lead to spurious SMPs and can render the inference inaccurate (see Sellinger et al 2021 for the effect of spurious SNPs). Hence, caution is advised when calling SMPs in other species and for different tissues.

    Finally, in some species methylated cytosines have mutation rates an order of magnitude higher than other nucleotides. The authors mention they assume independence, but how would violation of this assumption affect their inference?

    Answer: Indeed, we assume the mutation and epimutation process to be independent thus the probability for a SNP to occur does not depend on the local methylation state. If this was the case, the mutation rate use would indeed be wrong to a degree function of the dependency between the processes. We suggest that by ignoring this dependence, we are in the same situation as ignoring the variation of mutation rate along the genome. We have previously documented the effect of ignoring this biological feature of genomes in Strüt et al 2023 and Sellinger et al 2021. The variation in mutation rate along the genome if too extreme and not accounted for can lead to erroneous inference results. However, this problem could be easily solved (modelled) by adapting the emission matrix. To correctly model this dependency, additional knowledge is needed: either the mutation and epimutation rates must be known to quantify the dependency, or the dependency must be known to quantify the resulting rates. As far as we know, these data are at the moment not available, but could maybe be obtained using the MA lines of A. thaliana (used in Yao et al. 2023).

    Recommendations for the authors:

    All three reviewers liked this approach and found it a valuable contribution. I think it is important to address reviewer 1/3 concerns about quantifying the accuracy of inference (the TMRCA approach from reviewer 1 sounds pretty reasonable), and reviewer 1 also highlights an intriguing point about model accuracy being worse when the mutation rate is known. Additionally, I think some discussion is warranted about challenges dealing with empirical methylation data (points from Rev 2 and 3 as well as Rev 1's question about inferred vs published rates of epigenetic mutation).

    Answer : We have added tables containing the root mean square error (RMSE) of every demographic inference in the manuscript to better quantify accuracy. We have below given the explanation on why accuracy in presence of site and region epimutations can in some cases decrease when real rates are known (because methylation state at the region level needs to be first inferred). We added evidence that accounting for methylation can improve the accuracy when recovering the TMRCA along the genome when the rates are known. We also have enhanced the discussion on the challenges of dealing with epimutations data for inference. As is suggested, we hope this study will generate an interest in tackling these challenges by applying the methods to various methylome datasets from different species.

    Reviewer #1 (Recommendations For The Authors):

    Major comments:

    • For all of the simulated demographic inference results, only plots are presented. This allowsfor qualitative but not quantitative comparisons to be made across different methods. It is not easy to tell which result is actually better. For example, in Supp. Fig. 5, eSMC2 seems slightly better in the ancient past, and times the trough more effectively, while SMCm seems a bit better in the very recent past. For a more rigorous approach, it would be useful to have accompanying tables that measure e.g. mean-squared error (along with confidence intervals) for each of the different scenarios, similar to what is already done in Tables 1 and 2 for estimating $r$.

    Answer : We understand the concern of reviewer #1 for a more quantitative approach to compare the inference results. We agree that plots are not sufficient to fully grasp a method performance. To provide better supports to quantity approaches performance, we added Sup tables 1,6,8,9 and 10 containing the RMSE (in log10 for visibility) for all Figures. The root mean-squared error is calculated as in Sellinger 2021 and a description of how the root mean-squared error is calculated and now found in the method section lines 886-893.

    • 434: The discussion downplays the really odd result that inputting the true value of themutation rate, in some cases, produces much worse estimates than when they are learned from data (SFig. 6)! I can't think of any reason why this should happen other than some sort of mathematical error or software bug. I strongly encourage the authors to pin down the cause of this puzzling behaviour.

    Answer : There are unfortunately no errors in this plot and those results are perfectly normal and coherent, but we understand they can be confusing at first.

    As described in the method section and in the appendix, when accounting for regionlevel epimutations, our algorithm requires the regional methylation status which needs to be inferred as a first step from the data (real or simulated). Because region and single site epimutation events are occurring at similar rates in our simulated scenario, the methylation state of the region is very hard to correctly recover (e.g. there will be unmethylated site in methylated regions and methylated sites in unmethylated regions). In other words, the accuracy of the region estimation HMM procedure is decreased by the joint action of site and region epimutation processes.

    When subsequently applying the HMM for inference, as described in the appendix, the probabilities of two CG site being in the same or different methylation state depends on the methlylation state of the "region". Hence the mislabelling of the region methylation state is (to some extent) equivalent to spurious SMPs (or inaccurate SMP calling).

    If the true rates for site and region epimutations are given as input, the model forces the demography (and other inferred parameters) to fit the observed distribution of SMPs (given the inputted rates), resulting in the poor accuracy observed in the Figure (Now Supplementary Figure 7).

    Note: The estimated rates from real data in A. thaliana suffer from the same issue as the region and site epimutation rates are independently estimated, and the existence of regions first quantified using an independent HMM method (Denkena et al. 2022).

    However, when rates are freely inferred, they are inferred accordingly to the estimated methylation status of regions and SNPs. Therefore, even if the inferred rates are wrong, they are used by the SMC in a more consistent way.

    Note: When methylation rates violate the infinite site assumption, such as here, we first estimate the tree sequence along the genome using SNPs (i.e. DNA mutations). The algorithm then infers the epimutations rates given the inferred coalescent times and the observed methylation diversity.

    To summarise: when inputting rates to the model, if the model fails to correctly recover the region methylation status there will be conflicting information between SNPs and SMPs leading to accuracy loss. However if the rates are inferred this is realized with the help of SNPs, leading to less conflicting information and potentially smaller loss of accuracy. We apologize that the explanations were missing from the manuscript and have added them lines 449-460 and 702-716.

    A further argument is that if region and site epimutations occur at rates of at least two orders of magnitude difference, the inference results are better (and accurate) when the true rates are given. The reason is that one epimutational process overrides the other (see Supplementary Table 2). In that case one epimutation process is almost negligible and we fall back to results from Figure 5 or Supplementary Figure 6.

    • As noted at 580, all of the added power from integrating SMPs/DMRs should come fromimproved estimation of recent TMRCAs. So, another way to study how much improvement there is would be to look at the true vs. estimated/posterior TMRCAs. Although I agree that demographic inference is ultimately the most relevant task, comparing TMRCA inference would eliminate other sources of differences between the methods (different optimization schemes, algorithmic/numerical quirks, and so forth). This could be a useful addition, and may also give you more insight into why the augmented SMC methods do worse in some cases.

    Answer : We fully agree with reviewer 1. We have added a comparison in TMRCA inference as proof of principle between using or not using methylation sites. The results are written in Supplementary Table 7 and methodology is inspired by Schiffels 2014 and described at the end of the method section (line 894-907). Those results demonstrate the potential gain in accuracy when using methylation polymorphic. However, TMRCA (or ARG) inference is a very vast and complex subject in its own right. Therefore, we are developing a complete TMRCA/ARG inference investigation and an improve methodology than the one presented in this manuscript. To do so we are currently working on a manuscript focusing on this topic specifically. We hence consider further investigations of TMRCA/ARG inference beyond the scope of this current study.

    • A general remark on the derivations in Section 2 of the supplement: I checked theseformulas as best I could. But a cleaner, less tedious way of calculating these probabilities would be to express the mutation processes as continuous time Markov chains. Then all that is needed is to specify the rate matrices; computing the emission probabilities needed for the SMC methods reduces to manipulating the results of some matrix exponentials. In fact, because the processes are noninteracting, the rate matrix decomposes into a Kronecker sum of the individual rate matrices for each process, which is very easy to code up. And this structure can be exploited when computing the matrix exponential, if speed is an issue.

    Answer: We thank the reviewer for this very interesting suggestion! Unfortunately, it is a bit late to re-implement the algorithm and reshape the manuscript according to this suggestion. Speed is not yet an issue but will most likely become one in the future when integrating many different rates or when using a more complex SMC model. Hence, we added reviewer #1 suggestions to the discussion (line 648) and hope to be using it in our future projects.

    • Most (all?) of the SNP-only SMC methods allow for binning together consecutiveobservations to cut down on computation time. I did not see binning mentioned anywhere, did you consider it? If the method really processes every site, how long does it take to run?

    Answer: This is a very good question. We do the binning exactly as described in Mailund 2013 & Terhorst 2017, and added this information in the method section (lines 801-809). However, as described in Terhorst 2017, one can only bin observation of the same "type" (to compute the Baum-Welch algorithm). Therefore, the computation time gain by binning is reduced when different markers spread along the genome in high proportion. This is the approach we used throughout the study when facing multiple markers as it had the best speed performance. As for example, when the proportion of site with methylated information is 1% or less, computation time is only slightly affected (i.e. same order of magnitude).

    However, the binning method presented in Mailund 2013 can be extended to observation of different types, but parameters need to be estimated through a full likelihood approach (as presented in Figure 2). In our study this approach did not have the best speed performance. However, as our study is the first of its kind, it remains sub-optimal for now. Hence, we did not further investigate the performance of our approach in presence of many multiple different genomic marker (e.g. 5 different markers each representing ~20% of the genome each). Currently, with SMC approaches a high proportion of sites contain the information "No SNPs", making the Baum welch algorithm described in Terhorst 2017 very efficient. But when further developing our theoretical approach, we expect that most of the sites in a genome analysis will contain some "information", which could render the full likelihood approach computationally more tractable.

    • 486: The assumed site and region (de)methylation rates listed here are several OOMdifferent from what your method estimated (Supp. Tables 5-6). Yet, on simulated data your method is usually correct to within an order of magnitude (Supp. Table 4). How are we to interpret this much larger difference between the published estimates and yours? If the published estimates are not reliable, doesn't that call into question your interpretation of the blue line in Fig. 7 at 533?

    Answer: We thank the reviewer for asking this question. We believe answering this question is indeed the most interesting aspect of our study. Beyond demographic inference, our study has indeed unveiled a discrepancy between rates inferred through biological experiment and our study through the use of SNPs and branch length. There are several reasons which could explained the discrepancy between both approaches:

    • Firstly, our underlying HMM hypotheses are certainly violated. We ignoredpopulation structure, variation of mutations and recombination rate along the genome as well as the effect of selection. Hence, the branch lengths used for methylation rate estimations are to some extent inaccurate. We note that this is especially likely for the short branches of coalescent tree originating from background selection events in the coding regions and which are especially observable when using the methylation sites with a higher mutation rate than SNPs (Yao et al. 2023) at body methylated genes.

    • Secondly, calling single methylation site polymorphism is not 100 % reliable. If theerror rate is 0.1%, as the study was conducted on ~10 generations a minimum epimutation rate of 10-4 is to be expected. However, because our approach works at the evolutionary time scale, we expect that it suffers less from this bias as the proportion of diversity originating from actual epimutations, and not SMP calling error, should be greater.

    • Thirdly, as mentioned above, recovering the methylation status of a region is veryhard. Hence false region status inference could affect our inference accuracy as shown in Supplementary Figure 4.

    • Lastly and most importantly, the reason behind this discrepancy is the modelling ofepimutation and methylation between sites and regions. As we discuss, the current combination of rates and models is still limited to describe the observed diversity along the genome (as we intend in SMC methods). This is in contrast to the recent study by Yao et al. where very few regions of polymorphic SMPs are chosen, which implicitly avoids the influence of the methylation region effect. A study just published by Biffra et al. (Cell reports 2023) also uses a functional model of methylation modelling using a mix of region and site epimutation, albeit not tuned for evolutionary analyses. Thus we suggest, in line with functional studies, that epimutations are not independent from the local methylation context and may tend to stabilize the methylation state of a region. Therefore, the estimated methylation rates show a discrepancy to the previously measured ones. Indeed, the biological experiment would reveal a fast epimutation rate because epimutations can actually be tracked at sites which can mutate, while region mutation rate is much slower. However, because the methylation state of a region is rather stable through time it would reduce the methylation diversity over long time scale, and these rates would differ between methylated or unmethylated regions (i.e. the methylation rate is higher in methylated regions). Our results are thus in agreement with the observation by Biffra et al. that region methylation modelling is needed to explain patterns of methylation across the genome.

    To solve the discrepancy, one would need to develop a theoretical region + site epimutation model capable of describing the observed diversity at the evolutionary time scale (possibly based on the Biffra et al. model within an underlying population evolution model), and then use this model to reanalyse the sequence data from the biological experiment (i.e. in de Graaf et al. 2015 & Denkena et al. 2022) to re-estimate the methylation region sizes and epimutation rates.

    Minor comments:

    • 189: "SMCtheo" first occurs here, but it's not mentioned until 247 that this is the newmethod being presented.

    Answer : Fixed

    • 199: Are the estimates in this section from a single diploid sequence? Or is it n=5 (diploid) as mentioned in the earlier section?

    Answer : Yes, those results were obtained with 5 diploid individuals. We added it in the Table 1 description.

    • 336: I'm confused by the wording: it sounds like the test rejects the null if there is positivecorrelation in the methylation status across sites. But then, shouldn't 339 read "if the test is significant" (not non-significant)?

    Answer : We apologize for the confusion and rewrote the sentence line 339-348, the choice of word was indeed misleading .

    • Fig. 6: for some reason fewer simulations were run for 10Mb (panels C nad D) than for100Mb (A and B). Since it's very difficult to tell what's happening on average in the 10Mb case, I suggest running the same number of simulations.

    Answer : Yes we understand your concern. Actually, the same number of simulations were run but we plotted only the first 3 runs as it was less visually confusing. We now have added the missing lines to the plot C and D.

    Typos:

    • 104: "or or"
    • 292: build => built
    • 388: fulfil
    • 683: sample => samples

    Answer : Many thanks to reviewer 1 for pointing out the typos. They are all now fixed.

    Reviewer #2 (Recommendations For The Authors):

    The authors may find some valuable information in Pisupati et al (2023) "On the causes of gene-body methylation variation in Arabidopsis thaliana" on interpreting epimutation rates.

    Answer: Many thanks for the recommended manuscript. We add it to the cited literature as it strongly supports our use of heritability or methylation. We also added the recent Biffra et al. paper.

    Reviewer #3 (Recommendations For The Authors):

    There are many places throughout the manuscript with minor grammatical errors. Please review these. A few noted below as I read:

    L104: extra "or"

    L123: built not build

    L 160 "relies" instead of "do rely"

    L161 "events"

    L 336 "from methylation data"

    L 378 "exists"

    L 379 "regions are on average shorter" instead of "there are shorter"

    L 338 "a regional-level"

    L 349 "," instead of "but"

    L 394 DMRs

    Table 1 legend: parentheses not brackets?

    Answer : Many thanks to reviewer #3 for finding those mistakes. They are all now fixed.

    I think a paragraph in the discussion of considerations of when to use this approach might be helpful to readers. Comparison to e.g. increased sample size in MSMC2, while not necessary, might be helpful here. It may often be the case that doubling the number of haplotypes with SNP data may be easier and cheaper estimating methylation accurately.

    Answer : We discuss (lines 691-698) that our approach is always useful by design, but cannot always be used for the same purpose. If the evolutionary properties of the used marker used are not understood, we suggest that our approach can be used to investigate the marker heritability process through generations. This could help to correctly design experiments aiming to study the marker heritability through lineages. And if the properties of the marker are well understood and modelled, it can be integrated into the SMC framework to improve inference accuracy.

    Other minor notes:

    L 486 "known" is a stretch. empirically estimated seems appropriate.

    Answer : Fixed

    L 573 ARG? You are not estimating the full ARG here.

    Answer : We apologize for the wrong choice of word and have rephrased the sentence.

    Fig. 2 is not super useful and could be supplemental.

    Answer : We moved Figure 2 to the appendix (now sup fig 1)

  4. eLife

    eLife assessment

    This important study extends existing sequentially Markovian coalescent approaches to include the use of hypervariable loci such as epimutations. This is an intriguing addition 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 -- and thus the potential use of this approach -- there would be additional value in more extensive consideration of when and where these methods are best used.

  5. eLife

    Reviewer #1 (Public Review):

    The authors developed an extension to the pairwise sequentially Markov coalecent model that allows to simultaneously analyze multiple types of polymorphism data. In this paper, they focus on SNPs and DNA methylation data. Since methylation markers mutate at a much faster rate than SNPs, this potentially gives the method better power to infer size history in the recent past. Additionally, they explored a model where there are both local and regional epimutational processes.

    Integrating additional types of heritable markers into SMC is a nice idea which I like in principle. However, a major caveat to this approach seems to be a strong dependence on knowing the epimutation rate. In Fig. 6 it is seen that, when the epimutation rate is known, inferences do indeed look better; but this is not necessarily true when the rate is not known. (See also major comment #1 below about the interpretation of these plots.) A roughly similar pattern emerges in Supp. Figs. 4-7; in general, results when the rates have to be estimated don't seem that much better than when focusing on SNPs alone. This carries over to the real data analysis too: the interpretation in Fig. 7 appears to hinge on whether the rates are known or estimated, and the estimated rates differ by a large amount from earlier published ones.

    Overall, this is an interesting research direction, and I think the method may hold more promise as we get more and better epigenetic data, and in particular better knowledge of the epigenetic mutational process. At the same time, I would be careful about placing too much emphasis on new findings that emerge solely by switching to SNP+SMP analysis.

    Major comments:
    - For all of the simulated demographic inference results, only plots are presented. This allows for qualitative but not quantitative comparisons to be made across different methods. It is not easy to tell which result is actually better. For example, in Supp. Fig. 5, eSMC2 seems slightly better in the ancient past, and times the trough more effectively, while SMCm seems a bit better in the very recent past. For a more rigorous approach, it would be useful to have accompanying tables that measure e.g. mean-squared error (along with confidence intervals) for each of the different scenarios, similar to what is already done in Tables 1 and 2 for estimating $r$.

    - 434: The discussion downplays the really odd result that inputting the true value of the mutation rate, in some cases, produces much worse estimates than when they are learned from data (SFig. 6)! I can't think of any reason why this should happen other than some sort of mathematical error or software bug. I strongly encourage the authors to pin down the cause of this puzzling behaviour. (Comment addressed in revision. Still, I find the explanation added at 449ff to be somewhat puzzling -- shouldn't the results of the regional HMM scan only improve if the true mutation rate is given?)

    - As noted at 580, all of the added power from integrating SMPs/DMRs should come from improved estimation of recent TMRCAs. So, another way to study how much improvement there is would be to look at the true vs. estimated/posterior TMRCAs. Although I agree that demographic inference is ultimately the most relevant task, comparing TMRCA inference would eliminate other sources of differences between the methods (different optimization schemes, algorithmic/numerical quirks, and so forth). This could be a useful addition, and may also give you more insight into why the augmented SMC methods do worse in some cases. (Comment addressed in revision via Supp. Table 7.).

    - A general remark on the derivations in Section 2 of the supplement: I checked these formulas as best I could. But a cleaner, less tedious way of calculating these probabilities would be to express the mutation processes as continuous time Markov chains. Then all that is needed is to specify the rate matrices; computing the emission probabilities needed for the SMC methods reduces to manipulating the results of some matrix exponentials. In fact, because the processes are noninteracting, the rate matrix decomposes into a Kronecker sum of the individual rate matrices for each process, which is very easy to code up. And this structure can be exploited when computing the matrix exponential, if speed is an issue.

    - Most (all?) of the SNP-only SMC methods allow for binning together consecutive observations to cut down on computation time. I did not see binning mentioned anywhere, did you consider it? If the method really processes every site, how long does it take to run?

    - 486: The assumed site and region (de)methylation rates listed here are several OOM different from what your method estimated (Supp. Tables 5-6). Yet, on simulated data your method is usually correct to within an order of magnitude (Supp. Table 4). How are we to interpret this much larger difference between the published estimates and yours? If the published estimates are not reliable, doesn't that call into question your interpretation of the blue line in Fig. 7 at 533? (Comment addressed in revision.)

  6. eLife

    Reviewer #2 (Public Review):

    A limitation in using SNPs to understand recent histories of genomes is their low mutation frequency. Tellier et al. explore the possibility of adding hypermutable markers to SNP based methods for better resolution over short time frames. In particular, they hypothesize that epimutations (CG methylation and demethylation) could provide a useful marker for this purpose. Individual CGs in Arabidopsis tends to be either close to 100% methylated or close to 0%, and are inherited stably enough across generations that they can be treated as genetic markers. Small regions containing multiple CGs can also be treated as genetic markers based on their cumulative methylation level. In this manuscript, Tellier et al develop computational methods to use CG methylation as a hypermutable genetic marker and test them on theoretical and real data sets. They do this both for individual CGs and small regions. My review is limited to the simple question of whether using CG methylation for this purpose makes sense at a conceptual level, not at the level of evaluating specific details of the methods. I have a small concern in that it is not clear that CG methylation measurements are nearly as binary in other plants and other eukaryotes as they are in Arabidopsis. However, I see no reason why the concept of this work is not conceptually sound. Especially in the future as new sequencing technologies provide both base calling and methylating calling capabilities, using CG methylation in addition to SNPs could become a useful and feasible tool for population genetics in situations where SNPs are insufficient.

  7. eLife

    Reviewer #3 (Public Review):

    I very much like this approach and the idea of incorporating hypervariable markers. The method is intriguing, and the ability to e.g. estimate recombination rates, the size of DMRs, etc. is a really nice plus. I am not able to comment on the details of the statistical inference, but from what I can evaluate it seems reasonable and in principle the inclusion of highly mutable sties is a nice advance. This is an exciting new avenue for thinking about inference from genomic data. I remain a bit concerned about how well this will work in systems where much less is understood about methylation,

    The authors include some good caveats about applying this approach to other systems, but I think it would be helpful to empiricists outside of thaliana or perhaps mammalian systems to be given some indication of what to watch out for. In maize, for example, there is a non-bimodal distribution of CG methlyation (35% of sites are greater than 10% and less than 90%) but this may well be due to mapping issues. The authors solve many of the issues I had concerns with by using gene body methylation, but this is only briefly mentioned on line 659. I'm assuming the authors' hope is that this method will be widely used, and I think it worth providing some guidance to workers who might do so but who are not as familiar with these kind of data.

  8. Version published to 10.1101/2023.04.02.535252v4 on bioRxiv
  9. Version published to 10.1101/2023.04.02.535252v3 on bioRxiv
  10. Version published to 10.7554/elife.89470.1 on eLife
  11. eLife

    eLife assessment

    This important study advances existing approaches for demographic inference by incorporating rapidly mutating markers such as switches in methylation state. The authors provide a solid comparison of their approach to existing methods, although the work would benefit from some additional consideration of the challenges in the empirical use of methylation data. The work will be of broad interest to population geneticists, both in terms of the novel approach and the statistical inference proposed.

  12. eLife

    Reviewer #1 (Public Review):

    The authors developed an extension to the pairwise sequentially Markov coalecent model that allows to simultaneously analyze multiple types of polymorphism data. In this paper, they focus on SNPs and DNA methylation data. Since methylation markers mutate at a much faster rate than SNPs, this potentially gives the method better power to infer size history in the recent past. Additionally, they explored a model where there are both local and regional epimutational processes.

    Integrating additional types of heritable markers into SMC is a nice idea which I like in principle. However, a major caveat to this approach seems to be a strong dependence on knowing the epimutation rate. In Fig. 6 it is seen that, when the epimutation rate is known, inferences do indeed look better; but this is not necessarily true when the rate is not known. A roughly similar pattern emerges in Supp. Figs. 4-7; in general, results when the rates have to be estimated don't seem that much better than when focusing on SNPs alone. This carries over to the real data analysis too: the interpretation in Fig. 7 appears to hinge on whether the rates are known or estimated, and the estimated rates differ by a large amount from earlier published ones.

    Overall, this is an interesting research direction, and I think the method may hold more promise as we get more and better epigenetic data, and in particular better knowledge of the epigenetic mutational process. At the same time, I would be careful about placing too much emphasis on new findings that emerge solely by switching to SNP+SMP analysis.

  13. eLife

    Reviewer #2 (Public Review):

    A limitation in using SNPs to understand recent histories of genomes is their low mutation frequency. Tellier et al. explore the possibility of adding hypermutable markers to SNP based methods for better resolution over short time frames. In particular, they hypothesize that epimutations (CG methylation and demethylation) could provide a useful marker for this purpose. Individual CGs in Arabidopsis tends to be either close to 100% methylated or close to 0%, and are inherited stably enough across generations that they can be treated as genetic markers. Small regions containing multiple CGs can also be treated as genetic markers based on their cumulative methylation level. In this manuscript, Tellier et al develop computational methods to use CG methylation as a hypermutable genetic marker and test them on theoretical and real data sets. They do this both for individual CGs and small regions. My review is limited to the simple question of whether using CG methylation for this purpose makes sense at a conceptual level, not at the level of evaluating specific details of the methods. I have a small concern in that it is not clear that CG methylation measurements are nearly as binary in other plants and other eukaryotes as they are in Arabidopsis. However, I see no reason why the concept of this work is not conceptually sound. Especially in the future as new sequencing technologies provide both base calling and methylating calling capabilities, using CG methylation in addition to SNPs could become a useful and feasible tool for population genetics in situations where SNPs are insufficient.

  14. eLife

    Reviewer #3 (Public Review):

    I very much like this approach and the idea of incorporating hypervariable markers. The method is intriguing, and the ability to e.g. estimate recombination rates, the size of DMRs, etc. is a really nice plus. I am not able to comment on the details of the statistical inference, but from what I can evaluate it seems sound and reasonable. This is an exciting new avenue for thinking about inference from genomic data. I have a few concerns about the presentation and then also questions about the use of empirical methylation data sets.

    I think a more detailed description of demographic accuracy is warranted. For example, in L245 MSMC2 identifies the bottleneck (albeit smoothed) and only slightly overestimates recent size. In the same analysis the authors' approach with unknown mu infers a nonexistent population increase by an order of magnitude that is not mentioned.

    Similarly, it seems problematic that (L556) the approach requiring estimation of site and region parameters (as would presumably be needed in most empirical systems like endangered nonmodel species mentioned in the introduction) does no better than using only SNPs. Overall, I think a more objective and perhaps quantitative comparison of approaches is warranted.

    The authors simulate methylated markers at 2% (and in some places up to 20%). In many plant genomes a large proportion of cytosines are methylated (e.g. 70% in maize: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC8496265/). I don't know what % of these may be polymorphic, but this leads to an order of magnitude more methylated cytosines than there are SNPs. Couldn't this mean that any appreciable error in estimating methylation threatens to be of a similar order of magnitude to the SNP data? I would welcome the authors' thoughts here.

    A few points of discussion about the biology of methylation might be worth including. For example, methylation can differ among cell types or cells within a tissue, yet sequencing approaches evaluate a pool of cells. This results in a reasonable fraction of sites having methylation rates not clearly 0 or 1. How does this variation affect the method? Similarly, while the authors cite literature about the stable inheritance of methylation, a sentence or so more about the time scale over which this occurs would be helpful. Finally, in some species methylated cytosines have mutation rates an order of magnitude higher than other nucleotides. The authors mention they assume independence, but how would violation of this assumption affect their inference?

  15. Version published to 10.1101/2023.04.02.535252v2 on bioRxiv
  16. Version published to 10.1101/2023.04.02.535252v1 on bioRxiv