Inverse molecular docking reveals a novel function of thymol: Inhibition of fat deposition induced by high‐dose glucose in Caenorhabditis elegans

Abstract As a natural product isolated from thyme oil in thyme, thymol (2‐isopropyl‐5‐methylphenol) harbors antiviral, antioxidant, and other properties, and thus could be potentially used for the treatment of various diseases. However, the function of thymol has not been comprehensively studied. Here, we applied an inverse molecular docking approach to identify unappreciated functions of thymol. Potential targets of thymol in humans were identified by the server of DRAR‐CPI, and targets of interest were then assessed by GO and KEGG pathway analysis. Subsequently, homologous proteins of these targets in Caenorhabditis elegans were identified by Blast tool, and their three‐dimensional structures were achieved using Swiss‐Model workspace. Interaction between thymol and the targeted proteins in worms was verified using AutoDock 4.0. Analyses of the targets revealed that thymol could be potentially involved in the glycolysis/gluconeogenesis and fatty acid degradation pathways. To verify the activity of thymol on lipid deposition in vivo, the C. elegans model was established. The lipid content of nematodes induced by high‐dose glucose was determined by Oil Red O and Nile Red staining, and gene expression was assessed by qRT‐PCR. The results showed that thymol might lead to the acceleration of β‐oxidation by upregulating cpt‐1, aco, fabp, and tph‐1, causing the descent of lipid content in nematodes. Our findings indicated that thymol could be potentially used for the treatment of chronic metabolic diseases associated with increased fatty acid deposition.


| INTRODUC TI ON
Metabolic diseases, such as obesity and type-2 diabetes, have been becoming severe health problems. Obesity is considered to be the result of multiple genetic factors and environmental influences. In modern society, foods and drinks rich in high energy including fat and sugar can be seen everywhere. A high-fat and high-sugar diet has become the main diet structure at present and is considered as one of the main causes of obesity caused by the environment.
However, sugar-induced obesity is different from fat-induced obesity. On the one hand, sugar can inhibit fat degradation; on the other hand, it also can be oxidized and decarboxylated to acetyl-CoA that is the material for the synthesis of fatty acid. As we have known that glycolysis glucose is metabolized to acetyl-CoA from pyruvate by pyruvate dehydrogenase under aerobic conditions, and acetyl-CoA is the substrate of lipid biosynthesis. The rising blood glucose level triggers the release of insulin, leading to the dephosphorylation and activation of acetyl-CoA carboxylase (ACC) by insulindependent protein phosphatase. ACC catalyzes the formation of the first intermediate "malonyl-CoA" of fatty acid synthesis, which suppresses carnitine acyltransferase I, thereby preventing fatty acid entry into the mitochondrial matrix for degradation through βoxidation (Zheng & Greenway, 2012).
Recently, more national authorities and the World Health Organization (WHO) have provided advice on daily diets urging a limited consumption of rapidly digestible starches and sugars (Brouns, 2018). A low-energy diet and regular physical exercise are generally considered safe ways to prevent obesity. Besides, star molecular conjugated linoleic acid (CLA) has been widely applied to adjuvant therapy for chronic metabolic disease. CLA plays a role in the reduction of lipids storage, synthesis, and adipogenesis in adipocytes, and the enhanced fat utilization in muscle via fatty acid β-oxidation (Park & Pariza, 2007). However, some researchers have expressed concerns regarding the safety of CLA for humans (Kelley & Erickson, 2003;Larsen et al., 2003;Riserus et al., 2002;Tricon & Yaqoob, 2006). The main concerns were the involvement of CLA in fatty liver, lipodystrophy, oxidative stress, and glucose intolerance (Clement et al., 2002;Pariza, 2004;Poirier et al., 2005;Tricon & Yaqoob, 2006). Therefore, it is essential to identify alternatives to interfere with fat deposition.
Natural essential oil constituents present promising potential to treat infection , cancer (Xie et al., 2017), chronic inflammation (Yao et al., 2018), and even Alzheimer's disease (Postu et al., 2019). These constituents cover a group of biologically active molecules, including monoterpenes, sesquiterpenes, oxygenated monoterpenes, oxygenated sesquiterpenes, and phenolics among others. As one of the major constituents of essential oil in Monarda punctate, thymol (2-isopropyl-5-methylphenol) has antimicrobial (Marchese et al., 2016), antioxidant (Luna et al., 2017), antiviral (Lai et al., 2012), antitumor (Gunes-Bayir et al., 2019), and anti-inflamatory (Yao et al., 2018) properties. More importantly, evaluations from the Joint FAO/WHO Expert Committee on Food Additives (JECFA) have indicated that there are no safety concerns on thymol in terms of acceptable daily intake (ADI). Although multiple functions of thymol have been revealed, a comprehensive analysis of thymol's functionality has not been conducted. Computation methodologies for the repositioning of drug molecules onto their receptors have been extensively utilized to reveal novel functions of small chemical compounds. Through reverse screening, we found that thymol may be involved in the regulation of obesity. At present, studies on the effect of thymol on obesity are few and mainly focus on lipid-induced obesity. The effect of thymol on sugar-induced obesity has not been studied. This study will start from this point, and further through the analysis of the glucose metabolism pathway to find the potential target of the action of thymol.

| Strains and chemicals
The wild-type Caenorhabditis elegans strain N2 was obtained from the C. elegans Genetics Center (University of Minnesota, Minneapolis, MN, USA). The eggs of worms were incubated on nematode growth medium (NGM) agar plates at 20°C and seeded with a live Escherichia coli OP50 (E. coli OP50) as food resources.

| Inverse molecular docking of thymol with human proteins
DRAR-CPI server, a tool of inverse molecular docking, was developed to predict the drug target based on structural characteristics of protein (Luo et al., 2011). According to the algorithm in molecular docking program DOCK, the potential binding sites of receptor molecule in DRAR-CPI are defined by overlapping spheres, consequently, achieving the reverse docking of the drug with these binding sites. Affinity for interaction energy in ligand receptors is evaluated by scoring function and then ranked to predict the potential targets.
In this server, potential targets of a ligand can be obtained by submitting a mol2 formatted ligand file. In our study, the mol2-format thymol file was downloaded from the following ZINC online website: (https://zinc.docki ng.org/subst ance/1529592).

| Functional analysis for target genes
A PPI network was constructed to reveal the hub genes of the potential target proteins on STRING (https://strin g-db.org/cgi/input. pl), a web portal for undermining the integrated function of multiple genes (Szklarczyk et al., 2015). Besides, on STRING, Gene Ontology (GO) enrichment tool including biological processes (BPs), cellular components (CCs), and molecular functions (MFs) was used to perform functional analysis of target proteins. Similarly, the Kyoto Encyclopedia of Gene and Genomes (KEGG) enrichment tool was used to analyze biological pathway enrichment to identify those pathways that responded to thymol. It should be noted that the protein name achieved from DRAR-CPI needs to be changed to a discernible gene name by UniProtKB before submitting to STRING.
The correlation degree for the proteins of interests was analyzed by MCODE, a Cytoscape plugin that can be used to identify clusters (highly interconnected regions) in a network.

| Analysis of inverse docking and targeting sites prediction
Inverse docking was analyzed by Autodock 4.0, including AutoGrid and AutoDock. AutoGrid calculates related energy in a grid containing Van der Waals force, electrostatic force, and hydrogen bond force, while AutoDock conducts conformational search and evaluation using the Lamarckian Genetic Algorithm. For example, the docking of thymol with GAPDH could be realized by the following steps. The first step is the preparation of the system's coordinate files (PDBQT file) of GAPDH and thymol, including the removal of H 2 O (Eliminate the interference of solvation effect) and the addition of hydrogens (X-ray diffraction technology cannot obtain the coordinate data for H) in GAPDH molecule, and the setting of active keys number in thymol. The second step is the setting of search sites, mainly including the adjustment of grid size in "Grad Box." The third step is the establishment and operation of Autogrid parameter file with GPF format for thymol, generating GLG file and a series of MAP files of GAPDH that contain information of van der Waals forces, electrostatic forces, and dissolvent forces of different atom probe et al. The fourth step is the construction and running of Autodock parameter file with DPF format for thymol. Finally, we will achieve a DLG file of thymol. Docking files of GAPDH and thymol should be prepared for Pymol analysis, including the selection of thymol with the lowest energy, the sequence matching for GAPDH and thymol in Notepad++. In the end, a PDB file from Notepad++ needs to be opened in Pymol for the visualization and analysis of GAPDH-thymol docking. We can adjust the appearance of their conformation and search targeting site of thymol by observing its hydrogen bonding site in Pymol.

| Fat staining in Caenorhabditis elegans
Worms at the characteristic L4 larvae stage were chosen for synchronization. The procedure for glucose induction treatment was carried out as described previously (Zhao, 2015). Briefly, 1 mM and 5 mM glucose were added to the 50°C molten NGM medium. To determine the effect of thymol on lipid deposition, 1 μg/ml and 5 μg/ml of thymol were added to the NGM medium. The lipid content in the nematodes was evaluated by Oil Red O (ORO) (sigma) staining as described previously (Zhao, 2015). Caenorhabditis elegans at L4s were permeabilized through fixation, freezing and thawing, and dehydration process.
Permeabilized worms were stained with ORO in the dark overnight at 1.07 g, 28°C. Decolorization of nematodes was then performed with 1, 2-propylene glycol and washed with PBS buffer. The retention of ORO in nematodes could be observed with a 5× anatomical lens. The quantitative analysis was carried by the method of elution by ethanol.
Specifically, quantified nematodes were eluted with absolute ethyl alcohol after discarding PBS buffer by centrifugation, and the collected supernatant was transferred to a 96-well plate to measure the adsorption of ORO at 510 nm using a microreader.
Nile Red (Sigma) staining was also performed as described (Fouad et al., 2017). Briefly, 0.5 ml of the freshly diluted Nile Red (1 μg/ml) was added to plates seeded with E.coli OP50 and allowed to equilibrate for a minimum of 2 hr. Then, L1s obtained from bleaching were transferred to the plates and allowed to grow for the desired length of time. Stained nematodes were imaged under a 10× objective under a dark field and red fluorescence illumination with an inverted fluorescent microscope.

| RNA isolation and qRT-PCR
The total RNA of nematodes was isolated using TRIzol ® Reagent (Invitrogen) and purified with DNase I, RNase-free (Thermo Scientific™). cDNA was obtained using M-MuLV First Strand cDNA Synthesis Kit (Sangon). Real-time PCR was performed using SYBR GreenⅠprobes by QuantStudio ® 3 System (Peng et al., 2016). Primer sequences were listed in Table S1. Relative expression of genes was determined using the 2-ΔΔCt method and act-1 was chosen as the reference gene.

| Statistics analysis
All assays were repeated independently at least 3 times, and the experimental data were represented as mean ± SD. The data between two different treatments were compared statistically by one-way ANOVA in SPSS. The significance of differences was defined at the p < .05 level and p < .01.

| Protein targets prediction and analysis
Analysis by DRAR-CPI servers revealed that 392 human proteins have potential interaction with thymol. GO and KEGG analyses indicated that these targets of interest in glycolysis/gluconeogenesis and fatty acid degradation pathway were associated with carbon metabolism, tyrosine metabolism, pyruvate metabolism, retinol metabolism, drug metabolism-cytochrome p450, metabolism of xenobiotics by cytochrome p450, biosynthesis of amino acids and chemical carcinogenesis pathways (Figure 1a). Furthermore, in GO enrichment, these targets were mainly localized in the cytosol ( Figure 1b) and were associated with coenzyme binding, oxidoreductase activity, catalytic activity, oxidoreductase activity, acting on the CH-OH group of donors, NAD or NADP as acceptor and alcohol dehydrogenase activity, zinc-dependent (Figure 1c), and involved in organic substance catabolic process, single-organism biosynthetic process, small molecule metabolic process, xenobiotic process, ethanol metabolic process, pyruvate metabolic process, glucose metabolic process, oxidation-reduction process, and ethanol process (Figure 1d).
Our analysis also showed that module 1 (Figure 1e) was important in the whole PPI network, with the highest score of 19.714. GAPDH (Glyceraldehyde-3-phosphate dehydrogenase) is a node in module 1 with the MCODE-score of 12.72, which has the highest complexity in 12 proteins responding to thymol, predicting it may be the central protein in glycolysis. Independent of its glycolytic activity, the encoded protein has additionally been identified to have uracil DNA glycosylase activity in the nucleus (Meyer-Siegler et al., 1991). Also, GAPDH harbors a region that confers antimicrobial activity against Staphylococcus bacteria (Ong et al., 2019). Studies of a similar protein in mice have as-  Table 1.
Homology modeling is a technique that can build a 3D structure for proteins based on the 3D structure of a similar protein. A homology model can provide an alternative for subsequent receptor-ligand analyses such as molecular docking.  (Hao et al., 2015) No 1HSO Alcohol dehydrogenase 1A −23.3497 Catalyze the oxidation of alcohols to aldehydes. Variant confers susceptibility to esophageal squamous cell carcinoma (Cui et al., 2018). Polymorphisms in ADH1B associated with the increased risk of gastric cancer (Ghosh et al., 2017).  (Chen et al., 2015), and promotion of pancreatic cancer progression (Cui et al., 2015).

2VGB
Pyruvate kinase PKLR −24.1893 Catalyze the transphosphorylation of phohsphoenolpyruvate into pyruvate and ATP, which is the rate-limiting step of glycolysis. Correlate the promotion colorectal cancer liver colonization (Nguyen et al., 2016).

No
In our study, we constructed models of these homology proteins by SWISS-MOL and examined interactions between thymol and these proteins through Autodock 4.0 ( Figure 2). Docking analysis revealed that thymol potentially interacts with the different amino acids of these proteins (except 1DH-1) in nematodes (Table 2).

| Thymol reduces fat deposition due to high glucose exposure in nematodes
Fat accumulation induced by high-dose glucose (1 mM or 5 mM) in nematodes is inheritable across several generations (Zhao, 2015). So we chose 1mM or 5mM of glucose to induce fat deposition in nematodes. ORO allows for the quantification of total lipid levels, while Nile red facilitates the evaluation of lipid distribution among tissues such as the hypodermis, intestine, and the germline. When used in conjunction, ORO and Nile Red enable the researchers to determine how genotype and environmental changes affect lipid accumulation and wherein the worm these changes are occurring (Wilber et al., 2017).
The statistic result of ORO and Nile Red staining suggests both 1mM and 5mM of glucose are capable of inducing fat accumulation in the second generation of nematodes, and 5 mM-induced fat accumulation is more significant than 1mM-induced in the nematode (Figure 3b).
Total lipid levels and lipid droplets were visualized by staining with ORO ( Figure 3a) and Nile Red (Figure 3c), respectively. Further, the treatment of thymol can reduce lipid content. Importantly, 1 μg/ml and 5 μg/ml thymol treatment significantly reduced lipid deposition in glucose-induced nematodes, indicating thymol can reduce fat content.

F I G U R E 2
Homologous modeling of homologous proteins in nematodes at SWISS-MODEL network, evaluation of their 3D structure by SAVES, and molecular docking between target proteins and thymol through AutoDock 4.0

| Thymol restores the change of transcriptional level of genes in β-oxidation induced by highdose glucose
Our experimental data verified thymol could reduce lipid deposition in worms. We mentioned above that glucose-induced fat accumulation is on the one hand due to the synthesis of acetyl-CoA from glucose through glycolysis, which is the substrate of lipid biosynthesis, on the other hand, is due to the inhibition of fat degradation by the activation of acetyl-CoA carboxylase (ACC) triggered by the rising blood glucose level. Therefore, we suspected that the reduction of fat accumulation by thymol might be achieved by regulating the upstream glycolysis pathway and downstream of fatty acid β oxidation.
To confirm this idea, we determined the transcriptional level of key genes correlating with β-oxidation including carnitine palmitoyltransferase (cpt-1), acyl-CoA oxidase (aco), fatty acid-binding protein (fabp) in the process of β-oxidation and tryptophan hydroxylase (tph-1) in the upstream pathway of β-oxidation. The result of qRT-PCR suggested that the decreased level of cpt-1, aco, fabp, and tph-1 in nematodes exposed by 1 mM glucose was upregulated after the exposure of thymol ( Figure 4). These findings suggest the reduction of fat deposition caused by thymol is related to the oxidation of fatty acids.

| D ISCUSS I ON
In silico approaches can improve our ability to predict drug targets, thereby streamlining and accelerating the discovery process of a novel function of the drug. In this study, we found that thymol has a potential interaction with proteins associated with glycolysis/glu-   (Mak et al., 2006) and down-regulated after treatment with the hesperidin (Peng et al., 2016). In nematodes, four genes designated gpd-1 through gpd-4 encode glyceraldehyde-3-phosphate dehydrogenase, and △gpd-1 can cause greatly raised triglycerides (Hegele et al., 2014). However, there is no known link between the remaining genes (including pyk-2, tpi-1, 1dh1, to2g5.7, d2063.1, h24k24.3, and sodh-2) and obesity. These provide us with proof that this method has reliability for the prediction of thymol.
Our experimental results in vivo verified that thymol inhibited lipid deposition induced by high-dose glucose in nematodes. To understand how thymol makes a function in inhibiting fat accumulation, we detected the expression level of key genes correlation with β-oxidation by qRT-PCR. The findings suggested that thymol could recover the changes in the expression levels of these genes induced by high-dose glucose, including cpt-1, aco, fabp, and tph-1. So thymol probably accelerates β-oxidation according to the determination of qRT-PCR combined with bioinformatics prediction. Furthermore, we suppose that thymol may accelerate glucose metabolism by regulating the targets of GPD-1, GPD-2, GPD-3, GPD-4, PYK-2, 1DH-1, and 1MO0 in the glycolytic pathway, causing the reduction of blood glucose level and ACC inactivation. ACC inhibition resulted in the reduction of malonyl-CoA, which releases carnitine palmitoyltransferase (CPT) from inhibition and causes the transport of fatty acids into the mitochondria for β-oxidation. And the protein targets we predicted including H24K24.3, KAT-1, and T02G5.7 in the fatty acid degradation pathway may be located downstream of key genes correlating with β-oxidation ( Figure 5).

CO N FLI C T S O F I NTE R E S T
The authors declare no conflict of interest.

E TH I C A L S TATEM ENTS
This study does not involve any human or animal testing.

R E FE R E N C E S F I G U R E 5
Diagram of thymol alleviating the obesity induced by glucose