Computational Assignment of the Histidine Protonation State in (6-4) Photolyase Enzyme and Its E ﬀ ect on the Protonation Step

: The initial step to resolving the controversy regarding the repair mechanism of the pyrimidine (6-4) pyrimidone photoproduct DNA lesion is to determine the protonation state of the two conserved active site histidine residues (His365 and His369 in Drosophila melanogaster ) in the (6-4) photolyase enzyme ((6-4) PHR). Both His residues were experimentally determined to be crucial for catalysis. Most previous theoretical studies assumed the presence of protonated His365 in the active site, which would transfer its proton to N3 ′ of the lesion as the ﬁ rst step in repair; however, other empirical/semiempirical p K a calculations suggested the presence of two neutral His residues in the active site. Here, we conduct molecular dynamics simulations (MD) of the (6-4) PHR/DNA complex on the basis of the X-ray crystal structure to investigate all combinations of the three possible protonation states of each His. The MD results show that HIP365 (both ND1 and NE2 protonated) and HID369 (only ND1 protonated) are the most probable protonation states in the active site. Furthermore, we employ quantum mechanics-cluster and quantum mechanical/molecular mechanical (QM/MM) approaches to investigate the protonation of N3 ′ , starting with three plausible complexes resulting from MD (HIP365/HID369, HID365/HIE369, and HID365/HIP369). Surprisingly, protonation of the N3 ′ atom is found to be feasible starting from all three protonation states: protonated HIP365 transfers the proton via a barrierless step, and in the other cases of neutral HID365, the proton is transferred from O4 ′ of the lesion via His365 through an energy barrier of ∼ 80 kJ mol − 1 in both complexes. These results might explain the conservation of repair activity over a wide range of pH values.


■ INTRODUCTION
In DNA, adjacent pyrimidine bases can undergo harmful covalent modifications due to ultraviolet radiation, forming either cyclobutane pyrimidine dimers (CPD) or pyrimidine  pyrimidone photoproduct ((6-4) PP). 1,2 CPD is the most commonly formed DNA photolesion. 3 Fortunately, to conserve genetic information, these lesions are readily repaired either via nucleotide excision repair, by base excision repair, or by direct enzymatic cleavage of the modified covalent bonds in the pyrimidines. 1,4 Photolyase enzymes (PHRs) undo these covalent modifications using photorepair mechanisms. 4−6 On the basis of substrate specificity, PHRs are classified as CPD PHR or (6-4) PHR: 1,4 the former only repairs CPD, and the latter exclusively repairs  PP. The general repair mechanisms in both classes are quite similar. 2 First, PHRs absorb near the UV/blue light region via a photoantenna pigment such as 8-hydroxy-5-deazaflavin (8-HDF) or 5,10methenyltetrahydrofolate (MTHF), and then the excitation energy is transferred to the reduced flavin adenine dinucleotide (FADH) cofactor. Subsequently, FADH − transfers an electron to the damaged nucleotide, leading finally to its repair.
While the repair mechanism in CPD PHR is well understood, 7−9 the detailed repair pathway in (6-4) PHR remains elusive. 1 Many experimental and theoretical studies have proposed hypotheses explaining (6-4) PP repair, but the mechanism remains controversial. 4,5,10−15 This uncertainty mainly originates from the complexity of the (6-4) PP structure, which requires an intramolecular oxygen and hydrogen transfer from C5 of the 5′-base pyrimidine to C4′ and N3′ of the 3′-base pyrimidine, respectively.
Originally, the repair mechanism was proposed to occur via the formation of an oxetane intermediate upon O4′−C4′ bond formation, prior to electron transfer. 16 However, several theoretical studies have shown that the barrier for this step is energetically unfeasible, at ∼230 kJ mol −1 . 12,13 Alternatively, another study suggested the formation of an oxetane intermediate after electron transfer and protonation of N3′. 17 However, to date no oxetane intermediate has been experimentally identified. Faraji et al. 13 used quantum mechanical/molecular mechanical (QM/MM) calculations and proposed the direct transfer of the OH through an oxetane transition state via a 73.6−82.8 kJ mol −1 barrier. In contrast, Sadeghian et al. 11 suggested the formation of an oxetane intermediate through transfer of the electron back to FADH, followed by cleavage of the formed four-membered ring. Subsequently, a second photon excitation leads to a second electron transfer to the substrate, allowing pyrimidine restoration. It is important to mention that all the previous discussed mechanisms share a similar first step: namely, proton transfer from His365 to N3′ of the substrate. Recently, other oxetane-bypass mechanisms have been proposed, such as repair via a conical intersection of the electronic exited state as well as repair via the formation of a cationic radical (6-4) PP. 15,18 The active site in (6-4) PHR from Drosophila melanogaster (D. melanogaster) contains two indispensable histidine residues, His365 and His369. 10 A previous experimental study revealed that mutation of either His deactivates the enzyme, thus demonstrating their importance for repair. 19 As a result, several studies have proposed that His365 protonates the pyrimidine N3′ atom as the first step in the mechanism (see Scheme 1). 1,11,13 Notably, experimental studies also showed an ∼50% reduction of the repair catalytic rate in D 2 O, suggesting a protonation step in the mechanism. 4,19 Furthermore, a femtosecond fluorescence spectroscopy study by Li et al. 4 indicated that a proton transfer between an active site His and the substrate is a crucial initial step for the repair because it prevents fast back-electron transfer (50 ps), thus leading to catalysis. It is noteworthy that maximum repair activity was previously shown to occur at pH 8.5 and 7.8 in Xenopus laevis (X. laevis) and D. melanogaster (6-4) PHR, respectively. 10,19 More recently, the catalytic rate was shown to be unchanged over the range pH 7−9. 4 In general, (6-4) PHR repair activity was observed over the wide range of pH 6−10. 10 However, at pH 6, the catalytic rate was found to be significantly reduced. 10 Despite these proposed mechanisms, the actual protonation state of the His residues in the (6-4) PHR active site remains unidentified. Condic-Jurkic et al. 20 used several computational programs, including PROPKA, H++, and APBS, to estimate the His pK a values and suggested that both His residues are neutral, which disagrees with the aforementioned experimental observations. Furthermore, the authors also suggested that His369 is more likely to be protonated than His365. In contrast, Sadeghian et al. 11 suggested the presence of protonated His365 in the active site on the basis of a short (1 ns) molecular dynamics (MD) simulation, since the presence of two neutral His residues leads to a large rootmean-square deviation (RMSD) from the X-ray structure.
Clearly, determining the protonation states of the two His residues in the (6-4) PHR active site is the first step to resolving the dilemma regarding the repair mechanism. MD simulation was demonstrated to be a powerful tool to assess His protonation states in proteins, allowing measurement of the RMSD of the studied His and its neighboring residues with respect to the starting X-ray structures. 21 In this work, we conduct MD simulations for 100 ns for nine possible combinations of the two His residues in three protonation states in (6-4) PHR. Furthermore, quantum mechanical/ molecular mechanical (QM/MM) calculations and the QMcluster approach are employed to determine the previously proposed first step in the N3′ protonation reaction via the active site His from the best candidates resulting from the MD.

Molecular Dynamics
Simulations. The (6-4) PHR complex (PHR complexed with the T(6-4)T nucleotide and FADH − ) X-ray crystallographic structure of D. melanogaster (PDB ID: 3CVU) 10 was employed as the starting structure in the following simulations. This X-ray structure determined at 2.0 Å resolution was prepared at a pH value of 7.8, at which the (6-4) PHR of D. melanogaster has the maximum activity. Three protonation states were considered for both the His365 and the His369 residues: HID, HIE, and HIP, representing the protonation of only ND1, only NE2, and both ND1 and NE2, respectively. From the crystal structure, a total of nine complexes that represent all possible combinations of the protonation states were prepared and investigated using MD simulations (see Figure 1). MD simulations were performed Scheme 1. Illustration of the Proposed 1,11,13 First-Step Reaction in (6-4) PHR Repair Mechanism Figure 1. Possible combinations of His365 and His369 protonation states. HID, HIE, and HIP represent histidine residues with only ND1, only NE2, and both ND1 and NE2 protonated, respectively. The notations and colors used to represent these combinations are used in the remainder of the paper.

ACS Catalysis
Research Article using the AMBER14 software package. 22 The ff14SB force field was employed to parametrize the protein and DNA. Similarly, the generalized AMBER force field (GAFF) parameters were used for the damaged nucleotides and FADH − , for which AM1-BCC 23 was used to generate the corresponding atomic partial charges. The AMBER program xLEAP was used to prepare the MD starting structures by adding missing hydrogen atoms and sodium ions (Na + ) to neutralize the complexes. Then, all the complexes were solvated using the TIP3P water model in a box with a minimum buffer distance of 10 Å from the box edge. Ion parameters were obtained from the ff94/ff99 force field. It is important to mention that His protons were added using xLEAP and then examined for possible conflict in proton directions. Since no significant conflict was found, the MD simulations were started from the generated hydrogen-bonding network. To confirm the robustness of the result, we also performed a distinct initial MD procedure with positional restraints for heavy atoms while hydrogen atoms were relaxed freely up to the MD equilibration step. No significant difference was found between the results of the two equilibration procedures.
The following process was simulated using the PMEMD module in AMBER14. First, the solvated structures were energy minimized with positional restraints on the solute atoms for a total of 20000 steps: 5000 steps with steepest descent (SD), followed by 15000 steps using the conjugate gradient (CG) method. Each entire system was then freely minimized for a total of 30000 (10000 SD + 20,000 CG) steps. Next, the minimized solvated complexes were gradually heated from 0 to 300 K with positional restraints imposed on the solute with a force constant of 41.8 kJ mol −1 Å −2 , and then all the complexes were equilibrated at 300 K and 1 atm for 100 ps using the NPT ensemble. Finally, production runs were performed on the equilibrated complexes for 100 ns with the NPT ensemble (Langevin thermostat), using a time step of 2 fs. The SHAKE algorithm was employed for all of the bonds involving hydrogen atoms. Electrostatic interactions were treated using the particle mesh Ewald (PME) method within boundary conditions. 24 The VMD program was used to visualize the generated trajectories. 25 Similarly, the AMBER program CPPTRAJ was used to analyze and cluster all the trajectories using a default algorithm. 26 Five main clusters were created from each simulation on the basis of the RMSD of the active site residues. One representative structure (cluster centroid) of the most populated cluster was used as a starting point for the following QM-cluster and QM/MM calculations.
QM-Cluster Calculations. All calculations were performed using the Gaussian 09 suite of programs. 27 The unrestricted hybrid-meta-GGA density functional method UM06-2X was used for geometry optimization at the 6-31G(d) basis set. 28 The effect of polarization function addition to hydrogen atoms on the optimized geometries was also tested using the 6-31G(d,p) basis set, for which no significant difference was found to occur in the optimized geometries or relative energies in the DE complex (see Table S1 in the Supporting Information). Relative energies were obtained at the UM06-2X/6-311+G(2df,p) level of theory. Similarly, single-point calculations were performed using the UB3LYP-D3 functional 29−31 (see Table S2 in the Supporting Information). Furthermore, the effect of the protein environment was included via the integral equation formalism polarizable continuum model as implemented in Gaussian 09. 32 Specifi-cally, solvation corrections were obtained at the IEF-PCM-UM06-2X/6-31G(d) level of theory with a dielectric constant of 4, which mimics the protein environment surrounding the QM region. Gaussian09d program default parameter values were used for IEFPCM calculations.
On the basis of the MD results, three QM-cluster complexes, PD (HIP365/HID369), DE (HID365/HIE369), and DP (HID365/HIP369), generated from the average structure of the corresponding protonation states from the 100 ns MD simulations were investigated. The residues interacting with the active site substrate were included in the QM-cluster model, and the rest of the complex was deleted. In order to maintain the integrity of the truncated active site models, a minimal number of atoms (10 atoms) was kept fixed during the optimization process. The PD QM-cluster model (132 atoms) includes truncated models of the damaged T(6-4)T nucleotides, FADH − , Gln299, Tyr306, His365, His369, Tyr423, Asn406, and a water molecule. Similarly, the DE model (109 atoms) includes truncated models of the damaged T(6-4)T nucleotides, FADH − , Gln299, His365, His369, and Tyr423. Finally, the DP model (126 atoms) includes the damaged T(6-4)T nucleotides, FADH − , Gln299, His365, His369, Tyr423, Asn406, and a water molecule (see Figure S1 in the Supporting Information). In general, Tyr, His, and Gln/Asn were modeled as methanol, 4-ethylimidazole, and acetamide, respectively.
QM/MM Models and Calculations. Similarly to the QMcluster approach, all calculations were performed using Gaussian 09. 27 The ONIOM QM/MM subtractive scheme was used as implemented in Gaussian 09. 33 Similarly, UM062X 28 was used to optimize the high layer in the QM/ MM calculations, whereas the AMBER96 force field was used for the low-layer calculations. 34 The optimized geometries were obtained using the 6-31G(d) basis set in the mechanical embedding formalism. Subsequently, relative energies were calculated using a single-point calculation on all stationary points at the ONIOM UM06-2X/6-311+(2df,p) level of theory, using the electrostatic embedding formalism. All of the transition states were freely optimized and subsequently confirmed using frequency and intrinsic reaction coordinate (IRC) calculations. All QM/MM atoms were allowed to move freely in the optimization step.
The QM/MM starting structures for the PD, DE, and DP complexes were also obtained from the average structure of the corresponding MD simulations. The average structures from the MD simulations were first minimized at the MM level of theory. Then, all water molecules 5 Å or more from the solute were deleted to reduce the system's size. Finally, the QM/MM models for PD, DE, and DP were generated, in which the QM layers consist of the R group of the same residues included in the QM-cluster approach (see Figure 2 and Figure S2 in the Supporting Information).

■ RESULTS AND DISCUSSION
Protonation State of His365 and His369. In order to determine the protonation state of His365 and His369, the MD trajectories of the nine possible combinations (see Computational Methods) were analyzed using the following criteria. First, the RMSDs of the two active site His residues were calculated throughout the simulation with respect to the initial PDB structure. Second, the RMSDs of the active site (including Gln299, Trp302, His365, His369, Trp409, Tyr423, FADH − ,and (6-4) PP lesion) were also calculated. Third, active site interactions characterized by heavy-atom distances were

ACS Catalysis
Research Article calculated and compared to the X-ray distances. Finally, the changes in the dihedral angle defined by CA−CB−CG−ND1 of either His were also monitored. Note that the first 5 ns of each MD simulation was excluded from the analysis. Table 1 and Figures S3 and S4 in the Supporting Information show the calculated average RMSD values of His365 and His369 with respect to the X-ray structure. The lowest His RMSD was found to occur in the PD (HIP365 and HID369) protonation state, with a very low deviation of 0.27 Å. The DE (HID365 and HIP369) neutral protonation state was found to have the second lowest deviation, with an average value of 0.31 Å. In summary, the following RMSD order was obtained: PD (0. Similarly, the average RMSD of the active site residues shows that the PD protonation state has the lowest deviation (0.56 Å) and that the DE protonation state has the second lowest deviation (0.64 Å). The order of all protonation state RMSD values is as follows: PD (0.56) < DE (0.64) < DP (0.67) < ED (0.77) < EE (0.92) < DE (0.95) < PP (1.21) < PE (1.13) < EP (1.21). In conclusion, the RMSD analysis of the active site residues also suggests that the combination of protonated HIP365 and neutral HID369 in the active site is the most probable active site protonation state.
The only protonation states that show less than 0.5 Å deviation from all the X-ray distances were the HIP365/ HID369 (PD) and HID365/HIE369 (DE) protonation states. In addition, the PD trajectory also shows the presence of a water molecule in the active site. Notably, 90% of the simulation shows the presence of a water molecule in the active site (the water penetrates into the active site after 10 ns of the simulation). This water (W) was absent in the DE protonation state. It is important to note that we conducted a second MD simulation for the DE complex and confirmed that no water was stable in the active site. In our simulation, the crystal water was deleted prior to the MD and the detected water was found to enter the active site during the simulation. The other protonation state (HID365/HIE369 (DE)) consistently showed greater heavy-atom distances: only the HID365/HIP369 (DP) protonation state showed a D3 distance with less than 1.0 Å deviation from the X-ray structure, and the other deviations were smaller than 0.5 Å.
The distance fluctuation between the heavy-atom pairs was further analyzed, especially for the three most promising protonation state complexes (PD, DE, and DP). Figure 4 shows that the investigated interactions in the PD complex are maintained throughout the simulation, except for D3, which shows higher fluctuation with a standard deviation of 0.34 Å. Similar analyses were performed for the DE ( Figure S5 in the Supporting Information) and DP complexes ( Figure S6 in the Supporting Information). The DP complex shows higher fluctuation with respect to D2 ( Figure S6).
When the analyses of the RMSD and heavy-atom distances are combined, the PD (HIP365/HID369) protonation state was concluded to be the most probable in (6-4) PHR. The DE (HID365/HIE369) complex was ranked second in our analysis; however, as described above, it lacks the presence of a water molecule in the active site. The DP protonation state was ranked third. These results are in agreement with the previously proposed mechanisms in which His365 is expected to be protonated. 1,11,13 It is important to mention that the previous study on His pK a values using empirical programs 20 failed to

ACS Catalysis
Research Article predict the superiority of the PD protonation state in the active site.

First Reaction
Step in the (6-4) PHR Mechanism. Subsequently, we investigated the effect of the active site His protonation state on the first reaction step in the (6-4) PHR mechanism proposed previously. 4 As mentioned in the Introduction, the first step is suggested to be a proton transfer. Specifically, the protonation of N3′ of the (6-4) PP substrate was proposed in several studies. 1,11,13 Although the MD results suggested that the PD protonation state was the most probable protonation state, we investigated the mechanism using the top three candidates (PD, DE, and DP). Figure 5 shows that the first step of the reaction occurs through a proton transfer from the HIP365 residue to N3′ via a low barrier of 1.3 kJ mol −1 , as determined using the QM-cluster approach. The optimized reactive complex (RC) shows a strong hydrogen bond between HIP365 and N3′, with a length of 1.67 Å. As shown in Figure 5, the proton in the optimized transition state (TS1) is located a median distance between NE2 HIP365 and N3′: 1.23 and 1.35 Å, respectively. The structure of the following intermediate (I1) was found to be −34.5 kJ mol −1 more stable than the RC. The same step was also conducted using the QM/MM approach, and the proton was found to transfer freely during the optimization of the RC, forming I1.
Notably, the protonation of N3′ was also found to occur when the two His residues in the active site were neutral (DE protonation state), but protonation occurred at a higher energy cost in comparison to the PD complex. As can be seen in Figure 6, proton transfer occurs from O4′ to N3′ with the assistance of HID365. The energy barriers for this step were found to be 59.3 and 84.0 kJ mol −1 using the QM-cluster and QM/MM approaches, respectively. At the QM/MM level of theory, the optimized RC shows that HID365 is strongly

ACS Catalysis
Research Article hydrogen bonded to H O4′ , with a distance of 1.75 Å. In TS1, the proton is partially transferred to His365, with a H··· NE2 HID365 distance of 1.08 Å. In addition, the proton is located 1.59 and 2.28 Å from O4′ and N3′, respectively. Subsequently, the proton was found to be fully transferred to N3′, thus forming I1. The formed intermediate contains the oxyanion O4′, which is stabilized via hydrogen bonding with HIE369. I1 was found to be more stable than RC by −14.3 kJ mol −1 at the QM/MM level of theory. This intermediate was found to be less stable using the QM-cluster approach (see Figure 6).
Finally, we also investigated the protonation of N3′ in the DP complex. This step was quite similar to TS1 in the DE complex, in which the O4′ proton is transferred to N3′ with the assistance of HID365. However, the optimized TS1 showed that two other proton transfers occurred in the same TS1, from the protonated HIP369 to O4′ via a water molecule (Figure 7). The barriers for this step were found to be energy costs of 78.2 and 85.3 kJ mol −1 using the QM/MM and QM-cluster approaches, respectively. It should be noted that the mechanism was found to involve the MD-recovered crystallized water molecule (W) as it bridges the HIP369 interaction with O4′. The formed I1 was more stable than the RC by −142. In summary, our results show that the protonation of N3′ can occur in the presence of three protonation states of His365 and His369 in the (6-4) PHR complex. These results might

ACS Catalysis
Research Article explain the previous experimental results showing the conservation of PHR repair activity over a wide range of pH values, from pH 6−10, in Xenopus PHR. 19 These results might also explain the unchanged repair activity between pH 7 and 9 in the experimental study by Li et al. 4 Subsequent steps in the mechanism were also investigated in the three complexes (PD, DE, and DP); however, we could not freely optimize (not a scan) a TS that represents direct OH transfer in the presence of a 3′-base radical (data not shown). Instead, we were able to optimize a TS from the oxyanion O4′, but the energy barriers were found to be too high (see Figure  S7 in the Supporting Information). These results suggest that O4′ transfer occurs either via the previously proposed twophoton mechanism 11 or via an unknown alternative mechanism. Further investigations to elucidate these points are in progress.

■ CONCLUSION
In this study, the protonation states of the (6-4) PHR active site histidine residues (His365 and His369) were elucidated using MD simulations. All possible combinations of the two His protonation states were investigated for 100 ns with nine MD simulations of the (6-4) PHR/damaged DNA complex. Each protonation state was assessed by comparing the results of the MD trajectory to the starting X-ray structure using several criteria, including the His RMSD values, active site residue RMSD values, and active site heavy-atom distances, as well as the recovery of an X-ray water of crystallization during the simulation.
The PD complex comprising the HIP365/HID369 protonation state was found to be the closest to the X-ray structure. The two neutral His DE complexes (HID365/HIE369) were found to rank second; however, they lack the recovery of the crystal water molecule. It should be emphasized that the DE complex was previously suggested to be the most probable protonation state on the basis of His pK a calculations using empirical/semiempirical programs. 20 Our results suggest that the PD protonation state is the most probable protonation state in the active site. Our results also show that empirical/ semiempirical pK a calculations might fail in protein/nonprotein complexes.
Second, the first reaction step in the repair mechanism of (6-4) PHR was investigated using the QM-cluster and QM/MM approaches. Several previous studies suggested that the protonation of N3′ of the 3′-base of (6-4) PP is the first step in the mechanism, and therefore the proton transfer to N3′ was examined in the top three ranked protonation states (PD, DE, and DP).
In the DP complex, the proton was shown to be transferred from His365 to N3′ via a low barrier of 1.3 kJ mol −1 using the QM/MM and QM-cluster approaches, in agreement with previous theoretical results. Notably, the pronation of N3′ was also found to occur starting with two neutral His residues (DE complex), in which a proton is transferred from O4′ to N3′ via HID365 and feasible barriers of 59.3 and 84.0 kJ mol −1 as determined by the QM-cluster and QM/MM approaches, respectively.
A similar proton transfer mechanism was also observed starting with the HID365/HIP369 protonation state. However, this proton transfer was accompanied by another proton transfer from HIP369 to O4′ via a water molecule. The barrier for the three proton transfer TSs was found to be quite similar to that of the DE complex QM/MM energy, with values of 78.2 and 85.3 kJ mol −1 as determined by the QM/MM and QMcluster approaches, respectively.
As a result, the protonation of N3′ could occur starting with any of these protonation states (PD, DE, and DP). These results are consistent with a previous experimental study that showed the conservation of repair activity over a large range of pH values (pH 6−10). 19 In summary, our results determined the protonation state of the (6-4) PHR His, which is the first step to determining the repair mechanism. In addition, our results show the effectiveness of using MD simulations for determining His protonation states in proteins.