DNA Polymerase Î» Active Site Favors a Mutagenic Mispair Between

Subscriber access provided by UNIV OF NEWCASTLE

Article

DNA Polymerase # Active Site Favors a Mutagenic Mispair Between the Enol Form of Deoxyguanosine Triphosphate Substrate and the Keto Form of Thymidine Template. A Free Energy Perturbation Study Sergey N. Maximoff, Shina C. L. Kamerlin, and Jan Florian J. Phys. Chem. B, Just Accepted Manuscript • DOI: 10.1021/acs.jpcb.7b04874 • Publication Date (Web): 21 Jul 2017 Downloaded from http://pubs.acs.org on July 23, 2017

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 free 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 accessible to all readers and 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.

The Journal of Physical Chemistry B 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 40

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

The Journal of Physical Chemistry

DNA Polymerase λ Active Site Favors a Mutagenic Mispair between the Enol Form of Deoxyguanosine Triphosphate Substrate and the Keto Form of Thymidine Template. A Free Energy Perturbation Study. Sergey N. Maximoff,∗,† Shina Caroline Lynn Kamerlin,‡ and Jan Florián∗,† Department of Chemistry and Biochemistry, Loyola University, Chicago, Illinois, United States, and Science for Life Laboratory, Department of Cell and Molecular Biology, Uppsala University, Uppsala, Sweden E-mail: [email protected]; [email protected]

∗

To whom correspondence should be addressed Department of Chemistry and Biochemistry, Loyola University, Chicago, Illinois, United States ‡ Science for Life Laboratory, Department of Cell and Molecular Biology, Uppsala University, Uppsala, Sweden †

1

ACS Paragon Plus Environment


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

Abstract Human DNA polymerase λ is an intermediate fidelity member of the X family, which plays a role in DNA repair. Recent X-ray diffraction structures of a ternary complex of a loop-deletion mutant of polymerase λ, a deoxyguanosine triphosphate analog and a gapped DNA show that guanine and thymine form a mutagenic mispair with an unexpected Watson-Crick-like geometry rather than a wobble geometry. Hence, there is an intriguing possibility that either thymine in the DNA or guanine in the deoxyguanosine triphosphate analog may spend a substantial fraction of time in a deprotonated or enol form (both are minor species in aqueous solution) in the active site of the polymerase λ mutant. The experiments do not determine particular forms of the nucleobases that contribute to this mutagenic mispair. Thus, we investigate the thermodynamics of formation of various mispairs between guanine and thymine in the ternary complex at a neutral pH using classical molecular dynamics simulations and the free energy perturbation method. Our free energy calculations, as well as a comparison of the experimental and computed structures of mispairs, indicate that the Watson-Crick-like mispair between the enol tautomer of guanine and keto tautomer of thymine is dominant. The wobble mispair between the keto forms of guanine and thymine and the Watson-Crick-like mispair between the keto tautomer of guanine and enol tautomer of thymine are less prevalent, and mispairs that involve deprotonated guanine or thymine are thermodynamically unlikely. These findings are consistent with the experiment and relevant for understanding mechanisms of spontaneous mutagenesis.

Introduction Encoding of genetic information is robust because of the complementarity of purine and pyrimidine nucleobases in a double-stranded DNA. 1,2 Hydrogen bonding between nucleobases, their shapes and their interactions with their environment together interlace in a complex way to define the concept of nucleobase complementarity. 3–6 Adenine pairs with thymine and guanine pairs with cytosine in a Watson-Crick (W-C) geometry. On the other 2


Page 2 of 40

Page 3 of 40

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


hand, guanine mispairs with thymine (see fig. 1) and adenine mispairs with cytosine in a wobble geometry. 3,7–9 This picture applies to the nucleobases in their forms most prevalent thermodynamically, which are the keto tautomers of guanine (G) and thymine (T ) and the amino tautomers of adenine (A) and cytosine (C). However, the keto tautomers of guanine and thymine coexist with their enol tautomers (G∗ and T ∗ ) and deprotonated forms (G− and T − ) just as amino tautomers of adenine and cytosine coexist with their imino tautomers and protonated forms. These forms are remarkable since they may mispair in W-C geometries, which are sterically similar to those of canonical nucleobase pairs G · C and A · T . In particular, guanine and thymine may mispair in a W-C geometry when either guanine and thymine becomes an enol tautomer (G∗ · T and G · T ∗ ) or a deprotonated form (G− · T and G · T − ) in fig. 1. Alternatively, the imino tautomer of adenine and the amino tautomer of cytosine or the amino tautomer of adenine and the imino tautomer of cytosine may pair in a W-C geometry. 10 These mispairs with W-C geometry may further undergo transition mutations into canonical W-C pairs A · T or G · C under the action of a DNA polymerase. 3,11 Thus, it has been acknowledged early on that study of equilibria between tautomers 1,2,11 as well as protonated and deprotonated 12,13 nucleobases in polymerases is central to our understanding of spontaneous transition mutations. Tautomerism for individual nucleobases or nucleobase pairs in isolation has been a subject of experimental 14,15 and theoretical studies. 15–18 Structure and stability of nucleobase pairs and their tautomer forms in a DNA 19 and in isolation 20 have been studied using dispersion-corrected density functional theory and an implicit solvent model. Protonationdeprotonation equilibria in solution for nucleobases have also been investigated experimentally 21–24 and computationally. 25–27 Mispairs with W-C geometry due to ionized modified nucleobases appeared in an NMR study of a DNA in aqueous solution. 28 Transient mispairs in DNA in aqueous solution have also been resolved in a recent NMR study. 29 The environment in an enzyme active site, however, differs drastically from an aqueous solution or gas phase. 30–33 It is also harder to access experimentally. 34

3



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

Recently, a mispair of guanine and thymine with a W-C geometry has been reported by Bebenek et al. 35 in the active site of a loop deletion mutant of human DNA polymerase λ. More reports of experimental detection of transient W-C mispairs involving minor species within complexes with DNA polymerases 32,33,36,37 or within RNA or ribosomes 38,39 have followed. Kinetic study of modified nucleobase insertion by a polymerase suggest the formation of ionized nucleobase pairs. 40 Apart from an investigation of guanine and uracil mispair with W-C geometry in a ribosome, 41 we are not aware of theoretical work that addresses nucleobase tautomerism in an enzyme. In particular, the origin of W-C mispair in the active site of the loop deletion mutant of human polymerase λ 35 remains undetermined. Human polymerase λ (pol λ) is an intermediate fidelity DNA polymerase (error rate of around 10−4 42 ) that belongs to the X family of DNA polymerases and participates in both the base excision repair as well as in repairing double-stranded breaks in DNA. 42,43 The entire X family of enzymes is believed to have evolved from a common pol λ-like progenitor. 42,44 The mutant loop deletion type, pol λ/DL, arises from the wild-type pol λ by replacing nine amino acid residues in the loop 1 with four amino acid residues. The structure of the active site in the mutant and wild-type pol λ are nearly identical since the deleted loop resides far away from the active site, and the deletion does not alter the protein conformation. 45 Although this mutation does not affect the active site geometry locally, it does reduce the fidelity of the enzyme significantly. 45,46 Pol λ/DL forms a ternary complex with a DNA with a one-nucleotide gap in the primer strand and a deoxyribonucleoside triphosphate (dNTP) that the enzyme incorporates into the gap. Bebenek et al. made three observations. 35 First, their X-ray diffraction experiments reveal a mispair with a W-C geometry between guanine in a deoxyguanosine triphosphate (dGTP) analog and thymine in the template strand of the DNA in the active site of pol λ/DL for a structure obtained at a neutral pH (see fig. 2). Mispairs G∗ · T , G · T ∗ , G− · T , and G · T − are possible contributors to the W-C geometry. Second, the distance of 2.66 Å between the O6 atom in the guanine and the O4 atom in the thymine is consistent with a

4


Page 4 of 40

Page 5 of 40

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


hydrogen bond adjacent to these atoms. W-C mispairs G∗ · T or G · T ∗ but not G− · T or G · T − involve a hydrogen bond between the O6 and O4 atoms (see fig. 1). Third, kinetic experiments indicate that the rate of misinsertion of dGTP next to the template T increases with the pH, which is consistent with the increase in the incorporation of mispairs G− · T or G · T − at higher pH. 35 However, pH-dependent changes in the enzyme might also cause this trend. Thus, it remains unclear whether G∗ · T or G · T ∗ or G− · T or G · T − or a statistical mixture of the above appears in the experiment. In what follows, we use molecular dynamics simulations (MD) with a classical force field as well as the free energy perturbation method (FEP) to determine whether keto, enol, or deprotonated forms of guanine and thymine pair in the active site of pol λ/DL.

Methods Simulation strategy Differences ∆Gp (N1 · N1′ → N2 · N2′ ) in free energies for the mutation of a nucleobase pair N1 · N1′ into another nucleobase pair N2 · N2′ quantify the relative thermodynamic stability of N1 · N1′ vs. N2 · N2′ in the ternary complex p of the pol λ/DL with the DNA and dNTP. Here, N1 and N2 are G or G∗ or G− in the dNTP substrate, and N1′ and N2′ are T or T ∗ or T − in the nucleotide within the template strand of the gapped DNA opposite to the dNTP substrate (see fig. 2). ∆Gp (N1 · N1′ → N2 · N2′ ) includes contributions due to both intramolecular (e.g., due to changes in bonds and bond angles) and intermolecular interactions for the active site. A modified AMBER f14SB 47 molecular mechanics force field that we use here describes well the latter but not the former. Therefore, we apply a thermodynamic cycle in fig. 3(a), ∆Gp (N1 · N1′ → N2 · N2′ ) = ∆Gw (N1 · N1′ → N2 · N2′ ) + ∆∆Gw→p (N1 · N1′ → N2 · N2′ ). (1)

5



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 6 of 40

Here, ∆Gw (N1 · N1′ → N2 · N2′ ) = ∆Gw (N1 → N2 ) + ∆Gw (N1′ → N2′ ) is the free energy change for the mutation of a pair of individual nucleosides in aqueous solution in the limit of infinite dilution. This term includes the intramolecular contributions, and it should be evaluated using quantum chemistry calculations or inferred from the experimental data. On the other hand, the contribution

∆∆Gw→p (N1 · N1′ → N2 · N2′ ) = ∆Gw→p (N2 · N2′ ) − ∆Gw→p (N1 · N1′ ) = ∆GFp EP (N1 · N1′ → N2 · N2′ ) − ∆GFw EP (N1 · N1′ → N2 · N2′ )

(2)

is dominated by the intermolecular contributions since intramolecular contributions mostly cancel out. Hence, ∆∆Gw→p (N1 · N1′ → N2 · N2′ ) can be estimated accurately as the difference between ∆GFp EP and ∆GFw EP , which are computed with FEP MD using the molecular mechanics force field.

Estimates for ∆∆Gw→p (N1 · N1′ → N2 · N2′ ) FEP calculations We evaluate a free energy change using FEP along a path (parametrized by a coupling constant λ), which deforms the initial state to the final state. The free energy change 48 is

∆G = −RT

X i

Ui+1 − Ui ln exp − RT i

(3)

for the path subdivided into discrete windows. Here, Ui denotes the potential energy surface at a discreet value λi of the coupling constant, R is the gas constant, and T is the thermodynamic temperature. The modified AMBER ff14SB + parmbsc0_chiOL4_ezOL1 force field approximates potential energy surfaces. 47 Angular brackets h· · · ii indicate the average over an ensemble of configurations that a molecular dynamics simulation at the potential energy 6


Page 7 of 40

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


surface Ui generates. To parametrize the potential energy surfaces Ui along the path, we rely on the singletopology framework, 49 which describes the changes in the potential energy along the path in terms of the union of atoms in the initial and final states. The initial, intermediate and final potential energy surfaces then arise by gradually changing the force field (i.e., atomic charges and van der Waals (vdW) parameters as well as bond, bond angle, torsion angle, and improper torsion angle terms), which sequences of atoms in the single-topology framework label. Destruction or creation of atoms and their vdW potentials leads to singularities, which impact the numerical accuracy of FEP MD simulations. 50 Therefore, we employ non-singular soft-core vdW potentials to replace the infinite value of the hard-core vdW potential with a finite value of 20 kcal · mol−1 for those atoms that are created or destroyed. 51,52 A Q-region in fig. 2 encompasses all the atoms that contribute to changes in the potential energy along the path in this report. Table 1 describes the path that connects the initial and final states. This path contains five blocks, which are further subdivided in 17, 21, 321, 21, and 17 windows, respectively. Along the first block of 17 windows, the soft-core dummy atom types to be created are changed to the target atom types while soft-core potentials are still on and vdW parameters for all atoms that are not created or destroyed are adjusted. Along the second block of 21 windows, bond terms are changed, and the soft-core vdW potentials involving atoms that are created are transformed into the hardcore potentials. Along the third block of 321 windows, atomic charges in the Q-region are evolved, charges of the atoms to be destroyed are annihilated, and charges of the atoms to be created are turned on. Along the fourth block of 21 windows, bond angle terms are adjusted while the hard-core vdW potentials involving atoms that are destroyed are transformed into the soft-core potentials. Along the fifth block of 17 windows, torsion and improper torsion terms are adjusted while vdW terms involving atoms to be destroyed are further annihilated into dummy soft core vdW terms.

7



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 8 of 40

Force field parametrization The configuration ensembles in our free energy calculations are obtained from MD trajectories over the potential energy surfaces Ui that the modified AMBER ff14SB + parmbsc0_chiOL4_ezOL1 force field represents. The latter includes a series of enhancements 47,53 of the earlier AMBER force fields 54,55 for proteins and nucleic acids. Furthermore, α/γ, 56 ǫ/ζ, 57 and χ 58 torsional terms refinements for DNA are added on top of AMBER ff14SB. The force field parameters are ported from the AMBER format 47 to a format suitable for our simulations in the Q simulation package, version 5.10.7. 52 Parameters for deprotonated and enol forms of guanine and thymine nucleosides, DNA fragments, and nucleoside triphosphates are not parts of the standard AMBER ff14SB + parmbsc0_chiOL4_ezOL1 force field.

Generalized AMBER force field (GAFF) is a

parametrization, which applies to a broader range of organic compounds. 59 In our parametrization for the nucleosides, DNA fragments, and nucleoside triphosphates, we combine the AMBER and GAFF atom types as follows. Types of atoms within nucleobases other than those immediately adjacent to N1 in thymine or N9 in guanine (see fig. S1) are GAFF types (shown in the lower case). Meanwhile, types of the remaining atoms are AMBER types (shown in the upper case) as tables S1 to S6 list. Those bonding, angular, torsion and improper torsion parameters that are indexed by both GAFF and AMBER types are set to the respective AMBER values whenever possible or otherwise to the respective GAFF values according to a translation of AMBER atom types to GAFF atom types. We note that both AMBER and GAFF force fields are successful in describing solvation energies and water/chloroform partition coefficients of nucleobases. 60 The atom-centered RESP charges are fitted using a scheme, which is described for nucleosides and DNA fragments in Cieplak et al. 61 and for nucleoside triphosphates in Meagher et al. 62 Specifically, charge distributions in an isolated nucleoside and dimethyl phosphate or methyl triphosphate anions determine AMBER charges in DNA fragments and nucleoside triphosphates through gluing conditions that identify charges at 5′ and 3′ hydroxyls within 8


Page 9 of 40

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 nucleoside and those of a methyl group in the dimethyl phosphate or methyl triphosphate. Starting from HF/6-31G* geometries of G, G∗ , G− , T , T ∗ , T − nucleosides, we compute the electrostatic potential at HF/6-31G* level over a Mertz-Kollman grid (four spherical layers per atom) as implemented in the Gaussian suite of computer programs. 63 Then the potential is fitted to a potential of atom-centered point charges using RESP protocol 61 as implemented in AmberTools. 47 In all calculations, C2′ -endo conformations of the nucleosides are used. This particular level of model chemistry ensures the compatibility with the AMBER parametrization. 61 Tables S1 to S6 list atomic charges for various forms of guanine and thymine nucleosides, DNA fragments, and nucleoside triphosphates. We use the parameters of refs. 64,65 for Mg2+ . TIP3P parametrization 52,66 describes water molecules explicitly in our simulations. Boundary conditions and cut-offs To model systems in infinite aqueous solution, we employ the surface-constraint all-atom solvent (SCAAS) model 67,68 as implemented in the Q suite of computer programs. 52 Three concentric regions partition the system as fig. 2 shows. A sphere of 21 Å radius encloses the interior region. The surface region is a shell between the concentric spheres of 21 and 24.5 Å radii, respectively. The outer region encloses the remaining space. The solvent in the outer region is modeled implicitly, and the Born formula gives its effect on the inner region. −2

Harmonic restraints with a force constant of 200 kcal·mol−1 ·Å

localize solute atoms in the

outer region to their initial coordinates, and non-bonding interactions involving these atoms are turned off. The interior and surface regions contain solute and explicit solvent molecules. The harmonic polarization constraint with a force constant of 20 kcal · mol−1 · rad−2 ensures that orientation of solvent molecules in the surface region is consistent with being a part of the −2

bulk solution. Meanwhile, the radial constraint with a force constant of 60 kcal·mol−1 ·Å

in

the surface region controls the solvent density and prevents evaporation into the outer region. Nonbonding interactions between other atoms inside the simulation sphere are subjected to a 9



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 10 of 40

10 Å cut-off. Beyond the cut-off distance, the local reaction field (LRF) method 68,69 describes long-range electrostatic interactions. Nonbonding interactions of atoms in the Q-region are explicitly evaluated for all distances. Preparation of the initial structure of the simulated system Our FEP MD simulations of pol λ/DL ternary complexes rely on the experimental X-ray diffraction crystal structure (PDB code 3PML). 35 The crystal structure contains a dimer of pol λ/DL with monomers linked by a disulfide bridge. Since the active form of pol λ/DL is monomeric, we select a monomer unit in the unit cell that includes the W-C mispair between the guanine and the thymine in the ternary complex. The PDB structure from the databank is pre-processed using the amber4pdb and leap utilities from the AMBER14 suite of computer programs 47 to remove water molecules, to add missing heavy atoms, to remove hydrogen atoms, to identify terminal amino acids and nucleotides, and to rename residues and atoms per the AMBER force field conventions. Solvent molecules and missing hydrogen atoms are then generated using the qprep utility from the Q suite of computer programs. 52 Charges of amino acid and nucleotide residues inside the interior region correspond to their prevalent protonation states in aqueous solution at a neutral pH as determined by the Poisson-Boltzmann model. 70 To ensure the stability of FEP MD simulations, charge states of all residues outside the simulation region are set to zero, and charges of residues near the boundary region are also set to zero unless the residues form close ion pairs. Table S7 lists the charge states of residues in the ternary complex. Four sodium cations of the charge +1 and a sodium ion with a partial charge of +0.6921 are added to the pol λ/DL ternary complexes to reduce the total charge enclosed within the simulation sphere (−1 in systems that include a deprotonated nucleobase or 0 otherwise). To prevent diffusion of sodium ions towards the boundary of the simulation sphere, a flat−2

bottom harmonic potential at a force constant of 5.0 kcal·mol−1 ·Å

restrains their positions.

The potential vanishes for distances less than 20 Å from the center of the simulation sphere 10


Page 11 of 40

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 a 24.5 Å-radius with the center at the initial position of the C4′ atom in dGTP. Thermalization protocols First, solute atoms are restrained to their initial positions by a harmonic potential with a −2

force constant of 200 kcal · mol−1 · Å . Meanwhile, the distribution of solvent molecules is allowed to relax in MD runs at gradually increasing bath temperatures from 5 to 310 K, bath update times from 0.1 to 40 fs, and time steps from 0.01 to 1 fs. Then, the system is quenched back to 5 K and then gradually heated up to the bath temperature of 310 K as the harmonic restraints on positions of solute atoms are gradually lifted. The system is allowed to thermalize in an MD run for 1.3 ns at a 2 fs integration time step with the restraint on the solute switched off. Bonds between hydrogen atoms and atoms within the solute outside the Q-region are set to their equilibrium value by the SHAKE algorithm. 71 The result of this thermalization is an ensemble at the system temperature of 310 K. However, this is not the ensemble we are interested in. We use a Nose-Hoover chain thermostat at a chain length of 10 to ensure our MD simulations sample canonical ensemble. The thermostat couples to all degrees of freedom in a way that it thermalizes the system as a whole to the temperature that corresponds to the average kinetic energy of all atoms. However, the coupling of degrees of freedom within outer and interior regions is weak enough so that the outer region does not equilibrate with the rest of the system at given parameters of our MD simulations. This is to say that the two regions acquire two distinct steady state temperatures which differ from the temperature of 310 K for the whole system that the Nose-Hoover thermostat sets. To construct ensemble that corresponds to the temperature of 310 K for the interior region, we first measure during 100 ps following the thermalization the actual temperature that the interior region develops when the temperature of the bath is set to 310 K. Then we adjust the bath temperature so that the mean temperature within the interior region is 310 K. We then let the system thermalize at this new bath temperature for 650 ps. The bath temperatures of 254.4 K and 11



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 12 of 40

310 K have been used in our simulations of the pol λ/DL ternary complexes and nucleosides in aqueous solution, respectively. Simulation protocols To compute the free energy changes ∆GFp EP (N1 ·N1′ → N2 ·N2′ ) and ∆GFw EP (N1 ·N1′ → N2 ·N2′ ) in Eq. 2, we mutate G · T , G∗ · T , G · T ∗ , G− · T , and G · T − mispairs along the coupling constant paths in fig. 4 in the thermalized ternary complex or aqueous solution. The duration of a trajectory for a potential energy surface at a given coupling constant is 80 ps at the time step size of 1 fs. The bath update time is 40 fs. The SHAKE constraint is applied to all solvent bonds and to all covalent bonds that involve hydrogens outside the Q-region. The total duration of the simulation along the coupling constant path is 31.1 ns. The total energy is stored every 5 fs for further use in free energy calculations. For every coupling constant value, the initial 8.14 ps of data is disregarded in the free energy calculations. Our estimate for the free energy change is the average of the values, which are evaluated along the forward and reverse paths,

∆GF EP (N1 · N1′ → N2 · N2′ ) − ∆GF EP (N1 · N1′ ← N2 · N2′ ) /2. The total simulation time for this estimate is 62.2 ns. Meanwhile, the forward-reverse hysteresis error,

∆GF EP (N1 · N1′ → N2 · N2′ ) + ∆GF EP (N1 · N1′ ← N2 · N2′ ) /2,

(4)

sets a lower bound for the error in the free energy due to incomplete configuration space sampling. The hysteresis error in ∆∆Gw→p (see eq. 2) is the sum of the hysteresis errors for ∆GFp EP and ∆GFw EP .

12


Page 13 of 40

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


Estimates for ∆Gw (N1 · N1′ → N2 · N2′ ) Protonation-deprotonation equilibria in aqueous solution The free energy change ∆Gw (N → N − ) in Eq. 1 for N = G or T in aqueous solution at a given pH is related to the experimental pKa2 constant by ∆Gw (N → N − ) = 2.303 · RT · (pKa2 − pH).

(5)

We use experimental pKa2 values 21,22 in this work since existing quantum chemistry calculations estimate experimental pKa2 for G and T inaccurately. To estimate ∆Gw (N → N − ) at 310 K, we use the experimental pKa2 at the standard conditions and linearly extrapolate these values to 310 K using the experimental standard entropies for deprotonation in aqueous solution. For T at N3 − H and G at N1 − H, the experimental standard acidity constants pKa2 (or the standard deprotonation free energies) are 9.79 (13.4 kcal·mol−1 ) and 9.25 (12.6 kcal·mol−1 ). The respective standard deprotonation entropies are −18.8 cal · mol−1 · K−1 and −16.7 cal · mol−1 · K−1 . 21,22 Keto-enol equilibria in aqueous solution ∆Gw (N → N ∗ ) in Eq. 1 for keto-enol equilibria for N = G or T are not straightforward to access experimentally. Therefore, we compute these quantities. To estimate the free energy of tautomerization in nucleosides, we employ a model chemistry that evaluates tautomerization energies in nucleobases at a reliable ab initio QCISD(T) level of theory 72 and corrects for the transfer of the nucleobases to nucleosides as well as the thermal and solvent effects on a PBE0 level of theory. PBE0 is a balanced and well-tested parameter-free hybrid functional, which is one of the best in its class for a broad range of systems. 73 In the nucleosides, the

13



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 14 of 40

free energy changes upon tautomerization are given by

∆Gw (N → N ∗ ) = ∆EgQCISD(T),base + ∆EwPBE0,base − ∆EgPBE0,base +

(6)

(∆EwPBE0,nucleoside − ∆EwPBE0,base ) +

(∆GPBE0,nucleoside − ∆EwPBE0,nucleoside ) . w Here, ∆EgQCISD(T),base is the difference in the electronic energy upon tautomerization for guanine or thymine nucleobases in the gas phase for which accurate results at the QCISD(T)/TZVP(2df,2pd) level of theory are available. 72 In the second term, the difference between changes in the electronic energy upon tautomerization of the nucleobase in water, ∆EwPBE0,base , and in gas, ∆EgPBE0,base , accounts for solvation effects for the nucleobase at the PBE0/6-311++G(2df,p) level of theory. In the third term, the difference between changes in electronic energies upon tautomerization in water for the nucleoside, ∆EwPBE0,nucleoside , and for the nucleobase, ∆EwPBE0,base , accounts at the PBE0/6-311++G(2df,p) level of theory for the transfer of the nucleobase to the nucleoside. In the fourth term, the difference between the free energy change, ∆GPBE0,nucleoside , and the electronic energy change ∆EwPBE0,nucleoside in the nucleoside w at the PBE0/6-311++G(2df,p) level of theory is the thermal contributions to the free energy changes for the nucleoside. The free energy change ∆EwPBE0,nucleoside is computed from the ideal gas partition function for the nucleosides at 310 K and 1 atm in water with the PBE0/6-311++G(2df,p) geometries and harmonic frequencies. The solvent effects in water (subscript w in eq. 6) here are accounted for within the IEPCM polarizable continuum model as implemented in the Gaussian 09 suite of computer programs. 63 To corroborate our implicit solvent estimates, we independently evaluate ∆Gw (N → N ∗ ) using the FEP MD calculation with explicit solvent. A thermodynamic cycle in fig. 3(b) allows expressing the free energy change as

∆Gw (N → N ∗ ) = ∆Gg (N → N ∗ ) + ∆∆Gg→w (N → N ∗ ).

14


(7)

Page 15 of 40

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


Here, the free energy of tautomerization for an isolated nucleoside in gas is computed as

∆Gg (N → N ∗ ) = ∆EgQCISD(T),base + ∆EgPBE0,nucleoside − ∆EgPBE0,base + ∆GPBE0,nucleoside − ∆EgPBE0,nucleoside . g

(8)

As in eq. 6, the change in the electronic energy upon tautomerization for the nucleobase, ∆EgQCISD(T),base , is corrected for transfer of the nucleobase to the nucleoside in the second term as well as for the thermal effects in the third term at the PBE0/6-311++G(2df,p) level of theory. The term

∆∆Gg→w (N → N ∗ ) = ∆GFw EP (N → N ∗ ) − ∆GFg EP (N → N ∗ )

(9)

is the contribution that we infer from FEP MD simulations in gas and water. In order to compute ∆GFg EP (N → N ∗ ), we evolve an isolated nucleoside molecule whose center of mass is constrained to the center of a simulation sphere by an isotropic harmonic potential 2

with force constant of 50 kcal · mol−1 · Å to prevent the molecule from approaching the boundary of the simulation sphere. The simulation sphere of 24.5 Å radius fully encloses the simulated molecule. The boundary restraints are not included in these single-molecule FEP MD calculations. In our FEP MD calculations of ∆GFw EP (N → N ∗ ), we enclose a solvated nucleoside in the simulation sphere of 24.5 Å radius, which explicit water molecules fill this time. Furthermore, the SCAAS type boundary conditions are imposed in the surface shell. As in the case of the free nucleoside, the molecule’s center of mass is constrained to the center of the simulation sphere by the isotropic harmonic potential. The trajectories starting from the respective thermalized distribution of atomic position and velocities are propagated forward at a 2 fs time step during 12.56 ns, which translates to 0.08 ns per a λ interval in a FEP MD simulation. The bath coupling corresponds to an update time of 40 fs.

15



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

Results and discussion Nucleobases in aqueous solution To compare the thermodynamic stability of the various nucleobase pairs, we first determine the ∆Gw values in eq. 1. The experimental pKa2 for G and T at pH = 7 give rise through eq. 5 to the values of ∆Gw in table 3. Free energy changes upon tautomerization in water ∆Gw in table 3, which eq. 6 defines, are determined by individual contributions in table 2. We note that hybrid density functional theory model underestimates ∆Eg (relative to a more reliable QCISD(T)/TZVP(2df,2pd) 72 ) for guanine and thymine by 0.6 kcal · mol−1 and 0.5 kcal · mol−1 , respectively, which justifies our use of the higher level quantum chemistry for the electronic energies in eq. 6. Furthermore, different methods, the FEP MD with explicit solvent and the DFT with implicit solvent agree very well with each other as table 3 shows. The solvation effects stabilize T ∗ over T , and destabilize G∗ over G with respect to gas phase. Deprotonated forms G− and T − are more stable than the keto forms G and T . However, G and T remain more prevalent than G∗ , T ∗ , G− , and T − in aqueous solution at pH 7 and 310 K as the calculated values of ∆Gw in table 3 indicate.

Nucleobase pairs in pol λ/DL To estimate relative likelihood of G · T , G∗ · T , G · T ∗ , G− · T , G · T − mispairs to occur in the active site of the pol λ/DL, we use eq. 1 to compute values of ∆Gw and ∆∆Gw→p in table 4 along paths in fig. 4. The hysteresis errors for ∆GFw EP and ∆GFp EP (see eq. 4) are also shown in table 4. The hysteresis for ∆GFw EP is much smaller than that in ∆GFp EP as expected since the configurational landscape in the protein is rougher than in the enzyme. To reduce the hysteresis in ∆Gp , we evaluate the free energy difference between G · T − or G− · T and G · T along superposition paths G · T → G∗ · T → G− · T or G · T → G · T ∗ → G · T − . Our estimates of free energy changes in the active site of pol λ/DL in table 4 arrange the 16


Page 16 of 40

Page 17 of 40

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


different mispairs in the following order of decreasing stability: G∗ · T > G · T > G · T ∗ > G− · T > G− · T . The magnitude of the hysteresis error does not allow us to determine reliably whether G∗ · T or G · T is the dominant mispair. We further note that the free energy changes ∆∆Gw→p in table 4 indicate that both G∗ · T and G · T ∗ are significantly stabilized while G− · T and G · T − are destabilized in the active site of pol λ/DL in comparison to the aqueous solution. As table S8 elaborates, the largest contribution to ∆∆Gw→p along the path comes from the block III of the paths where electrostatic interactions change. Moreover, the electrostatic contributions to ∆∆Gw→p are comparable in the magnitude but are the opposite in sign for the tautomerization and deprotonation channels. Thus, electrostatic contributions along the paths are predominantly responsible for stabilization of the pairs that involve G∗ or T ∗ and destabilization of the pairs that involve G− or T − in the active site of pol λ/DL relative to aqueous solution. However, the relative importance of the remaining contributions due to changes vdW and intramolecular terms can be in the order of 20 % for paths that do not involve charged forms. It is further apparent from table S8 that the dominant contribution to the hysteresis error comes from changes in electrostatic contributions. Computed average structures over 1.3 ns trajectories for mispairs G · T , G∗ · T , G− · T , G · T ∗ , and G · T − in the active site of pol λ/DL are compared with the experimental structure 35 in fig. 5. The computed structure for the G · T mispair features an expected wobble geometry. This structure is quite different from the experimental W-C geometry. On the other hand, our calculated structure for G∗ · T has a W-C geometry, which aligns very well with the experimental structure. The average structure for G·T ∗ also arranges in a W-C geometry, but it is noticeably different from the experimental one. The computed structure for a deprotonated mispair G− · T features a W-C geometry, but the arrangement of the nucleobases is less planar than in the case of G∗ ·T and G·T ∗ . However, our structure for G·T − is a wobble geometry rather than a W-C one. Interatomic distances in table 5 and root-meansquare distances in fig. 5 arrange the computed structure from the most to the least similar

17



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

to the experimental mispair in the following order: G∗ · T > G · T ∗ > G− · T > G · T − > G · T . This ordering is parallel to the ordering based on the free energy except for the wobble mispair G · T . Although G∗ · T and G · T appear very similar thermodynamically, significant differences in their shapes are likely to result in significantly different rates of incorporation into the DNA of this mispairs by the polymerase because of the required realignment of the key protein residues that stabilize the chemical transition state. 74,75 Positively charged LYS-25, ARG-261, and ARG-264, as well as the negatively charged GLU-276, in pol λ/DL contribute to the stabilization of the arrangements of nucleobases in the mispairs as fig. 6 illustrates. Arrangements for ARG-264 and ARG-261 do not vary significantly among the mispairs. However, LYS-25, which points towards the guanine, appears relatively mobile and its orientation depends on whether the mispair is negatively charged or neutral. We note that the orientation of LYS-25 towards thymine in the experimental structure does not emerge in our simulations. The other amino acid residues in the computed average structures in fig. 5 reside near their experimental positions.

Conclusions Our free energy calculations place mispairs in the order of decreasing prevalence in the active site of pol λ/DL at a neutral pH: G∗ · T > G · T > G · T ∗ > G− · T > G · T − . In particular, mispairs involving deprotonated forms, G− · T and G · T − , are unlikely in pol λ/DL. We are unable to provide more detailed quantitative relative prevalence estimates for different mispairs since the free energy differences for some mispairs (G · T and G∗ · T ) are less than uncertainties due to configuration space sampling, errors in potential energy surfaces and quantum chemistry calculations. Major methodological advances outside of the scope of this report would be needed to improve upon these predictions. Therefore, as an independent evidence that G∗ · T rather than G · T is a dominant mispair, we rely on the comparison between the experimental 35 and our computed mispair geometries, which

18


Page 18 of 40

Page 19 of 40

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


indicates that the experimental and computed structures are the most similar in the case of the G∗ · T mispair. This report supports the notion that unusual forms of nucleobases can be significantly stabilized in the active site of pol λ/DL which offers environment significantly different from the gas phase or aqueous solution.

Supporting Information Available • Figure S1: Atom names in nucleosides, nucleotides, and nucleoside triphosphates. • Tables S1 to S6: Atom types and RESP charges. • Table S7: Charge states of residues in the pol λ/DL complex with the gapped DNA. • Table S8: Contributions to ∆∆Gw→p along different blocks of FEP paths. • Average structures for mispairs G · T , G∗ · T , G− · T , G · T ∗ , G · T − in the active site of pol λ/DL: GT.pdb, GsT.pdb, GmT.pdb, GTs.pdb, GTm.pdb, respectively. This material is available free of charge via the Internet at http://pubs.acs.org/.

Acknowledgement This work was supported by the National Cancer Institute under grant 5U19CA177547. The calculations in this work were performed under grant of computer time from the Swedish National Allocations Committee (SNAC) through SNAC projects SNIC 2015/16-12 and SNIC 2016/34-27.

19



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

References (1) Watson, J.; Crick, F. Genetical Implications of the Structure of Deoxyribonucleic Acid. Nature 1953, 171, 964–967. (2) Löwdin, P.-O. Proton Tunneling in DNA and its Biological Implications. Rev. Mod. Phys. 1963, 35, 724–732. (3) Saenger, W. Principles of Nucleic Acid Structure; Springer advanced texts in chemistry; Springer-Verlag, 1984. (4) Morales, J. C.; Kool, E. T. Efficient Replication between non-Hydrogen-Bonded Nucleoside Shape Analogs. Nature Structural Biology 1998, 5, 950–954. (5) Goodman, M. F. On the Wagon – DNA Polymerase Joins "H-Bonds Anonymous". Nature Biotechnology 1999, 17, 640–641. (6) Kool, E. T. Active Site Tightness and Substrate Fit in DNA Replication. Annu. Rev. Biochem. 2002, 71, 191–219. (7) Florián, J.; Šponer, J.; Warshel, A. Thermodynamic Parameters for Stacking and Hydrogen Bonding of Nucleic Acid Bases in Aqueous Solution: Ab Initio/Langevin Dipoles Study. The Journal of Physical Chemistry B 1999, 103, 884–892. (8) Hobza, P.; Šponer, J. Structure, Energetics, and Dynamics of the Nucleic Acid Base Pairs: Nonempirical Ab Initio Calculations. Chemical Reviews 1999, 99, 3247–3276. (9) Florián, J.; Warshel, A.; Goodman, M. F. Molecular Dynamics Free-Energy Simulations of the Binding Contribution to the Fidelity of T7 DNA Polymerase. The Journal of Physical Chemistry B 2002, 106, 5754–5760. (10) Goodman, M. F.; Creighton, S.; Bloom, L. B.; Petruska, J.; Kunkel, T. A. Biochemical Basis of DNA Replication Fidelity. Critical Reviews in Biochemistry and Molecular Biology 1993, 28, 83–126. 20


Page 20 of 40

Page 21 of 40

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


(11) Topal, M. D.; Fresco, J. R. Complementary Base Pairing and the Origin of Substitution Mutations. Nature 1976, 263, 285–289. (12) Lawley, P.; Brookes, P. Ionization of DNA Bases or Base Analogues as a Possible Explanation of Mutagenesis, with Special Reference to 5-Bromodeoxyuridine. Journal of Molecular Biology 1962, 4, 216–219. (13) Sowers, L. C.; Shaw, B. R.; Veigl, M. L.; Sedwick, W. D. DNA Base Modification: Ionized Base Pairs and Mutagenesis. Mutation Research/Fundamental and Molecular Mechanisms of Mutagenesis 1987, 177, 201–218. (14) Kennard, O. Structural Studies of DNA Fragments: The G·T Wobble Base Pair in A, B and Z DNA: The G·A Base Pair in B-DNA. Journal of Biomolecular Structure and Dynamics 1985, 3, 205–226. (15) de Vries, M. S. Gas-Phase IR Spectroscopy of Nucleobases. Topics in Current Chemistry 2014, 271–297. (16) Florián, J.; Leszczyński, J. Spontaneous DNA Mutations Induced by Proton Transfer in the Guanine·Cytosine Base Pairs: An Energetic Perspective. Journal of the American Chemical Society 1996, 118, 3010–3017. (17) de Vries, M. S.; Hobza, P. Gas-Phase Spectroscopy of Biomolecular Building Blocks. Annual Review of Physical Chemistry 2007, 58, 585–612. (18) Shukla, M. K.; Leszczynski, J. Tautomerism in Nucleic Acid Bases and Base Pairs: A brief Overview. Wiley Interdisciplinary Reviews: Computational Molecular Science 2013, 3, 637–649. (19) Poater, J.; Swart, M.; Bickelhaupt, F. M.; Fonseca Guerra, C. B-DNA Structure and Stability: The Role of Hydrogen Bonding, π–π Stacking Interactions, Twist-angle, and Solvation. Org. Biomol. Chem. 2014, 12, 4691–4700. 21



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

(20) Poater, J.; Swart, M.; Fonseca Guerra, C.; Bickelhaupt, F. M. Solvent Effects on Hydrogen Bonds in Watson–Crick, Mismatched, and Modified DNA Base Pairs. Computational and Theoretical Chemistry 2012, 998, 57–63. (21) Christensen, J. J.; Rytting, J. H.; Izatt, R. M. Thermodynamics of Proton Dissociation in Dilute Aqueous Solution. VII. pK, ∆H 0 , and ∆S 0 Values for Proton Ionization from Several Pyrimidines and Their Nucleosides at 250 . The Journal of Physical Chemistry 1967, 71, 2700–2705. (22) Christensen, J. J.; Rytting, J. H.; Izatt, R. M. Thermodynamic pK, ∆H 0 , ∆S0 , and ∆Cp 0 Values for Proton Dissociation from Several Purines and Their Nucleosides in Aqueous Solution. Biochemistry 1970, 9, 4907–4913. (23) Kampf, G.; Kapinos, L. E.; Griesser, R.; Lippert, B.; Sigel, H. Comparison of the Acid–Base Properties of Purine Derivatives in Aqueous Solution. Determination of Intrinsic Proton Affinities of Various Basic Sites. Journal of the Chemical Society, Perkin Transactions 2 2002, 1320–1327. (24) Mucha, A.; Knobloch, B.; Jezowska-Bojczuk, M.; Kozlowski, H.; Sigel, R. Comparison of the Acid-Base Properties of Ribose and 2′ -Deoxyribose Nucleotides. Chemistry - A European Journal 2008, 14, 6663–6671. (25) Jang, Y. H.; Goddard, W. A.; Noyes, K. T.; Sowers, L. C.; Hwang, S.; Chung, D. S. pKa Values of Guanine in Water: Density Functional Theory Calculations Combined with Poisson-Boltzmann Continuum-Solvation Model. The Journal of Physical Chemistry B 2003, 107, 344–357. (26) Thapa, B.; Schlegel, H. B. Calculations of pKa’s and Redox Potentials of Nucleobases with Explicit Waters and Polarizable Continuum Solvation. J. Phys. Chem. A 2015, 119, 5134–5144.

22


Page 22 of 40

Page 23 of 40

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


(27) Brovarets’, O. O.; Hovorun, D. M. The Nature of the Transition Mismatches with Watson–Crick Architecture: The G∗ ·T or G·T∗ DNA Base Mispair or Both? A QM/QTAIM Perspective for the Biological problem. Journal of Biomolecular Structure and Dynamics 2014, 33, 925–945. (28) Sowers, L. C.; Eritja, R.; Kaplan, B.; Goodman, M. F.; Fazakerly, G. V. Equilibrium between a Wobble and Ionized Base Pair Formed between Fluorouracil and Guanine in DNA as Studied by Proton and Fluorine NMR. Journal of Biological Chemistry 1988, 263, 14794–14801. (29) Kimsey, I. J.; Petzold, K.; Sathyamoorthy, B.; Stein, Z. W.; Al-Hashimi, H. M. Visualizing Transient Watson–Crick-like Mispairs in DNA and RNA Duplexes. Nature 2015, 519, 315–320. (30) Warshel, A. Energetics of Enzyme Catalysis. Proceedings of the National Academy of Sciences 1978, 75, 5250–5254. (31) Warshel, A.; Florian, J. Computer Simulations of Enzyme Catalysis: Finding Out What has been Optimized by Evolution. Proceedings of the National Academy of Sciences 1998, 95, 5950–5955. (32) Xia, S.; Konigsberg, W. H. Mispairs with Watson-Crick Base-Pair Geometry Observed in Ternary Complexes of an RB69 DNA Polymerase Variant. Protein Science 2014, 23, 508–513. (33) Xia, S.; Konigsberg, W. H. RB69 DNA Polymerase Structure, Kinetics, and Fidelity. Biochemistry 2014, 53, 2752–2767. (34) Kimsey, I.; Al-Hashimi, H. M. Increasing Occurrences and Functional Roles for High Energy Purine-Pyrimidine Base-Pairs in Nucleic Acids. Current Opinion in Structural Biology 2014, 24, 72–80.

23



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

(35) Bebenek, K.; Pedersen, L. C.; Kunkel, T. A. Replication Infidelity via a Mismatch with Watson–Crick Geometry. Proceedings of the National Academy of Sciences 2011, 108, 1862–1867. (36) Wang, W.; Hellinga, H. W.; Beese, L. S. Structural Evidence for the Rare Tautomer Hypothesis of Spontaneous Mutagenesis. Proceedings of the National Academy of Sciences 2011, 108, 17644–17648. (37) Koag, M.-C.; Nam, K.; Lee, S. The Spontaneous Replication Error and the Mismatch Discrimination Mechanisms of Human DNA Polymerase β. Nucleic Acids Research 2014, 42, 11233–11245. (38) Singh, V.; Fedeles, B. I.; Essigmann, J. M. Role of Tautomerism in RNA Biochemistry. RNA 2014, 21, 1–13. (39) Rozov, A.; Demeshkina, N.; Westhof, E.; Yusupov, M.; Yusupova, G. Structural Insights into the Translational Infidelity Mechanism. Nature Communications 2015, 6, 7251. (40) Yu, H.; Eritja, R.; Bloom, L. B.; Goodman, M. F. Ionization of Bromouracil and Fluorouracil Stimulates Base Mispairing Frequencies with Guanine. Journal of Biological Chemistry 1993, 268, 15935–43. (41) Satpati, P.; Åqvist, J. Why Base Tautomerization does not Cause Errors in mRNA Decoding on the Ribosome. Nucleic Acids Research 2014, 42, 12876–12884. (42) Bebenek, K.; Pedersen, L. C.; Kunkel, T. A. Structure–Function Studies of DNA Polymerase λ. Biochemistry 2014, 53, 2781–2792. (43) Kunkel, T. A.; Bebenek, K. DNA Replication Fidelity. Annu. Rev. Biochem. 2000, 69, 497–529. (44) Uchiyama, Y.; Takeuchi, R.; Kodera, H.; Sakaguchi, K. Distribution and Roles of Xfamily DNA Polymerases in Eukaryotes. Biochimie 2009, 91, 165–170. 24


Page 24 of 40

Page 25 of 40

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


(45) Bebenek, K.; Garcia-Diaz, M.; Zhou, R.-Z.; Povirk, L. F.; Kunkel, T. A. Loop 1 Modulates the Fidelity of DNA Polymerase λ. Nucleic Acids Research 2010, 38, 5419–5431. (46) Klvaňa, M.; Bren, U.; Florián, J. Uniform Free-Energy Profiles of the P–O Bond Formation and Cleavage Reactions Catalyzed by DNA Polymerases β and λ. The Journal of Physical Chemistry B 2016, 120, 13017–13030. (47) Case, D.; Babin, V.; Berryman, J.; Betz, R.; Cai, Q.; Cerutti, D.; Cheatham, T.; III,; Darden, T.; Duke, R. et al. AMBER 14. University of California, San Francisco 2014. (48) Zwanzig, R. W. High-Temperature Equation of State by a Perturbation Method. I. Nonpolar Gases. The Journal of Chemical Physics 1954, 22, 1420. (49) Jorgensen, W. L.; Ravimohan, C. Monte Carlo Simulation of Differences in Free Energies of Hydration. The Journal of Chemical Physics 1985, 83, 3050–3054. (50) Beutler, T. C.; Mark, A. E.; van Schaik, R. C.; Gerber, P. R.; van Gunsteren, W. F. Avoiding Singularities and Numerical Instabilities in Free Energy Calculations Based on Molecular Simulations. Chemical Physics Letters 1994, 222, 529–539. (51) Boukharta, L.; de Terán, H. G.; Åqvist, J. Computational Prediction of Alanine Scanning and Ligand Binding Energetics in G-Protein Coupled Receptors. PLoS Computational Biology 2014, 10, e1003585. (52) Marelius, J.; Kolmodin, K.; Feierberg, I.; Åqvist, J. Q: A Molecular Dynamics Program for Free Energy Calculations and Empirical Valence Bond Simulations in Biomolecular Systems. Journal of Molecular Graphics and Modelling 1998, 16, 213–225. (53) Hornak, V.; Abel, R.; Okur, A.; Strockbine, B.; Roitberg, A.; Simmerling, C. Comparison of Multiple Amber Force Fields and Development of Improved Protein Backbone Parameters. Proteins: Structure, Function, and Bioinformatics 2006, 65, 712–725.

25



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

(54) 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. Journal of the American Chemical Society 1995, 117, 5179–5197. (55) Cheatham, T. E.; Cieplak, P.; Kollman, P. A. A Modified Version of the Cornell et al. Force Field with Improved Sugar Pucker Phases and Helical Repeat. Journal of Biomolecular Structure and Dynamics 1999, 16, 845–862. (56) Pérez, A.; Marchán, I.; Svozil, D.; Sponer, J.; Cheatham, T. E.; Laughton, C. A.; Orozco, M. Refinement of the AMBER Force Field for Nucleic Acids: Improving the Description of α/γ Conformers. Biophysical Journal 2007, 92, 3817–3829. (57) Zgarbová, M.; Luque, F. J.; Šponer, J.; Cheatham, T. E.; Otyepka, M.; Jurečka, P. Toward Improved Description of DNA Backbone: Revisiting Epsilon and Zeta Torsion Force Field Parameters. Journal of Chemical Theory and Computation 2013, 9, 2339– 2354. (58) Krepl, M.; Zgarbová, M.; Stadlbauer, P.; Otyepka, M.; Banáš, P.; Koča, J.; Cheatham, T. E.; Jurečka, P.; Šponer, J. Reference Simulations of Noncanonical Nucleic Acids with Different χ Variants of the AMBER Force Field: Quadruplex DNA, Quadruplex RNA, and Z-DNA. Journal of Chemical Theory and Computation 2012, 8, 2506–2520. (59) Wang, J.; Wolf, R. M.; Caldwell, J. W.; Kollman, P. A.; Case, D. A. Development and Testing of a General Amber Force Field. Journal of Computational Chemistry 2004, 25, 1157–1174. (60) Wolf, M. G.; Groenhof, G. Evaluating Nonpolarizable Nucleic Acid Force Fields: A Systematic Comparison of the Nucleobases Hydration Free Energies and Chloroform-

26


Page 26 of 40

Page 27 of 40

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


to-Water Partition Coefficients. Journal of Computational Chemistry 2012, 33, 2225– 2232. (61) Cieplak, P.; Cornell, W. D.; Bayly, C.; Kollman, P. A. Application of the Multimolecule and Multiconformational RESP Methodology to Biopolymers: Charge Derivation for DNA, RNA, and Proteins. Journal of Computational Chemistry 1995, 16, 1357–1377. (62) Meagher, K. L.; Redman, L. T.; Carlson, H. A. Development of Polyphosphate Parameters for Use with the AMBER Force Field. Journal of Computational Chemistry 2003, 24, 1016–1025. (63) Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A.; Cheeseman, J. R.; Scalmani, G.; Barone, V.; Mennucci, B.; Petersson, G. A. et al. Gaussian 09 Revision E.01. Gaussian Inc. Wallingford CT 2009. (64) Florián, J.; Goodman, M. F.; Warshel, A. Theoretical Investigation of the Binding Free Energies and Key Substrate-Recognition Components of the Replication Fidelity of Human DNA Polymerase β. The Journal of Physical Chemistry B 2002, 106, 5739– 5753. (65) Florián, J.; Goodman, M. F.; Warshel, A. Computer Simulation of the Chemical Catalysis of DNA Polymerases:

Discriminating between Alternative Nucleotide Insertion

Mechanisms for T7 DNA Polymerase. Journal of the American Chemical Society 2003, 125, 8163–8177. (66) Jorgensen, W. L.; Chandrasekhar, J.; Madura, J. D.; Impey, R. W.; Klein, M. L. Comparison of Simple Potential Functions for Simulating Liquid Water. The Journal of Chemical Physics 1983, 79, 926–935. (67) King, G.; Warshel, A. A Surface Constrained All-atom Solvent Model for Effective Simulations of Polar Solutions. The Journal of Chemical Physics 1989, 91, 3647–3661.

27



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

(68) Sham, Y. Y.; Warshel, A. The Surface Constraint All Atom Model Provides Size Independent Results in Calculations of Hydration Free Energies. The Journal of Chemical Physics 1998, 109, 7940–7944. (69) Lee, F. S.; Warshel, A. A Local Reaction Field Method for Fast Evaluation of Longrange Electrostatic Interactions in Molecular Simulations. The Journal of Chemical Physics 1992, 97, 3100–3107. (70) Anandakrishnan, R.; Aguilar, B.; Onufriev, A. V. H++ 3.0: Automating pK Prediction and the Preparation of Biomolecular Structures for Atomistic Molecular Modeling and Simulations. Nucleic Acids Research 2012, 40, W537–W541. (71) Ryckaert, J.-P.; Ciccotti, G.; Berendsen, H. J. Numerical Integration of the Cartesian Equations of Motion of a System with Constraints: Molecular Dynamics of n-Alkanes. Journal of Computational Physics 1977, 23, 327–341. (72) Piacenza, M.; Grimme, S. Systematic Quantum Chemical Study of DNA-Base Tautomers. Journal of Computational Chemistry 2004, 25, 83–99. (73) Staroverov, V. N.; Scuseria, G. E.; Tao, J.; Perdew, J. P. Comparative Assessment of a New Nonempirical Density Functional: Molecules and Hydrogen-bonded Complexes. The Journal of Chemical Physics 2003, 119, 12129–12137. (74) Xiang, Y.; Oelschlaeger, P.; Florián, J.; Goodman, M. F.; Warshel, A. Simulating the Effect of DNA Polymerase Mutations on Transition-State Energetics and Fidelity: Evaluating Amino Acid Group Contribution and Allosteric Coupling for Ionized Residues in Human Pol β. Biochemistry 2006, 45, 7036–7048. (75) Klvaňa, M.; Jeřábek, P.; Goodman, M. F.; Florián, J. An Abridged Transition State Model To Derive Structure, Dynamics, and Energy Components of DNA Polymerase β Fidelity. Biochemistry 2011, 50, 7023–7032.

28


Page 28 of 40

Page 29 of 40


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 ACS Paragon Plus Environment



Page 30 of 40

Page 31 of 40




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: Paths connecting different mispairs in FEP MD calculations.

32


Page 32 of 40

Page 33 of 40

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


G · T (rmsd 0.73 Å)

G∗ · T (rmsd 0.22 Å)

G− · T (rmsd 0.33 Å)

G · T ∗ (rmsd 0.33 Å)

G · T − (rmsd 0.42 Å) Figure 5: Comparison of the experimental and computed average structures of nucleobase pairs in the active site of pol λ/DL. Only nucleobases are shown in the figure. The heavy atoms in the experimental and computed structures for the guanine moiety are aligned. The root-mean-square distances involving the nucleobase atoms as well as the C′1 atom of deoxyribose are indicated in the parenthesis. Experimental structures are shown in ochre. 33




Page 34 of 40

Page 35 of 40

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: Deformation of the potential along the path in the λ-space in FEP MD simulations. X denotes a linear adjustment in the parameter it corresponds to. D and C refer to the changes during destruction and creation of atoms, respectively. Block I II III IV V Number of Windows 17 21 321 21 17 Charges XCD vdW XC D C D Soft-core vdW Bonds X Angles X Torsions X X Improper torsions

35



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 36 of 40

Table 2: Electronic energy changes in gas (∆Eg ) and water (∆Ew ) as well as the free energy changes in gas (∆Gg ) and water (∆Gw ) for the nucleobases and nucleosides. ∆Eg a ∆Eg b ∆Gg b ∆Ew c ∆Gw c

G → G∗d 0.10 0.74 0.78 6.84 6.82

G → G∗e 0.50 0.60 6.77 6.79

T → T ∗d 12.10 12.63 12.58 11.27 11.18

T → T ∗e 13.16 12.79 12.52 11.94

∆EgQCISD(T),base (N → N ∗ ) are computed in the QCISD(T)/TZVP(2df,2pd) approximation. 72 b ∆EgPBE0,base and ∆EgPBE0,nucleoside are computed in the PBE0/6-311++G(2df,p) approximation. c ∆EwPBE0,base , ∆EwPBE0,nucleoside , ∆GPBE0,base , and ∆GPBE0,nucleoside are computed in the w w PBE0/6-311++G(2df,p) approximation with the IEPCM solvent model. d Values refer to a nucleobase. e Values refer to a nucleoside.

a

36


Page 37 of 40

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 3: Free energy changes ∆Gg and ∆Gw and relative free energy changes ∆∆Gg→w upon tautomerization or deprotonation of the nucleosides (in kcal · mol−1 ) at pH 7. dG → dG

∗

dT → dT ∗

∆Gag 0.0 12.3

dG → dG− dT → dT −

∆∆Gg→w 6.1b 6.2c −0.8b −0.4c

∆Gw 6.1 6.2 11.4 11.9 3.3d 4.0d

The values of free energies at 310 K and 1 atm are computed for free nucleosides at the PBE0/6-311++G(2df,p) level of theory. 63 b The values are ∆GIEPCM − ∆Gg , where ∆GIEPCM is the free energy change for the free w w nucleosides computed in the presence of the IEPCM self-consistent reaction field at the PBE0/6-311++G(2df,p) level of theory. 63 c The differences in tautomerization free energies for the free nucleosides and nucleosides in explicit water are computed using FEP MD with the AMBER force field. 47 d An estimate for the difference in the deprotonation free energies for the nucleosides at 310 K. a

37


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


38

∆∆Gw→p −6.27 −1.03 9.64 −9.57 19.86 5.04

Hysteresis in ∆∆Gw→p a ±1.22 ±2.98 ±1.18 ±0.62 ±0.58 ±0.61

The hysteresis error is computed using eq. 4.

G·T →G ·T G · T → G− · T G · T → G · T− G · T → G · T∗ G · T∗ → G · T− G∗ · T → G − · T

∗

∆Gw 6.15 3.27 4.04 11.42 −7.38 −2.88

Hysteresis in ∆GFw EP a ±0.04 ±0.06 ±0.02 ±0.03 ±0.02 ±0.07

Hysteresis in ∆GFp EP a ±1.18 ±2.92 ±1.16 ±0.59 ±0.55 ±0.54

∆Gp −0.17 2.24 13.68 1.83 12.50 3.91

Table 4: Free energy changes upon tautomerization or ionization of nucleobase pairs in aqueous solution and the active site of pol λ/DL (in kcal · mol−1 ) at pH 7 and 310 K.

The Journal of Physical Chemistry Page 38 of 40

Page 39 of 40

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 5: Structural parameters for G · T , G∗ · T , G− · T , G · T − , and G · T ∗ mispairs (see fig. 1) in the active site of pol λ/DL. Distances are in Å and the angles are in degrees. N1 − O2 O6 − N3 N1 − N3 O6 − O4 N2 − O2 C′1 − C′1 6 N9 C′ C′ 1 1 6 N1 C′ C′ 1 1

G·T 2.74 3.11 3.55 3.71 3.12 10.20 49.07 75.30

G− · T 3.62 3.80 3.02 3.11 2.88 10.64 55.41 59.35

G∗ · T 3.67 3.76 2.90 2.71 2.92 10.67 54.01 57.04

39

G · T− 4.41 4.55 3.15 3.22 2.90 10.65 65.12 54.29

G · T∗ 3.24 3.16 2.90 2.74 2.86 10.58 51.78 59.76


Experiment 4.01 3.83 2.97 2.66 3.06 10.80 54.30 50.50


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

Graphical TOC Entry

40


Page 40 of 40

DNA Polymerase Î» Active Site Favors a Mutagenic Mispair Between

Recommend Documents