Characterization, Docking and Molecular Dynamics Simulation of Gonadotropin-Inhibitory Hormone Receptor (GnIHR2) in Labeo Catla
Pravesh Kumara Mukesh Kumarb K. S. Wisdomb Gireesh-Babu Pathakotab
Sunil Kumar Nayakb Dhalongsaih Reangb N. S. Nagpureb Rupam Sharmab
aDepartment of Aquaculture, College of Fisheries, Dr. Rajendra Prasad Central Agricultural University, Pusa, Bihar, India, bICAR-Central Institute of Fisheries Education, Mumbai, India
Key Words
Labeo catla • Reproduction • RF313 • Molecular docking • MD simulation
Abstract
Background/Aims: GnIH receptors (GnIHRs) belong to the family of G-protein coupled receptors (GPCRs) and play a key role in the regulation of reproduction from fishes to mammals, either by inhibiting or stimulating the expression of gonadotropins. The aim of this study was to characterize GnIH receptor (GnIHR2) from Indian Major Carp, Labeo catla and its docking and simulation with GnIH antagonist RF313. Methods: The full length sequence of GnIHR2 was obtained with RACE PCR. The docking analysis of RF313 with GnIHR2 receptor was performed with AutoDock v. 4.2.6 and molecular dynamics (MD) simulation with GROMACS 5.0. Results: In the present study, we cloned full-length cDNA (1733 bp) of GnIHR2 from the brain of L. catla. The phylogenetic analysis showed clustering of catla GnIHR2 with goldfish and zebrafish in the GPR147 group. L. catla GnIHR2 receptor comprised seven transmembrane domains and the 3D-structure was predicted by I-TASSER tool. The docking analysis revealed high binding affinity (-11.6 kcal/mol) of GnIH antagonist, RF313 towards GnIHR2 receptor. The primary bonds involved were alkyl and hydrogen bonds while the amino acids participated were proline 43, 210, 339, cysteine 214, leucine 211, serine 213 and phenylalanine 338. The MD simulation analysis of docked complex for 100 nano-seconds (ns) in the lipid membrane environment showed the stability of the complex with time. Conclusion: Our study showed that GnIH antagonist, RF313 interact tightly with the GnIH receptor, GnIHR2 of L. catla. To the best of our knowledge, this is the first report on computational modelling and MD simulation of GnIH receptor in fishes. This will help in functional characterization studies of GnIH/GnIHR system in vertebrates.
Introduction
Gonadotropin-inhibitory hormone (GnIH), a hypothalamic neuropeptide, is the negative regulator of hypothalamus-pituitary-gonadal (HPG) axis in birds and mammals [1-2]. However, in fish, GnIH is reported to have negative and positive effects on reproduction depending upon the species, sex and stage of the reproductive cycle [3-10]. Recently, we have characterized the full-length cDNA sequence of the GnIH gene from the brain of Indian major carp, Labeo catla and showed that it stimulates the expression of gonadotropins and other essential genes in the brain [10]. To elucidate the molecular mechanism of action of GnIH on gonadotropin synthesis and release, it is essential to identify and characterize the GnIH receptor/s.
GnIH receptors (GnIHRs) are G-protein coupled receptors (GPCRs) with seven transmembrane domains, designated as GPR147 [11]. Bonini et al. [12] first identified two GPCRs for neuropeptide FF (NPFF) and labeled them as NPFF1/GPR147 and NPFF2/GPR74. GnIH and NPFF are paralogs with conserved LPXRFamide and PQRFamide motifs at C-terminal, respectively. NPFF is a pain modulatory neuropeptide, while GnIH has a role in reproduction [13]. GnIH showed higher affinity for the GPR147 receptor, whereas NPFF had higher activity for the GPR74 receptor [12, 14]. Further, Ikemoto and Park [15] also showed that chicken GnIH peptide inhibited the Gαi2 mRNA expression about 100-fold stronger in COS-7 cells transfected with GPR147 than GPR74. Also, in quail, the GnIH and GnIH related peptides bound specifically with GPR147 in transfected COS-7 cells [16]. These results suggest that GPR147 is the primary receptor for GnIH/LPXRF gene.
In birds and mammals only one type of GnIH receptor has been identified [15-19], while in fishes like goldfish, zebrafish and common carp up to three GnIH receptors (GnIHRs) were detected [8, 20-21]. However, in tilapia [3] and grouper [6], only one GnIH receptor was identified, indicating the structural and functional variability of GnIH/GnIHR system in fish. GnIH/GnIHRs may be involved in regulating HPG axis at all levels, as GnIHRs were observed from the hypothalamus to gonads on the whole axis [22, 23].
Selective receptor antagonists are generally used to block the ligand receptor interaction and study its downstream effect. Simonin et al. [24] developed an RFamide dipeptide, RF9, which showed antagonistic activity towards both GPR147 and GPR74 receptors. Intracerebral injection of RF9 induced robust LH surge in rats and mice; while in the ewe, it increased both LH and FSH secretion during the oestrus as well as anoestrus season [25-28]. However, it was later revealed that the RF9 is a kisspeptin agonist which increases the gonadotropin levels by activating the kisspeptin receptors and not by blocking the GnIH receptors [29-32]. Recently, Elhabazi et al. [33] developed a new GnIH antagonist, RF313, which displayed negligible affinity and no agonist activity (up to 100 μM) towards the kisspeptin receptor. In male hamster GnIH peptide exerts a stimulatory impact on gonadotropin release and it was showed that RF313 could block the LH surge entirely. However, there is no report on docking and molecular dynamic (MD) simulations of GnIHRs with agonists/antagonists in animals.
Materials and Methods
Animals and tissue collection
For gene characterization, adult catla fish (average weight 2.0 kg) was sampled from ICAR-CIFE regional research centre at Powarkheda, Madhya Pradesh, India, and brain tissue was dissected and stored in RNAlaterTM solution (Qiagen, Germany). The dissection was carried out after anesthetizing the animals with clove oil. Experimental procedures were conducted following the guidelines of the CPCSEA (Committee for Control and Supervision of Experiments on Animals), Ministry of Environment, Forest and Climate Change, Government of India on care and management of animals in scientific research.
Total RNA isolation and cDNA synthesis
Total RNA was isolated from 100 mg brain tissue using Trizol™ reagent (Invitrogen, USA) as per manufacturer's instructions. The integrity of the isolated RNA was checked on 1.5% agarose gel, while the purity and quantity of RNA were measured by Nanodrop 2000/2000c (ThermoFisher Scientific, USA). It was then treated with DNase I (ThermoFisher Scientific, USA) to remove the contamination of genomic DNA. First-strand cDNA for degenerate PCR was synthesized from 1 μg of total RNA using Affinity Script cDNA Synthesis Kit (Agilent Technologies, USA). For RACE-PCR, the first stand cDNA was synthesized using SMARTer® RACE 5’/3’ Kit (Clontech, USA) as per manufacturer’s instructions. The synthesized first-strand cDNAs were diluted and stored at −80 °C until use.
Molecular cloning of GnIHR2 receptor
A set of degenerate primers (DegF/R) for the partial amplification of catla GnIHR2 (Table 1) was designed based on the conserved amino acid (AA) motifs in closely related fish species. For this, GnIHR amino acid sequences of closely related fish species reported in NCBI GenBank database were obtained and multiple sequence alignment was performed to identify conserved domains. The degenerate primers thus designed were synthesized commercially by Eurofins Genomics Pvt. Ltd., Bangalore. PCR amplification of brain cDNA was performed in a 96-well Takara PCR System (Takara, Japan) using the following conditions: initial denaturation at 94 °C for 4 min, followed by 35 cycles of 94 °C for 30 sec, 62°C for 30 sec and 72 °C for 45 sec and a final extension of 7 min at 72 °C. The amplified cDNA fragments were resolved on 1.5% agarose gel, purified, cloned into pTZ57R/T vector (ThermoFisher Scientific, USA), sequenced and confirmed by NCBI BLAST.
To obtain the full-length cDNA sequence of GnIHR2, 5’ and 3’ RACE PCR was performed using gene-specific primers (Table 1) that were designed from the partial cDNA sequence. For the amplification of 5’ and 3’ cDNA ends, two PCR reactions were performed. The first one was touchdown PCR using 5’ and 3’ GSP1 primers (Table 1) and 10× universal primer mix (UPM) with following conditions: 94 °C (30 s), 72 °C (2 min) for 5 cycles, 94 °C (30 s), 70 °C (2 min) and 72 °C (2 min) for 5 cycles, 94 °C (30 s), 68 °C (30 s) and 72 °C (2 min) for 25 cycles. Following the first PCR reaction, nested PCR was performed on diluted first PCR product using 5’ and 3’ GSP2 primers (Table 1) and short universal primer. The nested PCR conditions were as follows: 94 °C (30 s), 65 °C (30 s), and 72 °C (1 min) for 35 cycles. The amplified PCR products were cloned in pTZ57R/T cloning vector (Thermo Scientific, USA), sequenced and confirmed by NCBI BLAST analysis. Full length cDNA sequence was obtained by aligning the partial sequences.
Sequence analysis of GnIHR2 receptor
The cDNA sequence of catla GnIHR2 was analyzed using the BLAST algorithm at NCBI (http://www.ncbi.nlm.nih.gov/blast). Open reading frame (ORF) analysis was performed using NCBI ORF finder (http://www.ncbi.nlm.nih.gov/gorf/gorf.html), and the nucleotide sequence of the ORF was translated into the respective AA sequence using ExPASY translation tool of EMBL (http://web.expasy.org/translate/). Multiple sequence alignment of the deduced AA sequence was carried out using Clustal Omega online software (https://www.ebi.ac.uk/Tools/msa/clustalo/). Calculated molecular weight (MW) and predicted isoelectric point was obtained by using the ExPASy online server (https://web.expasy.org/protparam/). The protein-protein interaction of the GnIHR2 protein with other proteins involved in the reproductive pathway was determined using online software STRING (http://string-db.org). The secondary structure of protein was predicted by SOPMA online server (https://npsa-prabi.ibcp.fr/cgi-bin/npsa_automat.pl?page=/NPSA/npsa_sopma.html). Phylogenetic analysis of GnIHR2 was conducted with MEGA 7.0 using the neighbor-joining method with bootstrap values of 1000.
Homology modeling and energy minimization of GnIHR2 receptor
For structural modeling of GnIHR2 receptor, first we used SWISS-MODEL, PHYRE2 Protein Fold Recognition Server and CABS-fold server, but due to very low homology with template proteins, the obtained models were not acceptable. After that, a threading-based GPCR- I-TASSER server was used [34]. This tool uses a computational method designed for 3D structure prediction of GPCRs. This is a hierarchical approach for protein structure and function prediction, which identifies structural templates from the PDB by multiple threading approache, with full-length atomic models constructed by iterative template fragment assembly simulations. For each receptor, ten models with different conformations of the loops were built. The best model was chosen based on C-score. The selected model was further subjected to energy minimization using the GROMACS 2018.1 package [35-36]. In the first step, the structure was minimized using steepest descent method for about 50000 steps. The conjugate gradient method of minimization was used to get the structure in a globally minimal energy state. RAMPAGE (http://mordred.bioc.cam.ac.uk/~ rapper/rampage.php), an online server, was used for Ramachandran plot development from the predicted model, which evaluates the quality of final GnIHR2 model.
Molecular docking of RF313 ligand with GnIHR2 receptor
Since the structure of GnIH antagonist, RF313 was not available in PubChem or any other databases, we designed the structure of RF313 in Avogadro molecule editor software. The structure was then converted to PDB file format by Open Babel GUI tool. AutoDock v. 4.2.6 was used for molecular docking of RF313 ligand with the GnIHR2 receptor to understand the interaction between them. In AutoDock software, the grid dimensions were set at 84 Å × 84 Å × 84 Å to cover the entire binding site of GnIHR2 with a spacing of 0.375 Å. The Gasteiger charges were added to the receptor and Lamarckian Genetic Algorithm was used as reported previously [37]. Total 5000 conformations were generated, which were then clustered and compared on the basis of binding energy, cluster size, intermolecular energy, electrostatic energy and van der Waals energy calculated by scoring function. The docked complexes were visualized using PyMOL [38]. The least energy conformation of RF313 with GnIHR2 was then selected for understanding the hydrogen-bond interactions and chosen as a starting structure for MD simulation in the lipid bilayer environment using GROMACS 2018.1.
MD simulation of complex in a lipid membrane
MD simulations were performed to know the stability of the GnIHR2:RF313 docked complexes using the CHARMM36m force field [39]. The CHARMM-GUI was used to embed GnIHR2-ligand complexes in a POPC lipid bilayer membrane as GnIHR2 is associated with the GPCR [40]. This receptor peptide complex embedded in POPC lipid bilayer was neutralized by counter ions, either Na+ or Cl-. After that, the complex was subjected to energy minimization using the steepest descent method followed by conjugate gradient method for 50000 steps each. By applying position restraint 500 ps simulations were performed using the NVT and NPT ensemble. Berendsen thermostat was used to maintain a constant temperature at 303.15 K and MD simulation was performed for 100 nano-seconds (ns) [41]. The constraint on the H-bond lengths is applied using the LINCS algorithm [42]. The binding energy was calculated using the ‘g_mmpbsa’ tool of GROMACS 5.0 [43], to understand the binding affinity of GnIHR2 with RF313. Here, we used the entire MD simulated trajectory (0 to100 ns) to calculate binding energy between the ligand and receptor. In this study, the contribution of entropy in the binding energy calculation was not considered as reported by other authors [44-45].
Ethics Statement
The care and treatment of animals used in this study were conducted based on the guidelines of the CPCSEA (Committee for Control and Supervision of Experiments on Animals), Ministry of Environment, Forest and Climate Change (Animal Welfare Division), Government of India on care and management of animals in scientific research. The Board of Studies and Authorities (BOSA) Fish Genetics and Biotechnology Division of ICAR- Central Institute of Fisheries Education, Mumbai, has approved the study.
Results
Cloning and characterization of GnIH receptor (GnIHR2)
In the present study, we cloned full-length cDNA of GnIHR2 containing 1733 nucleotides from the brain of L. catla by 5' and 3' RACE-PCR. The cDNA consisted of 1452 bp ORF that encoded a precursor protein of 483 AA with 5' and 3' untranslated region (UTR) of 65 bp and 216 bp, respectively. Predicted molecular weight and isoelectric point of the GnIHR2 have been recorded as 54.41 kDa and 8.74, respectively. The initiation codon (ATG) and stop codon (TAA) were identified at positions 66 bp and 1515 bp, respectively (Fig. 1). The sequence was submitted to NCBI GenBank (Acc. No. LC380835).
Sequence analysis of GnIHR2 gene
The NCBI protein BLAST of GnIHR2 AA sequence showed high similarity with Carassius auratus (91%) and Danio rerio (85%) followed by Astyanax mexicanus (79%), Scleropages formosus (74%), Seriola lalandi (72%) and Cynoglossus semilaevis (56%) (Supplementary Fig. 1). Further analysis of the predicted amino acid sequences of the GnIHR2 revealed seven putative transmembrane domains (TMD) (Fig. 2). The multiple sequence alignment of catla GnIHR2 AA sequence with other vertebrate species is presented in Fig. 3.
The secondary structure prediction of GnIHR2 showed 42% alpha helix and 14% extranded strand (Table 2). The phylogenetic analysis using the Neighbor-joining method showed that Catla GnIHR2 clustered with goldfish and zebrafish in the GPR147 group. All three species belong to the cyprinid group and formed a separate clade. The GPR147 of humans, mammals and birds formed different clades. Further, mammals, birds, and fish NPFFR2/GPR74 clustered in the GPR74 group (Fig. 4).
Homology modeling of GnIHR2
Out of ten models generated by GPCR I-Tasser server for catla GnIHR2, the best model with maximum C-score (-1.94) was selected for further analysis (Table 3). The structure of the selected model is represented schematically in Fig. 5 with seven transmembrane helices like other known GPCR structures. Ramachandran plot analysis showed that 64.4% (278) residues were present in the most favored region, while 25.7 % (111) residues were present in the additional allowed region and rest of the residues in general allowed and disallowed regions (Fig. 6).
Molecular docking of RF313 with GnIHR2
The molecular docking was performed to study the detailed molecular basis of interactions and to estimate the binding affinity of ligand, RF313 with the GnIHR2 receptor. The predicted structure of RF313 was docked into the active site of GnIHR2 and the results revealed high binding affinity of RF313 antagonist towards GnIHR2 protein with a binding energy of -11.6 kcal/mol. The 3D structures showing the interactions between the protein and ligand in different poses are shown in Fig. 7. Different types of bonds like chemical covalent, non-covalent, hydrogen and van der waals forces are found to be involved in the protein-ligand binding. The amino acids - leucine (LEU) at position 211 and phenylalanine (PHE) at 338 of GnIHR2 form the chemical covalent bond with ligand RF313, while proline (PRO) at positions 43, 210, 339 and cysteine (CYS) at 214 form the non-covalent bond. Hydrogen bond was formed by serine (SER) present at position 213 of GnIHR2 with RF313 (Fig. 8). Apart from these direct bonds, many other hydrophobic and van der waals forces are thought to play an essential role in stabilizing the complex.
MD Simulation of GnIHR2 and RF313
The simulation study was performed to check the interactions between ligand and receptor in the simulated lipid bilayer environment and the effect of ligand binding on conformational changes of receptor. We performed MD simulations of up to 100 nanoseconds, which is sufficient time for side-chain rearrangement and keeping the complex in the most stable binding confirmation. The snapshots of the complex during the entire simulation time are shown in Fig. 9, which depicts that the complex is stable and ligand binds with the receptor in the inner side. The stability of the GnIHR2:RF313 complex was assessed using all-atom explicit MD simulations. Root Mean Square deviation (RMSD) values plotted over the simulation time revealed the stable dynamics and equilibration after 60 ns for GnIHR2 protein and GnIHR2:RF313 complex, while ligand RF313 showed stable dynamics just after 30 ns. The stabilization of GnIHR2:RF313 complex at higher RMSD value of ~8Å (Fig. 10) suggests that this complex undergoes significant conformational changes during the simulation promoting tight binding between RF313 ligand and GnIHR2 receptor. We performed the Root Mean Square Fluctuation (RMSF), Radius of Gyration (Rg) and Solvent-accessible surface area (SASA) analysis for the GnIHR2-RF313 complex and results are depicted in Fig. 11. The RMSF value varies from 0.1 to 0.5 nm except for C and N terminals. Both the terminal ends of the protein, which lack native secondary structure showed an increased RMSF value of >1.5 nm. The overall residual fluctuation is less than 0.4 nm, which also revealed the stability of this protein during the simulation (Fig. 11A). The Rg is used to predict the compactness of protein structure. In this study, we plotted the Rg graph with simulated time for GnIHR2:RF313 complex, and the result suggests that it forms a compact globular shape after the equilibration period of 50 ns. Rg values first increased from 2.9 to 3.25 nm; after that it showed declining trend indicating the proper folding of the complexes to adopt compact globular shape during the simulation (Fig. 11B). SASA is the surface area of a biomolecule that is accessible to a solvent. The surface area of GnIHR2:RF313 complex increased from 265 to 300 nm2 in first 0 to 50 ns simulation time period, while in rest of simulation time the surface area was first decreased to 280 nm2 followed by increase to 290 nm2 by the end (Fig. 11C).
The binding energy calculated using the MM-PBSA tool of GROMACS 5.0 for the RF313:GnIHR2 complex was -181.08±60.41 kcal/mol (Table 4). We also plotted per residues interaction energies for the GnIHR2:RF313 complex (Fig. 12). Per residue decomposition energy provides an insight into the role of individual residues contribution in the binding affinity between the GnIHR2 and ligand.
Discussion
In fishes, GnIH receptors have been cloned from very few species [3, 6, 8, 20-21]. In the present study, we cloned full-length cDNA of GnIHR2 from Indian major carp, L. catla. Three GnIH receptors have been reported in common carp, goldfish and Zebrafish [8, 20-21]. However, in birds, humans and some other fish species like tilapia and orange-spotted grouper only one GnIH receptor has been detected [3, 6, 16, 19, 23]. We identified three GnIH receptors in catla and in this study the sequencing and analysis of only one GnIH receptor, GnIHR2 was performed. The other two receptors for GnIH; GnIHR1 and GnIHR3 are partially characterized and further work is going on in our laboratory (data not provided). The GnIHR2 cDNA comprised an ORF of 1452 bp, encoding a polypeptide of 483 AA. This is in concurrence with the study by Wang et al. [6] who reported an ORF of the same length encoding a protein of 483 AA for GnIHR in orange-spotted grouper. However, in the closely related species like goldfish and zebrafish GnIHR2 polypeptide of 484 amino acids were reported [20-21]. The full-length GnIHR2 receptor contained a seven-transmembrane domain, an extracellular N-terminus and a cytoplasmic C-terminus. This characteristic feature of transmembrane domain in GPCRs was detected in different studies [6, 8, 11]. In the phylogenetic analysis, GnIHR2 of L. catla and GPR147 of other vertebrates are clustered into a single clade discrete from GPR74. Zhang et al. [20] and Qi et al. [21] reported similar results in zebrafish and goldfish, respectively.
In this study, GPCR I-Tasser server was used to generate the 3D model of catla GnIHR2 protein due to low sequence similarity with template proteins available in PDB database. The best model was selected based on the c-score, stereochemical quality and the accuracy of the predicted protein models. The model's quality was checked using Ramachandran plot. The low percentage of residues in the most favored region may be due to the non-availability of template structure in the PDB database. However, this does not affect further analysis as the amino acids involved in the interaction with RF313 (Fig. 8) were present in the most favored region of the Ramachandran plot. Also, the results of MD simulation obtained in this study confirmed that the predicted model is stable with time. Kauffman et al. [46] showed that templates of low sequence similarity with the query sequence could also be used for developing a good model for docking analysis when the sequence alignment was handled correctly. Different studies for GPCRs also used the c-score and Ramachandran plot for analyzing the predicted protein models [47-51].
The docking results showed high binding affinity of GnIH antagonist, RF313 towards GnIHR2 receptor with the binding score of -11.6 kcal/mol. The 3D structure of complex showed the interactions between the receptor and ligand in different orientations. The primary bonds involved in the binding of ligand with receptor were hydrogen, pi-alkyl, pi-sigma bonds and Van der Waals (VdW) attractions. Previously, Rather et al. [47] docked the kisspeptin with kiss1 receptor, a GPCR receptor and reported a high binding affinity of -11.3 Kcal/mol. The interactions of hydrogen bonds, hydrophobic and VdW forces played an essential role in stabilizing the complex. Vijayakumar et al. [51] also studied the docking of Cannabinoid receptor 2 (CB2) protein with twenty-four (18 bioactive molecules, five agonists and one antagonist) small molecules and showed the docking score of -7.9 and above. The interactions were due to the hydrogen bond backbone and side-chain, hydrophobic interactions and pi-pi stacking. Manogar et al. [49] docked the eight cyanobacterial ligand molecules with cannabinoid receptor 1 (CNR1), a GPCR receptor with the help of the Schrodinger tool. The results showed a high docking score of -7.35 Kcal/mol for Symplocamide A ligand. The present study showed an excellent docking score (-11.3 Kcal/mol) compared to above studies, depicting the stability of the complex.
MD simulations is a technique that helps to understand molecular structure and biological interactions of proteins and macromolecules. In this study, the stability and conformational changes of GnIHR2 protein in the presence of RF313 ligand were investigated by performing MD simulations for upto 100 nanoseconds in lipid bilayer membrane enivironment. Simulation time of 100 nanoseconds used in the present study is adequate for the side chain reorganizations in order to facilitate the most stable binding conformation [47, 52]. Even, Feng et al. [53] showed that MD simulation for 20 nanoseconds between β2 adrenergic receptor and Gs protein results in stable conformational changes from the extracellular region to the intracellular region. The phospholipids which are most commonly used in membrane simulation studies are DLPC (1,2-dilaureoyl-sn-phosphatidylcholine), DMPC (1,2-dimyristoyl-sn-phosphatidylcholine), DPPC (1,2-dipalmitoylsn phosphatidylcholine), DOPC (1,2-dioleoyl-sn phosphatidylcholine), POPC (1-palmitoyl-2-oleoyl-sn-phosphatidylcholine), and POPE (1-palmitoyl-2-oleoyl-sn phosphatidylethanolamine). The large variety of possible acyl chain types of phospholipids in cell membrane would result in very complex and bulky systems for simulation, so to overcome this problem researchers generally use a palmitoyl and an oleyl chain in each phospholipid [54]. For the present study, POPC (1-palmitoyl-2-oleoyl-sn-phosphatidylcholine) has been used and this is the membrane of choice for several researchers in the previous studies [55-57].
RMSD values showed stable dynamics after 60 ns for GnIHR2 protein and GnIHR2:RF313 complex at ~8Å, while ligand RF313 showed stable dynamics just after 30 ns. The RMSD value depends on the binding interaction of protein and ligand. The RMSD plot will depict the ratio of fluctuation in residue level. Hence RMSD plots are important to transfer out the prediction of structural stability of protein [58-59]. Compactness of the protein structure can be predicted using radius of gyration (ROG). Proteins with high ROG are less tightly packed and vice versa [60]. Throughout the simulation time, GnIHR2 complex showed ROG in the range of 2.9 to 3.2 nm. The ROG of 3.2 nm was reached after 40 ns which remained nearly stable rest of the simulation period. The results suggest the stability of the complex without any major note worthy expansion/contraction. Binding free energy calculation permits the prediction of relative binding affinities between protein and peptides/ligands [61]. Similar to earlier studies, the entropic contribution in the binding energy was not considered in the present study, as its calculation is computationally costly [44, 62]. The binding free energy of complex is reasonably high (-181.08±60.41 kcal/mol) in this study, which showed that the complex is stable. The VdW (ΔEvdw) and electrostatic (ΔEele) interactions are important for the binding of protein-ligand complex. Here, electrostatic interactions make the highest contribution towards the binding free energy, while the polar solvation energy (ΔEpsol) is unfavorable for binding of ligand (Table 2). The total binding free energy is calculated by adding the ΔEvdw, ΔEele, ΔEpsol and SASA (Non-polar solvation) energies. Kumbhar et al. [62], calculated the binding free energy of indanocine with seven human tubulin isotypes αβI, αβIIa, αβIII, αβIVa, αβIVb, αβV, and αβVI and observed the highest (-50.70 kcal/mol) binding free energy with αβVI tubulin isotypes. Santoshi and Naik [63] used three ligands namely amino-noscapine, bromo-noscapine and noscapine for calculating the binding free energy with protein αβIII-tubulin isotype and found that amino-noscapine (-34.7 and -46.2 kcal/mol) bound more tightly than the other two ligands. Different researchers also showed that the high binding free-energy during the MD simulation increased the stability of the docked receptor-ligand complex [64-66].
Conclusion
In summary, the present study reports the full-length sequence of the GnIHR2 receptor from L. catla. The computational analysis showed seven transmembrane domains, and docking analysis revealed the high binding affinity of GnIH antagonist, RF313 with GnIHR2 receptor. The MD simulation analysis from RMSD, ROG, binding free-energy etc. showed the high stability of the GnIHR2:RF313 complex. The predicted ligand binding sites can be used for functional characterization of GnIH/GnIHR system in vertebrates. As per our knowledge, it is the first report of docking and MD simulation analysis of GnIHR2 in fishes.
Acknowledgements
This work is supported by Indian Council of Agricultural Research (ICAR), New Delhi, India.
Disclosure Statement
The authors declare they have no conflicts of interest.
References
|
1
Tsutsui K, Saigoh E, Ukena K, Teranishi H, Fujisawa Y, Kikuchi M, Ishii S,
Sharp PJ: A novel avian hypothalamic peptide inhibiting gonadotropin release.
Biochem Biophys Res Commun 2000;275:661-667. |
|
|
|
|
|
2 Yoshida H, Habata Y, Hosoya M, Kawamata Y,
Kitada C, Hinuma S: Molecular properties of endogenous RFamide-related
peptide-3 and its interaction with receptors. Biochim Biophys Acta
2003;1593:151-157. |
|
|
|
|
|
3 Biran J, Golan M, Mizrahi N, Ogawa S, Parhar
IS, Levavi-Sivan B: LPXRFa, the piscine ortholog of GnIH, and LPXRF receptor
positively regulate gonadotropin secretion in Tilapia (Oreochromis
niloticus). Endocrinology 2014;155:4391-4401. |
|
|
|
|
|
4 Moussavi M, Wlasichuk M, Chang JP, Habibi HR:
Seasonal effect of GnIH on gonadotrope functions in the pituitary of
goldfish. Mol Cell Endocrinol 2012;350:53-60. |
|
|
|
|
|
5 Moussavi M, Wlasichuk M, Chang JP, Habibi HR:
Seasonal effect of gonadotrophin inhibitory hormone on
gonadotrophin‐releasing hormone‐induced gonadotroph functions in
the goldfish pituitary. J Neuroendocrinol 2013;25:506-516. |
|
|
|
|
|
6 Wang Q, Qi X, Guo Y, Li S, Zhang Y, Liu X, Lin
H: Molecular identification of GnIH/GnIHR signal and its reproductive
function in protogynous hermaphroditic orange-spotted grouper (Epinephelus
coioides). Gen Comp Endocrinol 2015;216:9-23. |
|
|
|
|
|
7 Wang B, Yang G, Liu Q, Qin J, Xu Y, Li W, Liu
X, Shi B: Inhibitory action of tongue sole LPXRFa, the piscine ortholog of
gonadotropin-inhibitory hormone, on the signaling pathway induced by tongue
sole kisspeptin in COS-7 cells transfected with their cognate receptors.
Peptides 2017;95:62-67. |
|
|
|
|
|
8 Peng W, Cao M, Chen J, Li Y, Wang Y, Zhu Z, Hu
W: GnIH plays a negative role in regulating GtH expression in the common
carp, Cyprinus carpio L. Gen Comp Endocrinol 2016;235:18-28. |
|
|
|
|
|
9 Spicer OS, Zmora N, Wong TT, Golan M,
Levavi-Sivan B, Gothilf Y, Zohar Y: The gonadotropin-inhibitory hormone
(Lpxrfa) system's regulation of reproduction in the brain-pituitary axis of
the zebrafish (Danio rerio). Biol Reprod 2017;96:1031-1042. |
|
|
|
|
|
10 Kumar P, Wisdom KS, Bhat IA, Pathakota GB,
Nayak SK, Reang D, Nagpure NS, Sharma R: Molecular characterization of
gonadotropin-inhibitory hormone (GnIH) gene and effect of intramuscular
injection of GnIH peptide on the reproductive axis in Catla catla. Anim
Biotechnol 2020;31:335-349. |
|
|
|
|
|
11 Ubuka T, Son YL, Bentley GE, Millar RP,
Tsutsui K: Gonadotropin-inhibitory hormone (GnIH), GnIH receptor and cell
signaling. Gen Comp Endocrinol 2013;190:10-17. |
|
|
|
|
|
12 Bonini JA, Jones KA, Adham N, Forray C,
Artymyshyn R, Durkin MM, Smith KE, Tamm JA, Boteju LW, Lakhlani PP, Raddatz
R: Identification and characterization of two G protein-coupled receptors for
neuropeptide FF. J Biol Chem 2000;275:39324-39331. |
|
|
|
|
|
13 Osugi T, Daukss D, Gazda K, Ubuka T, Kosugi
T, Nozaki M, Sower SA, Tsutsui K: Evolutionary origin of the structure and
function of gonadotropin-inhibitory hormone: insights from lampreys.
Endocrinology 2012;153:2362-2374. |
|
|
|
|
|
14 Liu Q, Guan XM, Martin WJ, McDonald TP,
Clements MK, Jiang Q, Zeng Z, Jacobson M, Williams DL, Yu H, Bomford D:
Identification and characterization of novel mammalian neuropeptide FF-like
peptides that attenuate morphine-induced antinociception. J Biol Chem
2001;276:36961-36969. |
|
|
|
|
|
15 Ikemoto T, Park MK: Chicken RFamide-related
peptide (GnIH) and two distinct receptor subtypes: identification, molecular
characterization, and evolutionary considerations. J Reprod Dev
2005;51:359-377. |
|
|
|
|
|
16 Yin H, Ukena K, Ubuka T, Tsutsui K: A novel G
protein-coupled receptor for gonadotropin-inhibitory hormone in the Japanese
quail (Coturnix japonica): identification, expression and binding activity. J
Endocrinol 2005;184:257-266. |
|
|
|
|
|
17 Dardente H, Birnie M, Lincoln GA, Hazlerigg
DG: RFamide‐related peptide and its cognate receptor in the sheep: cDNA
cloning, mRNA distribution in the hypothalamus and the effect of photoperiod.
J Neuroendocrinol 2008;20:1252-1259. |
|
|
|
|
|
18 Ubuka T, Kim S, Huang YC, Reid J, Jiang J,
Osugi T, Chowdhury VS, Tsutsui K, Bentley GE: Gonadotropin-inhibitory hormone
neurons interact directly with gonadotropin-releasing hormone-I and-II neurons
in European starling brain. Endocrinology 2008;149:268-278. |
|
|
|
|
|
19 Ubuka T, Morgan K, Pawson AJ, Osugi T,
Chowdhury VS, Minakata H, Tsutsui K, Millar RP, Bentley GE: Identification of
human GnIH homologs, RFRP-1 and RFRP-3, and the cognate receptor, GPR147 in
the human hypothalamic pituitary axis. PloS One 2009;4:e8400. |
|
|
|
|
|
20 Zhang Y, Li S, Liu Y, Lu D, Chen H, Huang X,
Liu X, Meng Z, Lin H, Cheng CH: Structural diversity of the GnIH/GnIH
receptor system in teleost: its involvement in early development and the
negative control of LH release. Peptides 2010;31:1034-1043. |
|
|
|
|
|
21 Qi X, Zhou W, Li S, Lu D, Yi S, Xie R, Liu X,
Zhang Y, Lin H: Evidences for the regulation of GnRH and GTH expression by
GnIH in the goldfish, Carassius auratus. Mol Cell Endocrinol 2013;366:9-20. |
|
|
|
|
|
22 Bentley GE, Ubuka T, McGuire NL, Chowdhury
VS, Morita Y, Yano T, Hasunuma I, Binns M, Wingfield JC, Tsutsui K:
Gonadotropin-inhibitory hormone and its receptor in the avian reproductive
system. Gen Comp Endocrinol 2008;156:34-43. |
|
|
|
|
|
23 Maddineni S, Ocón‐Grove OM,
Krzysik‐Walker SM, Hendricks Iii GL, Proudman JA, Ramachandran R:
Gonadotrophin‐inhibitory hormone receptor expression in the chicken
pituitary gland: potential influence of sexual maturation and ovarian steroids.
J Neuroendocrinol 2008;20:1078-1088. |
|
|
|
|
|
24 Simonin F, Schmitt M, Laulin JP, Laboureyras
E, Jhamandas JH, MacTavish D, Matifas A, Mollereau C, Laurent P, Parmentier
M, Kieffer BL: RF9, a potent and selective neuropeptide FF receptor
antagonist, prevents opioid-induced tolerance associated with hyperalgesia.
Proc Natl Acad Sci 2006;103:466-471. |
|
|
|
|
|
25 Pineda R, Garcia-Galiano D, Sanchez-Garrido
MA, Romero M, Ruiz-Pino F, Aguilar E, Dijcks FA, Blomenröhr M, Pinilla L, Van
Noort PI, Tena-Sempere M: Characterization of the inhibitory roles of RFRP3,
the mammalian ortholog of GnIH, in the control of gonadotropin secretion in
the rat: in vivo and in vitro studies. Am J Physiol Endocrinol Metab
2010;299:39-46. |
|
|
|
|
|
26 Saito TH, Nakane R, Akazome Y, Abe H, Oka Y:
Electrophysiological analysis of the inhibitory effects of FMRFamide-like
peptides on the pacemaker activity of gonadotropin-releasing hormone neurons.
J Neurophysiol 2010;104:3518-3529. |
|
|
|
|
|
27 Caraty A, Blomenröhr M, Vogel GM, Lomet D,
Briant C, Beltramo M: RF9 powerfully stimulates gonadotrophin secretion in
the ewe: evidence for a seasonal threshold of sensitivity. J Neuroendocrinol
2012;24:725-736. |
|
|
|
|
|
28 García-Galiano D, van Ingen Schenau D, Leon S,
Krajnc-Franken MA, Manfredi-Lozano M, Romero-Ruiz A, Navarro VM, Gaytan F,
van Noort PI, Pinilla L, Blomenröhr M: Kisspeptin signaling is indispensable
for neurokinin B, but not glutamate, stimulation of gonadotropin secretion in
mice. Endocrinology 2012;153:316-328. |
|
|
|
|
|
29 Liu X, Herbison AE: RF9 excitation of GnRH
neurons is dependent upon Kiss1r in the adult male and female mouse.
Endocrinology 2014;155:4915-4924. |
|
|
|
|
|
30 Min L, Leon S, Li H, Pinilla L, Carroll RS,
Tena-Sempere M, Kaiser UB: RF9 acts as a KISS1R agonist in vivo and in vitro.
Endocrinology 2015;156:4639-4648. |
|
|
|
|
|
31 Korthanke CM, Thorson JF, Prezotto LD, Welsh
Jr TH, Cardoso RC, Williams GL: Secretion of Gonadotropins in Response to a Novel
Kiss-1 Receptor Agonist, RF9 in the Mare: Modulation by Estradiol-17β
and Half-Life of RF9 in the Peripheral Circulation. J Equine Vet Sci
2017;57:100-106. |
|
|
|
|
|
32 Ullah R, Batool A, Wazir M, Naz R, Rahman TU,
Wahab F, Shahab M, Fu J: Gonadotropin inhibitory hormone and RF9 stimulate
hypothalamic-pituitary-adrenal axis in adult male rhesus monkeys.
Neuropeptides 2017;66:1-7. |
|
|
|
|
|
33 Elhabazi K, Humbert JP, Bertin I, Quillet R,
Utard V, Schneider S, Schmitt M, Bourguignon JJ, Laboureyras E, Boujema MB,
Simonnet G: RF313, an orally bioavailable neuropeptide FF receptor
antagonist, opposes effects of RF-amide-related peptide-3 and opioid-induced
hyperalgesia in rodents. Neuropharmacology 2017;118: 188-198. |
|
|
|
|
|
34 Zhang J, Yang J, Jang R, Zhang Y:
GPCR-I-TASSER: a hybrid approach to G protein-coupled receptor structure
modeling and the application to the human genome. Structure
2015;23:1538-1549. |
|
|
|
|
|
35 Berendsen HJ, van der Spoel D, van Drunen R.
GROMACS: a message-passing parallel molecular dynamics implementation. Comput
Phys Commun 1995;91:43-56. |
|
|
|
|
|
36 Van Der Spoel D, Lindahl E, Hess B, Groenhof
G, Mark AE, Berendsen HJ: GROMACS: fast, flexible, and free. J Comput Chem
2005;26:1701-1718. |
|
|
|
|
|
37 Morris GM, Huey R, Lindstrom W, Sanner MF,
Belew RK, Goodsell DS, Olson AJ: AutoDock4 and AutoDockTools4: Automated
docking with selective receptor flexibility. J Comput Chem 2009;30:2785-2791. |
|
|
|
|
|
38 DeLano WL: The pymol molecular graphics system, 2002. URL: http://www. pymol.org. |
|
|
|
|
|
39 Huang J, Rauscher S, Nawrocki G, Ran T, Feig
M, de Groot BL, Grubmüller H, MacKerell AD: CHARMM36m: an improved force
field for folded and intrinsically disordered proteins. Nat Methods
2017;14:71-73. |
|
|
|
|
|
40 Jo S, Kim T, Iyer VG, Im W. CHARMM‐GUI:
a web‐based graphical user interface for CHARMM. J Comput Chem
2008;29:1859-1865. |
|
|
|
|
|
41 Bussi G, Donadio D, Parrinello M: Canonical
sampling through velocity rescaling. J Chem Phys 2007;126:014101. |
|
|
|
|
|
42 Hess B, Bekker H, Berendsen HJ, Fraaije JG.
LINCS: a linear constraint solver for molecular simulations. J Comput Chem
1997;18:1463-1472. |
|
|
|
|
|
43 Kumari R, Kumar R; Open Source Drug Discovery
Consortium, Lynn A: g_mmpbsa--a GROMACS tool for high-throughput MM-PBSA
calculations. J Chem Inf Model 2014;54:1951-1962. |
|
|
|
|
|
44 Kumbhar BV, Borogaon A, Panda D, Kunwar A:
Exploring the origin of differential binding affinities of human tubulin
isotypes αβII, αβIII and αβIV for
DAMA-colchicine using homology modelling, molecular docking and molecular
dynamics simulations. PloS One 2016;11:e0156048. |
|
|
|
|
|
45 Tripathi S, Srivastava G, Singh A, Prakasham
AP, Negi AS, Sharma A: Insight into microtubule destabilization mechanism of
3,4,5-trimethoxyphenyl indanone derivatives using molecular dynamics
simulation and conformational modes analysis. J Comput Aided Mol Des
2018;32:559-572. |
|
|
|
|
|
46 Kauffman C, Rangwala H, Karypis G: Improving
homology models for protein-ligand binding sites. Comput Syst Bioinformatics
Conf 2008;7:211-222. |
|
|
|
|
|
47 Rather MA, Basha SH, Bhat IA, Sharma N,
Nandanpawar P, Badhe M, Gireesh-Babu P, Chaudhari A, Sundaray JK, Sharma R:
Characterization, molecular docking, dynamics simulation and metadynamics of
kisspeptin receptor with kisspeptin. Int J Biol Macromol 2017;101:241-253. |
|
|
|
|
|
48 Jackson GE, Pavadai E, Gäde G, Timol Z,
Andersen NH: Interaction of the red pigment-concentrating hormone of the
crustacean Daphnia pulex, with its cognate receptor, Dappu-RPCHR: a nuclear
magnetic resonance and modeling study. Int J Biol Macromol 2018;106:969-978. |
|
|
|
|
|
49 Manogar P, Vijayakumar S, Rajalakshmi S,
Pugazhenthi M, Praseetha PK, Jayanthi S: In silico studies on CNR1 receptor
and effective cyanobacterial drugs: Homology modelling, molecular docking and
molecular dynamic simulations. Gene Rep 2019;17:100505. |
|
|
|
|
|
50 Rather MA, Dutta S, Guttula PK, Dhandare BC,
Yusufzai SI, Zafar MI: Structural analysis, molecular docking and molecular
dynamics simulations of G-protein-coupled receptor (kisspeptin) in fish. J
Biomol Struct Dyn 2020;38:2422-2439. |
|
|
|
|
|
51 Vijayakumar S, Manogar P, Prabhu S,
Pugazhenthi M, Praseetha PK: A pharmacoinformatic approach on Cannabinoid
receptor 2 (CB2) and different small molecules: Homology modelling, molecular
docking, MD simulations, drug designing and ADME analysis. Comput Biol Chem
2019;78:95-107. |
|
|
|
|
|
52 DuBay KH, Geissler PL: Calculation of
proteins' total side-chain torsional entropy and its influence on
protein-ligand interactions. J Mol Biol 2009; 391: 484-497. |
|
|
|
|
|
53 Feng Z, Hou T, Li Y. Studies on the
interactions between β2 adrenergic receptor and Gs protein by molecular
dynamics simulations. J Chem Inf Model 2012;52:1005-1014. |
|
|
|
|
|
54 Goossens K, De Winter H: Molecular Dynamics
Simulations of Membrane Proteins: An Overview. J Chem Inf Model
2018;58:2193-2202. |
|
|
|
|
|
55 Huber T, Botelho AV, Beyer K, Brown MF:
Membrane model for the G-protein-coupled receptor rhodopsin: hydrophobic
interface and dynamical structure. Biophys J 2004;86:2078-2100. |
|
|
|
|
|
56 Guixà-González R, Albasanz JL,
Rodriguez-Espigares I, Pastor M, Sanz F, Martí-Solano M, Manna M,
Martinez-Seara H, Hildebrand PW, Martín M, Selent J: Membrane cholesterol
access into a G-protein-coupled receptor. Nat Commun 2017;8:1-12. |
|
|
|
|
|
57 Saeedimasine M, Montanino A, Kleiven S, Villa
A: Role of lipid composition on the structural and mechanical features of
axonal membranes: a molecular simulation study. Sci Rep 2019;9:1-12. |
|
|
|
|
|
58 Hospital A, Goñi JR, Orozco M, Gelpí JL:
Molecular dynamics simulations: advances and applications. Adv Appl Bioinform
Chem 2015;8:37-47. |
|
|
|
|
|
59 Kufareva I, Abagyan R: Methods of protein
structure comparison; in Orry A, Abagyan R (eds): Homology Modeling, Methods
and Protocols. Humana Press, 2011, pp 231-257. |
|
|
|
|
|
60 Lobanov MY, Bogatyreva NS, Galzitskaya OV:
Radius of gyration as an indicator of protein structure compactness. Mol Biol
2008;42:623-628. |
|
|
|
|
|
61 Bhandare VV, Kumbhar BV, Kunwar A:
Differential binding affinity of tau repeat region R2 with neuronal-specific
β-tubulin isotypes. Sci Rep 2019;9:1-12. |
|
|
|
|
|
62 Kumbhar BV, Panda D, Kunwar A: Interaction of
microtubule depolymerizing agent indanocine with different human αβ
tubulin isotypes. PloS One 2018;13:e0194934. |
|
|
|
|
|
63 Santoshi S, Naik PK: Molecular insight of
isotypes specific β-tubulin interaction of tubulin heterodimer with
noscapinoids. J Comput Aided Mol Des 2014;28:751-763. |
|
|
|
|
|
64 Liao SY, Mo GQ, Chen JC, Zheng KC:
Exploration of the binding mode between (−)-zampanolide and tubulin
using docking and molecular dynamics simulation. J Mol Model 2014;20:2070. |
|
|
|
|
|
65 Ge Y, Wu J, Xia Y, Yang M, Xiao J, Yu J:
Molecular dynamics simulation of the complex PBP-2x with drug cefuroxime to
explore the drug resistance mechanism of Streptococcus suis R61. PLoS One
2012;7:e35941. |
|
|
|
|
|
66 Natarajan K, Senapati S: Understanding the
basis of drug resistance of the mutants of αβ-tubulin dimer via
molecular dynamics simulations. PLoS One 2012;7:e42351. |
|