Evaluation of potential molecular interaction between quorum sensing receptor, LuxP and grouper fatty acids: in-silico screening and simulation

Pathologically relevant behaviors of Vibrio, such as the expression of virulence factors, biofilm production, and swarming motility, have been shown to be controlled by quorum sensing. The autoinducer-2 quorum sensing receptor protein LuxP is one of the target proteins for drug development to suppress the virulence of Vibrio. Here, we reported the potential molecular interaction of fatty acids identified in vibriosis-resistant grouper with LuxP. Fatty acid, 4-oxodocosahexaenoic acid (4R8) showed significant binding affinity toward LuxP (−6.0 kcal/mol) based on molecular docking analysis. The dynamic behavior of the protein–ligand complex was illustrated by molecular dynamic simulations. The fluctuation of the protein backbone, the stability of ligand binding, and hydrogen bond interactions were assessed, suggesting 4R8 possesses potential interaction with LuxP, which was supported by the low binding free energy (−29.144 kJ/mol) calculated using the molecular mechanics Poisson–Boltzmann surface area.


INTRODUCTION
The aquaculture of brown-marbled grouper suffers from a high frequency of vibriosis outbreak, which often causes massive mortality. Thus, vibriosis-resistant grouper is of great interest to the aquaculture industry, as it could reduce economic losses and facilitate aquaculture management. Studies have been conducted to comprehend the disease etiology, and a marker-assisted selective breeding scheme has been developed to reproduce grouper offspring with greater disease resistance. Disease resistance in grouper has been extensively studied at the molecular level through transcriptomics (Huang et al., 2011;Low et al., 2014aLow et al., , 2015aMu et al., 2010), proteomics (Low et al., 2014b(Low et al., , 2015b, and metabolomics approaches (Johnson & Brown, 2011;Karakach et al., 2009). In addition to extensive studies of several metabolites with antibacterial properties (Dee & Gradle, 2011;Desbois & Smith, 2010;Heath & Rock, 2004;Ouattara et al., 1997;Zheng et al., 2005), a recently conducted study has identified highly abundant metabolites, such as icosapentaenoic acid, eicosa-8,11,14-trienoic acid, and linoleic acid in brown-marbled grouper, Epinephelus fuscoguttatus, which has resisted Vibrio vulnificus infection (Nurdalila, Mayalvanan & Baharum, 2019). Study by Zhao et al. (2015) on the global metabolic response of tilapia against streptococcosis showed the involvement of specific metabolites in fish defense system against bacterial infection where they have identified l-proline contributes to the increased survival rate. While study on the V. vulnificus resistance in grouper has identified several fatty acids that were highly abundant during infection (Nurdalila, Mayalvanan & Baharum, 2019), nine fatty acids were selected to evaluate their potential molecular interaction with quorum sensing receptor, LuxP through molecular docking and simulation analysis. The findings from this experiment would support the importance of specific metabolites in fish defense against bacterial infection.
Quorum sensing is a bacterial cell-to-cell communication that allows a population of pathogenic bacteria to coordinate their gene expression, achieving collective behavior to evade the host immune system, to express toxic virulence factors and form antibiotic-resistant biofilms (Annous, Fratamico & Smith, 2009;Kim, Lee & Choi, 2012;Kim et al., 2003;Liang et al., 2007;Liu et al., 2013). Quorum sensing system regulation in V. harveyi is well characterized (Chen et al., 2002;Liu et al., 2013;Zhu et al., 2012). In V. harveyi, the quorum sensing system is activated by a boron-containing signaling molecule, furanosyl borate diester, that binds to the periplasmic receptor protein LuxP (Chen et al., 2002), which then forms a complex with LuxQ, a membrane protein.
The activated LuxPQ complex dephosphorylates the downstream proteins LuxU and LuxO and subsequently activates the transcription of the luciferase targeted genes (Liu et al., 2013;Schauder et al., 2001;Zhu et al., 2002Zhu et al., , 2012, leading to the expression of bioluminescence, biofilm formation and siderophore and metalloprotease production (Guo et al., 2013). In the attempt to suppress the expression of the virulence genes regulated by the quorum sensing system (Amara et al., 2010;Hentzer et al., 2002;Kalia, 2013;Guo et al., 2013;Rasmussen & Givskov, 2006;Saedi et al., 2012), studies have identified a range of secondary metabolites from microorganism (Morohoshi et al., 2008;Park et al., 2005;Rasmussen et al., 2005;Swem et al., 2009;Teasdale et al., 2009) and plant species (Packiavathy et al., 2012) that possess potent inhibitory properties against quorum sensing. In addition, it has recently been reviewed that fatty acid, cis-2-decenoic acid possesses inhibitory properties against biofilm production that was regulated by quorum sensing (Marques, Davies & Sauer, 2015). Antibiofilm activity has also been demonstrated by mono-unsaturated chain fatty acids, palmitoleic, and myristoleic acids (Nicol et al., 2018). Even though mammalian enzymes that hydrolyze the quorum sensing signaling molecules have been identified and characterized, but these enzymes were reported to be absent in model fish species (Morohoshi et al., 2008;Yang et al., 2005). Therefore, it was proposed that fish species might possess metabolic strategies to suppress bacterial pathogenicity and prevent infections (Nurdalila, Mayalvanan & Baharum, 2019;(Zhao et al., 2014(Zhao et al., , 2015. The potential of previously identified fatty acids (Nurdalila, Mayalvanan & Baharum, 2019) to interact with autoinducer-2 (AI-2) quorum sensing receptor were assessed by docking and molecular dynamic simulations in this study. Through this computational approach, the list of fatty acids was screened by AutoDockVina (Trott & Olson, 2010) to select for the best docking position with the lowest binding affinity score toward LuxP receptor protein. The selected protein-ligand complexes were then examined by molecular dynamic simulations using GROMACS (Van der Spoel et al., 2013;Hess et al., 2008) to evaluate the stability of the structure.

Preparation of protein receptors, ligands, and reference compounds
Autoinducer 2-binding periplasmic LuxP protein sequence of V. vulnificus was obtained from the NCBI Protein Databank (accession number: WP_011152474). I-TASSER (http://zhanglab.ccmb.med.umich.edu/I-TASSER/) was used to predict the 3D structure of this protein. The top-ranked model was selected for further analysis and named as V. vulnificus LuxP (vvLuxP) throughout this paper. The crystallographic co-ordinates for LuxP V. harveyi (PDB ID: 1JX6), named as V. harveyi LuxP (vhLuxP) throughout this paper, were selected and included in the study to serve as a control model for vvLuxP. The model was prepared by removing the endogenous ligand (furanosyl borate diester) and water molecules using AutoDockTools-1.5.6, and the hydrogen atoms were added to the structure. The fatty acid structure coordinates were obtained from Ligand Expo in the Protein Data Bank (http://ligand-expo.rcsb.org/). Three-dimensional structures of the reference compounds were generated using an online open access tool (https://web.chemdoodle.com). These reference compounds were randomly selected based on the findings that they can inhibit quorum sensing in V. harveyi (Zhu et al., 2012). Molecules were converted into the PDBQT file format prior to molecular docking.

Docking and molecular dynamic simulations of protein-ligand complexes
A grid box of size 30 Â 30 Â 30 A 3 was generated to contain the LuxP protein, with the LuxP binding pocket set as a centroid using AutoDockTools-1.5.6. Fatty acids and reference compounds were docked in LuxP receptor binding pocket using AutoDockVina 1.1.2 (Trott & Olson, 2010). The fatty acids with closest approximate affinity toward LuxP receptor compared to the reference compounds were shortlisted and further analyzed by molecular dynamic simulations using GROMACS 4.6.5 ( Van der Spoel et al., 2013;Hess et al., 2008). Protonation and structure minimization was performed using the GROMOS 54A7 force field, where hydrogens were added for optimal hydrogen bond network by default. Topology files for molecules were generated using PRODRG server (http://davapc1.bioch.dundee.ac.uk/cgi-bin/prodrg). LuxP-ligand complexes were solvated and fully immersed in the center of a cubic box prior to electrostatic energy calculation. A default 3-point model of SPC water model of GROMACS was used to solvate the box. The number of water molecules adopted to solvate the complex were as follow: vhLuxP_4R8 (23,331 molecules); vhLuxP_EPA (23,333 molecules); vhLuxP_Lax (23,331 molecules); vhLuxP_C19 (23,336 molecules); vhLuxP_C31 (23,338 molecules); vvLuxP_4R8 (22,425 molecules); vvLuxP_EPA (22,427 molecules); vvLuxP_LAX (22,430 molecules); vvLuxP_LAX#2 (22,417 molecules); vvLuxP_C19 (22,429 molecules); and vvLuxP_C31 (22,428 molecules). Electrostatic energy was calculated using gromacs preprocessor, and the system was neutralized by adding in accordance Na + ions or Clions to create zero charged system, and subsequent energy minimization was performed. The energy minimization was performed using the steepest descent minimization of 5,000 steps (maximum number of minimization steps to perform). Energy minimization was stopped when the maximum force was less than 1.0 kJ/mol. The system was further equilibrated for 50 ps at constant volume and a temperature of 293 K. The molecular dynamic simulations were run for 10,000 ps for each protein-ligand complex, where the coordinates were saved every two ps interval. LINear Constraint Solver, LINCS algorithm was applied to constraint all bonds, including heavy atom-H bonds during the molecular dynamics (MD) simulations. Long-range electrostatic interactions were treated with the adoption of Particle Mesh Ewald method (Chen, Wang & Zhu, 2014), and the cut-off distances for the long range electrostatic and Van der Waals interactions were set at 1.0 nm. Lastly, the trajectories were saved for further analysis using the Xmgrace and UCSF Chimera (Pettersen et al., 2004).

Re-scoring of protein-ligand complexes using interaction energy and MM-PBSA approach
The interaction energy of the protein-ligand complexes was calculated using molecular mechanics Poisson-Boltzmann surface area (MMPBSA) (Kumari et al., 2014) tool in Gromacs. The binding energy components were calculated separately as the MM, PB, and SA energy. The binding free energy of each complex was calculated from 20 snapshots at time intervals of 0.5 from the 10 ns MD production run. In general, the binding free energy (DG bind ) of fatty acids/reference compounds at LuxP receptor protein was calculated as follows: where G complex , G protein , G ligand are the free energies of the complex, LuxP receptor protein, and fatty acid/reference compound, respectively. The free energy (G) of each state was calculated as follows: where E MM is the molecular mechanical energy, G PB and G SA are the polar and nonpolar terms of the free energy, and TS is the entropic contribution of the solute. The solvent accessible surface area (SASA) and solvent accessible volume (SAV) models were used to calculate its contribution to the binding free energy of the complexes.

Comparative modeling of LuxP protein
The V. vulnificus LuxP protein sequence was submitted to the I-TASSER online server for homology structure prediction and five top final models with C-scores ranging from -0.22 to -3.16 were generated. Higher C-score values correlate with higher confidence levels in the 3D model, and, the model with C-score -0.22 was selected for subsequent docking and MD simulation. The top-ranked template identified by LOMETS with the highest normalized Z-score of 6.94 is the crystal structure of V. harveyi LuxP (PDB: 1JX6), which is one of the closest species to V. vulnificus. conformational changes upon ligand binding (Fig. 2). Both protein models were docked with fatty acids and the reference compounds to analyze the stability of the complexes via MD simulations.

Docking and binding pose analysis of fatty acids and reference compounds
Fatty acids that were found to be differentially abundant in E. fuscoguttatus that demonstrated their tolerance to the experimental infection by V. vulnificus were docked in both vhLuxP and vvLuxP models using AutoDockVina 1.1.2. The results showed that 4-oxodocosahexaenoic acid (4R8), 5,8,11,14,17-icosapentaenoic acid (EPA), and eicosa-8,11,14-trienoic acid (LAX) have significant binding affinity toward both vhLuxP and vvLuxP model. The 2D conformation of the fatty acids and the reference compounds with their binding affinity scores are given in Table 1. In the vvLuxP model, LAX was found to present a second binding pose (LAX#2) posterior to the vvLuxP receptor binding pocket, as shown in Fig. 3, with binding affinity score -5.5 kcal/mol (Table 1). This binding pose of LAX#2 was included in the MD simulation to evaluate the stability of the  Eicosa-8,11,14-trienoic acid (LAX) -7.4 -5.9 -5.5 Icosanoic acid (DCR) -6.6 -4.4 Linoleic acid (EIC) -7.2 -4.8 9-Octadecenoic acid (ELA) -6.8 -4.6 Alpha-linolenic acid (LNL) -7.1 -4.9 Oleic acid (OLI) -6.9 -4.9 Palmitoleic acid (PAM) -6.7 -5.4 complex. A significant cut-off binding affinity score of -5.5 kcal/mol was selected; complexes with not significant higher binding affinity scores were excluded from further molecular dynamic simulations analysis. Only parallel models of vhLuxP were selected to serve as control models. Reference compounds showed binding affinity score ranging from -7.5 to -8.0 kcal/mol against both LuxP models. eicosapentaenoic acid (EPA)  presented the highest affinity toward both LuxP models, with affinity scores of -7.8 and -6.0 kcal/mol. To visualize the interaction pattern, the UCSF Chimera molecular visualization tool was used to generate graphical representations of the hydrogen bonds (H_bonds) formed between the amino acid residues. Details of the hydrogen bonds formed between the amino acid residues with the fatty acids were reported in Table 2. The EPA molecule was found to interact with three key binding residues in vhLuxP receptor protein, as shown in Fig. 4. All these complexes were further assessed by MD simulation for structural behavior and flexibility analysis.

Molecular dynamic simulations
The structural behavior and flexibility of vhLuxP and vvLuxP docked with fatty acids and reference compounds were assessed by 10 ns of MD simulation using Gromacs 4.6.5 for each complex. Preliminary simulation for 50 ns was assessed and the stabilized protein backbone was identified after 10 ns of simulation. Due to the limitation in computing power and the preliminary simulation result, the structural behavior and flexibility of LuxP complexes were assessed for 10 ns of MD simulation. The root mean square deviation values of both vhLuxP and vvLuxP were calculated against the initial structure in the protein-ligand complexes and plotted using the 3-D plotting tool Xmgrace to compare the protein backbone stability. The backbones of the vvLuxP-ligand complexes showed significant fluctuation compared to the vhLuxP-ligand complexes (Fig. 5), implying that the binding of fatty acids and the reference compounds in vhLuxP is more stable and does not affect the protein backbone stability. The LuxP-ligand complexes were snapshot at one ns intervals from the 10 ns MD production trajectory and were superimposed to assess the ligand binding stability as in Fig. 6, observed the stability of the ligand binding position in each of the complexes. Meanwhile, the stability of the LuxP-ligand complexes was examined by calculating the residual mobility. The Root Mean Square Fluctuation of the trajectory from the MD simulation for each complex was calculated, and the protein residual fluctuations in LuxP-ligand complexes are minimal for both LuxP models (Fig. 7).

Re-scoring of interaction and binding free energies
The complex stability is further assessed by calculating the binding free energy using the g_mmpbsa tool (Kumari et al., 2014). Polar and non-polar energy terms were calculated for each complex. Both vhLuxP and vvLuxP hindered the binding of all fatty acids and reference compounds in terms of the polar solvation energy, which was recorded between 108.944 (lowest polar solvation energy) and 494.049 kJ/mol (highest polar salvation energy) (Table 3). Various non-polar energy terms (Van der Waals, VdW; SASA; and SAV) are favorable for all fatty acids and the reference compounds that bind in both LuxP models. Although both reference compounds have the least hydrogen bond interaction in vhLuxP and vvLuxP (Fig. 8), Compound C19 has the lowest binding free energy recorded, which was mainly contributed by low electrostatic energy in both vhLuxP (496.516 kJ/mol) and vvLuxP (-336.767 kJ/mol)and stabilized the binding of C19 in LuxP.

DISCUSSION
It has been reported that LuxP undergoes conformational changes upon (AI-2) binding to form the AI-2-LuxPQ complex (Stock, 2006;Zhu et al., 2012). Template-model based alignment (Fig. 2) has revealed the potential helices and beta-sheets that might be involved in the major conformational changes upon ligand binding, as shown in Fig. 1C. Crystal structure of LuxP without endogenous ligand is not available due to the fact that in order to crystallize the target protein, the compositional and conformational stability of the protein are prerequisite (Deller, Kong & Rupp, 2016). The binding of endogenous ligand, furanosyl-borate diester (AI-2) as observed in vhLuxP (crystal structure) involving the H_bond interactions with three key residues, which are ASN159, ARG215, and ARG310. Although the binding affinity score of LAX was lower than the other molecules, it is still considered as a potential candidate to interact with LuxP because it forms H_bond with the key residue of ARG310 in vvLuxP. It also has the highest number of H_bonds with the most key residues in the binding pocket of vhLuxP (Table 2). Surprisingly, reference compound 19 shows no H_bond interaction in the binding pocket despite having been proven to inhibit AI-2 quorum sensing (Zhu et al., 2012). The structural behavior and flexibility of vhLuxP and vvLuxP were assessed and the results were not as anticipated as the protein backbone of vvLuxP is expected to be more flexible since it is known to undergo major conformational changes upon ligand binding (Stock, 2006;Zhu et al., 2012). This may be due to the stability of ligands that bind inside the binding pocket that is crucial to determine the efficiency of the ligand in inhibiting LuxP protein.
However, the relatively short MD simulation of 10 ns could be another limiting factor to study and conclude the protein flexibility. Snapshots on the LuxP-ligand complexes were taken at one ns intervals from the 10 ns MD production trajectory and were superimposed to assess the ligand binding stability (Fig. 6). It showed that the positions of the ligands bound in vhLuxP are more confined and stable compared to vvLuxP. The open binding pocket of vvLuxP possesses higher accessible volume, and therefore the interaction of the ligand inside the binding pocket must be stronger and more specific in order to stabilize. When the stability of the LuxP-ligand complexes was examined, no significant fluctuation on the residues involved in the interaction was recorded in both vvLuxP and vhLuxP model. The simulation results showed similar dynamics for vvLuxP model and the control vhLuxP model, further strengthening the proposed open conformation of LuxP which was depicted in Fig. 1B to represent the LuxP apo structure. Although the average number of H_bonds varies significantly among the different complexes, the binding of these molecules in the protein binding site is relatively stable, except for the vvLuxP-LAX#2 complex, in which a higher level of fluctuation in the ligand binding pose is observable (Fig. 6). Nevertheless, less H_bond interaction is observed in both reference compounds, implying that the binding stability is contributed by the electrostatic forces. The continuous Notes: Each value represents the average value calculated from 20 snapshots at 0.5 ns intervals of the 10 ns MD production run. SASA, solvent accessible surface area; SAV, solvent accessible volume; WCA, Weeks-Chandler-Andersen.
contribution of H_bond interactions in the binding pose analysis suggests that these fatty acids possess potential stable interaction with the LuxP protein.
The binding free energies calculated using the MMPBSA approach indicated that among the three fatty acids, 4R8 displayed the lowest binding free energy (-29.144 kJ/mol) in vvLuxP. The lower binding free energy suggests that the vvLuxP-4R8 complex is more stable than the vvLuxP-EPA and vvLuxP-LAX complexes, which have higher binding free energies of 31.739 and 12.389 kJ/mol, respectively. The results obtained from this study suggest that 4R8 is a potential candidate possesses molecular interaction with quorum sensing receptor LuxP protein. Nevertheless, further in vitro/in vivo experiments are recommended to evaluate the biological effect of 4R8 interaction with the quorum sensing receptor LuxP of quorum sensing signaling in Vibrio.

CONCLUSIONS
The MD of the LuxP-ligand complexes, as shown by MD simulation, identified 4-oxodocosahexaenoic acid (4R8) as a potential candidate to interact with the LuxP receptor protein. The binding of 4R8 in both vhLuxP and vvLuxP displayed a stable protein backbone conformation during the MD simulation and a stable ligand binding position with the highest hydrogen bond interactions within the 10 ns MD simulation. The computational analysis results suggest a molecular interaction of the fatty acid 4R8 with the quorum sensing receptor LuxP. This knowledge is highly valuable for further in vivo/in vitro analysis to evaluate the effect of this molecular interaction in the biological signaling of Vibrio quorum sensing system.