In silico docking of methyl isocyanate (MIC) and its hydrolytic product (1, 3-dimethylurea) shows significant interaction with DNA Methyltransferase 1 suggests cancer risk in Bhopal-Gas- Tragedy survivors.

DNA methyltransferase 1 (DNMT1) is a relatively large protein family responsible for maintenance of normal methylation, cell growth and survival in mammals. Toxic industrial chemical exposure associated methylation misregulation has been shown to have epigenetic influence. Such misregulation could effectively contribute to cancer development and progression. Methyl isocyanate (MIC) is a noxious industrial chemical used extensively in the production of carbamate pesticides. We here applied an in silico molecular docking approach to study the interaction of MIC with diverse domains of DNMT1, to predict cancer risk in the Bhopal population exposed to MIC during 1984. For the first time, we investigated the interaction of MIC and its hydrolytic product (1,3-dimethylurea) with DNMT1 interacting (such as DMAP1, RFTS, and CXXC) and catalytic (SAM, SAH, and Sinefungin) domains using computer simulations. The results of the present study showed a potential interaction of MIC and 1,3-dimethylurea with these domains. Obviously, strong binding of MIC with DNMT1 interrupting normal methylation will lead to epigenetic alterations in the exposed humans. We suggest therefore that the MIC- exposed individuals surviving after 1984 disaster have excess risk of cancer, which can be attributed to alterations in their epigenome. Our findings will help in better understanding the underlying epigenetic mechanisms in humans exposed to MIC.


Introduction
There is increasing evidence that exposure to toxic chemicals primarily affect the epigenome leading to the alteration of normal epigenetic process (Hou et al., 2011), In humans, chemically-induced epigenetic alterations has been reported to play an important role in cancer development and progression (Fraga et al., 2005;Baccarelli and Bollati, 2009;Shrivastava et al., 2013). Methyl isocyanate (MIC) is a toxic chemical used as an intermediate in the production of carbamate pesticides (HSFS, 2002;Senthilkumar et al., 2011;2012;2013). Accidental leakage of MIC from storage tank number E610 of Bhopal Union Carbide India Limited (UCIL) pesticide production plant on 3rd December, 1984 resulted in the mortality of 3600 inhabitants and about 2,00,000 were severely exposed. This catastrophe was caused by the introduction of water with 40 tonnes of MIC. Analysis conducted after Bhopal episode have concluded the 1
Almost three decades elapsed of Bhopal disaster, only a few efforts have been made to understand the cancer risk in the MIC-exposed survivors (Malla et al., 2011;Senthilkumar et al., 2011;2012;. Emerging epidemiological reports from Bhopal, India suggest diverse cancer pattern amongst the MIC-gas exposed survivors (Dikshit and Kanhere, 1999;Ganesh et al., 2005;Senthilkumar et al., 2011). It was therefore, worthwhile to investigate the underlying mechanisms of MIC on epigenetic alterations/DNA methylation using computational methods. In the present study, we aimed to understand the in silico interaction of MIC and its hydrolytic product (1,3-dimethylurea) with DNA methyltransferase (DNMT) domains, in order to reveal the cause for cancer risk in Bhopal population exposed to MIC during 1984.
DNA methylation is essential for the maintenance of normal epigenetic process (Miller et al., 1974). Aberrant methylation is the major cause of several cancers, shows both hyper-and hypo-methylation patterns (Robertson, 2005  contributes to oncogenesis (Dante et al., 1991;Jones and Martienssen, 2005). Global hypomethylation also results in genome wide instability that leads to oncogenesis (Guo et al., 2014). DNA Methyl Transferase (DNMT) is relatively large protein family that controls methylation and responsible for cell growth and its survival in mammals (Qin et al., 2011). Knockout of any regions of DNMT results in the aberrant development that causes cell death (Li et al., 1992;Okano et al., 1999). DNMT1 is abundantly found in the cell and maintains the methylation process during replication (Leonhardt et al., 1992). DNMT1 holds 1616 amino acids with N-terminal and C-terminal catalytic regions, linked together by seven lysyl-glycyl dipeptide repeats referred to as (KG)7 linker (Jurkowska et al., 2011). N-terminal region contains DNA methyl transferase associated protein 1 (DMAP1) interacting domain, proliferating cell nuclear antigen (PCNA) binding domain, nuclear localization signal (NLS), targeting sequence (TS) / Replication Foci Targeting Sequence (RFTS) domain, CXXC (Zinc Finger) domain, tandem bromo-adjacent homology (BAH) 1 and 2, and KG linker. N-terminal domain controls the activity of catalytic domain. The catalytic domain of DNMT1 forms a core complex that makes contact with the single side of DNA. DNMT 3a and b are involved in de novo methylation (Mortusewicz et al., 2005;Chen and Li, 2006). Recent studies have shown a significant association of DNMT3 in the maintenance of methylation and DNMT1 in de novo methylation (Jeltsch and Jurkowska, 2014). C-terminal region comprise a catalytic domain with different regions viz. I, IV, VI, and VIII-X, altogether constitutes methyl transferase motifs that helps in methyl group transfer from S-Adenosyl Methionine (SAM) to methyl cytosine through specialized base flipping mechanism (Jurkowska et al., 2011;Paganon et al., 2011;Clements et al., 2012). S-adenosyl homocysteine (SAH) (AdoHcy) is located at the active center of catalytic domain. CXXC-BAH1 domains are situated on different sides of the catalytic domain and connected by its long linker. BAH1 and BAH2 domains are located distantly from the bound DNA and separated from each other by α-helical linker (Figure 1) (Song et al., 2012).

Source, selection and preparation of protein
A total of four human DNMT1 proteins (viz. DMAP1, RFTS, DNMT1 in complex with DNA, and Sinefungin) and their three-dimensional (3-D) crystal structures were retrieved from the Protein Data Bank (PDB, 2014a-d) (ID: 4IEJ, 3EPZ, 3PTA, 3SWR) by X-Ray diffraction with the resolution range from 1.45 to 3.60 Å. The raw crystal structures were subjected to docking preparation using Protein Preparation Wizard (PPW). These raw crystal structures were scrutinized for excessive heavy atoms, improper bonding, missing of side chains, presence of co-crystallized ligand, water molecules, metal ions, and cofactors etc. Using PPW, hydrogen and appropriate bonds were added to the structures. Metals were treated by breaking bonds and assigned with formal charges. Disulfide bonds were added between two adjacent sulphur molecules. Overlaps were removed, water molecules were deleted, and missing loops were filled using the PRIME. N-Acetyl (ACE) and N-Methyl Amide (NMA) groups were added to uncap the N-and C-terminus. In addition, the hetro-groups were also detached. Protonation and metal charges states were generated. Hydrogen bonds were optimized using the PROT ASSIGN. Finally, the protein structure was subjected to restrained minimization in the Impref utility using the OPLS2005 force field. (Jacobson et al., 2004;Shelley et al., 2007;Shivakumar et al., 2010;Sastry et al., 2013).

Preparation of ligand
3-D ligand structures of MIC and 1,3-dimethylurea were retrieved from ZINC Database (Irwin and Shoichet, 2005). Ligand geometry was optimized by generating the structural variants, which reduced the structural problems in ligands and these structures were checked through LIGPREP. Ionization state was generated in the default pH 7.0±2.0 using ionizer. The structure was desalted to retain molecule with the largest number of atoms. Tautomers were generated for each neutralized/ionized molecule. Finally, around 32 stereoisomers were created with retained specific chiralities according to Greenwood et al. (2010).

Generation of receptor grid
The receptor grid is a 3-D boundary for ligand binding. The grid was generated specifically for various domains of the proteins. The residues of DMAP1 interacting with DNMT1 was determined using the structure (PDB ID: 4IEJ). The grid was created on the centroid region of selected residues in which DMAP1 and DNMT1 interact, in order to confine the ligands to that specific region during the docking process. Similarly, CXXC region interacting with DNA was analyzed using the structure (PDB ID: 3PTA) and grids were generated separately in three regions. RFTS Domain (PDB ID: 3EPZ) was also analyzed using the same procedure and the grid was generated for two different regions of DNMT1, where the interaction with replication foci occurs. With reference to catalytic domain of DNMT1, SAH is previously co-crystallized with the structure (PDB ID: 3PTA) and therefore the grid generation process was totally different. The grid was generated at the centroid of selected ligand with co-crystallized SAH. The grid was generated around this ligand and grid file also lack this ligand. Typically, using this procedure the grid was generated and ligands were docked only for this region. Similar procedure was adapted for Sinefungin (PDB ID: 3SWR). Van Der Waals scaling was made for all the grids using default scaling factor of 1.0 and partial cut-off charges were 0.25 as described earlier . In silico docking The docking sites for ligand were presented by grids. During the docking procedure of Schrodinger GLIDE, the ligands and grids were assessed and then subjected to induced fit flexible docking using Extra Precision (XP) mode. While docking the protein, the receptor was rigid, but the ligand was flexible. Van Der Waals radii scaling were done with the default scaling factor 0.8 and a partial charge cut-off as 0.15. Different conformations were created internally by GLIDE and subjected through a set of filters, viz. Euler angles, grid based force field evaluation and refinement, and Monte Carlo energy minimization according to Friesner et al. (2006). Finally, the docked conformers were estimated using GLIDE score and a single best position per ligand was generated as an output as described earlier (Halgren et al., 2004). For each ligand and protein site, the docked positions were determined accordingly to the docking profile (such as docking, XP GLIDE, and standard precision GLIDE scores). In addition, GLIDE-E model was calculated and it is directly proportional to the elevated negative value i.e. increased negative value is considered as a better docking score. This score is the measure of strength and stability of the bound-ligand-protein. Table 1 depicts the docking profile of MIC and 1,3-dimethylurea with DNMT1 (Figure 2A-H). Both the ligands showed significant interaction with amino acids and stable hydrogen bond formation with core-and sidechains on various regions of DNMT1. Hydrogen bonds formed by ligands with the protein residue were found between 1.5 to 2.5 Å. The minimum donor and acceptor angle is 120° and 90°, respectively. 1,3-dimethylurea has shown a strong interaction as compared to MIC. This interaction is due to the size and complexity of the molecule. Besides, the strongest binding was observed for SAH and Sinefungin region ( Figure 2 G and H).

Discussion and Conclusion
In silico docking reveals the effect of any chemicals and its binding affinity on 3-D conformational structures (Khamkar et al., 2013). With the use of computational methods, numerous studies have discovered that several cancers are markedly triggered by epigenetic alterations (Mahdavi et al., 2012;Wu et al., 2013;Jiang et al., 2014;Song et al., 2014;Zhuo e t al., 2014;Wang and Yan, 2015). The present study was also an attempt to determine the in silico interaction of MIC and 1,3-dimethylurea with several domains of DNMT1.
DMAP1 is a charge rich region that binds with MIC and 1,3-dimethylurea specifically interferes as a transcriptional repressor in the normal process. Interference of the repressive activity results in the aberrant expression of genes (Rountree et al., 2000). Furthermore, any impendence in the normal binding of DMAP1 to DNA at CpG sites may exhibit a lack of DNMT1 activity, which affects the gene expression (Bashtrykov et al., 2012). Besides, DMAP1 also form several complexes that participate in the replication process. DNMT1-DMAP1 complex is recruited in the replication foci during S-phase that facilitates the replication and also provides stability to DNMT1-histone deacetylase (HDAC) activity. Therefore, the binding of MIC and 1,3-dimethylurea with DMAP1 apparently hinders the normal methylation process that cause malfunction of transcription machinery (Rountree et al., 2000).
CXXC domain is a crescent shaped ZN+ binding domain, rich in cystine and contains 8 catalytically important cystine residues in the form CXXC motif. CXXC binds to CpG of DNA at both major-and minorgrooves and also links to the bases and phosphate groups (Pradhan et al., 2008;Risner et al., 2013). In the present study, we noticed the binding of both the ligands with various regions of CXXC, which may affects the normal methylation mechanism. RFTS domain targets DNMT1 in the Replication Foci during S-phase. Besides, it also induces DNMT1 to heterochromatin region from late S-phase to early G1-phase of the cell cycle. Therefore, the ligands interacting with RFTS may influence the DNA replication and cell cycle (Paganon et al., 2011). RFTS also contains Zn+ binding motifs, where its hydrophobic interactions cause dimerization of DNMT1. Ligand binding to RFTS inhibits the interaction, which results in to aberrant methylation process (Clements et al., 2012). SAM (Adomet) is considered as a key methyl donor in almost every methylation reactions of the cell. SAM is produced by single carbon metabolism network with the help of dietary folates (Lu, 2000). SAM acts a cofactor of DNMT1 and binds to its catalytic domain, which triggers the methylation reaction. During this process, the cytosine is flipped out from the DNA double helix, where the formation of covalent enzyme-substrate intermediates occurs as a result with the transfer of methyl group to the cytosine (Song et al., 2012). Thereafter, Adomet is converted into SAM, which is a potent inhibitor of methylation reactions (Finkelstein, 1990;Hermann et al., 2004). Both the MIC and 1,3-dimethylurea was found to bind strongly with this region, indicated the transfer of SAM-SAH may affect the normal methylation reaction.
Sinefungin, an antibiotic obtained from Streptomyces griseolus is an inhibitor of DNMT1. It acts antagonistically to SAM by blocking the methyl transfer reactions. Sinefungin binds strongly than SAM and SAH. Therefore, it is well-known to inhibit DNMT1 (Schluckebier et al., 1997;Kilgore et al., 2013). Both the ligands have shown significant interaction with this region, where as Sinefungin acted as an inhibitor. The preoccupied condition of MIC and 1,3-dimethylurea cause obstruction in the normal inhibitory mechanism of Sinefungin.
In agreement to our results, several studies proved that the global hypomethylation is due to many environmental and agricultural chemicals (Hou et al., 2011) have plausible role in cancer etiology. Interestingly, recent docking studies conducted by Shrivastava et al. (2010;2013) and Tripathi et al. (2015) also demonstrated a comprehensive understanding of MIC with different immunoproteins against Tuberculosis. To support our findings, recent epidemiological and experimental studies also indicate genetic instability in the MIC-exposed survivors, are highly susceptible to cancers (Malla et al., 2011;Senthilkumar et al., 2011;2013;. Almost 30 years elapsed after Bhopal disaster, owing to several bias and factors it is difficult to understand the direct or indirect effects of MIC in the exposed survivors. Hence, this kind of in silico docking analysis will be helpful to elucidate the underlying disease mechanism in the MIC-exposed survivors and their epigenetic status. To our knowledge, this is the first study to report the interaction of MIC and 1,3-dimethylurea with DNMT1 (DMAP1, RFTS, and CXXC) and catalytic (SAM/SAH and Sinefungin) domains.
In summary, we found a potential interaction of MIC and 1,3-dimethylurea with these domains. Apparently, the stronger binding of MIC with DNMT1 interrupts the normal methylation will lead to epigenetic alterations in the exposed survivors. We conclude therefore that the MIC-exposed individuals surviving after 1984 disaster have excess risk of cancer, which can be attributed to alterations in their epigenome. Our findings will help in better understanding the underlying epigenetic mechanisms in humans exposed to MIC.