Hoogsteen breathing in a

Subscriber access provided by University of Winnipeg Library

Biomolecular Systems

Computational probing of Watson-Crick/Hoogsteen breathing in a DNA duplex containing N1-methylated adenine Changwon Yang, Eunae Kim, Manho Lim, and Youngshang Pak J. Chem. Theory Comput., Just Accepted Manuscript • DOI: 10.1021/acs.jctc.8b00936 • Publication Date (Web): 02 Dec 2018 Downloaded from http://pubs.acs.org on December 2, 2018

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 43 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

Computational probing of Watson-Crick/Hoogsteen breathing in a DNA duplex containing N1methylated adenine Changwon Yang1, Eunae Kim2, Manho Lim1, and Youngshang Pak1,*

1. Department of Chemistry and Institute of Functional Materials, Pusan National University, Busan, 46241, South Korea

2. College of Pharmacy, Chosun University, Gwangju, 61452, South Korea

KEYWORDS: Molecular dynamics, WC/HG breathing, N1-methylated adenine, Methylated DNA, Free energy surface

ACS Paragon Plus Environment

1

Journal of Chemical Theory and Computation 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 2 of 43

ABSTRACT

DNA breathing is a local conformational fluctuation spontaneously occurring in double-stranded DNAs. In particular, the possibility of individual base pairs (bps) in duplex DNA to flip between alternate bp modes, i.e. Watson-Crick (WC) like and Hoogsteen (HG) like, at relevant time scales has impacted DNA research fields for many years. In this study, to computationally probe effects of chemical modification on the DNA breathing, we present a free energy landscape of spontaneous thermal transitions between WC and HG bps in a free DNA duplex containing N1methylated adenine (m1A). For the current free energy computation, a variant of well-tempered metadynamics simulation was extensively performed for a total of 40 μs to produce free energy surfaces. The free energy profile indicated that upon the chemical modification of adenine, the HG bp (m1A•T) was located in the most favourable conformation (96.7%), however, the canonical WC bp (m1A•T) was distorted into two WC-like bps of WC* (2.8%) and WC** (0.5%). The conformational exchange between these two minor WC-like bps occurs with the first hundred nanoseconds. The transition between WC-like and HG bp features multiple


2

Page 3 of 43 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


transition pathways displaying various extents of base flipping in combination with glycosidic rotation. Analysis of the simulated ensemble showed that the m1A-induced changes of the backbone and sugar pucker were in a reasonable agreement with previous results inferred from NMR experiments. Also, this study revealed that the formation of the stable HG bp upon the mutation alters the characteristics of dynamic fluctuations of the neighboring WC residues of m1A. We expect this simulation approach to be a robust computational scheme to complement and guide future high-resolution experiments on many outstanding issues of duplex DNA breathing.

INTRODUCTION The sequence of base pairs (bps) in duplex DNAs encode the genetic information. Beyond these genetic codes, they can undergo spontaneous thermal fluctuations called “DNA breathing”.1-3 Upon such a DNA breathing, target base pairs in duplex DNAs can adopt local conformations deviating from the most favourable structure. One example of DNA breathing is a thermal transition from the Watson-Crick (WC) to Hoogsteen (HG) bp.3-7 In this WC/HG breathing fluctuation, a target bp can flip between alternate bp modes, i.e. WC-like and HG-like at time scales of 𝜇s – ms. Such bp transition can alter the global and local geometrical features of duplex DNAs. These structural


3


Page 4 of 43

perturbations of DNAs can guide specific enzymes to recognize and access DNA templates to carry out DNA processing reactions. In normal biological environments, the methylation of DNAs typically occurs, leading to blocking of DNA replication, repression of gene transcription, aging, and carcinogenesis. To repair damaged DNAs via methylation, repair enzymes such as AlkB and ABH2 should detect a methylation site of damaged DNAs.8-10 Thus, it is of fundamental interest to investigate how this damaged site is recognized by such enzymes at their early stage of repairing processes. Therefore, the investigation of detailed molecular events of WC/HG breathing fluctuation of a free duplex DNA can shed light on initial stages of recognition by repair enzymes. As a result, the WC/HG breathing on a free methylated duplex DNA as well as an unmodified one is both of fundamental and practical importance, ranging from further resolving dynamics processes of duplex DNAs to possible therapeutic applications resulting from these processes. Recently, an NMR relaxation dispersion study on a free duplex DNA (A6DNA; see Figure 1A) indicated that the most favourable WC AT bp exists in thermal equilibrium with a transient HG bp in aqueous solution.4-5 More recently, a follow-up


4

Page 5 of 43 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


NMR study was performed for methylated duplex DNAs by introducing a single methylation of the adenine base at the N1 position (m1A).11 Upon this N1-methylation of the adenine base, an inter-base steric clash destabilizes the m1AT WC bp, and the m1AT HG bp is then adopted as a stable conformer.9-12 Previously, a solution NMR structure of the m1AT HG bp was determined for a mutated DNA duplex containing the m1A lesion (m1A-A6DNA).11 From this NMR study, in combination with a molecular dynamics (MD) simulation, the m1A lesion was shown to induce various local structural perturbations, such as kinking of DNA, conformational changes of backbone, sugarpuckering, and possibly melting of neighboring WC bp. This experimental study on the normal and modified duplex DNA led to the characterization of both WC and HG bp states, providing molecular insight into the WC/HG breathing of DNA duplexes and its potential impact on enzyme recognitions. Despite such experimental advances in WC/HG breathing studies, in the case of the chemically modified DNA duplex (m1AA6DNA), detailed atomistic views on the WC/HG breathing have not been reported. In this work, to probe this WC/HG breathing event at atomistic details, as an alternative approach, a fully atomistic computer simulation of free energy landscape was employed


5


Page 6 of 43

for in-depth interpretations of canonical and non-canonical bp pairing modes, transition pathways, and overall structural dynamics of the modified DNA duplex.13-14 Using a state-of-art molecular dynamics (MD) simulation method, we computed free energy profiles of WC/HG breathing fluctuations of both A6-DNA and m1A-A6DNA. For a proper free energy description of the WC/HG bp transition, as with a previous simulation study,13 two collective variables (CVs), glycosidic angle (𝜒) and pseudo-dihedral angle (𝜙), were used (see Figure 1B for the definition). In practice, this kind of free energy simulation is likely to encounter a limited sampling issue, since typical free energy barriers of the WC/HG transition have been previously reported to be 15 kcal/mol at ambient temperature.4-5 Therefore, the presence of such large free energy barriers makes free energy mapping a computationally difficult task. To circumvent this sampling issue by means of accelerating MD sampling, multiple-walker well-tempered metadynamics (mWTMetaD)15-16 was used. Computationally, this mWTMetaD is a more efficient version of WTMetaD for free energy mapping (see Methods). Herein, this mWTMetaD is combined with a modified version of the AMBER bsc1 force field (bsc1_vdW).17 In this bsc1_vdW, only van der Waals (vdW) radii of the polar base


6

Page 7 of 43 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


atoms (N and O) were previously fine-tuned from the standard AMBER bsc1 force field18 under the OPC water model,19 such that the experimental folding free energy of a G-quadruplex DNA was closely reached.17 As a benchmark test of the present protocol (mWTMetaD/bsc1_vdW), the 2-dimensional (2D) (𝜒, 𝜙) free energy simulation was conducted for the normal A6DNA, for which the experimental free energy difference between WC and HG bp (∆𝐺𝑊𝐶→𝐻𝐺) was available for comparison. This test simulation yields a ∆𝐺𝑊𝐶→𝐻𝐺 value (3.5 kcal/mol) in excellent agreement with the NMR value at ambient temperature, thereby well validating this computational scheme. Using the same simulation protocol, the 2D free energy surface was also computed for the damaged DNA (m1A-A6DNA). As expected, this resulting free energy surface at 298 K indicated that the HG bp state was most populated. On the other hand, the canonical WC bp state was distorted into two distorted WC-like bp states (WC* and WC**) in minor populations. According to the present free energy landscape of m1A-A6DNA, the overall transition from these distorted WCs to the HG bp states proceeds via multiple routes in the presence of a wide range of free energy barriers. Moreover, detailed structural perturbations due to the m1A lesion agree well with the NMR experiment.


7


Page 8 of 43

METHODS System setup The two duplex DNA structures were obtained from the NMR structures of A6DNA (PDB ID: 5UZF) and m1A-A6DNA (PDB ID: 5UZI). Each DNA structure was then solvated with 9235 OPC water molecules19 in a cubic box with 66.4 Å sides. For the DNA part, a variant of the AMBER bsc1 force field (bsc1_vdW)17 was used. For the charge neutralization of each system, 22 and 21 Na+ ions were inserted into the systems of A6DNA and m1A-A6DNA, respectively. For the Na+ ion, the Joung-Cheatham parameters optimized for the TIP4PEW20 were used. For N1-methyadenine (m1A), its partial charges were derived using the RED III program21 and the Gaussian 09 program.22 The RESP-A1 method (Connolly surface algorithm using HF/6-31G*)23 was used for the partial charge generation. All partial charges derived for m1A are listed in Table S1. All of the necessary files for the force field implementation are given in Supporting information. The glycosidic angle and other bonding parameters for m1A were taken from the standard AMBER99bsc1 and other missing parameters were borrowed from the AMBER RNA force field database (all_modrna08.frcmod and


8

Page 9 of 43 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


all_modrna08.lib).24 The simulation system was minimized using a steepest descent algorithm, followed by a 10 ns equilibration at 1 atm and 298 K. The pressure was controlled with the Berendsen barostat25 with a coupling constant of 1 ps, and the temperature was maintained with the modified Berendsen thermostat26 with a coupling constant of 0.1 ps. A MD time step of 4 fs was used in conjunction with the hydrogen mass repartitioning.27 The Lennard-Jones (LJ) interaction energy was calculated using a cut-off distance of 10 Å. The particle mesh Ewald (PME) method28 was used for longrange electrostatics interactions with a direct space cut-off of 10 Å. Enhanced sampling MD simulation For effective mapping of the free energy landscape, the mWTMetaD15-16 employs the aforementioned CVs: the glycosidic angle (𝜒) and base flipping angle (𝜙)29 (Figure 1B). For a single target nucleobase (A16 or m1A16), these 𝜒 and 𝜙 angles define the

syn/anti glycosidic rotation and the base flipping from the interior of the duplex DNA into the extra-helical space, respectively. In the framework of the standard WTMetaD, a history dependent bias potential is introduced to guide MD sampling. This external bias is constructed on the fly by a sum of Gaussians deposited on these two CV coordinates.


9


Page 10 of 43

The characteristic nature of this external potential impedes exploration of more frequently visited states, but promotes sampling of unexplored states. As a consequence, the conformational search on the 2D CV surface (𝜒, 𝜙) is accelerated via more frequent barrier-crossing events. In this study, the mWTMetaD was carried out at a constant temperature and pressure of 298 K and 1 atm, respectively. The system pressure was controlled with the Parrinello-Rahman barostat30 with a coupling constant of 2.0 ps. This mWTMetaD is a variant of WTMetaD with multiple replicas (walkers) running in parallel. Since all walkers in the mWTMetaD simultaneously contribute to the time evolution of a single bias potential, a linear scaling performance is achieved in terms of the number of walkers.15 For the present mWTMetaD, the initial height of the Gaussian hill was 0.076 kcal/mol and the width of the Gaussian was set to 5.7° for both CVs. To ensure a well-tempering scheme, the Gaussian Hill with a bias factor of 10 was deposited every 1.0 ps on the CV space on visit. In the mWTMetaD, a total of 20 independent walkers were employed to simultaneously explore the 2D conformational surface. Initial structures for the mWTMetaD run were generated from a 10 ns WTMetaD trajectory. Each of all 20 walkers ran for 2 μs, thereby producing a total of 40


10

Page 11 of 43 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


μs trajectory data combined. From these data, the free energy profile F(𝜒, 𝜙) was constructed via a reweighting process using the variational free energy profile (vFEP) method.31 The reweighting factor for the vFEP was computed using the TiwaryParrinello scheme.32 In this reweighting analysis, the plot of bias offset function C(t) vs. ln (time) indicated that a steady state was reached in 500 ns, as shown by a linear behaviour of C(t) in ln (time) (See Figure S1). Thus, the first 500 ns data per walker were discarded in the reweighting analysis. As shown in Figure S1, the estimated maximum standard deviation of the 1D reweighted free energy profile is about 1.4 kcal/mol. All simulations were performed using the GROMACS program33 in conjunction with the PLUMED software.34 The local kink angle (𝛽h) and global bending calculation Local kinking angles (𝛽h) were calculated using the ABG_calc program,35-37 and global bending angles were calculated using the Curves+ program.38-40 For the computation of average and standard deviation values of global/local kinking angles, a total of 800,000 - 1,500,000 data in the mWTMetaD-generated ensemble were reweighted to produce unbiased results. Details of the local kink angle definition are given in Figure S2. RESULTS


11


Page 12 of 43

Simulation reproduces experimental result in high accuracy. The 2D-free energy surfaces of the WC/HG bp transition for A6DNA and m1A-A6DNA are given in Figure 2A and B. Owing to the additional NMR data available for comparisons, the free energy result of A6DNA would serve as a benchmark test to assess the reliability of the current method. In the case of A6DNA, the general shape of the free energy profile and the presence of multiple pathways were similar to those previously obtained from a different computational protocol.13 In the present study of A6DNA, a noticeable improvement was observed for the free energy difference between the WC and HG bp (∆𝐺𝑊𝐶→𝐻𝐺). Compared to the previous simulation result (∆𝐺𝑊𝐶→𝐻𝐺 = 4.4 kcal/mol),13 the current result (3.5 kcal/mol) very closely reproduces the NMR value (3.3 kcal/mol).4-5 As shown by Figure 2A, various degrees of the base (A16) flipping into the extra-helical space shape up multiple routes connecting the WC to HG bp. In the previous simulation, the corresponding transition state (TS) barriers relative to the AT WC bp state were 10 - 14 kcal/mol. The present study produced a range of 10 - 18 kcal/mol, which was more in line with the NMR value (about 15 kcal/mol) estimated from


12

Page 13 of 43 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


a simple two-state kinetics model.4-5 These improved energetic features may well validate the current simulation protocol for studying DNA duplexes. Free energy surface of m1A-A6DNA discloses multiple pathways with three bp modes In the case of m1A-A6DNA, the three bp modes of HG, WC*, and WC** in the m1A16 T9 bp emerged in the simulated free energy surface (Figure 2B). The m1AT HG bp was unambiguously located as the most populated state (96.7%). But the canonical m1AT WC bp was no longer stable due to inter-base steric repulsion. Instead, the m1A16 at this WC position protruded into the minor and the major grooves to form two WC-like bp modes of WC* (2.8%) and WC** (0.5%), respectively. The occurrence of these two distorted WC pairs appears reminiscent of the wobble bps of TT mismatches.41 Moreover, the current WC* bp mode is similar to a distorted WC bp state of m1GC predicted in an earlier quantum-mechanical study.42 The WC* bp was marginally stable as compared to the WC** bp (∆𝐺𝑊𝐶 ∗ →𝑊𝐶 ∗∗ = 1.7 kcal/mol), and the free energy barrier from WC* to WC** was only 2.3 kcal/mol. Hence, we expected that interconversions between WC* and WC** would possibly be realized with a normal MD simulation. Our 5 μs ordinary MD simulation initiating from a WC* bp state


13


Page 14 of 43

demonstrated that the conformational exchange between these distorted WC took place in a time scale of several hundred nanoseconds (Figure 2C). Of note, this time scale was well-separated from a longer time scale of HG bp conformational exchange estimated from the prior NMR study.11 As shown in the free energy profile (Figure 2B), the 𝜒 angle rotation in either clockwise or anticlockwise direction in concert with the base flipping (𝜙) into either the major- and minor-grooves led to multiple transition pathways connecting the distorted WC-like states (WC* and WC**) to the HG state. Representative pathways and relative energetics of the three bp states and relevant transition state barriers are illustrated in Figure 3. Herein, a total of eight representative pathways connecting the WC* to HG bp state were identified and named as R1, R2, R3, R4, L1, L2, L3, and L4 (Figure 3A). In regard to these names, the first letter R or L denotes the clockwise or counter-clockwise rotation of the 𝜒 angle, respectively. Following this letter, the integer number from 1 to 4 represents an ascending order of free energy barriers relative to the HG bp state. As seen from Figure 3A, the R1 and L1 were relatively low barrier central pathways (with little or no base flipping) and all others were base flipping pathways in various extent into the major- (R4, L4) and minor-groove


14

Page 15 of 43 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(R2, R3, L2, L3). All free energy barriers relative to the HG bp state were in a range of 10 - 18 kcal/mol (Figure 3B). Of note, the larger the base flipping, the higher the free energy barrier. For each representative transition pathway identified, the corresponding movie file (in mp4 format) is provided in Supporting information. Structural changes due to m1A As indicated by the previous NMR study, the introduction of m1A into the A6DNA triggered diverse structural perturbations of the DNA duplex, such as global/local kinks, alteration of backbone and sugar-puckering, and 3'-end bp melting/dynamic changes. In direct comparison with the prior result, we computed the local kink of A6DNA (WC) and m1A-A6DNA (HG) (Figure 4A and B) by reweighting the biased simulation ensembles of the local free energy basin of each bp mode (N = 8×106 – 15×106). The computed local kink angle (𝛽h)35-37 along the junction position was in line with two previous results derived from residual dipolar coupling (RDC) and NOE data (Figure 4A and B). As shown in Figure 4C, the m1A-induced local kink was observed primarily at the m1A16T9 HG bp site and its 3'-neighboring WC bp site (A17T8), and was directed toward the major groove (Figure S2 for the definition of 𝛾h angle and Figure S3 for


15


Page 16 of 43

distribution of 𝛾h angle at m1A16•T9 and A17•T8 junction). For all other junction positions, no substantial change in 𝛽h was observed. At the m1A16T9 HG bp, the 𝛽h value increased by 4.5° relative to that of the A16T9 WC bp (control), but this increase was rather small compared to the RDC-driven value of 8°. Nevertheless, the m1Ainduced change of 𝛽h (∆𝛽h) across the junction bp position was comparable to the RDC result. In the case of A6DNA, the transient A16T9 HG bp was slightly more bent than the canonical A16T9 WC bp (by about 2°). As m1A is incorporated, the local kink at the stable m1A16T9 HG bp further increases by about 3° as compared to that at the unstable A16•T9 HG bp. In terms of the global bending, the global bending angle calculated from the biased ensemble was similarly arranged as m1A16T9 HG (19.1°) > A16T9 HG (16.4°) > A16T9 WC (11.0°) (Figure 4A and B). To help visualize overall structural variations of the DNA duplex for different bp modes, the free energy minimum predicted structures of A16T9 WC, A16T9 HG, and m1A16T9 HG bp state were superimposed in Figure 4D. As previously reported, the m1A lesion affected the BI/BII backbone conformations by perturbing the backbone angles of 𝜀 [C4'(i)-C3'(i)-O3'(i)-P(i+1)] and 𝜁 [C3'(i)-O3'(i)-P(i+1)-


16

Page 17 of 43 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


O5' (i+1)] near the m1A site (see Figure S4 for the angles definition in more details). In a similar manner to the previous analysis,11 the BI/BII conformations can be well resolved by 𝜀 – 𝜁 (BI: -160 ≤ 𝜀 – 𝜁 ≤ +20°, BII: +20 ≤ 𝜀 – 𝜁 ≤ +200°).43 From our biased ensembles for both A6DNA WC and m1A-A6DNA HG, these backbone angles were fully analyzed, and the resulting distributions as a function of 𝜀 − 𝜁 are given in Figure 5. Among eight sites examined, four sites yielded meaningful m1A-induced changes of the BI/BII population; the BI-population increased at G10G11, C15A16, and A16A17, while the BII-population increased only at A17A18. This finding was in agreement with the RDC ensemble result, but no changes were observed for T7T8, T8T9, and T9G10, in contrast to the implication from the RDC data. Besides these m1A-induced backbone changes, in the m1A-A6DNA, the sugar pucker was also shown to be perturbed. The sugar phase angle44 from the biased ensemble was computed via a reweighting scheme and the resulting distribution of sugar phase angle is given in Figure 6. In agreement with the RDC-driven result, the sugar pucker moved towards the C3'-endo (A-conformation, P-P distance 5.9 Å, Sugar phase angle = 0~36°, North) at the m1A site (m1A16 and T9), but the C2'-endo (B-conformation, P-P distance 7.0 Å, Sugar phase angle = 144~180°,


17


Page 18 of 43

South) at T8, G10, and C15. In contrast, the sugar pucker of A17 remains the same (C2'-endo) before and after the mutation. Compared to the P-P distance at A16T9 WC bp (6.8 Å), the P-P distance at m1A16T9 HG bp (6.5 Å) was reduced by 0.3 Å, indicative of a sugar pucker perturbation towards the C3'-endo position. To further justify the validity of the predicted m1A-induced sugar pucker changes (Figure 6), we compared these simulation results to the sugar pucker changes inferred from previous NMR chemical shift data of C3' and C4' (13CΔω) (See Table I).11, 45 Noting that positive (downfield shift) and negative (upfield shift) values of

13CΔω

indicate sugar pucker

changes toward C2'-endo and C3'-endo positions, respectively, the simulated sugar pucker changes were well supported by the NMR chemical shift data. HG bp increases dynamic fluctuations of nucleobase and backbone in the 3'-end It was previously suggested that the HG bp could induce partial melting of the neighboring WC bp at the 3'-end of the m1A.11 Unfortunately, there was no indication of such bp melting in the current ensemble of the HG bp state. Instead, we frequently observed this bp melting in barrier crossing regions of the bp transition pathways, as shown in movie files in Supporting information. This bp melting is more pronounced when


18

Page 19 of 43 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


crossing over large free energy barriers that are associated with more extensive base flipping. As an alternative to the HG-induced partial melting scenario, we looked for possible dynamic changes of the neighboring bp of C15, A16, A17, A18, and A19. The corresponding root mean square fluctuation (RMSF) values were computed separately for nucleobase and backbone (sugar + phosphates) parts (Figure 7). Interestingly, these RMSF results revealed that structural dynamics of the neighboring bp was altered. For the unmodified duplex (A6DNA), the dynamic fluctuations of the nucleobase and backbone of the A16 residue were most pronounced, but upon the mutation, the fluctuations of m1A16 were rather suppressed, but the fluctuations of its neighboring three adenine bases (A17 – A19) were enhanced as compared to those of m1A16. DISCUSSION To investigate the detailed structure and dynamics of WC/HG bp breathing, we performed an enhanced sampling MD simulation to predict the free energy profiles of the WC/HG bp transition of an unmodified duplex DNA (A6DNA) and its mutated DNA with a m1A lesion (m1A-A6DNA). As shown by the benchmark test on A6DNA, the mWTMetaD in conjunction with the bsc1_vdW/OPC force field enabled a reliable and


19


Page 20 of 43

accurate construction of the free energy profiles. Especially, compared to a prior simulation study, a substantial improvement was achieved for the free energy difference between the WC and HG bp states of A6DNA.13 In part, this improvement was ascribed to the use of the modified version of AMBER bsc1 force field (bsc1_vdW) in conjunction with the OPC water model. The main added feature in the bsc1_vdW was that the strength of bp H-bonds in the original bsc1 was enhanced effectively using a simple local correction of vdW radii.17 Notably, the current choice of the force field reinforces the importance of proper refinements of bp H-bond interactions in the standard AMBER bsc1 force field. Except for the distortion of the canonical WC bp and the population inversion of WC-like vs. HG bp state upon the N1-methylation, the 2D free energy profile of m1A-A6DNA shares similar topographic features with that of A6-DNA. Both free energy surfaces of the unmodified and modified DNA exhibited multiple transition pathways with a wide range of free energy barriers (10−18 kcal/mol) of the WC/HG bp transition. As mentioned above, the magnitude of these transition barriers depends on how much the target base is opened to extra-helical space during the transition process. The presence


20

Page 21 of 43 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


of multiple pathways undertaking various extents of base flipping implies that WC/HG breathing dynamics of both A6-DNA and m1A-A6DNA is governed by rather a complex kinetics beyond a simple two-state model. It was previously indicated that the stable m1AT HG bp and the transient AT HG bp show similar chemical shifts and HG H-bond type.11 But the absence of NMR structures of the transient AT HG bp prevented direct structural comparisons of these two HG bps. Since both bp structures became available from the present simulation, the two HG bp can be directly compared. Despite the aforementioned similarity of these HG bps, the simulated structure of m1A16T9 HG bp somewhat deviated from that of the A16T9 HG, as indicated by the root mean square deviation (RMSD) of 1.51 Å. This study finds that a major source of this structural deviation is the DNA kink angles. According to the global/local kink angles computed for these HG bps, apparently the m1AT HG bp state was more kinked than the AT HG bp state. Therefore, the m1A lesion facilitates bending of the DNA duplex by adopting a more stable version of HG bp and thereby providing possibly better access to the m1A site by repair-proteins.


21


Page 22 of 43

Notably, the present study produced the m1A-induced perturbation effects of the backbone and sugar pucker in reasonable agreement with the NMR experiment. As another interesting structural change due to the m1A lesion, the formation of HG bp was previously suggested to induce partial melting of neighboring WC bp at the 3'-end of the m1A.11 In contrast, our simulation could not find any evidence of such HG bp-induced partial bp melting. This absence may result from force field artifacts, but such force field issues need to be further investigated. Instead of the HG-induced partial bp melting, we find that the formation of HG bp launches dynamic changes of the neighboring WC bp at the 3'-end of m1A. In line with earlier observations,46-47 this change of residue dynamics at the 3'-end of m1A provides an interesting molecular insight into recognition process by repair proteins. All of the m1A-induced structural alterations mentioned here could be related to the creation of a “soft spot” for recognition by repair-enzymes.10 As with the HG bp state, the WC* bp state features substantial local bending of the duplex around and in the m1A16 (Figure S5). Moreover, the H-bond strength of the WC* bp is weakened (only a


22

Page 23 of 43 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


single bp H-bond: see Figure 2B). Therefore, we cannot rule out a possibility that this minor bp state can also be detected by repair-enzymes such as AlkB or ABH2.8, 48 In previous high-resolution NMR studies, 1 μs MD-generated ensemble obtained from the AMBER99 bsc049 and the SPC/E50 were sorted out to obey experimental RDC data.11, 51-53 In this study, however, no references to any experimental data were used in the reweighting analysis of the mWTMetaD ensemble. Therefore, it is encouraging that the pure simulation ensembles (a total of 40 μs trajectory) mimic most of the key structural features of the DNA duplexes probed from previous NMR experiments. To further improve the predicting power of the present simulation, one can use recently developed MD schemes by systematically combining experimental ensemble averages.54-56 Interestingly, these experimental guided MD techniques can also be used to improve empirical force fields.56 Overall, we expect that the current protocol (mWTMetaD/bsc1_vdW/OPC) can be a robust computational method for complementing and guiding ongoing DNA breathing experiments. Also, this scheme can be useful for tackling many other outstanding issues of DNA mismatches, DNA with rare bps, and DNA/protein interactions.


23


Page 24 of 43

ASSOCIATED CONTENT

Supporting Information. The following files are available free of charge on the ACS publications website.

Table S1, Figures S1-S5, and cluster analysis method details, force field parameter files for m1A used in this study (PDF)

Movies for transition pathways (mp4 format)

AUTHOR INFORMATION

Corresponding Author * Tel: +82-51-510-2235. Fax: +82-51-516-7421. Email: [email protected]

Funding Sources

This work was supported by the National Research Foundation [NRF2014R1A4A1001690,

NRF-2016R1A2B4009517] ACKNOWLEDGMENT


24

Page 25 of 43 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


C.Y. thanks Prof. Hashim Al-Hashimi and Mr. Honglue Shi for helpful comments on the ABG_calc code.

REFERENCES 1.

Frankkamenetskii, M., DNA Chemistry - How the Double Helix Breathes. Nature

1987, 328, 17-18. 2.

von Hippel, P. H.; Johnson, N. P.; Marcus, A. H., Fifty Years of DNA "Breathing":

Reflections on Old and New Approaches. Biopolymers 2013, 99, 923-954. 3.

Frank-Kamenetskii, M. D., DNA breathes Hoogsteen. Artif. DNA PNA XNA 2011,

2, 1-3. 4.

Nikolova, E. N.; Kim, E.; Wise, A. A.; O’Brien, P. J.; Andricioaei, I.; Al-Hashimi, H.

M., Transient Hoogsteen base pairs in canonical duplex DNA. Nature 2011, 470, 498502. 5.

Alvey, H. S.; Gottardo, F. L.; Nikolova, E. N.; Al-Hashimi, H. M., Widespread

transient Hoogsteen base pairs in canonical duplex DNA with variable energetics. Nat.

Commun. 2014, 5, 4786 6.

Watson, J. D.; Crick, F. H., Molecular structure of nucleic acids; a structure for

deoxyribose nucleic acid. Nature 1953, 171, 737-738.


25


7.

Page 26 of 43

Hoogsteen, K. R., The crystal and molecular structure of a hydrogen‐bonded

complex between 1‐methylthymine and 9‐methyladenine. Acta Crystallographica 1963,

16, 907-916. 8.

Yang, C. G.; Yi, C. Q.; Duguid, E. M.; Sullivan, C. T.; Jian, X.; Rice, P. A.; He, C.,

Crystal structures of DNA/RNA repair enzymes AlkB and ABH2 bound to dsDNA.

Nature 2008, 452, 961-965. 9.

Lu, L. H.; Yi, C. Q.; Jian, X.; Zheng, G. Q.; He, C. A., Structure determination of

DNA methylation lesions N-1-meA and N-3-meC in duplex DNA using a cross-linked protein-DNA system. Nucleic Acids Res. 2010, 38, 4415-4425. 10. Yang, C. G.; Garcia, K.; He, C., Damage Detection and Base Flipping in Direct DNA Alkylation Repair. Chembiochem 2009, 10, 417-423. 11. Sathyamoorthy, B.; Shi, H. L.; Zhou, H. Q.; Xue, Y.; Rangadurai, A.; Merriman, D. K.; Al-Hashimi, H. M., Insights into Watson-Crick/Hoogsteen breathing dynamics and damage repair from the solution structure and dynamic ensemble of DNA duplexes containing m1A. Nucleic Acids Res. 2017, 45, 5586-5601. 12. Yang, H.; Zhan, Y. Q.; Fenn, D.; Chi, L. M.; Lam, S. L., Effect of 1-methyladenine on double-helical DNA structures. FEBS Letters 2008, 582, 1629-1633.


26

Page 27 of 43 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


13. Yang, C.; Kim, E.; Pak, Y., Free energy landscape and transition pathways from Watson-Crick to Hoogsteen base pairing in free duplex DNA. Nucleic Acids Res. 2015,

43, 7769-7778. 14. Chakraborty, D.; Wales, D. J., Energy Landscape and Pathways for Transitions between Watson-Crick and Hoogsteen Base Pairing in DNA. J. Phys. Chem. Lett. 2018,

9, 229-241. 15. Raiteri, P.; Laio, A.; Gervasio, F. L.; Micheletti, C.; Parrinello, M., Efficient reconstruction of complex free energy landscapes by multiple walkers metadynamics. J.

Phys. Chem. B 2006, 110, 3533-3539. 16. Barducci, A.; Bussi, G.; Parrinello, M., Well-tempered metadynamics: A smoothly converging and tunable free-energy method. Phys. Rev. Lett. 2008, 100, 020603. 17. Yang, C.; Kulkarni, M.; Lim, M.; Pak, Y., In silico direct folding of thrombinbinding aptamer G-quadruplex at all-atom level. Nucleic Acids Res 2017, 45, 1264812656. 18. Ivani, I.; Dans, P. D.; Noy, A.; Perez, A.; Faustino, I.; Hospital, A.; Walther, J.; Andrio, P.; Goni, R.; Balaceanu, A.; Portella, G.; Battistini, F.; Gelpi, J. L.; Gonzalez, C.; Vendruscolo, M.; Laughton, C. A.; Harris, S. A.; Case, D. A.; Orozco, M., Parmbsc1: a refined force field for DNA simulations. Nat. Methods 2016, 13, 55-58.


27


Page 28 of 43

19. Izadi, S.; Anandakrishnan, R.; Onufriev, A. V., Building Water Models: A Different Approach. J. Phys. Chem. Lett. 2014, 5, 3863-3871. 20. Joung, I. S.; Cheatham, T. E., Determination of alkali and halide monovalent ion parameters for use in explicitly solvated biomolecular simulations. J. Phys. Chem. B 2008, 112, 9020-9041. 21. Dupradeau, F. Y.; Pigache, A.; Zaffran, T.; Savineau, C.; Lelong, R.; Grivel, N.; Lelong, D.; Rosanski, W.; Cieplak, P., The R.E.D. tools: advances in RESP and ESP charge derivation and force field library building. Phys. Chem. Chem. Phys. 2010, 12, 7821-7839. 22. Frisch, M. J.; Trucks, G.; Schlegel, H. B.; Scuseria, G.; Robb, M.; Cheeseman, J.; Scalmani, G.; Barone, V.; Mennucci, B.; Petersson, G., Gaussian 09, revision A. 1.

Gaussian Inc. Wallingford CT 2009, 27, 34. 23. Cornell, W. D.; Cieplak, P.; Bayly, C. I.; Gould, I. R.; Merz, K. M.; Ferguson, D. M.; Spellmeyer, D. C.; Fox, T.; Caldwell, J. W.; Kollman, P. A., A second generation force field for the simulation of proteins, nucleic acids, and organic molecules. J. Am.

Chem. Soc. 1995, 117, 5179-5197. 24. Aduri, R.; Psciuk, B. T.; Saro, P.; Taniga, H.; Schlegel, H. B.; SantaLucia, J., AMBER force field parameters for the naturally occurring modified nucleosides in RNA.

J. Chem. Theory Comput. 2007, 3, 1464-1475.


28

Page 29 of 43 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


25. Berendsen, H. J. C.; Postma, J. P. M.; Vangunsteren, W. F.; Dinola, A.; Haak, J. R., Molecular-Dynamics with Coupling to an External Bath. J. Chem. Phys. 1984, 81, 3684-3690. 26. Bussi, G.; Donadio, D.; Parrinello, M., Canonical sampling through velocity rescaling. J. Chem. Phys. 2007, 126, 014101. 27. Feenstra, K. A.; Hess, B.; Berendsen, H. J. C., Improving efficiency of large timescale molecular dynamics simulations of hydrogen-rich systems. J. Comput. Chem. 1999, 20, 786-798. 28. Darden, T.; York, D.; Pedersen, L., Particle Mesh Ewald - an N•Log(N) Method for Ewald Sums in Large Systems. J. Chem. Phys. 1993, 98, 10089-10092. 29. Song, K.; Campbell, A. J.; Bergonzo, C.; de los Santos, C.; Grollman, A. P.; Simmerling, C., An Improved Reaction Coordinate for Nucleic Acid Base Flipping Studies. J. Chem. Theory Comput. 2009, 5, 3105-3113. 30. Parrinello, M.; Rahman, A., Polymorphic Transitions in Single-Crystals - a New Molecular-Dynamics Method. J. Appl. Phys. 1981, 52, 7182-7190. 31. Lee, T. S.; Radak, B. K.; Pabis, A.; York, D. M., A New Maximum Likelihood Approach for Free Energy Profile Construction from Molecular Simulations. J. Chem.

Theory Comput. 2013, 9, 153-164.


29


Page 30 of 43

32. Tiwary, P.; Parrinello, M., A Time-Independent Free Energy Estimator for Metadynamics. J. Phys. Chem. B 2015, 119, 736-742. 33. Abraham, M. J.; Murtola, T.; Schulz, R.; Páll, S.; Smith, J. C.; Hess, B.; Lindahl, E., GROMACS: High performance molecular simulations through multi-level parallelism from laptops to supercomputers. SoftwareX 2015, 1, 19-25. 34. Tribello, G. A.; Bonomi, M.; Branduardi, D.; Camilloni, C.; Bussi, G., PLUMED 2: New feathers for an old bird. Comput. Phys. Commun. 2014, 185, 604-613. 35. Bailor, M. H.; Sun, X. Y.; Al-Hashimi, H. M., Topology Links RNA Secondary Structure with Global Conformation, Dynamics, and Adaptation. Science 2010, 327, 202-206. 36. Bailor, M. H.; Mustoe, A. M.; Brooks, C. L.; Al-Hashimi, H. M., 3D maps of RNA interhelical junctions. Nat. Protoc. 2011, 6, 1536-1545. 37. Zhou, H. Q.; Hintze, B. J.; Kimsey, I. J.; Sathyamoorthy, B.; Yang, S.; Richardson, J. S.; Al-Hashimi, H. M., New insights into Hoogsteen base pairs in DNA duplexes from a structure-based survey. Nucleic Acids Res. 2015, 43, 3420-3433. 38. Blanchet, C.; Pasi, M.; Zakrzewska, K.; Lavery, R., CURVES+ web server for analyzing and visualizing the helical, backbone and groove parameters of nucleic acid structures. Nucleic Acids Res. 2011, 39, W68-W73.


30

Page 31 of 43 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


39. Lavery, R.; Sklenar, H., Defining the Structure of Irregular Nucleic-Acids Conventions and Principles. J. Biomol. Struct. Dyn. 1989, 6, 655-667. 40. Lavery, R.; Moakher, M.; Maddocks, J. H.; Petkeviciute, D.; Zakrzewska, K., Conformational analysis of nucleic acids revisited: Curves+. Nucleic Acids Res. 2009,

37, 5917-5929. 41. He, G.; Kwok, C. K.; Lam, S. L., Preferential base pairing modes of T.T mismatches. FEBS Letters 2011, 585, 3953-3958. 42. Felske, L. R.; Lenz, S. A.; Wetmore, S. D., Quantum Chemical Studies of the Structure and Stability of N-Methylated DNA Nucleobase Dimers: Insights into the Mutagenic Base Pairing of Damaged DNA. J. Phys. Chem. A 2017, 122, 410-419. 43. Hartmann, B.; Piazzola, D.; Lavery, R., BI-BII transitions in B-DNA. Nucleic Acids

Res. 1993, 21, 561-568. 44. Westhof, E.; Sundaralingam, M., A method for the analysis of puckering disorder in five-membered rings: the relative mobilities of furanose and proline rings and their effects on polynucleotide and polypeptide backbone flexibility. J. Am. Chem. Soc. 1983,

105, 970-976. 45. Shi, H.; Clay, M. C.; Rangadurai, A.; Sathyamoorthy, B.; Case, D. A.; Al-Hashimi, H. M., Atomic structures of excited state A–T Hoogsteen base pairs in duplex DNA by


31


Page 32 of 43

combining NMR relaxation dispersion, mutagenesis, and chemical shift calculations. J.

Biomol. NMR 2018, 70, 229-244. 46. Nag, N.; Rao, B. J.; Krishnamoorthy, G., Altered dynamics of DNA bases adjacent to a mismatch: A cue for mismatch recognition by MutS. J. Mol. Biol. 2007,

374, 39-53. 47. Yang, W., Surviving the sun: Repair and bypass of DNA UV lesions. Protein Sci. 2011, 20, 1781-1789. 48. Yi, C. Q.; Chen, B. E.; Qi, B.; Zhang, W.; Jia, G. F.; Zhang, L.; Li, C. J.; Dinner, A. R.; Yang, C. G.; He, C., Duplex interrogation by a direct DNA repair protein in search of base damage. Nat. Struct. Mol. Biol. 2012, 19, 671- 676. 49. Perez, A.; Marchan, I.; Svozil, D.; Sponer, J.; Cheatham, T. E., 3rd; Laughton, C. A.; Orozco, M., Refinement of the AMBER force field for nucleic acids: improving the description of alpha/gamma conformers. Biophys. J. 2007, 92, 3817-3829. 50. Berendsen, H. J. C.; Grigera, J. R.; Straatsma, T. P., The Missing Term in Effective Pair Potentials. J. Phys. Chem. 1987, 91, 6269-6271. 51. Chen, Y.; Campbell, S. L.; Dokholyan, N. V., Deciphering protein dynamics from NMR data using explicit structure sampling and selection. Biophys. J. 2007, 93, 23002306.


32

Page 33 of 43 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


52. Frank, A. T.; Stelzer, A. C.; Al-Hashimi, H. M.; Andricioaei, I., Constructing RNA dynamical ensembles by combining MD and motionally decoupled NMR RDCs: new insights into RNA dynamics and adaptive ligand recognition. Nucleic Acids Res. 2009,

37, 3670-3679. 53. Salmon, L.; Bascom, G.; Andricioaei, I.; Al-Hashimi, H. M., A General Method for Constructing Atomic-Resolution RNA Ensembles using NMR Residual Dipolar Couplings: The Basis for Interhelical Motions Revealed. J. Am. Chem. Soc. 2013, 135, 5457-5466. 54. Camilloni, C.; Cavalli, A.; Vendruscolo, M., Replica-Averaged Metadynamics. J.

Chem. Theory Comput. 2013, 9, 5610-5617. 55. White, A. D.; Voth, G. A., Efficient and Minimal Method to Bias Molecular Simulations with Experimental Data. J. Chem. Theory Comput. 2014, 10, 3023-3030. 56. Cesari, A.; Gil-Ley, A.; Bussi, G., Combining Simulations and Solution Experiments as a Paradigm for RNA Force Field Refinement. J. Chem. Theory Comput. 2016, 12, 6192-6200.


33


This work

T8

Expt. 1

Expt. 2

Sathyamoorthy et al. 11

Shi et al. 45

shift toward C2'-endo

13CΔω

shift toward C3'-endo

13CΔω

(C4') = -0.5 ppm

G10 shift toward C2'-endo

13CΔω

(C4') = -0.2ppm

13CΔω

(C4') = 0.1 ppm

T9

A17

no change from C2'-endo

Page 34 of 43

A16 shift toward C3'-endo

13CΔω

C15 shift toward C2'-endo

13CΔω

(C4') = 0.3 ppm

(C4') = -0.7 ppm

(C4') = 0.5 ppm

13CΔω

(C3') = 0.5 ppm

13CΔω

(C4') = 0.3 ppm

13CΔω

(C3') = -1.7 ppm

13CΔω

(C4') = -0.4 ppm

13CΔω

(C3') = 0.6 ppm

13CΔω

(C4') = 0.4 ppm

13CΔω

(C3') = 0.6 ppm

13CΔω

(C4') = 0.4 ppm

13CΔω

(C3') = -3.5 ppm

13CΔω

(C4') = -1.3 ppm

13CΔω

(C3') = 2.5 ppm

13CΔω

(C4') = 0.6 ppm


34

Page 35 of 43 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Table 1. Comparison of m1A-induced sugar pucker changes from the current simulation and previous NMR chemical shift data. Chemical shift changes of C3'/C4' (13CΔω) show directions of sugar pucker due to the m1A mutation. Downfield (positive upfield (negative

13CΔω)

13CΔω)

and

shifts indicate m1A-induced sugar pucker change toward C2'-

endo and C3'-endo positions, respectively. For the definition of C2' and C3' atom types, see Figure S4.


35


Figure 1. (A) DNA sequences of A6DNA and m1A-A6DNA with the methylation site highlighted blue. Shown is the chemical structure of the m1A base (N1-methylated adenine) with its methyl group colored red. (B) Definition of collective variables used in this study: glycosidic angle (0 ≤ χ ≤180°) and base flipping angle (180 ≤ Φ ≤ 180°). The χ is a torsional angle consisting of O4'-C1'-N9-C4 atoms. The Φ is a torsional angle defined by P1-P2-P3-P4. P1 is the center of mass (CM) of heavy atoms of the upper and lower bp enclosed by red rectangles. P2 or P3 is the CM of the phosphate group and P4 is the CM of the 5-membered ring of A16 or m1A16. A positive or negative Φ value indicates base flipping into major- or minor-grooves. 73x64mm (300 x 300 DPI)


Page 36 of 43

Page 37 of 43 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure 2. Reweighted free energy profiles on (χ, Φ) for WC/HG bp transition. (A) A6DNA and (B) m1AA6DNA. Free energy differences between bp states are shown alongside the arrow. Simulated (black) and experimental (red) populations of each bp state are denoted in parenthesis. (C) Conformational exchange between the WC* and WC** bp states as shown by 5 μs MD simulation initiating from the WC* state. While the WC* bp is a small minor groove opening state (with a small negative Φ), the WC** bp is a major groove opening state (with a positive Φ angle). Inter-conversion between the WC* and WC** bp occurs within the first hundred nanoseconds. Schematic free energy profile governing this conformational exchange is shown in inset. 150x232mm (300 x 300 DPI)



Figure 3. (A) Detailed bp transition pathways shown on the 2D-free energy map. Clockwise χ rotation pathways are named R1, R2, R3, and R4. Anti-clockwise χ rotation are named L1, L2, L3 and L4. The integer numbers 1 – 4 indicate an ascending order of free energy barriers in the pathways. (B) Schematic free energy level diagram of bp states and transition states along each pathway. The transition state between WC* and WC** is coloured orange. 170x118mm (300 x 300 DPI)


Page 38 of 43

Page 39 of 43 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure 4. Local kinking angle (βh) as a function bp junction positions and global kinking angle calculated from the mWTMetaD-generated biased ensemble (A) A6DNA (A16•T9 WC bp and A16•T9 HG bp) and (B) m1A-A6DNA (mA16•T9 HG bp). These simulated values are compared with previous results inferred from NMR ensemble (XPLOR) and RDC-selected ensemble (RDC).11 (C) Local kinking angle changes (Δβh) as a function of bp junction positions from A6DNA (A16•T9 WC) to A6DNA (A16•T9 HG) and from A6DNA (A16•T9 WC) to m1A-A6DNA (mA16•T9 HG). Simulated values are compared to values driven from the RDC ensembles. (D) Overlapped structures of A16•T9 WC (grey), A16•T9 HG (red), and m1A16•T9 HG (blue). Each structure was taken from the corresponding local minimum of the free energy profiles of A6DNA and m1A-A6DNA. The heavy atom root mean square deviation (RMSD) values of A16•T9 WC and m1A16•T9 HG are 1.19 and 1.25 Å, respectively, relative to their reported NMR structures.11 The RMSD value between m1A16•T9 HG and A16•T9 HG is 1.51 Å. 137x104mm (300 x 300 DPI)



Figure 5. Reweighted distribution of (ε – ζ) calculated from biased ensemble for each bp of A6DNA (A16•T9 WC: red) and m1A-A6DNA (mA16•T9 HG: blue). Arrows indicate the direction of the backbone change (BI ⟷ BII) upon N1-methylation. The BI conformation is defined in the range of -160 ≤ ε – ζ ≤ +20° and the BII conformation is defined in the range of +20 ≤ ε – ζ ≤ +200°.43 230x324mm (300 x 300 DPI)


Page 40 of 43

Page 41 of 43 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure 6. Reweighted distributions of sugar phase angles44 calculated from biased ensemble for bps around the methylation site. Resulting distributions for A6DNA (A16•T9 WC: red) and m1A-A6DNA (mA16•T9 HG: blue). Arrows indicate the direction of sugar pucker change (C2'-endo ⟷ C3'-endo) upon N1-methylation. 175x193mm (300 x 300 DPI)



Figure 7. Root mean square fluctuation (RMSF) values of C15, A16, A17, A18, and A19 for A6DNA (A16•T9 WC) and m1A-A6DNA (mA16•T9 HG). (A) Base fluctuation (Å) and (B) Backbone fluctuation (Å). For each RMSF calculation of the biased ensembles of A6DNA and m1A-A6DNA, the most populated cluster-center structure was chosen as a reference structure. Details of the clustering method used in this analysis are given in Supporting information. 96x55mm (300 x 300 DPI)


Page 42 of 43

Page 43 of 43 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Table of Contents 35x18mm (300 x 300 DPI)


Hoogsteen breathing in a

Recommend Documents