Theoretical modelling of epigenetically modified DNA sequences

Gas-phase base pairs

To examine the effect of modifications on base pairing we examined the structure and energy of gas-phase CG pairs in both hydrogen bonded and stacked orientations, with results reported in Table 1 and Table 2 respectively. These data show that methylation has little effect on the geometry or stability of the Watson-Crick base pair. The presence of a hydroxymethyl slightly weakens the N₄-H₄…O₆ H-bond, perhaps due to the proximity of CH₂OH and NH₂ groups, reported as X…H₄ in Table 1. Formyl has a larger effect overall, lengthening N₃…H₁-N₁ and O₂…H₂-N₂ H-bonds and hence reducing binding by over 3 kcal/mol. The pattern of changes induced by carboxylate is different from all other modifications, lengthening the peripheral H-bonds N₄-H₄…O₆ and O₂…H₂-N₂ markedly, but shortening N₃…H₁-N₁. Despite this weakening, the carboxylate-substituted cytosine binds most strongly to guanine, presumably due to ion-dipole interactions within the anionic system. Both formyl and carboxylate contain close O…H₄ contacts, but overall the proximity of these groups does not appear to be related to strength or geometry of binding.

As well as the effect on H-bonding, epigenetic modifications can alter the stacking behaviour of DNA bases. Table 2 reports geometrical details, as well as binding energies, of the five modified cytosines considered here stacked with guanine. All such calculations started from the idealised B-DNA orientation (Cent…Cent = 4.390 Å, Dihedral = 4.9°), and overall this is retained in our gas-phase DFT optimisation. Table 2 shows that methylation leads to closer contact and greater stabilisation between bases, as might be expected due to the increased polarizability of this modified base. Hydroxymethylation leads to the most stable pair considered here, largely due to a strong H-bond between the H—O of hydroxymethyl and O6 of guanine (H…O = 1.770 Å), whereas formylation leads to longer, weaker interaction between bases. Carboxylate-substituted cytosine is the only case considered here that loses the approximately parallel orientation of bases. This appears to be driven as much by repulsion between the carboxylate group and C=O₆ of guanine as by H-bonding.

Double-stranded DNA trimers

While these gas-phase dimers give useful information on the intrinsic effect of modifications on cytosine’s ability to interact with guanine, environmental effects including the DNA sequence, sugar-phosphate backbone and solvent will play a major role in determining their effect in real systems. In order to better simulate the behaviour of modified cytosines in real systems, structures of double-stranded d(GCG) and d(ACA), as well as epigenetic modifications to the central cytosine were optimized using DFT in continuum solvent (PCM), and the resulting geometries of the local base pairs were analysed in the coordinate frame recommended by Olson et al.³⁶. Unlike the free dimers considered above, modifications have only subtle effects on this larger structure, which retains the overall canonical B-DNA shape of the unmodified WT structure.

Following Zubatiuk et al.³⁷, we summarise key aspects of trimer structure, which are displayed graphically in Figure 2 and Figure 3. The corresponding values are tabulated in Table S1 of the Supporting Information, with the base step and local helical parameters tabulated in Table S2. As with Zubatiuk et al.¹⁴, base pair step parameters are averaged over 3´ and 5´ directions. In the GCG oligomer, methylation has only a small effect on base pair distances, but does alter the propeller angle by over 4°. Hydroxymethylation has a larger effect on the GCG oligomer, especially on the stagger, buckle and propeller, whereas the stretch and opening parameters are much less affected. Formyl does not strongly affect base pair distances but does change angles substantially, especially buckle and propeller, which change by as much as 10°. In contrast, carboxylate induces a large change in stagger but only small changes in angular geometry. Base pair step parameters for d(GCG) in general are less affected than those for the base pair noted above, with the exception of formyl which exhibits smaller slide and less negative roll values than unmodified DNA.

Figure 2. Base pair (A) and step (B, C) parameters for central GC in d(GCG) and modifications (Å and °).

The corresponding data are provided in Table S1.

Figure 3. Base pair (A) and step (B, C) parameters for central GC in d(ACA) and modifications (Å and °).

The corresponding data are provided in Table S1.

Rather larger changes are evident on modification of d(ACA), as shown in Figure 3. In this case, even methylation induces significant changes in distances, especially stagger which increases by 0.1 Å, and angles (buckle and propeller change by 8 and 13°, respectively). At the base pair step level, methylation gives rise to substantial increase (0.9 and 1.5 Å) in slide and more negative roll in both 3´ and 5´ directions. Less apparent in Figure 3, but still notable, are changes in rise that are 0.1 and 0.3 Å smaller in the methylated structure, reflecting the greater stacking that results from addition of a methyl group. Other modifications induce different patterns of structural change: for the central base pair these changes are typically smaller than for methylation, but for base pair steps much larger changes are found in some parameters. Most notable of these are slide, which changes by over 3 Å and roll (up to 17°) in the 3´ direction, in a similar way to that reported previously for smaller systems^17,18. Other parameters such as the width of the DNA strand, measured as the distance between C1´ nuclei, and virtual angles λ_Y and λ_R, which describe the pivoting of complementary bases in the base-pair plane, vary only slightly from the idealised values for B-DNA.

QM/MM studies of double-stranded oligomers

The oligomers considered so far are close to the limit of our computational capabilities of current DFT methods (the largest structure, carboxylated d(ACA), has 962 electrons in 2743 basis functions), such that longer sequences cannot currently be routinely studied in this manner. However, they are too small to correctly represent how DNA behaves in a real system, where the conformations adopted by each base pair step depend on the neighbouring step. Moreover, simulations of nucleic acids are known to suffer problems due to greater elasticity of the terminal part of the structure (the so called “end-effect”³⁸). For these reasons, these small oligomers are inadequate models to probe the effects of epigenetic modifications on the structure of DNA. We therefore turn to hybrid QM/MM methods, in which a subset of the atoms in the system is treated with DFT, and the remainder of the system with much faster molecular mechanics methods. In order to test the validity of this approach, methylated GCG was optimized using either only two or six bases in the QM region (Figure 4). These tests show that including only two bases in the QM region leads to significant differences in geometry to that obtained from DFT, particularly in the stagger and buckle coordinates. In contrast, including six bases in the QM region reproduces the DFT structure reasonably well. Similar observations were made from analogous treatment of methylated ACA (data not shown).

Figure 4. Comparison of QM/MM with DFT geometry (Å and °).

The corresponding values are shown in Tables S3 and S4.

As a further test, we also compared DFT and QM/MM derived structures with those optimised using the force field parameters developed in our previous work. Figure 5 shows the base-pair parameter values of the methylated structure d(GC´G) for the different methods. The MM structures provide very close values to those obtained by both QM/MM and DFT approaches, showing slight difference only in the stagger and propeller angle. We can therefore conclude that for small DNA oligomers, DFT, QM/MM and MM methods can all produce almost equally adequate DNA structures, but that QM/MM and MM approaches are more similar to one another than those obtained from DFT alone.

Figure 5. Comparison of DFT, QM/MM and molecular mechanics energy minimised (EM) geometries of d(GC´G) (Å and °).

The corresponding values are shown in Tables S3 to S6.

QM/MM geometry optimization with six bases in the QM region was then applied to a set of larger DNA sequences. The experimental structure of Renčiuk et al.⁴ obtained using X-ray diffraction (PDB Entry 4GLG) was truncated to a sequence of 7 base pairs, i.e. 5´-ATT CGCG-3´, and the central 6 bases (TCG//CGA) assigned as QM atoms. The remaining atoms, including crystallographic water molecules and counterions, were assigned to the MM layer, and the entire system was geometry optimized. The resulting optimized structure of the system with methylated C in the central position is shown in Figure6. Base pair and base pair step geometries of wild type, methylated, hydroxymethylated structures optimised with QM/MM, along with experimental values for methylated C, are shown in Figure 7.

Figure 6. QM/MM optimised structure of 5´-ATTCGCG-3´ with 5-mC in central position, and the bases defined as QM atoms shown in CPK.

A purple sphere highlights the methylation position, and water molecules and counterions have been omitted for clarity.

Figure 7. Base pair step parameters for central GC in 5´-ATTCGCG-3´ and modifications (Å and °).

The corresponding values are shown in Table S2.

We find that the structural effect of methylation is larger in this longer sequence than in the trimers considered above. Particularly, the optimized values of shear, stagger and buckle of the central base pair differ markedly between the methylated and WT forms of DNA. In contrast, the base pair step parameters exhibit rather smaller changes. For the hydroxymethylated structures, we observe similar profiles to the methylated structures. Furthermore, our simulations also allow us to probe the preferred orientation of the hydroxymethyl group: our DFT calculations predict a slight preference for the OH group to point in 3´ over 5´ and an optimized of C6-C5-C5A-O5 torsion angle of 118.4°, while previous MD simulations show this torsion to vary between 85 and 120° over 100 ns of simulation¹⁹. This is in good agreement with the experimental and theoretical results of Renčiuk et al.³⁹, who reported values between 72 and 133° using X-ray diffraction methods.