- Email: [email protected]
Reaction mechanism of isoflavone O-methyltransferase: A theoretical investigation
Chemical Physics Letters 501 (2011) 502–507
Contents lists available at ScienceDirect
Chemical Physics Letters journal homepage: www.elsevier.com/locate/cplett
Reaction mechanism of isoflavone O-methyltransferase: A theoretical investigation Feng-Chao Cui, Xiao-Liang Pan, Jing-Yao Liu ⇑ State Key Laboratory of Theoretical and Computational Chemistry, Institute of Theoretical Chemistry, Jilin University, Changchun 130023, China
a r t i c l e
i n f o
Article history: Received 10 June 2010 In final form 14 November 2010 Available online 17 November 2010
a b s t r a c t The methyl-transfer reaction mechanism catalyzed by Isoflavone O-methyl-transferase (IOMT) and the roles of several residues around the active site are investigated by employing density functional method. The calculations confirm that the proton transfer from daidzein to His257 occurs barrierlessly, and the methyl group is transferred from S-adenosyl-L-methionine (SAM) to phenolate ion in a single SN 2 step with a barrier of 17.0 kcal/mol, consistent with experimental value. Glu318 and Asp288 play important roles in lowering the reaction barrier. Ó 2010 Published by Elsevier B.V.
1. Introduction Isoflavone O-methyltransferase (IOMT) belongs to S-adenosyl-Lmethionine (SAM) dependent plant natural product O-methyltransferase (OMT) involved in secondary metabolism [1] and is a key enzyme for the biosynthesis of the antimicrobial pterocarpan phytoalexin medicarpin, which is regarded as an important component of disease resistance response in alfalfa (Medicago sativa L.) [2,3]. As shown in Figure 1, in vivo, IOMT catalyzes the transfer of a methyl group from SAM to the B-ring 40 -hydroxyl group of the isoflavone daidzein, leading to the formation of S-adenosylhomocysteine (SAH) and formononetin, and the latter is an essential intermediate for synthesizing medicarpin in the isoflavonoid branch of the phenylpropanoid pathway [4]. However, in vitro, the A-ring 7-hydroxyl group of the isoflavone daidzein is preferentially methylated, yielding the isoformononetin, which is rarely found in plants and has no known physiological function [1,2]. The crystal structure of IOMT in complex with the products SAH and isoformononetin (40 -hydroxy-7-methoxyisoflavone) was first determined at 1.4 Å in 2001 by Noel et al. [1]. Figure 2 shows the major feature of the active site of this structure. Asn258, Asn310, Met168 and Met311 help in orienting the daidzein. Glu318 and Asp288 bracket the catalytic His257 residue. On the basis of this crystal structure, sequence alignments with the large family of plant OMTs and the mutational data, Noel et al. [1] put forward that similar to other methyltransferases [5–10], the putative catalytic base His257 is responsible for abstracting the proton from the 7-hydroxyl group of the A-ring of daidzein, and then the strong nucleophile (phenolate ion) generated on the deprotonated 7-hydroxyl of daidzein is now able to abstract the reactive methyl group from the positively charged SD of SAM in an SN 2 reaction [11,12], yielding the isoformononetin and SAH products [1]. They
⇑ Corresponding author. E-mail address: [email protected] (J.-Y. Liu). 0009-2614/$ - see front matter Ó 2010 Published by Elsevier B.V. doi:10.1016/j.cplett.2010.11.045
also speculated that the catalytic role of Glu318 is likely to proceed through the His–Glu charge-relay mechanism [13,14], that is, involving the proton transfer between Glu318 and His257, similar to protein arginine N-methyltransferase [8]. The purposes of this paper are to investigate the suggested mechanism and to examine the energetic feasibility of the methyl transfer reaction by employing the hybrid density functional theory (DFT) method B3LYP [15–20], which has been successfully used to study this kind of reaction [21–24]. The catalytic roles of several residues in the active site are also discussed by using several quantum chemical models with different sizes.
2. Computational details In the present study, quantum chemistry calculations were performed using the density functional theory (DFT) with Becke’s three-parameter (B3) [15–18] exchange functional along with the Lee–Yang–Parr (LYP) [19,20] nonlocal correlation functional (B3LYP). The geometries were optimized by using the 6-31 G(d,p) basis set. Hessians were calculated at the same level to confirm the nature of the stationary points, that is, each equilibrium species possesses all real frequencies, whereas the transition state possesses one and only one imaginary frequency. For each transition state, intrinsic reaction coordinate (IRC) calculations were performed to guarantee its correct connection to the designated local minima. To evaluate solvation effects of the missing protein environment, single-point calculations on the B3LYP/6-31 G (d,p) optimized geometries were performed using the conductor-like solvation model COSMO [25–28] at the same level. In the model, a cavity around the system is surrounded by polarizable dielectric continuum with the dielectric constant of 4, which is usually chosen for protein surrounding [29]. On the basis of these optimized geometries, single-point energies for all the stationary points were calculated at the B3LYP/6-311+G(2 d,2p) level to obtain more accurate energies. Some centers were fixed to preserve the positions of
F.-C. Cui et al. / Chemical Physics Letters 501 (2011) 502–507
HO 7
where pK cal a ðCPCMÞ value was estimated by using the thermodynamics cycles mentioned above. The above empirical correction formula was obtained by linear correlation relationship, and it has been tested that the calculated pK a values using Eq. (1) are very close to the corresponding experimental data. More computational details have been discussed in the literature [30]. In the present study, we performed the same-level calculations as those done in Ref.30 to estimate the pK cal a value. The gas-phase free energies of the acids and their conjugate bases were obtained by geometry optimizations and harmonic vibrational frequency calculations at the B3LYP/6-31+G* level of theory, and the solvation free energy of each species was calculated as a single point using conductor-like polarizable continuum model (CPCM) [31] at the HF/6-31+G* level based on HF/6-31+G*-optimized gas-phase structures. All computations were performed by using the GAUSSIAN03 package of programs.
O
A
daidzein O
B
4'
OH
In vivo HO
In vitro O
O CH3
O
O
formononetin
HO
CH3
O
O
isoformononetin
503
OH
O
3. Results and discussion O
medicarpin
O
CH3
Figure 1. O-methylation of the isoflavone daidzein catalyzed by Isoflavone Omethyltransferase (IOMT) in vivo and in vitro.
Figure 2. X-ray crystal structure of the active site of IOMT complex with SAH and isoformononetin.
the active site residues as similar as those in the crystal structure during the optimization. This constraint in geometry optimizations gives rise to a few imaginary frequencies, typically on the order of 10i—20i cm1 . These frequencies do not contribute significantly to the zero-point energies and can thus be ignored. In addition, it is found that the proton-transfer from daidzein to His257 occurs barrierlessly during the optimization process of reactant when Glu318 was included in the models. To provide the likely explanation for this barrierless proton-transfer process, the pK a values of the 7-hydroxyl of daidzein and the His257 N e in the presence or absence of Glu318 were estimated by employing thermodynamics cycles, which have been applied comprehensively for pK a calculations. This approach indirectly determines free-energy change for proton dissociation in aqueous via calculating free-energy change for proton dissociation in gas phase and the solvation free energies of the general acid and dissociated species. Here, we employed the following empirical correction equation proposed by Qian et al. [30] to evaluate pK a values:
pK a ¼ 0:418pK cal a ðCPCMÞ þ 2:27
ð1Þ
In the present study, several quantum chemical models of the IOMT active site with different sizes, ranging from 52 to 148 atoms, have been constructed based on the X-ray crystal structure of IOMT in complex with the products S-adenosyl-L-homocysteine (SAH) and isoformononetin (40 -hydroxy-7-methoxyisoflavone) (PDB code 1FP2) [1]. The substrate daidzein was built by displacing the methyl of 7-methoxy with hydrogen and SAH was converted to SAM by adding a methyl group at the S atom position. All hydrogen atoms were added manually. To keep the various groups in place to resemble the crystal structure as much as possible, the truncation points of the SAM cofactor and various active site residues were kept frozen to their X-ray positions during the geometry optimization. These fixed atoms are indicated by asterisks in the figures given later. In order to shed more light on the detailed methyl transfer process catalyzed by IOMT and evaluate the effects of the various residues at the active site of IOMT, we performed the calculations from smaller to larger models. The simplest model of the IOMT active site (called Model A) consists of 52 atoms and just contains the substrate daidzein, SAM cofactor and His257, which are essential for the catalytic reaction. The B-ring of the substrate daidzein was truncated and represented by a methyl. SAM was modeled by truncating the cofactor two carbons away on each side of the sulfur center. It has been substantiated by Himo et al. that this model of SAM is sufficient to reproduce the properties of the SD —C E bond [21–23]. His257 regarded as a base was represented by a methyl imidazole ring. The optimized structures of the reactant and transition state of Model are shown in Figure 3. The calculated potential energy profiles using all models with different sizes in the present investigation are displayed in Figure 4. As seen from TS-A in Figure 3, the distances of the breaking SD —C E bond and the forming C E —OD bond are 2.43 and 1.96 Å, respectively, and the transferred proton slightly moves away from the OD atom of the substrate daidzein toward the N e atom of His257, with OD —H and H—N e bond distances of 1.09 and 1.46 Å, respectively. This shows that the proton transfer from the substrate daidzein to His257 and the methyl cation transfer from SAM to the substrate daidzein occur in one concerted asynchronous step via an SN 2 reaction mechanism. The barrier of Model A is calculated to be 29.6 kcal/mol including solvation effects with ¼ 4, and the reaction is slightly endothermic by 1.5 kcal/mol. Clearly, the barrier of the simplest model with solvation effects is apparently higher than the experimental value of ca.19.8 kcal/mol estimated from the rate constantðkcat ¼ 0:03 s Þ [2] by using classical transition state theory. This poor agreement between theory and experiment indicates that some important residues for lowering the barrier of the catalytic reaction are missing and need to be included in the quantum chemical model.
504
F.-C. Cui et al. / Chemical Physics Letters 501 (2011) 502–507
Figure 3. Optimized structures of the reactant for Model A and transition states for Models A–C. Distances are in angstroms. Asterisks indicate that atoms are fixed to their Xray positions.
Relative energy (kcal/mol)
TS-C TS-A TS-B
30.1(31.0) 29.6(21.8) 29.5(31.0)
TS-D 23.8(17.6) TS-E 19.1(12.2) TS-F 17.0(14.6)
0.0
Pro-A
R-(A-F)
1.5(-2.8)
Pro-C -2.6(-2.9) Pro-B -4.1(-2.7) Pro-D -7.2(-18.7) Pro-E -7.4(-22.8) Pro-F -12.8(-25.6)
Reaction coordinate Figure 4. Computed potential energy profiles of all models in the solution and in gas phase (values in parentheses).
Two somewhat larger models (called Models B and C) of the IOMT active site were constructed. In Model B, the amide side chains of Asn258 and Asn310, which form hydrogen bonds to the substrate daidzein, were added based on Model A. In addition to parts included in Model B, Model C contains parts of side chains of Met168 and Met311, which sandwich the A ring of the substrate daidzein. The optimized TS structures for Models B and C (TS-B and TS-C) are also displayed in Figure 3. At both TSs, the critical SD —C E and C E —OD bond distances for the methyl transfer change a little compared to those of TS-A, while the OD —H and H—N e bond distances for proton transfer have larger variation, for example, H—N e bond distances are 1.08 and 1.07 Å in TS-B and TS-C (as shown in Figure 3), respectively. These results indicate that both models exhibit a single-step mechanism, which involves concerted the methyl-transfer and proton-transfer processes, similar to Model A, however, we find that the proton transfer takes place much earlier in these two models (almost finished at the two TSs) compared to that in Model A. On inclusion of solvent effects, the
calculated barriers of Models B and C are 29.5 and 30.1 kcal/mol, respectively, very close to that of Model A. Therefore, it appears that these four residues, which form hydrogen bonds or steric constraint to the substrate daidzein, may be helpful in preventing larger spatial reorganization of the substrate daidzein and the active site residues so as to keep them close to the crystal structure, while they don’t have effect on reducing the reaction barrier. Based on Model C, we constructed a larger model (called Model D), consisting of 104 atoms. This model contains a part of side chain of Glu318 (modeled by an acetate), which makes hydrogen bonding interaction to the N d of His257, and three water molecules around Glu318 presented in the crystal structure. The optimized structures of the reactant and transition state of Model D are shown in Figure 5. It is interesting to find that the proton of the substrate daidzein becomes unstable with inclusion of Glu318, and the optimization of the reactant always leads to the proton transfer to His257 (R-D in Figure 5). That is to say, in Model D, the barrierless proton-transfer from the substrate daidzein to
F.-C. Cui et al. / Chemical Physics Letters 501 (2011) 502–507
505
Figure 5. Optimized the structures of the reactants and transition-states for Models D and E.
His257 occurs firstly while the enzyme-substrate complex was formed, thereby generating two charged groups at the active site, i.e., the phenolate anion of the substrate daidzein and protonated His257, which is different from the behavior happening in the former three models. The likely reason of this phenomenon will be discussed in the later section of the text. Then, the methyl as a cation is transferred to the substrate daidzein in an SN 2 reaction, with the critical SD —C E and C E —OD bond distances being 2.29 and 2.13 Å, respectively. It is seen from Figure 4 that adding the residue Glu318 in the model results in a dramatic decrease of the barrier, about 6.3 kcal/mol lower than that of Model C, and an increase of the exothermicity by 4.6 kcal/mol, to 7.2 kcal/mol. It is thus conceivable to consider that Glu318 plays an important role in lowering the energetics of the catalytic reaction. However, the calculated barrier of Model D is still 4.0 kcal/mol higher than the experimental value (23.8 vs. 19.8 kcal/mol), implying more groups need to be added successively. In addition, because the hydrogen bonding interactions between the crystal water molecules and Glu318 may make the anionic Glu318 be stabilized, no proton-transfer from His257 to Glu318 occurs during the process of the catalytic reaction, implying that the catalytic role of Glu318 would proceed through the electrostatic interaction between carboxyl group of Glu318 and the incipient imidazolium cation of His257, rather than the charge-relay mechanism suggested by experiment [1]. A much larger model (Model E) was created with 116 atoms. In addition to parts included in Model D, in Model E, larger parts of His257 including backbone ðCa and a amino group) of His257 was kept, and Asp288, which forms hydrogen bond to the a amino group of His257, was included and represented by an acetate molecule. The optimized structures of the reactant and transition state are also shown in Figure 5. As seen from Figure 5, the geometrical
parameters of R-E and TS-E are very similar to those of R-D and TS-D. Comparing the energies for the two models, the reaction energy of 7.4 kcal/mol with solvent effects for Model E is almost identical with the value of 7.2 kcal/mol obtained in Model D, and the barrier for Model E, however, changes substantially, about 4.7 kcal/mol (in the protein environment) lower than that of Model D. It is found that the calculated value of 19.1 kcal/mol with solvation effects shows good agreement with the experimental one (19.8 kcal/mol). The calculations demonstrate that Asp288 is also very important for lowering the barrier and necessary to be included in the model. Although as more groups (especially Glu318 and Asp288) were included in the models, the reaction barriers with solvation are in reasonable agreement with the experimental value, we found that the relative solvent effects of Models D and E are slightly larger. Inclusion of solvation correction using ¼ 4, the barriers are raised by 6.2 kcal/mol for Model D and 6.9 kcal/mol for Model E. And the reaction energies are increased dramatically. As compared to quite small solvation effects in Model C, the larger relative solvation effects in Models D and E are likely to result from the adding of negatively charged Glu318 located at the edge of both models. To improve the solvation effects, we therefore created a considerably larger model (Model F), which consists of 148 atoms. In addition to parts used in Model E, Model F includes parts of side chains of two residues Met289 and Met307 around Glu318 so as to reduce the sensitivity of energies with solvation effects by effectively shielding the negatively charged Glu318 close to the edge of model; meanwhile, larger parts of Glu318 and SAM were kept to make the two groups more flexible and to keep the locked centers further away from the reaction moieties. The optimized structures of the reactant, transition state and product are shown in Figure 6. As
Figure 6. Optimized structures of the reactant, transition state and product for Model F. Some hydrogen atoms are omitted for clarity.
506
F.-C. Cui et al. / Chemical Physics Letters 501 (2011) 502–507
seen from Figure 6, the structural parameters in general do not differ much from those obtained for Model E. Compared to TS-E, the breaking SD —C E bond distance in TS-F is slightly shorter (2.18 vs. 2.26 Å) and the forming C E —OD bond distance becomes somewhat longer (2.21 vs. 2.10 Å). These changes of the bond lengths indicate that the transition state for Model F lies slightly earlier than that for Model E. The barrier for Model F is calculated to be 14.6 kcal/ mol without solvation, while it increases somewhat by 2.4 kcal/ mol to be 17.0 kcal/mol on inclusion of solvation effects, which is slightly lower than the experimental value (19.8 kcal/mol). This underestimation might be caused by some possible sources of error, such as the computational methods (B3LYP, which generally underestimates a few kilocalories per mole on the barrier [32], and the employed basis sets), the size of the active-site model, the positions of truncation, and the solvation model used, etc. On the other hand, it should be pointed out that, although the solvation effects on the barrier are improved using this larger-sized model, the difference on the reaction energies with and without solvation is still large, by up to 12.8 kcal/mol. The likely reason for this large solvation correction on the reaction energy is due to the charge-state changes of several groups from the reactant to the product. In the cases of Models D–F, the presence of Glu318 generates two isolated charged groups, i.e., the phenolated anion of the substrate and protonated His257 at the reactant, while the charge states of SAM and substrate change from ion to neutral at the product, particularly, the positive charge of SAM located at the edge of the models disappearing. This results in different solvation energies between the reactant and the product. Similar case can be found for the reaction of the methyl transfer in guanidinoacetate methyltransferase [21]. Despite this, Model F makes considerable improvement on the solvation effects of the barrier, and reproduces the experimental value well. Therefore, we believe that Model F with 148 atoms is appropriate to investigate the catalytic mechanism and could provide a proper description for the present system. Also, systematic build-up of models from small to larger could provide qualitatively identification of some important groups to the catalytic reaction for IOMT. In addition, it is known that the pK a value reflects the acidity or basicity of the species, i.e., the larger the pK a value, the higher the ability for the species combining with the proton. Since the experimental pK a values of Histone and 7-hydroxyl of daidzein are 6.0 and 7.6 [33], respectively, the proton-transfer from daidzein to His257 is unlikely to occur spontaneously. These are the cases in Models A–C. However, as Glu318 is added, a different case is found that the proton of the substrate moves to His257 at reactant with no barrier. In order to give a reasonable explanation for this process, we employ a simple model, where His257, Glu318 and three water molecules are included (as shown in Figure 7), to roughly estimate the pK a value for N e of His257 in the presence of
Glu318 by using the empirical correction equation proposed by Qian et al. [30]. For comparison, the pK a values of N e of His257 in the absence of Glu318 and the 7-hydroxyl of truncated substrate daidzein are also calculated by the same method, and the calculated values are 5.9 and 7.0, respectively, very close to the corresponding experimental pK a data. It is found that the estimated pK a value of His257 in the presence of Glu318 increases to 7.73, which makes His257 become a stronger base than the 7-hydroxyl group of substrate daidzein. Consequently, the proton transfer from the substrate daidzein to His257 take place spontaneously and barrierlessly in the process of the formation of enzyme–substrate complex, to generate the phenolate anion. This facilitates the abstraction of the methyl cation from SAM via an SN 2 reaction, thereby lowering the activation barrier significantly in the models with the presence of Glu318 compared to those without Glu318. 4. Conclusions In this work, the catalytic reaction of methyl transfer from SAM to daidzein by the IOMT has been investigated by applying the DFT approach, and the catalytic roles of residues in the active site have also been evaluated based on several quantum chemistry models. The calculations indicate that the 7-hydroxyl proton of the substrate daidzein should be transferred to the N e atom of His257 spontaneously with no barrier when the substrate daidzein enters into the active site, and then the phenolate ion (strong nucleophile) abstracts the methyl from the positively charged SD of SAM in a single SN 2 reaction. On inclusion of solvation effects, the activation barrier is calculated to be 17.0 kcal/mol with Model F, in good agreement with the experimental value, and the reaction was exothermic by 12.8 kcal/mol. The calculated results suggest that Glu318 and Asp288 play an important role in decreasing the reaction barrier, while residues Asn258, Asn310, Met168 and Met311 may be responsible for constraining the spatial orientation of the reactant. The catalytic role of Glu318 proceeds likely through the electrostatic interaction. Acknowledgment This work was supported by the National Natural Science Foundation of China (20973077), the Program for New Century Excellent Talents in University (NCET). References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11]
Figure 7. Optimized structures of the general acid for pK a calculation: (A) protonated His257; (B) truncated daidzein; (C) protonated His257 with inclusion of Glu318 and three waters.
[12] [13] [14] [15] [16] [17] [18] [19] [20] [21] [22]
C. Zubieta, X.Z. He, R.A. Dixon, J.P. Noel, Nat. Struct. Biol. 8 (2001) 271. X.Z. He, R.A. Dixon, Arch. Biochem. Biophys. 336 (1996) 121. J.L. Ingham, in: Phytoalexins, Wiley, New York, 1982, p. 21. P.M. Dewick, M. Martin, Phytochemistry 18 (1979) 597. J. Vidgren, L. Svensson, A. Liljas, Nature 368 (1994) 354. R. Pattanayek, M.E. Newcomer, C. Wagner, Protein Sci. 7 (1998) 1326. W. Gong, M. O’Gara, R.M. Blumenthal, X. Cheng, Nucl. Acids Res. 25 (1997) 2702. X. Zhang, L. Zhou, X. Cheng, EMBO J. 19 (2000) 3509. J. Komoto, T. Yamada, Y. Takata, K. Konishi, H. Ogawa, T. Gomi, M. Fujioka, F. Takusagawa, Biochemistry 43 (2004) 14385. C.L. Gee, J.D.A. Tyndall, G.L. Grunewald, Q. Wu, M.J. McLeish, J.L. Martin, Biochemistry 44 (2005) 16875. C. Walsh, Enzymatic Reaction Mechanisms, W.H. Freeman and Company, San Francisco, 1979. H.L. Schubert, R.M. Blumenthal, X. Cheng, Trends Biochem. Sci. 28 (2003) 329. D.M. Blow, J.J. Birktoft, B.S. Hartley, Nature 221 (1969) 337. W.W. Bachovchin, J.D. Roberts, J. Am. Chem. Soc. 100 (1978) 8041. A.D. Becke, Phys. Rev. A 38 (1988) 3098. A.D. Becke, J. Chem. Phys. 96 (1992) 2155. A.D. Becke, J. Chem. Phys. 97 (1992) 9173. A.D. Becke, J. Chem. Phys. 98 (1993) 5648. C. Lee, W. Yang, R.G. Parr, Phys. Rev. B 37 (1988) 785. B. Miehlich, A. Savin, H. Stoll, H. Preuss, Chem. Phys. Lett. 157 (1989) 200. P. Velichkova, F. Himo, J. Phys. Chem. B 110 (2006) 16. P. Velichkova, F. Himo, J. Phys. Chem. B 109 (2005) 8216.
F.-C. Cui et al. / Chemical Physics Letters 501 (2011) 502–507 [23] P. Georgieva, Q. Wu, M.J. McLeish, F. Himo, Biochim. Biophys. Acta 1794 (2009) 831. [24] P. Georgieva, F. Himo, Chem. Phys. Lett. 463 (2008) 214. [25] M. Cossi, N. Rega, G. Scalmani, V. Barone, J. Chem. Phys. 114 (2001) 5691. [26] M. Cossi, G. Scalmani, N. Rega, V. Barone, J. Chem. Phys. 117 (2002) 43. [27] R. Cammi, B. Mennucci, J. Tomasi, J. Phys. Chem. A 103 (1999) 9100.
[28] [29] [30] [31] [32] [33]
507
R. Cammi, B. Mennucci, J. Tomasi, J. Phys. Chem. A 104 (2000) 5631. L. Noodleman, T. Lovell, W.G. Han, J. Li, F. Himo, Chem. Rev. 10 (2004) 459. H.T. Dong, H.B. Du, X.H. Qian, J. Phys. Chem. A 112 (2008) 12687. M. Cossi, G. Scalmani, N. Rega, V. Barone, J. Comput. Chem. 24 (2003) 669. Y. Zhao, N. Gonzàlez-Garcı`a, D.G. Truhlar, J. Phys. Chem. A 109 (2005) 2012. G.Y. Gao, D.J. Li, W.M. Keung, J. Med. Chem. 44 (2001) 3320.
Contents lists available at ScienceDirect
Chemical Physics Letters journal homepage: www.elsevier.com/locate/cplett
Reaction mechanism of isoflavone O-methyltransferase: A theoretical investigation Feng-Chao Cui, Xiao-Liang Pan, Jing-Yao Liu ⇑ State Key Laboratory of Theoretical and Computational Chemistry, Institute of Theoretical Chemistry, Jilin University, Changchun 130023, China
a r t i c l e
i n f o
Article history: Received 10 June 2010 In final form 14 November 2010 Available online 17 November 2010
a b s t r a c t The methyl-transfer reaction mechanism catalyzed by Isoflavone O-methyl-transferase (IOMT) and the roles of several residues around the active site are investigated by employing density functional method. The calculations confirm that the proton transfer from daidzein to His257 occurs barrierlessly, and the methyl group is transferred from S-adenosyl-L-methionine (SAM) to phenolate ion in a single SN 2 step with a barrier of 17.0 kcal/mol, consistent with experimental value. Glu318 and Asp288 play important roles in lowering the reaction barrier. Ó 2010 Published by Elsevier B.V.
1. Introduction Isoflavone O-methyltransferase (IOMT) belongs to S-adenosyl-Lmethionine (SAM) dependent plant natural product O-methyltransferase (OMT) involved in secondary metabolism [1] and is a key enzyme for the biosynthesis of the antimicrobial pterocarpan phytoalexin medicarpin, which is regarded as an important component of disease resistance response in alfalfa (Medicago sativa L.) [2,3]. As shown in Figure 1, in vivo, IOMT catalyzes the transfer of a methyl group from SAM to the B-ring 40 -hydroxyl group of the isoflavone daidzein, leading to the formation of S-adenosylhomocysteine (SAH) and formononetin, and the latter is an essential intermediate for synthesizing medicarpin in the isoflavonoid branch of the phenylpropanoid pathway [4]. However, in vitro, the A-ring 7-hydroxyl group of the isoflavone daidzein is preferentially methylated, yielding the isoformononetin, which is rarely found in plants and has no known physiological function [1,2]. The crystal structure of IOMT in complex with the products SAH and isoformononetin (40 -hydroxy-7-methoxyisoflavone) was first determined at 1.4 Å in 2001 by Noel et al. [1]. Figure 2 shows the major feature of the active site of this structure. Asn258, Asn310, Met168 and Met311 help in orienting the daidzein. Glu318 and Asp288 bracket the catalytic His257 residue. On the basis of this crystal structure, sequence alignments with the large family of plant OMTs and the mutational data, Noel et al. [1] put forward that similar to other methyltransferases [5–10], the putative catalytic base His257 is responsible for abstracting the proton from the 7-hydroxyl group of the A-ring of daidzein, and then the strong nucleophile (phenolate ion) generated on the deprotonated 7-hydroxyl of daidzein is now able to abstract the reactive methyl group from the positively charged SD of SAM in an SN 2 reaction [11,12], yielding the isoformononetin and SAH products [1]. They
⇑ Corresponding author. E-mail address: [email protected] (J.-Y. Liu). 0009-2614/$ - see front matter Ó 2010 Published by Elsevier B.V. doi:10.1016/j.cplett.2010.11.045
also speculated that the catalytic role of Glu318 is likely to proceed through the His–Glu charge-relay mechanism [13,14], that is, involving the proton transfer between Glu318 and His257, similar to protein arginine N-methyltransferase [8]. The purposes of this paper are to investigate the suggested mechanism and to examine the energetic feasibility of the methyl transfer reaction by employing the hybrid density functional theory (DFT) method B3LYP [15–20], which has been successfully used to study this kind of reaction [21–24]. The catalytic roles of several residues in the active site are also discussed by using several quantum chemical models with different sizes.
2. Computational details In the present study, quantum chemistry calculations were performed using the density functional theory (DFT) with Becke’s three-parameter (B3) [15–18] exchange functional along with the Lee–Yang–Parr (LYP) [19,20] nonlocal correlation functional (B3LYP). The geometries were optimized by using the 6-31 G(d,p) basis set. Hessians were calculated at the same level to confirm the nature of the stationary points, that is, each equilibrium species possesses all real frequencies, whereas the transition state possesses one and only one imaginary frequency. For each transition state, intrinsic reaction coordinate (IRC) calculations were performed to guarantee its correct connection to the designated local minima. To evaluate solvation effects of the missing protein environment, single-point calculations on the B3LYP/6-31 G (d,p) optimized geometries were performed using the conductor-like solvation model COSMO [25–28] at the same level. In the model, a cavity around the system is surrounded by polarizable dielectric continuum with the dielectric constant of 4, which is usually chosen for protein surrounding [29]. On the basis of these optimized geometries, single-point energies for all the stationary points were calculated at the B3LYP/6-311+G(2 d,2p) level to obtain more accurate energies. Some centers were fixed to preserve the positions of
F.-C. Cui et al. / Chemical Physics Letters 501 (2011) 502–507
HO 7
where pK cal a ðCPCMÞ value was estimated by using the thermodynamics cycles mentioned above. The above empirical correction formula was obtained by linear correlation relationship, and it has been tested that the calculated pK a values using Eq. (1) are very close to the corresponding experimental data. More computational details have been discussed in the literature [30]. In the present study, we performed the same-level calculations as those done in Ref.30 to estimate the pK cal a value. The gas-phase free energies of the acids and their conjugate bases were obtained by geometry optimizations and harmonic vibrational frequency calculations at the B3LYP/6-31+G* level of theory, and the solvation free energy of each species was calculated as a single point using conductor-like polarizable continuum model (CPCM) [31] at the HF/6-31+G* level based on HF/6-31+G*-optimized gas-phase structures. All computations were performed by using the GAUSSIAN03 package of programs.
O
A
daidzein O
B
4'
OH
In vivo HO
In vitro O
O CH3
O
O
formononetin
HO
CH3
O
O
isoformononetin
503
OH
O
3. Results and discussion O
medicarpin
O
CH3
Figure 1. O-methylation of the isoflavone daidzein catalyzed by Isoflavone Omethyltransferase (IOMT) in vivo and in vitro.
Figure 2. X-ray crystal structure of the active site of IOMT complex with SAH and isoformononetin.
the active site residues as similar as those in the crystal structure during the optimization. This constraint in geometry optimizations gives rise to a few imaginary frequencies, typically on the order of 10i—20i cm1 . These frequencies do not contribute significantly to the zero-point energies and can thus be ignored. In addition, it is found that the proton-transfer from daidzein to His257 occurs barrierlessly during the optimization process of reactant when Glu318 was included in the models. To provide the likely explanation for this barrierless proton-transfer process, the pK a values of the 7-hydroxyl of daidzein and the His257 N e in the presence or absence of Glu318 were estimated by employing thermodynamics cycles, which have been applied comprehensively for pK a calculations. This approach indirectly determines free-energy change for proton dissociation in aqueous via calculating free-energy change for proton dissociation in gas phase and the solvation free energies of the general acid and dissociated species. Here, we employed the following empirical correction equation proposed by Qian et al. [30] to evaluate pK a values:
pK a ¼ 0:418pK cal a ðCPCMÞ þ 2:27
ð1Þ
In the present study, several quantum chemical models of the IOMT active site with different sizes, ranging from 52 to 148 atoms, have been constructed based on the X-ray crystal structure of IOMT in complex with the products S-adenosyl-L-homocysteine (SAH) and isoformononetin (40 -hydroxy-7-methoxyisoflavone) (PDB code 1FP2) [1]. The substrate daidzein was built by displacing the methyl of 7-methoxy with hydrogen and SAH was converted to SAM by adding a methyl group at the S atom position. All hydrogen atoms were added manually. To keep the various groups in place to resemble the crystal structure as much as possible, the truncation points of the SAM cofactor and various active site residues were kept frozen to their X-ray positions during the geometry optimization. These fixed atoms are indicated by asterisks in the figures given later. In order to shed more light on the detailed methyl transfer process catalyzed by IOMT and evaluate the effects of the various residues at the active site of IOMT, we performed the calculations from smaller to larger models. The simplest model of the IOMT active site (called Model A) consists of 52 atoms and just contains the substrate daidzein, SAM cofactor and His257, which are essential for the catalytic reaction. The B-ring of the substrate daidzein was truncated and represented by a methyl. SAM was modeled by truncating the cofactor two carbons away on each side of the sulfur center. It has been substantiated by Himo et al. that this model of SAM is sufficient to reproduce the properties of the SD —C E bond [21–23]. His257 regarded as a base was represented by a methyl imidazole ring. The optimized structures of the reactant and transition state of Model are shown in Figure 3. The calculated potential energy profiles using all models with different sizes in the present investigation are displayed in Figure 4. As seen from TS-A in Figure 3, the distances of the breaking SD —C E bond and the forming C E —OD bond are 2.43 and 1.96 Å, respectively, and the transferred proton slightly moves away from the OD atom of the substrate daidzein toward the N e atom of His257, with OD —H and H—N e bond distances of 1.09 and 1.46 Å, respectively. This shows that the proton transfer from the substrate daidzein to His257 and the methyl cation transfer from SAM to the substrate daidzein occur in one concerted asynchronous step via an SN 2 reaction mechanism. The barrier of Model A is calculated to be 29.6 kcal/mol including solvation effects with ¼ 4, and the reaction is slightly endothermic by 1.5 kcal/mol. Clearly, the barrier of the simplest model with solvation effects is apparently higher than the experimental value of ca.19.8 kcal/mol estimated from the rate constantðkcat ¼ 0:03 s Þ [2] by using classical transition state theory. This poor agreement between theory and experiment indicates that some important residues for lowering the barrier of the catalytic reaction are missing and need to be included in the quantum chemical model.
504
F.-C. Cui et al. / Chemical Physics Letters 501 (2011) 502–507
Figure 3. Optimized structures of the reactant for Model A and transition states for Models A–C. Distances are in angstroms. Asterisks indicate that atoms are fixed to their Xray positions.
Relative energy (kcal/mol)
TS-C TS-A TS-B
30.1(31.0) 29.6(21.8) 29.5(31.0)
TS-D 23.8(17.6) TS-E 19.1(12.2) TS-F 17.0(14.6)
0.0
Pro-A
R-(A-F)
1.5(-2.8)
Pro-C -2.6(-2.9) Pro-B -4.1(-2.7) Pro-D -7.2(-18.7) Pro-E -7.4(-22.8) Pro-F -12.8(-25.6)
Reaction coordinate Figure 4. Computed potential energy profiles of all models in the solution and in gas phase (values in parentheses).
Two somewhat larger models (called Models B and C) of the IOMT active site were constructed. In Model B, the amide side chains of Asn258 and Asn310, which form hydrogen bonds to the substrate daidzein, were added based on Model A. In addition to parts included in Model B, Model C contains parts of side chains of Met168 and Met311, which sandwich the A ring of the substrate daidzein. The optimized TS structures for Models B and C (TS-B and TS-C) are also displayed in Figure 3. At both TSs, the critical SD —C E and C E —OD bond distances for the methyl transfer change a little compared to those of TS-A, while the OD —H and H—N e bond distances for proton transfer have larger variation, for example, H—N e bond distances are 1.08 and 1.07 Å in TS-B and TS-C (as shown in Figure 3), respectively. These results indicate that both models exhibit a single-step mechanism, which involves concerted the methyl-transfer and proton-transfer processes, similar to Model A, however, we find that the proton transfer takes place much earlier in these two models (almost finished at the two TSs) compared to that in Model A. On inclusion of solvent effects, the
calculated barriers of Models B and C are 29.5 and 30.1 kcal/mol, respectively, very close to that of Model A. Therefore, it appears that these four residues, which form hydrogen bonds or steric constraint to the substrate daidzein, may be helpful in preventing larger spatial reorganization of the substrate daidzein and the active site residues so as to keep them close to the crystal structure, while they don’t have effect on reducing the reaction barrier. Based on Model C, we constructed a larger model (called Model D), consisting of 104 atoms. This model contains a part of side chain of Glu318 (modeled by an acetate), which makes hydrogen bonding interaction to the N d of His257, and three water molecules around Glu318 presented in the crystal structure. The optimized structures of the reactant and transition state of Model D are shown in Figure 5. It is interesting to find that the proton of the substrate daidzein becomes unstable with inclusion of Glu318, and the optimization of the reactant always leads to the proton transfer to His257 (R-D in Figure 5). That is to say, in Model D, the barrierless proton-transfer from the substrate daidzein to
F.-C. Cui et al. / Chemical Physics Letters 501 (2011) 502–507
505
Figure 5. Optimized the structures of the reactants and transition-states for Models D and E.
His257 occurs firstly while the enzyme-substrate complex was formed, thereby generating two charged groups at the active site, i.e., the phenolate anion of the substrate daidzein and protonated His257, which is different from the behavior happening in the former three models. The likely reason of this phenomenon will be discussed in the later section of the text. Then, the methyl as a cation is transferred to the substrate daidzein in an SN 2 reaction, with the critical SD —C E and C E —OD bond distances being 2.29 and 2.13 Å, respectively. It is seen from Figure 4 that adding the residue Glu318 in the model results in a dramatic decrease of the barrier, about 6.3 kcal/mol lower than that of Model C, and an increase of the exothermicity by 4.6 kcal/mol, to 7.2 kcal/mol. It is thus conceivable to consider that Glu318 plays an important role in lowering the energetics of the catalytic reaction. However, the calculated barrier of Model D is still 4.0 kcal/mol higher than the experimental value (23.8 vs. 19.8 kcal/mol), implying more groups need to be added successively. In addition, because the hydrogen bonding interactions between the crystal water molecules and Glu318 may make the anionic Glu318 be stabilized, no proton-transfer from His257 to Glu318 occurs during the process of the catalytic reaction, implying that the catalytic role of Glu318 would proceed through the electrostatic interaction between carboxyl group of Glu318 and the incipient imidazolium cation of His257, rather than the charge-relay mechanism suggested by experiment [1]. A much larger model (Model E) was created with 116 atoms. In addition to parts included in Model D, in Model E, larger parts of His257 including backbone ðCa and a amino group) of His257 was kept, and Asp288, which forms hydrogen bond to the a amino group of His257, was included and represented by an acetate molecule. The optimized structures of the reactant and transition state are also shown in Figure 5. As seen from Figure 5, the geometrical
parameters of R-E and TS-E are very similar to those of R-D and TS-D. Comparing the energies for the two models, the reaction energy of 7.4 kcal/mol with solvent effects for Model E is almost identical with the value of 7.2 kcal/mol obtained in Model D, and the barrier for Model E, however, changes substantially, about 4.7 kcal/mol (in the protein environment) lower than that of Model D. It is found that the calculated value of 19.1 kcal/mol with solvation effects shows good agreement with the experimental one (19.8 kcal/mol). The calculations demonstrate that Asp288 is also very important for lowering the barrier and necessary to be included in the model. Although as more groups (especially Glu318 and Asp288) were included in the models, the reaction barriers with solvation are in reasonable agreement with the experimental value, we found that the relative solvent effects of Models D and E are slightly larger. Inclusion of solvation correction using ¼ 4, the barriers are raised by 6.2 kcal/mol for Model D and 6.9 kcal/mol for Model E. And the reaction energies are increased dramatically. As compared to quite small solvation effects in Model C, the larger relative solvation effects in Models D and E are likely to result from the adding of negatively charged Glu318 located at the edge of both models. To improve the solvation effects, we therefore created a considerably larger model (Model F), which consists of 148 atoms. In addition to parts used in Model E, Model F includes parts of side chains of two residues Met289 and Met307 around Glu318 so as to reduce the sensitivity of energies with solvation effects by effectively shielding the negatively charged Glu318 close to the edge of model; meanwhile, larger parts of Glu318 and SAM were kept to make the two groups more flexible and to keep the locked centers further away from the reaction moieties. The optimized structures of the reactant, transition state and product are shown in Figure 6. As
Figure 6. Optimized structures of the reactant, transition state and product for Model F. Some hydrogen atoms are omitted for clarity.
506
F.-C. Cui et al. / Chemical Physics Letters 501 (2011) 502–507
seen from Figure 6, the structural parameters in general do not differ much from those obtained for Model E. Compared to TS-E, the breaking SD —C E bond distance in TS-F is slightly shorter (2.18 vs. 2.26 Å) and the forming C E —OD bond distance becomes somewhat longer (2.21 vs. 2.10 Å). These changes of the bond lengths indicate that the transition state for Model F lies slightly earlier than that for Model E. The barrier for Model F is calculated to be 14.6 kcal/ mol without solvation, while it increases somewhat by 2.4 kcal/ mol to be 17.0 kcal/mol on inclusion of solvation effects, which is slightly lower than the experimental value (19.8 kcal/mol). This underestimation might be caused by some possible sources of error, such as the computational methods (B3LYP, which generally underestimates a few kilocalories per mole on the barrier [32], and the employed basis sets), the size of the active-site model, the positions of truncation, and the solvation model used, etc. On the other hand, it should be pointed out that, although the solvation effects on the barrier are improved using this larger-sized model, the difference on the reaction energies with and without solvation is still large, by up to 12.8 kcal/mol. The likely reason for this large solvation correction on the reaction energy is due to the charge-state changes of several groups from the reactant to the product. In the cases of Models D–F, the presence of Glu318 generates two isolated charged groups, i.e., the phenolated anion of the substrate and protonated His257 at the reactant, while the charge states of SAM and substrate change from ion to neutral at the product, particularly, the positive charge of SAM located at the edge of the models disappearing. This results in different solvation energies between the reactant and the product. Similar case can be found for the reaction of the methyl transfer in guanidinoacetate methyltransferase [21]. Despite this, Model F makes considerable improvement on the solvation effects of the barrier, and reproduces the experimental value well. Therefore, we believe that Model F with 148 atoms is appropriate to investigate the catalytic mechanism and could provide a proper description for the present system. Also, systematic build-up of models from small to larger could provide qualitatively identification of some important groups to the catalytic reaction for IOMT. In addition, it is known that the pK a value reflects the acidity or basicity of the species, i.e., the larger the pK a value, the higher the ability for the species combining with the proton. Since the experimental pK a values of Histone and 7-hydroxyl of daidzein are 6.0 and 7.6 [33], respectively, the proton-transfer from daidzein to His257 is unlikely to occur spontaneously. These are the cases in Models A–C. However, as Glu318 is added, a different case is found that the proton of the substrate moves to His257 at reactant with no barrier. In order to give a reasonable explanation for this process, we employ a simple model, where His257, Glu318 and three water molecules are included (as shown in Figure 7), to roughly estimate the pK a value for N e of His257 in the presence of
Glu318 by using the empirical correction equation proposed by Qian et al. [30]. For comparison, the pK a values of N e of His257 in the absence of Glu318 and the 7-hydroxyl of truncated substrate daidzein are also calculated by the same method, and the calculated values are 5.9 and 7.0, respectively, very close to the corresponding experimental pK a data. It is found that the estimated pK a value of His257 in the presence of Glu318 increases to 7.73, which makes His257 become a stronger base than the 7-hydroxyl group of substrate daidzein. Consequently, the proton transfer from the substrate daidzein to His257 take place spontaneously and barrierlessly in the process of the formation of enzyme–substrate complex, to generate the phenolate anion. This facilitates the abstraction of the methyl cation from SAM via an SN 2 reaction, thereby lowering the activation barrier significantly in the models with the presence of Glu318 compared to those without Glu318. 4. Conclusions In this work, the catalytic reaction of methyl transfer from SAM to daidzein by the IOMT has been investigated by applying the DFT approach, and the catalytic roles of residues in the active site have also been evaluated based on several quantum chemistry models. The calculations indicate that the 7-hydroxyl proton of the substrate daidzein should be transferred to the N e atom of His257 spontaneously with no barrier when the substrate daidzein enters into the active site, and then the phenolate ion (strong nucleophile) abstracts the methyl from the positively charged SD of SAM in a single SN 2 reaction. On inclusion of solvation effects, the activation barrier is calculated to be 17.0 kcal/mol with Model F, in good agreement with the experimental value, and the reaction was exothermic by 12.8 kcal/mol. The calculated results suggest that Glu318 and Asp288 play an important role in decreasing the reaction barrier, while residues Asn258, Asn310, Met168 and Met311 may be responsible for constraining the spatial orientation of the reactant. The catalytic role of Glu318 proceeds likely through the electrostatic interaction. Acknowledgment This work was supported by the National Natural Science Foundation of China (20973077), the Program for New Century Excellent Talents in University (NCET). References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11]
Figure 7. Optimized structures of the general acid for pK a calculation: (A) protonated His257; (B) truncated daidzein; (C) protonated His257 with inclusion of Glu318 and three waters.
[12] [13] [14] [15] [16] [17] [18] [19] [20] [21] [22]
C. Zubieta, X.Z. He, R.A. Dixon, J.P. Noel, Nat. Struct. Biol. 8 (2001) 271. X.Z. He, R.A. Dixon, Arch. Biochem. Biophys. 336 (1996) 121. J.L. Ingham, in: Phytoalexins, Wiley, New York, 1982, p. 21. P.M. Dewick, M. Martin, Phytochemistry 18 (1979) 597. J. Vidgren, L. Svensson, A. Liljas, Nature 368 (1994) 354. R. Pattanayek, M.E. Newcomer, C. Wagner, Protein Sci. 7 (1998) 1326. W. Gong, M. O’Gara, R.M. Blumenthal, X. Cheng, Nucl. Acids Res. 25 (1997) 2702. X. Zhang, L. Zhou, X. Cheng, EMBO J. 19 (2000) 3509. J. Komoto, T. Yamada, Y. Takata, K. Konishi, H. Ogawa, T. Gomi, M. Fujioka, F. Takusagawa, Biochemistry 43 (2004) 14385. C.L. Gee, J.D.A. Tyndall, G.L. Grunewald, Q. Wu, M.J. McLeish, J.L. Martin, Biochemistry 44 (2005) 16875. C. Walsh, Enzymatic Reaction Mechanisms, W.H. Freeman and Company, San Francisco, 1979. H.L. Schubert, R.M. Blumenthal, X. Cheng, Trends Biochem. Sci. 28 (2003) 329. D.M. Blow, J.J. Birktoft, B.S. Hartley, Nature 221 (1969) 337. W.W. Bachovchin, J.D. Roberts, J. Am. Chem. Soc. 100 (1978) 8041. A.D. Becke, Phys. Rev. A 38 (1988) 3098. A.D. Becke, J. Chem. Phys. 96 (1992) 2155. A.D. Becke, J. Chem. Phys. 97 (1992) 9173. A.D. Becke, J. Chem. Phys. 98 (1993) 5648. C. Lee, W. Yang, R.G. Parr, Phys. Rev. B 37 (1988) 785. B. Miehlich, A. Savin, H. Stoll, H. Preuss, Chem. Phys. Lett. 157 (1989) 200. P. Velichkova, F. Himo, J. Phys. Chem. B 110 (2006) 16. P. Velichkova, F. Himo, J. Phys. Chem. B 109 (2005) 8216.
F.-C. Cui et al. / Chemical Physics Letters 501 (2011) 502–507 [23] P. Georgieva, Q. Wu, M.J. McLeish, F. Himo, Biochim. Biophys. Acta 1794 (2009) 831. [24] P. Georgieva, F. Himo, Chem. Phys. Lett. 463 (2008) 214. [25] M. Cossi, N. Rega, G. Scalmani, V. Barone, J. Chem. Phys. 114 (2001) 5691. [26] M. Cossi, G. Scalmani, N. Rega, V. Barone, J. Chem. Phys. 117 (2002) 43. [27] R. Cammi, B. Mennucci, J. Tomasi, J. Phys. Chem. A 103 (1999) 9100.
[28] [29] [30] [31] [32] [33]
507
R. Cammi, B. Mennucci, J. Tomasi, J. Phys. Chem. A 104 (2000) 5631. L. Noodleman, T. Lovell, W.G. Han, J. Li, F. Himo, Chem. Rev. 10 (2004) 459. H.T. Dong, H.B. Du, X.H. Qian, J. Phys. Chem. A 112 (2008) 12687. M. Cossi, G. Scalmani, N. Rega, V. Barone, J. Comput. Chem. 24 (2003) 669. Y. Zhao, N. Gonzàlez-Garcı`a, D.G. Truhlar, J. Phys. Chem. A 109 (2005) 2012. G.Y. Gao, D.J. Li, W.M. Keung, J. Med. Chem. 44 (2001) 3320.