Molecular Modeling and Docking Studies of 3'-Hydroxy-N-methylcoclaurine 4'-O-Methyltransferase from Coptis chinensis

Abstract

Coptis chinensis 3'-hydroxy-N-methylcoclaurine 4'-O-methyltransferase (HOMT), an essential enzyme in the berberine biosynthetic pathway, catalyzes the methylation of 3'-hydroxy-N-methylcoclaurine (HMC) producing reticuline. A 3D model of HOMT was constructed by homology modeling and further subjected to docking with its ligands and molecular dynamics simulations. The 3D structure of HOMT revealed unique structural features which permitted the methylation of HMC. The methylation of HMC was proposed to proceed by deprotonation of the 4'-hydroxyl group via His257 and Asp258 of HOMT, followed by a nucleophilic attack on the SAM-methyl group resulting in reticuline. HOMT showed high substrate specificity for methylation of HMC. The study evidenced that Gly117, Thr312 and Asp258 in HOMT might be the key residues for orienting substrate for specific catalysis.

Introduction

Berberine, a benzylisoquinoline alkaloid obtained from Coptis species, is an important pharmaceutical alkaloid with anti-bacterial, anti-inflammatory, hypotensive, hypoglycemic and anti-tumor activities.1 Both medical importance of berberine and interests in elucidating the biosynthesis of alkaloids had led to the complete characterization of berberine biosynthetic pathway2 (Figure 1). 3'-hydroxy-N-methylcoclaurine 4'-O-methyltransferase (HOMT),3 belonging to the S-adenosyl-L-methionine (SAM)-dependent O-methyltransferase family,4,5 is an essential enzyme in the biosynthetic pathway of berberine.6,7 HOMT possesses a high specificity to catalyze the transfer of a methyl group from SAM to the 4'-oxygen atom of 3'-hydroxy-N-methylcoclaurine (HMC), resulting in the formation S-adenosyl-Lhomocysteine (SAH) and reticuline8 (Figure 1).

The production of berberine is difficult to prepare by organic synthesis methods and the expression of a HOMT gene with high catalytic activity in a heterologous host is a challenging task that will lead to efficient biotransformation. One goal in HOMT research is to understand the specificity of that enzyme in stereo selective reactions for the production of berberine. However, the prerequisite to understand the specificity recognition of an enzyme is the knowledge of its 3D structure. To date, many attempts to crystallize HOMT have been unsuccessful. This prompted us to study the 3D structure of HOMT by homology modeling. In the present study, the homology model of HOMT had been created for deeper insights into its structure and function. After validating the consistency of the model, it was further subjected to docking with its ligands and molecular dynamics simulations to study the substrate recognition and catalytic mechanism of HOMT.

Figure 1.The methylation of 3'-hydroxy-N-methylcoclaurine catalyzed by HOMT in the biosynthetic pathway of berberine.8

 

Materials and Methods

Homology Modeling. The protein sequence of HOMT, which contained 350 amino acid residues, was obtained from NCBI database (Accession: ABY75613). For homology modeling of HOMT, the template Medicago truncatula isoflavanone O-methyltransferase (IOMT)9 (PDB: 2QYO_A) was identified from PDB10 by BLAST.11

MODELER12 module included in Discovery Studio (DS) 2.013,14 was used for homology modeling of HOMT. The sequence of the template protein (IOMT) was extracted and aligned with the target protein (HOMT) by Align Multiple Sequences module of DS. Based on the sequence alignment, the structures of conserved regions (SCRs) were determined. The structures of loops, variable regions and side chains not included in SCRs were built automatically by MODELER.

Molecular Dynamics Simulation. The initial model from MODELER was optimized by energy minimization (EM)14,15 and subsequent molecular dynamics (MD)16,17 simulation in the force field of CHARMm18 provided in DS. The initial model was solvated with the TIP3P model of water using the orthorhombic box with a minimum distance of at least 10 Å from any edge of the box to any protein atom. The solvated protein system was subsequently subjected to a thorough EM before MD simulation. Firstly, the water molecules of solvated protein system was energy-minimized by 2000 steps of steepest descent algorithm19 while holding the solute frozen. Then, 5000 steps of conjugate gradient energy-minimization20 of the whole system was carried out to remove close contacts and relax the system.

The system was then subjected to 2000 ps of heating from 50 to 300 K with a Harmonic constraint14 imposed on the solute to allow for relaxation of water molecules. The following 2000 ps equilibrium run for the system was performed with distance Harmonic constraints. The production run was next carried out without any constraints in the NPT ensemble14 with periodic boundary conditions for 10 ns at 1 atm and 300 K. During the simulations, all covalent bonds involving hydrogen were constrained using the SHAKE algorithm21 so that a time step was set to 2 fs. Long-range electrostatic interactions were treated with the particle mesh Ewald (PME) method22. Finally, the output conformers were collected at 10 ps interval and the optimized model was saved for validation.

Validation of Model. The stereo chemical quality of HOMT structure was checked by Ramachandran plot using PROCHECK.23 The PROSA24 test was applied for the final model to check for energy criteria in comparison with the potential of mean force derived from a large set of known protein structures. Further, the compatibility of the model with its sequence was measured by Verify-3D graph13,25 using Profile-3D module of DS. Finally, the root mean square deviation (RMSD) between the backbone atoms of HOMT and that of the template (IOMT) was calculated by Proteins Superimposing module of DS.

Docking of Ligands into the Protein. As mentioned above, the methylation of HMC catalyzed by HOMT needs a co-factor, SAM. Therefore, the ligands for docking consisted of HMC and SAM. The structures of HMC (PubChem ID: 25203738) and SAM (PubChem ID: 34756) were obtained from NCBI PubChem.26 The binding sites between the protein and the ligands were identified by Binding Site tool of DS. After determining the binding sites, SAM and HMC was docked into the protein using CDOCKER27 module of DS. The final docking model with the minimum binding energy and reasonable spatial conformation was chosen and then subjected to EM and MD simulation using the same method as that used for the initial model of HOMT.

In order to study the substrate specificity of HOMT, the compounds norcoclaurine (PubChem ID: 440927), coclaurine (PubChem ID: 46878403), N-methylcoclaurine (PubChem ID: 440591), reticuline (PubChem ID: 439653), scoulerine (PubChem ID: 439654) and tetrahydrocolumbamine (PubChem ID: 25203738) obtained from PubChem, structures of which were similar to HMC, were also docked with HOMT and further subjected to EM and MD simulation using the same method as that used for HMC.

 

Results

Homology Modeling of HOMT. Amino acid sequence alignment showed that HOMT shared 41% sequence identity and 63% sequence similarity with IOMT, suggesting that they may share a similar overall structure. Both HOMT and IOMT were SAM-dependent O-methyltransferases5,28 and the reactions they catalyzed were the methylation of similar compounds mediated by SAM as the co-factor. A number of potential active residues in IOMT were conserved in HOMT (Figure 2), including those interacting with SAM and substrate, as well as those involved in the methyl transfer system. Based on the above considerations, it was reasonable to choose IOMT as a template for constructing a 3D model of HOMT. IOMT consisted of 353 residues (5-357) visible in the crystal structure and the first four residues were not visible. The 3D structure of HOMT, which consisted of 347 residues (4-350) omitting the first three residues, was constructed based on the crystal structure of IOMT.

Figure 2.Sequence alignment of HOMT with IOMT by using Proteins Aligning module of DS. The secondary structures were demonstrated with “α” for α-helix, “β” for β-sheet and “c” for loop. The residues for substrate binding and SAM binding were shown in blue and pink boxes, respectively. The catalytic residues were shown in red boxes. These active residues were potential and identified from IOMT9 and other representative plant OMTs.28-30

Refinement of HOMT Model by MD Simulation. The constructed model was subjected to MD simulation, in order to assess the stability of the model and to find the energetically favorable structure for further docking study. The MD simulation trajectory-based analysis showed that the potential energy of the model gradually decreased from −241196.113 kcal/mol to −245847.711 kcal/mol and trended to balance. Structural stability of the constructed model during the MD simulation was examined using RMSD. As shown in Figure 3, the RMSD plot for the protein backbone atoms with reference to the initial structure, as a function of time. The first 2 ns were considered as an equilibration period followed by a plateau where the model became stable in the RMSD range about 2.2-2.7 Å.

Distinguishing the flexible regions of the protein could help in understanding the stability of the protein. To examine the flexible regions of the model, we had generated the average root mean square fluctuation (RMSF) plot for the Cα-atoms with respect to the residues. Those residues, which deviated more than 2.5 Å were considered as highly flexible elements of the protein.31 As shown in Figure 4, the high fluctuations occurred in the loop regions and the low fluctuations occurred in the α-helix and β-sheet regions. The N-terminal residues 4-8, loop-helix (23-26), loop-helix (63-65), loop (87-90), loop-helix (126-130), loop (145-147), helix (172-174), helix (180-183), loop (196-198), and loop (277-279) showed high flexibility. The regions with RMSF values lower than 2.0 clearly showed that the protein structural core is well-constructed and the model is stable in the course of MD simulation.

Figure 3.The RMSD plot calculated for backbone atoms of HOMT and ligands-bound complex, and all atoms of ligands during MD simulations.

Figure 4.The RMSF plot calculated for Cα-atoms of residues in HOMT and ligands-bound complex.

Validation of HOMT Model. The stereo chemical quality of the HOMT model was checked by PROCHECK based on the Ramachandran plot. The results, summarized in Figure S1, showed 99.7% of the residues in the favored or allowed regions and only 0.3% in disallowed region. The Goodness factors (G-factors), indicating the quality of covalent and bond/angle distance, were 0.05 for dihedrals, 0.19 for covalent, and 0.11 for overall. Then, the ERRAT score was calculated based on the non-bonded atomic interactions.14,32 In this case, the ERRAT score for the model was 94.69, which was within the range of a high quality model.

In order to investigate whether the interaction energy of each residue of the protein was negative, a second test was done using the energy criteria of PROSA plot. The PROSA analysis of the model showed that maximum residues of HOMT had negative interaction energy and very few residues displayed positive interaction energy (Figure S2).

Subsequently, the environmental profile of each residue was assessed by the Profile-3D program. The compatibility score above zero in the Verify-3D graph (Figure S3) corresponds to acceptability side-chain environments.33 About 90% of the residues had an average score more than 0.2, which suggested that the side chain environments were acceptable. Hence, the predicted model was determined to be good enough for further docking studies.

Structural Features of HOMT. The constructed 3D structure of HOMT was predominantly organized in 16 long α-helices (45.53% α-helical residues) and 12 short β-sheets (16.43% β-sheeted residues) (Figure 2). The sequences and structures of HOMT and IOMT were conserved due to their similar biological functions. Like IOMT, the HOMT activesites were composed of two domains (Figure 5). The substrate binding domain (residues 114-170 and 304-313) was composed of six α-helices (α7-12, α15) and five loops (c10-13, c24) (Figure 2). The SAM binding domain (residues 194-259) was composed of four β-sheets (β4-7), two α- helices (α12, α23) and six loop (c15-20) (Figure 2). The superimposition of HOMT model with its template (IOMT) resulted in RMSD value of 1.747 Å, indicating that the 3D structure of HOMT model was similar to that of IOMT (Figure 5). The substrate binding domain and SAM binding domain contributed significantly to the overall similarity of the two backbones. The high similarities of their putative active domains and their secondary structures attested to the high precision of the model.

Figure 5.The superimposition of the structure of HOMT with that of IOMT using Proteins Superimposing module of DS. The pink ribbon represented IOMT and the green one represented HOMT. The substrate binding domain and SAM binding domain were identified from IOMT.9

Active Sites of HOMT. The binding sites of HOMT suggested by Binding Site tool of DS were similar to those of IOMT determined from the X-ray crystal structure. The active sites of SAM-dependent methyltransferases, especially the SAM binding site, are highly conserved.34,35 As shown in Figure 2, the SAM binding site was composed of conserved residues and the substrate binding site was lined by mostly hydrophobic residues. Comparative substrate binding site analysis of IOMT and HOMT showed various differences in the respective substrate binding residues. For example, residues Gly117, Met118, Phe123, Met124, Trp127, Ile149, Met308, Leu309, and Thr312 in IOMT were replaced by Cys122, Val123, Lys128, Ser129, Phe132, Leu155, Ile316, Ser317, Asn320 in HOMT, respectively. These relative differences in the active site should undoubtedly affect the conformation and presumably play significant roles in the substrate specificity. 36,37 The residues His257, Asp258, Glu286, and Glu316 were considered to be involved in catalysis and completely conserved in these two proteins.

Docking and MD Simulation of HOMT with Ligands. Based on the active sites, SAM and HMC were docked into the protein with interaction energies of −65.482 kcal/mol and −61.741 kcal/mol, respectively (Table 1). Even though SAM and HMC were docked into the protein, the 3D structures of HOMT were not changed much. The comparison of the conformations before and after docking only gave less RMSD value than 2 Å. To confirm the stability of the model, HOMT-substrate-SAM complex was subjected to MD simulations. RMSD of backbone atoms of HOMT-substrate-SAM complex and all atoms of HMC and SAM were calculated from their initial structure as a function of simulation time. As shown in Figure 3, it was indicated that the structure of HOMT-substrate-SAM complex showed relatively low RMSD value (2.539 Å) when compared with the structure of HOMT (2.947 Å), suggesting that the presence of substrate and SAM stabilized HOMT-substrate-SAM complex. The overall RMSD of SAM was 1.194 Å and that of HMC was 1.245 Å, indicating that the position and orientation of both ligands were mostly maintained14 in the binding sites during simulations.

Table 1.aEvdW: van der Waals interaction energy. bEele: electrostatic interaction energy. cEtotal: total interaction energy;

A comparison of structures of ligands-bound complex and HOMT indicated that structural changes were in loops and active site regions (Figure 4). These regions showed lower RMSF value in the ligands-bound complex than in the HOMT, suggesting that the presence of ligands decreased the fluctuations of these regions. This is not surprising, because ligands were not bound with these regions of HOMT before docking. Comparison of the RMSF of the ligandsbound complex with that of HOMT also showed that there were many features common to both curves. However, several apparent differences were noted. For example, the RMSF values of loop (223-225), loop-sheet (231-233), sheet-loop (285-287) and loop-sheet (315-317) were significantly higher in the ligands-bound complex than in HOMT, indicating that these regions were relatively stable in HOMT but could be flexible with the presence of ligands.

The active site residues in loop regions showed high RMSF values than in helix or sheet regions. As a result, residues in loop region made the active site flexible. As shown in Figure 6, the active site residues were relatively flexible in HOMT than in the ligands-bound complex, indicating that the presence of SAM and HMC stabilized the active site in the ligands-bound complex. For example, the RMSF values of SAM binding residues (Gly196, Asp219, Asp239, Met240, and Lys253), HMC binding residues (Gly117 and Thr312 and Asp258) and catalytic residues (His257 and Asp258) were much lower in the ligands-bound complex than in HOMT, mostly due to their H-bonds with SAM and HMC. It was speculated from the results of RMSF value of active site residues that these residues might be involved in the catalytic mechanisms.14 In addition, the RMSF values of Leu254, Asp286 and Glu316 were relatively higher in the ligandsbound complex than in HOMT, suggesting that these residues were relatively stable with the absence of SAM and HMC but flexible in ligands-bound complex, which could contribute to the catalysis roles of His 257 and Asp258.

Figure 6.The RMSF plot for the Cα-atoms of the active site residues in HOMT and ligands-bound complex.

Binding mode of Ligands in HOMT. Like other SAMdependent O-methyltransferases,5,28 SAM binding in HOMT was mediated by an extensive network of van der Waals and H-bonds interactions (Figure 7). The heterocyclic adenine ring of SAM was sequestered by van der Waals interactions with a number of hydrophobic residues, in which Leu220 and Met240 bracketed the ring, and Asp219 and Phe241 were involved in a favorable edge-to-face interaction.30 The ribose of SAM was bracketed by Trp259 and Asp219. The tail of SAM was also sequestered by a number of hydrophobic residues (Gly196, Gly197, Thr202, Lys253, and Leu254). Furthermore, Gly196, Lys253, Asp219, Asp239, and Met240 formed H-bonds with the terminal amino group, carboxylate group, ribose hydroxyl group, exocyclic amino group and cyclic amino group in the adenine ring of SAM, respectively.

Figure 7.Binding mode of ligands in HOMT. (a) Ribbon diagram of the HOMT in complex with SAM and HMC. (b) Surface representation of the active site cavity for the HOMT-SAM-HMC complex. (c) The interactions of binding residues in HOMT with HMC and SAM. “Me” indicated methyl group and the dotted lines indicated H-bonds. Only key residues were shown in the figure. (d) Hypothetical mechanism for the methylation of HMC catalyzed by HOMT. The arrow indicated the transfer of electron. The dotted line indicated the H-bond between Asp258 and HMC.

HMC was stabilized by H-bonds and van der Waals interactions (Figure 7). The hydrophilic residues Thr312, Gly313, Gly314, and hydrophobic residue Leu309 sequestered the benzene ring linked by 4'-oxygen atom, which presented the reactive 4'-oxygen atom to the catalytic residues (His257 and Asp258) and SAM. The H-bond of Thr312 with the 3'-oxygen atom and that of Asp258 with the 4'-oxygen atom helped orient HMC for preferential methylation of the 4'-oxygen atom (Figure 7). The distances of the 4'-oxygen atom of HMC from Asp258 and the methyl group of SAM were 3.084 and 5.424 Å, respectively. Meanwhile, the 3'-oxygen atom of HMC was 6.830 and 10.334 Å far from Asp258 and the methyl group of SAM, respectively. These results hence pointed out that the distance constraints in the binding site of HOMT were important for the selectivity of the catalytic group. The closer the 3'-oxygen and 4'-oxygen atom of HMC approach SAM and Asp258, the easier they would be utilized by the methyl group of SAM and ultimately lead to potent catalysis. In this case, it could be assumed that the 4'-oxygen atom might be more catalytically potent as compared with the 3'-oxygen atom of HMC.

The environment for binding the heterocyclic fused ring of HMC consisted of a number of hydrophobic residues (Met118, Phe123, Met124, Trp127, Phe141, Ile149 and Met308) and the hydrophilic residue Gly149. The structures and relative hydrophobicity binding pocket for the heterocyclic fused ring, as well as the H-bond of Gly117 with the 7-oxygen atom, were able to constraint the heterocyclic fused ring and the 7-oxygen atom, which suggested the selectivity for the benzene ring linked by 4'-oxygen atom instead of the heterocyclic fused ring by this region property of HOMT. The specific pattern of favorable van der Waals interactions and H-bonds of HOMT with HMC most likely explained the kinetic preferences of HOMT for the substrate.

Proposed Catalytic Mechanism. As mentioned above, the potential catalytic residues of HOMT were composed of His257, Asp258, Glu286, and Glu316. The imidazole ring of His257 could be protonated in the acid environment provided by Asp258, Glu286, and Glu316. Asp258 served as a suitable active site residue because of its low pKa value, the dominant effect and stability of its negative charge, and its favorable charge interactions with different ligands.30 Like the scoulerine 9-O-methyltransferase,38 we could proposed that the methylation of HMC involved deprotonation of its 4'-hydroxyl group by the neighboring His257 and Asp258 residues (Figure 7(d)). An electron was transferred from water molecules within the reaction centre to the imidazole ring of His257, and further to Asp258 giving rise to 4'-hydroxyl group becoming an electron-enriched group.39 Then, the SAM-methyl group was nucleophilic attacked by 4'-hydroxyl group of HMC resulting in the formation of reticuline and SAH.

Figure 8.The interactions of different ligands with the catalysis residues and SAM-methyl group. The dotted line indicated the H-bond and the blue line indicated the distance between atoms.

Substrate Specificity of HOMT. In order to study the substrate specificity of HOMT, the compounds norcoclaurine, coclaurine, N-methylcoclaurine, reticuline, scoulerine and tetrahydrocolumbamine, structures of which were similar to HMC, were also docked with HOMT. The results were shown in Table 1 and Figure 8. The potential methylation oxygens of norcoclaurine, coclaurine, scoulerine and tetrahydrocolumbamine were too far from Asp258 and SAMmethyl group. Their interaction energies with HOMT were higer than that of HMC. In addition, the potential methylation oxygens of norcoclaurine and coclaurine were restrained by the H-bonds with Asn311. Therefore they could not be methylated by HOMT. N-methylcoclaurine and reticuline showed the similar conformation as that of HMC, but their potential methylation oxygens were relatively far from Asp258 and SAM-methyl group and restrained by H-bonds with Thr312. Their interaction energies with HOMT were relatively higer than that of HMC. It indicated that N-methylcoclaurine and reticuline were not a suitable substrate for the methylation. Thus, it was evidenced that HOMT showed high substrate specificity for methylation of HMC. Considering the H-bond of Gly117 with the 7-oxygen atom, which was able to constrain the 7-oxygen atom, a conclusion could be drawn that Gly117, Thr312 and Asp258 in HOMT might be the key residues for orienting substrate for specific catalysis.

 

Conclusions

The 3D structure of HOMT, predominantly α-helical in nature with some β-sheeted residues, was similar to that of IOMT but also revealed unique structural features which permitted the methylation of HMC. The SAM binding residues of HOMT were mostly identical to those of other SAM-dependent O-methyltransferases, while the substrate binding residues were variable. These relative differences in the substrate binding site should undoubtedly affect the conformation and presumably play significant roles in the substrate specificity. The structural changes of active sites before and after binding with ligands could contribute to the catalysis of HOMT. The methylation of HMC was proposed to proceed by deprotonation of the 4'-hydroxyl group via His257 and Asp258 of HOMT, followed by a nucleophilic attack on the SAM-methyl group resulting in reticuline. The study evidenced that Gly117, Thr312 and Asp258 in HOMT might be the key residues for orienting substrate for specific catalysis.