The adipokinetic hormones and their cognate receptor from the desert locust, Schistocerca gregaria: solution structure of endogenous peptides and models of their binding to the receptor

Peptide molecular dynamics

The three AKH peptides were built using Maestro (Schrödinger, Inc., New York, NY, USA) and energy minimized using a steepest descent algorithm. NMR restrained molecular dynamic (MD) simulations in vacuum, water and DPC were performed using GROMACS version 5.1.2 (Van Der Spoel et al., 2005). Although only the final results in DPC micelle solution are relevant to the current study, the three different environments were used to adequately search conformational space in a reasonable amount of time. Thus, the vacuum simulation was done at 600 K to allow the peptide to move faster and overcome any conformational barriers. All simulations were performed using the OPLS-AA all-atom force field with a time step of two fs. The LINCS algorithm was used to constraint all bonds. A cut-off of 1.0 nm was used for van der Waals interactions and electrostatic interactions for real space calculations. Vacuum simulations were first done to search conformational space by collecting 100 snapshots of the trajectory during a 10 ns simulation, at 600 K. Each conformation was then annealed to 300 K over 40 ps. Cluster analysis of the resulting structures, using the linkage algorithm of GROMACS and a cut-off of 0.1 nm on the backbone atoms, gave one single large cluster. The conformer in the cluster with the lowest energy was used for simulations in water. Using the tip4p water model, a box containing the peptide and ~7,000 water molecules was constructed. Following equilibration, MD was performed for 10 ns at 300 K under NPT conditions. In total, 200 structures were collected at 50 ps intervals. Cluster analysis was performed as before and the results used in the DPC/water simulations.

For simulations in a water/DPC mixture, the lowest energy structure obtained previously was placed in the center of a seven nm cubic box filled with ~10,000 SPC water molecules and a 52 DPC molecule micelle as obtained from Tieleman, Van Der Spoel & Berendsen (2000). The micelle was translated so that the center of the micelle was at the bottom edge of the box. This meant that, using periodic boundary conditions, half the micelle was at the bottom of the box and the other half was at the top. The peptide was then placed in the center of the box. Energy minimization was carried out using the steepest descent method to a tolerance of 10 kJ/mol or to machine precision. Two stages of system equilibration were performed to solvate the peptide and to achieve a steady state starting temperature, pressure, and density. The first stage of equilibration involved performing MD for 100 ps under NVT conditions at 300 K followed in the second stage by a further 1 ns MD under NPT conditions. The final MD simulation was for 10 ns during which 200 snapshots were collected. Cluster analysis was performed in the same manner as before.

Homology modeling

The primary sequence of the AKH receptor from S. gregaria, Schgr-AKHR, was obtained from the GenBank (GenBank ID: AVG47955.1). Transmembrane (TM) helix predictions were computed online using PSIPRED (http://bioinf.cs.ucl.ac.uk/psipred/) (Jones, 2007; Jones, Taylor & Thornton, 1994; Nugent & Jones, 2009). The results showed that this sequence had seven TM helices (Fig. S1).

The GPCR-ModSim Web server (http://gpcr-modsim.org/) (Rodríguez, Bello & Gutiérrez-De-Terán, 2012) was used for template selection and preliminary sequence alignment of Schgr-AKHR. This server accepts an amino acid sequence as input and searches templates by multiple sequence alignments with the query sequence. GPCR-ModSim identified the human kappa opioid receptor (hĸ-OR) (PDB ID: 4DJH; Wu et al., 2012) as the best template. Subsequently, the sequence alignment between Schgr-AKHR and hĸ-OR was manually edited to remove gaps in the TM domains without disrupting their conserved regions. Finally, homology models for Schgr-AKHR were constructed using Modeler 9v7. This is an automated homology modeling program that performs automated protein homology modeling and loop modeling for the receptor by satisfaction of spatial restraints (Sali, 1995). The quality of the selected model was evaluated for internal consistency and reliability, such as stereochemical quality, using PROCHECK (Laskowski et al., 1996). The quality of the non-bonded atom interactions was evaluated using ERRAT (Colovos & Yeates, 1993).

Docking studies

The validated Schgr-AKHR model and the structure of the DPC micellar structure of the three AKH peptides were prepared for docking simulations using Protein Preparation Wizard and LigPrep of the Schrödinger suite (Schrödinger Inc., New York, NY, USA). Site-directed mutagenesis studies (Kooistra et al., 2013), molecular modeling and structural analyses (Li, Hou & Goddard, 2010) suggest that most of the class A GPCRs share a similar binding pocket c.f. retinal bound to rhodopsin, carazolol bound to beta-2 adrenergic receptor (Vilar et al., 2011) and Anoga-HrTH bound to the AKHR of Anopheles gambiae (Mugumbate, Jackson & Van Der Spoel, 2011). Thus, the extracellular portion of the receptor was used for the docking simulations.

GLIDE docking (Trott & Olson, 2010) was used for peptide docking with a grid space of 72 × 72 × 72, which covered all extracellular loops and helices. The receptor grid was generated for peptide ligands and the docking precision was SP-Peptides. This setting automatically increases the number of poses collected.

MD of docked structure

The best poses from the docking studies were used as starting structures for a 2 μs MD simulation in a 1-palmitoyl-2-oleoyl-glycero-3-phosphocholine (POPC) membrane. Using the CHARMM-GUI (http://www.charmm-gui.org/) the docked complex was placed in a POPC membrane (128 POPC molecules) such that it spanned the membrane. The construct was then converted to an OPLS-AA, all–atom, force field. Using GROMACS, ~12,000 water molecules were added and the charge neutralized by adding 19 Cl⁻ ions. Several steps of equilibration were used, to pack the membrane around the receptor. This was followed by 1 μs of NPT simulation at 300 K with Berendsen pressure coupling (Berendsen et al., 1984) and a tau-p of 2.0. The free energy of binding of the final structures, from the dynamic simulations, was calculated using Prime-MM-GBSA (Schrödinger Inc., New York, NY, USA).

Spectral assignment

The chemical shift assignments of the three AKH ligands in DPC are given in Tables 2–4. Previously, Berjanskii & Wishart (2006) as well as Tremblay, Banks & Rainey (2010) have shown that the chemical shifts of peptides are sensitive to the secondary structure and flexibility of the polypeptide. Structuring-induced chemical shift changes (observed shifts minus random coil reference values) were analyzed using the CSDb algorithm available at https://andersenlab.chem.washington.edu/CSDb/ (Eidenschink et al., 2009; Fesinmeyer, Hudson & Andersen, 2004). Figures 1A–1C show such plots for the three ligands. For Schgr-AKH-II both H^N and H_α are shifted up-field, while for Locmi-AKH-I and Aedae-AKH, only the H^N deviations are up-field. The H_α deviations are small and random. It has been shown that polypeptides with a helical structure have, H^N and H_α chemical shifts which are, on average, −0.30 ppm less than their random coil values, while β-sheet structures have shifts of ca. 0.6 ppm (Szilágyi, 1995).

Thus, each of these peptides has a β-turn structure. Similar results were found for a number of other deca- and octapeptidic members of the AKH family, i.e., for Declu-CC, Melme-CC, and Dappu-RPCH (Jackson et al., 2018). However, downfield shifts were previously found for Anoga-HrTH from the Anopheles mosquito (Mugumbate, Jackson & Van Der Spoel, 2011; Mugumbate et al., 2013).

Using the Random Coil Index tool (Tremblay, Banks & Rainey, 2010), the chemical shifts were also used to estimate the model-free order parameter, S², of the peptides (see Fig. 1D). An order parameter of 1 means the peptide is rigid, while an order parameter of 0 means the peptide has no structure. Figure 1D shows that Locmi-AKH-I is very ordered, with a maximum order parameter of 0.9 around proline, whereas the C-terminal has less ordering (S² = 0.30). On the other hand, Schgr-AKH-II and Aedae-AKH are much more flexible. The order parameter for these two peptides range from 0.1 to 0.4, which is similar to that of Dappu-RPCH, an AKH peptide member from the crustacean water flea (Jackson et al., 2018).

MD simulation with DPC micelle

Figure 2 shows the solution structures of the three AKH ligands in DPC micelle solution. In each case, the MD was started with the peptide in water, but they rapidly diffused to interact with the DPC micelle. Depending on the starting orientation of the peptide relative to the micelle, the peptide would make contact with the phospholipid and move away until a stable orientation was established. This is shown in Fig. 3 where the peptide/DPC contact area is plotted as a function of time. For Locmi-AKH-I, contact between the DPC and micelle is established during the equilibration period, for Aedae-AKH it is established after 30 ns, while for Schgr-AKH-II, even after 60 ns the peptide is still not permanently attached to the micelle. It is interesting to note that the contact area between Locmi-AKH-I is much higher than that for the other two peptides and Locmi-AKH-I is much more rigid even though it is longer, a decapeptide versus an octapeptide. The interaction between the peptides and the lipid surface, as shown by the contact area, is important as it has been postulated that, before the ligand binds to its receptor, it first binds to the cell membrane surface. Thus, surface binding is an important step in receptor activation.

In order to show how the structure of the peptides changed during the MD, an overlay of each peptide is shown in Figs. 2A–2C. Cluster analysis of the trajectory gave a single large cluster for each peptide with a number of smaller, high energy, clusters. The root conformer and an overlay of each cluster is shown in Fig. 2. The predominant conformation of each peptide does have a turn feature but the details differ for each AKH. Locmi-AKH-I (Fig. 2A) has a clear β-turn around its proline residue; Aedae-AKH (Fig. 2B) has a more open structure compared to the other peptides and no marked turn around proline; Schgr-AKH-II (Fig. 2C) is tightly coiled in DPC solution.

Receptor construct

Use of the GPCR-ModSim server gave the crystal structure of the hĸ-OR (PDB ID: 4DJH), with 2.9 Å resolution, as a top template for Schgr-AKHR. The crystal structure has the selective antagonist JDTic complexed to the receptor (Wu et al., 2012). This template has the highest sequence identity (26.33%) compared to other templates as shown in Fig. S2. hĸ-OR belongs to the class A (rhodopsin-like), γ-subfamily of GPCRs and was selected to build models for Schgr-AKHR. The initial sequence alignment from the GPCR-ModSim server was aligned manually with the use of Chimera (Pettersen et al., 2004). The predicted TM helices of Schgr-AKHR and the PDB structural assignments of hĸ-OR were also used to confirm the alignment as shown in Fig. S1. The sequence analysis shows that the conserved residues in the seven TM helices (TM1–TM7) of hĸ-OR are also highly conserved within Schgr-AKHR (conserved residues highlighted by purple colored boxes in Fig. S1), indicating that they may be involved in functions, such as signaling and ligand binding, of Schgr-AKHR. Additionally, as in hĸ-OR, disulfide forming cysteine residues (Cys-131 and Cys-210) are conserved in Schgr-AKHR. The predicted TM helices of Schgr-AKHR (Fig. S3) are consistent with the TM helices of hĸ-OR: there are no gaps or insertions in these regions, signifying that the target sequence is correctly aligned with the template sequence and, hence, can be used for the modeling process.

The homology models of Schgr-AKHR were built using the Modeler 9v7 program. The input parameters of Modeler were set to generate 100 models with high structural optimization. A disulfide bridge, Cys-131–Cys-210, was defined as in the template structure. The best model was selected based on the lowest PDF (molecular probability density function) energies and DOPE score (discrete optimized protein energy) for the docking simulations. The selected models' qualities were subsequently assessed with structural evaluation programs such as PROCHECK and ERRAT. Ramachandran plot analyses of the Schgr-AKHR model from PROCHECK are shown in Fig. S4. All the residues were either in favored or in allowed regions, indicating that the backbone torsion angles (phi and psi) of this model are reliable. In addition, the ERRAT score, so-called overall quality factor, was computed on the Schgr-AKHR model to check the quality of its non-bonded atomic interactions. The normally accepted score range for a high-quality model is >50 (Colovos & Yeates, 1993). The ERRAT score for the Schgr-AKHR model was 78, showing that the model is within the high-quality range. All these validation methods demonstrated that the model is reliable and can be used for further studies. The final 3D structure of the Schgr-AKHR model, as shown in Fig. S5A, is similar in overall fold with the template protein hĸ-OR (Fig. S5B). The superposition of the Schgr-AKHR model with that of hĸ-OR displayed a 2.804 Å root mean squared deviation (RMSD) for 1,332 atoms pairs. This low RMSD demonstrates that the overall tertiary structure of the model is similar to the template structure. In addition, the superimpositions indicate that the seven TM helices are highly conserved with hĸ-OR. However, there are slight structural variations in the loop regions of the model.

Ligand docking

LigPrep was used to generate multiple conformers of the three AKH peptides, which were then docked to the Schgr-AKHR model. A receptor grid was generated for the extra-cellular half of the GPCR and the peptide docked using SP-Peptide precision. A total of 100 different poses were collected and scored for each peptide. An overlay of the highest scoring poses gave the same receptor binding site for all three peptides. This binding pocket consisted of a cleft running across the top of Schgr-AKHR, between helices 2, 6, and 7 and extra-cellular loops 2 and 4 (details shown in Figs. 4–6). The peptides lay along this cleft. It was found that, while the GLIDE protocol tried to dock many different conformations of the peptide, only those with some turn structure were successful. Many of the docked poses had similar GLIDE scores but different sets of peptide/receptor interactions. During the docking, it was found that the orientation of the peptide within the binding pocket did not change. For this reason, the docking was repeated with the peptide rotated through 180° around the axis of the TM helices, i.e., the direction of N-terminus to the C-terminus was reversed. Again, the peptide could be docked successfully, indicating that the binding pocket was quite promiscuous. Although either orientation of the peptide would bind to the receptor, the binding energies of the two orientations differed by some 50 kcal/mol and hence the original orientation was chosen for further study.

The docked structure of each AKH peptide with the highest binding energy, was used as the starting structure for MD of the complex in a POPC membrane. During the dynamics, the ligands were found to move, and individual side-chains rotate within the binding site, making and breaking H-bonds to various residues of the receptor. This was essentially the same as what was found during the GLIDE docking. While the POPC membrane added to the computational cost of the simulation, it was necessary to prevent the transmembrane helices of the receptor from moving apart. Snapshots of the simulation were transferred to Maestro, where Prime-MM-GBSA was used to calculate the binding energy.

Figures 4A and 4B shows an overlay of several snapshots of the dynamic simulation of the octapeptide, Schgr-AKH-II, in its receptor binding pocket. As can be seen, the overall conformation of Schgr-AKH-II remains the same (Fig. 4C) but the side-chains move within the binding site, forming and breaking intra- and inter-molecular H-bonds. There is also some movement of the receptor during the dynamics. The free energy of binding of the different snapshots were not significantly different and ranged from −94 to −116 kcal/mol over the 1 μs simulation.

The bound conformation of Schgr-AKH-II is shown in Fig. 4C, while Fig. 4D is an overlay of bound Schgr-AKH-II and its lowest energy conformation in DPC micelle solution. The agreement of these two conformers is remarkable, especially considering that the GLIDE—SP protocol generates some 100 different starting conformations for the peptide docking.

Figure 5A shows the bound conformation of Locmi-AKH-I, while Fig. 5B is an overlay of this bound conformer and the lowest energy conformer found in DPC micelle solution. The agreement here is not as close as that for Schgr-AKH-II, but the same turn structure is seen. Figure 5C shows Locmi-AKH-I in the receptor binding pocket, while Fig. 5D and 5E show the details of how Locmi-AKH-I fits into the receptor: the decapeptide stretches across the cleft in the receptor with the central portion of the peptide fitting into the binding pocket, but the two termini pointing outside the binding pocket. Locmi-AKH-I gave the most trouble during the docking stage as poses were frequently rejected. During the MD, the terminal amide of this decapeptide was sometimes found to H-bond to a POPC molecule, which is of course not present during the GLIDE docking. The final binding energy for Locmi-AKH-I was −98 kcal/mol.

Figures 6A and 6B show how Aedae-AKH fits into the binding pocket of the Schgr-AKHR model. The arrangement is similar to that of Locmi-AKH-I, except, in this case, the termini do not extend outside the receptor. In the case of Aedae-AKH, ECL4 folds over the top of the binding site, trapping the peptide inside. The conformation of bound Aedae-AKH is shown in Fig. 6C, and an overlay with the DPC micelle solution conformation is shown in Fig. 6D. Again, these two conformations are very similar, supporting the idea that the peptides are pre-arranged on the cell surface. The binding energy of Aedae-AKH was −88 kcal/mol.

Ligand interaction diagrams for the three ligands are shown in Fig. 7, while Table 5 lists the interactions between the ligands and Schgr-AKHR. From these data, it is clear that all three AKH ligands have very similar interactions with Schgr-AKHR. Both Schgr-AKH-II and Aedae-AKH, H-bond to His169 of the receptor, while W⁸ of Locmi-AKH-I π-stacks with this residue. In Locmi-AKH-I, the amide carbonyl, pE¹CO, H-bonds to Lys281, while in Aedae-AKH it is pE¹O_ɛ1 which H-bonds to Lys281. The terminal residue of both Schgr-AKH-II and Locmi-AKH-I H-bond to Lys288.

Analysis of molecular switches

A feature of class A GPCRs is the presence of highly conserved molecular switch motifs. These switches, which play key roles in the stabilization of the receptor in an inactive and active state, include a TM3–6 lock, a SPLF switch, a tyrosine toggle, and a DRY ionic lock. The breaking of these switches results in movement of the TM helices, which can activate the receptor (Trzaskowski et al., 2012). These switches are the same as those reported for the AKHR of Drosophila melanogaster, Tribolium castaneum, Anopheles gambiae, and Rhodnius prolixus (Rasmussen et al., 2015), suggesting that the activation mechanism of Schgr-AKHR may be the same.

The DRY ionic lock between arginine and tyrosine is postulated to open and close during receptor activation. This is shown in Fig. 8. In the inactive state, the DRY switch is closed (Fig. 8B) but upon ligand binding, TM6 and TM3 twist, opening this switch (Fig. 8C).

The TM3–6 lock involves two residues in the binding pocket, Arg¹⁰⁷ on TM3 and Tyr²⁶⁵ on TM6 (Fig. 9A). In the inactive state these two residues are far apart, but on ligand activation, they move closer together. Figure 9A also shows, in the active state, Arg¹⁰⁷ H-bonding with Glu¹⁹⁰ of ECL4. It was noted before that this loop closes over the binding pocket after ligand binding. In Fig. 9A one can also clearly see how TM6 and TM3 move together on the extra-cellular side but move away from each other on the intracellular side, upon activation.

The tyrosine toggle switch involves the NPxxY motif on TM7 and Tyr²¹³ on TM5 (Fig. 9B).