Preferred WMSA catalytic mechanism of the nucleotidyl transfer reaction in human DNA polymerase κ elucidates error-free bypass of a bulky DNA lesion

Structure preparation

The ternary crystal structure (PDB ID (46): 2OH2) of Pol κ_19-526 with DNA and incoming dTTP (32) was the basis for our prepared enzyme-substrate model. Details are given in Supplementary Data including Supplementary Scheme S1.

Classical MD

The initial model was subjected to 3 ns of MD simulations using AMBER 9 (47). We employed the Amber99SB (48–50) force field with modification for DNA by parmbsc0 (51), and previously calculated parameters for the dCTP (52). The structure was neutralized by 10 Na⁺ counterions and was solvated with a periodic rectangular box of TIP3P water (53,54) with 10 Å buffer around the enzyme-substrate complex. Details of the MD protocols are given in Supplementary Materials. A random snapshot from the last 0.5 ns of the 3 ns trajectory was utilized for the subsequent QM/MM–MD stage. In this snapshot the Mg⁺² coordination, Mg⁺²–Mg⁺² distance, and O3′–P_α distance were similar to those in the high resolution X-ray crystal structure of pol β (34), in our previous Pol κ model (30), and in the initial models utilized for our previous polymerase mechanism studies with Dpo4 (35) and T7 (36) (shown in Figure 2 and Supplementary Figure S2).

QM/MM and QM/MM–MD

The QM subsystem (81 atoms) was comprised of the dCTP, the primer 3′-nucleotide (excluding the C5′ phosphate group), the glutamic acid (E199), the catalytic and nucleotide-binding Mg⁺² ions (Mg⁺² A and Mg⁺² B, respectively), and two water molecules (Supplementary Figure S3). While the QM subsystem did not include all the ligands of the two Mg⁺² ions (Figure 2), the critical ones for the chemistry of the nucleotidyl transfer reaction are included: E199, the crystal water molecule, the attacking O3′ and the dNTP phosphate oxygens. The other ligands, which do not participate in the chemistry per se but act to maintain the active site structure in a reaction-ready organization, are simulated through the MM component of our computations. Hence, the overall treatment of the Mg⁺² ion ligands is sufficiently robust. This subsystem was treated by the hybrid density functional B3LYP (55–57) with a medium split valence basis set and polarization functions 6–31G* (58,59). All other residues and the surrounding water molecules (11 774 atoms), constituting the MM subsystem, were described by the TIP3P water model and the modified Amber99SB force field that was used in the classical MD simulations. Spherical boundary conditions were applied in all QM/MM calculations and all atoms outside the 20 Å radius centered at the primer 3′-end oxygen atom were frozen. A cutoff of 12.0 Å was used for van der Waals interactions, and a cutoff of 18.0 Å was employed for electrostatic interactions among the MM atoms. There was no cutoff for electrostatic interactions between the QM and MM atoms. Solvent molecules outside of a 27 Å sphere centered at the primer 3′-end oxygen atom were removed. The pseudobond approach (40,41) with improved parameters (60) was employed to treat the partition between QM and MM subsystems (Supplementary Figure S3). All QM/MM calculations were carried out utilizing the modified QChem (61) and Tinker (62) programs. The reactant was minimized by the efficient iterative optimization QM/MM approach. We mapped the minimum energy path for the investigated mechanisms utilizing the reaction coordinate driving method (42). For each determined structure along the path, a 500 ps MD simulation with the MM force field was carried out to equilibrate the MM sub-system. The resulting snapshot was used as the starting structure for ab initio QM/MM–MD simulations with umbrella sampling (38,39,63) that applies a harmonic potential to constrain the reaction coordinate (RC) at successive values. In order to ensure sufficient overlap between the successive windows, force constants in the range of 0 to 300 kcal/mol-Ų were employed; these were chosen based on the estimated slope of the free energy curve at the reaction coordinate. The potential of mean force (PMF) was calculated from the probability distributions along a reaction coordinate using WHAM (43–45).

Pol κ containing the B[a]P-dG lesion in the pentacovalent phosphorane transition state

For the studies involving the B[a]P-dG adduct we began with the pentacovalent phosphorane transition state structure of the characterized WMSA mechanism. We modeled the B[a]P-dG lesion by docking the B[a]P moiety to the amino group of the templating guanine, placing the B[a]P on the minor groove side in Pol κ, oriented in the 5′-direction of the modified strand as in the NMR solution structure (26). Because the transition state contains five covalent bonds about the α-phosphate, it was necessary to develop an Amber-compatible force field to simulate this transient chemical state for the unmodified Pol κ. A new residue comprised of the incoming dCTP and the 3′-terminal nucleotide of the primer DNA strand was defined. This residue contained four new atom types: one for the protonated γ-phosphate oxygen of the dCTP, a second for the γ protonating hydrogen, a third for the α-phosphate of the dCTP and the fourth atom type was defined to describe both the dCTP α,β-bridging oxygen and the primer O3′. Topology assignments, geometries and partial charge parameters were characterized on the basis of the QM/MM–MD calculations for the transition state (data are listed in Supplementary Tables S1 and S2). The force constants were chosen to be large enough to maintain the transition state bond lengths, bond angles and dihedral angles in the unmodified Pol κ. Partial charges were taken from the QM calculations. The force field parameters for the guanine lesion, B[a]P-dG, were taken from Jia and coworkers (30). The MD protocol was the same as described in Supplementary Data.

All initial model constructions were carried out using the INSIGHT II 2005 program (Accelrys Software, Inc.). PyMOL (Delano Scientific, LLC) (64) was employed to make molecular images and movies.

In the WMSA mechanism the γ-phosphate of the incoming dNTP serves as the base for the H(O3′) which shuttles through water molecules

We found that a variant of the WMSA mechanism, which was previously proposed for the bypass DNA polymerase Dpo4 (35) and the high-fidelity replicative T7 DNA polymerase (36), is preferred for the human bypass DNA polymerase κ over the many explored schemes. Two steps are involved in this mechanism. In the first step, the H of O3′ is transferred to the oxygen on the γ-phosphate of the dCTP via two mediating water molecules (Figures 3 and 4). Hydrogen bond analysis revealed stable interactions throughout the reaction pathway between the β-phosphate group of the dCTP with Thr138, Arg144 and Phe111, while the γ-phosphate group hydrogen bonds with Arg144 and Lys328 and with the Ala110 backbone (Figure 5 and Supplementary Figure S7). The γ-phosphate has the most negative partial charge among the phosphate groups of dCTP (52), and this is relieved by protonation (Figure 6). In addition, the α-phosphate maintained stable hydrogen bonds with the two water molecules and helped anchor the second water (not coordinated with Mg⁺²) throughout the reaction. The proton transfers in this shuttle are: from the O3′ to the crystal water molecule that participates in the Mg⁺² coordination; from this water to a solvent water, that is stably held by hydrogen bonds with oxygen atoms of the substrate, and from the second water to the γ-phosphate of the dCTP (Figures 3A and 4). The crystal water found in Pol κ (32) is also observed in a number of other polymerases (34,65–70). In the organized active site of Pol β (34) and in several other DNA polymerases (66,68–70) a second crystal water is located at the same position as the solvent water (WAT2) in our model (Supplementary Figure S1). A free energy profile was calculated for the first step of this mechanism by employing 37 umbrella windows along the reaction coordinate, each simulated for 30 ps of B3LYP(6–31G*) QM/MM–MD. An activation barrier of 11.2 ± 0.4 kcal/mol was determined for this free energy profile (Figure 3C). (Supplementary Movie S1).

In an associative mechanism a pentacovalent phosphorane is the transition state of the nucleotidyl transfer reaction

With a protonated γ-phosphate as a stable intermediate the reaction proceeds to the second step with a nucleophilic attack of the naked O3′ on the α-phosphate and leaving of the pyrophosphate (Figures 3B and 4). The transition state is a pentacovalent phosphorane structure, which results from an associative nucleotidyl transfer reaction; that is, the O3′–P_α bond forms before the P_α–O_α,β-bridge breaks (Supplementary Figure S8). The free energy barrier obtained for this step is 10.6 ± 0.3 kcal/mol, calculated using 50 umbrella sampling windows along the reaction coordinate each for 35 ps of B3LYP(6–31G*) QM/MM–MD (Figure 3C) (Supplementary Movie S2). As the reaction progresses the distance between the catalytic Mg⁺² A and the nucleotide-binding Mg⁺² B becomes shorter and achieves a minimum of 3.59 ± 0.09 Å at the transition state (Supplementary Figure S2A) which promotes the formation of the O3′–P_α bond; it subsequently elongates again to facilitate the breaking of the P_α–O_α,β-bridge bond. By approaching metal B, metal A brings the nucleophile O3′ within attack distance of the P_α, as was suggested by Yang et al. (71).

Stabilization of the pentacovalent phosphorane transition state involves interactions with nearby amino acid residues

We wished to elucidate the roles of the amino acid residues in the active site in facilitating the WMSA mechanism. We analyzed hydrogen bonding interactions, and calculated intermolecular interactions (including van der Waals and electrostatic energies). We also analyzed group partial charges of each of the three phosphates in the dCTP along the reaction pathway for the nucleophilic attack. We compared the transition state values with those for the intermediate in which the γ-phosphate was protonated but phosphoryl transfer had not been initiated, in order to reveal the specific interactions that stabilize the transition state.

Key hydrogen bonding interactions present in the transition state are shown in Figure 5A. The partial charge analysis showed that negative charge builds up on the α and also to some extent on the β phosphates as the pentacovalent phosphorane forms during the nucleophilic attack (Figure 6). This motivates the proton on the γ-phosphate to migrate toward the α,β-bridge. Consequently the γ-phosphate group becomes highly negatively charged (Figure 6). This negative charge provokes stronger electrostatic interactions between Arg144, Tyr141 and Lys328 and the γ-phosphate group compared to their interactions in the γ-protonated intermediate (Figure 7, Supplementary Figures S9 and S10). In addition there is a water mediated hydrogen bond involving Tyr141 with the γ-phosphate throughout both steps of the reaction. In the transition state the hydrogen bond is direct without water mediation (Supplementary Figure S11).

Transition state alignment is facilitated by C3′-endo sugar pucker

The conformation of the sugar pucker of the primer terminus through the progression of the reaction is shown in Supplementary Figure S12. In the original crystal structure, which had a hydrogen in place of the 3′OH, the sugar pucker was C2′-endo, but upon remodeling and MD simulation (see ‘Computational Methods’ section) the pucker changed to C3′-endo. Supplementary Figure S12 shows the pseudorotation parameter P (72,73) during the QM/MM–MD as a function of the reaction coordinate. We see that the pucker remains near C3′-endo until the neighborhood of the intermediate (Figure 4), where it oscillates around C1′-exo but returns to C3′-endo in the transition state, providing the proper alignment for the transition state geometry. Crystal structures of pre-reactant polymerase ternary complexes have been predominantly C3′-endo (34,66,68–70,74–76) except for Y-Family polymerases (2).

Spontaneous proton transfer from the γ-phosphate to the α,β-bridging oxygen occurs during leaving of the pyrophosphate

Under the influence of the QM/MM–MD the dynamics revealed spontaneous proton transfer from the γ-phosphate to the α,β phosphate bridge during the leaving of the pyrophosphate (Figure 4). This protonation accompanied the nucleotidyl transfer reaction step after formation of the pentacovalent phosphorane transition state and before the pyrophosphate left completely. The leaving of the pyrophosphate is assisted by this proton transfer, which relieves the accumulating negative charge on the α,β-bridging oxygen during the elongation of the P_α–O_α,β-bridge bond (Figure 6).

Other mechanisms provide unstable intermediates or higher energy barriers

Other pathways that we explored for the deprotonation of the H(O3′) included: (i) Amino acids as the general base for H(O3 ′). There are two carboxylate-containing amino acids within a short distance of the primer (E199 and D198). Both are part of the conserved amino acid triad in polymerase active sites (4). We considered the closer one, E199 as a natural candidate for the base that absorbs the O3′ proton. While the potential energy for this transfer revealed a promising activation barrier of 30.1 kcal/mol, the E199 protonated intermediate was not stable and the hydrogen quickly returned to the O3′ upon an unrestrained QM/MM–MD of 0.5ps (Supplementary Figure S13). This contrasts with the case for the WMSA mechanism discussed above where the O3′ proton remained stably on the γ-phosphate during 35 ps of unrestrained QM/MM–MD. We next tried a one-step mechanism in which the O3′ donates its proton to E199 while it simultaneously approaches the α-phosphate for the nucleophilic attack. This mechanism always produced an uphill reaction energy profile (Supplementary Figure S14).

(ii) Amino Acids as donor to protonate the dNTP. Experimental evidence indicates that the incoming dNTP must become protonated to ease the phosphoryl transfer and the leaving of the pyrophosphate, and it has been suggested that an amino acid could be the donor in the case of certain other polymerases (77). In Pol κ a suitable donor amino acid could be Arg144 because a guanidino group proton is 2.8 Å from the γ phosphate oxygen (Figure 5 and Supplementary Figure S7). Therefore, we explored a path that involves protonation of the γ-phosphate by this Arg144 (Supplementary Figure S15A). Another path we explored for the γ-phosphate protonation combined proton donation by Arg144 with E199 serving as the base for the H(O3′) (Supplementary Figure S15B). In addition we considered the possibility that each of these two paths occurred simultaneously with the nucleophilic attack. All produced uphill reaction energy profiles (Supplementary Figure S15C).

(iii) The α-phosphate of the dNTP as the base or an intermediate in the path of H(O3′) to the γ-phosphate. Mechanisms involving protonation of the α-phosphate were examined: either a direct protonation or protonation via a water were considered (Supplementary Figure S6). The protonated α-phosphate did not provide a stable intermediate as found for the Dpo4 WMSA mechanism (35); upon a short unrestrained QM/MM–MD the hydrogen went back to the primer O3′ to produce the reactant structure. By contrast the protonated γ-phosphate intermediate in the WMSA mechanism was stable during 35 ps of unrestrained QM/MM–MD simulation.

Disfavored mechanisms remain so irrespective of H(O3′) orientation

We made an extensive investigation of the H(O3′) orientation and its effect on the reaction mechanisms that we studied. We carried out a survey of the free energy profile of the C2′–C3′–O3′–H torsion angle which governs the orientation of the hydrogen of the O3′. This was carried out with AMBER 9 (47) by restraining this torsion angle at 10° intervals from 0° to 360° with 2 ns MD simulations for each window, utilizing a force constant of 30 kcal/mol-rad². The WHAM method was then utilized to obtain the free energy profile (78). Figure 8 shows minima at 45°, 177° and 309° with relative energies of 2.6, 7.8 and 0 kcal/mol, respectively. Our results showed that for the WMSA mechanism it was necessary for the H(O3′) to be near the 60° energy well, which allows the proton to shuttle through the water molecule, as in our earlier work (35,36); however, the disfavored mechanisms remained disfavored irrespective of the H(O3′) orientation (e.g. Supplementary Figure S14).

The pentacovalent phosphorane transition state structure is preserved in the presence of the B[a]P-dG adduct

We wished to evaluate whether the pentacovalent phosphorane transition state in our characterized WMSA mechanism is preserved in the presence of the B[a]P-dG adduct, modeled on the minor groove side of the templating base in Pol κ, as in the NMR solution structure (26). For this purpose, we developed a new Amber-compatible force field for the transition state structure in the unmodified Pol κ (see ‘Computational Methods’ section and Supplementary Tables S1 and S2). Transition states are extremely transient, with a lifetime of the order of 0.01 to 0.1 ps, in the same time frame as a single bond vibration (79). We carried out an MD simulation of the transition state structure in the presence of the B[a]P-dG adduct, and for an unmodified control transition state structure. We wished to determine if the geometric properties of the QM/MM–MD determined transition state for the unmodified Pol κ are preserved when the lesion is present. Supplementary Figure S16 and Supplementary Table S3 shows that the characteristic phosphorane transition state geometry, after 1 ns of classical MD, is preserved as in the QM/MM–MD transition state for both unmodified and modified systems. Figure 5 shows the QM/MM–MD transition state and that of the B[a]P-containing structure simulated by MD. These data reveal that the transition state is maintained in the presence of the lesion as well as it is in the unmodified case, and as found in the QM/MM–MD computations for the unmodified system. Figures 9 and 10, and Supplementary Movie S3 show the adduct structures prior to reaction and in the transition state after the 1 ns MD, revealing that they are very similar. Supplementary Figure S17 shows the linkage site torsion angles χ, α’ and β’; 242 ± 12°, 204 ± 39° and 58 ± 9°, respectively (as defined in Figure 1) that govern the carcinogen orientation on the minor groove side. They are very close to those obtained previously by Jia and coworkers, who performed MD simulations (30) for this adduct positioned in Pol κ on the templating guanine opposite dCTP prior to the nucleotidyl transfer reaction (χ, α′ and β′; 225 ± 13.0°, 187 ± 8.1° and 64.4 ± 9.4°, respectively). These results indicate that the lesion does not interfere with the formation of the transition state as the reaction progresses from its initial state.