Disagreement Between the Structure of the dTpT Thymine Pair

Article pubs.acs.org/JPCB

Disagreement Between the Structure of the dTpT Thymine Pair Determined by NMR and Molecular Dynamics Simulations Using Amber 14 Force Fields Collins Nganou,† Scott D. Kennedy,‡ and David W. McCamant†,* †

Department of Chemistry, University of Rochester, Rochester, New York 14627, United States Department of Biochemistry and Biophysics, University of Rochester, Rochester, New York 14642, United States

‡

S Supporting Information *

ABSTRACT: We report a disagreement between the predicted structures of the dTpT thymine pair (thymidylyl(3′ → 5′)thymidine) using nuclear magnetic resonance (NMR) spectroscopy and molecular dynamics (MD) simulations using the AMBER ff14SB and ff14 + ε/ζOL1 + χOL4 force fields for DNA. The NMR structure was determined using NOE couplings to thymine’s H6 and JHH couplings between sugar protons. The MD simulation used replica exchange methods to produce converged statistics in a 500 ns trajectory. NMR data indicate that both thymine nucleotides in the pair display an anti conformation of B-DNA, while the MD simulations predict a structure in which the 5′-thymine is flipped into a syn conformation and the 3′-thymine is in an anti conformation. The syn conformation of the 5′-thymine predicted by MD appears by a ∼ 180-deg flip of the glycosidic angle in comparison to the B-form anti structure. Differences in the distortion of the sugar pucker between 5′-thymine and 3′-thymine further highlighted the surprisingly different conformation of the 5′- and 3′-ends. While both MD and NMR indicate the deoxyribose sugars to be primarily in the 2′-endo conformation typical of B-form DNA, the MD simulations predict a more twisted conformation (2′-endo/1′-exo) for the 5′sugar and significant flexibility of C3′ of the 3′-sugar. We conclude that the current AMBER force field does not accurately predict the conformation of single-stranded thymine, in agreement with previous work investigating single-stranded DNA. reaction.5 In other words, the low quantum-yield is really a measure of the low probability that the reactive bonds are in close proximity at the moment a UV photon is absorbed. The hypothesis that base stacking geometry in the ground state influences the excited state behavior of DNA has triggered an intensive collection of computational studies of DNA.7,9−12 Based on MD simulations, Kohler and co-workers hypothesized that two parameters need to be measured to predict the QY of the (2 + 2) cycloaddition reaction in TT: the C5C6:C5C6 adjacent bond distance, d, and the dihedral angle of those bonds, η.7 Utilizing 25 6 ns trajectories of dTpT in the AMBER-99 force field, they found a dominant probability distribution of ∼4 Å and a dihedral angle of ∼30°.7 Further, they found that the QY can be predicted by integrating the probability of forming conformations with d < 3.63 Å and |η| < 48°. Later, McCullagh et al. examined short double-stranded hairpin DNA structures with ten 6 ns production runs using the CHARMM27 force field.10 They showed that only one parameter, the distance between adjacent C5C6 bonds, was sufficient to determine the probability of CPD production and that the reactive conformations were those with d < 3.52

1. INTRODUCTION In this report, we discuss the structure and structural fluctuations of a single-stranded thymine pair (thymidylyl(3′ → 5′)thymidine or dTpT). The dTpT DNA sequence is the minimum strand necessary to simulate the damaging effects of ultraviolet light on DNA, in particular the formation of cyclobutane pyrimidine dimers (CPD) and pyrimidinepyrimidone dimers (6−4 adduct) shown in Figure 1. Interest in the UV photochemistry of DNA has long been motivated by interest in the events leading to photocarcinogenesis.1−3 Thymine dimerization generates two potential photolesions: a CPD formed by a 2 + 2 cycloaddition reaction of the adjacent C5C6 double bonds and a pyrimidine-pyrimidone formed by a 6−4 cycloaddition reaction between the CO double bond of the 3′-thymine (3′-end) and the CC double bond of the 5′-thymine (5′-end).4 Recently, researchers have used ultrafast spectroscopy to show that the dimerization of the thymine pair occurs in less than 1 ps following absorption of UV light.5 Other steady-state experiments have quantified the dimerization quantum yield (QY) as ∼0.01. 6−8 The rapidity of the excited-state dimerization reaction and its low yield suggest that the yield may be controlled by the probability of finding the ground-state DNA in a conformation that promotes the dimerization © XXXX American Chemical Society

Received: January 7, 2016

A

DOI: 10.1021/acs.jpcb.6b00191 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B

Figure 2. Schematic diagrams of (a) the glycosidic rotation, χ, and (b) the sugar pucker in nucleic acids. In part a, the two stable minima are shown at χ = −120°, termed “anti”, and χ = 60°, termed “syn”. Ribose schematics show the possible deoxyribose sugar puckers in both 3′endo and 2′-endo DNA. Images in part b use the X-ray crystallography parameters in NAB (Ambertools 2014) and 3DNA and (a) follows the IUPAC drawing scheme.22,23 In part b, oxygen atoms are red, nitrogens of the nucleic acid are blue and carbons are gray; hydrogens are not shown.

NMR is also uniquely capable of characterizing the solution conformation of the backbone ribose in DNA.17 The sugar puckering can result in a variety of out-of-plane deformations of the sugar ring. These sugar conformations are elegantly described by the sugar’s pseudorotation angle, P, and puckering amplitude, νm, first introduced by Altona and Sundaralingam and reviewed in multiple sources.18−21 In nucleic acids, these are typically constrained to two broad categories as shown in Figure 2b. In B-form DNA, the deformation is a 2′-endo conformation with only the C2′ out of plane, puckered up toward the nucleic acid and characterized by a pseudorotation angle from 120 to 180°.20 In A-form DNA, the sugar is in a 3′endo conformation, with the C3′ carbon tilted up toward the nucleic acid and a pseudorotation angle from 0 to 30°. Recent optimization of DNA force fields have focused on the potential energy for rotation about the sugar-base bond, χ, and the backbone sugar−phosphate dihedral angles. Here, we note that although dTpT is small enough to be treated by small molecule force fields, the specific interactions that make it unique have been optimized significantly from the foundational small molecule forces, so that the latest models perform much better at predicting both the structure of the monomeric nucleotides and double stranded DNA. The Amber force fields are reviewed concisely and clearly in the most recent user manual but will be reviewed here as well.24 The ff99 force field improved upon the original force fields for nucleic acids by improving the atomic charge parameters and, around the same time, the handling of long-range electrostatic interactions was significantly improved.25 These parameters were further optimized to correct the α/γ dihedral angle potentials in the “bsc0” improvements in 2008,26,27 making the “ff14SB” force field, that has been broadly and successfully used in a variety of nucleic acid studies. About the same time, concerns arose about the glycosidic potential energy, which governs the base flipping about the nucleic acid-sugar bond, i.e. the χ angle. The rotation about χ in ff99 was optimized to ff99χ using quantum chemical calculations by Ode et al. but these parameters were not implemented in a broadly distributed version of Amber.28 Later

Figure 1. Photoreaction of two stacked thymines causes two possible cycloadducts ((2 + 2) and (6−4)) after 250−290 nm irradiation.

Å.10 However, in 2013 Cheatham and co-workers showed that single stranded RNA tetramers are poorly converged even after >5 μs of direct MD simulation time and that structural convergence only occurred after extensive replica exchange MD was performed.13,14 Hence, the validity of the previous dTpT simulations are called into question since the sampled structures are biased toward the initial, B-form structure, and not the equilibrium structure. Although some benchmarking of double-stranded nucleic acid structures has occurred,15,16 there has been little to no benchmarking of the dTpT structures predicted by MD simulations to experimental measures of solution phase structure. NMR is particularly well suited for generating testable predictions of structure. For instance, the work of Chattopadhyaya and co-workers on the nuclear Overhauser effect (NOE) of single-stranded DNA (ssDNA) in solution has revealed a high population of ssDNA with a right-handed helical structure, which is evocative of the right-handed helical duplex, the B-DNA form of DNA.17 Modeling the NOE from the overall MD simulation and comparing it with NMR data can provide an effective test that shows how well the current state of the force field can predict the conformation of the thymine pair. In particular, the NOE coupling to the H6 of thymine can measure the configuration of the base-flipping motion as the nucleic acid rotates about the glycosidic bond (see Figure 2a). Rotation about the glycosidic bond is characterized by the dihedral angle χ with two stable minima at χ = −120°, termed “anti”, and χ = 60°, termed “syn”. Stable B-form pyrimidines predominately adopting the anti configuration. B


Article

The Journal of Physical Chemistry B studies discovered that separate χ optimizations were required for DNA and RNA.29−31 In their effort to improve the accuracy of the ff14SB force field (i.e., ff99 + bsc0), Sponer and co-workers optimized the χ parameters for DNA specifically, generating the χοL4 parameter set.32 The optimization compared energies of high-level quantum mechanical calculations with energies from the classical force fields in which both methods used implicit solvent models, a fixed χ angle and relaxed structures during the torsional scan.29 Relative to ff99, the new parameters of χOL4 exhibits a smaller barrier at χ = ∼120°, as the C2O of thymine passes over C2′ of the sugar, and a much higher barrier at χ = 0°, as the C2O passes over the ribose oxygen, O4′. Additionally, χOL4 reduces the energy difference between the syn minimum (χ = ∼ 60°) and the anti minimum (χ = 190− 250°), making them roughly isoenergetic in deoxyguanosine.32 The DNA backbone force field was also recently updated by Zgarbova et al., who improved the ε (C4′−C3′−O3′−P dihedral) and ξ (C3′−O3′−P−O5′ dihedral) using a MP2/ CBS quantum mechanical simulation.33 The refinement was designed to correct the relative populations of structural subensembles in duplex DNA that were incorrectly predicted in the ff14SB force field. They showed that the optimized force field improved the agreement between the MD simulated structure from a ∼1-μs run and X-ray and solution NMR structures of B-form duplex DNA. Importantly, all of the new refinements of the AMBER force field (ff14 + e/zOL133 + qOL432) were optimized and benchmarked with double stranded DNA.33 In this work, we will compare the ff14SB force field to the latest force field incorporating the bsc0, χ oL4, and ε and ξ improvements (ε/ξOL1), which we term “ff14SB/ chi/ez”.24 Because these force fields have been used to predict the proximity of photoreactive thymines to each other, we must test whether the current AMBER force field can indeed predict the structure of ssDNA.7,10,34−36 To achieve this goal, here we compare the conformation of dTpT obtained from solution NMR with that predicted from MD simulations using the ff14SB and ff14SB/chi/ez force field parameters. Because of the conformational flexibility of short single-stranded DNA strands, we have utilized replica exchange MD (REMD37,38) and long simulation times (500 ns) to ensure accurate sampling of the dTpT structure during the simulations.

NOE cross-peak intensities by assuming 1/r6 scaling and using intraresidue H1′:H4′ cross-peaks as reference intensities. The H1′ to H4′ distance was taken to be 3.2 Å in the C2′-endo conformation and 3.4 Å in the C3′-endo conformation. A scaling factor was calculated from these distances and sugarpucker populations of each residue determined from spin−spin couplings (see above). Populations of anti/syn conformations of the glycosidic bonds were determined from measured H6:H1′ NOE cross-peak intensities and assuming that this distance is 3.7 Å in the anti conformation and 2.4 Å in the syn conformation. Molecular Dynamics. The B-form ssDNA dTpT structure was built using UCSF Chimera40 to edit a B-form dsDNA created using Nucleic Acid Builder (NAB).41 The dTpT was then solvated in a periodic box filled with TIP3P water using TLeap. The minimum distance from the box wall to all the DNA atoms was 1 nm. A single sodium counterion was added to neutralize the system. The total system size was 1050 atoms. We equilibrated and ran REMD with the Amber14 package with the force fields “ff14SB” and ff14SB/chi/ez (i.e., parmbsc0_chiOL4_ezOL1) and with the Gromacs package with the AMBER force field port “AMBER99SB-ILDN″, which uses the ff99 force field for DNA.15,42 During the equilibration, we applied a holonomic constraint via LINCS in Gromacs and via SHAKE in AMBER14 to constrain all bonds involving a hydrogen. The first equilibration step was the relaxation of the solvent around the DNA shell in 100 ps steepest-descent minimization in Gromacs and 100 steps steepest-descend in Amber14, followed by 20 ps NVT dynamics at each replica exchange temperature, T. The set of eight temperatures (300.00, 304.96, 310.00, 315.10, 320.27, 325.51, 330.83, and 336.22 K) were those that produced the desired exchange probability of 14% during the REMD. The initial temperatures were generated using the method of Patriksson.43 The REMD time step was set to 2 fs. The exchange attempted for REMD was every 1 ps in Amber and Gromacs. To match the 14% exchange probability with AMBER force field ff14SB, the gap between adjacent temperatures was adjusted by hand after 50 ns of production REMD. The same temperatures were used in Gromacs with Amber. To achieve NVT equilibration, we initially coupled the DNA and the solvent to baths at 10 K and T, respectively. Temperatures were maintained using a modified Berendsen thermostat with τ= 0.5 ps. During the first step of each NVT, we heated the DNA oligomers linearly from 10 to T. During the next 20 ps NVT, we coupled both DNA and solvent to the same thermostat at T. The third equilibration step was to perform 100 ps NPT dynamics at T and 1 bar, using the modified Berendsen barostat in Gromacs and the Monte Carlo barostat in Amber14. We initiated REMD production after NPT equilibration. The REMD production was performed in Amber and Gromacs under PME with a 10 Å cutoff. We monitored the ground state (300 K) MD trajectories for 500 ns. After 50 ns, the exchange probability reached 14%. Structures were saved every 2 ps and the statistical distributions were analyzed with MATLAB (Mathworks) and IGOR PRO (Wavemetrics).

2. METHODS Sample Preparation and NMR. Nonmodified thymidylyl(3′-5′)-thymidine (dTpT) oligomers were purchased from Sigma-Aldrich. A ∼200 nmol sample was dissolved without further purification in ∼0.7 mL of a 25-mM phosphate buffer solution at a pH of 7.0 to allow for the fixation of phosphate and then dried with argon. The residual sample was redissolved in deuterium oxide (99.9% D2O), dried, and redissolved in D2O to minimize the presence of water in the system with final dimer concentration of ∼0.3 mM. Spectra were collected using a Varian 600 MHz NMR spectrometer. We assigned the dTpT nonexchangeable protons from NOESY and DQF−COSY spectra acquired at 25 and 16 °C (see Figures S1 and S2 in the Supporting Information). NOESY spectra acquired at 16 °C were used for quantification of NOE interactions with mixing time (τm) set to 400 ms. Decoupling of 13P was applied in both dimensions. Sugar pucker was characterized using spin−spin coupling constants determined from peak splittings in these spectra.39 Internuclear distances (r) were determined from

3. RESULTS AND DISCUSSION NMR Results. The NMR spectra were used to determine three-bond scalar coupling constants (3JHH) between sugar protons by measuring the resulting peak splittings (Table 1).20,21 The 3JHH scalar-coupling is dependent on the dihedral angle defining the orientation of the three bonds and, C


Article

The Journal of Physical Chemistry B Table 1. Spin−Spin Constants, 3JHH, in dTpTa

Table 2. Predicted Sugar Pseudorotation Anglesa,b and base anti/syn distributionc from NMR with predictions from MDd in parentheses

constant, Hz

a

sugar couplings

DT1(5′) base

DT2(3′) base

H1′−H2′ H1′−H2″ H2′−H2″ H3′−H2′ H3′−H2″ H3′−H4′ H4′−H5′ H4′−H5″

7.1 6.3 14 6.5 3.6 3.6 3.5 4.3

6.8 6.7 14 6.3 4.3 4.3 2.4 3.2

PS νm,S XS P Nb νm, Nb X Nb % anti/sync

DT1 and DT2 denote the 5′- and 3′-end bases, respectively.

DT1 NMR (MD)

DT2 NMR (MD)

141° (146°) 26° (42°) 72.5% −19° 46° 27.5% NMR 86/14 (ff14SB 17/82) (ff14SB/chi/ez 38/62)

149° (133°) 33° (42°) 59.5% 18° 35° 40.5% NMR 90/10 (ff14SB 87/1) (ff14SB/chi/ez 96/4)

a

Calculated using the NMR 3JHH couplings with the Matlab GUI from Hendrickx, 2008.21 Parameters are for the “South” and “North” populations corresponding to the 2′-endo (“S” subscript) and 3′-endo (“N” subscript) structures, including: the pseudorotation angle, P; the puckering amplitude, νm; and mole fraction of each population, X. b Note that Altona has warned against the validity of the calculated parameters for the minor population, in this case, the North (A-form) population, which is poorly determined by the NMR couplings.19 Here, we report them for completeness though without further discussion. cNMR percentage of anti and syn conformations determined from H6:H1′ NOE cross-peak intensities. See text for details. The sum of MD percentages does not equal 100 in the ff14SB simulations because of the percentage in the high-anti region (χ = −35° to −90°). dBoth the ff14SB and ff14SB/chi/ez force fields produced the same pseudorotation distributions.

consequently, scalar couplings between sugar protons are useful for determining sugar pucker. The coupling constants result in our assignment of both sugars to a 2′-endo envelope conformation. The H1′−H2′ and H1′−H2″ splitting produces strong coupling constants of 7.1 and 6.3 Hz for DT1 and 6.8 and 6.7 Hz for DT2, respectively. This denotes a high C2′-endo population, in which the similarity of the values is likely due to motion of H2′ and H2″.44 The H3′-H2′ coupling was similarly strong with coupling constants 6.5 and 6.3 Hz; however, the coupling of H3′−H2″ was relatively weak with constants of 3.6 and 4.3 Hz. The coupling of H3′−H4′ is also weak, at 3.6 and 4.3 Hz, further confirming that the C3′ cannot be endo. Thus, the sugar conformation of dTpT is primarily 2′-endo, the conformation found in B-form DNA.45 A more detailed calculation of the sugar pucker can be calculated using the NMR coupling constants with freely available software from Hendrickx and Martins, following the method originally developed by De Leeuw and Altona.21,46 Using the Hendrickx program, the sugar 3JHH couplings predict the dominant population to be that of B-form DNA, with a 2′endo/1′-exo population with pseudorotation angle, P, of 141° for DT1 and 149° for DT2 (Table 2). However, the puckering amplitude, νm, is quite small for each ribose at 26° for DT1 and 33° for DT2 and the dominant structures only consist of 72.5% (DT1) and 59.5% (DT2) of the total population, indicating particularly flexible sugar conformations (Table 2). These values of P and νm place the sugars well within the normal distribution of B-form DNA structures.20 Plots of RMSD vs P and νm are shown in the Supporting Information, Figure S3, which also shows the minor contribution from A-form, 3′-endo structures. A detailed comparison of the NMR and MD predicted structures of the furanose in dTpT can be made based on the HCCH dihedral angles of each sugar. Table 3 presents the predicted HCCH dihedral angles calculated from P and νm, which were determined from the NMR 3JHH couplings. Dihedral angles were calculated using parameters from Wijk.19 We now turn to the orientation of the thymine about the glycosidic bond. The measurable proton distances of the dominant B-form conformation are depicted in Figure 3. The 2D NOESY spectrum containing cross-peaks from thymine’s H6 proton is shown in Figure 4 and the complete spectrum is in Figure S2. Both deoxyribose H1′ protons display a medium NOE cross peak to their own pyrimidine H6 protons (H1′:H6) advocating for an anti conformation about χ(Table 2). The strength of the H6:H1′ cross peaks predict an 86:14 anti/syn distribution for DT1 and a 90:10 anti/syn distribution for DT2. The H6:H2′ NOE cross peak on each base is strong, consistent with the anti conformation and the C2′ being out of plane in an

Table 3. Sugar HH Dihedral Angles (deg) from NMR and MD DT1 a

DT2 b

NMR

MDb

155 34 −28 93 −110

166 43 −29 94 −117

dihedral

NMR

MD

H1′C1′C2′H2′ H1′C1′C2′H2″ H2′C2′C3′H3′ H2″C2′C3′H3′ H3′C3′C4′H4′

148 27 −19 101 −116

167 44 −34 90 −105

a

a

NMR values were calculated using the pseudorotation values and the empirical parameters of Wijk.19 bMD values are determined by fitting the histograms of Figure 7 and S9. The histograms are identical for the ff14SB and the ff14SB/chi/ez force fields.

Figure 3. Proton couplings apparent via 1H−1H NOEs from H6.

endo conformation. On DT1, the deoxyribose H3′ is weakly coupled to its own pyrimidine H6, whereas on DT2, there is D


Article


Figure 4. NOESY spectrum with mixing time of 800 ms at 16 °C. The spectral region shown includes cross-peaks from aromatic protons (DT1 H6 at ω1 = 7.65 ppm and DT2 H6at ω1 = 7.67 ppm) to sugar protons. Note that for brevity, we label the 5′-thymine protons with the “T1” prefix and the 3′-thymine protons with “T2”.

Figure 5. Trajectories monitoring the glycosidic angle of DT1 (a and c) and DT2 (b and d) during REMD runs. Parts a and b were generated using ff14SB, and parts c and d were generated using ff14SB/chi/ez. Note: DT1 and DT2 denote the thymine at the 5′- and 3′-ends, respectively.

Figure 6. Conformational distribution of the glycosidic angles, χ, of DT1 and DT2 using (a) the ff14SB force field or (b) the latest ff14SB/chi/ez Amber force field. (c) Structure after initial equilibration, showing the anti−anti glycosidic conformation at the beginning of the simulation, and (d) after 500 ns REMD, showing the resultant syn-anti conformation.

B-form and indicates deviations of the phosphate backbone angles from those typical of B-form and a lack of π-stacking between the bases. In particular, the NMR spectra eliminate the possibility that either thymine takes on a syn conformation. If they did, one would expect strong H6:H1′ NOE coupling, regardless of sugar

weak to moderate H6:H3′ coupling. The weaker NOE cross peak we observe between H6 and H3′ on DT1 might support a twist deformation with C3′ pushed exo on the 5′-sugar, or it may indicate a rotation of the DT1 thymine to larger χ angles, moving H6 away from C3′. Lastly, the NOE between H6 of DT2 and H2″ of DT1 is substantially weaker than expected for E


Article


DT2 the REMD accurately predicts that the anti conformer will dominate the distribution. The syn or anti conformation of each thymine in the MD simulations is also apparent in the average distances from each H6, which are directly comparable to the NOE-predicted results. In particular, the histograms of H6:H1′ and H6:H3′ distances (Figure 7) can be used to calculate weighted average

pucker, and relatively weak H6:H2′ NOE coupling, none of which are observed. Molecular Dynamics Simulations. In this section, we present the results of the REMD simulation using both the ff14SB and ff14SB/chi/ez force fields. We concentrate on the sugar pucker and glycosidic angle using the H6 distance to all sugar protons to obtain information comparable to the NMR findings. The trajectories of the glycosidic angle during the 500 ns REMD runs are shown in Figure 5. During the 500 ns run with ff14SB, DT1 flips 40 times and DT2 flips 26 times, indicating relatively good sampling. Hence, with ff14SB, the modest barriers of 4−5 kcal/mol calculated by Krepl allow the production of good statistics for base flipping using REMD.32 However, with ff14SB/chi/ez, DT1 undergoes just 3 flips (Figure 5c) and DT2 only a single flip event (Figure 5d). The elimination of frequent base flipping events in ff14SB/chi/ez indicates quite clearly that, in dTpT, the barriers have increased significantly relative to ff14SB. This is consistent with the calculations of Krepl, who showed that ff14SB/chi/ez produces a large (7−8 kcal/mol) barrier to flip deoxyguanosine through χ = 0° and a modest 2−3 kcal/mol barrier around χ = 180°.32 The trajectories indicate that the ff14SB simulation has converged but that the statistics from the ff14SB/chi/ez simulation should be considered preliminary until longer simulations or multidimensional REMD can be employed.13 The base flipping events in ff14SB always occur by crossing gradually through χ = 0° and occur slowly in 20−40 ps (see Figure S4), consistent with a low and gradual potential energy barrier at that configuration as calculated by Krepl.32 In contrast, in ff14SB/chi/ez, the base flipping always occurs via rotation through the smaller barrier at χ = 180° with rapid relaxation to the anti or syn minima in a few ps (Figure S4). This indicates that the barrier at χ = 0° has increased enough to completely eliminate that pathway in ff14SB/chi/ez. The new force field, if accurate, would indicate that base flipping never occurs via passage of the thymine C2 carbonyl over the O4′ of the ribose and, instead, occurs only by passing the carbonyl over C2′. For DT2, flipping this direction requires passage of the C5 methyl group over O4′ where there is likely significant crowding from DT1, which may contribute to the infrequency of base flipping for DT2. However, when DT1 rotates through χ = 180°, the methyl group passes through open solvent without crowding from the adjacent thymine, producing the more facile flipping of DT1 in Figure 5c. Despite starting from the equilibrated B-form all-anti configuration, the REMD consistently predicts a tendency of the 5′-DT1 base to flip to a syn conformation. This is apparent in the histogram of χ angles (Figure 6, parts a and b). As shown in Table 2, when using ff14SB, DT1 is in the syn conformation 82% of the time and anti 17% of the time. (The remaining ∼1% of the distribution is due to “high-anti” conformations with χ angles between −35° and −90°.) However, DT2 is predicted to have a 87:1 anti:syn distribution in ff14SB, with an additional 12% in the high-anti conformer. This propensity of DT1 to flip to the syn configuration is reproduced in the ff14SB/chi/ez simulations, in which DT1 has an anti:syn ratio of 38:62, while DT2 has a ratio of 96:4. To clarify the nature of this base flip, Figure 6c shows the anti−anti conformation when the simulations are initiated and Figure 6d shows the flipped structure with the syn-anti conformation that is prevalent later in the run. Clearly, for DT1, the REMD predicts a prevalent syn conformer that is not consistent with the NMR results, but for

Figure 7. Comparison of population histograms determined using the old (ff14SB) and new (ff14SB/chi/ez) force fields in Amber. (a) H6:H1′ in ff14SB, (b) H6:H3′ in ff14SB, (c) H6:H1′ in ff14SB/chi/ ez, and (d) H6:H3′ in ff14SB/chi/ez. Also indicated are distances derived from NOESY NMR data (arrows), as well as REMD-predicted distances (vertical bars). The vertical axis, “Probability Distribution”, is in units of probability per Angstrom.

distances,⟨r−6⟩−1/6,for the REMD simulations that can be directly compared with the NOE calculated distances. The ff14SB REMD histogram of H6:H1′ distances (Figure 7a) predicts an average r−6 distance of 2.35 Å for DT1 and 3.57 Å for DT2, consistent with the syn(DT1), anti(DT2) conformation visible in the glycosidic angle distributions. The REMD distance for DT1 is barely changed in the ff14SB/chi/ez simulation (Figure 7c), but the DT2 distance is slightly smaller (3.38 Å) due to the slightly increased probability of the syn conformation for DT2 in ff14SB/chi/ez. With either force field, the REMD distance for DT1 is much too short compared to the predicted NMR H6:H1′ distance of 3.14 Å but the DT2 REMD distance is in rough agreement with that from NMR, 3.23 Å. The syn(DT1), anti(DT2) conformation also manifests itself in the average H6:H3′ distance (Figure 7b and d). With either force field, the DT1 H6:H3′ distance is quite large (4.83 Å with ff14SB (Figure 7b) and 4.59 Å with ff14SB/chi/ez (Figure 7d)) and much longer than the NMR distance of 3.44 Å. The DT2 H6:H3′ distance is significantly shorter (3.54 Å for ff14SB and 3.49 Å for ff14SB/chi/ez) and is much closer to agreement with the NMR distance of 2.96 Å. Similarly, the H6:H2′ distance is quite short for DT2 but significantly longer for DT1 (Figure S6). Lastly, the anti(DT2) MD conformation is evident in the close proximity predicted between the H6 of DT2 and the H2″ of DT1 (Figure S7). The REMD-predicted structure of dTpT clearly contradicts the NMR data that showed the anti conformation for both ends. F


Article


but have nothing directly comparable to NMR. Dihedral angles to H2″ are shown only in Figure S10 but agree with the shown distributions of H2′, shifted by ∼120°.) The sugar dihedral angle histograms were exactly the same in both force fields, as is expected given that they differ only in the χ and the ε/ζ potentials. Lastly, the most-probable dihedral angles from the MD simulations can be used to back calculate the P and νm from the MD simulations (Table 2), for comparison with those determined by NMR. For the deoxyribose, we find that the 5′-sugar conformation (DT1) predicted by MD is significantly different from that measured by NMR, but that the 3′-sugar (DT2) is fairly well reproduced from MD and NMR (see Figure 8 and Tables 1 and 2). For DT2, the most-probable pseudorotation angles, puckers and dihedral angles from MD agree very well with those calculated by NMR. The H1′C1′C2′H2′ dihedral (Figure 8a) histogram has a relatively narrow distribution for both sugars, while the histograms for the H2′C2′C3′H3′ (Figure 8b) and H3′C3′C4′H4′ (Figure 8c) are narrow distributions for the 5′-sugar (DT1, black) but have long tails for the 3′-sugar (DT2, red). Hence, the MD dihedral histograms indicate a significant fraction of the population of DT2 samples a distorted C3′ geometry, evident in the tail of the distribution in Figure 8, parts b and c. In contrast, MD predicts only a small variation in the structure of the DT1 sugar and the MD structure differs significantly from that of the NMR. It is surprising that the DT2 dihedral angles are so close, given that the pseudorotation angles from MD and NMR are significantly different. Most importantly, on each base, we see strong correlation between the agreement of the glycosidic angle and that of the sugar pucker. Given the strong coupling between the glycosidic angle and the sugar pucker, it is likely that the poor agreement in the DT1 sugar angles are induced by the syn flip of DT1 or vice versa. Similarly, the good agreement of DT2 sugar angles are likely made possible by, or induce, the correct anti orientation of DT2. The NMR data shows that the medium NOE interaction between H6 and H1′ is consistent with both thymine residues of dTpT being in an anti conformation. However, the MD structure obtained with both ff14SB and ff14SB/chi/ez predicts a syn χ angle for DT1 (Figure 6a). Concerning the possible origin of this glycosidic angle discrepancy, Zgarbova and coworkers proposed that it may be due to the short-ranged van der Waals and coulombic interactions, which take place during dihedral rotation.29 We have shown that the effect does not seem to be affected by the choice of water model but can be affected by changes in the seemingly peripheral backbone forces, which should allow force-field developers to focus on the intramolecular nucleic acid parameters. Prior predictions of thymine photodimerization have been obtained with both the current and the old Amber force field. From the work of Law7 it was proposed that MD can predict which conformations of the ground state lead to photodimerization when excited by UV light. To predict the conformation of dTpT in various cosolvents, which can be classified as dimerizable or nondimerizable, the Amber-99 force field was used in Gromacs.7 Prior to that work, Johnson used an MD simulation of 50 ns to explore the reactive conformers that exist in T18 ssDNA.9 We now know that the force field used was not able to reproduce accurately the distribution about χ, as was later shown by its reparameterization by Ode and Krepl.28,32 Other REMD simulations have shown that the current nucleic acid force fields have an excessive propensity for

To check whether or not the water model may play a role through the coupling between dTpT and TIP3P water, we repeated the REMD with the SPC/E water model (Figure S8b) and found that the glycosidic angle distribution exactly matches that produced with TIP3P, indicating that the parametrization of χ is likely to be independent of the water model. Additionally, to test the dependence of the glycosidic angle distribution on changes to peripheral backbone forces, we repeated the REMD simulations with the ff99 potential, which lacks the updates to the α (O3′−P−O5′−C5′ dihedral) and γ (C5′−C4′−C3′−O3′ dihedral) potentials in the DNA backbone that were implemented with the bsc0 updates.26 Curiously, the ff99 calculations (Figure S8c) do predict a larger percentage of the DT1 anti conformer (30%) compared to ff14SB (17%) despite the fact that the DT1 glycosidic bond is so far from where the changes in the backbone potential occurred. This suggests that the errors we are observing in χ may be due to subtle effects in the overall structure of the sugar and backbone of dTpT and not necessarily due to errors in the χ potential. A comparison of the NMR- and MD-predicted structures of the deoxyribose in dTpT can be made based on the HCCH dihedral angles of each sugar. Table 3 presents the average HH dihedral angles calculated from P and νm, which were determined from the NMR JHH couplings using parameters from Wijk.19 These can be compared with the distributions of dihedral angles observed in the MD simulations, shown in Figure 8. (The complete set of sugar torsions are in Figure S9,

Figure 8. HCCH dihedral angles of the deoxyribose sugar determined by MD and NMR: (a) H1′C1′C2′H2′, (b) H2′C2′C3′H3′, and (c) H3′C3′C4′H4′. Both force fields, ff14SB and ff14SB/chi/ez, gave exactly the same distributions. NMR dihedral angles (vertical lines) were calculated from the pseudorotation angles determined by JHH couplings. G


Article

The Journal of Physical Chemistry B base-stacking in single stranded nucleic acids.13 Given the significant problems in using MD to predict ssDNA χ angles, sugar pseudorotations, and base stacking, it brings into question how well previous dTpT studies measured reasonable statistics for the C5C6:C5C6 distances. Also, we have shown that dTpT does not converge in 100 ns of direct MD and therefore neither of the previous simulations sampled a converged distribution. Ironically, the fact that none of the previous studies sampled a converged distribution may have improved their accuracy; the sampled structures would have been biased toward structures similar to the initially built B-form structure, which is closer to the actual solution structure than the equilibrated, converged structure would have been. To proceed with confidence, the conformational studies of dTpT7,9 and dTpT flanked by purines11,12 should be subsequently revisited. In his work with cytidine and uridine for glycosidic reparameterization, Yildirim30 suggested that one cannot assess the glycosidic (χ) behavior in a force field with dsDNA because it is restricted by the hydrogen bonding network. Now that REMD has proven to be a method to obtain structural convergence in ssDNA systems, future optimizations should focus on parametrizing the force field so that it produces structures consistent with the few ssDNA experimental structures that are in the literature.

Integrated Research Computing at the University of Rochester, to the Grossfield, Turner, Mathews, and Roitberg laboratories for assistance with the molecular dynamics simulations and to Danielle Barnett for assistance with the samples and manuscript.

■

(1) Beukers, R.; Berends, W. Isolation and Identification of the Irradiation Product of Thymine. Biochim. Biophys. Acta 1960, 41, 550− 551. (2) Bernerd, F.; Sarasin, A.; Magnaldo, T. Galectin-7 Overexpression Is Associated with the Apoptotic Process in UVB-Induced Sunburn Keratinocytes. Proc. Natl. Acad. Sci. U. S. A. 1999, 96, 11329−11334. (3) Wang, Y.; DiGiovanna, J. J.; Stern, J. B.; Hornyak, T. J.; Raffeld, M.; Khan, S. G.; Oh, K.-S.; Hollander, M. C.; Dennis, P. A.; Kraemer, K. H. Evidence of Ultraviolet Type Mutations in Xeroderma Pigmentosum Melanomas. Proc. Natl. Acad. Sci. U. S. A. 2009, 106, 6279−6284. (4) Pfeifer, G. P. Formation and Processing of Uv Photoproducts: Effects of DNA Sequence and Chromatin Environment. Photochem. Photobiol. 1997, 65, 270−283. (5) Schreier, W. J.; Schrader, T. E.; Koller, F. O.; Gilch, P.; CrespoHernández, C. E.; Swaminathan, V. N.; Carell, T.; Zinth, W.; Kohler, B. Thymine Dimerization in DNA is an Ultrafast Photoreaction. Science 2007, 315, 625−629. (6) Crespo-Hernandez, C. E.; Cohen, B.; Kohler, B. Base Stacking Controls Excited-State Dynamics in A•T DNA. Nature 2005, 436, 1141−1144. (7) Law, Y. K.; Azadi, J.; Crespo-Hernández, C. E.; Olmon, E.; Kohler, B. Predicting Thymine Dimerization Yields from Molecular Dynamics Simulations. Biophys. J. 2008, 94, 3590−3600. (8) Johns, H. E.; Pearson, M. L.; LeBlanc, J. C.; Helleiner, C. W. The Ultraviolet Photochemistry of Thymidylyl-(3′→5′)-Thymidine. J. Mol. Biol. 1964, 9, 503−IN501. (9) Johnson, A. T.; Wiest, O. Structure and Dynamics of Poly(T) Single-Strand DNA: Implications toward CPD Formation. J. Phys. Chem. B 2007, 111, 14398−14404. (10) McCullagh, M.; Hariharan, M.; Lewis, F. D.; Markovitsi, D.; Douki, T.; Schatz, G. C. Conformational Control of TT Dimerization in DNA Conjugates. A Molecular Dynamics Study. J. Phys. Chem. B 2010, 114, 5215−5221. (11) Pan, Z.; McCullagh, M.; Schatz, G. C.; Lewis, F. D. Conformational Control of Thymine Photodimerization in PurineContaining Trinucleotides. J. Phys. Chem. Lett. 2011, 2, 1432−1438. (12) Pan, Z.; Hariharan, M.; Arkin, J. D.; Jalilov, A. S.; McCullagh, M.; Schatz, G. C.; Lewis, F. D. Electron Donor−Acceptor Interactions with Flanking Purines Influence the Efficiency of Thymine Photodimerization. J. Am. Chem. Soc. 2011, 133, 20793−20798. (13) Bergonzo, C.; Henriksen, N. M.; Roe, D. R.; Swails, J. M.; Roitberg, A. E.; Cheatham, T. E. Multidimensional Replica Exchange Molecular Dynamics Yields a Converged Ensemble of an RNA Tetranucleotide. J. Chem. Theory Comput. 2014, 10, 492−499. (14) Henriksen, N. M.; Roe, D. R.; Cheatham, T. E. Reliable Oligonucleotide Conformational Ensemble Generation in Explicit Solvent for Force Field Assessment Using Reservoir Replica Exchange Molecular Dynamics Simulations. J. Phys. Chem. B 2013, 117, 4014− 4027. (15) Sorin, E. J.; Rhee, Y. M.; Pande, V. S. Does Water Play a Structural Role in the Folding of Small Nucleic Acids? Biophys. J. 2005, 88, 2516−2524. (16) Sen, S.; Andreatta, D.; Ponomarev, S. Y.; Beveridge, D. L.; Berg, M. A. Dynamics of Water and Ions near DNA: Comparison of Simulation to Time-Resolved Stokes-Shift Experiments. J. Am. Chem. Soc. 2009, 131, 1724−1735. (17) Isaksson, J.; Acharya, S.; Barman, J.; Cheruku, P.; Chattopadhyaya, J. Single-Stranded Adenine-Rich DNA and RNA Retain Structural Characteristics of Their Respective Double-Stranded

■

CONCLUSION We have shown that the current parameters of the AMBER force field 14 (the ff14SB and ff14 + e/zOL1 + qOL4) fail to accurately reproduce the conformation of dTpT due to the stabilization of the 5′-thymine in a syn conformation. Smaller discrepancies are found in the sugar pucker of the 5′deoxyribose. This work shows the importance of benchmarking new nucleic acid force fields with single-stranded DNA oligomers, which exhibit structures that are extremely sensitive to the force-field parameters. Clearly, parameters that faithfully reproduce structural characteristics of double stranded DNA have difficulties when simulating these small ssDNA systems.

■

ASSOCIATED CONTENT

S Supporting Information *

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jpcb.6b00191. Sums of NMR coupling constants, dTpT DQF−COSY and dTpT NOESY spectra, pseudorotation wheel, examples of base flipping events, schematic syn nucleotide, statistics of the H6:H2′ and H6(DTn):H2″(DTm) distances from REMD, comparison of glycosidic and delta angles, conformation of the five pseudorotation torsions, dihedral angles to H2″ (PDF)

■

REFERENCES

AUTHOR INFORMATION

Corresponding Author

*(D.W.M.) E-mail: [email protected]. Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS This project was supported by grants from the National Center for Research Resources (5R21RR025344-03) and the National Institute of General Medical Sciences (8R21GM103438-03) from the National Institutes of Health. D.W.M. was supported as an Alfred P. Sloan Fellow. We are grateful to the Center for H


Article

The Journal of Physical Chemistry B Conformations and Show Directional Differences in Stacking Pattern. Biochemistry 2004, 43, 15996−16010. (18) Altona, C.; Sundaralingam, M. Conformational Analysis of the Sugar Ring in Nucleosides and Nucleotides. New Description Using the Concept of Pseudorotation. J. Am. Chem. Soc. 1972, 94, 8205− 8212. (19) Van Wijk, J.; Huckriede, B. D.; Ippel, J. H.; Altona, C. Furanose Sugar Conformations in DNA from Nmr Coupling Constants. Methods Enzymol. 1992, 211, 286−306. (20) Sun, G.; Voigt, J. H.; Filippov, I. V.; Marquez, V. E.; Nicklaus, M. C. Prosit: Pseudo-Rotational Online Service and Interactive Tool, Applied to a Conformational Survey of Nucleosides and Nucleotides. J. Chem. Inf. Model. 2004, 44, 1752−1762. (21) Hendrickx, P.; Martins, J. A User-Friendly Matlab Program and GUI for the Pseudorotation Analysis of Saturated Five-Membered Ring Systems Based on Scalar Coupling Constants. Chem. Cent. J. 2008, 2, 20. (22) Lu, X.-J.; Olson, W. K. 3DNA: A Versatile, Integrated Software System for the Analysis, Rebuilding and Visualization of ThreeDimensional Nucleic-Acid Structures. Nat. Protoc. 2008, 3, 1213− 1227. (23) IUPAC. Abbreviations and Symbols for the Description of Conformations of Polynucleotide Chains. Eur. J. Biochem. 1983, 131, 9−15. (24) Amber 2015 Reference Manual, 2015. http://ambermd.org/ doc12/Amber15.pdf. (25) Cheatham, T. E., III; Kollman, P. A. Molecular Dynamics Simulations of Nucleic Acids. Annu. Rev. Phys. Chem. 2000, 51, 435− 471. (26) Pérez, A.; Marchán, I.; Svozil, D.; Sponer, J.; Cheatham, T. E., III; Laughton, C. A.; Orozco, M. Refinement of the Amber Force Field for Nucleic Acids: Improving the Description of α/γ Conformers. Biophys. J. 2007, 92, 3817−3829. (27) Svozil, D.; Šponer, J. E.; Marchan, I.; Pérez, A.; Cheatham, T. E.; Forti, F.; Luque, F. J.; Orozco, M.; Šponer, J. Geometrical and Electronic Structure Variability of the Sugar−Phosphate Backbone in Nucleic Acids. J. Phys. Chem. B 2008, 112, 8188−8197. (28) Ode, H.; Matsuo, Y.; Neya, S.; Hoshino, T. Force Field Parameters for Rotation around X Torsion Axis in Nucleic Acids. J. Comput. Chem. 2008, 29, 2531−2542. (29) Zgarbová, M.; Otyepka, M.; Sponer, J. i.; Mládek, A. t.; Banás, P.; Cheatham, T. E.; Jurecka, P. Refinement of the Cornell et al. Nucleic Acids Force Field Based on Reference Quantum Chemical Calculations of Glycosidic Torsion Profiles. J. Chem. Theory Comput. 2011, 7, 2886−2902. (30) Yildirim, I.; Stern, H. A.; Kennedy, S. D.; Tubbs, J. D.; Turner, D. H. Reparameterization of RNA X Torsion Parameters for the Amber Force Field and Comparison to NMR Spectra for Cytidine and Uridine. J. Chem. Theory Comput. 2010, 6, 1520−1531. (31) Yildirim, I.; Stern, H. A.; Tubbs, J. D.; Kennedy, S. D.; Turner, D. H. Benchmarking Amber Force Fields for RNA: Comparisons to NMR Spectra for Single-Stranded r(GACC) Are Improved by Revised X Torsions. J. Phys. Chem. B 2011, 115, 9261−9270. (32) Krepl, M.; Zgarbová, M.; Stadlbauer, P.; Otyepka, M.; Banás,̌ P.; Koča, J.; Cheatham, T. E.; Jurečka, P.; Šponer, J. Reference Simulations of Noncanonical Nucleic Acids with Different X Variants of the Amber Force Field: Quadruplex DNA, Quadruplex RNA, and ZDNA. J. Chem. Theory Comput. 2012, 8, 2506−2520. (33) Zgarbová, M.; Luque, F. J.; Šponer, J.; Cheatham, T. E.; Otyepka, M.; Jurečka, P. Toward Improved Description of DNA Backbone: Revisiting Epsilon and Zeta Torsion Force Field Parameters. J. Chem. Theory Comput. 2013, 9, 2339−2354. (34) Neelakandan, P. P.; Pan, Z.; Hariharan, M.; Lewis, F. D. Facially-Selective Thymine-Thymine Photodimerization in TTT Triads. Photochem. Photobiol. Sci. 2012, 11, 889−892. (35) Pan, Z.; Chen, J.; Schreier, W. J.; Kohler, B.; Lewis, F. D. Thymine Dimer Photoreversal in Purine-Containing Trinucleotides. J. Phys. Chem. B 2012, 116, 698−704.

(36) Pan, Z.; McCullagh, M.; Schatz, G. C.; Lewis, F. D. Conformational Control of Thymine Photodimerization in PurineContaining Trinucleotides. J. Phys. Chem. Lett. 2011, 2, 1432−1438. (37) Roe, D. R.; Bergonzo, C.; Cheatham, T. E. Evaluation of Enhanced Sampling Provided by Accelerated Molecular Dynamics with Hamiltonian Replica Exchange Methods. J. Phys. Chem. B 2014, 118, 3543−3552. (38) Romo, Tod D.; Grossfield, A. Unknown Unknowns: The Challenge of Systematic and Statistical Error in Molecular Dynamics Simulations. Biophys. J. 2014, 106, 1553−1554. (39) Kruger, K.; Grabowski, P. J.; Zaug, A. J.; Sands, J.; Gottschling, D. E.; Cech, T. R. Self-Splicing RNA: Autoexcision and Autocyclization of the Ribosomal RNA Intervening Sequence of Tetrahymena. Cell 1982, 31, 147−157. (40) Pettersen, E. F.; Goddard, T. D.; Huang, C. C.; Couch, G. S.; Greenblatt, D. M.; Meng, E. C.; Ferrin, T. E. UCSF Chimeraa Visualization System for Exploratory Research and Analysis. J. Comput. Chem. 2004, 25, 1605−1612. (41) Macke, T. J.; Case, D. A. In Molecular Modeling of Nucleic Acids; American Chemical Society: 1997; Vol. 682; pp 379−393. (42) DePaul, A. J.; Thompson, E. J.; Patel, S. S.; Haldeman, K.; Sorin, E. J. Equilibrium Conformational Dynamics in an RNA Tetraloop from Massively Parallel Molecular Dynamics. Nucleic Acids Res. 2010, 38, 4856−4867. (43) Patriksson, A.; van der Spoel, D. A Temperature Predictor for Parallel Tempering Simulations. Phys. Chem. Chem. Phys. 2008, 10, 2073−2077. (44) Li, H.; Huang, F.; Shaw, B. R. Conformational Studies of Dithymidine Boranomonophosphate Diastereoisomers. Bioorg. Med. Chem. 1997, 5, 787−795. (45) Znosko, B. M.; Barnes, T. W.; Krugh, T. R.; Turner, D. H. NMR Studies of DNA Single Strands and DNA:RNA Hybrids with and without 1-Propynylation at C5 of Oligopyrimidines. J. Am. Chem. Soc. 2003, 125, 6090−6097. (46) De Leeuw, F. A. A. M.; Altona, C. Computer-Assisted Pseudorotation Analysis of Five-Membered Rings by Means of Proton Spin−Spin Coupling Constants: Program Pseurot. J. Comput. Chem. 1983, 4, 428−437.

I


Disagreement Between the Structure of the dTpT Thymine Pair

Recommend Documents