Computational Interaction Analysis of Sirex noctilio Odorant-Binding Protein (SnocOBP7) Combined with Female Sex Pheromones and Symbiotic Fungal Volatiles

: Sirex noctilio , a major forestry quarantine pest, has spread rapidly and caused serious harm. However, existing methods still need to be improved because its olfactory interaction mechanisms are poorly understood. In order to study the role of male-speciﬁc protein SnocOBP7 in the protein–ligand interactions, we selected it as the object of computational simulation and analysis. By docking it with 11 ligands and evaluating free binding energy decomposition, the three best binding ligands were found to be female sex pheromones (( Z )-7-heptacosene and ( Z )-7-nonacosene) and symbiotic fungal volatiles (( − )-globulol). Binding mode analysis and computational alanine scanning suggested that ﬁve residues play key roles in the binding of each female sex pheromone to SnocOBP7, whereas two residues play key roles in ( − )-globulol binding. Phe108 and Leu36 may be the crucial sites via which SnocOBP7 binds female sex pheromones, whereas Met40 may regulate the courtship behavior of males, and Leu61 may be related to mating and host ﬁnding. Our studies predicted the function of SnocOBP7 and found that the interaction between SnocOBP7 and pheromone is a complex process, and we successfully predicted its binding key amino-acid sites, providing a basis for the development of new prevention and control methods relying on female sex pheromones and symbiotic fungi.


Introduction
As one of the main mechanisms through which insects perceive the external environment, the olfactory system plays a very important role in foraging, defense, mating, reproduction, information exchange, and habitat selection [1]. Studying how odor molecules in the environment act on insect sensilla and induce insects to produce behavioral responses will illuminate the molecular interaction mechanisms underlying insect host identification, interspecific interaction, and intraspecific communication, allowing the formulation of corresponding management and protection strategies with advantages over currently available methods.
Signal transduction in the peripheral olfactory system of insects can be summarized as follows: odor molecules enter the antennal sensilla lymph through the stratum corneum, where they bind to odorant-binding proteins (OBPs) or chemosensory proteins (CSPs) [2]. Soluble proteins transport hydrophobic odorants to membrane-bound odorant receptors (ORs) or ionotropic receptors (IRs) located on the dendritic membrane of olfactory neurons through the sensory lymph, after which ORs or IRs convert chemical signals into electrical signals and transmit them to the central nervous system, thereby causing an olfactory response [3,4]. Finally, odor molecules are degraded by odor-degrading enzymes (ODEs) [5]. In general, after entering the organism through a series of olfactory proteins, the odor molecules are recognized and degraded, which is regarded as a vital ecological interaction process between organisms and the environment.
The combination of OBPs and odor molecules is the initial biochemical reaction through which insects specifically recognize external odorants [6]. A typical OBP has a stable and compact three-dimensional binding hydrophobic cavity, which is likely to be the key region involved in binding with external odor molecules [7,8]. However, OBPs are specific to the selection of environmental odor molecules based on their expression levels in different tissues. Correspondingly, through the regulation of the concentration or presence of certain odor molecules, organisms will affect the response of the biological olfactory system from the molecular perspective, and then manifest in the behavioral interactions. Thus, exploring how OBPs bind ligands in their binding cavity and revealing the binding interaction mechanisms underlying the effects of chemical volatiles on insects are important research subjects.
With the development of computer technology, computer simulations have been widely used in interdisciplinary studies of protein-ligand binding methods. Without getting the actual structure of the protein, homology modeling has become a commonly used method of predicting protein structure because proteins with high homology generally have similar three-dimensional structures [9]. For example, MvitOBP3 in Maruca vitrata was specifically predicted to bind to several host plant volatiles from legumes through homology modeling and molecular docking analysis [10]. However, because molecular docking only models transient and stable binding, the binding modes of protein and odorant molecules are not fully simulated. Therefore, molecular dynamics (MD) simulations must be included to determine whether the binding between a protein and odorant is stable, as well as to predict the major driving forces. For example, using a MD simulation, the seventh α-helix of Agrilus mali OBP8 quickly formed a loop structure upon binding with geranyl formate, which indicated that insect OBPs may be able to modify their secondary structures to increase the range of proteins with which they may bind [11]. In work aimed at deciphering the binding mechanism of Athetis lepigone GOBP2 to Chlorpyrifos and Phoxim, MD simulation was used to analyze the forces driving these interactions, revealing that the main driving forces were alkyl-π and hydrophobic interactions [12].
Sirex noctilio Fabricius (wood wasp; Hymenoptera: Siricidae: Sirex) is a major international forestry quarantine pest that is native to Eurasia and North Africa and has a wide range of hosts, mainly Pinus, and a few species in Picea, Abies, and Larix) [13][14][15][16]. A wide range of host tree species are conducive to the invasion and spread of wood wasps, and they increase the difficulty of its prevention and control, which has led to its continuous expansion worldwide for many years. In the last 100 years, with increased human activity and international trade, S. noctilio has invaded Oceania (New Zealand and Australia), South America (Uruguay, Argentina, Brazil, and Chile), North America (Canada and the United States), and Asia (China) [17][18][19][20][21][22][23], thus becoming a global invasive pest. Sun [24] predicted that every continent except Antarctica has areas suitable for S. noctilio. As a result, effective control of S. noctilio requires extreme vigilance, because this invasive pest can rapidly multiply and spread in new areas lacking competing species and natural enemies, causing great damage to host plants and serious economic losses [25]. Until now, S. noctilio in China has been distributed mainly in the northeast and parts of Inner Mongolia, which is consistent with the highly suitable area predicted by the CLIMEX model [26]. In countries where invasive S. noctilio colonization has occurred, parasitic wasps and nematodes have achieved good control effects [27]. However, the implementation of biological control is greatly affected by environmental factors. To achieve real-time monitoring of low-density populations, new methods of controlling S. noctilio from the perspective of chemical ecology through attractants are necessary.
Generally, odor molecules with potential utility as wood wasp attractants are divided into three categories: host plant volatiles, symbiotic fungal volatiles, and sex pheromones. One or more odor molecules are reasonably formulated into attractants to regulate the interactions of wood wasps on the environment and achieve the purpose of monitoring or prevention and control. Plant-derived attractants have long been used in forested environments for prevention and control of wood wasps, but their effects are relatively limited. Liu [28] used lure cores with four host plant volatile components to conduct forest experiments in Heilongjiang Province for four consecutive years, but this method was not found to be suitable for monitoring wood wasp populations. Although the use of bait wood to prevent and control wood wasps was in line with their habit of damaging weak host plants, and although it was also suitable for monitoring insect populations with low density, this method is costly and difficult to apply. After discovering the gap in trapping work, there is an urgent need to identify new attractants with improved specificity and optimize the attractant formulation [29]. Whether there is a new perspective that can effectively regulate the interactions between male adults and odor molecules and achieve effective prevention and control effects needs to be studied.
Pheromone control is becoming an important measure in comprehensive pest management due to its high efficiency, nontoxicity, and strong specificity. However, according to existing work [30,31], the use of S. noctilio pheromones alone was found to be an unsatisfactory method of forest trapping; thus, pheromones and plant volatiles were often used in combination. Previous studies have mostly focused on the development of male pheromone attractants. For example, Sarvary et al. [32] developed a lure consisting of three male pheromones ((Z)-3-decen-1-ol, (Z)-4-decen-1-ol, and (E,E)-2, 4-decadienal), which were found to be effective attractant ingredients in an indoor attractant activity study. However, a study of female pheromones derived from cuticular washes identified the pheromones (Z)-7-heptacosene, (Z)-7-nonacosene, and (Z)-9-nonacosene, and these pheromones were found to be active in laboratory assay experiments and posited to be short-distance contact pheromones, which meant that the actual effective distance must be evaluated in field experiments [33]. In addition, the female adult carries and spreads the symbiotic bacterium Amylostereum areolatum [34], which is injected into the tree trunk with a phytotoxin when eggs are laid, causing the host to become debilitated and eventually die. It can be seen that the symbiotic fungus is involved in a crucial part of the reproductive oviposition process of the wood wasps, and its volatiles also perhaps play an important role in the interactions between S. noctilio and the environment. Therefore, we combined knowledge from previous experiments and the established habits of wood wasps to develop a new attractant based on female sex pheromones and symbiotic fungal volatiles. In previous studies using antennae transcriptome data and qPCR analysis of tissue expression to identify olfactory-related genes, a total of 16 SnocOBPs with different expression levels were obtained. SnocOBP7 was chosen for this study because of its genderspecific expression pattern (it is expressed only in male genitalia) [35], which allowed us to identify gender-specific attractants through reverse chemical ecology and clarify its biological interaction mechanism.
Is SnocOBP7 related to the recognition and interaction of odor molecules? How does it play a regulatory role in protein-ligand interactions? Protein-ligand interactions (PLI) are important processes in organisms. In pest control, they are of great significance for the development and design of attractants and understanding the molecular mechanism of interactions, aiming to clarify the molecular interaction that occurs in the specific binding region. On the basis of the specific expression of SnocOBP7 in the olfactory system of male wood wasps, this study explored the binding interaction relationship between SnocOBP7 and four types of odor molecules: male aggregation pheromones, female sex pheromones, host plant volatiles, and symbiotic fungal volatiles. In order to dynamically characterize the binding between SnocOBP7 and different chemical volatiles and identify the key residues involved in the binding process, we performed 50 ns MD simulations of the ligand-connected complexes. After obtaining stable representative conformations and identifying the ligands with the best binding to SnocOBP7 [36], we performed binding free energy calculations, pairwise per-residue free energy decomposition, and computational alanine scanning (CAS) to determine their key amino-acid residues. The prediction of key protein binding sites can guide the structural analysis and function prediction of protein complexes, as well as design molecules that can regulate biological functions at the system level [37]. The results of this study demonstrate the presence of functional OBPs with gender-specific expression in S. noctilio and provide a foundation for the development of improved methods for prevention and control through studying its protein-ligand interaction mechanism.

RNA Extraction and cDNA Synthesis
Adult samples of S. noctilio were collected from five pieces of damaged wood from Xindian Forest Station in Durbert Mongolian Autonomous County, Daqing City, Heilongjiang Province. After the adults emerged and flew for the first time, the external genitalia of the male adults were collected and dissected, immediately frozen in liquid nitrogen, and stored at −80 • C for later use.
RNA was extracted from the frozen male genitalia using the RNeasy Mini Kit (QIA-GEN, Hilden, Germany). The quality and concentration of RNA were analyzed using an ultramicro-spectrophotometer (Nanodrop 8000, Thermo, Waltham, MA, USA). First-strand complementary DNA (cDNA) was synthesized using the Prime Script™ RT Reagent Kit with gDNA Eraser Kit (Takara, Dalian, China) following the manufacturer's instructions. The cDNA was subjected to PCR amplification, and all PCR products were stored at −20 • C.

Gene Cloning and Verification of the SnocOBP7 Sequence
The validity of the putative OBPs was confirmed by cloning and sequencing ORFs using corresponding primer pairs. On the basis of the sequence of SnocOBP7, primers were designed using Primer3plus (http://www.primer3plus.com/cgi-bin/dev/primer3plus.cgi; accessed on 10 January 2021) to clone the entire ORF of SnocOBP7 (Table 1). Table 1. Primers used for cloning of SnocOBP7.

Name Sequence
SnocOBP7F (5 to 3 ) GGCGGACATTAGAAGAGACTGT SnocOBP7R (3 to 5 ) TGCTTCAAGATCTCGGGCTG ExTaq DNA polymerase (TransGen, Beijing, China) was used to amplify individual sequences under the following conditions: initial denaturation at 95 • C for 5 min, 94 • C for 30 s, 55 • C for 30 s, then 35 cycles at 72 • C for 30 s, and a final extension at 72 • C for 10 min. The PCR products were gel-purified using an AxyPrep DNA Gel Extraction Kit (Axygen, Union City, CA, USA) and then subcloned into the pEasy ® -T3 vector (TransGen, Beijing, China). Finally, the product was sequenced to identify whether the target fragment was inserted correctly.

Homology Modeling and Molecular Docking
The program SWISS-MODEL (http://swissmodel.expasy.org/; accessed on 13 November 2020) was used to search for the best suitable model template for 3D structure construction [39]. In the results, Locusta migratoria OBP1 (PDB ID: 4PTI) was selected for homology modeling because it had the highest identity (32.69%) with SnocOBP7. From 20 models built by Modeller 9.25 [40][41][42][43], the best initial model was selected according to the lowest discrete optimized protein energy (DOPE) score and refined by UCSF Chimera1.14 [44,45]. The residue compatibilities and stereochemical rationalities of the model were examined using the online program SAVEs v6.0 (https://saves.mbi.ucla.edu/, accessed on 16 November 2020). Molecular docking was carried out using the AutoDock Tools 1.5.6 package [46]. The required ligand files were downloaded from Pubchem (https://pubchem.ncbi.nlm.nih.gov/; accessed on 25 November 2020) and were divided into three categories: host plant volatiles, male and female pheromones, and symbiotic fungal volatiles. The size of the Autogrid box was determined by Protein plus DoGSiteScore (https://protein.plus/; accessed on 25 November 2020). The optimal conformation of the SnocOBP7-ligand complex was determined by the lowest estimated free energy of binding and subjected to visual analysis using PyMol2.3.0 [47].

Molecular Dynamics (MD) Simulations of SnocOBP7 and 11 Ligands
After obtaining protein conformations as described above, all MD simulations of SnocOBP7-ligand complexes were carried out using the GROMACS2019.6 software package [48] with the AMBER-ff99sb-ildn force field [49] for protein SnocOBP7 and the Am-berTools 18 GAFF force field [50], which was optimized with the ACPYPE script [51] for ligands. Ligand information is listed in Table 2 and Figure 1. First, Na + or Cl − ions were added according to the total charge of the system. For example, the SnocOBP7-(Z)-7heptacosene complex contained 6082 TIP3P water molecules and the total charge was +8; hence, eight Cl − ions were added for charge neutralization. System energy minimization was achieved by a conjugated gradient (CG) method, after which NVT and NPT (N means the number of particles, P means pressure, T means temperature, and any letter in the name means the value is constant) ensembles for position-restricted MD simulations were performed to relax water solvents. MD simulation was performed for 50 ns with a time step of 2 fs with constant temperature (298.15K) and pressure (1 bar) coupling using a V-rescale thermostat and a Parrinello-Rahman barostat, respectively. The particle mesh Ewald (PME) algorithm was used to calculate the long-range electrostatic interactions, while the LINCS algorithm was applied to constrain all covalent bonds with hydrogen atoms. The cutoff radii of both Coulomb and van der Waals interactions were set to 10 Å. The coordinate information of all systems was recorded every 2 ps, and 25,000 configurations were dumped into the corresponding trajectories for an analysis of the correlations between the conformations. The equilibrium of SnocOBP7 with 11 ligands was analyzed according to the degree of curve fluctuation of the root-mean-square deviation (RMSD). The gmx_cluster module was used for cluster analysis to distinguish different structures and define the representative dominant conformations [52].

Binding Free Energy Calculation and Per-Residue Free Energy Decomposition
The molecular mechanics/Poisson-Boltzmann surface area (MM/PBSA) method [53] and g_mmpbsa script [54] were selected to calculate the free binding energy of SnocOBP7 with 11 ligands. The energy contribution of each residue was further decomposed into van der Waals and electrostatic energy, polar solvation free energy, and nonpolar solvation free energy. The binding free energy (∆G bind ) was defined as follows: The binding free energy (∆G bind ) was equal to the difference in energy of the complex (∆G complex ), receptor protein(∆G receptor ), and the ligand (∆G ligand ); the binding free energy (∆G bind ) was evaluated by calculating the sum of the molecular mechanics energy (∆E MM ), the solvation free energy (∆G sol ), and the conformational entropy effect on binding (−T∆S) in the gas phase; ∆E MM was the sum of the potential energy of the internal bond energy (∆G internal ), electrostatic energy (∆G ele ), and van der Waals energy (∆G vdw ); ∆G PB and ∆G SA were the solvation free energies of polar and nonpolar solvents, respectively. In addition, −T∆S was neglected because of high computational cost [55], and the residues of proteins with free energy contributions greater than −1.00 kcal·mol −1 can be considered as key contributors to binding affinity [56].   Table 2 for details.

Computational Alanine Scanning (CAS)
Computational alanine scanning (CAS) is an energy-based method of identifying key amino-acid residues involved in protein-ligand binding [57], and it can be used to illustrate the contributions of key contributors to binding based on the results of MM/PBSA. All alanine mutation calculations were performed using Discovery Studio Client 2019 (DS, Accelrys Inc., USA) according to the publisher's protocol. The calculation principle was as follows: In this formula, ∆∆G mut is the binding free energy difference between the wild-type and mutants, ∆G bind (mutant) is the binding free energy of the mutants, and ∆G bind (wildtype) is the energy of the wild-type. From the value of ∆∆G mut , the effect of the mutant residue on the binding affinity can be determined. If ∆∆G mut is between −0.5 and 0.5, the mutant residue has no significant effect on the affinity; if it is greater than 0.5, the mutant residue weakens the interaction and reduces affinity; if it is lower than −0.5, the mutant residue enhances the binding affinity for the protein-ligand interaction. Finally, we selected residues with ∆∆G mut greater than 0.5 and analyzed their effects on the protein-ligand system [58].

Sequence Analysis and Homology Modeling
After obtaining the complete sequence of SnocOBP7, ORFs were identified and BLASTP analysis was performed. Five OBPs from other Hymenoptera species were selected for multiple sequence alignment with SnocOBP7 (specific information is shown in Figure 2). The results showed that SnocOBP7 contains 131 amino acids with six conserved cysteine residues. Amino acids 1-22 at the N-terminus are the signal peptide region (Figures 2a and 3). According to the results of the bioinformatics analysis, SnocOBP7 is a small-molecule hydrophobin with a molecular weight of 15.29 ku, which performs secretory functions outside the cell. In summary, SnocOBP7 conforms to the basic characteristics of OBPs. In order to find the most suitable template to create the 3D structure of SnocOBP7, we first identified Locusta migratoria OBP1 (PDB ID: 4PTI) (Figure 2b) as the currently available template with the highest homology (identity = 32.69%), and it was used as the initial template for optimization. Finally, we obtained an optimized model of SnocOBP7 (Figure 2c).  Before molecular docking, the quality of the model must be assessed to ensure its usability. As shown in the Ramachandran plot ( Figure S1a), 92.91% of the amino-acid residues of the SnocOBP7 model were located in the most favored regions (A, B, L, the red region in the figure), and all non-Gly residues were located in the allowed regions. In addition, the overall quality factor of ERRAT was 98.0198, and 88.99% of the residues were found to meet the Verify_3D standard, which strongly suggested that the model was of high quality ( Figure S1b,c). By observing the 3D protein structure of SnocOBP7, the model was found to possess seven α-helices and a hydrophobic cavity. Classical OBPs generally have several common structural characteristics, including six or seven α-helices, three disulfide bonds, and a hydrophobic binding cavity [59]. In addition, in classical OBPs, the α7 helix at the C-terminus forms the wall of the hydrophobic binding cavity.

Stability of SnocOBP7-Ligand Complexes in MD Simulation
The putative binding pocket of SnocOBP7 was predicted in the 3D model (shown in the yellow, green, and purple grid area) ( Figure S2). The yellow area in the hydrophobic cavity was identified as the docking area where the SnocOBP7-ligand complexes were generated. According to the time-evolution root-mean-square deviation (RMSD) curves of all 11 systems (Figure 4), most of the SnocOBP7-ligand complex systems could eliminate spatial conflicts and reach equilibrium at 20 ns, and all systems converged within 50 ns, indicating that 50 ns MD simulations have practical significance and can be used as a basis for exploring the global conformation of SnocOBP7. The average RMSD value of the 11 systems ranged from 2.7 Å to 3.5 Å. The RMSD curve stabilized after 20 ns; the maximum value of the standard deviation of the RMSD was 0.23 Å, and the minimum value was only 0.09 Å. The root-mean-square fluctuation value (RMSF) reflects the local motion characteristics and degree of freedom of the protein secondary structure elements when combined with the ligand ( Figure S3). The RMSF curve showed that the largest fluctuation was in α4 (amino acids between 54 and 60) and α5 (amino acids between 66 and 76), and there was a common peak in the 60-65 residue region (loop and part of α5); the fluctuations in the N-terminus and C-terminus were also large. The maximum RMSF values of S5, S6, S9, and S11 were all located at the 71st amino acid (arginine).
Cluster analysis of the complexes was performed, and the complex with the highest rate of occurrence was selected as the representative conformation for free binding energy decomposition [60,61]. The free binding energy (∆G bind ) of the SnocOBP7 complexes was decomposed into the van der Waals energy (∆G vdw ), the electrostatic energy (∆G ele ), the polar solvation energy (∆G PB ), and the apolar solvation energy (∆G SA ) ( Table 3). The ∆G bind values of S1 (−55.656 ± 0.351 kcal·mol −1 ) and S2 (−56.783 ± 0.260 kcal·mol −1 ) were significantly greater than those of the other complexes, and the ∆G bind of S10 was also relatively high (−31.057 ± 0.154 kcal·mol −1 ). Among the four types of calculated energy, ∆G vdw had the highest value (as high as −62.877 ± 0.220 kcal·mol −1 in S2), indicating that it was the main force promoting the formation of SnocOBP7 complexes [62]. ∆G ele and ∆G SA were weak, but they still played a positive role in SnocOBP7 complex formation; in general, the effect of ∆G SA was slightly stronger than that of ∆G ele . ∆G PB was greater than 0, which showed that it had a negative effect and inhibited the formation of SnocOBP7 complexes. When the side-chains of nonpolar amino acid residues form the three-dimensional structure of the protein, they twist and fold together to form the nonpolar region of the active site, which is known as the hydrophobic pocket. Polar residues were generally located in the opening or bottom of the pocket due to their hydrophilicity, which was why ∆G PB inhibited binding between SnocOBP7 and ligands.  Table 1 for specific information regarding the 11 systems and ligands shown here.

Binding Mode Analysis of SnocOBP7-Ligand Complexes
According to the results described above, we selected the three systems with the lowest free binding energy (S1, S2, and S10) for subsequent analysis. Using cluster analysis, several clusters with different frequencies that appeared during the process of protein binding with ligands were obtained. The most frequently appearing cluster was regarded as the typical conformation of the system. S1, S2, and S10 generated nine, eight, and five clusters, respectively, and the appearance rates of the most common clusters were 38.21%, 39.53%, and 66.11%, respectively (Table S1). The most commonly appearing cluster (cluster 1) of each system was selected for binding mode analysis and combined with the key amino-acid residues predicted by MM-PBSA for visual display (Table 4, Figure 5). As shown in Figure 5, the hydrophobic cavity of S1 was composed of eight amino-acid residues (Leu36, Phe39, Met40, Ile45, Leu61, Ile102, Ile106, and Phe108), whereas that of S2 was composed of five amino acid residues (Leu36, Phe39, Leu61, Leu103, and Phe108), and that of S10 was composed of four amino acid residues (Met40, Ile45, Leu61, and Ile106). Leu, Phe, Ile, and Met were all nonpolar amino acids. These hydrophobic residues were located inside the protein and formed a hydrophobic core facing all directions, maintaining the tight structure of SnocOBP7. The ligands were firmly anchored in the hydrophobic binding cavity; in other words, SnocOBP7 and the ligands could bind stably in it. The side chains of nonpolar amino-acid residues were the structural basis for the formation of hydrophobic interactions, and the hydrophobic force was the main force for ligand binding. Therefore, follow-up studies on these nonpolar residues were conducted. By observing the local motion characteristics of amino-acid residues, it was found that the RMSF values of the abovementioned binding residues were generally low in the 50 ns MD simulation. The minimum RMSF value of S1 was only 0.356 Å (Leu36), whereas that of S2 was 0.361 Å (Leu36), and that of S10 was 0.424 Å (Met40); each of these values was calculated for the common depression formed by residues 36 to 40. The other key residues were also located in the concave area of the RMSF curve shown in Figure 6. The RMSF values at other sites, such as the N-terminus, the C-terminus, and a small peak near Lys64, were higher and unstable. In summary, these key amino-acid residues could interact with small ligand molecules and form stable complexes. Figure 6. The root-mean-square fluctuation (RMSF) curve of SnocOBP7-ligand complexes (S1, (a); S2, (b); S10, (c)).

Per-Residue Free Energy Decomposition
The MM-PBSA, which has been widely used to evaluate the binding free energy of complexes, was used to calculate the contribution of each amino-acid residue of the SnocOBP7-ligand complex to its binding free energy. As shown in Figure 7 and Table 4, there were eight (Leu36, Phe39, Met40, Ile45, Leu61, and Ile102), five (Leu36, Phe39, Leu61, Leu103, and Phe108), and four (Met40, Ile45, Leu61, and Ile106) key amino-acid residues that contributed more than −1 kcal·mol −1 to the binding free energy for S1, S2, and S10, respectively. Therefore, these key amino acids were selected for energy decomposition analysis. Figure 7. Per-residue contribution to the binding free energy of SnocOBP7-ligand complexes (S1, (a); S2, (b); S10, (c)) calculated from the equilibrated conformations. The residues contributing more than −1.00 kcal/mol to the binding free energy are marked by the red line.
The highest value of ∆G bind was −3.374 kcal·mol −1 and the lowest was −1.08 kcal·mol −1 . As shown in Table 4, one amino-acid residue (Ile106 of S1) had a ∆G bind value that exceeded −3 kcal·mol −1 because of its high ∆E MM (−3.076 kcal·mol −1 ). The ∆E MM of Phe108 in S2 also exceeded −3 kcal·mol −1 ; however, due to the strong inhibitory ∆G PB (2.188 kcal·mol −1 ), the final ∆G bind was only −1.744 kcal·mol −1 . ∆G PB was mostly positive, which was unfavorable for binding, while ∆G SA was negative and had a positive, but weak effect on binding; ∆E MM was the main driving force for binding.

Computational Alanine Scanning (CAS)
Computational alanine scanning is an effective tool for evaluating the free energy change caused by the substitution of amino acids with alanine (Ala) residues [63]. Therefore, we used the key residues of S1, S2, and S10 obtained above as the target for CAS to determine the potential for these residues to participate in binding ligands. As shown in Table 5, replacement of five residues of S1 (Leu36, Met40, Leu61, Leu73, and Phe108five 5 residues of S2 (Leu36, Phe39, Leu61, Leu103, and Phe108), and two residues of S10 (Met40 and Leu61) with Ala resulted in a mutation energy (∆∆G mut ) change greater than 0.5 kcal·mol −1 . These residues were also identified by the per-residue free energy decomposition analysis, demonstrating that the results were consistent and reliable. The ∆∆G mut values of other residues from the per-residue free energy decomposition analysis results were less than 0.5 kcal·mol −1 , indicating that replacement of these residues with alanine had little effect on the binding affinity of the SnocOBP7 protein with ligands. According to the DS assessment criteria [64,65], the residues identified via CAS may have important effects in stabilizing SnocOBP7-ligand complexes, and amino-acid changes at these positions would be expected to change the active conformation of SnocOBP7, ultimately changing the ability of the protein to bind ligands.

Discussion
It is well established that differential expression of OBPs is associated with differences in physiological structure and physiological function, and gender-specific OBPs have significant promise in interactions with odor molecules and pest control and monitoring applications. In insects, OBPs are most abundant in the antennae, but specific expression of OBPs also occurs in the abdomen, head, and external genitalia. Several studies have shown that differential expression of OBPs is widespread in Hymenoptera insects. For example, in Cotesia vestalis, six OBPs were found to be highly expressed in the antennae of females, while three different OBPs were expressed in the antennae of males [66]. However, the three OBPs of Sclerodermus alternatusi were found to be expressed only in the abdomen of female adults; hence, it was speculated that they played a specific role in oviposition behavior [67]. In this study, due to the gender and tissue specificity of SnocOBP7, we speculated that it may play an important role in guiding male adults to complete reproductive mating. Therefore, with the aim of establishing new directions for the development of compounds to control and monitor S. noctilio, we used our previous antenna transcriptome data to identify OBP7 as a gender-specific OBP expressed only in the external genitalia of male adults. Subsequently, bioinformatics analysis and molecular dynamics simulations were carried out to reveal the interaction mechanism between SnocOBP7 and various odor molecules.
To explore interaction avenues for improving S. noctilio attractants, binding between SnocOBP7 and 11 types of odor molecules (including host plant volatiles, symbiotic fungal volatiles, and sex pheromones) was simulated. Subsequently, energy decomposition analysis of the SnocOBP7-ligand complexes was carried out. Firstly, the screening results showed that (Z)-7-heptacosene and (Z)-7-nonacosene had the highest binding affinity for SnocOBP7 among female wood wasp sex pheromones. The strategies used by S. noctilio to locate mates during courtship are not well understood, mainly because mating generally occurs in the high canopy layer and is, thus, difficult to observe. Most studies have indicated that male and female wood wasps are attracted by host plant volatiles [68] or light [69] and, therefore, congregate in the same areas in the upper canopy during mating. After meeting in the upper canopy, male and female wood wasps must locate each other through various signals to increase the possibility of successful mating. Before mating, the male always approaches from behind the female and touches the end of the female's body with his genitalia [70], which suggests that the male wood wasps use their genitalia to recognize females. This is a specific behavioral interaction that odor molecules regulate insect physiology and then reflect behavior. According to current research progress, female sex pheromones were predicted to be short-distance pheromones. To explore the scope of its effectiveness in the actual application process, it was necessary to carefully adjust the attractant formula and additive dosage in the follow-up work. If the action distance was short, it was compounded with plant volatiles to observe whether it could achieve a close-range attracting effect and enhance the overall effect. Moreover, previous studies have also indicated that female sex pheromones were obtained from female adult body surface extraction. Our findings suggest that (Z)-7-heptacosene and (Z)-7-nonacosene might be used by male S. noctilio to locate females; therefore, these female pheromones could guide the development of gender-specific attractants.
Secondly, among the tested symbiotic fungal volatiles, the binding affinity of (−)globulol far exceeded that of the volatiles, which may be linked to the mating effect of symbiotic bacteria and male adults and an important factor in the regulation of biological interactions. The symbiotic fungus of S. noctilio is A. areolatum, which has weak pathogenicity but can degrade the lignocellulose of host plant and facilitate larval feeding [71]. (−)-Globulol is notable for its attractiveness to female adult S. noctilio [72], but relevant behavioral experiments have not been performed using male adults. Our computer simulations demonstrate that (−)-globulol is capable of binding strongly with male-specific SnocOBP7. Female adults inject eggs and symbiotic fungi into the host plant at the same time when laying eggs. Therefore, it is likely that SnocOBP7, which is expressed only in male genitalia, binds to symbiotic fungal volatiles such as (−)-globulol to signal to adult males that females could be present. This reminds us that, in the follow-up control work, it may be possible to influence the behavioral interaction between male adults and the environment by adjusting the expression of (−)-globulol, so as to achieve the purpose of attracting and trapping. Males may be more inclined to fly to areas with a high density of symbiotic fungal volatiles, because such areas are likely to contain weakened host plants and more mating resources. Symbiotic fungal volatiles are released from the fungal mycangium of females or at oviposition sites on damaged trees [73]. Therefore, the fungal volatile (−)-globulol may also be used in the courtship process of wood wasps as an attractant. S. noctilio and its symbiotic fungi have gradually formed a highly evolved and strict obligate dependency mutualism [74]. Because S. noctilio locates suitable egglaying locations through symbiotic fungi, it is possible to mix symbiotic fungal volatiles into attractants formulated with plant-derived volatiles and pheromones to improve their control efficiency. (−)-Globulol, which has a strong binding affinity with males and has been proven to be effective in attracting females, may be an important regulatory odor factor that has been neglected and has not been used in field traps. Therefore, analysis of simulations of SnocOBP7-ligand binding interactions could lead to the development of new strategies to limit the reproduction of S. noctilio via intervention at the courtship stage.
Lastly, in our simulation analysis, SnocOBP7 did not bind strongly with the selected male aggregation pheromone, presumably because the antennae of adults generally perceive this type of pheromone. Adult males can cluster in the wild [75], mainly because they produce aggregation pheromones, which cause individuals of the same sex to gather. These male pheromones should play similar role to host volatiles that attract both sexes to gather in the high canopy layer. Different OBPs have different interactions with the environment and have different sensitivity to different odor molecules, which are all related to their expression levels in different tissues.
Different types of free binding energy have different numerical values, which represent different contributions to binding affinity. Among the four energies we analyzed, the value of ∆G vdw was obviously the highest, indicating that it was the main force driving binding [76]. Indeed, ∆G vdw was higher than the final ∆G bind because the negative offset effect of ∆G PB on the binding could not be ignored. In contrast, although ∆G ele and ∆G SA played a role in promoting binding, their contributions were relatively weak. The decomposition of free binding energy also showed that stable SnocOBP7-ligand binding was attributed to the hydrophobic interaction between SnocOBP7 and ligands. OBPs have a 3D hydrophobic binding cavity, and hydrophobic odor molecules are usually anchored in it to achieve relatively stable binding with proteins. For example, analysis of the binding energy of A. lepigone AlepPBP1 showed that the main driving force promoting its binding was the van der Waals force, demonstrating the key role of hydrophobicity [77]. In the CpomPBP2-Dod system of Cydia pomonella, in addition to the dominance of the van der Waals force, the electrostatic force also plays a significant role due to the formation of strong salt bridges and H-bonds [78].
Clarifying the binding mechanism of OBP and odor molecules and deciphering their biological interactions with the environment has always been a hot topic in the field of OBP research. The per-residue free energy decomposition showed that each system had several key amino-acid residues, which together formed a hydrophobic binding cavity. For example, the eight amino-acid residues of S1 that formed a hydrophobic binding cavity were Leu36, Phe39, Met40, Ile45, Leu61, Ile102, Ile106, and Phe108. As shown in Figure 5, the ligand was surrounded by these residues, allowing it to be pulled in different directions to achieve balanced binding. Among these amino-acid residues, only the ∆G bind of Ile106 was greater than −3 kcal·mol −1 , and ∆E MM was the main contributing force. From the perspective of overlapping key amino-acid positions, Phe108 and Leu36 appeared as key amino acids in S1 and S2. Therefore, it could be inferred that they were the key sites involved in binding of OBP7 with female sex pheromones. Met40 appeared in both S1 and S10, suggesting that it plays a key role in binding to female sex pheromones and symbiotic fungal volatiles; therefore, Met40 could play an important role in regulating the courtship behavior of males. Leu61 was identified in all three systems, suggesting that it could play key roles in mating and host plant identification. More importantly, after obtaining the prediction results of the key amino-acid binding sites, it is possible to induce mutations at these sites from the perspective of the protein SnocOBP7 to achieve the reverse regulation of the biological interactions. Therefore, CAS was performed on the key residues derived from the results described above. If a particular amino-acid residue was replaced by alanine and had a significant effect on the binding energy, then it was determined to play a key role in ligand binding. The key residues, Met40, Leu61, and Ile106, were found to be in the same region, which contained a depression according to the RMSF analysis, which is consistent with our simulation and prediction results. Similarly, in a study of Anopheles gambiae AgOBP1, molecular dynamics analysis confirmed that eugenyl acetate was a better insect repellent than DEET and also revealed the main features of the binding site of AgOBP1 [79].
According to reverse chemical ecology, after obtaining the key amino-acid residues and odor ligands for binding, OBPs can be used as a major entry point for the development of new and environmentally friendly attractants and as a powerful tool for interfering with or regulating biological interactions. However, although female pheromones were obtained from the body wall of female adults through chemical leaching in previous studies, actual molecular experiments are needed to determine their potential utility in control and monitoring applications. For example, in a study of the olfactory function of OBP2 in Apis cerana, site-directed mutagenesis was used to verify the key amino-acid binding sites predicted by molecular docking [80]. In this study, the results of our computer simulation provide new ideas for explaining the protein-binding interaction mechanism and the subsequent pest control and prevention from the perspective of female sex pheromones and symbiotic fungi, which could lead to the development of improved attractants via the principle of chemical ecology.