Docking studies for TBM and AF models of MhOR5:
The three experimentally resolved structures available for MhOR5 in PDB are the apo-structure, 7LIC (PDBID: 7LIC), and 7LID (PDBID: 7LID) and 7LIG (PDBID: 7LIG), where it is bound to eugenol (EOL) and N,N-diethyl-m-toluamide (DEET), respectively12. EOL is a natural product insecticide and DEET a synthetic repellant35,36. We compared our pre-docking models against 7LID (minus its ligand) because it was solved with higher resolution than 7LIC (2.9 vs 3.3 Å).
As expected from the preliminary study25, both models accurately predicted the MhOR5 fold but TBM consistently outperformed AF across a range of measures. Template modeling scores were 0.55 and 0.54 for MhOR5_TBM and MhOR5_AF, respectively, but the RMSD between 7LID and MhOR5_AF was higher (2.37 Å) than that between 7LID and MhOR5_TBM (0.31Å) (Fig. 1(a), (b)). Kinks in helices 2 and 7 of 7LID were accurately represented in both models (Fig. 1(c)) but there were prominent differences between 7LID and MhOR5_AF, but not MhOR5_TBM, in extracellular regions, specifically helices 3, 5, 6 and 7 (Fig. 1(c)). These differences were also reflected in structure-based sequence alignments at the end of helices 3, 5 and 7 and the beginning of helix 6 (Fig. 1(l)). Moreover, helices 4 and 5 were longer in MhOR5_AF than 7LID (Fig. 1(c), (d)), which might be because intracellular loop 2 (ICL2) connecting helices 4 and 5 was not resolved in the empirical MhOR5 structure.
We then compared the binding conformations and interactions of the two modeled structures with EOL and DEET with those in the corresponding 7LID and 7LIG empirical structures. To do this we first downloaded the structures of EOL and DEET (Fig. 2(a)) from the PubChem database and re-docked them to 7LID and 7LIG, respectively. This was done to validate the recovery of RLIs after docking and, as such, to then provide controls for the comparison of the ICM docking scores39 and contact residues obtained with the models. The re-docking recovered 100% of the experimental RLIs in both the 7LID_EOL and 7LIG_DEET complexes, with no additional contacts, and yielded control ICM docking scores of -24.33 for 7LID with EOL and − 20.41 for 7LIG with DEET (Table S1).
The corresponding ICM docking values for MhOR5_TBM and MhOR5_AF with the two ligands were all higher than the control values but those for MhOR5_TBM, -19.20 and − 12.96 respectively, were both closer to the control values than those for MhOR5_AF, -13.69 and − 8.34 respectively (Table S1). Inspection of the structures confirmed that the binding conformations of both EOL and DEET inside the binding pocket of MhOR5_TBM were closer to the experimental structures than those inside MhOR5_AF (Figure S1). Consistent with the docking scores, the greatest discrepancy was between DEET in MhOR5_AF and 7LIG_DEET.
The ligand binding pocket of MhOR5 involves helices 2, 3, 4, and 6 and, using our generic residue numbering system for each helix (see Methods), the contact sites for 7LID with EOL are 2.3, 2.6, 2.7, 3.26, 3.29, 3.30, 3.33, 4.13, 4.17, 6.7 and 6.10, and those for 7LIG with DEET are 2.3, 2.6, 2.7, 3.25, 3.26, 3.29, 3.30, 3.33, 4.17, 6.6, 6.7 and 6.10. These residues are mostly aromatic and hydrophobic (Figure S2), and all are non-polar. Functional assays of site-directed mutants have also implicated homologs of six of these contact sites, 2.3, 2.7, 3.29, 3.33, 4.13 and 4.17, in ligand binding in ORs from other insect species11,17–22,40,41. Some other sites close to these, namely 3.34, 4.13, 4.17, ECL3, and 6.3, have also been implicated in these studies (Table S2).
All the contact residues predicted when EOL and DEET were docked into MhOR5_TBM were also contact residues in the empirical structure but MhOR5_AF only captured 54% of the contact residues for EOL in the empirical structure, and just 16% of those for DEET in that structure (Figure S2). The low values for the MhOR5_AF modelare consistent with its high ICM scores and failure to properly orient either ligand, particularly DEET, in the binding pocket. Notably, the MhOR5_TBM model also predicted three other contact residues not in the empirical data (2.2, 4.16, and 6.3) while the corresponding number for MhOR5_AF was five (3.37, 3.38, 3.39, ECL3, and 6.4).
Molecular dynamics and simulations for TBM and AF models of MhOR5:
MD simulations were executed for all four combinations of the two models and two ligands. We also executed MD simulations for the empirical 7LID_EOL and 7LIG_DEET complexes to validate our MD simulation protocol (through recovery of the experimental RLIs) and to provide positive control values for post-dynamics energy calculations. As detailed below, we found the simulations generally improved the performance of both models but TBM continued to perform better than AF.
Post-dynamics binding energy values for the control 7LIG_DEET (GB: -32, PB: -9.02) and 7LID_EOL (GB: -35.75, PB: -9.62) complexes were notably low, as expected, and those for the simulations using the two TBM complexes (GB: -26.6, PB: -5.49 for TBM_MhOR5_EOL and GB:-34.13, PB: -12.59 for TBM_MhOR5_DEET ) were generally closer to the corresponding controls than those based on the AF models (GB: -24.06, PB: -4.2 for AF_MhOR5_EOL and GB:-26.27.13, PB: -3.97 for AF_MhOR5_DEET ) (Fig. 2(b)).
Post-dynamics stability for both the complex and ligand was assessed from the standard deviations (CSD and LSD, respectively) for the last 10ns of the standard 50 ns simulations. Values for the control empirical structures (0.09 for 7LID_EOL and 0.10 for 7LIG_DEET) were again low, as expected. Both the TBM and AF complexes likewise generated low CSD values, with the TBM values lower than the AF values (0.07 and 0.16 respectively) for MhOR5_EOL but no difference between them (0.15 for both) for MhOR5_DEET. However, the LSD values for the modelled complexes were generally somewhat higher than their experimental counterparts, with the AF generating lower values than TBM in this case (0.33 and 0.13 for TBM and AF, respectively for MhOR5_EOL and 0.28 and 0.19 for TBM and AF, respectively for MhOR5_DEET).
The post-simulation position and orientation of EOL and DEET in the TBM complexes largely overlapped with their positions in the empirical structures, whereas the position of EOL in the AF complex differed slightly from its experimental position (Fig. 2(c)) and that of DEET in the AF complex was quite distinct from its experimental position (Fig. 2(c)). These differences were also reflected in the post-simulation RLIs of the various complexes (Fig. 2(d)). Overall, TBM and AF recovered 100% and 45%, respectively, of the experimental RLIs for EOL, and 91% but just 0.08%, respectively, of them for DEET (Figs. 2(d), S3). The only differences in the post-simulation RLI sets predicted by TBM from the empirical set were the absence of 3.25 and the presence of 2.10, 3.32, 4.16, and 6.3 in the respective predicted sets. Notably too, recovery rates of the empirical RLIs were generally higher post- than pre-simulation (100, 54, 100 and 0.16% c.f. the 100, 45, 91 and 0.08%) respectively, above (Figure S4).
Hydrophobicity-based MSAs of key helices for the TBM
The ability of TBM to predict the known functionalities of our panel of five B. dorsalis and three B. minax ORs for which no empirical structural data were available was dependent on optimizing the alignment of the sequences of their key helices to those in the known structures of the template1. This is challenging for insect ORs for two reasons. One is that sequence similarities between any two ORs within a species can be as low as 17%9. The other is that, unlike the situation for mammals, very few highly conserved residues to use as anchor-points for the alignments are known for insect ORs; the only well-known motif in insect ORs is tryptophan at the end of predicted helix 69. Exploiting the crucial importance of hydrophobicity to the structure and function of the key helices, and drawing on multiple transcriptome datasets for B. minax14,42–44 and B. dorsalis17,31,45–48 (Table S3, Supporting information Sheet 1), we therefore developed an extensive database of hydrophobicity-guided alignments of those helices in MhOR5 and the 130 ORs from those two species for which high confidence sequence data were available.
As expected, a MSA of the 84 BdORs and 46 BmORs generated through ClustalW showed very little conservation (data not shown) but a hydrophobicity-based, helix-wise MSA of them to 7LID produced by iterative application of AlignME revealed several conserved motifs in most of the seven shared helices (Supporting information: Sheet 2). Conservation was greatest in helix 5 and lowest in helices 1 and 3, although some conservation of hydrophobic residues was evident even in the latter (see e.g., sites 1.6, 1.8, 1.11, 1.12, and 1.15). Highly conserved residues or short motifs identified in individual helices of this MSA were then used as anchor points to guide further manual curation of the alignments of the helices containing them (Fig. 3, Table 1).
Table 1
Conserved sites/motifs identified in four helices. Anchoring sites used to facilitate pairwise alignments between the Bactrocera ORs and MhOR5 are bold. Minimum levels of conservation for the respective sites/motifs in our MSA are also indicated.
Helix | Anchoring site(s) | Conservation (%) |
2 | K/R (2.13) | 97 |
4 | Y/F (4.2), D (4.24), L(4.41) | 75 |
5 | LxxxI/VxxH (5.4), F/LxQF (5.31), H (5.11) | 80, 84, 85 |
6 | Y/F(6.7), W(6.41) | 80, 90 |
Only one anchor point was identified for helix 2, site K/R (2.13), which occurred in 97% of the sequences in the dataset. However, three such sites were identified for helix 4, namely Y/F (4.2), D (4.24) and L (4.41), which all occurred in over 75% of the sequences. Helix 5 contained three conserved multi-site motifs: LxxxI/VxxH at sites 5.4–5.11 which was conserved in more than 80% of the sequences; F/LxQF at sites 31–34, with F/L (5.31) more highly conserved than QF (84 cf 45%); and L/IC/G at sites 40–41, which was conserved in 65% of the sequences. Since H (5.11) was the most conserved site of all these residues (85%), it was used as the primary anchor point for adjusting the alignment of helix 5. Three conserved sites or motifs were also identified in helix 6; site Y/F (6.7), which was conserved in about 80% of the sequences; CYxG/A, with 78% C (6.19), 62% Y (6.20) and 95% G/A (6.22) conservation across the dataset; and Y/FxxxW at sites 6.37–6.41, which showed over 90% conservation across the dataset. Note that Y/F (6.7) is a contact site for EOL and DEET in the binding pockets of 7LID and 7LIG respectively. W (6.41) alone has previously been reported as a binding motif in D. melanogaster and Anopheles gambiae9 and was used here as the primary anchor point for this helix. All the above conserved sites except Y/F (4.2) and D (4.24) were also detected as motifs by MEME suite (Figure S5).
We could not identify any anchoring sites for helices 1, 3, and 7 and therefore used the alignment generated through AlignME for helices 1 and 7. However, as noted, the available structural and mutational data indicates that helix 3, along with helices 2, 4 and 6, contributes to ligand binding (Fig. 2(f), Table S2), so, following the approach of Jabeen et al49, we manually adjusted its alignment in order to optimise the hydrophobicity correspondence (HC) per site, using the sum of squared differences (SSD) calculator embedded in Bio-GATS50. The approach is detailed and illustrated for helix 3 for BdOR88a and BmOR7a1 in Supporting information: note 1 and Figures S6 and S7, respectively. The adjustments were then validated through simulations; for instance, complexes based on adjusted helix 3 models yielded lower SSDs per residue than AlignME for both BdOR88a_ME and BmOR7a1_linalool (BmOR7a1_LIN) for both GB (-27.62 vs -14.9 and − 28.34 vs -22.65, respectively) and PB (-7.5 vs -4.8 and − 8.78 vs -3.93, respectively). The adjusted, validated helix 3 alignments were integrated into the full length pairwise alignments, serving as the basis for used in subsequent model building.
Testing the abilities of TBM and AF to capture empirical functional evidence
We identified a panel of 63 odorants which have been tested for binding against one or more B. dorsalis or B. minax OR (Fig. 4, Table S4). Five of these complexes have EC50 values in the µM range which is widely regarded as potentially relevant biologically; biokinetic studies of ORs typically define high affinity binding as within the nM to µM range8,51,52. These five complexes plus the structurally validated 7LID-EOL and 7LIG-DEET were selected as a binder cohort for further analysis. A complementary cohort of non-binders was then selected from the many non-binding OR-ligand combinations in the panel. Selection of these was based on maximizing the overlap with the binder cohort in respect of the ligands tested but was otherwise random. In total, the 14 combinations selected involved eight ORs and four odorants (1-octen-3-ol (Oct-3OL), LIN, ME, and EOL) (Figure S8). Two of the ORs, BmOR13a1 and BdOR13a1, were orthologs in the two species53.
TBM and AF modeling were then carried out on the eight ORs in these combinations, using 7LID as a template for the TBM because it was bound to the agonist EOL and most of the other ligands in our binder cohort were also agonists of the receptors. The TBM also used final pairwise alignments of the bactroceran ORs in the cohorts tested with 7LID (Figures S9, S10). Sequence identity between those ORs and 7LID varied from 6.9 to 20% (Table S5).
The docking results for the ORs in the binder cohort modelled using either TBM or AF with their respective ligands did not reflect the corresponding empirical data. The ICM docking scores (Table S6) were higher than those of the MhOR5_EOL control (-24.33) and not significantly different from those in the non-binder cohort (p-values of 0.76 for the TBM and 0.22 for the AF models). Therefore, we proceeded with MD simulations, which, for each model, generated GB, PB, CSD and LSD values as defined above, and a total of 64 RLIs. Thirty RLIs were retained after filtering out several with low variances and high correlations (r2 ≥ 0.7) with other RLIs (Figure S11).
Linear discriminant analysis (LDA) of the energy and stability variables showed that TBM successfully distinguished between the binders and non-binders (Fig. 5(a)), whereas AF left some overlap between the two cohorts (Fig. 5(b)). To probe these results more closely we then used the distributions of normalized values predicted by each model for each energy and stability variable across all the complexes in the two cohorts (Fig. 5(c)).
Considering first the two key energy variables, TBM generated the clearly negative normalized values (<= -0.25) expected of true binders for both GB and PB for five of the seven binder complexes (MhOR5a_DEET, BdOR13a1_Oct-3OL, BdOR7a6_Oct_3OL, BmOR7a1_LIN and BdOR88a_ME), whereas AF only did so for two (BmOR13a1_Oct-3OL and BdOR13a1_Oct-3OL). Conversely, TBM generated the clearly positive normalized values ( > = 0.25) expected of true non-binders for both GB and PB for four of the seven non-binder complexes (BmOR7a1_Oct-3OL, BdOR82a_EOL, BdOR82a_Oct-3OL, and BdOR94b2_ME), whereas AF only did so for one (BdOR82a_EOL). Moreover, TBM only produced the counter-expectation, clearly positive values for both GB and PB for one binder (BmOR13a1_Oct-3OL) and clearly negative values for one non-binder complex (BmOR13a1_LIN), while AF produced clearly positive pairs for four binders (MhOR5_EOL, MhOR5_DEET, BdOR7a6_Oct-3OL, and BmOR7a1_LIN) and three clearly negative pairs for three non-binders (BmOR13a1_LIN, BdOR82a_ME, and BdOR94b2_ME). Thus, while not in agreement with the empirical data in all individual cases, the energy estimates from TBM clearly outperformed those from AF in predicting both binders and non-binders, with AF having little predictive value in either case.
The corresponding analysis for the two stability variables showed neither TBM nor AF consistently predicted either the empirical binder or non-binder cohorts. This was largely because there was much less correlation between the two stability variables than there was between the two energy variables generated by either TBM or AF (Figure S12). Consequently, there were relatively few cases where CSD and LSD were both <= -0.25 or both > = 0.25. However, there was a general trend for TBM to produce the more negative CSD and LSD values indicating higher stability in the empirical binder than non-binder complexes. This was particularly the case for MhOR5_EOL, MhOR5_DEET, BdOR7a6_Oct-3OL, and BdOR88a_ME, respectively, for CSD and BmOR13a1_Oct-3OL, BdOR13a1_Oct-3OL, BmOR7a1_LIN, and BdOR88a_ME, respectively, for LSD. By contrast, the opposite trend was evident in the AF results (BdOR13a1_Oct-3OL, and BdOR88a_ME, respectively for CSD; MhOR5_EOL, MhOR5_DEET, and BdOR7a6_Oct-3OL, respectively for LSD).
Most of the experimentally validated RLIs in the 7LIG and 7LID complexes (2.3, 2.6, 2.7, 3.26, 3.29, 3.30, 3.33, 4.13, 4.17, 6.7, and 6.1012) were frequently found in either or both of the TBM and AF complexes (Figs. 5(d), S13 and S14) and most of the others predicted by the two models also lay in the four helices containing the 7LIG and 7LID RLIs. However, a few helix 1 and helix 5 RLIs were also predicted, the former mainly by AF and the latter mainly by TBM. Helix 5 has been implicated in binding in in silico docking of the ligand VUAA1(N-(4-ethylphenyl)-2-{[4-ethyl-5-(pyridine-3-yl)-4H-1,2,4-triazol-3-yl]sulfanyl}acetamide) with coreceptor Orco in D. melanogaster54 but there is no other evidence as yet to implicate helix 1 in ligand binding by insect ORs. Empirical evidence from site-directed mutagenesis of DmelOR85b has implicated its site equivalent to 3.34 in its binding to 2-heptanone19.
TBM and AF each generated several differences in the frequencies of particular RLIs between the binder and non-binder cohorts but there was little overlap between TBM and AF in these RLIs (Fig. 5(d)). The TBM and AF models also differed in the orientation of the ligands in the binding pockets (Figures S13, S14). At the helix-level, the RLIs in helix 5 implicated by TBM and those in helix 1 implicated by AF were much more frequent in the binder complexes. To confirm the involvement of helix 5 and helix 1 in the TBM and AF binder models respectively, we extended the respective simulation times to 200 ns for one of the strongest binding complexes, namely BdOR13a1 with Oct-3OL. Helix 5 in TBM_BdOR13a1_Oct-3OL and helix 1 in AF_BdOR13a1_Oct-3OL remained consistently present in the simulation trajectories throughout the 200ns. Nor did values for GB, PB, CSD, and LSD change significantly throughout that time (Table S7).
Eleven individual RLIs were also discriminating in so much as each was markedly more prevalent in just one of the four TBM/AF/binder/non-binder cohorts than it was in the other three. These were 4.20, 5.35, 5.36 and 5.39 in the TBM binder cohort, 2.2 in the TBM non-binder cohort, 2.7 and ECL2 in the AF binder cohort, and 4.20, 6.10, 6.11, and 6.14 in the AF non-binder cohort. Thus, inspection of the individual discriminating RLIs again showed prominent representation of helix 5 in the TBM binding dataset, although it did not do so for helix 1 in the AF non-binder dataset.
Direct side-by-side comparisons of predicted parameter values for binder and non-binder complexes were possible for two of the ORs in our panels, BmOR7a1 and BmOR13a1, the latter being of particular interest because it had been the most refractory to TBM modelling (Figs. 5(c), S13 and S14). The binding of BmOR7a1 with LIN and its non-binding with Oct-3OL were both successfully predicted by TBM (GB -0.46 vs 1.07; PB -0.69 vs 0.73; CSD − 0.11 vs 2.48; LSD − 0.74 vs 2.51. Its RLIs with LIN lay in three of the four helices expected above (2–4) plus helix 5 but only three of the four expected (3, 4, 6) with Oct-3OL. However, the AF models for BmOR7a1 yielded predictions opposite to the empirical data for both ligands (GB 1.32 vs -0.31; PB -0.97 vs -0.1; CSD 0.8 vs 0.46; LSD 0.12 vs 0.02), and its RLIs with LIN were in all four of the helices expected plus helix 1 and ECL2, and all four expected plus 5 with Oct-3OL.
On the other hand, the interactions of BmOR13a1 with its binder and non-binder ligands proved problematic to predict for both AF and TBM. AF yielded GB and PB estimates (-1.06 and − 1.47, respectively) correctly predicting its binding to Oct-3OL, although its CSD and LSD estimate (0.28 and 0.59, respectively) indicated that both the complex and Oct-3OL were unstable in the binding pocket (Figures S15(a), S16(a)). However, the estimates of GB, PB and LSD (-1.44, -1.59, -0.46, respectively), albeit not CSD (2.71), which AF produced also incorrectly suggested it bound to LIN. In contrast, TBM yielded GB, PB and CSD estimates (0.69, 1.89, and 1.01, respectively) incorrectly predicting non-binding to Oct-3OL, although the LSD (-0.64) suggested ligand stability in the binding pocket. TBM also failed to predict the non-binding of BmOR13a1 to LIN, although its GB and PB values (-0.24 and − 0.4, respectively) were less extreme than those produced by AF (-1.43 and − 1.59, respectively) and its CSD estimate (0.45) suggested complex instability. Both models also failed to produce RLIs which distinguished the binder and non-binder interactions; AF predicted RLIs in helix 1 and ECL2 as well as the four helices expected with both ligands, while TBM predicted RLIs in three of the four expected helices (2, 4 and 6) plus 5 with both. Importantly also, none of the results for any of the parameters changed materially when the simulations of BminOR13a1 with LIN and Oct-3OL were extended to 200ns (Table S7). Their problems with this OR emphasise the scope for further refinement of both TBM and AF.