Anti-Methanogenic Effect of Phytochemicals on Methyl-Coenzyme M Reductase—Potential: In Silico and Molecular Docking Studies for Environmental Protection

Methane is a greenhouse gas which poses a great threat to life on earth as its emissions directly contribute to global warming and methane has a 28-fold higher warming potential over that of carbon dioxide. Ruminants have been identified as a major source of methane emission as a result of methanogenesis by their respective gut microbiomes. Various plants produce highly bioactive compounds which can be investigated to find a potential inhibitor of methyl-coenzyme M reductase (the target protein for methanogenesis). To speed up the process and to limit the use of laboratory resources, the present study uses an in-silico molecular docking approach to explore the anti-methanogenic properties of phytochemicals from Cymbopogon citratus, Origanum vulgare, Lavandula officinalis, Cinnamomum zeylanicum, Piper betle, Cuminum cyminum, Ocimum gratissimum, Salvia sclarea, Allium sativum, Rosmarinus officinalis and Thymus vulgaris. A total of 168 compounds from 11 plants were virtually screened. Finally, 25 scrutinized compounds were evaluated against methyl-coenzyme M reductase (MCR) protein using the AutoDock 4.0 program. In conclusion, the study identified 21 out of 25 compounds against inhibition of the MCR protein. Particularly, five compounds: rosmarinic acid (−10.71 kcal/mol), biotin (−9.38 kcal/mol), α-cadinol (−8.16 kcal/mol), (3R,3aS,6R,6aR)-3-(2H-1,3-benzodioxol-4-yl)-6-(2H-1,3-benzodioxol-5-yl)-hexahydrofuro[3,4-c]furan-1-one (−12.21 kcal/mol), and 2,4,7,9-tetramethyl-5decyn4,7diol (−9.02 kcal/mol) showed higher binding energy towards the MCR protein. In turn, these compounds have potential utility as rumen methanogenic inhibitors in the proposed methane inhibitor program. Ultimately, molecular dynamics simulations of rosmarinic acid and (3R,3aS,6R,6aR)-3-(2H-1,3-benzodioxol-4-yl)-6-(2H-1,3-benzodioxol-5-yl)-hexahydrofuro[3,4-c]furan-1-one yielded the best possible interaction and stability with the active site of 5A8K protein for 20 ns.


Introduction
Global warming has become an increasing concern to sustaining life on our planet. Various factors have led to global warming, especially the emission of greenhouse gases 2 of 14 from various biotic and abiotic sources. One of the most concerning greenhouse gases is methane (CH 4 ) which has a 28-fold greater warming potential over carbon dioxide (CO 2 ) [1]. Though anthropogenic emission of greenhouse gases of 7% to 18% is generally accepted, the increasing trend of emissions has raised global concern [2,3]. When investigating the sources of methane emissions, a large contribution of emissions trace back to livestock, and especially to ruminants such as sheep and cattle. The production of methane by ruminants is referred to as ruminant methanogenesis. This methanogenesis is also a factor for loss of energy in ruminants which could have been otherwise used for growth development or that of milk production [4]. The process of methanogenesis in the ruminants is outlined as the conversion of CO 2 and H 2 into CH 4 which is carried out by the ruminant microbiome which includes archaea, bacteria, fungi, and protozoa [5,6]. Studies on the ruminant microbiome reveal that bacteria inhabit 95% of the microbiome and 2-4% of archaea and 1% of protozoa and fungi [7]. For methanogenesis, rumen archaea bacteria require the enzyme methyl coenzyme M reductase (MCR). Therefore, the MCR protein was used as a marker of methanogenesis in diverse environments [8][9][10]. The most common methanogens (hydrogenotropicarchae) are from the genus Methanobrevibacter, closely related to methane emissions [11,12]. Some of the microorganisms involved in ruminant methane production are Methanobacterium bryantii [8], Methanomicrobium mobile [9], Methanobrevibacter olleyae, Methanobrevibacter millerae [10], Methanobacterium formicicum [11], Methanobrevibacter ruminantium [12], and Methanobrevibacter gottschalkii [4]. These microorganisms possess MCR, which is involved in the methanogenesis metabolism and referred to as a hydrogenotrophic process [13]. MCR is a critical enzyme involved in the final step of nearly all methanogenesis-based metabolism, along with coenzyme M (CoM) which is 2-merceptoethane sulfonic acid. MCR aids in the reduction of methylated CoM, which results in the production of methane. This shows that methyl-CoM reductase is an important enzyme which plays a significant role in methanogenesis [14]. This has led researchers to find a significant target molecule for inhibiting this enzyme so as to reduce methane production in ruminants. A significant source of compounds for various applications have been provided by nature's best progeny: plants and their related phytochemicals. Globally, various plants have been studied and documented for various biological activity and other applications. Selected phytochemicals from various plants can be studied for their inhibitory activity against this critical enzyme as a solution for attenuating ruminant methane production. However, every plant possesses numerous phytochemicals, each with their respective functions and methane emission. To address the increasing global emissions in an effective manner, there is a need for promising and rapid testing methods to save time and resources for in-vitro or laboratory screening of these phytoactive compounds against MCR. Thus, in-silico studies can afford screening of numerous phytochemicals against this enzyme so as to find a potential inhibitor of MCR without the requirement of in-vitro laboratory resources. Molecular docking is one of the most used in-silico screening activities of any biomolecules or ligands with the target protein. Thus, the current study involves the process of selecting highly bioactive phytochemicals based on previous literature and screening of compounds for potential inhibitory activity against methyl-CoM reductase. The target phytochemicals with promising results can be used for in-vitro studies that can be further used as phytogenic feed additives for ruminants to ensure reduced methane emissions.

Plant Selection
An extensive literature study was carried out to identify suitable plants with highly bioactive phytochemicals that can be utilized as a potential inhibitor for rumen methane production. Initially, phytochemical databases such as Dr. Duke's Phytochemical and Ethnobotanical Dictionary of Natural Products 27.2 [15], IMPPAT: Indian Medicinal Plants, Phytochemistry And Therapeutics [16], Super Natural Database V2 [17] were taken into consideration and based on the collected literature data, 11 indigenous plants were chosen, as follows: Cymbopogon citratus, Origanum vulgare, Lavandula officinalis, Cinnamomum zeylanicum, Piper betle, Cuminum cyminum, Ocimum gratissimum, Salvia sclarea, Allium sativum, Rosmarinus officinalis and Thymus vulgaris. These plants were selected for rumen methanogenesis analysis.

Virtual Screening of Ligands
A virtual screening of ligands was carried out to scrutinize the database for the identification of compounds with higher bioactivity and low toxicity. The ADME properties of the compounds were predicted using the Swiss ADME web server and Lipinski rule of five along with "drug likeliness" predictions that were used to filter out the best compounds with required activity [29].

Molecular Docking
Phytochemical compounds which showed nil violation for Lipinski rule of five and moderate "drug likeliness" were selected for the process of molecular docking. Docking was carried out to evaluate the real time interaction between the phytochemical compounds with the target protein of interest (methyl-coenzyme M reductase, catalyst for rumen methane production) by measuring the binding energy of the complex (drug with target protein). AutoDock 4.0 (AD) and AutoDock Vina 1.5.6 (ADV) (Centre for Computational Structural Biology, La Jolla, CA, USA) were used for this process and the results were visualized using the Chimera UCSF visualization tool (Resource for Biocomputing, Visualization, and Informatics (RBVI), San Francisco, CA, USA).

Ligand Preparation and Target Protein
The three-dimensional structure of the target protein (methyl-coenzyme M reductase) was retrieved from the protein data bank using its PDB ID and it was visualized using discovery studio visualizer (DSV) software. Water molecules, ions, and standard inhibitors if present were also removed using the discovery studio visualization tool [30]. Marvin sketch software was used to construct the three-dimensional structures of the ligands to be studied. Avogadro (molecular modelling tool) was used for energy minimization of ligand molecules to obtain the best pose with low energy (Most stable 3D form). Both local charge and torsion (rotational motion) was provided to three-dimensional ligand structure for obtaining an effective binding complex with target protein. Similarly, Kollmann charge and polar hydrogen bonds were added to the protein [31]. Three dimensional PDBQT files of both the ligand and the protein were generated using MGL tools.

Protocol of Docking Studies
Auto dock: Auto dock version 4.0 software was used for automated docking studies. AutoGrid grid, a component of auto dock, was used to compute the grid maps with the interaction energies depending upon the macromolecule target of the docking study. The grid center was placed on the active target site of the protein (predicted using prankweb.cz online tool) [32]. Then, the binding free energy of the inhibitors was evaluated using automated docking studies. The best conformations search was done by adopting a genetic algorithm with local search (GA-LS) method. The docking parameters were set at default values with 100 independent docking runs using the software ADT (AutoDock Tool Kit). Root mean square (RMS) tolerance of 2.0 Å was performed using structures generated after completion of docking via cluster analysis. Molecular graphics and visualization were performed with the UCSF Chimera package [33,34].
Autodock vina: By using an Autodock vina software, binding energies were reported and determined in kcal/mol unit for ligand-protein. The grid box was generated by targeting the active site with a size of (x = 40, y = 40, z = 40) and a center of (x = 54.17, y = 9.53, z = 37.09) was set to cover the binding site for 5A8K protein, which was adopted and the conformations simulated. Virtual molecular docking and analysis were performed using the AutoDock Vina 1.5.6 tool. All the ligand molecules were docked to the active site of 5A8K protein. Finally, protein-ligand interactions were analyzed by using Biovia discovery studio visualizer. The top-ranked ligands with the most negative, favorable interactions were carried forward for further analysis of molecular dynamics simulations.

Molecular Dynamics Simulation Studies
The molecular dynamic investigation was performed for the top two compounds (3R,3aS,6R,6aR)-3-(2H-1,3-benzodioxol-4-yl)-6-(2H-1,3-benzodioxol-5-yl)-hexahydrofuro-[3,4-c]furan-1-one and rosmarinic acid, which were obtained from the molecular docking. This approach was used to analyze the ligand consistency and the binding mode stability in the binding pocket of protein [35]. To maintain the system's neutral conditions, the minimum quantity of Na + ions with salt ions were added. Once the compound attained energy minimization, it was put through equilibrium using NPT gathering for 2 ns [36]. The relaxed system was then submitted to 20 ns MD simulations, which were performed using Marlyna-Tobias-Klein barostats at 1 bar pressure and a Nosé-Hoover thermostat set to 300 K under an NPT ensemble. In each course of the box, the protein complexes had a minimum of 7 Å buffer to allow for significant fluctuations along the MD simulation. By using desmond modules of the Schrodinger 2020-2 suite software (New York, NY, USA), the molecular dynamics simulation study was carried out [37].

Plant Selection and Database Preparation
Phytochemical compounds recognized from GCMS reports of existing literature and collected from various databases were used to prepare a primary database for the study. Initially, a total of 168 potential chemical entities (data not shown) were identified as bioactive phyto-ingredients present in 11 indigenous medicinal plants. Two dimensional structures of the compounds were downloaded from pubchem and were compiled into a single structural database as shown in the Supplementary Material (cf. Figure S1).

Virtual Screening of Ligands
The virtual screening of ligands was carried out using the Swiss ADME server and parameters like Lipinski's rule of five and "drug likeliness" were used as a filter for narrowing down the number of compounds to be taken up for docking. Initially, among 168 compounds present in the database, a total of 51 compounds with less than three violations were filtered out. Then these compounds were subjected to the same Lipinski filter (without any violation) along with "drug likeliness" and finally 25 compounds were found to have nil violations in Lipinski rule and also zero to one violation for "drug likeliness", proving it as a lead target molecule for docking studies. The compounds selected through virtual screening are listed in Table 1.

Molecular Docking of Compounds
The three-dimensional structure of methyl-coenzyme M reductase (Figure 1) was retrieved from the protein data bank using the PDB ID 5A8K (https://www.rcsb.org/ structure/5A8K, accessed on 17 October 2021), which was prepared for docking analysis using discovery studio visualizer software. Additional data is presented in the Supplementary Materials (cf. Table S1) related to the hydrogen bond interactions of ligands with Methyl co-enzyme M reductase. The protein consists of six chains (chains A-F), whereas chain A, B, C represent a dimer unit of chain D, E, F. The 25 ligands taken for the study were also preprocessed by the addition of charge and torsion with the help of AutoDock tools 4.0. The active site of the protein was predicted using PrankWeb server and the selected residues were used as coordinates for generating the AutoGrid. The active site that was predicted by the PrankWeb server is tabulated in Table 2 and the active pocket is illustrated in Figure 1.

Sample Ligand Ki
Upon considering the various parameters of docking such as binding energy, the number of hydrogen bonds, and inhibition constant, it was clearly evident that the selected phytochemicals have greater specificity towards the methyl-CoM reductase binding site and could serve as potent anti-methanogen inhibitors. Micromachines 2021, 12, x 9 of 14

Conclusions
With a steady increase in greenhouse gas production and severe natural calamities due to climatic change arising because of global warming, there is a need to find green solutions to address this issue, hence efforts were made to correlate the methanogenic activity of bioactive phytochemicals present in various indigenous medicinal plants. This study was carried out to determine the methanogenic activity of 11 medicinal plants. Using an extensive literature survey, a total of 168 bioactive phytochemicals were taken into consideration and virtual screening was carried out with ADME and Lipinski rules of five as filters. Finally, 25 compounds were screened, and molecular docking was carried out for the same against methyl co-enzyme M reductase as the target protein. Most of the compounds showed good results, where five compounds had exceptional binding affinity with low binding energy and better hydrogen bonding. The results support that phytochemicals such as (3R,3aS,6R,6aR)-3-(2H-1,3-benzodioxol-4-yl)-6-(2H-1,3-benzodioxol-5-yl)hexahydrofuro [3,4-c]furan-1-one, rosmarinic-acid, biotin, 2,4,7,9-tetramethyl-5decyn-4,7diol, and α-cadinol can be utilized for real time methanogenic applications. MD simulations were carried out for two compounds ((3R,3aS,6R,6aR)-3-(2H-1,3-benzodioxol-4-yl)-6-(2H-1,3-benzodioxol-5-yl)-hexahydrofuro [3,4-c]furan-1-one and rosmarinic-acid). Improved docking energies revealed the stability of protein-ligand complexes and the existence of hydrogen bonds. However, further in-vitro and in-vivo studies are encouraged to confirm the anti-methanogenic ability of the selected plant compounds.
Supplementary Materials: The following are available online at https://www.mdpi.com/article/ 10.3390/mi12111425/s1, Figure S1. An illustrative example of the developed database for various ligand systems examined in this study. Table S1. List of ligands (n = 166) selected for the study.

Data Availability Statement:
The data presented in this study are available on request from the co-corresponding coauthor (R.J.). The data are not publicly available due to the raw/processed data required to reproduce these findings that cannot be shared at this time as the data also forms part of an ongoing study.