Uniform Free-Energy Profiles of the P-O Bond Formation and Cleavage Reactions Catalyzed by DNA Polymerases β and λ.

Human X-family DNA polymerases β (Polβ) and λ (Polλ) catalyze the nucleotidyl-transfer reaction in the base excision repair pathway of the cellular DNA damage response. Using empirical valence bond and free-energy perturbation simulations, we explore the feasibility of various mechanisms for the deprotonation of the 3'-OH group of the primer DNA strand, and the subsequent formation and cleavage of P-O bonds in four Polβ, two truncated Polλ (tPolλ), and two tPolλ Loop1 mutant (tPolλΔL1) systems differing in the initial X-ray crystal structure and nascent base pair. The average calculated activation free energies of 14, 18, and 22 kcal mol-1 for Polβ, tPolλ, and tPolλΔL1, respectively, reproduce the trend in the observed catalytic rate constants. The most feasible reaction pathway consists of two successive steps: specific base (SB) proton transfer followed by rate-limiting concerted formation and cleavage of the P-O bonds. We identify linear free-energy relationships (LFERs) which show that the differences in the overall activation and reaction free energies among the eight studied systems are determined by the reaction free energy of the SB proton transfer. We discuss the implications of the LFERs and suggest pKa of the 3'-OH group as a predictor of the catalytic rate of X-family DNA polymerases.

Introduction

Preserving genetic and epigenetic information in DNA yet allowing for mutation and recombination, two essential properties of life, to occur, requires a proper balance between DNA mutagenesis and DNA repair processes.[1] Two X-family DNA polymerases,[2] β (Polβ) and λ (Polλ), have been implicated in numerous DNA repair pathways including base excision repair (BER) of oxidative and alkylation damage.[3−16] Both Polβ and Polλ possess 5′-deoxyribosephosphate (dRP) lyase and DNA polymerase activities[8,9,17,18] conferred by two different domains.[19−25] The lyase domain creates a single-nucleotide gap by removing dRP and the polymerase domain inserts a specific 2′-deoxyribonucleotide 5′-triphosphate (dNTP) into the gap according to the identity of the template nucleotide. The fidelity, determined by the ratio between correct and incorrect dNTP insertion events, depends on the proof-reading capacity and reflects the specific biological role of each DNA polymerase. Both Polβ and Polλ are low-fidelity DNA polymerases, in part because they lack an intrinsic 3′ → 5′ exonuclease activity.[26−31] To advance our understanding of the structural basis of the catalytic power and the fidelity determinants of Polβ and Polλ, we need to explore the reaction mechanisms of these enzymes at the atomic level of detail, which can be achieved by complementing experimental (structural and kinetic) observations with the insight obtained using theoretical calculations.[32] Further progress in this field will eventually lead to a thorough realization of the principles that govern the balance between the protective and harmful involvement of DNA polymerases in carcinogenesis.[33,34]

The lyase and polymerase domains of human Polβ and Polλ are evolutionarily homologous, sharing 35.9% amino acid sequence identity (Clustal Omega[35]). Polλ contains an additional, N-terminally located, nuclear localization motif (Met1–Glu35), followed by a BRCT domain (Glu36–Arg132) and a proline-rich domain (Arg133–Pro244),[36] none of which is, however, essential for the lyase and polymerase activities of Polλ.[22,29,37] X-ray crystal structures of the ternary complexes of Polβ (Thr10–Glu335; PDB code: 2FMP) and truncated Polλ (tPolλ; comprising the lyase and polymerase domains, Ser245–Trp575; PDB code: 2PFO) with single-nucleotide gapped dsDNA, dNTP, and two active site cations superimpose with a root-mean-square deviation (RMSD)Cα of 1.4 Å (Figure ). The structures differ in amino acid sequence, length, and geometry of the surface loop between β-strands 3 and 4 (Loop 1). This loop modulates the fidelity of tPolλ by controlling the transition from the tPolλ·DNA binary complex to the tPolλ·DNA·dNTP·2 Mg2+ ternary complex, which is accompanied by a conformational change in the template DNA strand.[38,39] A chimeric tPolλ that contains a shorter, Polβ-like Loop1 (tPolλΔL1) assumes the tPolλ-like catalytically competent conformation for the template DNA strand in the ternary complex (Figure ) but differs from tPolλ in the inactive conformation of the binary complex. In contrast, the transition from the inactive to the active conformation in Polβ involves an open-to-closed conformational change in the polymerase domain, whereas the template strand remains in the reaction-ready conformation in both binary and ternary complexes.[21]

Figure 1

Crystal structures of the ternary complexes of human DNA polymerases β and λ. White cartoon, protein; black cartoon, DNA; orange spheres, template nucleotide; green spheres, 3′-terminal primer nucleotide; red spheres, dNTP; yellow spheres, metal ions; magenta spheres, catalytic triad (Cα atoms); cyan spheres, Loop 1 (Cα atoms).

The rate-limiting nucleotidyl-transfer reaction[19,29,30,40,41] (Figure ) is initiated in the active ternary complex by deprotonation of the 3′-OH group of the primer DNA strand. The loss of proton is facilitated by catalytic Mg2+, which lowers pKa of the hydroxyl group and stabilizes the evolving negative charge of the nucleophilic oxygen (Onuc). The proton acceptor can be either an active site group via a general base (GB) mechanism or a hydroxide anion in bulk water via a specific base (SB) mechanism. In the latter scenario, the proton transfer can be described by a pH-dependent equilibrium between the protonated and deprotonated Onuc. The activated Onuc forms a covalent bond with Pα of dNTP and thus extends the primer DNA strand by one nucleotide. Concomitantly, the cleavage of the Pα–Olg bond in the triphosphate group of dNTP generates an inorganic pyrophosphate (PPi), the byproduct of the reaction. Despite broad consensus on the overall mechanism of the nucleotidyl-transfer reaction,[19] three important questions remain: (1) What is the preferred mechanism for the deprotonation of the 3′-OH group of the primer DNA strand? (2) What is the degree of concertedness of the P–O bond formation and cleavage reactions? (3) Which of the three reaction steps limits the overall reaction rate?

Figure 2

Nucleotidyl-transfer reaction catalyzed by Polβ, tPolλ, and tPolλΔL1. Atoms and ions are colored by element: light green, carbon; light red, oxygen; white, hydrogen; light blue, nitrogen; light orange, phosphorus; light gray, magnesium ion. A, catalytic Mg2+; B, structural Mg2+; W, oxygen atom of water molecule; N, nitrogen atom of the nucleobase of the primer nucleotide; N′, nitrogen atom of the nucleobase of dNTP; α1, Cα of Asp256/490/485 in Polβ/tPolλ/tPolλΔL1; α2, Cα of Asp192/429/429 in Polβ/tPolλ/tPolλΔL1; α3, Cα of Asp190/427/427 in Polβ/tPolλ/tPolλΔL1; a–f, proton transfer options. The atoms assigned to the reactive region in MD simulations of the reaction in this study are labeled by atom identifiers (1–15); see the Methods section for details.

Several theoretical studies have already addressed these challenging questions by employing hybrid quantum mechanics/molecular mechanics (QM/MM), free-energy perturbation (FEP), or empirical valence bond (EVB) simulation methods to Polβ[30,42−48] and Polλ,[49−51] but definitive answers are yet to be reached: (1) An ordered active site water molecule, two active site catalytic aspartates, a nonbridging oxygen of the Pα-group, and bulk water have been proposed as the proton accceptors.[30,42−44,49,51] (2) Associative and concerted P–O bond formation and cleavage reactions have been suggested.[30,43,48,49] (3) Either the proton transfer or the P–O bond formation/cleavage reaction has been reported as the rate-limiting step.[30,43−45,49]

This article is a contribution to this ongoing debate. Using FEP[52] and EVB[53,54] methods, we examine the nucleotidyl-transfer reaction in water, Polβ, tPolλ, and tPolλΔL1 with the aim of (1) obtaining the free-energy profiles of the reaction; (2) calculating the catalytic rate constants (kpol) using the transition-state (TS) theory approximation,[55] comparing them with available experimental values (Polβ and tPolλ), and offering them as genuine predictions (tPolλΔL1); (3) assessing the relative feasibility of the proton transfer to active site aspartates (Asp256 and Asp192 in Polβ and their corresponding Asp490/485 and Asp429/429 in tPolλ/tPolλΔL1) and to bulk water; (4) exploring the concertedness and tightness of the P–O bond formation and cleavage reactions by employing a trigonal-bypiramidal geometry around the Pα atom of dNTP in the TS model;[56,57] and (5) determining the dependence of the calculated energetics on the initial atomic coordinates and on the identity of the nascent base pair (Table ). Our free-energy calculations suggest the same reaction mechanism of the nucleotidyl-transfer reaction and similar energetics of the P–O bond formation and cleavage reactions in Polβ and Polλ and indicate linear free-energy relationships (LFERs) as potentially useful concept for further theoretical and experimental studies aimed at deepening our insight into DNA polymerase catalysis and fidelity.

Table 1

Model Systems

enzyme	system	XRCa	dNpb	dNTP·dNtc	Ad	Bd
Polβ	Polβ1	2FMP	dC	dCTP·dG	Mg2+	Mg2+
	Polβ2	2PXI	dC	dCTP·dG	Mg2+	Mg2+
	Polβ3	2FMP	dC	dGTP·dC	Mg2+	Mg2+
	Polβ4	2PXI	dC	dGTP·dC	Mg2+	Mg2+
tPolλ	Polλ1	2PFO	dC	dCTP·dG	Mg2+	Mg2+
	Polλ2	2PFP	dC	dCTP·dG	Mg2+	Mg2+
tPolλΔL1	Polλ3	3PML	dC	dGTP·dC	Mg2+	Mg2+
	Polλ4	3PML	dC	dGTP·dT	Mg2+	Mg2+

XRC, X-ray crystal structure (PDB code).

dNp, 3′-terminal nucleotide of the primer DNA strand.

dNt, templating nucleotide of the template DNA strand.

A and B, active site metal ions (see Figure ).

Methods

Structural Models

The structures of 12 water systems and 24 enzyme systems (Tables S1–S3) were prepared using PyMOL 0.99[58] and Coot 0.6.2[59] molecular graphics programs. The Qprep 5.04 module of the molecular dynamics (MD) package Q[60] was then used to immerse the structures in a sphere of TIP3P[61] water molecules. The sphere, which measured 48 Å in diameter, was centered at the C4′ atom of dNTP. Amino acid residues and nucleotides located outside the sphere (Figure S1) were excluded from nonbonded interactions and restrained to their initial coordinates by a harmonic potential with a force constant of 200 kcal mol–1 Å–2. Each solvated structure (Supporting Files 1 and 2) was equilibrated by employing the ff94 AMBER force field[62] in a series of 12 MD simulations using the Qdyn 5.06 module of Q (Table S4). The equilibrating MD simulation protocol was identical to that of our previous studies:[63,64] Water molecules were subjected to the surface constraint all-atom solvent (SCAAS)[65] boundary condition. The SHAKE algorithm[66] was applied to all hydrogen atoms. Nonbonded interactions were evaluated explicitly and using the local reaction field method[67] within and beyond 10 Å cutoff distance, respectively. The values of stepsize, temperature, and external bath coupling[68] parameters were gradually increased from 0.01 to 2 fs, 5 to 298 K, and 0.1 to 40 fs, respectively. The equilibrated structures were subsequently used in (1) EVB simulations of the GB proton transfer, Onuc–Pα bond formation, and Pα–Olg bond cleavage reactions, and (2) FEP simulations of the SB proton transfer (Table S5).

EVB Simulations

The EVB method (eqs S1–S19) is based on the assumption that the most important effect of the environment (water vs enzyme) is the change in the relative energies of the covalent and ionic valence bond (VB) states.[53] The application of the EVB method thus comprises five steps: (1) Defining VB states. (2) Sampling the potential energy surface between adjacent VB states in water and enzyme systems by FEP simulations. (3) Obtaining free-energy profile (energy gap vs EVB free energy) of the reaction in the water system by superimposing and mixing the potential energy curves of adjacent VB states and reproducing the experimental or ab initio activation and reaction free energies. (4) Applying the water-derived superimposing and mixing parameters to VB potential energy curves of the enzyme system. (5) Calculating the quantitative catalytic effect of the enzyme by comparing its activation and reaction free energies to those of the water system.

Note that the FEP simulations in step 2 are practically performed on the linear combination of molecular mechanical potentials of the two states involved in the given reaction step. EVB free energies in step 3 are generated via the umbrella sampling technique. This technique combines a statistical ensemble, obtained in step 2, with EVB energies obtained by the diagonalization of the EVB Hamiltonian. Although the umbrella sampling step could be, in principle, avoided if the FEP simulations were performed using the EVB forces rather than the classical forces, such simulations would be slower. Perhaps more importantly, they would have to be repeated many times for the correct values of the mixing and superimposing terms in the EVB Hamiltonian to be obtained, so as to reproduce the experimental activation and reaction free energies of the reference reaction.

To model the nucleotidyl-transfer reaction in the water, Polβ, and Polλ systems, a reactive region was designated as consisting of 15 and 11 atoms in the GB and SB pathways, respectively (Figure ). Four VB states of the reactive region were defined (Tables and S6–S24): (1) the reactant state (RS) represented by the catalytically competent Polβ and Polλ systems, and the corresponding reference water systems; (2) the proton transfer intermediate state (IS1); (3) an artificial intermediate state (IS) along the concerted P–O bond formation and cleavage reactions (IS2); and (4) the product state (PS) with the primer DNA strand extended by one nucleotide and PPi fully formed. The artificial IS2 allowed us to explore the TS region of the free-energy surface in a thorough and controlled way.

Table 2

EVB/FEP States for Modeling of the GB and SB Pathways of the Nucleotidyl-Transfer Reaction in the Water, Polβ, and Polλ Systemsa

state	GB	SB
RS	Oδ(−0.8514) + H3t–Onuc + Pα–Olg	H3t–Onuc + Pα–Olg
IS1	Oδ–H3t + Onuc(−0.8549) + Pα–Olg	Onuc(−0.8549) + Pα–Olg
IS2a	Oδ–H3t + Onuc···(1.90 Å)···Pα···(1.90 Å)···Olg	Onuc···(1.90 Å)···Pα···(1.90 Å)···Olg
IS2b	Oδ–H3t + Onuc···(2.28 Å)···Pα···(2.25 Å)···Olg	Onuc···(2.28 Å)···Pα···(2.25 Å)···Olg
PS	Oδ–H3t + Onuc–Pα + Olg(−0.8100)	Onuc–Pα + Olg(−0.8100)

Each state is described schematically by the atoms or atomic groups most directly involved in the formation and cleavage (···) of the covalent bonds (−). Atomic charges of selected atoms are indicated in superscripts, and distances between the oxygen and phosphorus atoms in IS2 and IS2b are given in parentheses. See Tables S6–S28 for the full description of the EVB/FEP states. The missing proton in the IS1, IS2, and PS states in the SB pathway was compensated by sodium and/or chloride ions in the water systems, and by a charged arginine or a sodium ion in the Polβ and Polλ systems (see Tables S1 and S3).

To model the compact and loose associative P–O bond formation and cleavage reactions, two IS2 VB states, IS2a and IS2b, differing in Onuc···Pα···Olg distances and, consequently, in the charge distribution were defined. The atomic charges of RS, IS1, IS2a, and PS were obtained from Florián et al.[56] These charges are consistent with the atom-centered point charge model of the ff94 AMBER force field, which is based on molecular electrostatic potentials obtained using the Hartree–Fock method and the 6-31G* basis set. The atomic charges of IS2b were based on IS2a and adjusted to account for the difference between PCM/B3LYP/TZVP//HF/6-31G* Mulliken atomic charges in QM models of the compact and loose associative-like TS between IS1 and IS2.[69] To achieve net zero charge of the reactive region in IS1, IS2, and PS in the SB pathway, the atomic charge of C3′ in 2′-deoxyribose (dRib; water systems) and dC (Polβ and Polλ systems) was further adjusted (cf. Tables S6 and S7).

The FEP method was employed to sample the transition between adjacent VB states in both directions using Qdyn: RS ↔ IS1 ↔ (IS2a or IS2b) ↔ PS in the GB pathway and IS1 ↔ (IS2a or IS2b) ↔ PS in the SB pathway. The transition was carried out in 51 k steps by gradually changing the value of the FEP mapping parameter λ from 1 to 0; Δλ = λ – λ = −0.02. Each FEP step consisted of 10 ps MD simulation at 298 K with a stepsize of 1 fs; MD parameters were identical, except for the stepsize, to the equilibrating MD simulation XII (see Table S4). Mapping and perturbation potential energies of the reactive region with the surrounding environment were recorded every 10 fs with no cutoff applied to nonbonded interactions. The recorded potential energies, 51 × 999 = 50 949 for each FEP transition, were processed using the Qfep 5.01 module of Q: EVB free energies (y coordinate) were calculated as a function of the energy gap (x coordinate, consisting of 100 and 60 bins in the RS ↔ IS1 and IS1 ↔ IS2 ↔ PS reactions, respectively) using the umbrella sampling method. The coupling of adjacent VB and VB potential surfaces and difference in their minima were modulated using EVB parameters H and Δα to set the activation and reaction free energies in the water systems in accordance with Florián et al.[57] See eqs S1–S19 for a detailed description of the EVB method.

FEP Simulations

The SB proton transfer free energy, which corresponds to the free energy of IS1, was calculated from pKa = 12.67 for the 3′-OH group of dRib[70] and from FEP simulations of the hydrogen annihilation and creation in the water, Polβ, and Polλ systems. The reactive region consisted of the C3′, Onuc, and H3t atoms of dRib in the water systems or dC in the Polβ and Polλ systems (Tables S25–S28). Two FEP states were defined: (1) RS with the 3′-OH group and (2) intermediate state 1 (IS1) with Onuc– (Table ). The transition between RS and IS1 was carried out in 51 k FEP steps, each consisting of 10 ps MD simulation, using Qdyn. The sampling was performed in both directions (RS ↔ IS1); the hydrogen creation was initiated from the final coordinates of the hydrogen annihilation MD simulation. All MD parameters were identical to the sampling MD simulations in the EVB approach. The reaction free energy was calculated as the cumulative sum of average differences between mapping and perturbation potential energies using Zwanzig’s FEP formula.[52] See eqs S1–S5 and S20 for a detailed description of the FEP method.

Calibration of the Free-Energy Profile of the RS → IS1 Reaction in Solution

The EVB reaction free energy of the GB proton transfer in the reference water system was set to 10.8 kcal mol–1, according to the difference in pKa between the hydrogen donor (dRib) and acceptor (acetate)where R = 1.9872041 × 10–3 kcal mol–1 K–1 is the gas constant, T = 298 K is the simulation temperature, pKa,dRib = 12.67,[70] and pKa,acetate = 4.74. The EVB activation free energy of the GB proton transfer was set to 12.6 kcal mol–1, that is, 1.8 kcal mol–1 above the reaction free energy,[57] which is in accordance with the assumption of a low activation barrier for proton transfer between two oxygen atoms.[71]

The FEP reaction free energy of the SB proton transfer in the reference water system was set to 7.7 kcal mol–1, according to the difference between pKa of dRib and pH of 7The differences in pKa between the enzyme and water systems were calculated from FEP free energies (see eq S20)and used for the calculation of the FEP reaction free energy of the RS → IS1 step in the enzyme systemsThe value of ΔgFEP,water in eq is a weighted average, based on Mg2+ occupancy, of the FEP free energies obtained from FEP simulations of the water system with/without Mg2+ ions. The Mg2+ occupancy was set to 0.7, which yielded ΔpKa for Polβ3 and Polβ4 in accordance with the difference between the approximate inflection point in the kpol versus pH plot of 8.3 for dGTP·dC nascent base pair in Polβ[41] and the pKa of dRib.

Calibration of the Free-Energy Profile of the IS1 → IS2 Reaction in Solution

The EVB activation free energy of the IS1 → IS2 step in the reference water system was set to 23.1 kcal mol–1,[57] which was calculated as the sum of three contributionswhere Δgcpmp‡ = 27.5 kcal mol–1 is the activation free energy for the hydrolysis of (3-chlorophenyl) methyl phosphate, whose pKa = 9.02[72] is similar to pKa4 = 8.93 of PPi,[73] by hydroxide anion in solution;[74,75] Δgcage‡ = −2.4 kcal mol–1 is the solvent cage effect;[57] and Δgmg‡ = −2.0 kcal mol–1 is the correction for the presence of Mg2+ ions in the reference water system.[76] The EVB reaction free energy of 19.5 kcal mol–1 for the IS1 → IS2 step was obtained as followswhere ΔgEVB,water‡ = 23.1 kcal mol–1 is the activation free energy of the IS1 → IS2 step (eq ), ΔΔg0‡ = −2.6 kcal mol–1 is the difference between the reaction and activation free energies used by Florián et al.,[57] and c = −1.0 kcal mol–1 is the correction introduced in this study to make the intermediate observable on the EVB free-energy surface and to change the rate-limiting step in aqueous solution from the IS2 → PS to the IS1 → IS2 step, in accordance with ref (69).

Calibration of the Free-Energy Profile of the IS2 → PS Reaction in Solution

The EVB activation free energy of the IS2 → PS step in the reference water system was set to 2.6 kcal mol–1, the value used by Florián et al.[57] The EVB reaction free energy of the IS2 → PS step in the reference water system was set to −32.8 kcal mol–1 in the GB pathway and −29.7 kcal mol–1 in the SB pathway (−32.8 + 10.8 – 7.7 kcal mol–1; cf. eqs and 2), which yielded the intended cumulative RS → PS reaction free energy of −2.5 kcal mol–1 in both pathways, a value in accordance with the equilibrium constant of ∼102 observed for the extension of a single-stranded DNA by one nucleotide, dA + dATP ↔ dA + PPi, catalyzed by terminal deoxynucleotidyltransferase.[77]

Calculation of the Catalytic Effect

The free-energy profiles of the entire nucleotidyl-transfer reaction were created by combining the free-energy profiles of individual steps: RS ↔ IS1, IS1 ↔ IS2, and IS2 ↔ PS, where IS2 was either IS2a or IS2b. The value of IS1 was obtained from EVB (in the GB pathway) or FEP (in the SB pathway). The IS1 → PS part in the SB pathway was calculated as the average of two versions, SB1 and SB2, which differed in the composition of the simulated system (see Tables S1, S3, and S5). The catalytic effect was then calculated as the difference in activation free energies between reactions in the enzyme and water systems. To facilitate comparison of the calculated and observed reaction energetics, the calculated activation free energies (kcal mol–1) were converted into the catalytic rate constants (s–1; eq ) and the observed rate constants were converted into activation free energies (eq )[55]where T (K) is the simulation (eq ) or experimental temperature (eq ), kB = 1.3806488 × 10–23 m2 kg s–2 K–1 is the Boltzmann constant, and h = 6.62606957 × 10–34 m2 kg s–1 is the Planck constant.

Structural Analysis of MD Trajectories

VMD[78] and PyMOL molecular graphics programs were used for the visualization of MD trajectories and individual structures, respectively. The geometrical parameters, including RMSD, interatomic distances, and angles, were extracted from MD trajectories using VMD. PyMOL was used for creating all molecular structure images for this article. Onuc–Pα and Pα–Olg distances in TSs were calculated as weighted averages of all FEP k steps that contributed to the TS bin. (Note that configurations from different k steps may contribute to one EVB bin, and, vice versa, configurations from the same k step may belong to different EVB bins; see eq S16.) Animations of the nucleotidyl-transfer reaction were generated from the final snapshots of the FEP steps that corresponded to each state (FEP simulations) or contributed by the highest number of energy points to each state (EVB simulations).

Results

Free-Energy Profiles of the Proton Transfer Reaction

Two mechanisms for the deprotonation of the 3′-OH group of the primer DNA strand were considered: (1) GB proton transfer to the proximal and distal carboxylic oxygens of Asp256/490/485 (labeled Oδ2 and Oδ1, respectively, in Figure ) and to the proximal carboxylic oxygen (Oδ2) of Asp192/429/429 in the Polβ/tPolλ/tPolλΔL1 systems and (2) SB proton transfer to bulk water.

GB Proton Transfer (RS → IS1)

The proton transfer to the active site aspartates in enzyme systems and to the corresponding acetate ions in the reference water system (Figure S2; Table S29) was modeled using the EVB approach: The reactants and products were defined using RS and IS1 VB states, and the RS → IS1 transformation was carried out using the FEP method. The initial set of six enzyme systems comprised Polβ1, Polβ3, Polλ1, Polλ2, Polλ3, and Polλ4. Oδ2 of Asp192/429/429 was the least feasible proton acceptor with minimum and maximum activation free energy (Δg‡) values of 24.5 and 32.1 kcal mol–1, respectively, among the six enzyme systems. Of the two carboxylic oxygens of Asp256/490/485, Oδ1, which accepted the proton with the minimum, maximum, and average Δg‡ of 16.7, 23.4, and 19.1 kcal mol–1, respectively, was surpassed in feasibility by Oδ2, which accepted the proton with the minimum, maximum, and average Δg‡ values of 10.6, 13.7, and 12.0 kcal mol–1, respectively. The proton transfer to Oδ2 of Asp256/490/485 also showed the lowest reaction free energies (Δg0) with an average of 6.5 kcal mol–1 among the six enzyme systems, which was 10.5 and 7.5 kcal mol–1 lower than the average Δg0 for the two less favorable proton acceptors, Oδ2 of Asp192/429/429 and Oδ1 of Asp256/490/485, respectively.

To increase the robustness of the calculated Δg‡ and Δg0 for the preferred proton transfer to Oδ2 of Asp256/490/485, EVB simulations were also performed in the reverse direction (RS ← IS1), and Polβ2 and Polβ4 were added to the set of enzyme systems. The reverse reaction failed to proceed in Polβ4, because the Morse potential for the distance between Onuc and the proton did not provide the sufficient driving force needed for overcoming the energy barrier caused by a strong electrostatic attraction between the proton and Oδ2 of Asp190. In other enzyme systems, the values of both Δg‡ and Δg0 were on average 1.2 and 1.8 kcal mol–1 lower than the corresponding free energies in the forward direction. The final values of Δg‡ and Δg0 were obtained by averaging the free energies from the forward and reverse proton transfer (RS ↔ IS1). The GB proton transfer in the Polβ and Polλ systems showed average catalytic effects of 2.7 and 0.2 kcal mol–1, respectively, when compared to that of the reference water system calibrated to Δg‡ = 12.6 kcal mol–1. The proton transfer product state (IS1) in the Polβ and Polλ systems was on average 6.0 and 4.4 kcal mol–1 more stable than the corresponding IS1 in the water system calibrated to Δg0 = 10.8 kcal mol–1. Thus, the GB proton transfer reaction was endergonic with average free-energy cost of 4.8 and 6.4 kcal mol–1 in the Polβ and Polλ systems, respectively.

Specific Base (SB) Proton Transfer (RS → IS1)

The proton transfer to bulk water was modeled using the FEP approach: The reactants and products were defined by the RS and IS1 FEP states; the RS ↔ IS1 transformation was carried out in the forward (proton annihilation) and reverse (proton creation) directions. The obtained FEP free energies were averaged and converted into the SB proton transfer reaction free energies on the basis of the differences in pKa of the 3′-OH group between the enzyme and reference water systems (Table S30). The SB proton transfer reaction was endergonic with average reaction free energies of 1.3, 6.6, and 10.2 kcal mol–1 in Polβ, tPolλ, and tPolλΔL1 systems, respectively. Compared to the reference water system, Δg0 values were 6.4 kcal mol–1 lower in Polβ, 1.1 kcal mol–1 lower in tPolλ, and 2.5 kcal mol–1 higher in tPolλΔL1. Compared to the GB proton transfer, Δg0 values were 3.5 kcal mol–1 lower in Polβ, 0.2 kcal mol–1 higher in tPolλ, and 3.9 kcal mol–1 higher in tPolλΔL1.

Free-Energy Profiles of the Concerted Nucleophilic Attack and PPi Formation Reactions

The concerted Onuc–Pα formation and Pα–Olg cleavage reactions were split into two successive steps by introducing an artificial VB state (IS2) between the proton transfer product (IS1) and the final product of the nucleotidyl-transfer reaction (PS): IS1 → IS2 → PS. Two variants of IS2 were employed for exploring the tightness of the TS of the overall reaction: IS2a for the compact TS and IS2b for the loose TS. IS2a and IS2b differed in the atomic charges in the reactive region (−0.65e, 1.16e, and −0.55e for Onuc, Pα, and Olg in IS2a vs −0.75e, 1.22e, and −0.65e for the corresponding atoms in IS2b; see Tables S6 and S7) and in the equilibrium P–O distances in the Morse potential energy function (1.90 Å for Onuc–Pα and Pα–Olg distances in IS2a vs 2.28 Å for Onuc–Pα and 2.25 Å for Pα–Olg distances in IS2b; see Tables S13–S15).

Part I of the P–O Bond Formation and Cleavage Reactions (IS1 → IS2)

In the GB pathway, the IS1 → IS2a reaction proceeded with an average Δg‡ of 17.6 kcal mol–1 among Polβ1, Polβ2, and Polβ4, and 15.7 kcal mol–1 among Polλ. Polβ3 was an outlier with an average Δg‡ of 21.9 kcal mol–1 (squares and circles in Figure S3; Table S31). The IS1 → IS2b reaction proceeded with average Δg‡ values of 18.5 and 16.2 kcal mol–1 in Polβ and Polλ systems, respectively (triangles and pentagons in Figure S3; Table S32). IS2b was less stable (by 1.5 and 2.2 kcal mol–1 in Polβ and Polλ systems, respectively) than IS2a. However, considering both Δg‡ and Δg0, the energetics of the IS1 → IS2a and IS1 → IS2b reactions were too similar to allow for a confident designation of a preferred IS2 geometry.

The IS1 → IS2 → PS reaction steps following the SB proton transfer were modeled using two different approaches to compensate for the total charge of the systems for the lost proton: (1) setting the net charge of the remote Arg152/389/389, located >14 Å from the Pα atom of dNTP in Polβ/tPolλ/tPolλΔL1, to +1e (SB1; Figure S4A; Tables S33 and S34) or (2) adding one unrestrained Na+ ion into the active site, about 4 Å from the Pα atom of dNTP (SB2; Figure S4B; Tables S35 and S36). The energetics of the IS1 → IS2 → PS reaction steps obtained with the SB1 and SB2 structures could not be confidently distinguished (Tables S33–S36 and S41–S44); therefore, the Δg‡ and Δg0 values were averaged (Tables S37, S38, S45, and S46).

The IS1 → IS2a reaction in the SB pathway (squares and circles in Figure S4; Table S37) proceeded with an average Δg‡ of 12.8 kcal mol–1 in Polβ systems and 11.8 kcal mol–1 in Polλ systems. Among the eight enzyme systems, the maximum difference in Δg‡ was only 2.2 kcal mol–1, and the standard deviation of the average Δg‡ of 12.3 kcal mol–1 was only 0.3 kcal mol–1. The average catalytic effect among the enzyme systems of 10.8 kcal mol–1 was 4.3 kcal mol–1 larger than that for the corresponding reaction in the GB pathway. Similarly to Δg‡, the values of Δg0 showed very little variation among the enzyme systems, averaging 8.4 kcal mol–1 in Polβ systems and 7.5 kcal mol–1 in Polλ systems. Thus, the free-energy difference between Δg0 and Δg‡ was, on average, nearly identical in the water (−4.6 kcal mol–1), Polβ (−4.4 kcal mol–1), and Polλ (−4.3 kcal mol–1) systems.

The IS1 → IS2b reaction in the SB pathway (triangles and pentagons in Figure S4; Table S38) proceeded with an average Δg‡ of 11.8 kcal mol–1 in Polβ systems and 11.6 kcal mol–1 in Polλ systems. Among the eight enzyme systems, the maximum difference in Δg‡ was only 2.3 kcal mol–1, and the standard deviation of the average Δg‡ of 11.7 kcal mol–1 was only 0.3 kcal mol–1; these parameters are nearly identical to those of the reaction involving IS2a. The average catalytic effect among the enzyme systems of 11.4 kcal mol–1 was 0.6 kcal mol–1 larger compared to that of the IS1 → IS2a reaction and 4.9 kcal mol–1 larger compared to that of the IS1 → IS2b reaction in the GB pathway. Similarly to Δg‡, the values of Δg0 also showed very little variation among the enzyme systems, averaging 5.7 kcal mol–1 in Polβ systems and 5.6 kcal mol–1 in Polλ systems. However, the free-energy difference between Δg0 and Δg‡ of −6.0 kcal mol–1 in the enzyme systems was 1.4 kcal mol–1 larger compared to that of the reference water system and 1.6 kcal mol–1 larger compared to that of the reaction involving IS2a. Thus, the IS1 → IS2b reaction was more favorable than the corresponding reaction involving IS2a in both Δg‡ and Δg0.

Part II of the P–O Bond Formation and Cleavage Reactions (IS2 → PS)

In the GB pathway, the IS2a → PS and IS2b → PS reactions proceeded with average Δg‡ values of 3.3 and 3.0 kcal mol–1 among Polβ and Polλ systems, that is, 0.7 and 0.4 kcal mol–1 above Δg‡ of the reference water system (Figure S5; Tables S39 and S40). The reactions were exergonic with average Δg0 values in Polβ and Polλ systems of −34.0 and −29.7 kcal mol–1 in the reactions involving IS2a and IS2b, respectively, that is, 1.2 kcal mol–1 below and 3.1 kcal mol–1 above Δg0 of the reference water system (−32.8 kcal mol–1). In the SB pathway, the IS2a → PS and IS2b → PS reactions proceeded with average Δg‡ values of 3.3 and 2.3 kcal mol–1 among Polβ and Polλ systems (Figure S6; Tables S45 and S46); these Δg‡ values were very similar to those of the GB pathway. The IS2a → PS and IS2b → PS reactions in the SB pathway were exergonic with an average Δg0 in Polβ and Polλ systems of −38.7 kcal mol–1 in the reactions involving IS2a or IS2b, that is, 9.0 kcal mol–1 below Δg0 of the reference water system (−29.7 kcal mol–1).

Cumulative Free-Energy Profiles of the Nucleotidyl-Transfer Reaction

The average activation and reaction free energies of the three successive steps, RS → IS1, IS1 → IS2, and IS2 → PS, were adjoined, and the cumulative Δg‡ and Δg0, relative to RS, were obtained for each step. The artificial IS2 in the water systems was removed by shifting the free energy of IS2 in both GB and SB pathways (IS2a by +3.17 and +3.25 kcal mol–1; IS2b by +3.16 and +3.27 kcal mol–1); these free-energy shifts were then applied to the corresponding IS2 in the enzyme systems (Tables and 4). The reaction coordinate of each state (Tables S47–S50) was obtained as the cumulative sum of energy gaps of the individual reaction steps (RS to PS in the GB pathway and IS1 to PS in the SB pathway) or by adoption from the GB pathway (RS → IS1 in the SB pathway).

Table 3

Cumulative Activation and Reaction Free Energies in the GB Pathway of the Nucleotidyl-Transfer Reaction in the Water, Polβ, and Polλ Systemsa

	TS1	IS1	TS2	IS2	TS3	PS
system	kcal mol–1
water	12.6	10.8	33.9	33.8	32.9	–2.5
Polβ1	10.2	5.2	23.1	23.8	24.1	–13.7
			22.3	19.7	19.5	–13.7
Polβ2	9.7	4.4	21.7	22.2	22.5	–14.0
			20.9	18.1	18.2	–13.3
Polβ3	8.5	3.2	25.2	28.0	27.4	–11.8
			23.3	23.1	22.4	–8.9
Polβ4	11.2	6.2	23.8	23.3	23.2	–12.2
			26.5	23.2	23.4	–9.7
Polλ1	12.4	6.8	23.2	23.7	23.5	–14.5
			23.7	21.1	20.6	–12.0
Polλ2	13.4	6.0	22.3	23.0	22.5	–15.1
			22.5	21.1	20.8	–12.9
Polλ3	11.5	5.9	20.9	21.0	19.8	–17.3
			22.0	19.3	18.1	–14.5
Polλ4	12.4	6.8	21.6	21.4	22.0	–15.1
			21.9	18.7	17.5	–16.6

TS1, transition state 1; IS1, intermediate state 1; TS2, transition state 2; IS2, intermediate state 2a or 2b (see Table ); TS3, transition state 3; PS, product state. Upper rows for the enzyme systems: IS1 → IS2a → PS; lower rows for the enzyme systems: IS1 → IS2b → PS. Numbers in bold denote the lowest rate-determining activation free energies among the GB (this table) and SB pathways (Table ).

Table 4

Cumulative Activation and Reaction Free Energies in the SB Pathway of the Nucleotidyl-Transfer Reaction in the Water, Polβ, and Polλ Systemsa

	IS1	TS2	IS2	TS3	PS
system	kcal mol–1
water	7.7	30.8	30.7	29.8	–2.5
Polβ1	1.6	14.6	13.7	13.5	–28.3
		13.6	11.4	10.2	–30.6
Polβ2	0.0	11.8	11.1	11.5	–30.3
		11.5	8.1	6.9	–33.3
Polβ3	2.2	15.2	13.9	13.8	–29.4
		13.8	11.4	9.9	–32.4
Polβ4	1.2	14.7	13.9	13.0	–29.0
		13.2	10.9	9.6	–31.3
Polλ1	7.4	19.3	19.4	19.8	–23.9
		18.5	16.3	16.2	–24.4
Polλ2	5.8	17.9	16.5	16.4	–24.0
		17.2	14.1	12.6	–27.8
Polλ3	11.6	23.1	22.8	21.9	–20.0
		22.4	20.5	19.0	–22.4
Polλ4	8.8	20.9	19.1	18.7	–22.3
		21.9	19.6	17.9	–22.7

IS1, intermediate state 1; TS2, transition state 2; IS2, intermediate state 2a or 2b (see Table ); TS3, transition state 3; PS, product state. Upper rows for the enzyme systems: IS1 → IS2a → PS; lower rows for the enzyme systems: IS1 → IS2b → PS. Numbers in bold denote the lowest rate-determining activation free energies among the SB (this table) and GB pathways (Table ).

Activation Free Energies of the Nucleotidyl-Transfer Reaction

The highest point in the free-energy profiles of the enzyme-catalyzed reactions (ΔG‡) was TS2 (as in the water system), IS2a, or TS3 (Figure S7). TS2 was the overall TS in 24 cases, including all Polβ and Polλ systems in the GB and SB pathways involving IS2b. IS2a in the GB pathway was the overall TS in Polβ3, Polλ1, Polλ2, and Polλ3. TS3 in the GB pathway involving IS2a was the overall TS in Polβ1, Polβ2, and Polλ4; TS3 in the SB pathway involving IS2a was the overall TS in Polλ1. The maximum difference in the free energy of TS2, IS2, and TS3 ranged from 0.5 to 4.4 kcal mol–1 in the GB pathway, and from 0.5 to 4.6 kcal mol–1 in the SB pathway.

No significant differences were observed in the values of ΔG‡ in the reactions involving IS2a and IS2b. In the GB pathway, the average ΔG‡ was 23.9, 23.2, and 21.7 kcal mol–1 in Polβ, tPolλ, and tPolλΔL1 systems. In the SB pathway, the average ΔG‡ was 13.6, 18.4, and 22.1 kcal mol–1 in Polβ, tPolλ, and tPolλΔL1 systems. The average catalytic effect in the GB pathway was 10.0, 10.7, and 12.2 kcal mol–1 in Polβ, tPolλ, and tPolλΔL1 systems. The average catalytic effect in the SB pathway was 17.2, 12.4, and 8.7 kcal mol–1 in Polβ, tPolλ, and tPolλΔL1 systems.

LFERs

LFERs with correlation coefficients >0.80 were detected for the free energies of TS3 versus PS in the GB pathway involving IS2b (R2 of 0.85) and for the free energies of IS1 versus PS, TS2 versus PS, and TS3 versus PS in the SB pathway involving IS2a or IS2b with R2 ranging from 0.94 to 0.98 (Figure ). In the SB pathway, the slope of the linear regression lines in IS1 versus PS, TS2 versus PS, and TS3 versus PS relationships for the reactions involving IS2a and IS2b ranged from 0.92 to 1.08. Averaging over the compact (via IS2a) and loose (via IS2b) variants of the SB pathway yielded the following LFER equations for relating the free-energy changes in IS1, TS2, TS3, and PS

Figure 3

Linear free energy relationships (LFERs, kcal/mol) among IS1, TS2, TS3, and PS in the GB and SB pathways of the nucleotidyl-transfer reaction in Polβ and Polλ systems. Gray, Polβ1; red, Polβ2; blue, Polβ3; green, Polβ4; cyan, Polλ1; magenta, Polλ2; yellow, Polλ3; orange, Polλ4. Squares and triangles, the reactions involving IS2a; circles and pentagons, the reactions involving IS2b. Solid lines, linear regressions for the reactions involving IS2a; dashed lines, linear regressions for the reactions involving IS2b; correlation coefficients (R2) of the linear regressions are given in parentheses (solid line; dashed line).

Comparison of the Calculated Energetics with the Available Experimental and Theoretical Data

The lowest ΔG‡ values for Polβ were obtained in the SB pathway: 12.5 kcal mol–1 for dCTP·dG and 13.5 kcal mol–1 for dGTP·dC nascent base pairs, which was 2.9–4.3 kcal mol–1 below the experimental and 1.5–5.6 kcal mol–1 below the theoretical values obtained from the literature (Table S51). The lowest ΔG‡ values for tPolλ were obtained in the SB pathway: 17.9 kcal mol–1 for the dCTP·dG nascent base pair, which was 0.4 kcal mol–1 above the experimental and 0.3–0.5 kcal mol–1 above the theoretical (albeit for a different nascent base pair) values obtained from the literature (Table S52). The lowest ΔG‡ values for tPolλΔL1 were 21.0 kcal mol–1 for dGTP·dC in the GB pathway and 20.9 kcal mol–1 for dGTP·dT in the SB pathway; the latter was 1.7 kcal mol–1 above the experimental value for tPolλ obtained from the literature and represented the first reported prediction for tPolλΔL1 for which no experimental pre-steady-state kinetic data exist. The values of ΔG0 for Polβ and Polλ were 12.8–30.9 kcal mol–1 below the corresponding values that had been reported in other theoretical studies.

Geometries of the RS, IS, TS, and PS of the Nucleotidyl-Transfer Reaction

Animations of the nucleotidyl-transfer reaction generated from the representative structures of each state showed similar atomic motions in the active site for all enzyme systems (Supporting Files S3). The geometric determinants of ΔG‡ and ΔG0 were searched for in the structures of IS1, TS2, TS3, and PS.

TS Structures

Pα–Olg and Onuc–Pα distances in TS2 and TS3 differed in the pathways involving IS2a and IS2b (Figure ; Tables S53–56). TS2a and TS3a featured shorter Pα–Olg (by 0.11 and 0.17 Å on average) and Onuc–Pα (by 0.09 and 0.12 Å on average) distances than TS2b and TS3b among all enzyme systems in the GB and SB pathways. Pα–Olg distances in TS2 and TS3 were only marginally shorter (by 0.01 and 0.03 Å on average) and Onuc–Pα distances were longer (by 0.06 and 0.03 Å on average) in the SB than the GB pathway. The water systems showed no difference in Pα–Olg distances in TS2 between the GB and SB pathways but showed longer Onuc–Pα in TS2 of the SB (than the GB) pathway (by 0.13 Å). Compared to water systems, enzymes showed: (1) similar Pα–Olg distance (shorter by 0.02 Å on average); (2) longer Onuc–Pα distance in TS2 (by 0.12 Å on average) and TS3 (by 0.06 Å on average) in the GB and SB pathways, with TS2a of Polβ3 in the GB pathway being an outlier, having water-like Onuc–Pα distance and also the highest ΔG‡ (28.0 kcal mol–1) among all enzyme systems.

Figure 4

P–O distances (Å) in TS2 (lower left) and TS3 (upper right) in the GB (squares and circles) and SB (triangles and pentagons) pathways of the nucleotidyl-transfer reaction in the water, Polβ, and Polλ systems. Squares and triangles, TSs in the reactions involving IS2a; circles and pentagons, TSs in the reactions involving IS2b.

TS2 structures in all enzyme systems were stabilized by a network of hydrogen-bonding interactions (Figures and S8–S15). The Watson–Crick hydrogen-bonding interactions between complementary nucleobases were always present, as were the interactions between the O3β atom of dNTP and an amino group of Arg183/420/420 in Polβ/tPolλ/tPolλΔL1 and between the hydroxyl group of Tyr271/505/500 in Polβ/tPolλ/tPolλΔL1 and the carbonyl-group oxygen of the 3′-terminal dC of the primer DNA strand. The arginine could also interact with the hydroxyl group of dNTP, which could also form alternative interactions with the O1β atom of dNTP and with the carbonyl-group oxygen of the tyrosine which could simultaneously form a hydrogen bond with the amino group of Asn279/513/508 in Polβ/tPolλ/tPolλΔL1. The asparagine could also be involved in alternative hydrogen-bonding interactions with the nucleobases of the dNTP·dN nascent base pair.

Figure 5

Representative TS2 structures in Polβ, tPolλ, and tPolλΔL1 systems. Polβ, Polβ4 in SB2b pathway; tPolλ, Polλ2 in SB2b pathway; tPolλΔL1, Polλ4 in SB2b pathway. Atoms are colored by element: green, carbon in amino acid residues; gray, carbon in nucleotides; red, oxygen; blue, nitrogen; orange, phosphorus; white, hydrogen. Hydrogens bound to carbon atoms are omitted for clarity. Dotted lines indicate hydrogen-bonding interactions and the partial Onuc–Pα bond. See Figures S4–S11 for the full set of representative TS2 structures.

PS Structures

PS structures of all enzyme systems showed both similarities to and differences from the 4BU3 X-ray crystal structure of the Polβ product complex (Figures S16–S23). The similarities in the atomic positions included: (1) the nucleotides; (2) the PPi product; (3) the MgB2+ ion coordinating Asp190/427/427 and Asp192/429/429 in Polβ/tPolλ/tPolλΔL1 with PPi; and (4) Arg183/420/420 in Polβ/tPolλ/tPolλΔL1 interacting with the O1β atom of the inserted nucleotide. The differences in atomic positions included: (1) closer interaction between the amino groups of Arg149/386/386 and oxygen atoms of PPi in the PS structures; and (2) the MgA2+ ion coordinating the Onuc, Asp256/490/485, Asp190/427/427, and Asp192/429/429 in the PS structures of Polβ/tPolλ/tPolλΔL1 versus coordinating Olg of PPi and a phosphate-group oxygen of the inserted dNMP in the X-ray crystal structure; in the PS structures, the new crystallographic position of the MgA2+ ion was occupied by a sodium ion (Polλ1 in the SB2 pathway involving IS2b) or (in all other cases) by water molecules.

Cumulative Charge Profiles

The simple geometrical analysis failed to explain the observed LFERs. Thus, this phenomenon had to be due to a complex combination of short- and long-range electrostatic interactions. But an attempt to quantify these interactions, as the charge distribution around Onuc in IS1 and around the O3β atom of PPi in PS, did not provide any hint about the source of the LFER (Figure S24): No obvious link between the cumulative charge profiles and Δg0 could be established.

Discussion

X-family DNA polymerases β and λ are involved in DNA damage repair and DNA recombination. In vitro and in vivo knockout experiments have demonstrated their vital roles: Polβ and Polλ double null mouse embryonic fibroblasts are hypersensitive to alkylating and oxidizing DNA-damaging agents;[15] Polβ knockout mice exhibit apoptosis in post-mitotic neuronal cells and die at birth;[79] heterozygous Polβ+/– mice develop normally, but have increased incidence of lymphomas;[80] and Polλ knockout mice are viable and fertile but display reduced immunoglobulin heavy-chain junction variability due to impaired nonhomologous end-joining.[81] Polβ and Polλ protect the integrity of the genome by being the principal nucleotidyltransferases in single-nucleotide gap-filling DNA synthesis during BER;[15] Polλ also contributes to generating the diverse repertoire of immunoglobulins and T-cell receptors via the gap-filling activity during V(D)J recombination.[81] When overexpressed in tumors, Polβ and Polλ can act as oncogenes by increasing the mutation rate and enhancing the capability of DNA repair.[82] Polβ and Polλ play their biological roles by dRP lyase and nucleotidyltransferase activities; we explore the latter.

Polβ and Polλ share the common mechanism of the nucleotidyl-transfer reaction:[19]In this study, we modeled the chemical transformation of RS to PS (reaction II in eq ) in Polβ, truncated Polλ (tPolλ), and Loop1 mutant of tPolλ (tPolλΔL1) using FEP and EVB simulations in three stepswhere IS1 and IS2 are the intermediate states, RS → IS1 is the proton transfer step (DNA-3′-OH → DNA-Onuc− + H+), and IS1 → IS2 → PS comprises the Onuc–Pα bond formation and the Pα–Olg bond cleavage reactions. Recent ab initio QM calculations of a model reaction between tetrahydrofuranol and methyl triphosphate have suggested that IS2 does not exist,[69] which means that the reaction scheme of the uncatalyzed reaction between 2′-deoxyribose (dRib) and dNTP can be simplified intoWe removed IS2 in the FEP/EVB free-energy profiles by increasing the reaction free energy (Δg0) of the IS1 → IS2 step by about 3 kcal mol–1, which resulted in the free energy of IS2 being almost equal to TS2 and TS3. When we applied the same free-energy shift to IS2 in the free-energy profiles of the reactions catalyzed by Polβ and Polλ, we also obtained a flat free-energy surface from TS2 to TS3The absence of IS2 in the water, Polβ, and Polλ systems means that, using the terminology offered by Lassila et al.,[83] the Onuc–Pα bond formation and the Pα–Olg bond cleavage reactions are concerted,[30,44,48,49,75,84] and, unlike T7 DNA polymerase,[57] the studied enzymes do not change the concertedness of the reactions. The broad, flat TS region of the IS1 → PS reaction in the water and enzyme systems accommodates a range of associative, pentavalent-intermediate-like geometries: from the “early” TS (with less Onuc–Pα than Pα–Olg bonding) to the “late” TS (with more Onuc–Pα than Pα–Olg bonding), and with the total bonding at the Pα atom between 1.2 (compact TS) and 0.9 (loose TS) relative to RS.

IS1, the proton transfer intermediate, might also be artificial as it is absent in QM/MM calculations of the nucleotidyl-transfer reaction catalyzed by tPolλ,[49] which (if true) simplifies the reaction mechanism into the most simple form possibleBut if IS1 does not exist, then the rate-limiting step includes all chemical transformations; can the proton transfer to Asp490 be distinguished as the rate-limiting step in such a scenario?[49] (No, it cannot.) We assumed that IS1 does exist (eq ),[43−48] which allowed us to view the proton transfer as a distinct step. We could thus explore various proton transfer mechanisms and assess their relative feasibility from the calculated activation (Δg‡) and reaction free energies (Δg0).

In the RS → IS1 reaction, the 3′-OH group of the 3′-terminal nucleotide of the primer DNA strand loses its proton; the oxygen becomes Onuc, poised for the nucleophilic attack on the Pα atom of dNTP. The proton transfer is facilitated by the catalytic magnesium ion (MgA2+) and by the availability of a suitable proton acceptor. What are the requirements for the initial proton acceptor? It must (1) be thermodynamically “willing” to accept the proton, (2) be located in the vicinity of the 3′-OH group (or be able to migrate there), and (3) its protonation should not impede the subsequent reaction steps. In Polβ, tPolλ, and tPolλΔL1, there are six candidates that might satisfy these conditions (see Figure ): (1) oxygen of a water molecule migrating to the active site from the bulk solution,[30,46,48,51] (2) the distal carboxylic oxygen of Asp256/490/485 (Oδ1), (3) the proximal carboxylic oxygen of Asp256/490/485 (Oδ2),[43,45,47,49,51] (4) the proximal carboxylic oxygen of Asp192/429/429 (Oδ2),[45,49] (5) oxygen of the water molecule that coordinates the MgA2+,[44,45,49] and (6) a nonbridging oxygen of the Pα-group of dNTP (Oα2).[42,45,49]

Two QM/MM studies compared multiple proton transfer pathways: In Polβ, Alberts et al.[45] found the e → c proton transfer to be more feasible than the direct proton transfer to acceptors c, d, or f; in tPolλ, Cisneros et al.[49] found the proton transfer to acceptor c to be more feasible than the proton transfer to acceptors d, e, or f. Neither of these two studies, however, assessed the possibility of the proton transfer to the bulk solution (acceptor a), which requires calculating pKa of the 3′-OH group. The proton transfer to acceptor a has so far been reported only for Polβ, using FEP or protein-dipoles-Langevin-dipoles linear-response approximation (PDLD/S-LRA) simulations, by Sucato et al.,[30] Xiang et al.,[46] Ram Prasad and Warshel,[48] and Matute et al.:[51] acceptor a (Δg0 of 4 kcal mol–1) seems to be at least as good as c (Δg0 of 6 kcal mol–1).[51] How about the remaining proton acceptors? We did not find any report about the feasibility of acceptor b. Acceptor e was an intermediate for the transfer to d(44) or c(45) in QM/MM calculations with Polβ, or scored badly (ΔE‡ ∼ 35 kcal mol–1) in QM/MM calculations with tPolλ.[49] And acceptor f was first suggested for Polβ, based on pioneering QM calculations, by Abashkin et al.,[42] but the proton transfer to acceptor f failed to provide reasonable activation potential energies in two subsequent QM/MM studies on Polβ[45] and tPolλ.[49] We therefore excluded e and f from our consideration and modeled the proton transfer to the four remaining acceptors: the SB proton transfer to the bulk solution (a) using FEP simulations, and the GB proton transfer to b, c, and d using EVB simulations.

We found the acceptor d unlikely (Δg‡ > 24 kcal mol–1) because of the long distance between Onuc and Oδ2 of Asp192/429/429 in Polβ/tPolλ/tPolλΔL1 and because of steric hindrance and electrostatic repulsion by MgA2+. The proton transfer to acceptor b proceeded with better energetics (Δg‡ > 16 kcal mol–1) but required a conformational change of Asp256/490/485 in Polβ/tPolλ/tPolλΔL1 for reducing the distance between the Onuc and the Oδ1 atom of the aspartate; the conformational change perturbed the interaction of the aspartate with MgA2+. The proton transfer to acceptor c proceeded smoothly, over the shortest distance, and without disturbing the interaction between Oδ2 of the aspartate and MgA2+, which contributed to a better energetics: Δg‡ (Δg0) was 10 (5), 13 (6), and 12 (6) kcal mol–1 in Polβ, tPolλ, and tPolλΔL1, respectively. But we obtained similar, or even more favorable, energetics for the proton transfer to the bulk solution. This mechanism is structurally supported by the existence of the half-open active site, replete with water molecules, which form a continuous chain all the way to the bulk solution in crystal structures of ternary complexes in the so-called “closed conformation”. Proton transfer reactions along such proton-acceptor chains are very fast,[71] which means that the corresponding activation free-energy barriers are primarily determined by their reaction free energy. If we assume a small <2 kcal mol–1 free-energy difference between Δg‡ and Δg0 for the proton transfer to the bulk solution at pH 7, this pathway yields Δg‡ (Δg0) of 3 (1), 9 (7), and 12 (10) kcal mol–1 in Polβ, tPolλ, and tPolλΔL1, respectively. Thus, the proton transfer to the bulk solution (via Asp256/490/485, an active site water molecule, or a hydroxide anion) seems to be the proton transfer pathway of the least resistance.[51] Moreover, the subsequent IS1 → IS2 reaction step, which also contributes into the overall free-energy barrier, proceeds with lower Δg‡ when preceded by the SB proton transfer. Two questions remain unresolved: (1) What makes the proton transfer to the bulk solution superior over the proton transfer to the active site aspartate in Polβ and tPolλ, but not in tPolλΔL1? (2) What is the source of the large differences in Δg0 of the proton transfer to the bulk solution among the three enzyme systems?

Although the seemingly simple proton transfer to the bulk water yielded large differences between the studied systems, the subsequent concerted Onuc–Pα bond formation and Pα–Olg bond cleavage reactions (IS1 → PS), modeled using three VB states and involving large changes in atomic positions and charge distributions, proceeded uniformly in Polβ, tPolλ, and tPolλΔL1, regardless of the initial coordinates (X-ray crystal structures: 2FMP, 2PXI, 2PFO, 2PFP, or 3PML) and the nascent base pair (dCTP·dG, dGTP·dC, or dGTP·dT): Δg‡ of 12 kcal mol–1 and Δg0 of −32 kcal mol–1 for the reactions involving compact or loose TS geometries. The activation free-energy barrier of this step is 9, 6, and 0 kcal mol–1 larger than the assumed Δg‡ of the proton transfer in Polβ, tPolλ, and tPolλΔL1. Thus, the P–O bond formation/cleavage reaction is the rate-limiting step in Polβ (in agreement with refs (30, 43, 46, 48, 51) and in disagreement with ref (44)) and tPolλ (in disagreement with ref (49)), and rate-colimiting (with the proton transfer) in tPolλΔL1where k1 is the rate constant of the RS → IS1 reaction, k–1 is the rate constant of the RS ← IS1 reaction, and k2 is the rate constant of the IS1 → PS reaction. Because k2 is a constant in Polβ, tPolλ, and tPolλΔL1, the relative reaction rates among these enzymes depend solely on the pKa of the OnucEquation describes an intermolecular LFER (not to be confused with the intramolecular Brønsted LFER) between Δg0 of the proton transfer step and ΔG‡ of the overall reaction with a slope of ∼1. This idea can be tested experimentally using nonhydrolyzable dNTP analogues with the Olg atom replaced by a CH2, CF2, or NH group, which will trap the system in the RS ⇌ IS1 reaction and thus allow the evaluation of the Keq = k1/k–1 term in eq . If the values of Keq for the dNTP analogues correlate with the values of catalytic rate constants kpol, measured for the corresponding native dNTPs, the proposed intermolecular LFER will be validated. The values of ΔΔG‡ could then be calculated using linear-response simulations on the reactant and either TS[63,85] or IS1 structures.

The uniformly depressed Δg‡ of the rate-limiting P–O bond formation/cleavage step is also the main contribution (−11 kcal mol–1) to the average catalytic effect of −17, −12, and −9 kcal mol–1 in Polβ, tPolλ, and tPolλΔL1, which we obtained by comparing ΔG‡ for the enzymes and for the reference reaction between dRib and dNTP in solution. This catalytic contribution appears to be very robust with respect to small perturbations of the protein structure and even with respect to mispairing between the templating and incoming nucleobases.[86] When the closed, “catalytically competent”, ternary complexes of DNA polymerases are used for the simulations, the fidelity control originates from the preceding steps: thermodynamics of dNTP binding[44,86−88] or the deprotonation kinetics of the 3′-OH group of the primer DNA strand[44] because these steps are more sensitive to “small” structural variations in the ternary complexes. When the enzyme carries out the reaction while in the (partially-)open, “catalytically incompetent”, conformation,[89,90] the catalysis of the P–O bond formation/breakage reactions is greatly diminished.[46,85]

In search for the main source of the catalytic effect of Polβ, tPolλ, and tPolλΔL1, we explored the representative structures of the RS, IS, TS, and PS and asked what interactions are responsible for the higher stability of the TS in the enzyme compared to that in the solution. The characteristic feature of the active sites of Polβ, tPolλ, and tPolλΔL1 is a network of favorable electrostatic (including hydrogen-bonding) interactions, which interconnects the three apparent binding sites for the nucleobase, dRib, and triphosphate groups of dNTP.[40,63,90,91] These interactions, present already in the RS and preserved during the entire nucleotidyl-transfer reaction, involve the sidechains of Arg149/386/386, Arg183/420/420, and Asn279/513/508, the backbone and sidechain of Tyr271/505/500, and the nucleobase of the templating nucleotide in Polβ/tPolλ/tPolλΔL1. In the solution reaction, the role of the amino acid residues and the templating nucleotide is fulfilled only suboptimally by water molecules, which are incapable of establishing an electrostatically preorganized reaction cage.[32,87,99]

Assuming the validity of the TS theory (eq ), we can make a direct comparison between the calculated and available experimental ΔG‡ of the nucleotidyl-transfer reaction catalyzed by Polβ and tPolλ: Our FEP/EVB simulations overestimate the catalytic effect of Polβ by 3–4 kcal mol–1[30,31,41,92−94] but give an accurate prediction of the catalytic effect of tPolλ.[29] Note that all our calculated ΔG‡ values are genuine predictions because the only experimental inputs were the X-ray crystal structures of ternary complexes of Polβ, tPolλ, and tPolλΔL1, pKa of dRib, and pH of maximum kpol in Polβ. In this light, the low single digit difference in ΔG‡ might seem small except that kpol scales exponentially with ΔG‡: the difference of 4 kcal mol–1 in ΔG‡ translates into >300-fold difference in kpol. In addition, our ΔG‡ values carry a cumulative uncertainty of 6 kcal mol–1. The standard deviations of individual reaction steps that add up to this total uncertainty have been assessed by simulating the same reaction in the forward and reverse directions. Thus, the activation and reaction free energies for each reaction step, each reaction mechanism, and each enzyme were determined by at least four (and up to 16) separate 510 ps FEP or FEP/EVB simulations. Of course, larger number of longer simulations (at least by 1 order of magnitude) might be a way to improve the conformational sampling along the reaction coordinate and thus, a way to increase the accuracy and reduce the uncertainty of ΔG‡ values. However, in studies aimed at exploring, quantitatively, complex enzymatic reaction mechanisms, the strategy of running multiple short simulations and averaging over very different simulation setups[95] is more cost-effective than running fewer but longer simulations for the same total simulation time. Previous QM/MM and FEP/EVB predictions for Polβ were more accurate than ours, for unknown reasons, deviating only by 1 to 2 kcal mol–1 from the experimental ΔG‡.[30,44,46] The only computational prediction for tPolλ was reported by Cisneros et al.,[49] who obtained an activation potential energy of 17.6 kcal mol–1 using QM/MM calculations, that is, in perfect agreement with our FEP/EVB activation free energy of 17.9 kcal mol–1.

The FEP/EVB approach, which was developed more than 35 years ago, remains underappreciated and underused by the molecular modeling community at large, despite multiple case studies demonstrating its robustness and versatility in modeling of enzyme catalysis. No molecular modeling tool can stand on its own and FEP/EVB is not an exception. But when complemented by QM and QM/MM calculations, linear-response simulations and, most importantly, by experiments (X-ray crystallography, NMR, pre-steady-state kinetics, kinetic isotope effect, Brønsted LFER, and our proposed intermolecular LFER), the FEP/EVB approach can provide a unique insight into the structure–dynamics–function–evolution relationships of DNA polymerases β and λ and into their (ab)normal roles in DNA repair and recombination.

Conclusions

Assuming our computational model is a correct imitation of reality, the three introductory questions can be answered as follows: (1) What is the preferred mechanism for the deprotonation of the 3′-OH group of the primer DNA strand? The preferred mechanism for the deprotonation of the 3′-OH group of the primer DNA strand is the transfer of the proton to the bulk aqueous solution via a multitude of possible pathways, involving active site amino acid residues (including, but not limited to, Asp256/490/485 in Polβ/tPolλ/tPolλΔL1), hydroxide anions accessing the active site, and/or chains of water molecules connecting the active site to the bulk solution. (2) What is the degree of concertedness of the P–O bond formation and cleavage reactions? The Onuc–Pα bond formation and Pα–Olg bond cleavage reactions are concerted, and proceed via an early, loose TS. (3) Which of the three reaction steps limits the overall reaction rate? Only two steps can be distinguished: proton transfer followed by Onuc–Pα bond formation and Pα–Olg bond cleavage; and whereas both these steps contribute to the overall activation barrier of the nucleotidyl-transfer reaction, the P–O bond formation/cleavage represents the rate-limiting step. The overall activation barrier (and hence kpol), which can be perturbed by small alterations in the structures of the catalytically competent ternary complexes, is linearly proportional to the stability of the deprotonated nucleophilic oxygen, that is, to the pKa of the 3′-OH group of the primer DNA strand in the enzyme active site. The observed intermolecular LFER can be exploited, computationally or experimentally, for predicting the catalytic effects of nondeleterious mutations in DNA polymerases.