Insight into the catalytic mechanism of meta-cleavage product hydrolase BphD: a quantum mechanics/molecular mechanics study

Abstract

The catalytic mechanism of BphD (meta-cleavage product hydrolase, the fourth enzyme of the biphenyl catabolic pathway) toward its natural substrate 2-hydroxy-6-oxo-6-phenylhexa-2,4-dienoic acid (HOPDA) was investigated in atomistic detail by the QM/MM approach. The calculated Boltzmann-weighted average barriers favor a substrate-assisted acylation mechanism, and the most feasible acylation pathway involves a catalytic triad (Ser-His-Asp). The product (2-hydroxypenta-2,4-dienoic acid) of the acylation process is replaced by three water molecules, and one of which is involved in the deacylation process. The established acylation and deacylation mechanism may shed light on investigating the degradation processes of wt BphD toward hundreds of other differently chlorinated HOPDA. The roles of seventeen residues during the catalytic process of wt BphD toward HOPDA were also reported in search of new promising experimental mutation targets for the improvement of BphD catalytic efficiency.

---

1. Introduction

Polychlorinated biphenyls (PCBs), known as persistent aromatic pollutants, were widely manufactured and used during 1920s to 1970s in a variety of industrial applications.¹ Although banned worldwide by the Stockholm Convention in 2004, there are still about 280 million kilograms of PCBs in mobile environmental reservoirs, causing severe environmental concerns.^2–4 It is thus critically desiderated to set up strategies to minimize continued exposure of the carcinogenic bioaccumulative PCBs to humans. Environmental biotransformation, one of the most promising strategies with the lowest power input, has provided some encouraging results and received worldwide acceptance in dealing with PCBs.^5–7

It has been evidenced that higher chlorinated PCB congeners are firstly transformed to lower chlorinated PCB congeners through reductive dechlorination pathways by anaerobic bacteria, while lower chlorinated PCB congeners are degraded by aerobic bacteria.^8,9 Among all the PCB-degrading aerobic bacteria, Burkholderia xenovorans LB400 displays the broadest substrate specificity and thus serves as the model organism for PCB aerobic degradation studies.¹⁰ The gene cluster that is responsible for PCB aerobic degradation was named as bph because biphenyl is the primary substrate for these PCB-degrading bacteria (including Burkholderia xenovorans LB400).^11,12 The four enzymes involved in bph upper pathway are thus named as BphA, BphB, BphC, and BphD, respectively. They can successively degrade biphenyl to cis-2,3-dihydro-2,3-dihydroxybiphenyl,2,3-dihydroxybiphenyl,2-hydroxy-6-oxo-6-phenylhexa-2,4-dienoic acid (HOPDA), and finally to 2-hydroxypenta-2,4-dienoic acid (HPD) plus benzoic acid.⁹ The conversion of HOPDA to HPD and benzoic acid by enzyme BphD is extremely important in the bph pathway and thus gets the most attention.^13–20

BphD (EC 3.7.1.8) belongs to meta-cleavage product (MCP) hydrolases which possesses an alpha/beta-hydrolase fold and utilizes a Ser-His-Asp triad during catalysis.^14,21 The crystallization data show that the active site of BphD has a nonpolar and a polar sub-site which can accommodate the phenyl and dienoate moiety of its natural substrate HOPDA.^16,17 Recent crystallization data, EI/GC/MS studies, and site-directed mutation analysis have evidenced an acylation and deacylation mechanism for BphD.²² Several following experimental investigations have led to the proposal of a substrate-assisted nucleophilic mechanism for BphD acylation.^23–25 These experimental explorations have deepened our understanding on BphD catalytic processes. Although the importance of substrate HOPDA during BphD acylation was spotted, the involvement of His265 is still an open question.^18–22 In addition, there are still many short-lived intermediates cannot be experimentally characterized. Understanding the details of the catalytic processes of BphD can be of great importance in designing new promising catalysts to accelerate the degradation of PCBs.

Combined quantum mechanics/molecular mechanics (QM/MM) method has become an increasingly powerful tool to complement experimental enzyme chemistry.²⁶ In the present work, the detailed catalytic mechanisms of wild type (wt) BphD at atomistic level were performed with the aid of QM/MM method. In order to check the role of the catalytic triad (Ser112-His265-Asp237) during the acylation process of wt BphD, the acylation mechanism of BphD–H265A mutant was also investigated. In addition, electrostatic analysis was performed to investigate the influence of seventeen active-site-constituting residues on the rate-determining step to look for new promising targets for further mutation studies.

2. Methods

The initial models for the acylation and deacylation simulations were built on the basis of the X-ray crystal structure of BphD–HOPDA complex (PDB code 2PUH, resolution 1.82 Å) and acylated BphD (PDB code 3V1N, resolution 1.59 Å) obtained from the Protein Data Bank (http://www.rcsb.org).^18,22 The mutated residues (S112A in 2PUH and H265Q in 3V1N) presented in the crystal structures were manually transformed back into their natural forms. The protonation states of ionizable residues were determined on the basis of the pK_a values obtained via the PROPKA procedure and were manually verified through visual inspection.²⁷ The missing hydrogen atoms of BphD were added through CHARMM22 force field in the HBUILD module of CHARMM package.^28,29 MolProbity software was used to check the flipped Asn/Gln/His residues.³⁰ The entire enzyme was dissolved in a water droplet (TIP3P model³¹) with a radius of 32 angstroms. Then, the enzyme–water system was neutralized by sodium ions via random substitution of solvent water molecules before being relaxed through energy minimizations. The whole system was firstly heated from 0 K to 298.15 K in 50 ps (1 fs per step) and equilibrated thermally for 500 ps (1 fs per step) to reach the equilibration state. After that, a 12 ns stochastic boundary molecular dynamics (SBMD) simulation was performed at 298.15 K by using NVT ensemble for conformational sampling.³² During the SBMD simulations, the whole system moves freely except the substrate (HOPDA for 2PUH, acylated Ser112 for 3V1N). Since both of the crystal structures are determined with high resolution (1.82 Å for 2PUH, 1.59 Å for 3V1N), it is reasonable to restrain the coordinates of the substrates and keep their positions consistent with that in the crystal structures. The leap-frog algorithm and Langevin temperature coupling method implemented in CHARMM program were applied during the simulations.

The QM/MM calculations were performed by using ChemShell³³ platform, which can integrate programs Turbomole³⁴ and DL_POLY.³⁵ The charge shift model³⁶ and electrostatic embedding method³⁷ were used during the QM/MM calculations. The geometries of the intermediates were optimized by using hybrid delocalized internal coordinates optimizer while transition state searches were done by using microiterative TS optimizer at the B3LYP/6-31G(d,p)//CHARMM22 level.³⁸ Microiterative TS optimizer can accelerate TS searches with less computing cost through splitting the system into a reaction core containing with few atoms and the environment (Scheme S1, ESI†). Frequency calculations were performed to validate the one imaginary frequency character of transition state structures, and the suitability of the transition vector was also confirmed. Additional single point energy calculations were carried out at the B3LYP/cc-pVTZ//CHARMM22 level for better describing of the energy profiles.

The QM-region for wt BphD acylation system contains the side chains of residues Ser112, Asp237, His265, and substrate HPODA while the QM-region for BphD–H265A acylation system contains the side chains of residues Ser112, Asp237, Ala265, and substrate HPODA. For deacylation system, the QM-region contains the side chains of residues Asn111, His265, Asp237, acly-Ser112 and a water molecule. The truncated C–C bonds for the three systems were complemented by hydrogen link atoms. There are totally 50, 42, and 55 atoms included in the QM region of wt BphD acylation system, BphD–H265A acylation system, and deacylation system, respectively. All the atoms within 20 angstrom of N1 atom (Scheme 1) from His265 were selected to be the active region (about 3700 movable atoms). Atoms that lie beyond 20 angstrom of N1 were fixed during the QM/MM calculation. Further details of the QM/MM setup can be found in ESI.†

Scheme 1 Substrate-assisted acylation (Ser112-His265-Asp237 involved) and deacylation pathways of wild type (wt) BphD. Atoms (1–7 for C, 1–5 for O, 1–2 for N, and 1–7 for H) are labeled in blue, distances (D1–D20) are labeled in red.

3. Results and discussion

Recent room-temperature single molecule experiments have shown that the rate constant of a single enzyme molecule exhibits large fluctuations with a broad range (1 ms to 100 s), while the experimentally determined overall k_cat is believed to be the integrated results by considering all the rate constant fluctuations.^39–42 By assuming that each snapshot from the dynamics trajectory corresponds to a local rate constant, many excellent studies have shown that the calculated Boltzmann-weighted average barrier correspond well to k_cat.^43,44 Thus, three snapshots extracted at 6, 9 and 12 ns from the trajectory were used as the start points in the present study. The differences of these snapshots were illustrated by some key distances, as shown in Table S1, ESI.† The QM/MM results indicate that different snapshots can be associated with different barriers. To analyze the barrier difference, Boltzmann-weighted averaging method on the rate constant was used:^45–48where, ΔE_ea is the average barrier, R is gas constant, n is the number of snapshots, ΔE_i is the energy barrier of path i, and T is the temperature. If not specially noted, the “average barrier” mentioned in the following paragraphs refers to Boltzmann-weighted average barrier.

3.1 Substrate-assisted acylation and Ser112-His265-Asp237 triad

It should be noted here that the basic idea of substrate-assisted acylation for BphD was proposed by Eltis and co-workers.^16,18,22 In our QM/MM models, three investigated substrate-assisted acylation pathways (based on structures extracted at 6, 9, and 12 ns) are defined as acy-6, acy-9, acy-12, respectively. For each of the three pathways, three major steps were found: enol-to-keto tautomerization (R-1 to IM-3), acylation (IM-3 to IM-4), and HPD formation (IM-4 to IM-5), as indicated in Scheme 1. Energy profiles are provided in Fig. 1 (details in Table S2, ESI†). The step of enol-to-keto tautomerization consists of two elementary reactions. The first elementary reaction is the deprotonation of HOPDA (R-1 to IM-1, average barrier 0.8 kcal mol⁻¹). By ignoring the role of Ser112, a second elementary reaction (IM-1 to IM-2, His265-Asp237 based rearrangement) with a high average barrier (25.9 kcal mol⁻¹) was established. This average barrier agree well with the experimental findings that IM-2 decays slowly in BphD–S112A mutant with a half-life 4.4 h (corresponds to a k_cat of 4.4 × 10⁻⁵ s⁻¹ and a barrier of 23.3 kcal mol⁻¹).¹⁸ However, by considering the role of Ser112, a second elementary reaction (IM-1 to IM-2, Ser112-His265-Asp237 based rearrangement) with a much lower average barrier (8.0 kcal mol⁻¹) can be established. Our calculated results combined with previous experimental findings suggest the non-substitutable character of Ser112 during substrate-assisted acylation. The detailed electron flows of the triad based and the His265-Asp237 based rearrangements are indicated in Scheme 1 in red and blue arrows.

Fig. 1 Energy profiles of three individual pathways along substrate-assisted acylation (Ser112-His265-Asp237 involved) and deacylation pathways calculated at B3LYP/cc-pVTZ//CHARMM22 level. The Boltzmann-weighted average barriers of each elementary step are provided in the braces.

After enol-to-keto tautomerization, the energies of IM-3 are about 6.2, 8.3, and 4.9 kcal mol⁻¹ lower than that of R-1 in pathways acy-6, acy-9, and acy-12. The calculated average barrier for acylation step (IM-3 to IM-4) is 14.1 kcal mol⁻¹, which is higher than triad based rearrangement (8.0 kcal mol⁻¹). IM-4 (acylated BphD) is unstable and it quickly turns into HPD through C–C bond cleavage with an average barrier 4.2 kcal mol⁻¹. As a result, acylation step (IM-3 to IM-4) is rate-determining during the whole substrate-assisted acylation process of BphD (R-1 to IM-5), and the corresponding average barrier (14.1 kcal mol⁻¹) is in good consistency with the pre-steady-state kinetics result (∼15 kcal mol⁻¹).²²

The possibility that only Ser112 but not Ser112-His265-Asp237 is involved in substrate-assisted acylation was also considered. The corresponding reaction pathways and detailed energies were provided in ESI (Scheme S2 and Table S3†). The calculated average barrier of the rate-determining step (21.6 kcal mol⁻¹) is much higher than that in Ser112-His265-Asp237 involved substrate-assisted acylation (14.1 kcal mol⁻¹). This indicates the importance of residue His265. To confirm this conclusion, QM/MM studies on the acylation process of the BphD–H265A mutant were performed, and the calculated average barrier of the rate-determining step is 20.2 kcal mol⁻¹. These data suggest that although substrate HOPDA can be solely attacked by Ser112, the mutation of His265 to Ala265 is accompanied with a significant increase of the average barrier (from 14.1 to 20.2 kcal mol⁻¹). This conclusion was also supported by the experimental findings that the acylation rate of BphD is dramatically reduced after H265A mutation.^18,22

3.2 Deacylation and residue electrostatic analysis

The proposed deacylation mechanism and the corresponding energies are represented in Scheme 1 and Fig. 1. The three water molecules presented in our QM/MM model root in the crystal water molecules (Wat440, Wat455, Wat465) found in X-ray data of BphD–H265Q mutant (PDB code 3V1N, resolution 1.59 Å), as shown in Fig. S1, ESI.† The average barriers for His265 assisted attack and benzoic acid formation are 19.6 and 1.9 kcal mol⁻¹. This makes His265 assisted attack step (R-5 to IM-6) not only the rate-limiting step of the deacylation process (R-5 to P-7), but also the rate-limiting step of the whole acylation and deacylation process. Thus, the electrostatic analysis of seventeen residues on this step was performed, and the results were provided in Fig. 2. All these residues are spatially close to the active site, as indicated in Fig. S2, ESI.† It should be noted here that His265 assisted attack step is rate-limiting for the degradation process of BphD toward substrate HOPDA, but may not be the rate-limiting step toward other differently chlorinated HOPDA, such as 3-Cl-HOPDA.¹⁹ However, the current work can serve as a model study and provide basic concepts on the acylation and deacylation mechanism of BphD. Much more future work aimed at studying the degradation processes toward hundreds of other chlorinated HOPDA is needed. The activation energy difference caused by residue “i” can be described as:where ΔEⁱ⁻⁰ is the energy changes of the barrier, ΔEⁱ is the energy barrier with charges on each atom of residue i set to 0, and ΔE⁰ is the original values of the energy barrier. During all these energy calculations, the geometry structures of the stationary points were kept unchanged. We mention here that a positive value of ΔEⁱ⁻⁰ means that residue “i” facilitates the reaction.

Fig. 2 Electrostatic influences of seventeen residues on the rate-determining step (His265 assisted attack step, R-5 to TS-5-6) of pathway dea-6.

As represented in Fig. 2, the residue influence on the rate-determining step ranges from −3.8 to 2.2 kcal mol⁻¹. In general, the positive or negative value of a residue is due to the different stabilizing ability toward the transition state and the reactant. If a residue stabilizes the transition state more than the reactant, it has a positive value. On the other hand, if a residue stabilizes the reactant more than the transition state, it has a negative value. For instance, large positive value of residue Gly42 is likely due to the stronger hydrogen bond between residue Gly42 and acylated Ser112 in TS-5-6 (1.75 Å) than in R-5 (1.78 Å), as shown in Fig. S3.† Similar situation was found for Trp266: positive value and stronger hydrogen bonds in TS-5-6 than in R-5. This is in agreement with the experimental findings that substitution of Trp266 will reduce the catalytic activity.^15,49 In addition, it is likely that the distance (average distance between each atom of Phe239 and oxygen atom in the water molecule) increase from 4.24 Å (R-5) to 5.00 Å (TS-5-6) is responsible for the large negative value (3.8 kcal mol⁻¹) of Phe239 toward the rate-determining step (Fig. S3†). This indicates that proper mutation of Phe239 may improve the catalytic efficiency of BphD toward HOPDA.

3.3 Catalytic itinerary of wt BphD

The catalytic itinerary of wt BphD was represented by atom–atom distance and NPA (Natural Population Analysis) charge variations along pathway acy-6 and dea-6 (Fig. 3). More detailed information of the fourteen critical points were provided in ESI† (Tables S4 and S5; coordinates of the QM atoms). The catalytic itinerary was listed as follows: deprotonation of HOPDA (R-1, TS-1-2, IM-2), triad based rearrangement (IM-2, TS-2-3-T, IM-3), acylation (IM-3, TS-3-4, IM-4), HPD formation (IM-4, TS-4-5, IM-5), HPD release (IM-5, R-5), His265 assisted water attack (R-5, TS-5-6, IM-6), and benzoic acid formation (IM-6, TS-6-7, P-7).

Fig. 3 Key distance and NPA charge variations along substrate-assisted acylation (Ser112-His265-Asp237 involved) and deacylation pathways (acy-6 and dea-6).

The deprotonation step (R-1 to IM-2) is monitored by distance variation of D7 and D8. The NPA charge of substrate HOPDA varies from −1.02 to −1.86 while the NPA charge of its phenyl group varies from −0.01 to −0.05. This indicates that the newly generated negative charge in substrate HOPDA is hardly distributed to phenyl group. The following triad based rearrangement (IM-2 to IM-3) can be represented by distance increase of D8, D10 and distance decrease of D12, D14. It is worth mentioning that Mulliken charge of atom C5 varies from −0.49 (IM-2) to −0.65 (TS-2-3-T) and −0.58 (IM-3), this suggests that electron sinks are generated in both TS-2-3-T and IM-3 during the acylation of wt BphD, which confirms the previous assumptions.²² Charge of atom O2 drops to −0.87 during the acylation step (IM-3 to IM-4), while it increases back to −0.69 in the following HPD formation step (IM4 to IM5), suggesting the C6=O2 double bond formation. After HPD release, the empty space of the active pocket is fulfilled by three water molecules (IM5 to R-5). Those water molecules are directly or indirectly anchored by surrounding residues (Asn111 and His265) through hydrogen bonds. This agrees with the experimental results that substitution of residue Asn111 reduces the reaction rate.¹⁵ The replacement of HPD by three water molecules increases the charge of C6, and thus makes it more electrophilic. A water molecule activated by His265 attacks C6 atom and thus results in an unstable intermediate IM-6, as represented by D17, D18, and D19 variations (R-5 to IM-6). Distance variation of D11 indicates that the benzoic acid is formed during the last step (IM-6 to P-7).

4. Conclusions

In the present study we investigate the acylation mechanism of wt BphD and BphD–H265A mutant in atomistic detail by using the QM/MM approach. A Ser112-His265-Asp237 triad is involved in the substrate-assisted acylation process of wt BphD. Three water molecules are found in the active pocket while one of them is involved in the deacylation process. Residues Asn111 and His265 are found to regulate the water molecules through hydrogen bonds. Significant roles of residues Phe175, Arg190, and Trp266 reported by experimental method on the rate-determining step were confirmed. Possible mutation targets (such as Phe239) in accelerating the degradation rates of BphD were also highlighted. The established degradation mechanism and residue influence analysis may shed light on investigating the degradation processes of BphD toward hundreds of other differently chlorinated HOPDA.

Acknowledgements

The work was financially supported by NSFC (National Natural Science Foundation of China, project no. 21337001 and 21177077), FRFSU (Fundamental Research Funds of Shandong University, project no. 2015GN007), and GFGCPSF (General Financial Grant from the China Postdoctoral Science Foundation, project no. 2015M570594).

References

1. D. L. Bedard, Polychlorinated biphenyls in aquatic sediments: environmental fate and outlook for biological treatment, in Dehalogenation: Microbial Processes and Environmental Applications, ed.M. M. Haggblom and I. D. Bossert, Kluwer Academic Publishers, Boston, MA, 2003, pp. 446–465.
2. O. M. Faroon, L. S. Keith, C. Smith-Simon, and C. T. de Rose, Polychlorinated Biphenyls: Human Health Aspects. World Health Organization, Geneva, 2003.
3. L. A. Rodenburg, J. Guo, S. Du and G. J. Cavallo, Environ. Sci. Technol., 2010, 44, 2816–2821.
4. C. Steinlin, C. Bogdal, M. Scheringer, P. A. Pavlova, M. Schwikowski, P. Schmid and K. Hungerbühler, Environ. Sci. Technol., 2014, 48, 7849–7857.
5. D. D. Focht, Curr. Opin. Biotechnol., 1995, 6, 341–346.
6. J. Borja, D. M. Taleon, J. Auresenia and S. Gallardo, Process Biochem., 2005, 40, 1999–2013.
7. D. H. Pieper, Appl. Microbiol. Biotechnol., 2005, 67, 170–191.
8. J. F. Brown Jr., D. L. Bedard, M. J. Brennan, J. C. Carnahan, H. Feng and R. E. Wagner, Science, 1987, 236, 709–712.
9. J. A. Field and R. Sierra-Alvarez, Environ. Pollut., 2008, 155, 1–12.
10. P. S. Chain, V. J. Denef, K. T. Konstantinidis, L. M. Vergez, L. Agulló, V. L. Reyes, L. Hauser, M. Córdova, L. Gómez, M. González, M. Land, V. Lao, F. Larimer, J. J. LiPuma, E. Mahenthiralingam, S. A. Malfatti, C. J. Marx, J. J. Parnell, A. Ramette, P. Richardson, M. Seeger, D. Smith, T. Spilker, W. J. Sul, T. V. Tsoi, L. E. Ulrich, I. B. Zhulin and J. M. Tiedje, Proc. Natl. Acad. Sci. U. S. A., 2006, 103, 15280–15287.
11. D. A. Abramowicz, Crit. Rev. Biotechnol., 1990, 10, 241–249.
12. K. Furukawa, J. Gen. Appl. Microbiol., 2000, 46, 283–296.
13. S. Y. Seah, G. Labbé, S. Nerdinger, M. R. Johnson, V. Snieckus and L. D. Eltis, J. Biol. Chem., 2000, 275, 15701–15708.
14. B. Hofer, L. D. Eltis, D. N. Dowling and K. N. Timmis, Gene, 1993, 130, 47–55.
15. C. Li, J. J. Li, M. G. Montgomery, S. P. Wood and T. D. Bugg, Biochemistry, 2006, 45, 12470–12479.
16. G. P. Horsman, J. Ke, S. Dai, S. Y. Seah, J. T. Bolin and L. D. Eltis, Biochemistry, 2006, 45, 11071–11086.
17. S. Y. Seah, J. Ke, G. Denis, G. P. Horsman, P. D. Fortin, C. J. Whiting and L. D. Eltis, J. Bacteriol., 2007, 189, 4038–4045.
18. G. P. Horsman, S. Bhowmik, S. Y. Seah, P. Kumar, J. T. Bolin and L. D. Eltis, J. Biol. Chem., 2007, 282, 19894–19904.
19. S. Bhowmik, G. P. Horsman, J. T. Bolin and L. D. Eltis, J. Biol. Chem., 2007, 282, 36377–36385.
20. H. Zhou, Y. Qu, C. Kong, E. Shen, J. Wang, X. Zhang, Q. Ma and J. Zhou, Appl. Microbiol. Biotechnol., 2013, 97, 10399–10411.
21. D. L. Ollis, E. Cheah, M. Cygler, B. Dijkstra, F. Frolow, S. M. Franken, M. Harel, S. J. Remington, I. Silman, J. Schrag, J. L. Sussman, K. H. G. Verschueren and A. Goldman, Protein Eng., 1992, 5, 197–211.
22. A. C. Ruzzini, S. Ghosh, G. P. Horsman, L. J. Foster, J. T. Bolin and L. D. Eltis, J. Am. Chem. Soc., 2012, 134, 4615–4624.
23. A. C. Ruzzini, G. P. Horsman and L. D. Eltis, Biochemistry, 2012, 51, 5831–5840.
24. A. C. Ruzzini, S. Bhowmik, S. Ghosh, K. C. Yam, J. T. Bolin and L. D. Eltis, Biochemistry, 2013, 52, 7428–7438.
25. A. C. Ruzzini, S. Bhowmik, K. C. Yam, S. Ghosh, J. T. Bolin and L. D. Eltis, Biochemistry, 2013, 52, 5685–5695.
26. W. Thiel and G. Hummer, Nature, 2013, 504, 96–97.
27. H. Li, A. D. Robertson and J. H. Jensen, Proteins, 2005, 61, 704–721.
28. B. R. Brooks, R. E. Bruccoleri, B. D. Olafson, D. J. States, S. Swaminathan and M. Karplus, J. Comput. Chem., 1983, 4, 187–217.
29. B. R. Brooks, C. L. Brooks III, A. D. Mackerell, L. Nilsson, R. J. Petrella, B. Roux, Y. Won, G. Archontis, C. Bartels, S. Boresch, A. Caflisch, L. Caves, Q. Cui, A. R. Dinner, M. Feig, S. Fischer, J. Gao, M. Hodoscek, W. Im, K. Kuczera, T. Lazaridis, J. Ma, V. Ovchinnikov, E. Paci, R. W. Pastor, C. B. Post, J. Z. Pu, M. Schaefer, B. Tidor, R. M. Venable, H. L. Woodcock, X. Wu, W. Yang, D. M. York and M. Karplus, J. Comput. Chem., 2009, 30, 1545–1615.
30. V. B. Chen, W. B. Arendall III, J. J. Headd, D. A. Keedy, R. M. Immormino, G. J. Kapral, L. W. Murray, J. S. Richardson and D. C. Richardson, Acta Crystallogr., Sect. D: Biol. Crystallogr., 2010, 66, 12–21.
31. W. L. Jorgensen, J. Chandrasekhar, J. D. Madura, R. W. Impey and M. L. Klein, J. Chem. Phys., 1983, 79, 926–935.
32. C. L. Brooks and M. Karplus, J. Chem. Phys., 1983, 79, 6312–6325.
33. P. Sherwood, A. H. D. Vries, M. F. Guest, G. Schreckenbach, C. R. A. Catlow, S. A. French, A. A. Sokol, S. T. Bromley, W. Thiel, A. J. Turner, S. Billeter, F. Terstegen, S. Thiel, J. Kendrick, S. C. Rogers, J. Casci, M. Watson, F. King, E. Karlsen, M. Sjovoll, A. Fahmi, A. Schafer and C. Lennartz, J. Mol. Struct.: THEOCHEM, 2003, 632, 1–28.
34. R. Ahlrichs, M. Bär, M. Häser, H. Horn and C. Kölmel, Chem. Phys. Lett., 1989, 162, 165–169.
35. W. Smith and T. R. Forester, J. Mol. Graphics, 1996, 14, 136–141.
36. A. H. de Vries, P. Sherwood, S. J. Collins, A. M. Rigby, M. Rigutto and G. J. Kramer, J. Phys. Chem. B, 1999, 103, 6133–6141.
37. D. Bakowies and W. J. Thiel, J. Phys. Chem., 1996, 100, 10580–10594.
38. S. R. Billeter, A. J. Turner and W. Thiel, Phys. Chem. Chem. Phys., 2000, 2, 2177–2186.
39. H. Yang, G. Luo, P. Karnchanaphanurach, T. M. Louie, I. Rech, S. Cova, L. Xun and X. S. Xie, Science, 2003, 302, 262–266.
40. X. S. Xie, Science, 2013, 342, 1457–1459.
41. H. P. Lu, L. Xun and X. S. Xie, Science, 1998, 282, 1877–1882.
42. W. Min, B. P. English, G. Luo, B. J. Cherayil, S. C. Kou and X. S. Xie, Acc. Chem. Res., 2005, 38, 923–931.
43. P. Saura, R. Suardíaz, L. Masgrau, J. M. Lluch and À. González-Lafont, ACS Catal., 2014, 4, 4351–4363.
44. R. Lonsdale, K. T. Houghton, J. Żurek, C. M. Bathelt, N. Foloppe, M. J. de Groot, J. N. Harvey and A. J. Mulholland, J. Am. Chem. Soc., 2013, 135, 8001–8015.
45. R. Lonsdale, J. N. Harvey and A. J. Mulholland, Chem. Soc. Rev., 2012, 41, 3025–3038.
46. R. Lonsdale, J. N. Harvey and A. J. Mulholland, J. Phys. Chem. B, 2010, 114, 1156–1162.
47. Y. Li, X. Shi, Q. Zhang, J. Hu, J. Chen and W. Wang, Environ. Sci. Technol., 2014, 48, 5008–5016.
48. Y. Li, R. Zhang, L. Du, Q. Zhang and W. Wang, RSC Adv., 2015, 5, 13871–13877.
49. J. J. Li, C. C. Li, A. Blindauer and T. D. H. Bugg, Biochemistry, 2006, 45, 12461–12469.

---

Footnote

† Electronic supplementary information (ESI) available: Key distances of the reactants during the substrate-assisted acylation (R-1) and deacylation pathway (R-5) at 6, 9, 12 ns (Table S1); relative energies of substrate-assisted acylation and deacylation pathway (Tables S2 and S3); key distances and charges along the substrate-assisted acylation and deacylation pathway (Tables S4 and S5); core atoms used during microiterative TS optimizations (Scheme S1); itinerary of the substrate-assisted acylation (Ser112 involved) for both wt BphD and BphD–H265A mutant (Scheme S2); additional details on the methods; coordinates and figures of structures R-1, TS-1-2, IM-2, TS-2-3-T, IM-3, TS-3-4, IM-4, TS-4-5, IM-5, R-5, TS-5-6, IM-6, TS-6-7, and P-7. See DOI: 10.1039/c5ra09939k

---