Author response:
The following is the authors’ response to the original reviews.
Reviewer #1 (Public Review):
Summary:
In this manuscript, the authors used a coarse-grained DNA model (cgNA+) to explore how DNA sequences and CpG methylation/hydroxymethylation influence nucleosome wrapping energy and the probability density of optimal nucleosomal configuration. Their findings indicate that both methylated and hydroxymethylated cytosines lead to increased nucleosome wrapping energy. Additionally, the study demonstrates that methylation of CpG islands increases the probability of nucleosome formation.
Strengths:
The major strength of this method is that the model explicitly includes elastic constraints on the positions of phosphate groups facing a histone octamer, as DNA-histone binding site constraints. The authors claim that their model enhances the accuracy and computational efficiency and allows comprehensive calculations of DNA mechanical properties and deformation energies.
Weaknesses:
A significant limitation of this study is that the parameter sets for the methylated and hydroxymethylated CpG steps in the cgNA+ model are derived from all-atom molecular dynamics (MD) simulations that suggest that both methylated and hydroxymethylated cytosines increase DNA stiffness and nucleosome wrapping energy (P´erez A, et al. Biophys J. 2012; Battistini, et al. PLOS Comput Biol. 2021). It could predispose the coarse-grained model to replicate these findings. Notably, conflicting results from other all-atom MD simulations, such as those by Ngo T in Nat. Commun. 2016, shows that hydroxymethylated cytosines increase DNA flexibility, contrary to methylated cytosines. If the cgNA+ model was trained on these later parameters or other all-atom force fields, different conclusions might be obtained regarding the effects of methylated and hydroxymethylation on nucleosome formation.
Despite the training parameters of the cgNA+ model, the results presented in the manuscript indicate that methylated cytosines increase both DNA stiffness and nucleosome wrapping energy. However, when comparing nucleosome occupancy scores with predicted nucleosome wrapping energies and optimal configurations, the authors find that methylated CGIs exhibit higher nucleosome occupancies than unmethylated ones, which seems to contradict their findings from the same paper which showed that increased stiffness should reduce nucleosome formation affinity. In the manuscript, the authors also admit that these conclusions “apparently runs counter to the (perhaps naive) intuition that high nucleosome forming affinity should arise for fragments with low wrapping energy”. Previous all-atom MD simulations (P´erez A, et al. Biophys J. 2012; Battistini, et al. PLOS Comput Biol. 202; Ngo T, et al. Nat. Commun. 20161) show that the stiffer DNA upon CpG methylation reduces the affinity of DNA to assemble into nucleosomes or destabilizes nucleosomes. Given these findings, the authors need to address and reconcile these seemingly contradictory results, as the influence of epigenetic modifications on DNA mechanical properties and nucleosome formation are critical aspects of their study. Understanding the influence of sequence-dependent and epigenetic modifications of DNA on mechanical properties and nucleosome formation is crucial for comprehending various cellular processes. The authors’ study, focusing on these aspects, will definitely garner interest from the DNA methylation research community.
Training the cgNA+ model on alternative MD simulation datasets is certainly of interest to us. However, due to the significant computational cost, this remains a goal for future work. The relationship between nucleosome occupancy scores and nucleosome wrapping energy is still debated, with conflicting findings reported in the literature, as noted in our Discussion section. Interestingly, we find that our predicted log probability density of DNA spontaneously acquiring a nucleosomal configuration is a better indicator of nucleosome occupancy than our predicted DNA nucleosome wrapping energy.
Reviewer #2 (Public Review):
Summary:
This study uses a coarse-grained model for double-stranded DNA, cgNA+, to assess nucleosome sequence affinity. cgNA+ coarse-grains DNA on the level of bases and accounts also explicitly for the positions of the backbone phosphates. It has been proven to reproduce all-atom MD data very accurately. It is also ideally suited to be incorporated into a nucleosome model because it is known that DNA is bound to the protein core of the nucleosome via the phosphates.
It is still unclear whether this harmonic model parametrized for unbound DNA is accurate in describing DNA inside the nucleosome. Previous models by other authors, using more coarse-grained models of DNA, have been rather successful in predicting base pair sequence-dependent nucleosome behavior. This is at least the case as far as DNA shape is concerned whereas assessing the role of DNA bendability (something this paper focuses on) has been consistently challenging in all nucleosome models, to my knowledge.
It is thus of major interest whether this more sophisticated model is also more successful in handling this issue. As far as I can tell the work is technically sound and properly accounts for not only the energy required in wrapping DNA but also entropic effects, namely the change in entropy that DNA experiences when going from the free state to the bound state. The authors make an approximation here which seems to me to be a reasonable first step.
Of interest is also that the authors have the parameters at hand to study the effect of methylation of CpG-steps. This is especially interesting as it allows us to study a scenario where changes in the physical properties of base pair steps via methylation might influence nucleosome positioning and stability in a cell-type-specific way.
Overall, this is an important contribution to the question of how the sequence affects nucleosome positioning and affinity. The findings suggest that cgNA+ has something new to offer. But the problem is complex, also on the experimental side, so many questions remain open.
Strengths:
The authors use their state-of-the-art coarse-grained DNA model which seems ideally suited to be applied to nucleosomes as it accounts explicitly for the backbone phosphates.
Weaknesses:
(1) According to the abstract the authors consider two “scalar measures of the sequence-dependent propensity of DNA to wrap into nucleosomes”. One is the bending energy and the other, is the free energy. Specifically in the latter, the authors take the difference between the free energies of the wrapped and the free DNA. Whereas the entropy of the latter can be calculated exactly, they assume that the bound DNA always has the same entropy (independent of sequence) in its more confined state. The problem is the way in which this is written (e.g. below Eq. 6) which is hard to understand. The authors should mention that the negative of Eq. 6 is what physicists call free energy, namely especially the free energy difference between bound and free DNA.
We have included the necessary clarifications in the revised manuscript, below Eq. 6.
(2) In Eq. 5 the authors introduce penalty coefficients ci. They write that values are “set by numerical experiment to keep distances ... within the ranges observed in the PDB structure, while avoiding sterical clashes in DNA.” This is rather vague, especially since it is unclear to me what type of sterical clashes might occur. Figure 1 shows then a comparison between crystal structures and simulated structures. They are reasonably similar but standard deviations in the fluctuations of the simulation are smaller than in the experiments. Why did the authors not choose smaller ci-values to have a better fit? Do smaller values lead to unwanted large fluctuations that would lead to steric clashes between the two DNA turns? I also wonder what side views of the nucleosomes look like (experiments and simulations) and whether in this side view larger fluctuations of the phosphates can be observed in the simulation that would eventually lead to turn-turn clashes for smaller ci-values.
The side view plots of the experimental and predicted nucleosome structures are now added to Supplementary material (Figure S8). Indeed, smaller ci values lead to steric clashes between the two turns of DNA – this is now specified in the Methods section. A possible improvement of our optimisation method and a direction of future work would be adding a penalty which prevents steric clashes to the objective function. Then the ci values could be reduced to have bigger fluctuations that are even closer to the experimental structures. We added this explanation to the Results section.
Reviewer #3 (Public Review):
Summary:
In this study, the authors utilize biophysical modeling to investigate differences in free energies and nucleosomal configuration probability density of CpG islands and nonmethylated regions in the genome. Toward this goal, they develop and apply the cgNA+ coarse-grained model, an extension of their prior molecular modeling framework.
Strengths:
The study utilizes biophysical modeling to gain mechanistic insight into nucleosomal occupancy differences in CpG and nonmethylated regions in the genome.
Weaknesses:
Although the overall study is interesting, the manuscripts need more clarity in places. Moreover, the rationale and conclusion for some of the analyses are not well described.
We edited the manuscript according to the reviewer’s suggestions and hopefully improved its readability.
Reviewer #1 (Recommendations For The Authors):
(1) The cgNA+ model parameters are derived from all-atom molecular dynamics (MD) simulations, yet there is no consensus within all-atom MD simulations regarding the impact of CpG methylation on DNA mechanical properties. The authors could consider fitting the coarsegrained model with a different all-atom force field to verify whether the conclusions regarding the effects of methylation and hydroxymethylation on DNA nucleosome wrapping energies still hold. For further details on MD simulations related to CpG methylation effects, the authors are advised to consult the review paper by Li et al. (2022) titled “DNA methylation: Precise modulation of chromatin structure and dynamics” published in Current Opinion in Structural Biology.
Parametrizing the cgNA+ model using MD simulations with various force fields is certainly of interest to us. However, due to the computational cost involved, it remains a goal for future work.
(2) Beyond DNA mechanical properties, which are directly linked to nucleosome wrapping energies in this study, the authors might also consider other factors such as geometric properties that could influence nucleosome formation. This approach might help the authors to reconcile the observed higher nucleosome occupancy scores for methylated CpGs. The authors are encouraged to review the aforementioned paper for additional experimental and MD simulation studies that could support this perspective.
Geometric properties of DNA are directly incorporated into our method through the cgNA+ model equilibrium shape prediction µ. We compute the mechanical energy needed deform µ to a nucleosomal configuration. Notably, the equilibrium shape µ is sensitive to methylation, as demonstrated in Figure 3.
(3) There are some issues with citation accuracy in the manuscript. For instance, in the Discussion section, the authors attribute a statement to Collings et al. and Anderson (2017), claiming that “methylated regions, known to have high wrapping energy, are among the highest nucleosome occupied elements in the genome.” However, upon reviewing this paper, it appears that it does not make any claims about the high wrapping energy of methylated regions.
The paragraph is now edited and a separate citation, P´erez et al. (2012), is given for the statement that methylation regions have high wrapping energy.
Reviewer #2 (Recommendations For The Authors):
Please improve the readability by:
(1) making clear that -ln ρ in Eq. 6 on page 4 is actually the free energy. Also, the word entropy comes too late (on page 7) where the best explanation of Eq. 6 is presented.
We added a comment about -ln ρ being the free energy after Eq. 6 and also included an equation, relating ln ρ and entropy.
(2) page 12 and 13 show two sets of experimental data. They are quite different from each other. When reading this, I wondered why there is this difference. But only on page 16, you explain that these are different cell types. The difference should be explained already when the papers are introduced on page 12.
A corresponding sentence already appeared in page 12: “The observations about nucleosome occupancy should be regarded as preliminary, and be treated with caution, as they are based on experimental data obtained for the cancerous HeLa cells Schwartz et al. (2019) and human genome embryonic stem cells Yazdi et al. (2015)”. Now we also added this information to the first paragraph of the subsection for clarity.
Finally, I add here some general thoughts that came up when reading the paper, comparing your findings with earlier findings in the field. This is not a strict one-to-one comparison and thus does not have to find its way into this manuscript but might give ideas for future studies. Experiments suggest that nucleosomes prefer DNA with a high content of C’s and G’s. Figure 2 does not look at the GC content but at the number of CpG’s. But in any case, let’s use this as a proxy for GC content. Figure 2a suggests that there is not a strong dependence of the bending energy on the number CpG steps. This is consistent with earlier work with the rigid basepair model which shows the same behavior for GC content (for both MD and crystal parametrizations). Figure 2c (related to the negative free energy) shows that with an increasing number of CpG steps the propensity to bind goes down. This suggests that the entropic cost to confine CpG-rich DNA increases, which in turn reflects that these DNA stretches are softer. This is rather interesting since in the case of the rigid basepair model this effect is observed only when stiffnesses are extracted from crystal data not MD data (however, this refers again to CG content). This might indicate a difference between the rigid bp model and cgNA+ which will be interesting to study in the future. Interesting is also the effect of CpG methylation. The stiffer methylated steps lead to an increase in the energy with the number of such steps (Figure 2a). The entropic cost for binding is thus expected to be smaller and this is indeed observed in Figure 2c when compared to the non-methylated steps.
We thank the reviewer for this comment. As for the GC content, the energy and lnp plots are indeed very similar to those in Figure 2.
Reviewer #3 (Recommendations For The Authors):
(1) The formulation of the cgNA+ model in the method section was not easy to follow and can be described better to improve clarity.
We have revised the model description and hope that its clarity has been improved
(2) The authors mention utilizing 100 human genome sequences with 100 configurations from DB. It would be helpful to clarify the source of these 100 human genome sequences. Are these 100 distinct regions on the human reference genome, or are they from a specific dataset or database?
We now include an explanation about the origin of sequences: “The human genome sequences are a random subset of our sequence sample for the CGI and NMI intersection in the Chromosome 1, but the following observations remain unchanged for sequence samples from different genomic regions.”
(3) The authors mention the lack of tail unwrapping in their model. It would be beneficial to understand the magnitude of this issue and its potential impact on the overall results. How significant is the lack of unwrapping events in their current model?
We observed the unwrapping of approximately five base-pairs at each end of our predicted nucleosome configurations, in comparison to the experimental configurations (Figure 1). This issue could be solved by adding additional constraints at the ends of the 147 bp sequence. The wrapping energy would increase marginally, as only about 10 of 147 bp would be affected. We added this remark to the main text.
(4) Observations from Figure 3 are not described properly. Are these differences statistically significant? Why is twist higher for CpG sites but lower for a roll?
We added an explanation of how the statistics was computed into the caption of Figure 3. In fact, we didn’t use statistical estimates here, but generated all the possible cases and computed the exact statistics (for the given set of our model parameters). Regarding the changes in twist and roll, we have added the following comment on page 7: “The ground state changes resulting from cytosine modifications – primarily characterized by an average increase in roll and a decrease in twist – may be linked to steric hindrance caused by the cytosine 5-substituent (Battistini et al. (2021)). Notably, the negative coupling between twist and roll has already been observed in X-ray crystallography data (Olson et al. (1998)).”
(5) Figure 4 does not clarify the authors’ conclusion of higher stiffness for ApT and TpA dinucleotides. The authors should provide further explanation for this observation.
We revised the text to clarify that the statement regarding ApT and TpA being the most stiff and the most flexible dinucleotides is not a conclusion derived from Figure 4, but rather from earlier work that we cite.
(6) In Figure 7, the authors note that methylated CGIs have higher nucleosome occupancy on average than unmethylated sequences. Is this observation statistically significant?
We observe that methylated sequences have a higher average occupancy than unmethylated sequences in Yazdi et al. data, when the CpG count falls into the intervals from 5 to 14 and from 15 to 24. For each of the two intervals this difference is statistically significant: the permutation test, used due to the lack of normality, yields a p-value of 0.0001 for both cases. The differences in mean scores shown in Figure 8 are also statistically significant. Such test results are expected, given the large sample sizes and the observed differences in means, therefore we prefer not to include this discussion in main text.
(7) The authors note that their analyses to correlate nucleosome occupancy profile with the methylation state of underlying sequences are preliminary, as different cell lines were used to perform these analyses. Given this inconsistency, it needs to be clarified why this analysis was performed and what the takeaway is.
We added the following comment at the end of the Results section: “Although comparing data from different cell lines is not optimal, to the best of our knowledge, no publicly available methylation and nucleosome occupancy data exist for the entire human genome within the same cell type. Nevertheless, since the lowest log probability densities in the human genome are predicted for CpG-rich sequences regardless of their methylation state (Figure 2d), and the same holds for both sets of the nucleosome occupancy scores (Figure 7), we conclude that the lowest occupancies occur for sequences with the lowest log probability densities.”
