Anti-methanogenic effect of rhubarb (Rheum spp.) - An in silico docking studies on methyl-coenzyme M reductase (MCR).

The present study explored anti-methanogenic properties of rhubarb compounds using in silico analysis on methyl-coenzyme M reductase (MCR) for identifying its anti-methanogen mechanism. To identify pharmacokinetics of 35 compounds from rhubarb, molecular docking and ADME analysis were performed against MCR using AutoDockVina, FAFDrugs3 and PROTOX programs. Docking results successfully indicated three possible candidate compounds 9,10-anthracenedione, 1,8-dihydroxy-3-methyl (-6.92 kcal/mol); phthalic acid isobutyl octadecyl ester (-5.26 kcal/mol); and diisooctyl phthalate (-5.61 kcal/mol) showed minimum binding energy (kcal/mol) with the target protein MCR which catalyze the biosynthesis of rumen methane. In conclusion, the identified compounds showed the most docking fitness score against the target methyl-coenzyme M reductase and the decrease in ruminal methane emission by rhubarb might be a result of these compounds by inhibition of methanogenesis.

Introduction

Ruminal methane emissions from livestock are responsible for a significant Green House Gas (GHG) production in the agriculture sector (Hristov et al., 2013), of which 90% of the results from microbial methanogenesis (McAllister et al., 2015). Ruminal methane productions not only affect the environment, but also serve as energy loss to animals by affecting the growth and milk production (Johnson and Johnson, 1995). Therefore, decreasing ruminal CH4 emissions would increase animal productivity and benefit the environment. Several methane mitigation strategies have been carried worldwide such as changes in management practices (Beaucheminet al., 2008, Morgaviet al., 2010); use of feed additives (Jayanegara et al., 2018) and plant secondary metabolite’s inhibiting the rumen methane emission (Bodaset al., 2008, Goel and Makkar, 2012, Demirtaset al., 2018).

In ruminants, rumen methanogens convert the H2 and CO2 (produced by bacteria, protozoa, and anaerobic fungi) into CH4 through methanogenesis pathway (Hydrogenotrophic) (Patraet al., 2010, Patra and Saxena, 2010, Cieslaket al., 2013). For this process (methanogenesis), methanogenic archaea require the methyl-coenzyme M reductase (MCR) for the formation of methane. The MCR used as a reliable marker of methanogenesis in diverse environments (Lutonet al., 2002, Palacio-Molina et al., 2013). The most common methanogens (hydrogenotropicarchae) are from the genus Methanobrevibacter, closely related to methane emissions (Danielssonet al., 2012, Shiet al., 2014).

Bioinformatics tools (CADD) ensure a great potential in not only reducing the cost but also the proficiency with which they can be designed. A number of novel tools and techniques have supported in the speeding up of drug discovery processes such as molecular docking, QSAR, and pharmacophore designing (Stalin et al., 2016). Docking analysis grants the scientist to virtually screen a database of compounds and envisions the strongest binders based on several scoring functions. It discovers ways in which two molecules such as drugs and a receptor protein (MCR) competent and dock to each other well. Similarly, a previous in silico study reported, the compound 3-nitrooxypropanol was found to be natural ligand with methyl-coenzyme M and able to decrease the rumen methane production (Duin et al., 2016).

According to our previous study, rhubarb decreased ruminal methane emission in vivo by reducing Methano brevibacter population but the mode of mechanism not reported (Kim et al., 2016) and in another study, the chemical composition (35 compounds) of rhubarb was reported (Arokiyaraj et al., 2017). To extend our research in methane mitigation strategies, we made a new approach to find the interaction between the phytochemical compounds and MCR for its anti-methanogenic mechanism using molecular docking techniques. Therefore, we investigated the in silico docking analysis of methyl-coenzyme M reductase with the Rhubarb compounds for its anti-methanogen mechanism.

Materials and methods

Ligand preparation and target protein structure

We selected 35 compounds as ligand molecules and their names were listed in the supplementary table 1 (Arokiyaraj et al., 2017). The energy of these compounds was minimized using open babel in PyRx 0.8 as a ligand for virtual screening analysis and in silico binding studies of the identified compounds with the receptor MCR (RCSB PDB- 1MRO) were determined. The retrieved protein structure was further used for active site predictions and ligand docking analysis. CASTp tool was used to predict the active site of the selected target proteins (Tian et al., 2018).

Molecular docking and virtual screening

Molecular docking simulation was performed using virtual screening tools such as AutoDockVina in PyRx 0.8 to find the potent drug-like molecules based on the energy scores as per the method and parameters (Morriset al., 2009, Trott and Olson, 2010, Dallakyan and Olson, 2015). Scoring function was calculated using the standard protocol of lamarckian genetic algorithm (Morris et al., 1998). The grid map for docking calculations was centered on the target proteins. From virtual screening analysis, the best successive hits of drug-like compounds were selected on the basis of higher scoring function and the interaction of ligand with all selected protein models was evaluated. The finalized selected molecules were again docked using Auto dock tools for the confirmation of the ligand-protein interaction sites and visualized by PyMol molecular graphics system (http://www.pymol.org). The overall studies were performed in Corei5-6200U, Intel processor CPU @ 2.3. GHz with 8 GB DDR3 RAM bundled with Windows 10 operating system.

In silico ADME prediction

For ADME property analysis, the final compound hits were used for prediction of the drug-likeness and pharmacokinetic properties. FAF Drugs-3 web server was used for evaluating ADME parameters under logP computation program XLOgP3 (Lagorce et al., 2008) using the Lipinski rule of five (LROF) physchem filter (Lipinski, 2004). Additionally, using the ProTox server, oral toxicity and drug-likeness were checked for the finalized compounds (Drwalet al., 2014).

Results

Ligand-Protein interaction

In this study, three ligands such as 9,10-anthracenedione 1,8-dihydroxy-3-methyl- (ligand 29), phthalic acid isobutyl octadecyl ester (ligand 31) and diisooctyl phthalate (ligand 33) obtained from the set of 35 compounds exhibited higher least energy (minimum binding energy) than other molecules (data not shown) to bind with the target protein MCR (Table 1). Hydrophobic interactions between the ligands and target protein showed in Fig. 1, Fig. 2, Fig. 3. The binding affinity values for the top three hits (ligands 29, 31 and 33) varied from −5.26 to −6.92 kcal/mol against the target protein MCR. The ligand 31 showed a docking score of −5.26 kcal/mol. The oxygen atom of the carbonyl group and the phenolic ring has been bind with ASN‘481 (Asparagine). The ligand 33 attainment the score of about −5.61 kcal/mol and one of the carbonyl group of the ester side chain bind with VAL‘482 (Valine) and ASN‘481; the oxygen atom attached to the carbonyl in the other side chain also bind with ASN‘481. Among 35 compounds, the compound 9,10-anthracenedione 1,8-dihydroxy-3-methyl- (ligand 29) showed a higher score (−6.92 kcal/mol) against MCR than others. The di-chelated carbonyl at C-10 bound with GLY‘397 (Glycine) and ARG‘401 (Arginine) and the phenolic hydroxyl at C-4 binds with SER‘399 (Serine) of methyl-coenzyme M reductase. Additionally, the center phenolic ring shows π-π interaction with TYR‘333 (Tyrosine) and PHE‘396 (Phenylalanine).

Table 1

Molecular docking of selected GCMS compounds with MCR protein.

Ligand no.	Compound name	Molecular formula	ProteinPDB ID	No of H bonds	Binding amino acid Residues	Binding Energy(kcal/mol)	Inhibition Constant uM	VDW_HB desolv_energy(kcal/mol)	Reff. RMSD	Ligand efficiency
29.	9,10-Anthracenedione, 1,8-dihydroxy-3-methyl-	C15H10O4	MCR (PDB ID: 1MRO)	3	GLY‘397/O, SER‘399/HN, ARG‘401/HN1TYR‘333, PHE‘396 (π-π interaction)	−6.92	8.53 (uM)	−6.62	31.66	0.36
31.	Phthalic acid, isobutyl octadecyl ester	C30H50O4	”	1	ASN‘481/2HD2	−5.26	140.46 (uM)	−9.33	53.83	0.15
33.	Diisooctyl phthalate	C24H38O4	”	2	ASN‘481/2HD2, VAL‘482/HN	−5.61	77.07(uM)	−7.86	41.85	0.20

9,10-Anthracenedione, 1,8-dihydroxy-3-methyl- (ligand 29) showed least energy value and also good ligand efficiency.

Fig. 1

The ligand 29 (9,10-Anthracenedione, 1,8-dihydroxy-3-methyl) with corresponding amino acid residues of methyl-coenzyme M reductase (PDB ID: 1MRO) (A); hydrophobic interactions between the ligand 29 and methyl-coenzyme M reductase (PDB ID: 1MRO) (B). The yellow dotted lines indicated the hydrogen bond interaction.

Fig. 2

The ligand 31 (Phthalic acid, isobutyl octadecyl ester) with corresponding amino acid residues of methyl-coenzyme M reductase (PDB ID: 1MRO) (A); hydrophobic interactions between the ligand 33 and methyl-coenzyme M reductase (PDB ID: 1MRO) (B). The yellow dotted lines indicated the hydrogen bond interaction.

Fig. 3

The ligand 33 (Diisooctyl phthalate) with corresponding amino acid residues of methyl-coenzyme M reductase (PDB ID: 1MRO) (A); hydrophobic interactions between the ligand 31 and methyl-coenzyme M reductase (PDB ID: 1MRO) (B). The yellow dotted lines indicated the hydrogen bond interaction.

ADME properties

Further, the three ligands were selected based on their binding affinity. Molecular properties such as total polar surface area, rotatable bonds, rigid bonds, octanol-water partition coefficient (LogP), molecular weight, hydrogen bond acceptor (HBA), hydrogen bond donor (HBD), water solubility (LogS) and stereocenters for the three ligands were presented (Table 2). The ligand 29 has lower rotatable bonds whereas ligands 31 and 33 have high rotatable bonds.

Table 2

Physio-chemical and ADME properties.

Ligand no.	Compound name	Molecular weight	Oral bioavailability	Rotatable bonds	Flexibility	Log P	HBD	HBA	Rings	Lipinski violations	Solubility (mg/l)	Stereo Centers
29	9,10-Anthracenedione, 1,8-dihydroxy-3-methyl-	254.24	Good	0	0	3.53	2	4	1	0	4913.02	0
31	Phthalic acid, isobutyl octadecyl ester	474.72	Good	23	0.74	10.21	0	4	1	1	204.40	0
33	Diisooctyl phthalate	390.56	Good	16	0.67	8.41	0	4	1	1	546.92	0

9,10-Anthracenedione, 1,8-dihydroxy-3-methyl- (ligand 29) showed good in oral bioavailability and satisfied the Lipinski rule of 5. HBD – Hydrogen Bond Donor, HBA – Hydrogen Bond Acceptor.

Higher LogP (partition coefficient) values of ligand 31 and 33 suggest hydrophobic nature than ligand 29 which meets one of the Lipinski criteria of LogP of <5. However, ligand 29 has a lower molecular weight than others with 2 HBD and 4 HBA without Lipinski violations. Toxicity evaluation by ProTox server which tests compounds through 1 and 6 for the most and the least toxicity levels respectively suggest that the ligand 29, with LD50 of 5000 mg/kg, is falling in toxicity class 5, and ligands 31 and 33 with LD50 of 1340 mg/kg fall in the toxicity class of 4 indicating their least toxic nature.

Discussion

The molecular docking of ligands onto the selected protein has turned out to be an effective method to analyze their docking patterns that give us a view about their binding affinity and corresponding inhibitory effect (Perola et al., 2004). Among three selected ligands, the ligand 29 expressed high minimal binding energy (−6.92 kcal/mol) and ligand efficiency with the active sites on MCR protein. This could be due to more number of potential hydrogen bonds. It is observed that hydrogen bonding plays a vital role as functional determinants of protein-ligand interactions especially in inhibition of a complex. It is noteworthy that these compounds show non-covalent interaction with target proteins and can serve as a new class of non-covalent inhibitors. The results are agreeing with Siwek et al., (2012) that topoisomerase IV showed strong hydrogen binding affinity with the ligand thio-semicarbazide.

Results of the molecular properties indicated that the ligand 31 and 33 had high rotatable bonds suggesting huge conformation space and flexibility and ligand 29 with lower rotatable bonds. Even though ADME and oral toxicity analyses by FAFDrugs3 and ProTox virtual tools on the selected ligands suggest their drug-like properties (acceptable bioavailability and solubility (LogS)), the ligand 29 has a higher solubility than the other two. In a previous study by Lu et al. (2004) reported that oral bioavailability is affected by the compound’s flexibility and by the number of rotatable bonds (<15).

Overall, the ligand 29 showed good compatibility as they satisfied Lipinski rule of 5 indicating that these compounds might exhibit orally active drug-likeness (Table 2) features such as permissible number of HBD (<5), acceptors (<10) and acceptable molecular weight (not more than 500 g/M) (Lipinski et al., 2001). Therefore, compounds possessing values for Lipinski’s rule of 5 in acceptable ranges can be observed possibilities to ensure the good intestinal absorption or permeation over the gut-blood barrier (Artursson et al., 2001).

Further, the study revealed that the ligand 29 showed hydrogen bond interactions with GLY‘397, SER‘399, ARG‘401 and π-π interaction with TYR‘ 333 and PHE‘396 which are the most active amino acids present in the MCR protein. The ligand molecules were significantly bound with lower energy value, good Inhibition Constant (uM) and ligand efficiency. Additionally, the ligand 29 has been found to be the very similar region of the active sites of MCR to that of 1-Thioethanesulfonic acid which was present in the crystal structure of the MCR (PDB ID: 1MRO).

Based on our docking result, ligand 29 (from Rheum sp. root) showed good binding interaction with MCR thereby reducing the methanogenesis activity, which has been directly proportional to increasing animal performance via distributing metabolic hydrogen to fermentation pathways further dynamically favorable to the animals (Guyader et al., 2017). In a previous study, Kung et al. (2003) have been reported that the supplementation of 9,10-anthraquinone in sheep has triggered a reduction in methane (CH4) production.

Anthraquinones are the major groups in quinones, and naturally present in huge amounts on the plant-based ethno-pharmacological activity (Thomson, 1971, Muelleret al., 1999, Odom, 1997), which as good antibacterial molecules through inhibiting the function of bacterial protein synthesis (Anke et al., 1980) and membranes (Chan et al., 2011). Similar to our result, 3-nitrooxypropanol has been reported for the specific reduction of methane emissions from ruminants and bind in the active site region as like as 9,10-Anthracenedione 1,8-dihydroxy-3-methyl- (ligand 29) (Duin et al., 2016).

Conclusion

Molecular binding interaction of an in silico analysis demonstrated that the Rhubarb compounds 9,10-Anthracenedione, 1,8-dihydroxy-3-methyl (−6.92 kcal/mol); phthalic acid isobutyl octadecyl ester (−5.26 kcal/mol); and diisooctyl phthalate (−5.61 kcal/mol) has more specificity towards the methyl-coenzyme M reductase binding site and could be a potent anti-methanogen inhibitor. This study concludes that three candidates have the potential for developing an anti-methanogenic drug among the compounds derived from the Rhubarb.