Metadynamics Simulation Study on the Conformational Transformation of HhaI Methyltransferase: An Induced-Fit Base-Flipping Hypothesis

DNA methyltransferases play crucial roles in establishing and maintenance of DNA methylation, which is an important epigenetic mark. Flipping the target cytosine out of the DNA helical stack and into the active site of protein provides DNA methyltransferases with an opportunity to access and modify the genetic information hidden in DNA. To investigate the conversion process of base flipping in the HhaI methyltransferase (M.HhaI), we performed different molecular simulation approaches on M.HhaI-DNA-S-adenosylhomocysteine ternary complex. The results demonstrate that the nonspecific binding of DNA to M.HhaI is initially induced by electrostatic interactions. Differences in chemical environment between the major and minor grooves determine the orientation of DNA. Gln237 at the target recognition loop recognizes the GCGC base pair from the major groove side by hydrogen bonds. In addition, catalytic loop motion is a key factor during this process. Our study indicates that base flipping is likely to be an “induced-fit” process. This study provides a solid foundation for future studies on the discovery and development of mechanism-based DNA methyltransferases regulators.


Post processing protocol
Trajectory analysis is performed using different tools. We use CatDCD converted original DCD type trajectory into TRR and PDB format. RMSD of catalytic loop and Distance map was plotted by g_rms and g_dist tool in Gromacs 4.5.5 4 . Hydrogen bond existence map and Hydrogen bond number map was plotted using g_hbond tools in Gromacs. We employ Curves+ 5 to monitor the groove width. Elctrostatic surface is generated by using APBS 6 plug-ins in PyMol. The secondary structure profile of M.HhaI was generated using do_dssp in Gromacs 4.6. Then, we use xpm2ps converted xpm format into PostScript file.
Three dimensional structures were made by PyMol. Free energy profile was generated using the RGL module in the R program 7 . The other plots are drawn using the XmGrace (http://plasma-gate.weizmann.ac.il/Grace/) plotting tool. The VMD program 8 was used for trajectory visualization.

Preliminary Metadynamics Simulations
Choosing proper collective variables (CVs) is vital to metadynamics simulation. Because catalytic loop motion was a very complex process, we tried six different CVs and their combination: (1) The distance between the center of mass (COM) of the catalytic loop and two target recognition loops; (2) The distance between amino group in the side chain of Gln237 and carboxyl group of Ile86; (3) Different psi and phi angle combinations of residues from 81 to 89, in order to monitor the loop formation process; (4) Alpha helix formation between residue 82 and 84; (5) Contact map 9 between GCGC motif and Target recognition loop of M.HhaI.
We also use TMD trajectory to test these CVs and their combinations (in Table S1). Plotting the value of these CVs along TMD trajectory, we find they can separate initiation and ending successfully.
However, no obvious conformational change of catalytic loop can be found in following metadynamic simulation. All these runs are considered unsuccessful. Since the complexity of the simulation system, it requires a CV to describe the motion of catalytic loop. Thus, we tried DMSD as CV to define the transition path. We also tried different ways to describe this CV: using atoms in the backbone of catalytic loop or heavy atoms in the catalytic loop. Finally, the result shows that using all non-hydrogen atoms could successfully drive the catalytic loop moving from open to closed form only.

Flipped base in semi-closed form
Base flipping take place at about 230ns, and this progress is rapid. After target base flipped, hydrogen bonds between Arg163, Arg165, Phe79 and target cytosine is observed. Overlapping the snapshot of 265ns with the ternary complex (PDBID: 2HR1), the RMSD value of the catalytic loop is 2.703 Å, and the RMSD of overall structure is 1.521 Å (shown in Figure S2). As a result, our starting model and following simulation is reasonable. Whiel there are also some deficiency in it: since the ending of this simulation, we have not observed the hydrogen bonds between Arg240 and Guanine 3' to target cytosine, and the formation of catalytic pocket. On the other hand, target cytosine is not stable (shown in Figure S3). Target cytosine rotate around an axis which connecting N1 and N6 90 degrees, then this cytosine tilts about 60 degrees, and the N6 atom hydrogen bonded to phosphor atom of nucleoside 5' to cytosine. Additionally, Amide group of Gln237 have not hydrogen bonded to Ser85. As a result, the catalytic loop has not proceeded to closed form, and the catalytic pocket of M.HhaI has not formed. We consider these insufficiency is caused by limited time scale and force field parameter accuracy 11,12 . Figure S2. A Comparison between the 265ns snapshot (the ending of our simulation) and the crystal structure of ternary complex. Structure derived from simulation is colored cyan, while the crystal structure is colored magenta. The residues that located around catalytic pocket are drawn using stick.