The Inhibition Mechanism of Emodin to the Core Protein of Aeromonas Hydrophila at the Molecular Level


 Aeromonas hydrophila is a main pathogen and is of great harm in aquaculture, but the effective methods to inhibit it is not clear. This study reveals that emodin inhibits the activity of A. hydrophila through experimental methods. Further proteomics analysis showed that atpE was the core protein of the inhibition. To disclose how emodin affects the inhibition mechanism, molecular docking and molecular dynamics simulations were carried out for A. hydrophila core protein binding with emodin. Comparing with these simulations of free atpE, this study uncovered that the inhibition of emodin involves emodin directly binding to core protein atpE with this simulations. The emodin engages with residues TYR10, ALA13, ALA14, MET17, LEU70 and TYR73 of atpE, which eventually leads to the inhibition of the activity of A. hydrophila. Additionally, this phenomenon was confirmed by the experimental results. The finding provides a novel insight for the study of the underlying mechanisms of emodin on antimicrobial activity of A. hydrophila at the molecular level.


Introduction
Aeromonas hydrophila widely exists in water and soil [1]. It can not only cause the outbreak of motile aeromonad septicemia in aquatic animals, resulting in a large number of deaths, but also cause human gastroenteritis and hemorrhagic disease, threatening human health. Previous researches report that A. hydrophila can induce sepsis, a serious disease in the sh farming industry, which has caused severe economic losses to sheries [2,3,4,5,6]. Many factors such as cytotoxins, adhesins, and hemolysins as well as metabolic pathways, bio lms formation speci c, and virulence factor expression regulation have been con rmed to be related to the virulence of A. hydrophila [7,8,9]. However, the inhibition mechanism of A. hydrophila at molecular level is still uncertain, it is an urgent need to probe the underlying mechanism to restrain pathogens.
Emodin (1, 3, 8-trihydroxy-6-methyl-anthraquinone) is a kind of monomer ( Fig. 1) isolated from rhubarb, Polygonum cuspidatum, Polygonum multi orum and other plants [10,11]. To date, emodin has growthinhibitory and cytotoxic effects on a variety of bacteria such as A. hydrophila [12,13]. At the same time, many studies have shown that emodin as a feed additive can enhance the immunity of sh to A. hydrophila infection [14,15]. However, most of the current research focuses on the changes in physical structures and bacteriostatic concentration, and lack of research on key targets. Therefore, it is important to understand the exact molecular mechanism by which emodin inhibits A. hydrophila and nd effective targets for emodin.
In recent years, computer-assisted drug design methods have developed rapidly, and these methods have achieved high accuracy that allows them to be routinely applied in the in-depth interpretation of experiments and drug discovery campaigns [16,17]. Protein structure-based methods are useful for predicting the binding modes of small molecules and their relative a nities. High-throughput molecular docking of up to 10 6 small molecules and scoring based on implicit solvent force elds can reliably identify micro-polar binding partners using rigid protein targets. Explicit solvent molecular dynamics (MD) simulation is a low-throughput technique used to characterize exible binding sites and accurately evaluate binding pathways, kinetics, and thermodynamics. This method can be used to guide experiments and explain experimental phenomena [17,18,19,20]. MD simulations have long been thought to provide insights into protein dynamics, not just insights into available crystal structures, and reveal new cryptic binding sites, extending the drug properties of targets. One of the rst methods described is the relaxation complex scheme, which combines all-atomic nanosecond-long molecular dynamics simulations of protein targets to describe its conformational exibility, while small molecules quickly dock to proteins stored in molecular dynamics snapshot [21,22]. This method has been successfully applied to a variety of cancer-related MDM2/MDMx-p53 interactions [23] or HIV integrase [24]. It is gratifying that in the latter case, the hidden pocket rst proposed in silico has been experimentally validated and approved for the development of clinical enzyme inhibitors against HIV infection [25].
Therefore, we performed molecular docking and MD simulations to assess the inhibition mechanism of A. hydrophila under the in uence of emodin on the basis of the experiments. The mode of binding of emodin to A. hydrophila core protein and interaction with residues was analyzed combing with experimental data. This study provides the underlying mechanism of inhibition behavior of the core protein of A. hydrophila under emodin at the molecular level.

RT-PCR analysis ad RNA extraction
The total RNA of A. hydrophila NJ-35 from the two groups (9 samples per group) were extracted by RNAiso Plus (TaKaRa, Japan) and tested by Nanodrop 2000 (Thermo Fisher Scienti c, USA). The RNA in each sample was diluted to 10 ng/ml. After that the quantitative analysis were performed on 2 µg of the total RNA using the Two Steps SYBR® PrimeScript® Plus RT-PCR Kit (TaKaRa, Dalian). Based on the method of Liu et al. [27], the Real-time quantitative PCR (RT-PCR) was adopted. ABI 7500 real-time PCR system was conducted for detection and the 2 −△△CT technique was used for further analysis.

Molecular docking calculations
The docking simulations were implemented by the molecular docking algorithm [28]. The atomic structures of emodin for the docking calculations were obtained from the protein data bank (http://www.rcsb.org/). Meanwhile, the homology model of atpE was obtained through the SWISS MODEL (https://swissmodel.expasy.org/). The potential CHRMM36 force elds were adopted, which was accurately adopt in different simulations based on experimental results [29].

Molecular Dynamics Simulations
MD simulations were performed with the GROMACS 5.1.3 through the potential of GROMOS96 54a7 force eld [30]. The lowest binding energy (the most negative) docking conformation generated by the software Autodock was taken as the initial conformation of MD simulation. The topological structure parameters of proteins were established by GROMACS program. The topological parameters of emodin were generated by Dundee PRODRG server. The simulation models were immersed in a tetrahedral box full of water. The solvation system (emodin, atpE and water) was neutralized by adding 2 Na ions. In order to release the collision contact, 50000 steps steepest descent method and 50000 steps conjugate gradient calculations were used to minimize the energy. MD simulation study includes equilibration stage and production stage. In order to balance the system, a 500 ps position constrained kinetic simulation (NVT and NPT) of solute (protein, anti-ion and emodin) was carried out at 300 K. Finally, the MD production of 100000 ps was carried out at 300 K and 1 bar pressure. The software of VMD is used for visualization [31].

MM/PBSA calculations
The binding free energies of atpE-emodin complex were obtained by the Molecular Mechanics Generalized Born Surface Area (MM/PBSA) methodology [32,33] based on where, ∆E MM is the energy of molecular mechanics interactions, including van der Waals energy and electrostatic energy. Where ∆G sol is the solvation free energy calculated via continuum solvent methods, which can be divided into nonpolar (∆G nonpolar ) and electrostatic (∆G polar ) contributions. The nonpolar component (∆G nonpolar ) was expressed as ∆G nonpolar = λSASA + b, where, b = 0.92 kcal mol − 1 and λ = 0.00542 kcal Å −2 [34], this component is a function of solvent accessible surface area (SASA). ∆G polar is calculated through solving the Poisson-Boltzmann equations. Where -T∆S is the entropy change, which was ignored here due to the low prediction accuracy and large computational overhead [35,36]. The entropy term also was ignored. Thus, the free energy obtained is not true binding energy. Besides, a perresidue decomposition of the total energy by the MM/GBSA algorithm was also carried out to assess the contribution of each residue to the binding. The polarity contribution of desolvation was calculated by GB model [37]. A more detailed theory of energy decomposition is described elsewhere [38].

Results And Discussions
3.1 Emodin inhibits the activity of A. hydrophila As shown in Fig. 2, with the increasing of the emodin concentration, the concentration of A. hydrophila decreased in a dose-dependent manner. The OD600 under the concentrations of 0 and 100 were almost the same. However, the concentrations of 200, 300, 400, and 500 µm/ml emodin groups were signi cantly different from the control group, and there were also signi cant differences between them, indicating that emodin can inhibit the growth of A. hydrophila. According to Fig. 2, the relative expression of atpE was at high level, but the expression was signi cantly down-regulated under the action of emodin, which indicates that atpE is a target site of emodin binding.

The binding mechanism of atpE with emodin
Computer-aided drug design was further performed to characterize the detailed mechanisms of this inhibition. Yet, many researches have shown that atpE at the protein was down-regulated under emodin treatment, the molecular mechanisms of inhibition expression remain elusive. Here, we hope to determine the potential binding of emodin to atpE by molecular docking and MD simulation. In order to study the mechanism of emodin inhibiting atpE, the preferred binding mode of atpE and emodin was determined by 10 ns MD simulation based on the docking results of emodin with atpE. In this calculations, emodin represents a kind of ligand, which binds to atpE through molecular mechanical interaction. During the simulation, emodin was localized to the region of atpE. VMD software was used to predict the binding mode of emodin with atpE and the binding energy of each residue [31] as illustrated in Fig. 3a. Figure. 3b showed that the key residues around emodin. Specially, the binding models (Fig. 3b) revealed that the TYR10, ALA13, ALA14, MET17, TYR73 and LEU70 anchor emodin. As shown in Fig. 3b, the chains of TYR10 and LEU70 were close to the emodin, indicating that TYR10 and LEU70 have strong interactions with emodin. Moreover, Fig. 3b shows that the benzene ring plane of emodin in atpE complex is parallel to that of residue. TYR10 and TYR10 were also close to emodin. Therefore, a strong interaction between emodin and these residues exists. These calculations will be further proved by the analysis of energy decomposition.
The root meant square uctuation (RMSF) of all residues of free atpE and atpE-emodin complex (residues 1-79) was shown to uncover the exibility of all the residues. As shown in Fig. 4, the RMSF of these residues clearly depicted the different exibility of atpE binding sites in the complex. In the case of RMSF less than 0.60 nm in the complex, all residues in the atpE binding sites (residues 10-20 and 60-75) of emodin exhibited small degree of exibility, showing that these residues became more rigid due to their binding to emodin. The regions are colored with yellow that binds with emodin show small degree of exibility with low RMSF (less than 0.6 nm) when compared with free atpE. This indicates that these residues seem to be more rigid because of binding to atpE. Figure. 3b shows that the stability of atpE binding cavity in the complex was mainly due to residues TYR10, ALA13, ALA14, MET17, LEU70 and TYR73.
The root-mean-square deviation (RMSD) of the backbone atoms with regard to the initial con rmation was obtained as a function of simulation time. The value of RMSD is to evaluate the conformational stability of the complex during the simulations. Figure. 5a showed that the mean value of RMSD pro le of the free atpE is 0.32 nm, which is higher than complex of 0.29 nm for the simulation process. The result showed that the conformation of the MD simulations after energy minimization was reliable for analyses. This indicates that the atpE in the complex is very rigid due to the constraint of emodin. On the other hand, Fig. 5b showed that the radius of gyration of complex bound conformation was found less than free atpE, indicating that the complex with the low radius of gyration and more stable.

The binding site in the atpE-emodin complex
In order to access the energetic effects on the contributions of key residues, the interaction energy between atpE and emodin was modeled based on the energy decomposition of the residues (Fig. 6a). The main energy contributions from − 8 to -1.5 kJ mol − 1 mainly came from TYR10, ALA13, ALA14, MET17, LEU70 and TYR73. These data agree with the described in Fig. 3 of the MD simulation results. Figure 6b showed that the binding mode of emodin with atpE. In this binding mode, color from red to blue represented total free energy from low to high. The results showed that emodin and atpE had a good binding effect.
In order to obtain underlying mechanism about the residues around the binding sites and their contribution to the complex, the molecular mechanics, polar, nonpolar and total energy contribution of the residues to the binding free energy were calculated by MM/PBSA algorithm. The analysis was performed using 1000 snapshots obtained from the 10 ns simulations. The sum of the each residue free energies were separated into molecular mechanics (∆E MM ), polar (∆E Polar ), nonpolar (∆E APolar ), and total contribution (∆E Total ). The energy decompositions from the key residues were showed in Fig. 7. We can see that the atpE-emodin complex, TYR73 has an appreciable molecular mechanics (∆E MM ) contribution, with a ∆E MM of -12 kJ/mol. In fact, TYR73 was close to the ligand of emodin, and more molecular mechanics interactions exist between atpE and emodin. Besides, residue TYR10, with a ∆E MM of -8.5 kJ/mol, has a strong molecular mechanics interaction with the ligand due to the close proximity between the emodin and residue TYR10, which facilitates a strong pep interaction in the atpE-emodin complex. For all the residues, the major decomposed energy interaction originated from molecular mechanics interactions, the second was by the polar interactions, while nonpolar contribution showed a minor effect on the formation of key residues in atpE-emodin complex. With the summation of the other energy term, an estimated ∆E bind of -0.271 kJ/mol was showed for emodin, indicating that emodin can strongly interact and bind with the binding site of atpE. These simulation results well explain the mechanism of emodin in inhibiting the activity of A. hydrophila.

Conclusions
In this paper, the growth curves of A. hydrophila at 0, 100, 200, 300, 400, and 500 µm/ml emodin were determined, and the gene expression level of atpE was measured. The results showed that emodin has a good inhibition effect of A. hydrophila, and atpE was the core protein of the inhibition. In order to understand the effect of emodin on the inhibition mechanism, MD simulations were performed for A. hydrophila core protein blinding with emodin. By this method, the analysis showed that the inhibitory of emodin involves the direct binding of emodin to the core protein atpE. The emodin engages with residues TYR10, ALA13, ALA14, MET17, LEU70 and TYR73 of atpE, which eventually led to inhibition of the activity. This mechanism was veri ed by the experimental results. The nding provides a novel insight for the study of mechanism of emodin on antimicrobial activity of A. hydrophila at molecular level.

Competing Financial Interests Statement
The authors declare no competing nancial interests. Figure 1 Chemical structure of emodin with the atom and moiety names used in this study    The RMSF of all residues (1-79) within the 10 ns simulation with regard to the initial positions of the atpEemodin complex. The region of the protein backbone appeared to depend on local environment uctuations (bound to emodin).