Mechanism of the Glycosidic Bond Cleavage of Mismatched Thymine

Aug 31, 2015 - (3) A detailed comprehension of the BER mechanism will thus help us to understand important disorders related to DNA damage and may hel...
0 downloads 15 Views 4MB Size
Article pubs.acs.org/JPCB

Mechanism of the Glycosidic Bond Cleavage of Mismatched Thymine in Human Thymine DNA Glycosylase Revealed by Classical Molecular Dynamics and Quantum Mechanical/Molecular Mechanical Calculations Natalia Kanaan,† Ramon Crehuet,‡ and Petra Imhof*,† †

Institute of Theoretical Physics, Free University Berlin, 14195, Berlin, Germany Institute of Advanced Chemistry of Catalonia (IQAC), CSIC, c/Jordi Girona 18-26, Barcelona 08034, Spain



S Supporting Information *

ABSTRACT: Base excision of mismatched or damaged nucleotides catalyzed by glycosylase enzymes is the first step of the base excision repair system, a machinery preserving the integrity of DNA. Thymine DNA glycosylase recognizes and removes mismatched thymine by cleaving the C1′−N1 bond between the base and the sugar ring. Our quantum mechanical/molecular mechanical calculations of this reaction in human thymine DNA glycosylase reveal a requirement for a positive charge in the active site to facilitate C1′−N1 bond scission: protonation of His151 significantly lowers the free energy barrier for C1′−N1 bond dissociation compared to the situation with neutral His151. Shuttling a proton from His151 to the thymine base further reduces the activation free energy for glycosidic bond cleavage. Classical molecular dynamics simulations of the H151A mutant suggest that the mutation to the smaller, neutral, residue increases the water accessibility of the thymine base, rendering direct proton transfer from the bulk feasible. Quantum mechanical/molecular mechanical calculations of the glycosidic bond cleavage reaction in the H151A mutant show that the activation free energy is slightly lower than in the wild-type enzyme, explaining the experimentally observed higher reaction rates in this mutant.



INTRODUCTION The base excision repair (BER) system is one of the machineries protecting the genome stability by removing damaged or mispaired DNA,1 arising from deaminations, oxidations, or alkylations.2 In the first step of the BER pathway, a DNA glycosylase recognizes the damaged or mispaired nucleobase and flips it out of the DNA strand, burying it into its active site.1 Then, the glycosidic C1′−N1 bond is cleaved, and the base is expelled.1 After that, according to the currently accepted model of the core BER pathway, the AP endonuclease APE1 cleaves the DNA backbone at the abasic site.1 In a fourth step, polymerase β fills the existing gap with the correct nucleotide, and finally, a DNA ligase reanneals the broken DNA strand.1 Damaged or mismatched bases can cause mutations and give rise to genetic diseases and cancer.3 A detailed comprehension of the BER mechanism will thus help us to understand © 2015 American Chemical Society

important disorders related to DNA damage and may help in developing new ways to selectively improve its performance. In this paper, we focus on the glycosidic bond cleavage performed by the monofunctional pyrimidine specific glycosylase thymine DNA glycosylase (TDG). TDG belongs to the same superfamily as uracil DNA glycosylase (UDG) and mismatch-specific uracil DNA glycosylase (MUG).4 TDG has a broad range of substrates that includes oxidation, alkylation and deamination products of methylcytosine, thymine, and adenine, implicating a rather general function in the repair of DNA base damage.5 In particular, it can excise uridine, thymine, 5-hydroxymethyluridine, 5-formylcytosine, and 5-carboxylcytosine when paired to guanine.5,6 The primary Received: June 9, 2015 Revised: August 27, 2015 Published: August 31, 2015 12365

DOI: 10.1021/acs.jpcb.5b05496 J. Phys. Chem. B 2015, 119, 12365−12380

Article

The Journal of Physical Chemistry B Scheme 1. Proposed Mechanism of the Glycosidic Bond Cleavage in the TDG (Adapted from Reference 8)a

a

Note that the gray arrow indicating the water attack applies only to the concerted mechanism shown in the upper route.

to discrimination of G:T and G:C) is a possible explanation for TDG’s weak activity compared to other glycosylases.4,16 The active site of TDG lacks a residue that could activate the nucleophile similar to the role of Asp145 in UDG that has been proposed to be the activator in that enzyme.4 In TDG, Asn140, is located at approximately the same site as Asp145 in UDG. According to Drohat et al.,14 Asn140 plays a major role in catalysis, possibly by positioning the attacking water molecule. Furthermore, based on a crystal structure of TDG with a substrate analog Thr197 has been suggested to help positioning Asn140, explaining the 32-fold lower activity of the T197A mutant compared to the native TDG.16 Tyr152 has been discussed to help in the positioning of the flipped nucleotide.17 The mutation of Asn191 to Ala causes a 15-fold loss in G:T activity. On the contrary, the mutations H151A and A145G have shown to improve the catalytic activity of TDG on G:T in 13-fold. Moreover, these mutants demonstrate an increased activity on A:T, a fact that could explain the balance reached by this enzyme: it has a weak G:T activity, but it avoids the aberrant activity on A:T.16 To shed more light on this important but still not fully understood enzyme, we embarked on exploring the molecular mechanism by which TDG excises the mispaired thymine in DNA. This paper reports molecular dynamics simulations of wild-type and mutant TDG in complex with G:T mismatched DNA, and quantum mechanical/molecular mechanical studies of the glycosidic bond cleavage reaction, elucidating mechanistic details and the catalytic role of various active site residues.

target of TDG is mismatched thymine from guanine−thymine (G:T) mispairs, which occurs as a result of methylcytosine deamination.3 Two mechanisms have been discussed for glycosidic bond cleavage reaction in the literature:4,7 an associative SN2 mechanism involving a concerted addition of a water nucleophile and expulsion of the leaving group (Scheme 1, upper intermediate4); or a stepwise mechanism involving the formation of a discrete oxacarbenium ion intermediate that is subsequently trapped by a water molecule (Scheme 1, lower intermediate7). The first glycosylase to be discovered9 was uracil DNA glycosylase (UDG), and as of today it is the most studied glycosylase. UDG has been extensively studied by molecular simulations (small models using quantum mechanical computations8,10 and quantum mechanical/molecular mechanical (QM/MM) calculations of the enzymatic reaction11). According to a QM/MM study,11 UDG follows a stepwise mechanism in which two residues, Asp145 and His268, stabilize the oxacarbenium ion and the anionic base, respectively. In this study another histidine residue, His148, was found to destabilize the transition state. Based on experimental kinetic isotope effects,4 and in agreement with the computational studies,11 a highly dissociative reaction mechanism involving an oxacarbenium cation/uracil anion pair intermediate has been suggested for UDG.8 His148 is reported to contribute in helping to position the flipped uracil or thymine base, then together with Asp145, orient the nucleophilic water molecule, and finally electrostatically stabilize the glycosyl cation. Moreover, it is conceivable that His148 changes its protonation state from neutral in the free UDG enzyme, as NMR data shows,12,13 to protonated, as indicated by the crystal structure of the complex with trapped glycosyl cation intermediate.14,15 For the removal of the “normal” DNA base thymine from a mismatch, thymine DNA glycosylase must recognize the target base taking into account the complementary base on the other strand. This extra burden on specificity, discriminating thymine in G:T mismatches over thymine in A:T matches (in addition



METHODS Model Setup. The initial structure chosen for carrying out this study is the X-ray crystal structure of the human thymine DNA glycosylase bound to a substrate analogue 2′-fluoro-2′deoxyuridine already flipped (but not cleaved) in the active site of the enzyme (PDB ID code 3UFJ, with a resolution of 2.97 Å).16 The substrate analogue was replaced in silico by thymine. In the computational model, a protein chain of 185 amino acids and a DNA strand of 11 bp length were included. The hydrogen atoms were placed by psfgen (VMD)18−20 and the 12366

DOI: 10.1021/acs.jpcb.5b05496 J. Phys. Chem. B 2015, 119, 12365−12380

Article

The Journal of Physical Chemistry B

Figure 1. Snapshot of the MD simulation of the wild-type (a) and H151A mutant (b) TDG-DNA complex used in subsequent QM/MM calculations. The QM part is depicted with sticks; the protein (gray) and the DNA (orange) are drawn with ribbons.

ensemble in which the pressure was maintained using the Nosé−Hoover Langevin piston with a decay period of 500 fs. Finally, three MD production runs of 100 ns each (time step 2 fs) were carried out started with different initial velocities. In order to reduce the impact of the initial overlay of water molecules with the solute atoms (protein and DNA) on the hydration of the active site, we performed additional MD simulations where the solute has been “resolvated” after an initial equilibration phase of 100 ps. “Resolvation” means that all water and ion atoms have been deleted and newly placed with the above conditions around the new, equilibrated conformation of the solute. Then, another round of minimization, heating, and equilibration was conducted. Those “resolvated” models were then also subjected to 100 ns long MD production runs. This procedure has been repeated three times. The model for the H151A TDG mutant was obtained following the same steps as for the native model, but replacing histidine 151 by alanine (before solvating and energy minimizing the system). To set up the N140A mutant model we replaced Asn140 by alanine. The T197A mutant model was built replacing Thr197 by alanine. The subsequent protocol of solvation, minimization, heating and relaxation was the same as the one used for the models of the native TDG-DNA complex and for the H151A mutant. For all mutant models, H151A, N140A, and T197A, we again performed three MD production runs with different initial velocities. pKa Calculations. The protonation state of His151 has been further investigated by pKa calculations. To explore the conformation dependence of the histidine’s protonation probabilities, we took ten snapshots of each of the classical molecular dynamics simulations (see above) of wild-type TDG with neutral (HSD and HSE) and protonated His151 (HSP). For each snapshot, the difference in the electrostatic contribution of the solvation free energy between a neutral and a charged histidine residue was evaluated by solving the Poisson−Boltzmann equation for both forms, neutral and charged. The same way the difference in the electrostatic component of the solvation free energy between protonated and neutral form was calculated for the isolated residue. The difference between the solvation free energy differences of

protonation state of all the titratable amino acids was computed by means of the propKa program of Jensen et al.21 All of the amino acids, except for His151, are treated in their standard protonation state. His151 points toward bulk water, and it is therefore likely to fluctuate between neutral forms (one proton either at Nδ (HSD), or one proton at Nϵ (HSE)), and the protonated form (one proton at Nδ and a second proton at Nϵ), see Figure S23. Classical Molecular Dynamics Simulations. We have conducted molecular dynamics simulations of wild-type TDG with neutral histidine His151 (named “HSD” in the model with the hydrogen atom bound to the Nδ atom and “HSE” in the model with the hydrogen bond atom bound to the Nϵ atom, see Figure S23 and Figure 2) and protonated His151 (HSP) and of the H151A mutant. Further classical molecular dynamics simulations were performed for the N140A mutant and the T197A mutant. After the construction of the model, a short conjugate gradient minimization (1000 steps) was performed using NAMD2.922 to avoid the hydrogen clashes. The ends of the DNA strand were restrained such that the centers of mass of the first and last DNA bases were kept around a distance of 3 Å. This restraint was maintained for all subsequent calculations. Next, the system was solvated in a cubic box of water molecules (90 Å side length), with a total of about 22 000 water molecules from which all the water molecules closer than 2.4 Å from any solute atom were removed. To obtain a neutral charge of the system, a total of 12 sodium ions were placed randomly with a minimum distance of 10.5 Å from the solute and with a distance of 5 Å between sodium ions. Charmm27 parameters and the TIP3P water model were used to describe the force field.23,24 Nonbonding interactions were treated by means of periodic boundary conditions (switch function with a cutoff of 14−12 Å and a pair list of 16 Å). Particle-mesh Ewald (PME25) was used to calculate the periodic electrostatic interactions. The solvated (and neutralized) system was first energy minimized (5000 steps of conjugate gradient with an energy tolerance of 10−4 kcal/mol) followed by a molecular dynamics (MD) simulation of 30 ps (time step 1 fs) to heat up the system to 295 K by velocity scaling. After that, a relaxation MD of 100 ps at 295 K (time step 1 fs) was computed for an NPT 12367

DOI: 10.1021/acs.jpcb.5b05496 J. Phys. Chem. B 2015, 119, 12365−12380

Article

The Journal of Physical Chemistry B

(Scheme S6). All stationary points (reactants, products, and transition states) were located and characterized using an analogous version of the micro/macroiteration algorithm29 with convergence criterion of 0.01 kcal·mol−1·Å−1. The final Hessian was computed to check the existence of a single imaginary frequency (saddle point) or no imaginary frequency (minimum). Steepest descent path calculations were carried out for each step so as to verify that the transition states connect the respective reactant and product states. The different reaction coordinates used for the individual steps are listed together with the reaction and transition state energies in Tables S5, S6 in the associated content for the wild-type enzyme, and in Table S7 for the H151A mutant, respectively. The PES scans were carried out and analyzed with the pDynamo-1.8.0 program30 and the GTK visual interface.31 QM/MM Free Energy Calculations. Starting from each of the structures obtained for the different reaction coordinate values along every PES scan of each step in the mechanisms (Schemes S4−S6), we performed free energy calculations for all the steps following the same coordinates as those employed for the scans (described in Tables S5−S7). Those free energy profiles are the potentials of mean force (PMF) along the coordinate followed (termed from now on “reaction coordinate”). For each step which had been computed as a two-dimensional PES scan, we calculated the two-dimensional potentials of mean force (see Figures S3 and S4). An exception is the C1′−N1 bond dissociation which according to the PES scan follows a one-dimensional reaction coordinate (see Figures S1, S2, S10, and S11). We have also computed two-dimensional PMFs with the C1′−N1 distance and C1′−Ow1 distance as reaction coordinates, for both models with protonated His151, the direct mechanism (the proton is on His151), and the histidine mechanism where the proton has been transferred to the O4 atom of the thymine base. These PMFs, presented in Figures S1 and S2, clearly show a stepwise mechanism with C1′−N1 bond cleavage prior to nucleophilic attack. Since only one reaction coordinate is followed in each step, we recomputed the C1′−N1 bond dissociation step in all models in a onedimensional PMF so as to avoid putting unnecessary restraints on the system. The reaction coordinate of every step was binned into windows of 0.05 Å width (0.1 Å in the twodimensional PMFs). A harmonic biasing potential with a force constant of 180 kcal·mol−1·Å−2 was used to restrain the system around the desired value of the reaction coordinate corresponding to the center of the window.32 We computed 5 ps of relaxation and 50 ps of production (with a time step of 0.5 fs) for each window, using Langevin dynamics at 295 K to generate a canonical ensemble (NVT). The weighted histogram analysis method (WHAM)33 was used to reconstruct the complete distribution function from the free energy profiles that have been computed separately for each window. All PMF calculations were performed with the pDynamo-1.8.0 program30 and the GTK visual interface.31 We computed all the PMFs three times using different starting velocities, drawn from a Maxwell−Boltzmann distribution at the simulation temperature T = 295 K. The PMFs presented here are averages over these three replicas with error bars computed from the standard deviation from the mean. Analysis. All classical MD simulations have been analyzed for hydrogen bonds and important distances between active site residues. Hydrogen bonds between two atoms (donor and acceptor) have been analyzed based on the distance of these

isolated histidine and the histidine residue in the protein−DNA complex is then used to calculate the pKa shift from isolated histidine in water (pKa = 6.0) to histidine in the full protein− DNA complex environment. The dielectric constant of the solvent water was set to 80, that of the protein interior to 10. We used a grid spacing of 0.75 Å. Poisson−Boltzmann calculations have been performed using Charmm.26 QM/MM Model. From the classical MD simulations of wild-type and H151A TDG we obtained relaxed conformations serving as starting points for the subsequent QM/MM study.27 These starting points were chosen as a conformation in which those hydrogen bonds in the active site of the enzyme that are most probable (∼70% occupancy in the three replicas) are formed. In the QM/MM model all the residues likely to be involved in the chemical reaction of this enzyme, as indicated by mutation studies16 or based on the proximity of the thymine, and all water molecules close enough to make hydrogen bonds with those residues were treated quantum-mechanically, the remainder of the system was treated with molecular mechanics. Specifically, the QM part consists of the side chains of Asn140, His151 (Ala151), Asn191, and thymine, and two (five in the H151A mutant) water molecules, making a total of 66 QM atoms (68 QM atoms in the H151A mutant), see Figure 1. We also performed QM/MM calculations of the C1′−N1 bond dissociation step with Tyr152 additionally included in the QM part. However, comparison of free energy profiles of this reaction step, with and without Tyr152 treated quantum mechanically, showed only negligible differences. We therefore used the smaller QM part throughout the QM/MM calculations so as to enable more sampling. Since a recent computational study showed no impact of an adjacent phosphate on the reaction barriers,10 the phosphate group of the thymine nucleotide was not included in the QM part. To check whether the glycosidic bond cleavage reaction can be initiated by proton transfer from a H3O+ ion to the thymine base we have set up a model of the H151A mutant with a protonated water molecule, chosen based on the proximity to the thymine base. We have then computed the cleavage reaction in the H151A mutant starting from that state with a protonated water molecule (wat2) close to the thymine base. Five link atoms were placed between the Cα−Cβ atoms of the amino acids and between the C4′−C5′ atoms of the thymine nucleotide.27 The QM region was treated with the semiempirical Hamiltonian PDDGPM3.28 The rest of the system (i.e., the rest of the enzyme and DNA strand, the water molecules, and the counterions, a total of about 69 000 atoms) was treated classically with the Charmm27 and TIP3P force fields.23,24 All atoms farther than 12 Å from the active site were kept stationary, leaving about 3200 atoms free to move in this QM/MM study. Additional small molecule and active site model calculations that evaluate the semiempirical PDDGPM3 method against calculations with density functional theory are presented in the Supporting Information. All energies reported in this main text have been obtained with PDDGPM3 or PDDGPM3/Charmm. QM/MM Potential Energy Surface Scans. The mechanistic study of the cleavage reaction in the native TDG and the H151A mutant was commenced by exploring reaction mechanisms on the QM/MM potential energy surface. For each step in the reaction catalyzed by TDG, we performed relaxed potential energy surface (PES) scans, from which we were able to reconstruct the reaction mechanisms for wild-type TDG (Schemes S4 and S5), and for the H151A mutant 12368

DOI: 10.1021/acs.jpcb.5b05496 J. Phys. Chem. B 2015, 119, 12365−12380

Article

The Journal of Physical Chemistry B

Figure 2. (a) Sketch of the state before C1′−N1 bond dissociation in wild-type TDG highlighting the proton positions in the models with neutral and protonated His151. The positions of the protons in the different models are highlighted with the following colors: HSE is yellow, HSD is orange, protonated His151 in the direct mechanism is light blue, and protonated His151 (HSP) in the histidine mechanism is dark blue. The C1′−N1 bond is shown with a dashed line. (b) Free energy profile of the C1′−N1 bond dissociation for neutral histidine: HSE (yellow) and HSD (orange), and protonated histidine: direct mechanism (light blue) and histidine mechanism (HSP, dark blue).

(HSD, HSE, and HSP). In the following, the energetic and mechanistic results of the glycosidic bond cleavage reaction calculated for wild-type TDG are presented. After a comparison of water distribution in the active site of wild-type TDG and the H151A mutant, we present the results of the mechanistic study of the glycosidic bond cleavage reaction in the H151A mutant. Finally, the implications of structure and charge on the role of the catalytic residues Asn140, Asn191, Tyr152, and Thr197 are analyzed. H151 Can Be Protonated at Physiological pH, Depending on Conformation. In the crystal structure, reported in ref 16, the distance between His151 Nϵ and the backbone oxygen atom of Pro125 is short, suggesting a hydrogen bond between them. This hydrogen bond is observed with only low to moderate occupancy. The simulation with protonated His151 (HSP) shows 21% occupancy in one of the three replicas and 59% occupancy in the simulation of one of the neutral forms of His151 (HSE). The closest distance to Pro125 fluctuates around 3.0 ± 0.6 Å in the HSE, 4.2 ± 0.4 Å in the HSD simulation, and 4.3 ± 0.3 Å in the HSP simulation, respectively. Fluctuations in the His151-Asp178 distance are large for the HSE simulation, whereas this distance shows a small decrease in the simulations of the HSD model and an almost constant value in the HSP simulation (see Figures S25, S26, and S27 for the HSP, HSE, and HSD simulation, respectively). Hence, the protonation state of His151 cannot be determined from comparison of its conformations alone. The pKa values computed from the Poisson−Boltzmann calculations of several snapshots from the MD simulations of protonated His151 (HSP) and the two neutral forms of His151 (HSD and HSE) are listed in Tables S18−S20. His151 shows a

two atoms (less than 3 Å) and the acceptor-hydrogen-donor angle (more than 138°) using a python coded program. Hydrogen bond occupancies were computed for all the models (wild-type, H151A, N140A, and T197A mutants) and averaged over the three replicas of each. For wild-type TDG (HSP) and the H151A mutant, we also analyzed the water density around the O4 atom of the thymine base. These simulations together with the simulation of the N140A mutant, and the T197A mutant were further analyzed for the water density distribution around residue 140. These analyses have been performed using the radial distribution function tool from Gromacs34 and averaged over three replicas, each with different rounds of “resolvation”, except for the T197A mutant where no “resolvation” has been performed. Errors have been estimated as the standard deviation of these replicas. In order to quantify the effect of polarization and changes in electron density along the simulated reaction pathways, we have computed the Mulliken charges of some atoms of the active site residues involved in the reaction. Again, errors have been estimated as the standard deviation of the three replicas performed for each simulation. Graphics. For the graphical analysis, we used two different softwares: PyMOL,35 a package for rendering three-dimensional structures, and VMD,19,20 which also provides tools to animate the molecular dynamics trajectories and to analyze them.



RESULTS We report first the results of the classical MD simulations for wild-type TDG in three different protonation states of His151 12369

DOI: 10.1021/acs.jpcb.5b05496 J. Phys. Chem. B 2015, 119, 12365−12380

Article

The Journal of Physical Chemistry B

Figure 3. (a) Free energy profile of all steps in the direct mechanism of the glycosidic bond cleavage by the enzyme TDG computed with QM/MM methods. For all the intermediate structures the part of the active site that is relevant in the respective step is shown. Solid arrows indicate electron movement, and dotted arrows indicate a conformational transition. All free energies are shown relative to the reactant state. Reaction coordinates used are listed in Table S5. (b) Network of computed reaction pathways for glycosidic bond cleavage in wild-type TDG with protonated His151. The boxes correspond to reactant, product and intermediate states, labeled according to subfigures (a) and (c), respectively. Arrows between boxes indicate calculated transitions. Green italic numbers within the boxes are free energies of that state relative to the previous state, i.e. the starting point of the arrow. Blue numbers are free energies relative to the reactant state R. For states C, E, and F, the second, inner set of green and blue numbers refers to free energies relative to, or accumulated by, pathways coming from the respective counterpart in the upper (d) or lower (h) route. Brown italic numbers at the arrows represent free energy barriers with respect to the previous state; red numbers are free energy barriers relative to the reactant state. All energies are given in kcal/mol and are rounded for clarity. For higher precision numbers see Tables S8−S10. All pathways connecting R and P must pass a point of 12−14 kcal/mol free energy. The pathways in that energy range are highlighted with bold arrows. (c) Free energy profile of all steps in the histidine mechanism of the glycosidic bond cleavage by the enzyme TDG computed with QM/MM methods. For all the intermediate structures the part of the active site that is relevant in the respective step is shown. Solid arrows indicate electron movement and dotted arrows indicate a conformational transition. All free energies are shown relative to the reactant state. Reaction coordinates used are listed in Table S6.

rendering protonation of His151 possible at neutral pH. Thus, both Poisson−Boltzmann pKa calculations and molecular dynamics trajectories indicate that a significant population of His151 can be protonated at physiological pH. Glycosidic Bond Cleavage in Wild-Type TDG Is Stepwise and Acid-Catalyzed. Free energy profiles for the complete glycosidic bond cleavage mechanism have been computed in five variants: (i) wild-type TDG with neutral His151 (HSD), (ii) wild-type TDG with neutral His151

broad range of pKa values, depending on the MD simulation that has been used to sample the conformations. For those conformations, that have been sampled in the HSP simulation, the protonated form is clearly the more probable one, as manifested by pKa values from 7.9 up to 8.7. Conformations that have been sampled with neutral His151, however, are less likely to be protonated, and pKa values fluctuating around 7.1 (HSD) and 7.5 (HSE) are observed. In the case of the HSE simulation, however, values as high as 7.9 pKa are calculated, 12370

DOI: 10.1021/acs.jpcb.5b05496 J. Phys. Chem. B 2015, 119, 12365−12380

Article

The Journal of Physical Chemistry B

dissociation step has a free energy barrier of 11.0 ± 0.6 kcal/ mol. Like in the direct mechanism, the subsequent steps are the attack of the nucleophilic water molecule wat1 (Ch to Dh), rearrangement of the base (Dh to Eh), and proton transfer from water molecule wat1 to the N1 atom of the thymine base (Eh to Fh). The corresponding free energy barriers are 2.8 ± 1.1, 5.6 ± 0.7, and 5.9 ± 0.4 kcal/mol, respectively. The reaction is completed by back transfer of the proton from thymine oxygen atom O4 to His151, again via the water molecule wat2 (steps Fh-Gh-P). The free energies of the transition states in these proton transfer steps are 5.4 ± 0.5 and 3.2 ± 0.1 kcal/mol relative to their initial states, Fh and Gh, respectively. The overall reaction free energy for this histidine mechanism, calculated by adding the reaction free energies of the individual steps, is −10.4 ± 0.1 kcal/mol. This overall reaction free energy differs by 4 kcal/mol from the one calculated for the direct mechanism. In addition to the intrinsic error due to the summation of individually obtained relative free energy differences, this deviation can be attributed to errors that arise from limited sampling such as hysteresis. For example, the difference in relative free energies computed from forward and backward transitions in the same step is ∼1−2 kcal/mol (see Figures S6−S8). This is in the same range as the error estimated from the standard deviation of the different replicas in each free energy calculation. The C1′−N1 bond dissociation (Bh to Ch) has a transition free energy of 12.0 ± 1.5 kcal/mol, relative to the reactant state, which is the highest free energy on the histidine mechanism, albeit only marginally higher than the transition free energy for proton transfer from the nucleophilic water molecule to the N1 atom of the leaving group (step Eh to Fh: 10.9 ± 2.2 kcal/mol relative to reactant). The two mechanisms computed for the wild-type TDG with protonated His151 both consist of the same essential steps, C1′−N1 bond dissociation, water attack, base rearrangement, and proton transfer from the water nucleophile to the N1 atom of the leaving base. The differences between the two mechanisms are the initial proton transfer steps from His151 to the O4 atom of the base and the respective back transfer in the final steps of the histidine mechanism. In order to explore whether those proton transfer steps are conceivable also at different points along the two mechanisms, we computed free energy profiles for proton transfer between His151 and the O4 atom of the thymine base (via water molecule wat2) in both directions, i.e., to and from the O4 atom for different points of the mechanism: after C1′−N1 bond dissociation (between states Cd and Ch), after water attack (between states Dd and Dh), and after base rearrangement (states Ed and Eh). The respective free energy profiles are shown in the Supporting Information (Figures S6−S8). These additionally computed proton transfer steps can be understood as connections between the direct mechanism and the histidine mechanism, resulting in a number of different possible combinations. A network connecting all computed transitions within and between the two mechanisms is presented in Figure 3b). The upper route R−Cd−Dd−Ed−P corresponds to the direct (d) mechanism, the lower route R−Ah−Bh−Ch−Dh−Eh−Fh−Gh−P corresponds to the histidine (h) mechanism. Transitions between states C, D, and E connect from d to h and from h to d. The free energy barriers of the proton transfers for these “connections” are between ∼10 and ∼17 kcal/mol.38 Taking into account all possible pathways between reactant state R and product state P, all pathways via state Cd, i.e., a state

(HSE), (iii) wild-type TDG with protonated His151 (HSP) termed “histidine mechanism” due to the active role this residue, His151, plays in that mechanism, (iv) wild-type TDG with protonated His151 (HSP), termed “direct mechanism”, and (v) the H151A mutant. The two-dimensional PES scans and PMFs with the C1′−N1 distance and the C1′−Ow1 distance as reaction coordinates (Figures S1 and S2) clearly show a stepwise mechanism with C1′−N1 bond cleavage prior to nucleophilic attack for the neutral models HSD and HSE (see Figures S10 and S11) and for both the “direct mechanism” and the “histidine mechanism” in the model with protonated His151. An associative glycosidic bond cleavage mechanism is therefore unlikely in wild-type TDG. Figure 2 shows the free energy profile of the C1′−N1 bond cleavage step in the different models and mechanisms i−iv of wild-type TDG. In the case of the models with neutral His151 (HSE and HSD), the free energy barrier for the C1′−N1 bond dissociation is about 35 kcal/mol. In contrast to the lower free energy barriers for this step with protonated His151 (∼18 kcal/ mol in the direct mechanism and ∼12 kcal/mol for C1′−N1 bond cleavage in the histidine mechanism where the proton has been transferred to the thymine). The free energy barrier of ∼35 kcal/mol of the C1′−N1 bond cleavage step already renders the reaction unlikely to occur in TDG with neutral histidine. Moreover, all subsequent steps (water attack after C1′−N1 bond dissociation, base rearrangement, and proton transfer from the nucleophile to the N1 atom of the thymine base) show transition free energies that are even higher (41 and 44 kcal/mol for HSD, and 39 and 45 kcal/mol for HSE, respectively; see Figures S12 (HSD) and S13 (HSE)). The free energy profile of all steps in the direct mechanism in the model with protonated His151 is depicted in Figure 3a), together with the corresponding structures R−Cd−Dd−Ed−P (see also Scheme S4) . The first step of the direct mechanism is the dissociation of the glycosidic C1′−N1 bond, leading to a negatively charged base and an oxacarbenium ion at the sugar (intermediate Cd). This step is followed by the attack of the nucleophilic water molecule, wat1, (Cd to Dd) and a subsequent rearrangement of the thymine base relative to the sugar (Dd to Ed) thereby reducing the distance between the nitrogen atom N1 of the leaving group and one of the water hydrogen atoms (Hw1), facilitating the last step which is the proton transfer from water molecule wat1 to the N1 atom of the thymine base (Ed to P). The rate-determining step of this mechanism is the C1′−N1 bond cleavage in the first step with a free energy barrier of 17.9 ± 0.9 kcal/mol. The subsequent steps of water attack, rearrangement and proton transfer have comparably low free energy barriers of 3.4 ± 0.5, 3.2 ± 0.2, and 1.5 ± 0.3 kcal/ mol, respectively. The sum of the reaction free energies of the individual steps gives an overall reaction free energy of −6.4 ± 0.1 kcal/mol for the direct mechanism.36 The free energy profile of the alternative histidine mechanism and the corresponding structures R−Ah−Bh−Ch− Dh−Eh−Fh−Gh−P are shown in Figure 3c (see also Scheme S5). This mechanism involves a number of additional proton transfer steps compared to the direct mechanism (Figure 3a). The initial steps are a proton transfer from the protonated His151 to the O4 atom of the thymine base, via water molecule wat2 (R to Ah, and Ah to Bh). The free energies of those two steps are 7.2 ± 0.1 and 7.5 ± 0.3 kcal/mol, respectively.37 The proton transfer is followed by the departure of the protonated (and thus neutral) thymine base. This C1′−N1 bond 12371

DOI: 10.1021/acs.jpcb.5b05496 J. Phys. Chem. B 2015, 119, 12365−12380

Article

The Journal of Physical Chemistry B

Figure 4. Distribution of the water molecules close to thymine and His151 (Ala151) are depicted for: (a) native TDG and (b) H151A mutant TDG. (c) Radial distribution function of the water molecules around the O4 atom of the thymine is shown for the native system (blue line) and for the mutant system (green line).

wild-type (Figure 4a) and H151A mutant (Figure 4b) of TDG, respectively. In the H151A mutant, bulk water can easily enter the active site and is free to move close to the nucleotide. Figure S28 highlights the positions visited by one water molecule in the course of the simulations. The radial distribution functions of water oxygen atoms around the O4 atom of thymine (Figure 4c) show a peak at ∼2.8 Å in both the native system and the H151A mutant, indicating a high probability for a water molecule to be hydrogen bonded to the O4 atom of thymine. The smaller peak at a larger distance in the radial distribution function near 4.0 Å can be interpreted as arising from a second solvation shell that is less well-ordered in the H151A mutant than in the wild-type indicating higher mobility of the water molecules in the active site. As shown in Figure S29, the number of water molecules in the active site of the mutant, measured as those between Cβ atom of residue 151 and the O4 atom of the thymine base, is significantly larger in the mutant than in wild type TDG. It is therefore conceivable that a H3O+ ion enters the active site or a proton is shuttled along several water molecules resulting in a protonated thymine. Glycosidic Bond Cleavage in H151A TDG Has a Lower Free Energy Barrier than in Wild Type. The mechanism and the corresponding free energy profiles for the glycosidic bond cleavage reaction in the H151A mutant of TDG are shown in Figure 5. In the H151A mutant of TDG, the reaction is initiated by a proton transfer from the hydronium ion to the O4 atom of the thymine base (Rm to Bm). The protonated thymine intermediate has a free energy that is 9 kcal/mol lower than the reactant state. The following C1′−N1 bond dissociation (Bm to Cm) has a free energy barrier of 5.5 ± 0.2 kcal/mol and leads to a dissociative intermediate (Cm). The subsequent steps are water attack (Cm to Dm), followed by

with cleaved C1′−N1 bond and unprotonated thymine, must cross a transition state of 16.7 ± 3.2 kcal/mol (via connection Ch to Cd) or 17.9 ± 0.9 kcal/mol (the R to Cd step in the direct mechanism), the highest free energies of the entire network of pathways. Cleavage of the C1′−N1 bond in a protonated thymine has a lower free energy barrier (12.0 ± 1.5 kcal/mol) than in neutral thymine. All other transitions in the network require less or similar free energies. All pathways leading from state Ch to product (avoiding state Cd) are feasible with a free energy of ∼14 kcal/mol. These pathways are highlighted with thick arrows in Figure 3b). In addition to the histidine (h) mechanism, possible pathways are R−Ah−Bh−Ch−Dh−Dd− Ed−P, corresponding to the following steps: first the thymine base is protonated, then the C1′−N1 bond cleavage and water attack take place with protonated thymine. The deprotonation of the departed thymine base (Dh−Dd) is the crossover to the direct mechanism. Its transition state has a free energy of 13.0 ± 0.6 kcal/mol or 9.2 ± 1.7 kcal/mol, depending on the direction in which the transition has been calculated (d to h or h to d, respectively). Base rearrangement and final proton transfer between the nucleophile and N1 atom (Ed−P) take then place with the base remaining unprotonated, passing only points with again similar free energies (12.9 ± 1.2 and 10.8 ± 1.5 kcal/mol, respectively) . Alternatively, base rearrangement can take place with a protonated base (Dh−Eh), followed then by base deprotonation (Eh−Ed transition) and the final proton transfer between wat1 and the N1 atom of the thymine (Ed−P). The highest free energy along this pathway is in the C1′−N1 bond scission step (12.0 kcal/mol). Larger Number of Water Molecules in H151A TDG Renders Entrance of Catalytic Proton as H3O+ Feasible. Figure 4 compares the distribution of water molecules as observed in classical molecular dynamics simulations of the 12372

DOI: 10.1021/acs.jpcb.5b05496 J. Phys. Chem. B 2015, 119, 12365−12380

Article

The Journal of Physical Chemistry B

free energy of the C1′−N1 bond dissociation of protonated thymine is only slightly lower in the H151A mutant (−0.6 kcal/ mol) than in the wild-type histidine mechanism (+2.5 kcal/ mol). In contrast, the free energy barrier of the C1′−N1 bond dissociation is significantly lowered (5.5 kcal/mol) in the mutant mechanism compared to the wild-type histidine mechanism (11.0 kcal/mol). Catalytic Residue Asn140 Positions Nucleophilic Water Molecule. Active site residue Asn140 has been shown by mutagenesis experiments17 to be crucial for catalysis. According to our simulations, Asn140 positions the nucleophilic water molecule that attacks the oxacarbenium ion (see Figure 1). Figure 6 highlights this positioning role of Asn140, showing snapshots of state C in the direct mechanism, Cd, the histidine mechanism, Ch, and in the mutant mechanism, Cm. In Figure 7 the region around Asn140 is compared for the wildtype TDG (Figure 7a), the H151A mutant (Figure 7b), and the N140A mutant (Figure 7c) TDG in the final frames of 100 ns MD simulation of each system, clearly showing at least one water molecule close to Asn140 in the wild-type and the H151A mutant. The radial distribution function of water oxygen atoms around the carbonyl oxygen atom in Asn140 shows a prominent peak at ∼2.8 Å indicating a well-preserved hydrogen bond between Asn140 and a water molecule that is a likely candidate for the nucleophile. In contrast, there is no evidence of a plausible nucleophilic water molecule positioned close enough to attack the thymine sugar in the N140A mutant system. Instead, the space that is filled by water molecules in the presence of Asn140 is partially occupied by Ile139 that is closer to the thymine sugar in the N140A mutant than in wildtype TDG (see Figure 7c). Catalytic Residue Asn191 Helps To Position and (De)Polarize the Thymine Base. Analysis of the classical molecular dynamics simulations with respect to the hydrogen bonds in the active site show a rather strong hydrogen bond between Asn191 and the thymine base (more than 90% occupancy in the classical MD simulations) in the wild-type TDG and both mutants, H151A and N140A (see Table S16), confirming the role of this residue in positioning the thymine base prior to the cleavage reaction. This initial hydrogen bond is not maintained through the entire reaction, but, the OD1 atom of Asn191 still remains rather close to the nucleotide (see Figures S16, S17, and S18). Analysis of the changes in electron density, expressed as Mulliken charges, in the course of the reaction, reveals that Asn191 is involved in the reaction by accepting and donating charge (Figures S20−S22). This effect is most pronounced in the histidine mechanism, in which the proton transfer to the thymine base leads to a significant change

Figure 5. Free energy profile of the H151A mechanism of the glycosidic bond cleavage by the enzyme TDG, computed with QM/ MM methods. For all the intermediate structures the part of the active site that is relevant in the respective step is shown. Solid arrows indicate electron movement, and dotted arrows indicate a conformational transition. All free energies are shown relative to state Bm. Reaction coordinates used are listed in Table S10.

rearrangement of the thymine base (Dm to Em), and proton transfer from nucleophilic water wat1 to the N1 atom of the base (Em to Fm), as in the mechanisms computed for wild-type TDG. The final step is a back transfer of the proton from the O4 atom of thymine to the water molecule wat2. This last step is the transition with the highest free energy barrier (12.0 ± 0.1 kcal/mol) relative to its initial state (Fm). The significantly lower free energies of the states with protonated thymine Bm and Fm compared to their corresponding states with neutral thymine and protonated water wat1, Rm and Pm, respectively, suggest that states Bm and Fm can be regarded as the actual reactant and product state of the glycosidic bond cleavage in the H151A mutant, respectively. The free energy profiles in Figure 5 are therefore shown relative to state Bm.39 Like in the histidine mechanism computed for the wild-type enzyme, the attack of the water nucleophile leads to a state (Dm) with lower free energy than the previous C1′−N1 bond cleaved intermediate (state C). The subsequent rearrangement steps (Dm to Em), and the proton transfer from the nucleophilic water molecule wat1 to the N1 atom of the base (Em to Fm) show higher transition free energies in the mutant than in the histidine mechanism (see Figure 3b), Figure 5, and Tables S9 and S10 for free energies in the wild-type histidine mechanism and the mutant mechanism, respectively). This is likely a consequence of the larger C1′−N1 distance in the Dm state of the mutant (3.5 Å) than in the Dh state of the wild-type histidine mechanism (2.8 Å) (see Figure S19). The reaction

Figure 6. Close-up view of the conformation prior to the nucleophilic water attack in (a) the direct mechanism, Cd, (b) the histidine mechanism, Ch, and (c) the mutant mechanism, Cm (cf. Figures 3a, 3c, and 5). The nulceophilic water molecule, wat1, that attacks the sugar ring is positioned and oriented by residue N140. 12373

DOI: 10.1021/acs.jpcb.5b05496 J. Phys. Chem. B 2015, 119, 12365−12380

Article

The Journal of Physical Chemistry B

Figure 7. Close-up view of the (a) wild-type, (b) H151A mutant, and (c) N140A mutant TDG after 100 ns MD. The clouds show the spatial distribution of water molecules close to C1′ and N140. For clarity only the distribution in the last 10 ns of the simulation is shown: in the N140A mutant water is absent; instead Ile139 occupies the space filled by water in the other two models, indicating that a nucleophilic water-attack at C1′ is unlikely. (d) Radial distribution function of the water molecule oxygen atoms around the carbonyl oxygen atom of N140 in the native system (blue), the H151A mutant (green), and around the Cβ atom of A140 in the N140A mutant (magenta).

Table 1. Computed Free Energies of the C1′−N1 Bond Scission Step and the Highest Free Energy Point along the Different Mechanisms Together with the Corresponding Reaction Rates (According to Transition State Theory) in Comparison to Experimentally Obtained16 Reaction Rates (in Italics) and Corresponding Free Energy Barriers for Wild-Type and H151A TDG mechanism HSD HSE HSP direct HSP histidine H151A mutant exp. WT exp. H151A

barrier C1′−N1 bond scission (kcal/mol) 34.2 35.3 17.9 12.0 5.5

rate (s−1) 2.82 4.32 0.34 79.1 51.7

highest free energy (kcal/mol)

−13

× 10 × 10−14 × 102 × 107

41.4 44.9 17.9 12.0 9.8 20.7 19.3

rate (s−1) 1.31 3.34 0.34 79.1 33.7 2.74 3.30

× 10−18 × 10−21 × × × ×

102 104 10−3 10−2

the main chain of Thr197 and the side chain of Asn140 in either of the wild-type or T197A mutant simulations (cf. Table S17). Moreover, the water distribution around Asn140, analyzed as radial distribution function (cf. Figure S30) is similar in both systems with a small reduction of water molecules observed in hydrogen bond distance to Asn140 in the T197A simulation. The position of Asn140 with respect to the thymine base, however, shows no significant difference between the simulations of the wild-type and the T197A mutant. These findings suggest the role of Thr197 in positioning Asn140 and hence the nucleophile to be only minor, at least after formation of the complex with the flippedout thymine base accommodated in the active site.

of charge density in the active site and both the base itself as well as Asn191, show a significant change in the Mulliken charges. Active Site Residue Tyr152 Helps To Position Thymine Base. According to our classical molecular dynamics simulations, Tyr152 is in a stacked conformation to the thymine base throughout the simulation time with an average distance of 4.5 Å and an average angle between planes of 28 degrees (see Figure S24). This relative orientation positions the thymine base in a flipped-out conformation that is further stabilized by interaction to Asn191. Our QM/MM calculations with and without Tyr152 in the QM part (cf. Figures S14 and S15) show that this residue does not contribute to the reaction in other than a sterical manner. Thr197 Does Not Play an Active Role in Prereaction Complex. Thr197 and Asn140 form hydrogen bonds between their side chains about half of the simulation time in two out of three replicas of the wild-type of TDG. Only a fraction of the simulation time (11−28%) a hydrogen bond is formed between



DISCUSSION

We have computed reaction pathways of the catalytic mechanism of glycosidic bond cleavage in wild-type TDG for two models with neutral His151 (HSD and HSE, respectively) and two mechanisms in a model with protonated His151: a 12374

DOI: 10.1021/acs.jpcb.5b05496 J. Phys. Chem. B 2015, 119, 12365−12380

Article

The Journal of Physical Chemistry B

between pH 6.5 and 7.5 (and between 6 and 9 for G:U). The G:T cleavage rate increases at pH values below 6.5 and decreases above pH 7.5. These data are consistent with either acid catalysis (by e.g. a protonated Histidine) which does not play a role for pH > 6.5 and His has a pKa < 6; or the acid (the protonated His151) is present throughout the pH range 6.5 to 7.5 and deprotonated at pH > 7.5. The calculations show a pKa > 7 for all snapshots from all three simulations, except for two snapshots whose pKa is also larger than 6. It is therefore more likely that pH independence of the reaction rate below pH = 7.5 is due to a protonated His151 with high pKa value that can act as a catalytic acid. Most Likely Mechanism in a Network of Reaction Pathways for Wild-Type TDG. For the model with protonated His151 we have computed a whole network of reaction pathways, consisting of the direct mechanism and the histidine mechanism connected to each other at various steps along the mechanism. Aside from the proton transfer to and from the thymine O4 atom, the direct and the histidine mechanism have the same sequence of steps: (1) C1′−N1 bond dissociation, (2) water attack, (3) base rearrangement, and (4) proton transfer from nucleophile to thymine base. The only difference is the fact that in the histidine mechanism, all steps are carried out with a protonated thymine base, whereas in the direct mechanism the proton remains on the His151. A comparison of the steps that are similar in both mechanisms shows that the activation free energy for the C1′−N1 bond dissociation step is lower with the protonated thymine, i.e., in the histidine mechanism than in the direct mechanism. On the other hand, the free energy barrier for the proton transfer from the nucleophile to the N1 atom of the (already neutral) thymine base is higher in the histidine mechanism than in the direct mechanism in which the cleaved base is negatively charged. In fact, in the histidine mechanism this proton transfer from the attacking nucleophilic water molecule wat1 to the N1 atom of the base has a free energy that is comparable to the one for the C1′−N1 bond dissociation step. It is therefore also possible that the system circumvents this second high free energy step in the histidine mechanism by reverting to the direct mechanism and carrying out the last steps with an unprotonated thymine base. This crossover to the other mechanism involves a proton transfer from the thymine base back to the His151 (via water molecule wat2). Taking together the two mechanisms computed for the wildtype TDG with protonated His151 and all possible branching points, the point with the highest free energy (relative to the reactant state) of that network of reaction steps is 17.9 ± 0.9 kcal/mol. This value is very similar to the free energy of 20.7 kcal/mol calculated from the experimentally measured rate (kobs at 295 K is 0.16 ± 0.03 min−1).16 The most probable pathway is the one with the lowest free energy in its ratedetermining step. On the network shown in Figure 3b), this pathway involves C1′−N1 bond dissociation with protonated thymine (R−Ah−Bh−Ch). The proton can then be shuffled back from O4 to His151 after water attack, after the base rearrangement, or after the final proton transfer between nucleophile and the N1 atom of the base (as in the histidine mechanism). Note that the highest free energies along all of these pathways are ∼12−14 kcal/mol, rendering all of these variants equally probable. However, all of these options are feasible only in the presence of an extra proton, i.e., a protonated His151 and not with neutral histidine.

direct one which we term direct mechanism and one assisted by His151 (termed histidine mechanism). All computed pathways are stepwise, dissociative mechanisms. Computed free energies of the rate-determining steps of all mechanisms are summarized in Table 1, and the experimentally determined rates16 in wildtype and H151A mutant TDG are given for comparison. The free energy profiles computed for the reaction pathways in the models with neutral His151 show significantly larger barriers for the C1′−N1 bond scission step than in the two pathways computed for the model with protonated His151. Already the presence of a proton on His151, even without explicit transfer to the thymine base (direct mechanism), lowers the activation free energy for the C1′−N1 bond dissociation. Proton transfer to the O4 atom of the thymine base prior to C1′−N1 bond dissociation leads to a further lowering of the activation free energy. The subsequent steps in the pathways of the HSD and HSE model show even higher free energies than the C1′−N1 bond scission step, rendering glycosidic bond cleavage in the presence of a protonated His151 more likely. The computed free energy barriers suggest an enhancement of C1′−N1 bond dissociation as a consequence of base protonation. This effect is more pronounced in the H151A mutant than in wild-type TDG. Also the highest free energy point along the mutant mechanism is lower than in the histidine mechanism, in agreement with the experimentally determined higher reaction rates for the mutant. The computed free energies, however, overestimate the reaction rate compared to the experimental values. Taking into account that in the computed mechanisms an additional proton is required in the active site (protonation of His151 in wild-type TDG or protonation of an active site water molecule in the H151A mutant), it is likely that only the fraction of enzyme-DNA complexes in which this proton is present react at a significant rate. The experimentally observed rate can thus be explained as an apparent rate, lowered by the amount of nonreactive complexes. Protonation State of His151 and pH Dependence of Reaction Rate. There are several indicators for His151 in TDG to be protonated at least in a fraction of the ensemble. Our classical MD simulations show a higher similarity to the crystal structure in the case of protonated His151 than for the neutral models (HSD and HSE). Moreover, the pK a calculations show that, depending on the conformation, His151 reaches pKa values of 8 and above. This is an unusually high pKa for histidine, however, there are a number of examples in the literature that report high pKa values for histidine residues. Based on histidine hydrogen−deuterium exchange mass spectrometry exchange rates and pKa values have been measured for different forms of dihydrofolatereductase (DHFR), H124 “is surrounded by the backbone carbonyl oxygen atoms of Asp11, Arg12, Gly121, and Asp122 and has a pKa = 7.6. And H141 is close to Glu139 and has a pKa = 8.04”.40 Older NMR work on hen egg white and human lysozyme shows that “the pK, of His residues are largely affected by surrounding ionized and polar groups”.41 The observed pKa values range from 5.32 to 7.49. The broad range of histidine pKa values, between 4.6 and 9.2, has further been demonstrated in ref 42. The rather strong conformation dependence of the calculated pKa values, however, suggests, that His151 is protonated only in the fraction of protonation-supporting conformations. Reference 43 reports the pH dependence of enzymatic activity of TDG on G:T to be rather constant 12375

DOI: 10.1021/acs.jpcb.5b05496 J. Phys. Chem. B 2015, 119, 12365−12380

Article

The Journal of Physical Chemistry B

ion has entered the active site, proton transfer to the O4 atom thymine base is very likely to occur. This entrance of hydronium ions into the active site and thus the reaction rate in the H151A mutant of TDG are anticipated to be largely pHdependent. To our knowledge, no corresponding experiments, probing the pH-dependence of the glycosidic bond cleavage reaction rate in H151A TDG have been reported, yet. Effect of Protonation on Leaving Group Ability. Taking together the results of the mechanistic studies on wild-type and H151A mutant TDG, the major catalytic effect is the presence of a proton in the active site. Both the proximity of a proton (on His151) to the O4 atom of the thymine base and even more the actual transfer to O4 (from His151 or a hydronium ion in the case of the H151A mutant), significantly reduce the barrier for the rate-determining C1′−N1 bond scission step. Kinetic experiments on the TDG activity that were conducted with different nucleobases indicate a strong dependence on the leaving group pKaN1 (the microscopic ionization constant for the N1 site).44 These findings suggest that leaving group departure, upon which the N1 site accumulates negative charge, is rate determining. This is in agreement with our simulations that find the C1′−N1 bond scission as the step with the highest transition free energy, albeit comparable to the transition free energies of subsequent steps. The importance of leaving group ability of the base, and hence the ability to carry extra negative charge as a consequence of C1′−N1 bond scission, has been shown in previous studies of nonenzymatic pyrimidine hydrolysis reactions.45,46 The importance of leaving group protonation has also been discussed for other base-excision enzymes. For example, in a calcium-dependent purine-specific nucleoside hydrolase, N7 protonation strongly facilitates cleavage of the guanidine nucleoside.47,48 The N7 position of guanine is the most nucleophilic site among the DNA bases and 2′-deoxyguanine has been found to be particularly susceptible to acid-catalyzed depurination.49 According to calculations with density functional theory protonation at N7 of the guanine nucleobase greatly catalyzes the hydrolysis of the N-glycosydic bond. Moreover, the reaction also occurs via a stepwise mechanism with a discrete oxocarbenium ion intermediate, similar to the mechanism we observe in TDG.50,51 According to our simulations, the amount of negative charge on N1 accumulated upon bond scission is the same with anionic and neutral leaving group in the direct and histidine mechanism respectively (see Supporting Information, Figures S20 and S21). The N1 pKa of the leaving group therefore does not determine whether the leaving group is anionic or neutral. However, a stabilization of the negative charge accumulating on the departing base in the transition state, e.g., by a nearby proton, transferred or close to the O4 atom, likely facilitates the glycosidic bond cleavage. In ref 44 the fact that a proton transfer is feasible to thymine (and uracil), whereas cytosine lacks the accepting O4 atom, is suggested as an explanation for the “enhanced (20 000-fold) difference in hTDG activity for G:T versus G:C pairs, as compared to the 150-fold difference in intrinsic reactivity between dT and dC”. The reactivity of the N-glycosidic bond has also been discussed to account for TDG’s activity against oxidized forms of 5methyl-cytosine (5hydroxymethyl-cytosine, 5hmC, 5formyl-cytosine, 5fC, and 5carboxyl-cytosine, 5caC). TDG excises 5fC and 5caC, whereas it is not active against 5hmC (and

Computed overall activation free energies that do not take into account that only a fraction of His151 is indeed protonated; i.e., the probability for the reaction to take place (computed as free energies) is a conditional probability given that His151 is protonated, are likely to underestimate the activation free energy obtained from experiments (cf. Table 1). As the protonated form is much more reactive, even a small fraction of HSP will contribute much more to the reactive flux than the neutral histidine tautomers. Mechanism in the H151A Mutant and Role of His151. At first sight this crucial role of His151 seems to contradict the experimental result of an increased activity of TDG upon H151A mutation. The simulations of the mechanism in the mutant enzyme, however, reveal that mutation of His151 does not conflict with an increased activity. Once the protonated thymine has been formed, the mechanism in the H151A mutant is almost equivalent to the histidine mechanism in the wild-type. The energetics of the individual steps, however, are different due to the changed protein environment. The overall activation free energy is on the same order of magnitude as the one obtained for the native enzyme, although the ratedetermining step is the proton transfer from the nucleophile to the thymine base (see Figure 5). The classical MD simulations show that, in the mutant system, the water (and the protonated water) molecules are clearly more free to move in the mutant system than in the wild-type. Moreover, entrance of bulk water (and hydronium ions) is less hindered by the smaller alanine side chain than in the wild-type where the (neutral) histidine is likely to “catch” the entering protons. The higher mobility of water molecules and of the excised base also renders other pathways for the reprotonation steps possible. For example, the proton from water molecule wat1 could be passed on to another water molecule after the base has left and moved farther apart. Various ways of proton transfer between the O4 or N1 atom of the base and surrounding water molecules are likewise conceivable along the release of the excised base product. The smaller size of A151 compared to His151 facilitates water movement and probably also base release, possibly resulting in an overall rate acceleration compared to wild-type TDG. In the kinetic experiments reported in ref 16, the detailed rates of individual steps are not resolved. According to Morgan et al.,6 the measured rate corresponds to the overall rate for all steps from DNA binding over base flip up to a state with cleaved C1′−N1 bond, but not base release. It is likely that nucleophilic attack has also taken place as part of this chemical step. Any subsequent reprotonation steps are, however, not contained or not resolved in the reported rate.16 It is interesting to note that, according to our network of simulated pathways in the wild-type enzyme, there is no clear preference for the exact order of those reprotonation steps. A comparison of the free energy computed for the ratedetermining step up to water attack in the H151A mutant with the histidine mechanism shows that the free energy of C1′−N1 bond scission is significantly lower in the H151A mutant than in the wild-type, a result in agreement with the experimentally measured rate acceleration upon H151A mutation. Again, the reaction probability, computed as transition free energy, does not take into account that the reaction occurs only given the thymine base is protonated. The computed reaction rate is thus overestimated. In a nonenzymatic environment, thymine is unlikely to accept a proton at pH > 2. However, our calculations show that in the enzyme, once the hydronium 12376

DOI: 10.1021/acs.jpcb.5b05496 J. Phys. Chem. B 2015, 119, 12365−12380

Article

The Journal of Physical Chemistry B

reaction, however, Asn191 remains close to the thymine base. This proximity enables Asn191 to transfer electron density to and from the nucleotide and thereby facilitating the charge transfer associated with protonation of the base. Tyr152 does not directly participate in the reaction but rather supports the position of the thymine base. Comparison to Other Glycosylases. The conserved protein motifs in the members of the UDG superfamily and sequences in the MUG family 2 suggest these enzymes to interact similarly with their admittedly large range of substrates, and to possibly also employ similar cleavage mechanisms.5,54 Asn140, for example is highly conserved in members of the MUG family, and crystal structures of MUG55 and UNG17,56 show a putative water nucleophile close to that residue. Moreover, UNG forms hydrogen bonds between Asn204 and its substrate uracil57,58 similar to the interactions between Asn191 and the thymine base in TDG. The K68N mutant of MUG shows enhanced UDG activity,59 possibly by similar interactions with the base. The mechanistically best-studied enzyme of the superfamily UDG has been reported to follow a stepwise mechanism.4 According to our simulations, the mechanism followed by TDG is also stepwise. In contrast to the many similarities between TDG and UDG, the two enzymes have different active site architectures. TDG lacks the aspartate residue that can serve as a nucleophile activator in UDG. Moreover, in UDG there is no histidine at the position of His151 in TDG, and UDG has been reported not to employ acid catalysis.14,60 However, the catalytic His148 in UDG closely interacts with the O2 atom of the uracil base61 and is suggested to thereby also enhance leaving group departure. His151 is not conserved in the MUG family. According to our simulations, the role of His151 is to provide a proton, and this role can also be played by a hydronium ion as the simulations of the H151A mutant show. Therefore, it is conceivable that the presence of His151 in TDG additionally provides specificity (for thymine) rather than rate acceleration only. The ability to discriminate thymine from mispairs over well-paired thymine is important to avoid aberrant base excision, a feature presumably less important for glycosylases processing other substrates. Mismatched thymine is a product of methyl-cytosine deamination, and therefore removal of mismatched thymine is important only in organisms with a significant degree of cytosine methylation, e.g., in mammals but not in yeast. Accordingly, thymine is not processed significantly in enzymes of bacteria or yeast.62,63 Human TDG, in contrast, does recognize and process thymine with evolution allowing it to develop a high level of specificity.

5methyl-cytosine, 5mC). In particular, the electron-withdrawing C5-substituent in 5fC has been attributed to increase the leaving group ability. For 5caC only the form with a protonated carboxyl group has a similar effect, rendering leaving group protonation catalytically important. Indeed, the cleavage rate for 5caC is enhanced at low pH.52 The probability for a proton to be transferred to the carboxyl group of 5caC or to the formyl oxygen atom of 5fC is likely to be higher than a proton transfer to the hydroxyl group of 5hmC. Hence, a glycosidic bond cleavage mechanism that requires proton transfer to the base, like the one reported here for thymine excision, may be a possible explanation of TDG’s activity against modified cytosine bases. Exploring the exact mechanism of glycosidic bond cleavage in oxidized methyl-cytosines clearly remains a subject of future studies. Role of Active Site Residues. According to our QM/MM calculations and the classical MD studies, the role of Asn140 is to properly position the nucleophilic water molecule prior to attack (see Figures 6 and 7) in agreement with earlier suggestions.17 Crystal structures of the pre- and postreactive complex of the N140A mutant with DNA carrying an oxidized methyl-cytosine (5caC)53 indeed do not reveal a putative water nucleophile. In the postreactive complex, this absence can be explained by the nucleophile having already attacked the cleaved sugar. In the prereactive complex, the lack of a putative nucleophile could also be due to the resolution of the crystal structure (3.0 Å). It is interesting to note in this context, that in the crystal structure of wild-type TDG, which has a similar resolution (2.97 Å), the putative nucleophilic water molecule is the only water molecule resolved. This could be regarded as an indicator of its tight positioning. The clear need for a properly positioned nucleophile and the lack of a likely nucleophile water molecule in the N140A mutant explain the experimental result of a 24300-fold decrease in the reaction rate when N140 is mutated to A140.17 Without a properly positioned water molecule, nucleophilic attack is likely to be the rate-determining step as opposed to the C1′−N1 bond scission that is ratedetermining in the presence of a well-poised nucleophilic water molecule. Thr197 has been suggested to in turn help positioning Asn140. Our MD simulations of wild-type and T197A mutant TDG show only moderate interaction between these two residues as manifested by the rather low population of hydrogen bonds between them and a similar position of Asn140 in both wild-type and mutant simulations, whereas the crystal structure of the wild-type TDG complex with the substrate analogue fluoro-uracil indicates a hydrogen bond between the hydroxyl group of Thr197 and Asn140, which maybe attributed to less flexibility in the crystal than in the solution. The simulations hence cannot confirm the suggested role of Thr197, at least after formation of the complex with the thymine base flipped into the active site and properly positioned for base excision. That is, Thr197 may help positioning Asn140 during the complex formation and base flip but appears not to contribute to the stabilization of the active site conformation. The catalytic role of residue Asn191 that has also been reported to be essential16 is mainly structural, helping to keep and probably also to position the flipped-out base in the active site prior to attack as can be seen from the rather strong hydrogen bonds formed between Asn191 and the thymine base (see Figure 1 and Table S16). This hydrogen bond between Asn191 and the thymine base is not maintained throughout the



CONCLUSION We have investigated the glycosidic bond cleavage reaction in human thymine DNA glycosylase by QM/MM calculations of small molecule models and QM/MM calculations of the enzymatic reaction supported by classical MD simulations of the TDG−DNA complex in wild type and some important mutants. The results of our calculations strongly indicate a dissociative mechanism in which the C1′−N1 bond is cleaved prior to the attack of a water nucleophile to be favorable. Our QM/MM calculations of the enzymatic mechanism show that C1′−N1 bond dissociation is significantly facilitated by protonation of the thymine base. Moreover, even in the case of the direct mechanism where the proton (at His151) is close to, but not explicitly transferred to the thymine base, the activation free energy for C1′−N1 bond dissociation is 12377

DOI: 10.1021/acs.jpcb.5b05496 J. Phys. Chem. B 2015, 119, 12365−12380

Article

The Journal of Physical Chemistry B significantly lowered, and thus the corresponding transition rate is significantly increased, compared to the enzymatic models with neutral His151. The probability for His151 to be protonated is conformation dependent but feasible also at neutral or slightly elevated pH. Computation of a whole network of reaction pathways that contains all steps along the mechanism with and without protonated thymine base, reveals proton shuttling between His151 and the O4 atom of thymine to be feasible at several stages along the mechanism. Exploiting this fact, in the energetically most feasible pathway, thymine is protonated for C1′−N1 bond scission, providing a neutral leaving group and a more positive electrophile at the sugar. The (catalytic) proton is equally likely to be shuttled back to His151 before and after water attack, base rearrangement or the final proton transfer between the nucleophile and the N1 atom of the base, giving rise to several possible routes for completing the catalytic reaction (and even more so for subsequent base release). Simulations of the H151A mutant confirm that the positive charge in the form of an extra proton is essential for catalysis rather than the His151 residue itself. In the wild-type TDG, His151 carries the catalytic proton whereas in the H151A mutant the proton is provided by a hydronium ion. Bulk water entering the active site of the H151A mutant is conceivable to allow for a direct proton transfer to the thymine base without a shuttle mechanism involving His151. The activation free energy computed for glycosidic bond cleavage in the H151A mutant is significantly lower, corresponding to a significantly higher transition rate than in the wild-type enzyme, explaining the experimentally observed increased activity of the H151A mutant. Our mechanistic studies further reveal the role of other active site residues that have been shown to be catalytically important. Asn191 is strongly hydrogen-bonded to the thymine base in the reactant state indicating that it plays a role in positioning the base. A similar positioning role can be attributed to Tyr152. The nucleophile is a water molecule positioned close to the thymine sugar by Asn140. This important role of residue Asn140 is further supported by our MD simulations of N140A TDG where, in contrast to the simulations of the wild-type, no water molecule is observed in a position that renders it likely to be the nucleophile. Taken together, by exploring a variety of reaction pathways, our simulations afforded us a more precise description of the most likely mechanism by which thymine DNA glycosylase achieves base excision of a mismatched thymine together with alternative routes that are energetically feasible. We could identify the role of a number of essential active site residues and reveal the catalytic effect, in addition to positioning the nucleotide and the nucleophile, to be mainly due to an additional proton that facilitates the leaving group departure.





tests for convergence and for the size of the QM region, Mulliken charges along the reaction, active site geometry information and pKa calculations, and water distribution analyses (PDF)

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Notes

The authors declare no competing financial interest.



ACKNOWLEDGMENTS This work has been funded by the German Science foundation (DFG) through grant IM141/1-1. R.C. acknowledges financial support from the Spanish MINECO (CTQ2012-33324). We would like to acknowledge the pDynamo group for the useful suggestions. Furthermore, we express our gratitude for the support given by GTK developers, especially Jose Fernando R. Bachega. We thank Jens Dreger for his constant help and support with the computer cluster of the department. We are grateful to Nadia Elghobashi-Meinhardt and Guillermo PérezHernández for carefully reading the manuscript and to them and to Frank Noé for helpful discussions. Computational resources provided by the North-German Supercomputing Alliance (HLRN) and the Zedat cluster “Soroban” of the Freie Universität Berlin are gratefully acknowledged.



REFERENCES

(1) Jiricny, J. DNA repair: An APE that proofreads. Nature 2002, 415, 593−594. (2) Brooks, S.; Adhikary, S.; Rubinson, E.; Eichman, B. Recent advances in the structural mechanisms of DNA glycosylases. Biochim. Biophys. Acta, Proteins Proteomics 2013, 1834, 247−271. (3) Maiti, A.; Morgan, M. T.; Pozharski, E.; Drohat, A. Crystal structure of human thymine DNA glycosylase bound to DNA elucidates sequence-specific mismatch recognition. Proc. Natl. Acad. Sci. U. S. A. 2008, 105, 8890−8895. (4) Stivers, J.; Jiang, Y. A mechanistic perspective on the chemistry of DNA repair glycosylases. Chem. Rev. 2003, 103, 2729−2759. (5) Cortazar, D.; Kunz, C.; Saito, Y.; Steinacher, R.; Schar, P. The enigmatic thymine DNA glycosylase. DNA Repair 2007, 6, 489−504. (6) Morgan, M.; Bennett, M.; Drohat, A. Excision of 5-halogenated uracils by human thymine DNA glycosylase: robust activity for DNA contexts other than CpG. J. Biol. Chem. 2007, 282, 27578−27586. (7) Berti, P.; McCann, J. Toward a detailed understanding of base excision repair enzymes: Transition state and mechanistic analyses of N-glycoside hydrolysis and N-glycoside transfer. Chem. Rev. 2006, 106, 506−555. (8) Werner, R.; Stivers, J. Kinetic isotope effect studies of the reaction catalyzed by uracil DNA glycosylase: Evidence for an oxocarbenium ion-uracil anion intermediate. Biochemistry 2000, 39, 14054−14064. (9) Lindahl, T. N-Glycosidase from Escherichia coli releases free uracil from DNA containing deaminated cytosine residues. Proc. Natl. Acad. Sci. U. S. A. 1974, 71, 3649−3653. (10) Lenz, S.; Kellie, J.; Wetmore, S. Glycosidic bond cleavage in deoxynucleotides: Effects of solvent and the DNA phosphate backbone in the computational model. J. Phys. Chem. B 2012, 116, 14275−14284. (11) Dinner, A.; Karplus, G. B. M. Uracil-DNA glycosylase acts by substrate autocatalysis. Nature 2001, 413, 752−755. (12) Drohat, A.; Stivers, T. Escherichia coli uracil DNA glycosylase: NMR characterization of the short hydrogen bond from His187 to uracil O2. Biochemistry 2000, 39, 11865−11875.

ASSOCIATED CONTENT

S Supporting Information *

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jpcb.5b05496. Small molecule model calculations, active site model calculations, data obtained from the potential energy surface studies (schemes and energies), bidimensional free energy plots, free energy profiles of connection transitions between direct and histidine mechanism, tables with calculated free energies for all mechanisms, 12378

DOI: 10.1021/acs.jpcb.5b05496 J. Phys. Chem. B 2015, 119, 12365−12380

Article

The Journal of Physical Chemistry B (13) Drohat, A. C.; Stivers, J. T. NMR Evidence for an unusually low N1 pKa for uracil bound to uracil DNA glycosylase:Implications for catalysis. J. Am. Chem. Soc. 2000, 122, 1840−1841. (14) Drohat, A.; Jagadeesh, J.; Ferguson, E.; Stivers, J. Role of electrophilic and general base catalysis in the mechanism of Escherichia coli uracil DNA glycosylase. Biochemistry 1999, 38, 11866−11875. (15) Parikh, S. S.; Walcher, G.; Jones, G. D.; Slupphaug, G.; Krokan, H. E.; Blackburn, G. M.; Tainer, J. A. Uracil-DNA glycosylase-DNA substrate and product structures: Conformational strain promotes catalytic efficiency by coupled stereoelectronic effects. Proc. Natl. Acad. Sci. U. S. A. 2000, 97, 5083−088. (16) Maiti, A.; Noon, M.; MacKerell, A. D.; Pozharski, E.; Drohat, A. Lesion processing by a repair enzyme is severely curtailed by residues needed to prevent aberrant activity on undamaged DNA. Proc. Natl. Acad. Sci. U. S. A. 2012, 109, 8091−8096. (17) Maiti, A.; Morgan, M.; Drohat, A. Role of two strictly conserved residues in nucleotide flipping and N-glycosylic bond cleavage by human thymine DNA glycosylase. J. Biol. Chem. 2009, 284, 36680− 36688. (18) Gullingsrud, J.; Saam, J.; Phillips, J. psfgen User’s Guide, 2012; http://www.ks.uiuc.edu/Research/vmd/plugins/psfgen/ug.pdf (accessed November 2013). (19) Humphrey, W.; Dalke, A.; Schulten, K. VMD: Visual molecular dynamics. J. Mol. Graphics 1996, 14, 33−38. (20) Humphrey, W., Dalke, A., Schulten, K. University of Illinois at Urbana-Champaign, http://www.ks.uiuc.edu/Research/vmd/ (accessed November 2013). (21) Li, H.; Robertson, A.; Jensen, J. Very fast empirical prediction and rationalization of protein pKa values. Proteins: Struct., Funct., Genet. 2005, 61, 704−721. (22) Phillips, J.; Braun, R.; Wang, W.; Gumbart, J.; Tajkhorshid, E.; Villa, E.; Chipot, C.; Skeel, R.; Kale, L.; Schulten, K. Scalable molecular dynamics with NAMD. J. Comput. Chem. 2005, 26, 1781−1802. (23) Vanommeslaeghe, K.; Hatcher, E.; Acharya, C.; Kundu, S.; Zhong, S.; Shim, J.; Darian, E.; Guvench, O.; Lopes, P.; Vorobyov, I.; et al. CHARMM general force field: A force field for drug-Like molecules compatible with the CHARMM all-atom additive biological force fields. J. Comput. Chem. 2010, 31, 671. (24) Jorgensen, W. L.; Chandrasekhar, J.; Madura, J. D.; Impey, R. W.; Klein, M. L. Comparison of simple potential functions for simulating liquid water. J. Chem. Phys. 1983, 79, 926−935. (25) Darden, T.; York, D.; Pedersen, L. G. Particle mesh Ewald: an Nlog(N) method for Ewald sums in large systems. J. Chem. Phys. 1993, 98, 10089−10092. (26) Brooks, B.; Brooks, C. L.; MacKerell, A. D.; Nilsson, L.; Petrella, R.; Roux, B.; Won, Y.; Archontis, G.; Bartels, C.; Boresch, S.; et al. CHARMM: The biomolecular simulation program. J. Comput. Chem. 2009, 30, 1545−1614. (27) Field, M.; Bash, P.; Karplus, M. A combined quantummechanical and molecular mechanical potential for moleculardynamics simulations. J. Comput. Chem. 1990, 11, 700−733. (28) Repasky, M.; Chandrasekhar, J.; Jorgensen, W. PDDG/PM3 and PDDG/MNDO: Improved semiempirical methods. J. Comput. Chem. 2002, 23, 1601−1622. (29) Martí, S.; Moliner, V.; Tuñoń , I. Improving the QM/MM description of chemical processes: A dual level strategy to explore the potential energy surface in very large systems. J. Chem. Theory Comput. 2005, 1, 1008−1016. (30) Field, M. The pDynamo program for molecular simulations using hybrid quantum chemical and molecular mechanical potentials. J. Chem. Theory Comput. 2008, 4, 1151−1161. (31) Bachega, J.; Timmers, L.; Assirati, L.; Bachega, L.; Field, M.; Wymore, T. GTKDynamo: A PyMOL plug-in for QC/MM hybrid potential simulations. J. Comput. Chem. 2013, 34, 2190−2196. (32) Torrie, G.; Valleau, J. Non-physical sampling distributions in Monte-Carlo free-energy estimation: Umbrella sampling. J. Comput. Phys. 1977, 23, 187−199.

(33) Kumar, S.; Bouzida, D.; Swendsen, R.; Kollman, P.; Rosenberg, J. The weighted histogram analysis method for free-energy calculations on biomolecules 1. The method. J. Comput. Chem. 1992, 13, 1011− 1021. (34) van der Spoel, D.; Lindahl, E.; Hess, B.; Groenhof, G.; Mark, A.; Berendsen, H. GROMACS: Fast, flexible, and free. J. Comput. Chem. 2005, 26, 1701−1718. (35) PyMOL Molecular Graphics System, Schrödinger LLC. PyMOL Molecular Graphics System; http://www.pymol.org (accessed December 2013). (36) This summation of individual free energy differences assumes an integration over the same phase space volume for all states. However, this is only the case for those states that have been sampled by the same free energy calculation and by following the same reaction coordinates in the simulations with biasing potentials. (37) These two transitions have been sampled with the same reaction coordinate, and their free energies can therefore be compared directly. (38) Note that the free energies relative to the reactant state (blue and red numbers in Figure 3b) for intermediate and transition states, respectively) should be unique for one state, regardless of the way that led to this state. The existing differences can be explained by rounding, limited sampling leading to hysteresis effects, and the intrinsic error made by adding free energies that have been obtained from simulations sampling different volumes in phase space. (39) Note that the corresponding states in the histidine mechanism of wild-type TDG, Bh and Fh, respectively; both have free energies comparable to the wild-type reactant state R which serves as a reference state for the direct and the histidine mechanism. (40) Miyagi, M.; Wan, Q.; Ahmad, M.; Gokulrangan, G.; Tomechko, S.; Bennett, B.; Dealwis, C. Histidine hydrogen-deuterium exchange mass spectrometry for probing the microenvironment of histidine residues in dihydrofolate reductase. PLoS One 2011, 6, e17055. (41) Takahashi, T.; Nakamura, H.; Wada, A. Electrostatic forces in two lysozymes: Calculations and measurements of histidine pKa values. Biopolymers 1992, 32, 897−909. (42) Edgcomb, S. P.; Murphy, K. P. Variability in the pKa of histidine side-chains correlates with burial within proteins. Proteins: Struct., Funct., Genet. 2002, 49, 1−6. (43) Maiti, A.; Drohat, A. Dependence of substrate binding and catalysis on pH, ionic strength, and temperature for thymine DNA glycosylase: Insights into recognition and processing of G:T mispairs. DNA Repair 2011, 10, 545−553. (44) Bennett, M.; Rodgers, M.; Hebert, A.; Ruslander, L.; Eisele, L.; Drohat, A. Specificity of human thymine DNA glycosylase depends on N-glycosidic bond stability. J. Am. Chem. Soc. 2006, 128, 12510− 12519. (45) Shapiro, R.; Kang, S. Uncatalyzed hydrolysis of deoxyuridine, thymidine and 5-bromodeoxyuridine. Biochemistry 1969, 8, 1806− 1810. (46) Shapiro, R.; Danzig, M. Acidic hydrolysis of deoxycytidine and deoxyuridine derivatives. General mechanism of deoxyribonucleoside hydrolysis. Biochemistry 1972, 11, 23−29. (47) Wu, R.; Gong, W.; Liu, T.; Zhang, Y.; Cao, Z. QM/MM molecular dynamics study of purine-specific nucleoside hydrolase. J. Phys. Chem. B 2012, 116, 1984−1991. (48) Chen, N.; Zhao, Y.; Lu, J.; Wu, R.; Cao, Z. Mechanistic insights into the rate-limiting step in purine-specific nucleoside hydrolase. J. Chem. Theory Comput. 2015, 11, 3180−3188. (49) Das, R. S.; Samaraweera, M.; Morton, M.; Gascon, J. A.; Basu, A. K. Stability of N-glycosidic bond of (5′S)-8,5′-cyclo-2′-deoxyguanosine. Chem. Res. Toxicol. 2012, 25, 2451−2461. (50) Rios-Font, R.; Rodriguez-Santiago, L.; Bertran, J.; Sodupe, M. Influence of N7 protonation on the mechanism of the N-glycosidic bond hydrolysis in 2prime-deoxyguanosine. A theoretical study. J. Phys. Chem. B 2007, 111, 6071−6077. (51) Rios-Font, R.; Bertran, J.; Sodupe, M.; Rodriguez-Santiago, L. On the mechanism of the N-glycosydic bond hydrolysis of 2′deoxyguanosine: insights from first principles calculations. Theor. Chem. Acc. 2011, 128, 619−626. 12379

DOI: 10.1021/acs.jpcb.5b05496 J. Phys. Chem. B 2015, 119, 12365−12380

Article

The Journal of Physical Chemistry B (52) Hashimoto, H.; Hong, S.; Bhagwat, A. S.; Zhang, X.; Cheng, X. Excision of 5-hydroxymethyluracil and 5-carboxylcytosine by the thymine DNA glycosylase domain: its structural basis and implications for active DNA demethylation. Nucleic Acids Res. 2012, 40, 10203− 10214. (53) Hashimoto, H.; Zhang, X.; Cheng, X. Selective excision of 5carboxylcytosine by a thymine DNA glycosylase mutant. J. Mol. Biol. 2013, 425, 971−976. (54) Hardeland, U.; Bentele, M.; Jiricny, J.; Schär, P. Separating substrate recognition from base hydrolysis in human thymine DNA glycosylase by mutational analysis. J. Biol. Chem. 2000, 275, 33449− 33456. (55) Barrett, T.; Savva, R.; Panayotou, G.; Barlow, T.; Brown, T.; Jiricny, J.; Pearl, L. Crystal structure of a G:T/U mismatch-specific DNA glycosylase: mismatch recognition by complementary-strand interactions. Cell 1998, 92, 117−129. (56) Xiao, G.; Tordova, M.; Jagadeesh, J.; Drohat, A. C.; Stivers, J. T.; Gilliland, G. L. Crystal structure of Escherichia coli uracil DNA glycosylase and its complexes with uracil and glycerol: Structure and glycosylase mechanism revisited. Proteins: Struct., Funct., Genet. 1999, 35, 13−24. (57) Slupphaug, G.; Mol, C. D.; Kavli, B.; Arvai, A. S.; Krokan, H. E.; Tainer, J. A. A nucleotide-flipping mechanism from the structure of human uracil-DNA glycosylase bound to DNA. Nature 1996, 384, 87− 92. (58) Kavli, B.; Slupphaug, G.; Mol, C.; Arvai, A. S.; Peterson, S. B.; Tainer, J.; Krokan, H. Excision of cytosine and thymine from DNA by mutants of human uracil-DNA glycosylase. EMBO J. 1996, 15, 3442− 3447. (59) Lee, D.-H.; Liu, Y.; Lee, H.-W.; Xia, B.; Brice, A. R.; Park, S.-H.; Balduf, H.; Dominy, B. N.; Cao, W. A structural determinant in the uracil DNA glycosylase superfamily for the removal of uracil from adenine/uracil base pairs. Nucl. Acid. Res. 2015, 43, 1081−1089. (60) Stivers, J.; Drohat, A. Uracil DNA glycosylase: Insights from a master catalyst. Arch. Biochem. Biophys. 2001, 396, 1−9. (61) Parker, J.; Stivers, J. Uracil DNA glycosylase: revisiting substrate-assisted catalysis by DNA phosphate anions. Biochemistry 2008, 47, 8614−8622. (62) Zhu, B. B.; Zheng, Y.; Hess, D.; Angliker, H.; Schwarz, S.; Siegmann, M.; Thiry, S.; Jost, J. 5-methylcytosine-DNA glycosylase activity is present in a cloned G/T mismatch DNA glycosylase associated with the chicken embryo DNA demethylation complex. Proc. Natl. Acad. Sci. U. S. A. 2000, 97, 5135−5139. (63) Hardeland, U.; Bentele, M.; Jiricny, J.; Schär, P. The versatile thymine DNA-glycosylase: a comparative characterization of the human, Drosophila and fission yeast orthologs. Nucl. Acid. Res. 2003, 31, 2261−2271.

12380

DOI: 10.1021/acs.jpcb.5b05496 J. Phys. Chem. B 2015, 119, 12365−12380