Multiscale Simulations on the Catalytic Plasticity of CYP76AH1

The catalytic promiscuity and fidelity of cytochrome P450 enzymes are widespread in the skeletal modification of terpenoid natural products and have attracted much attention. CYP76AH1 is involved in key modification reactions in the biosynthetic pathway of tanshinone, a well-known medicinal norditerpenoid. In this work, classical molecular dynamic simulations, metadynamics, and DFT calculations were performed to investigate the protein conformational dynamics, ligand binding poses, and catalytic reaction mechanism in wide-type and mutant CYP76AH1. Our results not only reveal a plausible enzymatic mechanism for mutant CYP76AH1 leading to various products but also provide valuable guidance for rational protein engineering of the CYP76 family.

Introduction

Terpenoids, such as the well-known artemisinin, taxol, tanshinones, and ginsenosides, have received increasing attention for their multiple pharmaceutical values (Yang W. et al., 2020; Zheng et al., 2019). It is well-known that the diverse biological activities of terpenoids are strongly related to the oxidative modification of the cyclic molecular skeletons of terpenoids (Sawai and Saito, 2011). Generally, terpene synthases (TPS) make up the major cyclic carbon skeleton, while its subsequent oxidation reactions of the molecular scaffold are involved by various cytochrome P450s (Bathe and Tissier, 2019). Tanshinones, a kind of abietane norditerpenoid present in the plant Salvia miltiorrhiza, have been widely used for the treatment of cerebrovascular and cardiovascular-related diseases (Zhou et al., 2005), and its potential antibacterial, antioxidant, anti-inflammatory, and anticancer biological activities have also been reported (Dong et al., 2011; Robertson et al., 2014). As shown in Figure 1, the cyclization of linear GGPP (geranylgeranyl pyrophosphate) into cyclic miltiradiene is the common step in the biosynthesis of tanshinones (Gao et al., 2009). And P450s are responsible for the subsequent oxidative modification of miltiradiene. So far, several P450s in the CYP76 family were identified to be involved in various oxidative modifications of tanshinones (Guo et al., 2016; Guo et al., 2013). Previous experiments revealed that CYP76AH1 (Guo et al., 2013) could transform miltiradiene to ferruginol via hydroxylation at the C12 position (see Figure 1). It was further demonstrated that the substrate of CYP76AH1 is readily oxidized to the aromatic intermediate abietatriene, which can be spontaneously converted from miltiradiene without enzymatic catalysis (Zi and Peters, 2013). Later, CYP76AH3 was found to be involved in further C7- or (and) C11-oxidation of ferruginol to produce sugiol, 11-hydroxy ferruginol, and 11-hydroxy sugiol (Guo et al., 2016). Recently, it was found that CYP76AH1^D301E,V479F efficiently produces a variety of products, including 11-hydroxy sugiol (major), 11-hydroxy ferruginol, sugiol, and ferruginol (minor), with the two mutant residues in the same position as the corresponding positions of CYP76AH3 (Mao et al., 2020), as shown in Supplementary Figure S1, that is, CYP76AH1 would exhibit promiscuous catalytic function through the D301E and V479F double mutations, whereas its wild type produces ferruginol with high fidelity (as summarized in Figure 1).

FIGURE 1

FIGURE 1. Product fidelity and promiscuity in wild-type and mutant CYP76AH1.

Regarding the catalytic fidelity and promiscuity of P450s, it has attracted extensive computational and experimental studies. In experimental studies, site-directed mutations have been widely used to identify the plastic residues responsible for fidelity and promiscuity, aiming at more efficient production of desired products (Scheler et al., 2016; Sun et al., 2020; Zhang et al., 2020). In the computational field, various molecular dynamic (MD) simulations, DFT, and QM/MM calculations have been employed to probe key features of P450s, including regioselectivity, stereoselectivity, pocket accessibility, substrate specificity, reactivity, and degree of oxidation (Lonsdale et al., 2013; Guo et al., 2016; Bonomo et al., 2017; Leth et al., 2019; Yang S. et al., 2020). In this study on CYP76AH1, classical and metadynamic MD simulations were performed to examine the differences in pocket dynamics between CYP76AH1 and CYP76AH1^D301E,V479F, in order to detect the key factors regulating its catalytic fidelity and promiscuity. A rational reaction pathway for the oxidative modification was then proposed and finally demonstrated by DFT calculations. Our results not only reveal the chemical logic of the enzymatic mechanism corresponding to the active pocket plasticity of CYP76AH1 but also provide valuable guidance for rational protein engineering of CYP76AH1 to enhance its catalytic efficiency or selectivity.

Computational Details

Preparation of Simulation Systems

The crystal structure of CYP76AH1 (PDB ID: 5YLW) (Gu et al., 2019) was utilized to build eight modeling systems, as summarized in Table1. As well-known, the resting state and Compound I (Cpd I) state (see Figure 1) are two critical states in the P450 catalytic cycle (Meunier et al., 2004; Dubey et al., 2016). The systems with resting state are aiming to analyze the binding pocket and the plausible binding poses of ferruginol. The systems of Compound I state are the reactive states prepared to detect reactivity and possible reaction pathways. Structures of the mutant enzyme and Cpd I were prepared by Maestro of Schrödinger package (version 11.1) (Schrödinger, LLC: New York, NY, 2017). Since the ligand is missing in crystal structures, molecular docking of abietatriene (abi for short) and ferruginol (fer for short) was performed with Glide module. Two molecules were generated and prepared by LigPrep module. The rules for selecting optimal docking pose were defined as below: 1) the distance of Fe-C12 should be less than 5 Å (systems wt/mut-abi-Rest); 2) distance of O-C11 should be less than 4 Å (systems wt/mut-abi-Cpd I); and 3) distance of Fe-C11/C7 should be less than 5 Å (systems mut-fer-C11/C7-Rest). The protonation states of all residues were confirmed by H++ Web server (Gordon et al., 2005), and the pocket residues were further carefully checked.

TABLE 1

TABLE 1. Illustration of simulation models.

Molecular Dynamics Simulation

Partial atomic charges of substrates were obtained from a restrained electrostatic potential (RESP) method at the HF/6-31G* level with Gaussian09 package (G09 A.01) (Frisch et al., 2009). The General Amber force field (GAFF) (Wang et al., 2005) was used to generate the force field parameters of small-molecule ligands (the abi and fer shown in Figure 1) by leap modules. Partial atomic charge and parameters of heme parameters of the resting state and Cpd I state were both from Shahrokh’s study (Shahrokh et al., 2012). Proteins were described with Amber ff14SB force field (Case et al., 2005). All systems were solvated by cuboid TIP3P (Jorgensen et al., 1983) water boxes with a 10 Å extended distance to the boundary of proteins. Sodium ions were utilized to neutralize all systems.

Three steps of minimization via steepest descent and conjugate gradient methods were taken to prepare the reasonable initial state. First, only water molecules and sodium ions were minimized with a restraint of 100 kcal⁻¹ mol Å⁻². Then only backbone atoms (C, CA, and N) of amino residues were restrained with a smaller force of 50 kcal⁻¹ mol Å⁻², and finally, all atoms were minimized without restraint. Subsequently, systems were heated from 0 to 300 K slowly for 80 ps and then kept at a temperature of 300 K for 20 ps under an NVT ensemble with restraint of 25 kcal⁻¹ mol Å⁻². And additional 300 ps of density equilibration simulations under the NPT ensemble were performed at a constant temperature of 300 K and a pressure of 1.0 atm. At last, all eight systems were implemented MD product and run for at least 400 ns in the NVT ensemble with a time step of 2 fs. Langevin dynamics was used to control the temperature with a collision frequency of 1.0 ps (Izaguirre et al., 2001). Long-range electrostatic interactions were computed with the particle mesh Ewald (PME) method (Essmann and Perera, 1995). Non-bounded cutoff distance was set to 10.0 Å. The SHAKE algorithm (Ryckaert et al., 1977) was applied to constrain hydrogen atoms. All systems were simulated under GPU-accelerated PMEMD program in Amber16 package (Salomon-Ferrer et al., 2013). The MD trajectories were analyzed by several programs including CPPTRAJ (Roe and Cheatham, 2013), Visual Molecular Dynamics (VMD) 1.9.3 (Humphrey et al., 1996), and POVME3.0 (Wagner et al., 2017). The root-measure-square deviations (RMSDs) and root-measure-square fluctuations (RMSFs) involved backbone atoms (C, CA, and N) of proteins were calculated. All the RMSD profiles of eight modeling systems were provided in Supplementary Figure S2. The K-means method was used for the cluster analysis. Principal component analysis (PCA), which recombines the original variables into a new set of unrelated variables to reflect principal information of original variables (Prompers and Bruschweiler, 2002), was performed to analyze the backbone atoms of pocket residues in wt/mut-abi-Rest systems. Pocket volumes of wt/mut-abi-Rest along the MD simulations were calculated by POVME3.0 (Wagner et al., 2017) and provided in Supplementary Figure S3A, and abietatriene was chosen to define the pocket with keyword DefinePocketByLigand, and the parameter GridSpacing was set to 1.0 Å. Pocket shape–based clustering was executed with respect to calculation results of each frame. In total, 800 frames were extracted from 400 ns MD trajectory at regular intervals for calculation of pocket volumes.

Besides, as a powerful technique of enhanced sampling, metadynamics (MTD) (Laio and Parrinello, 2002) was employed in the two systems (mut-abi-C11-Rest and mut-abi-C7-Rest) to map out free energy profiles for estimating the transformation of different ligand binding states, using NAMD 2.15 (Phillips et al., 2005) and Colvars module (Fiorin et al., 2013). The collective variables (CVs) of both systems were set to Fe-C11 distance and Fe-C7 distance, respectively. The Gaussian bias was used to fill the potential energy surface. The hills and widths were set to 0.1 kJ/mol and 0.05 Å, respectively. The Gaussian bias potential was added to the potential energy surface at every 2 ps. A well-tempered algorithm (Barducci et al., 2008) was utilized to decrease the Gaussian hills gradually with a bias factor of 6 to assist free energy profiles converging to definite values. The time step of MTD was set to 2 fs. In total, MTD simulation time of mut-abi-C11-Rest system was 190 and 550 ns for the mut-abi-C7-Rest system.

DFT Calculations

The reactivity of abietatriene was probed by DFT calculation on cluster models. In detail, the representative structures of cluster analysis from MD simulation trajectory of systems wt-abi-Cpd I and mut-abi-Cpd I were used to generate cluster models (Supplementary Figure S4). The model of Compound I was simplified by removing the side chains of porphyrin. The cysteine connected with porphyrin was replaced by the SCH₃^- group referring to previous modeling (Shaik et al., 2010). Amino residues were truncated with the boundary of α-C atoms which were then saturated with hydrogen atoms. All α-C atoms, hydrogen atoms linked to α-C, and methyl atoms linked to sulfur were frozen during calculations. All DFT calculations were performed with Gaussian 16 package (G16 A.03) (Frisch et al., 2016). All complexes were calculated with doublet states at the level of B3LYP (Lee et al., 1988; Becke, 1993), and empirical dispersion was corrected with D3BJ (Grimme et al., 2010). Geometry optimization and frequency calculations were performed at the B3LYP-D3/def2-SV(P) level (Weigend and Ahlrichs, 2005). The SMD model was applied with water (ε = 78.4) (Marenich et al., 2009). The transition-state (TS) structures were searched by Berny algorithm and verified with IRC calculations. Single-point energy calculation was performed at the level of B3LYP-D3/def2-TZVP, and zero-point energy (ZPE) correction was considered. The above theoretical methods and basis set size were used in previously computational studies of P450s (Dubey et al., 2016; Dubey et al., 2017; Su et al., 2019).

Results and Discussion

Different Pocket Dynamics of Wild-Type and Mutant CYP76AH1

The RMSF values of apo-(ligand-free) pocket residues showed minor differences between wild-type and mutant CYP76AH1 (Figure 2A). Differently, the corresponding RMSF values of holo-(abietatriene-bound) pocket residues in CYP76AH1^D301E,V479F (see Figure 2B) were much larger, whereas they became smaller in holo-state of wild-type CYP76AH1, suggesting that the active pocket becomes more flexible due to the ligand-induced fit effect and the D301E/V479F double mutant. Regarding the two mutant residues, as shown in Figure 2C, D301 and E301 overlap well, while the conformations of V479 and F479 do not overlap anymore. SRS1 is one of the six well-known substrate recognition sites (SRS) that serve as a key factor for substrate specificity and protein flexibility in previous studies (Gricman et al., 2015). Apparently, here, SRS1 was also found to be the most flexible region of the mutant CYP76AH1 based on the apparent RMSF value (see Figures 2B,D). All these are also confirmed by PCA of the abietatriene-bound active pocket, as shown in Figure 3. The PC1 and PC2 of wild-type system almost overlapped (Figure 3A), while the mutant system appeared to have a scattered distribution, mainly along PC1 (Figure 3B), similar to the PCA of SRS1 (Figure 3C), confirming that the flexibility of SRS1 is enhanced with the binding of the ligand (abietatriene) to the pocket of mutant CYP76AH1 (Figures 2C,D).

FIGURE 2

FIGURE 2. (A) Comparison of RMSF values of pocket residues between wt-apo-Rest and mut-apo-Rest systems. (B) Comparison of RMSF values of pocket residues between wt-abi-Rest and mut-abi-Rest systems. (C) Illustration of the active pocket and the two mutant residues in CYP76AH1. (D) Illustration of the flexible SRS1 region and the reactive center Fe.

FIGURE 3

FIGURE 3. PCA of pocket residues of systems wt-abi-Rest (A), mut-abi-Rest (B), and SRS1 residues of mut-abi-Rest system (C).

Based on pocket volume calculations and further cluster analysis, Figures 4A,B show the most representative active pockets of wild-type and mutant CYP76AH1 bound to abietatriene, respectively. Due to the double mutation (D301E and V479F), the pockets are much larger (showing larger volume and width), bringing enhanced adaptability to different ligands, which is very similar to the pocket analysis in P450s regarding the drug metabolism (Hendrychova et al., 2012). Therefore, we hypothesized that its novel catalytic function to generate multiple products might be due to a change in the pocket shape caused by the D301E/V479F mutant. To determine the exact effect of the mutant, we performed a comprehensive analysis of the binding sites (see infra).

FIGURE 4

FIGURE 4. Representative pocket shape of abietatriene-bound CYP76AH1: (A) wild type (wt-abi-Rest) and (B) the mutant (mut-abi-Rest). The pocket volume and its standard deviation are provided. The distance between two α-C atoms of L367 and V479 (F479 in mutant) is selected to characterize the width of pockets; the details are provided in Supplementary Figure S3B.

The Induced-Fit Binding Poses Change due to Mutants

As shown in Figure 1, in order to decipher the transition from abietatriene (abi) to ferruginol (fer, C12-hydroxy) in the wild type or multiproducts (C7-, C11-, and C12-oxidation) in the mutant, one needs to identify the near-attack binding poses as they trigger oxidative modifications on different carbon atoms. Therefore, MD trajectories for the wt-abi-Rest and mut-abi-Rest systems tracked the evolution of representative binding sites, as summarized in Figure 5. For wild-type CYP76AH1, it implies a potential reactive site on C12 of abietatriene, with a stable Fe-C12 distance of about 5 Å (Supplementary Figure S5A). Interestingly, as shown in Figure 5, not only the C12 site of abietatriene but also an additional reactive site on the C7 site of abietatriene was observed in CYP76AH1^D301E,V479F. Two attack sites of abietatriene were detected, one tending to Fe-C12 linkage and the other tending to Fe-C7 attack (see Supplementary Figure S5B) for a detailed distance profiles). As can be seen in Figure 5, the first is the C12-near-attack pose (i.e., state A in Figure 5), where the Fe-C12 distance is around 4.4 Å until 130 ns, maintaining an apparent “face-to-face” aromatic stacking interaction between the A ring of abietatriene and the benzene ring of Phe479. Then an intermediate state B was detected by flipping 90°, shifting to an “edge-to-face” aromatic stacking interaction. After that, it attained a new reactive binding pose (i.e., state C in Figure 5) with an Fe-C7 distance of about 5.1 Å. In conclusion, the induced ligand flip and conformational change of Phe479 would result in two reactive binding poses (C12-near-attack and C7-near-attack) in the mutant CYP76AH1, making the modification of the C7 position feasible. It should be noted that the value of the dihedral angle is ranging from 150° to 160° in the A state, and this value becomes more extensive between 100° and 160°, implying that both “face-to-face” and “edge-to-face” aromatic stacking interactions between the A ring and Phe479 would be present in the C state due to ligand flip and dynamic flip of the Phe479 aromatic ring (see Figure 5).

FIGURE 5

FIGURE 5. Dihedral angle between A-ring of the ligand and aromatic side chain of F479 along the MD trajectories in mut-abi-Rest and mut-fer-Rest systems. The corresponding representative binding poses are also given. The distance evolutions are supplemented in Supplementary Figures S5, S6.

In addition, an edge-to-face aromatic stacking interaction between the A-ring and Phe479 was also detected in the mut-fer-C7-Rest system (see Figure 5), where ferruginol adopted a stable C7 site near-attack binding pose with a stable Fe-C7 distance of about 4.8 Å (see Supplementary Figure S6A for details). This suggests that C12-hydroxylation from abietatriene to ferruginol may have occurred first in the mutant CYP76AH1, followed by further oxidation on C7, which is consistent with previous experiments that sugiol can be obtained in CYP76AH1^D301E,V479F. If 11-hydroxy-ferruginol is converted from ferruginol, the C11-near-attack binding pose of ferruginol is required. However, the C11-near-attack binding pose of ferruginol is not present in mutant CYP76AH1, which is confirmed by the expanded distance of Fe-C11 (∼7.5 Å) after 100 ns MD simulations of mut-fer-C11-Rest system (see Supplementary Figure S6B). To further confirm the above analysis of the binding poses, free energy profiles for the plausible binding pose transition were obtained by performing metadynamics simulations (Figure 6). In the mut-fer-C11-Rest system, the C11-near-attack binding pose (CV = 5.0 Å) lies at a metastable minimum, which is thermodynamically and kinetically facile to the other two more stable minimums (CV = 7.4 Å and 11.2 Å). In contrast, in the mut-fer-C7-Rest system, the dominant pose (CV = 4.7 Å) is located in a deep energy basin and is energetically difficult to transfer to the other minimum pose (CV = 11.4 Å).

FIGURE 6

FIGURE 6. Free energy profiles of binding poses transformation for systems (A) mut-fer-C11-Rest and (B) mut-fer-C7-Rest. The distance of Fe-C7 and Fe-C11 is defined as CVs which range from 2 to 15 Å.

In conclusion, the predominant binding poses of abietatriene and ferruginol in wild-type and mutant CYP76AH1 were comprehensively determined by detecting the different aromatic stacking patterns seen at the near-attack distance base on our classical MD simulations, and the spontaneous binding pose transition in mutant CYP76AH1 was clearly revealed by MTD simulations. Furthermore, our simulations demonstrate that the V479F mutation is responsible for the induced-fit binding pose transition. We then propose a plausible oxidative modification at the C7 position of abietatriene, which was not clear in previous studies. Furthermore, in the mutant CYP76AH1, C7- and C12(11)-near-attack binding poses were alternatively detected, and their conversion was also captured in the abietatriene bound state. The analysis of all these binding poses explains well the experimental results that the C7-oxidation products are detected in either the V479F single mutant or the D031E/V479F double mutant to CYP76AH1. Considering that only the C7-near-attack binding pose was observed in the mutant CYP76AH1 and no stable C11-near-attack binding pose was observed (see Figure 6), we inferred that all C11-hydroxylation products (see Figure 1) might not be derived from ferruginol but from the precursor abietatriene, which will be further demonstrated by DFT calculations (see infra).

The Plausible Reaction Pathways of Abietatriene in Wild-Type and Mutant CYP76AH1

Based on the above elaborative analysis of pocket dynamics and ligand (abietatriene and ferruginol) binding poses, we put forward a plausible reaction pathway for wild-type and mutant CYP76AH1 (D301E, V479F) to produce distinctive products, as shown in Figure 7, where both wild-type and mutant enzymes can catalyze the production of ferruginol from abietatriene. Sugiol can be obtained directly from ferruginol, and the formation of 11-hydroxy-ferruginol in this manner is prohibited. In contrast, according to recent studies on the mechanism of 1,2-naphalenediol production from naphthalene, 11-hydroxy-ferruginol may be obtained directly from abietatriene by epoxidation (Bao et al., 2019). Therefore, DFT calculations were performed to probe the reactivity of abietatriene in wild-type and mutant systems based on the truncated protein–ligand structure validated in the aforementioned MD simulations.

FIGURE 7

FIGURE 7. Plausible pathways of CYP76AH1 and CYP76AH1^D301E,V479F.

First, Figure 8 shows the electrophilic addition reaction of abietatriene, which is the rate-limiting step for aromatic oxidation. The energy barrier of 11.8 kcal/mol for the wild type is much lower than that of the mutant type, which has a reaction barrier of 16.8 kcal/mol. Considering the C12-near-attack binding poses of abietatriene and the critical distances (O-C12, Fe-O) and angles (Fe-O-C12) of the two transition state structures (wt-TS1 and mut-TS1) are very similar (Figures 9A,B), we deduce that the lower reaction barrier of the wild type may be due to the presence of two additional water molecules around the oxygen in the Cpd I reaction center (see Supplementary Figure S4 for details), whereas in our long-time MD simulations, the mutant model did not detect water molecules in this region because the larger side chain of Phe479 prevents water molecules from entering the pocket. According to previous studies (Altun et al., 2006a; Altun et al., 2006b; Kumar et al., 2011), water molecules near the oxygen of Cpd I stabilize the transition state by forming H-bonds and increasing the negative charge of the oxygen. This is consistent with our Mullikan charge prediction (see Figure 8), where the charge of oxygen in wt-TS1 (−0.646) is more negative than that in mut-TS1 (−0.437).

FIGURE 8

FIGURE 8. Reaction diagram [B3LYP-D3/def2TZVP//B3LYP-D3/def2SV(P)] for (A) wt-abi complex and (B) mut-abi complex. The values in brackets represent Mullikan charge. The cluster models are prepared from MD simulations of systems wt-abi-Cpd I and mut-abi-Cpd I as shown in Supplementary Figure S4. The distance of Fe-C12 and Fe-O-C12 angles is reasonable, since the distance is less than 4 Å, and the angle is fluctuating between 110° and 130° (Supplementary Figure S7). The corresponding optimized structures are given in Figure 9.

FIGURE 9

FIGURE 9. DFT optimized structures of along reaction for the two systems (A) wt-abi and (B) mut-abi.

The electrophilic addition reaction would then lead to a tetrahedral σ-complex (IM1 in Figure 8). Considering the energetically favorable character of the conversion from abietatriene to wt-IM1 in the wild type, ferruginol can be further efficiently and specifically yielded through a well-established NIH pathway (see Figure 8) (de Visser and Shaik, 2003). For the mutant system, the electrophilic addition reaction barrier becomes higher (16.8 kcal/mol) but is still kinetically reasonable, and the mut-IM1 can be yielded in view of the release heat (−7.5 kcal/mol). Interestingly, as shown in Supplementary Figure S1, the mutated residues (D301, F479) are indeed conserved and located at the same site of CYP76AH3, which also produces ferruginol but with a lower catalytic efficiency than wild-type CYP76AH1 (Mao et al., 2020). In this sense, these computational results explain the experimental results on the formation of the minor product ferruginol in the mutant CYP76AH1.

Finally, regarding the generation of C11-hydroxy products in the mutant enzyme, the reaction pathway from ferruginol was ruled out since ferruginol maintains an unfavorable attack distance (greater than 7 Å, i.e., fer-C11 binding pose), as shown in Figure 7. Here, we refer to a previously established reaction mechanism (Bao et al., 2019) that naphthalene can be converted to 1,2-naphalenediol, a double ortho-hydroxy product similar to 11-hydroxyl ferruginol, via epoxidation and then opening by OH radicals and water molecules (see Supplementary Figure S8 for details). Since a reactive C12-near-attack binding pose of abietatriene was detected in our MD simulations (provided in Figure 7), we propose that oxygen will attack C11 to form the epoxidation product (mut-IM2). As shown in Figure 8B, DFT calculations show that this process is very facile with a low activation barrier (10.1 kcal/mol) and heat release (2 kcal/mol). Then mut-IM2 can be converted to the final product 11-hydroxy-ferruginol via a ring-opening reaction according to a previous study (Bao et al., 2019), as shown in Figure 8B, Supplementary Figure S8.

Conclusion

In this work, the catalytic plasticity of the CYP76AH1 enzyme was investigated by multiscale simulations. The reasons for the different reaction pathways of the wild-type and D301E/V479F mutants were revealed, with the former being highly specific for a single product and the latter being very promiscuous for multiple products. It was found that the product promiscuity of mutant CYP76AH1 was mainly attributed to the presence of an alternative (C12-vs C7-near-attack) binding pose for the key intermediate abietatriene, whereas a unique C12-binding pose for abietatriene exists in the wild-type enzyme. Meanwhile, the strong pocket plasticity of CYP76AH1 was confirmed by highly adaptive pocket dynamics in response to ligand-induced conformational changes, which was also observed in many P450s, especially those involved in drug metabolism. All these findings provide clues not only for understanding oxidative modifications of tanshinone-like norditerpenoid natural products but also for redesigning CYP76 family enzyme function.

Data Availability Statement

The original contributions presented in the study are included in the article/Supplementary Material; further inquiries can be directed to the corresponding author.

Author Contributions

All authors listed have made a substantial, direct, and intellectual contribution to the work and approved it for publication. RW supervised the whole project; YQ performed the most calculation; HD assisted the DFT calculation; and YZ assisted the sequence/structure alignment and manuscript writing.

Funding

This work was supported by the National Natural Science Foundation of China (21773313) and Sun Yat-sen Univerisity (20ykzd13).

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

We thank the Guangzhou and Shenzhen Supercomputer Center for providing computational source.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fchem.2021.689731/full#supplementary-material

References

Altun, A., Guallar, V., Friesner, R. A., Shaik, S., and Thiel, W. (2006a). The Effect of Heme Environment on the Hydrogen Abstraction Reaction of Camphor in P450camCatalysis: A QM/MM Study. J. Am. Chem. Soc. 128, 3924–3925. doi:10.1021/ja058196w

PubMed Abstract | CrossRef Full Text | Google Scholar

Altun, A., Shaik, S., and Thiel, W. (2006b). Systematic QM/MM Investigation of Factors that Affect the Cytochrome P450-Catalyzed Hydrogen Abstraction of Camphor. J. Comput. Chem. 27, 1324–1337. doi:10.1002/jcc.20398

PubMed Abstract | CrossRef Full Text | Google Scholar

Bao, L., Liu, W., Li, Y., Wang, X., Xu, F., Yang, Z., et al. (2019). Carcinogenic Metabolic Activation Process of Naphthalene by the Cytochrome P450 Enzyme 1B1: A Computational Study. Chem. Res. Toxicol. 32, 603–612. doi:10.1021/acs.chemrestox.8b00297

PubMed Abstract | CrossRef Full Text | Google Scholar

Barducci, A., Bussi, G., and Parrinello, M. (2008). Well-tempered Metadynamics: A Smoothly Converging and Tunable Free-Energy Method. Phys. Rev. Lett. 100. doi:10.1103/PhysRevLett.100.020603

CrossRef Full Text | Google Scholar

Bathe, U., and Tissier, A. (2019). Cytochrome P450 Enzymes: A Driving Force of Plant Diterpene Diversity. Phytochemistry 161, 149–162. doi:10.1016/j.phytochem.2018.12.003

PubMed Abstract | CrossRef Full Text | Google Scholar

Becke, A. D. (1993). Density‐functional Thermochemistry. III. The Role of Exact Exchange. J. Chem. Phys. 98, 5648–5652. doi:10.1063/1.464913

CrossRef Full Text | Google Scholar

Bonomo, S., Jørgensen, F. S., and Olsen, L. (2017). Mechanism of Cytochrome P450 17A1-Catalyzed Hydroxylase and Lyase Reactions. J. Chem. Inf. Model. 57, 1123–1133. doi:10.1021/acs.jcim.6b00759

PubMed Abstract | CrossRef Full Text | Google Scholar

Case, D. A., Cheatham, T. E., Darden, T., Gohlke, H., Luo, R., Merz, K. M., et al. (2005). The Amber Biomolecular Simulation Programs. J. Comput. Chem. 26, 1668–1688. doi:10.1002/jcc.20290

PubMed Abstract | CrossRef Full Text | Google Scholar

De Visser, S. P., and Shaik, S. (2003). A Proton-Shuttle Mechanism Mediated by the Porphyrin in Benzene Hydroxylation by Cytochrome P450 Enzymes. J. Am. Chem. Soc. 125, 7413–7424. doi:10.1021/ja034142f

PubMed Abstract | CrossRef Full Text | Google Scholar

Dong, Y., Morris-Natschke, S. L., and Lee, K.-H. (2011). Biosynthesis, Total Syntheses, and Antitumor Activity of Tanshinones and Their Analogs as Potential Therapeutic Agents. Nat. Prod. Rep. 28, 529–542. doi:10.1039/c0np00035c

PubMed Abstract | CrossRef Full Text | Google Scholar

Dubey, K. D., Wang, B., and Shaik, S. (2016). Molecular Dynamics and QM/MM Calculations Predict the Substrate-Induced Gating of Cytochrome P450 BM3 and the Regio- and Stereoselectivity of Fatty Acid Hydroxylation. J. Am. Chem. Soc. 138, 837–845. doi:10.1021/jacs.5b08737

PubMed Abstract | CrossRef Full Text | Google Scholar

Dubey, K. D., Wang, B., Vajpai, M., and Shaik, S. (2017). MD Simulations and QM/MM Calculations Show that Single-Site Mutations of Cytochrome P450BM3 Alter the Active Site's Complexity and the Chemoselectivity of Oxidation without Changing the Active Species. Chem. Sci. 8, 5335–5344. doi:10.1039/c7sc01932g

PubMed Abstract | CrossRef Full Text | Google Scholar

Essmann, U., Perera, L., Berkowitz, M. L., Darden, T., Lee, H., and Pedersen, L. G. (1995). A Smooth Particle Mesh Ewald Method. J. Chem. Phys. 103, 8577–8593. doi:10.1063/1.470117

CrossRef Full Text | Google Scholar

Fiorin, G., Klein, M. L., and Hénin, J. (2013). Using Collective Variables to Drive Molecular Dynamics Simulations. Mol. Phys. 111, 3345–3362. doi:10.1080/00268976.2013.813594

CrossRef Full Text | Google Scholar

Frisch, M., Trucks, G., Schlegel, H., Scuseria, G., Robb, M., Cheeseman, J., et al. (2009). GAUSSIAN09. Revision A. 01. Wallingford, CT: Gaussian Inc.

Frisch, M., Trucks, G., Schlegel, H., Scuseria, G., Robb, M., Cheeseman, J., et al. (2016). GAUSSIAN16. Revision A. 03. Wallingford, CT: Gaussian Inc.

Gao, W., Hillwig, M. L., Huang, L., Cui, G., Wang, X., Kong, J., et al. (2009). A Functional Genomics Approach to Tanshinone Biosynthesis Provides Stereochemical Insights. Org. Lett. 11, 5170–5173. doi:10.1021/ol902051v

PubMed Abstract | CrossRef Full Text | Google Scholar

Gordon, J. C., Myers, J. B., Folta, T., Shoja, V., Heath, L. S., and Onufriev, A. (2005). H++: a Server for Estimating pKas and Adding Missing Hydrogens to Macromolecules. Nucleic Acids Res. 33, W368–W371. doi:10.1093/nar/gki464

PubMed Abstract | CrossRef Full Text | Google Scholar

Gricman, Ł., Vogel, C., and Pleiss, J. (2015). Identification of Universal Selectivity-Determining Positions in Cytochrome P450 Monooxygenases by Systematic Sequence-Based Literature Mining. Proteins 83, 1593–1603. doi:10.1002/prot.24840

PubMed Abstract | CrossRef Full Text | Google Scholar

Grimme, S., Antony, J., Ehrlich, S., and Krieg, H. (2010). A Consistent and Accurate Ab Initio Parametrization of Density Functional Dispersion Correction (DFT-D) for the 94 Elements H-Pu. J. Chem. Phys. 132, 154104. doi:10.1063/1.3382344

PubMed Abstract | CrossRef Full Text | Google Scholar

Gu, M., Wang, M., Guo, J., Shi, C., Deng, J., Huang, L., et al. (2019). Crystal Structure of CYP76AH1 in 4-PI-Bound State from Salvia Miltiorrhiza. Biochem. Biophysical Res. Commun. 511, 813–819. doi:10.1016/j.bbrc.2019.02.103

PubMed Abstract | CrossRef Full Text | Google Scholar

Guo, J., Ma, X., Cai, Y., Ma, Y., Zhan, Z., Zhou, Y. J., et al. (2016). Cytochrome P450 Promiscuity Leads to a Bifurcating Biosynthetic Pathway for Tanshinones. New Phytol. 210, 525–534. doi:10.1111/nph.13790

PubMed Abstract | CrossRef Full Text | Google Scholar

Guo, J., Zhou, Y. J., Hillwig, M. L., Shen, Y., Yang, L., Wang, Y., et al. (2013). CYP76AH1 Catalyzes Turnover of Miltiradiene in Tanshinones Biosynthesis and Enables Heterologous Production of Ferruginol in Yeasts. Proc. Natl. Acad. Sci. 110, 12108–12113. doi:10.1073/pnas.1218061110

PubMed Abstract | CrossRef Full Text | Google Scholar

Hendrychova, T., Berka, K., Navratilova, V., Anzenbacher, P., and Otyepka, M. (2012). Dynamics and Hydration of the Active Sites of Mammalian Cytochromes P450 Probed by Molecular Dynamics Simulations. Cdm 13, 177–189. doi:10.2174/138920012798918408

PubMed Abstract | CrossRef Full Text | Google Scholar

Humphrey, W., Dalke, A., and Schulten, K. (1996). VMD: Visual Molecular Dynamics. J. Mol. Graphics 14, 33–38. doi:10.1016/0263-7855(96)00018-5

PubMed Abstract | CrossRef Full Text | Google Scholar

Izaguirre, J. A., Catarello, D. P., Wozniak, J. M., and Skeel, R. D. (2001). Langevin Stabilization of Molecular Dynamics. J. Chem. Phys. 114, 2090–2098. doi:10.1063/1.1332996

CrossRef Full Text | Google Scholar

Jorgensen, W. L., Chandrasekhar, J., Madura, J. D., Impey, R. W., and Klein, M. L. (1983). Comparison of Simple Potential Functions for Simulating Liquid Water. J. Chem. Phys. 79, 926–935. doi:10.1063/1.445869

CrossRef Full Text | Google Scholar

Kumar, D., Altun, A., Shaik, S., and Thiel, W. (2011). Water as Biocatalyst in Cytochrome P450. Faraday Discuss. 148, 373–383. doi:10.1039/c004950f

PubMed Abstract | CrossRef Full Text | Google Scholar

Laio, A., and Parrinello, M. (2002). Escaping Free-Energy Minima. Proc. Natl. Acad. Sci. 99, 12562–12566. doi:10.1073/pnas.202427399

PubMed Abstract | CrossRef Full Text | Google Scholar

Lee, C., Yang, W., and Parr, R. G. (1988). Development of the Colle-Salvetti Correlation-Energy Formula into a Functional of the Electron Density. Phys. Rev. B 37, 785–789. doi:10.1103/PhysRevB.37.785

PubMed Abstract | CrossRef Full Text | Google Scholar

Leth, R., Ercig, B., Olsen, L., and Jørgensen, F. S. (2019). Both Reactivity and Accessibility Are Important in Cytochrome P450 Metabolism: A Combined DFT and MD Study of Fenamic Acids in BM3 Mutants. J. Chem. Inf. Model 59, 743–753. doi:10.1021/acs.jcim.8b00750

PubMed Abstract | CrossRef Full Text | Google Scholar

Lonsdale, R., Houghton, K. T., Żurek, J., Bathelt, C. M., Foloppe, N., De Groot, M. J., et al. (2013). Quantum Mechanics/molecular Mechanics Modeling of Regioselectivity of Drug Metabolism in Cytochrome P450 2C9. J. Am. Chem. Soc. 135, 8001–8015. doi:10.1021/ja402016p

PubMed Abstract | CrossRef Full Text | Google Scholar

Mao, Y., Ma, Y., Chen, T., Ma, X., Xu, Y., Bu, J., et al. (2020). Functional Integration of Two CYP450 Genes Involved in Biosynthesis of Tanshinones for Improved Diterpenoid Production by Synthetic Biology. ACS Synth. Biol. 9, 1763–1770. doi:10.1021/acssynbio.0c00136

PubMed Abstract | CrossRef Full Text | Google Scholar

Marenich, A. V., Cramer, C. J., and Truhlar, D. G. (2009). Universal Solvation Model Based on Solute Electron Density and on a Continuum Model of the Solvent Defined by the Bulk Dielectric Constant and Atomic Surface Tensions. J. Phys. Chem. B 113, 6378–6396. doi:10.1021/jp810292n

PubMed Abstract | CrossRef Full Text | Google Scholar

Meunier, B., De Visser, S. P., and Shaik, S. (2004). Mechanism of Oxidation Reactions Catalyzed by Cytochrome P450 Enzymes. Chem. Rev. 104, 3947–3980. doi:10.1021/cr020443g

PubMed Abstract | CrossRef Full Text | Google Scholar

Phillips, J. C., Braun, R., Wang, W., Gumbart, J., Tajkhorshid, E., Villa, E., et al. (2005). Scalable Molecular Dynamics with NAMD. J. Comput. Chem. 26, 1781–1802. doi:10.1002/jcc.20289

PubMed Abstract | CrossRef Full Text | Google Scholar

Prompers, J. J., and Brüschweiler, R. (2002). Dynamic and Structural Analysis of Isotropically Distributed Molecular Ensembles. Proteins 46, 177–189. doi:10.1002/prot.10025

PubMed Abstract | CrossRef Full Text | Google Scholar

Robertson, A. L., Holmes, G. R., Bojarczuk, A. N., Burgon, J., Loynes, C. A., Chimen, M., et al. (2014). A Zebrafish Compound Screen Reveals Modulation of Neutrophil Reverse Migration as an Anti-inflammatory Mechanism. Sci. Translational Med. 6, 225ra29. doi:10.1126/scitranslmed.3007672

PubMed Abstract | CrossRef Full Text | Google Scholar

Roe, D. R., and Cheatham, T. E. (2013). PTRAJ and CPPTRAJ: Software for Processing and Analysis of Molecular Dynamics Trajectory Data. J. Chem. Theor. Comput. 9, 3084–3095. doi:10.1021/ct400341p

PubMed Abstract | CrossRef Full Text | Google Scholar

Ryckaert, J.-P., Ciccotti, G., and Berendsen, H. J. C. (1977). Numerical Integration of the Cartesian Equations of Motion of a System with Constraints: Molecular Dynamics of N-Alkanes. J. Comput. Phys. 23, 327–341. doi:10.1016/0021-9991(77)90098-5

CrossRef Full Text | Google Scholar

Salomon-Ferrer, R., Götz, A. W., Poole, D., Le Grand, S., and Walker, R. C. (2013). Routine Microsecond Molecular Dynamics Simulations with AMBER on GPUs. 2. Explicit Solvent Particle Mesh Ewald. J. Chem. Theor. Comput. 9, 3878–3888. doi:10.1021/ct400314y

PubMed Abstract | CrossRef Full Text | Google Scholar

Sawai, S., and Saito, K. (2011). Triterpenoid Biosynthesis and Engineering in Plants. Front. Plant Sci. 2, 25. doi:10.3389/fpls.2011.00025

PubMed Abstract | CrossRef Full Text | Google Scholar

Scheler, U., Brandt, W., Porzel, A., Rothe, K., Manzano, D., Božić, D., et al. (2016). Elucidation of the Biosynthesis of Carnosic Acid and its Reconstitution in Yeast. Nat. Commun. 7, 12942. doi:10.1038/ncomms12942

PubMed Abstract | CrossRef Full Text | Google Scholar

Shahrokh, K., Orendt, A., Yost, G. S., and Cheatham, T. E. (2012). Quantum Mechanically Derived AMBER-Compatible Heme Parameters for Various States of the Cytochrome P450 Catalytic Cycle. J. Comput. Chem. 33, 119–133. doi:10.1002/jcc.21922

PubMed Abstract | CrossRef Full Text | Google Scholar

Shaik, S., Cohen, S., Wang, Y., Chen, H., Kumar, D., and Thiel, W. (2010). P450 Enzymes: Their Structure, Reactivity, and Selectivity-Modeled by QM/MM Calculations. Chem. Rev. 110, 949–1017. doi:10.1021/cr900121s

PubMed Abstract | CrossRef Full Text | Google Scholar

Su, H., Wang, B., and Shaik, S. (2019). Quantum-Mechanical/Molecular-Mechanical Studies of CYP11A1-Catalyzed Biosynthesis of Pregnenolone from Cholesterol Reveal a C-C Bond Cleavage Reaction that Occurs by a Compound I-Mediated Electron Transfer. J. Am. Chem. Soc. 141, 20079–20088. doi:10.1021/jacs.9b08561

PubMed Abstract | CrossRef Full Text | Google Scholar

Sun, W., Xue, H., Liu, H., Lv, B., Yu, Y., Wang, Y., et al. (2020). Controlling Chemo- and Regioselectivity of a Plant P450 in Yeast Cell toward Rare Licorice Triterpenoid Biosynthesis. ACS Catal. 10, 4253–4260. doi:10.1021/acscatal.0c00128

CrossRef Full Text | Google Scholar

Wagner, J. R., Sørensen, J., Hensley, N., Wong, C., Zhu, C., Perison, T., et al. (2017). POVME 3.0: Software for Mapping Binding Pocket Flexibility. J. Chem. Theor. Comput. 13, 4584–4592. doi:10.1021/acs.jctc.7b00500

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, J., Wolf, R. M., Caldwell, J. W., Kollman, P. A., and Case, D. A. (2005). Junmei Wang, Romain M. Wolf, James W. Caldwell, Peter A. Kollman, and David A. Case, "Development and Testing of a General Amber Force field"Journal of Computational Chemistry (2004) 25(9) 1157-1174. J. Comput. Chem. 26, 114. doi:10.1002/jcc.20145

CrossRef Full Text | Google Scholar

Weigend, F., and Ahlrichs, R. (2005). Balanced Basis Sets of Split Valence, Triple Zeta Valence and Quadruple Zeta Valence Quality for H to Rn: Design and Assessment of Accuracy. Phys. Chem. Chem. Phys. 7, 3297–3305. doi:10.1039/b508541a

PubMed Abstract | CrossRef Full Text | Google Scholar

Yang, S., Demars, M. D., Grandner, J. M., Olson, N. M., Anzai, Y., Sherman, D. H., et al. (2020a). Computational-Based Mechanistic Study and Engineering of Cytochrome P450 MycG for Selective Oxidation of 16-Membered Macrolide Antibiotics. J. Am. Chem. Soc. 142, 17981–17988. doi:10.1021/jacs.0c04388

PubMed Abstract | CrossRef Full Text | Google Scholar

Yang, W., Chen, X., Li, Y., Guo, S., Wang, Z., and Yu, X. (2020b). Advances in Pharmacological Activities of Terpenoids. Nat. Product. Commun. 15, 1934578X2090355. doi:10.1177/1934578x20903555

CrossRef Full Text | Google Scholar

Zhang, X., Peng, Y., Zhao, J., Li, Q., Yu, X., Acevedo-Rocha, C. G., et al. (2020). Bacterial Cytochrome P450-Catalyzed Regio- and Stereoselective Steroid Hydroxylation Enabled by Directed Evolution and Rational Design. Bioresour. Bioproc. 7. doi:10.1186/s40643-019-0290-4

CrossRef Full Text | Google Scholar

Zheng, X., Li, P., and Lu, X. (2019). Research Advances in Cytochrome P450-Catalysed Pharmaceutical Terpenoid Biosynthesis in Plants. J. Exp. Bot. 70, 4619–4630. doi:10.1093/jxb/erz203

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhou, L., Zuo, Z., and Chow, M. S. S. (2005). Danshen: An Overview of its Chemistry, Pharmacology, Pharmacokinetics, and Clinical Use. J. Clin. Pharmacol. 45, 1345–1359. doi:10.1177/0091270005282630

PubMed Abstract | CrossRef Full Text | Google Scholar

Zi, J., and Peters, R. J. (2013). Characterization of CYP76AH4 Clarifies Phenolic Diterpenoid Biosynthesis in the Lamiaceae. Org. Biomol. Chem. 11, 7650–7652. doi:10.1039/c3ob41885e

PubMed Abstract | CrossRef Full Text | Google Scholar