Interaction Analysis of Odorant-Binding Protein 12 from Sirex noctilio and Volatiles from Host Plants and Symbiotic Fungi Based on Molecule Dynamics Simulation

As a quarantine pest of conifer, Sirex noctilio has caused widespread harm around the world. It is expected that the molecular mechanism of protein–ligand binding can be elucidated to carry out the pest control. Through studies of SnocOBP12–ligand hydrophobic binding and dynamics and responsible amino acid residues identification, we got some promising results. SnocOBP12 had a general and excellent affinity for host plant volatiles, and may be a key protein for S. noctilio to find host plants. Among the many odor molecules that are bound to SnocOBP12, (−)-α-cedrene and (E)β-farnesene from host plants and (−)-globuol from the symbiotic fungi of Sirex noctilio stood out and formed highly stable complexes with SnocOBP12. By the molecular mechanics Poisson–Boltzmann surface area (MM-PBSA) method, the calculated free binding energy of the three complexes was −30.572 ± 0.101 kcal/mol, −28.349 ± 0.119 kcal/mol and −25.244 ± 0.152 kcal/mol, respectively. It was found that the van der Waals energy contributed to the stability of the complexes. Some key amino acid residues were also found: LEU74 and TYR109 were very important for SnocOBP12 to stably bind (−)-α-cedrene, while for (E)-β-farnesene, ILE6, MET10, and LEU74 were very important for the stable binding system. We discovered three potential ligands and analyzed the interaction pattern of the protein with them, this paper provides a favorable molecular basis for optimizing the attractant formulation. Investigation of the binding characteristics in the olfactory system at the molecular level is helpful to understand the behavior of S. noctilio and develop new methods for more effective and environmentally friendly pest control.


Introduction
Insects use their olfactory system to sense semiochemicals in the environment and subsequently adjust behaviors such as mating, host positioning, food selection, inter-species recognition, intraspecies communication and habitat selection [1,2]. As selective signal filters, odorant-binding proteins (OBPs) can transport these hydrophobic semiochemicals penetrating into the antennal olfactory sensilla, a tiny hairy porous epidermal structure, through hemolymph, thereby activating specific olfactory receptors (ORs) or ionotropic receptors (IRs) expressed by olfactory sensory neurons (OSNs), thereby initiating signaling transduction [3]. In this process, the selective binding of OBPs to odor molecules lays foundations for selective identification of odors and plays a crucial role in the peripheral olfactory system. OBPs of insects, which was first reported in the antennae of Antheraea polyphemus (Lepidoptera: Saturniidae) [4], is a water-soluble protein composed of a single polypeptide chain, around 140 residues [5]. Classical OBPs usually rely on three interlocking disulfide bonds paired from six conserved cysteines to maintain the stability of protein tertiary structure formed by alpha-helical domains and have a hydrophobic cavity inside used for ligand binding [6][7][8]. Along with these, there are also three categories of OBPs classified according to the number of conserved cysteines in their primary protein sequences and other sequence characteristics, i.e., Plus-C OBPs (more than six Cys residues and a highly conserved Pro), Minus-C OBPs (less than six Cys residues) and Atypical OBPs (greater than or equal to six Cys residues and a long C-terminus) [9]. Benefiting from the advent of nextgeneration sequencing (NGS) technologies, more and more insect odorant-binding proteins have been identified and their functions have been studied in recent years. An example is that 28 classic OBPs and five Plus-C OBPs were first identified in Tessaratoma papillosa antennae based on NGS data, which provided a comprehensive resource for in-depth analysis on olfactory proteins functions of stink bugs [10].
Expression of this class of genes has also been found in other organs. Meanwhile, gene expression differences have been observed in the different tissues of conspecific insects. This means that the function of OBPs is not limited to transporting odor molecules [11]. Activation of the olfactory receptor DmelOR67d requires a complex comprised of LUSH (DmelOBP76a) and cis-vaccenyl acetate (cVA), the male-specific pheromone of Drosophila melanogaster, rather than cVA alone [12]. Involved in the uptake of nutrients, DmelOBP19b, highly expressed in D. melanogaster taste sensilla, is necessary for detection of L-phenylalanine and L-glutamine that it cannot synthesize [13]. The regulation of release of pheromones may be contributed by OBP10 in seminal fluid of Helicoverpa armigera and H. assulta [14]. Apart from the above, OBPs may be also involved in the regulation of insecticides resistance. The expression of PxylOBP13 increases in Plutella xylostella after treatment with permethrin at low heterogeneous concentrations [15].
The function of OBPs is currently considered to be diverse, especially when they have no antenna-specific expression. However, most OBPs are mainly expressed in the antennae, which is the essential peripheral sensory structure of insects; most relevant studies confirmed that OBPs were more inclined to reflect its odorant-binding function [16]. Furthermore, how OBPs bind with various ligands has gradually drawn widespread attention. To date, a large number of OBPs have been discovered and their binding affinity to odorant ligands has been investigated via fluorescence competitive binding experiments [17]. However, the mechanisms underlying the interactions between OBPs and ligands are still unclear. To cope with the lack of crystallographic models of OBPs, computational simulation has become another important method in deeper studies of protein-ligand interactions after experiment and theory [18][19][20]. For example, using simulation of molecular docking and site-directed mutagenesis, what was revealed was nine key residues for the ligand binding specificity and their effect on binding activity in the binding modes of multiple chemicals to SlitOBP1 in Spodoptera litura [21]. Two dominant sites of Pheromone-Binding Protein 3 in Plutella xylostella were also found to distinguish (Z)-11-hexadecenyl acetate from (Z)-11-hexadecenal by molecular dynamics (MD) simulation and computational alanine scanning (CAS) [22]. In S. noctilio, computational interaction of SnocOBP7 and female sex pheromones and symbiotic fungal volatiles were analyzed. The three best binding ligands were found to be female sex pheromones ((Z)-7-heptacosene and (Z)-7-nonacosene) and symbiotic fungal volatiles ((−)-globulol). Five residues played key roles in the binding of each female sex pheromone to SnocOBP7, whereas two residues play key roles in (−)-globulol binding [23].
As a high-risk forestry invasive alien species, Sirex noctilio (Hymenoptera: Siricidae), originating in parts of Eurasia and North Africa, has been transmitted to and colonized in multiple regions of the world, owing to human factors such as frequent international trade for almost a century [24][25][26]. The lack of competing organisms and natural enemies in new areas leads to heavy ecological and economic losses caused by this wood wasp [27]. It is S. noctilio, the wood-boring insect, along with its venom and symbiotic fungi, that harm the pine-dominated hosts together, aggravating the host tree weakening and accelerating death [28,29]. In northeast China, this non-native insect was first observed in 2013 and has shown a tendency to spread [30]. Thus, the development of effective means of monitoring and controlling pest populations in China is of tremendous importance to prevent the spread of S. noctilio and reduce the damage they cause. A control strategy based on semiochemicals provides a new direction for the pest management. As an environmentally friendly method, attractants are broadly utilized to trap and kill adults.
Without a doubt, pheromones are pivotal for chemical communication in insects. To monitor insect populations, pheromones are widely used. At present, some pheromones of S. noctilio have been identified [31][32][33]. There are (Z)-dec-3-en-1-ol, (Z)-dec-4-en-1-ol, (2E,4E)-deca-2,4-dienal, (Z)-oct-3-en-1-ol and (Z)-dodec-3-en-1-ol from the male and (Z)-7heptacosene, (Z)-7-nonacosene and (Z)-9-nonacosene from the female, which can provoke strong responses from S. noctilio in gas chromatography-electroantennographic detection (GC-EAD) assays, and Y-tube olfactometer and wind tunnel assays. Besides, kairomone (plant volatiles) lure traps were found to be the most effective method for trapping S. noctilio individuals in areas with large populations, and (+)-limonene, (−)-limonene, α-pinene, β-pinene, 3-carene and camphene were attractive host plant volatiles in early reports [34]. However, it is discovered that the effective trapping in forests is primarily dependent on the attractants of plant origin, such as the mixture of α-pinene, β-pinene, 3-carene and camphene, but not pheromone lure cores [35]. At high population densities of S. noctilio, the trapping method using host plant volatiles as lures can achieve awesome results, but it is unfavorable and unsuitable at low densities. Fortunately, the addition of volatiles released from fungi symbiotic with S. noctilio may improve the effectiveness of attractants at low population densities [36,37]. Accordingly, the three promising chemical cues pheromones, host plant volatiles and symbiotic fungi volatiles are chosen as the focus of this study on the basis of the above findings.
Two major factors affecting olfactory binding protein function are expression abundance and binding affinity. In Sirex noctilio, 16 OBPs were identified from transcriptome analysis by our team, one of which is our present subject, SnocOBP12, highly expressed in antennae based on quantitative polymerase chain reaction (qPCR) analysis of the tissue [38]. It suggests that SnocOBP12 may be involved in the process of chemoreception in insects. However, its specific role in the olfactory system is uncertain and the molecular interaction mechanism between SnocOBP12 and ligands is still not shown.
Given the high expression of SnocOBP12, we focus primarily on computational methods to investigate its ability to bind the three semiochemicals, which is essential to explore its function. In this paper, our study was divided into the following two parts: (1) The biological characteristics of SnocOBP12 were obtained through the systematic analysis of bioinformatics. (2) There is a three-dimensional model of SnocOBP12 predicted by homology modeling in this study. Thereafter, with computational simulation, its binding preference was analyzed and the critical amino acid residues for ligand binding was further identified. This is beneficial to better explain the molecular mechanism of SnocOBP12 action.
It is expected that we can find the potential lure molecules with high specificity via the research on the SnocOBP12-ligand binding mode, to optimize the attractant formula for better results.

Verification of the SnocOBP12 Sequence
There were Sirex noctilio eggs or larvae in the wood segment felled in Durbert Mongolian Autonomous County, Daqing City, Heilongjiang Province. After natural emergence, adult antennae were collected, immediately frozen in liquid nitrogen, and stored at −80 • C. Total RNA was extracted according to the instructions of Trizol and RNA extraction kits. After 1% agarose gel electrophoresis and ultra-trace UV spectrophotometer testing to meet the quality requirements, the first strand of cDNA was synthesized according to the reverse transcription kit of TaKaRa company, and stored at −20 • C for future use. To verify the SnocOBP12 sequence by cloning and sequencing, the specific primers were designed with Primer3plus (http://www.primer3plus.com/cgi-bin/dev/primer3plus.cgi; accessed on 12 April 2021) ( Table 1). PCR amplification was performed using cDNA as template. Pre-denaturation was conducted at 95 • C for 2 min, 95 • C for 20 s, 53 • C for 5 s; and 72 • C for 5 s; then, 20 cycles, at 72 • C, for 5 min. The amplified products were detected by 1% agarose gel electrophoresis and recovered. The target gene SnocOBP12 obtained by double digestion was ligated with the expression vector pET28A with T4 ligase (NEB). To ensure the correctness of the sequence, the product was sequenced after completing the transformation of the OBP12-pET28A expression vector.

Homology Modeling and Molecular Docking
The SnocOBP12 sequences were aligned with the standard protein blast (blastp) in NCBI (https://blast.ncbi.nlm.nih.gov/Blast.cgi; accessed on 25 August 2021), and the protein sequences with better homology and scores were chosen as candidate templates for modeling. The 3D model of SnocOBP12 was calculated with MODELLER and the model quality assessment was completed with SAVES 6.0 (https://saves.mbi.ucla.edu/; accessed on 26 September 2021) [41]. Comparison of amino acid sequence identity between SnocOBP12 and template proteins were performed using ClustalW (www.genome.jp/tools/ clustalw/; accessed on 28 September 2021) and the result was depicted by ENDScript (https: //espript.ibcp.fr/ESPript/cgi-bin/ESPript.cgi; accessed on 28 September 2021). In this paper, we used Autodock Tool 1.5.6 to perform semi-flexible molecular docking between SnocOBP12 and 24 ligands selected from S. noctilio pheromones, host plant volatiles and symbiotic fungi volatiles [42]. The Grid Box center was located at (−28.147, −11.763, −16.449). The information of all ligands was searched in PubChem (https://pubchem.ncbi. nlm.nih.gov/; accessed on 21 May 2021) ( Table 2). Evaluating and sorting according to the reasonableness of the complex conformation and the free binding energy and other parameters, the conformation with a reasonable binding site and the lowest free binding energy was selected.

Molecular Dynamics Simulation of the SnocOBP12−Ligand Complexes
A 50 ns molecular dynamics (MD) simulation with a 2-fs time step was performed to analyze conformational changes of SnocOBP12-ligand complexes. The model obtained by molecular docking served as the initial conformation for the simulation. The production runs were carried out with the GROMACS package (version 2019.6) using the amber99sb-ildn force field for protein and the tip3p water model [43]. For ligands, they were parametrized using the AMBER force field (GAFF) 18, through use of the ACPYPE script [44]. The system energy was minimized with the conjugate gradient method, and the whole system was balanced at 298.15 K with a V-rescale thermostat and 1 bar with the Parrinello-Rahman barostat. Simulation experiments were repeated under the same conditions to examine reproducibility. In this paper, the typical result was selected for subsequent analysis.
During the 50 ns simulation, the long-range electrostatic interactions had been treated with the particle mesh Ewald (PME) method, and the covalent bonds with hydrogen atoms had been constrained with the Linear Constraint Solver (LINCS) algorithm. Meanwhile, the cutoff radius for non-bonded interactions was set to 10 Å in all stages of the simulation. The root-mean-square deviation (RMSD) was obtained as the important basis to measure the stability of the SnocOBP12-ligand system, while the root-mean-square fluctuation (RMSF) was used to determine the flexibility of a region of SnocOBP12. Dominant conformations were obtained to observe the interaction between SnocOBP12 and ligands by cluster analysis, and visualized with the Discovery Studio 2019 (DS, Accelrys Inc., San Diego, CA, USA).

Binding Free Energy Calculations and Per-Residue Free Energy Decomposition
For energy calculations, the relatively stable part of the 50 ns simulation time was intercepted. Calculating the free binding energy of SnocOBP12 with ligands through the molecular mechanics Poisson-Boltzmann surface area (MM-PBSA) method [45], it was well known how much energy each residue of SnocOBP12-ligand complexes contributes. In addition, the binding free energy was divided into polar and nonpolar solvation energy, van der Waals, and electrostatic energy.

Computational Alanine Scanning (CAS)
In this paper, the prediction of key amino acid residues was carried out with Discovery Studio 2019. The alanine replaced key amino acid residues with lower than −1.00 kcal/mol in per-residue free energy decomposition, and the resulting impact was measured by mutation energy (∆E mut ) (∆E mut < −0.5 kcal/mol, stabilizing; −0.5 kcal/mol ≥ ∆E mut ≤ 0.5 kcal/mol, neutral; ∆E mut > 0.5 kcal/mol, destabilizing). Combining the results of MM-PBSA and CAS, the key amino acid residues, it was evident which key amino acid residues affected the binding affinity of SnocOBP12 to different ligands.

Sequence Analysis and Homology Modeling
After a preliminary bioinformatics analysis, it is shown that SnocOBP12 is an extracellular protein (molecular weight: 13.66 ku) with a theoretical pI of 6.57. The full-length protein is composed of 140 amino acids including the signal peptide region of the N-terminal amino acids 1-22. As a classical odorant-binding protein, the stability of the 3D structure of SnocOBP12 comprised of six α-helices is maintained by three disulfide bonds formed by six conserved cysteines (CYS17-CYS48, CYS44-CYS99, CYS89-CYS108) ( Figure 1C). The target hydrophobic binding cavity is mainly composed of α1 helix, α4 helix and α5 helix. Additionally, a few amino acid residues in the α3 helix and α6 helix were also involved in the formation of the binding cavity. To construct the valid model, the crystal structure of OBP 4 in Chrysopa pallens (PDB ID: 6JPM, chain A) ( Figure 1B) was selected as the most satisfactory template, which had a sufficient homology for SnocOBP12 (sequence identity > 30.0%) ( Figure 1C). After a simple optimization, the final model of SnocOBP12 ( Figure 1A) satisfying the rationality verification was obtained for in-depth experiments. In the Ramachandran plot ( Figure S1A), 92.7% (>90%) of non-glycine and non-proline residues were located in the most favored regions, implying good stereochemical quality of the constructed SnocOBP12 model. Meanwhile, 89.17% (>80%) of amino acids had a proper score (≥0.2) for compatibility between the three-dimensional structure and primary structure, and the overall quality factor of ERRAT was 94.643 ( Figure S1B,C). The above exactly proved the reliability of the model.

Molecular Docking and Binding Affinities for SnocOBP12 with Ligands
The binding energy between the receptor and the ligand reflects the binding situation between the two. The smaller the binding energy is, the more stable the receptor-ligand complex is [46]. In this study, the binding energy of SnocOBP12 to most ligands was concentrated between −1.0 kcal/mol and −1.4 kcal/mol, and for symbiotic fungi volatiles, the binding energy was uneven ( Figure 2). The binding energy of (E)-hex-3-enyl acetate was −1.01 kcal/mol, while (−)-globulol was −1.45 kcal/mol. Specifically, (Z)-oct-3-en-1-ol had the highest binding energy of −0.93 kcal/mol, while (−)-α-cedrene had the lowest binding energy of −1.73 kcal/mol. In the docking simulation of α-pinene and β-pinene with SnocOBP12, their levoisomers were always more advantageous.
Considering the source of the ligand, the docking result ( Figure 2) suggested that S. noctilio pheromones generally had higher binding energy than the other two types of substances. By contrast, host plant volatiles and some symbiotic fungi volatiles had a stronger affinity to SnocOBP12.

MD and Conformational Stability of SnocOBP12-Ligand Complexes
With binding energy below −1.2 kcal/mol, the complexes of host plant volatiles and some symbiotic fungi volatiles with SnocOBP12 were selected for 50 ns molecular dynamics simulations. From the RMSD diagram (Figure 3), it can be seen that most of the proteinligand systems reached a relative equilibrium at around 10 ns, while the three complexes, C H -7, C H -9, and C F -1, reached equilibrium at around 20 ns. The average RMSD value of the C H -1 complex was about 3.9 Å, and the other 15 complexes were concentrated between 2.6 Å and 3.5 Å. Verifying the docking result, the 16 complex systems all showed obvious stability during the simulation process, and there was the excellent interaction between the ligands and SnocOBP12.
The local motility properties of amino acids residues when SnocOBP12 complexed with ligands were further determined by the root-mean-square fluctuation (RMSF). In the 50 ns simulation, the 16 complexes showed similar motion trends (Figure 4). Among them, the maximum RMSF value (6.193 Å) appeared on LYS92 with a high degree of freedom in the C F -3 complex. On the whole, there were three main areas with general violent fluctuations in the RMSF image, including the α4 helix, α5 helix and N-terminus, which were very close to the SnocOBP12-ligand binding site. In this region, the stability of SnocOBP12-ligand complex interactions was not high. There were also two relatively stable regions, residues 42-51 and 99-110. It was worth noting that the residues 99-110 in the α5 helix were still very close to the binding site. It was speculated that their existence was one of the key factors in maintaining binding stability.

Energy Calculation and Binding Modes Analysis
Based on the relatively stable orbital data intercepted during the 50 ns molecular dynamics simulation, the free binding energy (∆G bind ) of 16 protein-ligand complexes was calculated by the MM-PBSA method and subdivided into four parts-van der Waals energy (∆G vdw ), electrostatic energy (∆Gele), polar solvation energy (∆G PB ), and non-polar solvation energy (∆G SA )-for specific analysis ( Table 3). The binding free energy of most complexes was around −19 kcal/mol, while the values of C H -11 (−30.572 ± 0.101 kcal/mol), C H -12 (−28.349 ± 0.119 kcal/mol) and C F -2 (−25.244 ± 0.152 kcal/mol) were prominent in all subjects. The major contributors to the C H -11 complex were van der Waals energy and non-polar solvation (SASA) energy. For C H -12 and C F -2, van der Waals energy, electrostatic energy and SASA energy provided the main guarantee for the stability of the system. According to Table 3, van der Waals energy, which dominates the good binding, had the most prominent contribution to the binding of proteins and ligands among all the simulated objects. SASA energy was also beneficial to the combination, although its contribution was weak. Polar solvation energy was greater than 0 kcal/mol, which was the main obstacle to the binding of protein and ligand, especially the binding of SnocOBP12 and symbiotic fungi volatiles. Having two sides, the functionality of electrostatic energy may vary depending on the ligand. The values of electrostatic energy only in the C H -2 (0.001 ± 0.011 kcal/mol), C H -3 (0.027 ± 0.013 kcal/mol), C H -6 (0.357 ± 0.015 kcal/mol), C H -7 (0.005 ± 0.017 kcal /mol) and C H -11 complexes (0.047 ± 0.008 kcal/mol) were positive, inhibiting binding, while all other complexes were negative, promoting binding.
After calculating the binding energy of the complex as a whole, the per-residue free energy decomposition of the 16 complexes was also completed by the MM-PBSA method, and the specific data were visualized ( Figure 5). The residues with binding energy below −1.0 kcal/mol were individually annotated in Figure 5, including ILE6, MET10, MET53, ILE70, ILE73, LEU74 and TYR109. It can be seen that TYR109 had the highest frequency, followed by LEU74 and MET53. Combined with the local movement of amino acid residues, the RMSF values of ILE6, MET10 and TYR109 fluctuated around 1.0 Å, while the RMSF value of MET53 was lower. They were all located in the trough region of the RMSF curve, which was positive for stably binding the ligand. The region where LEU74 was located had large fluctuations during the binding process, which might cause the interaction between LEU74 and ligands to be unstable.  After having an overall understanding of the binding situation of the two types of ligands, host plant volatiles and symbiotic fungi volatiles, the three systems with the lower free binding energy were selected for the following canonical conformational analysis. The complex conformations were clustered by constraining the range of motion of the three complexes, and the corresponding conformations were taken as canonical conformations. Considering the corresponding energy calculation results and the specific positions of the residues in the protein, the specific states and interactions of the three highly active ligands when they bound to SnocOBP12 were highlighted ( Figure 6). As shown in Figure 6, it was LEU74 and TYR109 that bound (−)-α-Cedrene tightly to SnocOBP12. ILE6, MET10, MET53 and LEU74 bound (E)-β-farnesene tightly and MET53, ILE70 and ILE73 bound (−)-globuol tightly to the protein. Among them, TYR109 belonged to polar amino acids, and the rest are non-polar amino acids. These hydrophobic residues tightly controlled the ligands in the binding cavity. At the same time, these amino acid residues also had different interaction forms with ligands (Figure 7), van der Waals and pi-alkyl being the most common. There are two special points worth noting. The first is that, due to the existence of TYR109, there was a hyperconjugation effect in the C H -11 complex, which favored the stability of this system. The second is that MET53 and the hydroxyl of (−)-globuol formed a conventional hydrogen bond to help stabilize the C F -2 complex system.

Computational Alanine Scanning
Based on the per-residue decomposition results, the key residues in the dominant system C H -11, C H -12 and C F -2 were targeted to perform the computational alanine scan, and the resulting mutation energy was calculated to measure the impact (Table 4). The results showed that the binding of (−)-α-cedrene to SnocOBP12 was destabilizing after the mutation of LEU74 and TYR109 in the C H -11 complex to Ala, and for C H -11, the mutation of ILE6, MET10 and LEU74 destabilized the SnocOBP12-(E)-β-farnesene complex. Unusually, the mutation of MET53, ILE70, and ILE73 in C F -2 complex had no obvious effect on the stability of the whole system. In other words, in the C H -11 complex, two amino acid residues were identified as key residues, and in the C H -12 complex, there were four key residues. These residues affected the binding ability of SnocOBP12 to ligands. However, none of the three residues in the C F -2 complex were key residues, although they contributed a lot of energy in the binding process. It was worth noting that LEU74 played an important role in both systems. In contrast, MET53 also had a prominent contribution in both systems, but its mutation did not affect binding stability in either binding system.

Discussion
Various odor molecules scattered in the air have a certain guiding effect on the life activities of Sirex noctilio. Under the situation that environmental problems are getting worse and the 3R problem (resistance, resurgence and residue) of pesticides is becoming more and more prominent, in-depth research on the mechanism of insect peripheral olfactory system recognition of odor molecules, and the development of drugs with certain tropism or targeting effect on insects will help to develop an efficient, green, safe, broad-spectrum pest control agent and play a positive role in the management of harmful insects.
There are two methods of biological control based on insect olfaction: trapping and poisoning. It is cost-effective and predictive to profile the selectivity and affinity of proteinligand binding, and computational methods have been to explore the potential mosquito OBP1 inhibitors from natural products in essential oils so as to facilitate further development of novel repellents [47]. In this study, we also used a similar method to identify three ligands with higher binding activity that have great potential to optimize existing attractant formulations to combine traps for more efficient killing. In addition, according to their binding characteristics and organic structure, specific insecticide molecules can be designed to inhibit or block the signal transmission of SnocOBP12 and interfere with normal life activities.
Stable complexes are more favorable for OBP transport of ligands [12]. Through homology modeling, molecular docking and molecular dynamics simulation, the binding characteristics of SnocOBP12 and 24 ligands have been analyzed. It suggests that Sno-cOBP12 is a classic odorant-binding protein mainly expressed in the antennae, which may be involved in the life activities of S. noctilio in finding and locating the host plants. It is expected that the existing attractant formulations can be optimized based on the discovery of potential ligands such as (−)-α-cedrene. The same method was used to analyze and reveal the molecular chemical mechanism of AlepPBP2 and AlepPBP3 in Athetis lepigone and the phoxim, and key amino acids have been identified, which would be of great help in the development of highly effective insecticides [48]. In contrast, when studying the mechanism of AgamOBP1-ligand interaction, the hierarchical virtual screening was used to screen potential ligands on a large scale and prioritize them [49]. For Anopheles gambiae, ligands are recognized by the AgamOBPs homodimer and AgamOBPs can bind to multiple ligands. The conformational changes when AgamOBP1 bound both N,N-diethyl-3-toluamide and 6-methyl-5-hepten-2-one were investigated by molecular dynamics [50]. According to the existing research, both SnocOBP12 and SnocOBP7 [23] can bind (−)-globuol well. The interaction between SnocOBPs remains to be explored.
OBPs have different binding affinities with different odor molecules, although they have a broad spectrum of ligand binding [51]. The ligand preference also exists in Sno-cOBP12. The protein prefers host plant volatiles, followed by symbiotic fungi volatiles, and the tendency towards pheromones is the weakest. The overall amount of volatiles released by the stressed host pine trees was significantly higher than that of the healthy pine trees. Plant-derived attractants consisting of stressed host plant volatiles were studied and used in some trapping experiments. For example, eight kinds of substances ((+)-α-Pinene (12.5%), (−)-α-Pinene (12.5%), (−)-β-Pinene (25.0%), 3-Carene (30.0%), Camphene (5.0%), β-Myrcene (10.0%), (−)-Limonene (2.5%), (+)-Limonene (2.5%)) were mixed as the main ingredients of lure core [52]. The forest trapping experiments using the volatiles from the ring-cut of Pinus sylvestris var. mongolica confirmed that α-pinene, β-pinene and 3-carene and others have strong attraction to S. noctilio. This is why they are always used in current attractant formulations. This is consistent with the results of molecular docking and molecular dynamics simulations in this paper. The binding energy of the complexes formed by the above volatiles combined with SnocOBP12 is relatively low [53].
What is more, three ligands with higher binding affinity were successfully screened out, two of which were host plant volatiles and the other was from symbiotic fungi. Currently, it is unclear what message (−)-α-cedrene conveys to S. noctilio, but it cannot be ignored that the affinity of SnocOBP12 for this odorant molecule is so pronounced. It may be a potential key to optimizing formulations and increasing lure rates. Thus, it is worth studying how (−)-α-cedrene will affect the physiological behavior of S. noctilio. For another host plant volatile, (E)-β-farnesene, there are few reports on its attraction to insects, but it has been reported that (E)-β-farnesene is an alarm signal in aphids [54]. It is speculated that (E)-β-farnesene might have an avoidance effect on S. noctilio. When the presence of (E)-β-farnesene is detected, they would choose to stay away. The volatile (−)-Globuol, of the symbiotic fungus, Amylostereum areolatum, may convey to S. noctilio that eggs have been laid on the plants and the hosts have been damaged, which is helpful for the female to locate the host plants with low growth vigor and lay eggs. Relevant experiments have confirmed that it has a significant attracting effect on the female. While for the male, the discovery of eucalyptol means that there may be female S. noctilio nearby, so as to carry out mating behavior [55].
From a microscopic point of view, the binding of ligands to proteins occurs in a hydrophobic binding cavity mainly surrounded by three α-helices, in which some residues, such as LEU74 and TYR109, contribute indispensably to complex formation strength.
Making an outstanding contribution, TYR109 is located on the α6 helix, which may be the key to locking the ligand into the binding cavity. This is also one of the reasons why (−)-α-cedrene binds best to SnocOBP12. In both C H -11 and C H -12 complexes, LEU74 was identified as a key amino acid residue and it is speculated that the other key residues, such as ILE6 and MET10, might change the active conformation of the protein to distinguish (−)-α-cedrene from (E)-β-farnesene. It is reported that the strong hydrophobic environment formed by multiple leucine residues represented by LEU79 lays a very favorable foundation for the stable binding of AmalOBP8 from Agrilus mali to (Z)-3-hexenyl hexanoate [56]. When analyzing protein-ligand interactions (Figure 7), multiple leucines can be found. For SnocOBP12, the leucine residues represented by LEU79 may be essential for the formation of a hydrophobic environment. In the binding process of SnocOBP12 and ligands, van der Waals force is the main driving force. Due to its existence, the ligands can stably bind to the protein and reach the receptor through the sensory lymph, and the olfactory nerve begins to interpret the information contained in the molecule. Similar to SnocOBP12, AmalOBP8 also had a larger ligand binding profile. However, it was hydrogen bonds rather than van der Waals forces that stabilized the complex between AmalOBP8 and geranyl formate [56].
Although the research on the binding characteristics of SnocOBP12 in this paper is fruitful, there is still room for improvement. First, this study only studies and analyzes the binding of proteins to parts of substance from S. noctilio pheromones, host plant volatiles and symbiotic fungi volatiles. The choice of ligands has certain limitations, and other ligands that can bind well to proteins may be available. Second, this study only analyzes the binding characteristics from the perspective of computational simulation to help guide and explain physiological experiments. The widely used amber99sb-ildn force field was chosen by our team for molecular dynamics simulations. However, different force fields (amber99sb force field and residue-specific force field) were used during simulations of AlepOBP2 binding to chlorpyrifos and phoxim. It was shown that the residue-specific force field can make calculations more efficient and reliable, providing favorable conditions for the subsequent systematic analysis of the binding of AlepOBPs to various molecules [57]. Besides, the fluorescence competitive binding assay was also widely used to analyze the binding capacity of proteins to ligands from an experimental point of view, such as BminOBP9 to 11 citrus volatiles [58]. These targeted changes and experimental verification that are beneficial to experiments are worth learning and exploring in investigations.
To sum up, as a classic OBP, SnocOBP12 highly expressed in the antennae takes on the responsibility of transporting molecules, and is closely involved in locating hosts and finding mates, etc. The relevant findings deserve further exploration. As a new strategy, the reverse chemical ecology approach was used to find an effective oviposition attractant acetaldehyde to Culex quinquefasciatus through the studies on the structure and function of olfactory proteins [59]. From the perspective of this idea, two research directions are proposed: (1) to optimize attractant formula, more potential ligands can be found with the help of hierarchical virtual screening or other methods, and their attractant effect can be tested in forests. (2) The amino acid residues in insect OBPs can be occupied by antagonists so that the protein cannot bind suitable ligands, thereby disrupting the normal behavior of insects to achieve control effects. In the future, fluorescence competition binding experiments, electrophysiology, gene silencing and other means are needed to verify our results, and to comprehensively and systematically understand the role of this protein in the life activities of S. noctilio in order to better control their population and reduce economic and ecological losses.

Conclusions
In this paper, the binding properties of SnocOBP12 were explored and analyzed by computational simulation methods including homology modeling, molecular docking, molecular dynamics and computational alanine scanning. We also determined the basic biological characteristics of SnocOBP12 through bioinformatics analysis to increase understanding of this protein. The discovery of (−)-α-cedrene, (E)-β-farnesene and (−)globuol will facilitate the optimization of existing attractants for enhanced control and monitoring of the S. noctilio population. Key amino acid residues discovered will also facilitate the antagonist development. As mentioned above, this study provides effective insights for explaining the binding properties and interaction mechanisms of OBPs and odorant molecules, and provides a theoretical basis and reliable reference for integrated pest management and biological control.
Supplementary Materials: The following are available online at https://www.mdpi.com/article/ 10.3390/agronomy12040861/s1, Figure S1: Quality evaluation of the OBP12 model.  The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.

Data Availability Statement:
The authors confirm that the data supporting the findings of this study are available within the article and its Supplementary Materials.