Accepted Manuscript Title: Characterization, molecular docking, dynamics simulation and metadynamics of kisspeptin receptor with kisspeptin Authors: Mohd Ashraf Rather, Syed Hussain Basha, Irfan Ahmad Bhat, Niti Sharma, Priyanka Nandanpawar, Mohan Badhe, Gireesh-Babu P, Aparna Chaudhari, Jitendra Kumar Sundaray, Rupam Sharma PII: DOI: Reference:
S0141-8130(17)30544-5 http://dx.doi.org/doi:10.1016/j.ijbiomac.2017.03.102 BIOMAC 7268
To appear in:
International Journal of Biological Macromolecules
Received date: Revised date: Accepted date:
12-2-2017 16-3-2017 18-3-2017
Please cite this article as: Mohd Ashraf Rather, Syed Hussain Basha, Irfan Ahmad Bhat, Niti Sharma, Priyanka Nandanpawar, Mohan Badhe, GireeshBabu P, Aparna Chaudhari, Jitendra Kumar Sundaray, Rupam Sharma, Characterization, molecular docking, dynamics simulation and metadynamics of kisspeptin receptor with kisspeptin, International Journal of Biological Macromoleculeshttp://dx.doi.org/10.1016/j.ijbiomac.2017.03.102 This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.
Characterization, molecular docking, dynamics simulation and metadynamics of kisspeptin receptor with kisspeptin Mohd Ashraf Rather1*, Syed Hussain Basha2, Irfan Ahmad Bhat3, Niti Sharma4, Priyanka Nandanpawar5, Mohan Badhe5, Gireesh-Babu P3, Aparna Chaudhari3, Jitendra Kumar Sundaray5 and Rupam Sharma3 1. Department of Fisheries Biology, College of Fisheries Shirgaon, Rantagiri Maharashtra, India. 2. Innovative Informatica Technologies, Mayurinagar, Miyapur, Hyderabad - 500 049, India. 3. Division of Fish Genetics and Biotechnology, Central Institute of Fisheries Education, Mumbai-61, India 4. Central Inland Fisheries Research Institute, Regional Centre, Guwahati, Assam - 781 006, India. 5. Division of Fish Genetics and Biotechnology, Central Institute of Freshwater Aquaculture, Odisha – 751 002, India.
*Correspondent authors: Mohd Ashraf Rather: Department of Fisheries Biology, College of Fisheries Shirgaon, Rantagiri, Maharashtra, India. Tel: +91 9177247605; Email:
[email protected]/
[email protected]
1
Highlights
The full-length cDNA sequence of 1786 bp kiss1r was characterised
Normal tissue expression pattern of kiss1r mRNA revealed mainly expressed in brain and testis.
ASN4, SER5, GLY7, ARG9 and PHE10 of kisspeptin and TRP7, ASN8, GLU11, ILE17, ASN19 and TYR183 of kiss1r could be crucial role players.
First paper presenting molecular docking, molecular dynamics and metadynamics of kiss1r
Abstract We report molecular characterization of the kisspeptin receptor (kiss1r), an essential gatekeeper for reproduction and onset of puberty in vertebrates. The full-length cDNA sequence of kiss1r is 1786 bp which consist of 5' UTR (untranslated region) 261 bp, 3' UTR of 424 bp and open reading frame of 1101 encoding a putative protein of 366 amino acids. Basal tissue expression pattern of kiss1r mRNA revealed that it is mainly expressed in the brain and testis. We also report the structure of the kiss1r, along with plausible activation mechanism of this receptor by kisspeptin using computational modelling and dynamic simulation approach of multiple 100 ns of timescale. A present modelling and simulations studies shed light on the molecular level of interaction, suggesting that direct hydrogen bonds between ASN4, SER5, GLY7, ARG9 and PHE10 of kisspeptin and TRP7, ASN8, GLU11, ILE17, ASN19 and TYR183 of kiss1r could be crucial role players in initial binding of receptor and the kisspeptin towards allosteric modulatory effects of kisspeptin on the receptor. To the best our knowledge, this is the first report on computational modelling and molecular dynamic simulations of kiss1r in animals shedding light on its possible mode of activation.
Keywords:
Kisspeptin receptor, ORF, Homology Modelling, Docking, Molecular
Dynamics. 2
Introduction Kisspeptin, a member of the RF-amide related peptide family, has emerged recently as an indispensable gatekeeper of various reproductive processes via its ability to activate kisspeptin receptors (kiss1rs) at puberty [1, 2]. In mice, with loss-of-function mutations of kiss1r, leads to hypogonadotropic and hypogonadism has added entirely novel view regarding the regulation of puberty in vertebrates [3, 4]. In mammalian model organisms kisspeptin produced in the brain acts through the kiss1r, or GPR54 to stimulate gonadotropin-releasing hormone (GnRH) secretion by GnRH neurons in the hypothalamus, which in turn stimulates synthesis and secretion of gonadotropins in the pituitary gland [5]. Kisspeptin signalling in the brain has been involved in many reproduction processes, including sexual differentiation of the brain itself, onset of puberty, seasonal reproduction and metabolic control of fertility in humans [6, 7]. Kisspeptin has also a non-hormonal role and was originally named metastin after its ability suppressing gene in malignant melanoma (metastasis) cells [8 ]. Here we have followed the nomenclature of kisspeptin gene and kiss1r based on Human Genome Organization Gene Nomenclature Committee (HGNC) recommendations. Agonists of kiss1r with increased half-life and its efficacy to kisspeptin in vitro and in vivo have provided beneficial applications in breeding management and gonadal maturation of various species [9,10, 11]. Whitlock et al., [9] have reported agonist of kiss1r (FTM080; 4fluorobenzoyl-Phe-Gly-Leu-Arg-Trp-NH2) arouse the gonadotropic axis of ruminants in vivo and increased half-life. Scott et al., 2013 showing excellent gonadotropin releasing activity in vivo at low doses. But, in vitro and in vivo efficacy of kiss1r agonists will not be always same [12, 13]. Decourt et al., [10] designed kisspeptin analogs and injected intramuscularly to ewes showed improved potency and dramatically enhanced pharmacodynamics, which ultimately leads synchronized ovulations in both breeding and non-breeding seasons. However, this is the no report on computational modelling and molecular dynamic simulations of kiss1r in animals Kisspeptin in non-mammalian vertebrates, like teleost, may play similar roles as in the mammals [14]. After the discovery of the leading role of kiss1r in mammalian reproductive system, in 2004 existence of the fish ortholog was reported in Nile tilapia (Oreochromis niloticus) and its co-expression in GnRH neurons was the first confirmation of the existence of a kisspeptin system in fish [15]. Due to gene duplication events and subsequent gene loss, unlike higher vertebrates, fish often express multiple kisspeptin genes and receptor encoding 3
proteins that form one or more receptor/ligand pairs [16,17]. In teleost, neuroanatomical studies of brain have revealed that the ligand/receptor pairs are structurally segregated, suggesting these systems are not unnecessary and may have separate, yet not fully distinct functions [18,19]. Till date, kiss1rs have been cloned and sequenced in different teleost fish species [20]. Moreover, two distinct kiss1r and kiss2r subtypes have been found in the fish [21], medaka (Oryzias latipes) [22] sablefish (Anoplopoma fimbria) [23] among others. Molecular endocrinology mechanisms responsible for puberty in fish are not fully understood as in other vertebrates [25]. A recent study in fish suggested distinct roles between kiss1r and kiss2r with kiss2r being pertinent to reproductive functions while kiss1r has purposes is still ambiguous [23]. However, at present there is no information available on structural features of kiss1r even at least by using computational modelling approaches. Hence inspite of kiss1rs importance, so far, its structure has not been elucidated in vertebrates. In this work, we hereby report the identification and molecular characterisation of the kiss1r in Catla calta, a commercial important aquaculture fish. Furthermore, quantitative real-time PCR (qRT-PCR) was used to detect the mRNA expression level of the
kiss1r in different tissues of C. catla. On the other hand, we have tried to elucidate its structure using computational modelling approaches and studied the impact of kisspeptin on its binding with kiss1r using molecular dynamic simulations. The components of the receptor have been deduced using state-of-the-art molecular modeling, docking, molecular dynamics (MD) and metadynamic simulations which have enabled us to study the molecular motions of the receptor at atomic level. Furthermore, we present a three-dimensional (3-D) model for kiss1r which was subjected to dynamic simulations (100 ns) in its apo state along with in the presence of its substrate kisspeptin. Besides, docking methods were applied to explore the kiss1r/kisspeptin binding pocket with respect to its ligand recognition properties and to investigate the crucial interactions responsible for this complex formation leads to the receptor activation/deactivation mechanism.
Materials and methods Sample collection: Adult C. catla were collected from a local fish farm at Maharashtra, India. They were maintained for 30 days to acclimatize in the wet laboratory of ICAR (Indian Council of Agricultural Research) Central Institute of Fisheries Education, Mumbai, India, before the start of the experiment. The fishes were fed twice daily with a diet containing 35% crude
4
protein. For molecular cloning of kiss1r, whole brain was collected from fish and stored at 800C for futher work. To study the normal expression (tissue distribution) of kiss1r mRNA level in different tissues of C. catla brain, pituitary, testis, ovary, gill, intestine, liver and kidney were taken (n=4) and stored at -80 0C. RNA extraction and cDNA synthesis of kiss1r Total RNA was extracted from all the tissues using Trizol™ reagent (Invitrogen, USA). Total RNA was treated with DNase I (Thermofisher scientific, USA) to remove any impurity of genomic DNA. First strand cDNA was synthesized from total RNA using RevertAidTM cDNA Synthesis kit (Thermofisher scientific USA) primed with OligodT primer as per the manufacturer's instructions. Molecular cloning of of kiss1r To do PCR amplification of kiss1r, oligonucleotide primers (Supplementary Table 1) were designed for kiss1r mRNA based on the conserved regions of kiss1r sequences from the closely related fish species available at National Center for Biotechnology Information (NCBI) using the Gene Runner soft (v. 4.0.9.62). For partial cds, PCR amplification kiss1r performed in a total volume of 20μL containing 2.5μL of 10X Taq buffer, 0.5μL of 10mM dNTP mix 2.5μL of 25mM MgCl2, 1μL each of 10 pmol sense and antisense primers, 1μL of cDNA, 0.25μL of 5U/uLTaq DNA polymerase and 12.25μL of nuclease-free Amplification of kiss1r was performed in a 96-well Takara PCR System (Takara, Japan). The PCR conditions included an initial denaturation at 94°C for 4 min, followed by 35 cycles of 93°C for 30 sec min, 56°C for 1 min and 72°C for 45 sec and a final extension of 10 min at 72°C. The amplified cDNA fragments were purified from the gel and then cloned into pTZ57R/T vector (Thermo Scientific, USA), sequenced and confirmed BLAST algorithm at NCBI web site (http://www.ncbi.nlm.nih.gov/blast). The full-length cDNA sequence of kiss1r consists of open reading frame (ORF), 3' and 5' UTR regions was obtained using SMARTer™ RACE (Rapid Amplification of Complimentary cDNA Ends) cDNA Amplification Kit (Clonetech, USA) following the manufacturer’s protocol. Gene specific primers (Table 1 Supplementary file A) for 3' and 5' RACE-PCR were designed from the identified partial sequence. Tissue specific distribution study using real-time quantitative PCR (qRT-PCR) Normal expression of C. catla kiss1r mRNA level was analyzed in the tissue samples mentioned above by qRT-PCR. The 25 uL reaction mix was prepared with 12.5 uL of Maxima™ SYBR Green qPCR master mix (Thermo Scientific, USA), 0.5 uL of (10 pmol) each gene-specific primer (Supplementary file A) and 1uL of cDNA. The primer efficiency for qRT-PCR was determined by investigation of a serial dilution of cDNA of different 5
samples. Among three reference/housekeeping genes, β-actin gene was the most stable gene to normalize mRNA levels for tissue distribution, rest detail of validation of housekeeping genes are given in our previous report [26]. The relative expression level of target gene to housekeeping gene was determined using 2-ΔΔCt method recommended by Livak and Schmittgen, [27] Sequence analysis The cDNA sequences of C. catla kiss1r were analyzed using the BLAST algorithm at NCBI web site (http://www.ncbi.nlm.nih.gov/blast). ORF analysis of kiss1rperformed using NCBI ORF finder. Bayesian phylogenetic analysis was performed across many representative species from different taxonomic groups with kiss1r sequence available. The sequences were initially aligned on Guidance2 server [28]. The columns, residues and sequences with unreliable alignment were filtered and the alignment was rendered suitable for the phylogenetic analysis. The phylogenetic analysis was performed in BEAST v 1.8(Bayesian Evolutionary Analysis Sampling Trees) [29]. The analysis was run incorporating a lognormal relaxed molecular clock model [30], a yule speciation tree prior [31]. The best fitting amino acid substitution model was obtained through MEGA v.6 by maximum likelihood. The BEAST MCMC chains were run twice for 10 million generations each and sampling every 5000th generation. The convergence of the MCMC chains was tested using TRACER v1.6 (http://tree.bio.ed.ac.uk/software/tracer/) after removing 10% of the samples and it was noticed that the effective sample sizes were more than 200 for all the parameters for both the runs and the combined log. Maximum clade credibility tree was created from single run in Tree Annotator v1.8 [28 ] for the 3600 trees sampled. Homology Modelling: Homology modelling of kiss1r was performed with the program Modeller software with the X-ray crystal structure of P2Y1 receptor (4XNV.pdb acquired from RCSB Protein Data Bank) with a 2.2-Å resolution as a structural template as per the protocol followed details are given elsewhere in detail [32]. For each receptor, five models with different conformations of the loops were built. The crude model with least energy coming from Modeller was capped with an acetyl group at the amino-terminus and with an N-methyl group at the carboxyl-terminus to prevent electrostatic interactions and was optimized by short 10ns of molecular dynamic simulation using the default protocol as briefly explained below and elsewhere in detail [33].
6
Molecular docking: Autodockver4.0 [34] is used to predict the binding mode between the kiss1r and kisspeptin along with the aim to reveal the amount of binding energy involved in the thermodynamically stable complex formation. The Protocol followed for carrying out the docking studies using autodock to predict binding energies and mode of binding for receptorpeptin complex formation is of default parameters similar to the protocol followed elsewhere [35]. However, special care was taken during the preparation of the kisspeptin as ligand as it contains 54 rotatable bonds in total allowing which, it is extremely difficult to predict the close to accurate binding poses. To overcome this problem with high number of rotatable bonds, we have applied a simple yet effective method to decrease the number of allowed rotations without losing the flexibility of the peptide to form best possible binding poses via allowing the rotation of amide only in the peptide as shown in figure 10. The energy scoring grid box was set to 36.48; 37.10 and 35.36 Å (x, y, and z) centered at X = -10.22; Y = 16.63; and Z = −11.22 with 0.375 angstroms grid points spacing assigned with default atomic salvation parameters. The grid box was designed such that the extracellular portion of the kiss1rwas surrounded by the three-dimensional grid box centered. Lamarckian Genetic Algorithm (LGA) [36 ] was selected as a docking engine, with all the docking parameters set to default. After each LGA run, Autodock reports the best docking mode of binding along with binding energy values for the docked complex in K.cal/mol units and the results are reported based on the cluster analysis. Binding Gibbs free energy (ΔG) is calculated as a sum of six energy terms of dispersion/repulsion, electrostatic interactions, hydrogen bonding, deviation from covalent geometry, desolvation effects, and internal ligand torsional constraints. From a total of 10 docking modes represented by LGA cluster analysis, the lowest energy docking mode with respective binding energy prediction was selected from the docking simulation.
MD simulations in aqueous solution: ‘Desmond v 3.6 Package’ [37, 38] was used to study the thermodynamic stability of the kiss1r in its apo state and in complex with kisspeptin separately. Predefined TIP3P water model [38] was used to simulate water molecules using the OPLS2005 force field. Orthorhombic periodic boundary conditions were set up to specify the shape and size of the repeating unit buffered at 10 Å distances. Boundary conditions box volume was initially calculated as 450,000 Å3 (for kiss1r in apo form); 630,000 Å3 (for kiss1rin complex with 7
kisspeptin). In order to neutralize the system electrically, appropriate counter ions were added to balance the system charge and were placed randomly in the solvated system. After building the solvated system containing protein receptor and in complex with the kisspeptin, the system was minimized and relaxed using default protocol integrated within Desmond module using OPLS 2005 force field parameters. Molecular dynamic simulations were carried out with the periodic boundary conditions in the NPT ensemble [39]. The temperature and pressure were kept at 300 K and one atmospheric pressure using Nose–Hoover temperature coupling and isotropic scaling [40] the operation was followed by running the 100ns NPT production simulation and saving the configurations thus obtained at 5ps intervals. All the docking and molecular dynamics simulation snapshots were rendered using Autodocking tools and Schrodinger’s maestro interface v9.5 respectively. Metadynamics: Metadynamics is a valuable tool which enhances sampling of the underlying free energy space by biasing against previously-visited values via making use of the specified collective variables (CVs). The biasing is achieved by periodically dropping repulsive kernels of Gaussian shape at the current location of the simulation in the phase space of the collective variables. This history-dependent potential encourages the system to explore new values of the collective variables and the accumulation of potential allows the system to cross potential barriers much more quickly than would occur in standard molecular dynamics. The periodic boundary conditions and other simulation parameters were set same as molecular dynamic simulations discussed above. Metadynamic simulation was run for 100ns with selective variables between the tip of the peptide ligands’ first amino acid and bottom of the protein receptor. Ethics Statement The guidelines of the CPCSEA (Committee for the Purpose of Control and Supervision of Experiments on Animals), Ministry of Environment and Forests (Animal Welfare Division), Govt of India was followed for maintenance and treatment of animals used in the present study. Methods used in this study were approved by the Board of Studies and authorities of the Central Institute of Fisheries Education (Deemed University), Mumbai61. Statistical Analysis The differences in kiss1r mRNA transcript levels between different tissues were tested for statistical significance using the statistical package SPSS 16.0 (USA). Data was subjected to one-way analysis of variance (ANOVA), followed by Duncan’s Multiple Range 8
Test. P < 0.05 was considered as statistically significant.
Results and discussion: Cloning and characterization ofkiss1r The full-length cDNA sequence of kiss1r was 1786 bp. The sequence analysis of kiss1r showed 1101 bp of ORF region, 261 and 424 bp 5'- and 3'-terminal UTR respectively. The ORF encoded a putative protein of 366 amino acids. Further analysis of the predicted amino acid sequences of the kiss1r revealed seven putative transmembrane domains (TMD) (Figure 1). The multiple sequence alignments of deduced amino acid sequence of C. catla revealed greater than 80% identity with all other teleost. The C. catla kiss1r amino acid sequence revealed high similarity to Cyprinus carpio (common carp) kiss1r (97%), followed by D. rerio (zebrafish) (96%). The sequence of kiss1r is available on GenBank with accession number of KM114877. We isolated cDNA encoding kiss1r from the brain of C. catla. The kiss1r ORF region has been identified from 262 to1362 bp (1101bp) in the cDNA sequecne of C. catla. In D. rerio and C. auratus (gold fish) the kiss1r has an ORF of 1101bp encoding a putative protein of 366 amino acid residues [21, 41]. Analysis of the predicted amino acid sequence of C. catla kiss1r revealed seven TMDs (Trans-membrane domain), a characteristic feature of membrane-bound G-protein coupled receptors [42]. High sequence resemblance among fish species confirmed C. catla kiss1r one of the teleost kiss1r ortholog. Phylogenetic and molecular evolutionary analysis of C. catla kiss1r showed divergence from Chinese carp, C. carpio, Gobiocypris rarus and C. auratus (Figure 2) because C. catla is Indian major carp like Lebeo rohita (Rohu) and Cirrhinus cirrhosus (Mrigal) are
different species from
Chinese carps, but both Chinese and Indian carp belong to same family i.e., cyprinidae.
Overall architecture of kiss1r The overall structure of kiss1r is represented schematically in Figures 3 & 4. Kiss1r is 366 residues length canonical shape architecture with seven-transmembrane helices like other known GPCR structures. Like some other GPCRs, the kiss1r structure lacks helix VIII. Initial 30 residues at the N-terminal end were appeared to be disordered. Small beta-sheet between ASP179-TYR181 and LYS187-CYS189 in loop 4 along with a small beta sheet between VAL237-ASN239; HIS244-LEU246 and PHE281-SER283 in loop 5 was observed. Current model of kiss1r can act as a template of GPCR based on its secondary structural elements which are well within permissible range as per the Ramachandran plot analysis (figure 4b),
9
moreover its molecular size of 366 amino acids, is very close to bovine rhodopsin 348 amino acid size and thus can feature most of the essential parts of functional importance in Gprotein activation. The lengths of the three extracellular loops and seven transmembrane helices and many of the features are expected to be nearly the same for most of the GPCR family [43] with a difference in the arrangement of the seven helices. The variations in the arrangement of the helices along with their corresponding amino acid position are speculated to be the reasons behind its functionality specificity. Expression pattern of kiss1r Quantitative real-time polymerase chain reaction (qRT-PCR) analysis revealed that kiss1r mRNA transcript were expressed more in brain and testes while low expression level was detected in the pituitary, ovary, kidney, liver and gills. Interestingly, there was significant difference in mRNA transcript of kiss1r expression in testes and ovaries (figure 5). This expression pattern of kiss1r shows similarity with zebrafish and goldfish [21, 41]. Tariq et al., [44] reported in rhesus monkey Kiss1 gene and Kiss1r expression in adult testes, which was earlier validated in human testes, and immunoreactivity for Kiss1 and Kiss1r in spermatocytes, spermatids and sertoli cells respectively. Furthermore, the study suggested that the kisspeptin system could regulate spermatogenesis and play a role in sperm motility. A recent study in sablefish also showed ovaries and testes exhibited opposite Kiss1r expression profiles [23]. Molecular docking We have performed the molecular docking studies in order to study the detailed molecular basis of interactions and to estimate the binding affinity of the kisspeptin with kiss1r active site. A prepared structure of kisspeptin was docked into the active site of kiss1r located at the extracellular side of the protein. The docking results of the kisspeptin with kiss1r as per Autodock calculations have been reported as -11.3 Kcal/mol. The results of docking simulations demonstrated high binding affinity along with good interactions with plausibly critical residues present at the extracellular opening channel of the kiss1r. As per the docking simulations, kisspeptin was successfully docked into kiss1r with −11.3 kcal/mol of binding energy and shown to be forming direct hydrogen bond between ASN4, SER5, GLY7, ARG9 and PHE10 of kisspeptin and TRP7, ASN8, GLU11, ILE17, ASN19 and TYR183 of kiss receptor (Figure 6). Apart from direct hydrogen bonds many other hydrophobic and VdW forces are thought to play an important role towards stabilizing this complex. Present docking simulation results and the interaction revealed has shown significant evidence for reports of Gutierrez et al., [13] who reported that positions 6 and 10 10
are essential for KP-10 action and the C-terminal amino acids of KP-10 are vital for efficient kiss1r
binding
[45].
The
main
sequence,
kisspeptin-10
(Kiss-10)
decapeptide
(FNY/FNPFGLRF) reveals conserved substitution at positions 1 and 10 with Phenylalanine (F) amino acid among cyprinid fishes and 80-90% with other animals [46]. It is (Kiss-10) functional peptide whose administration in fishes is having stimulatory influence at the hypothalamic–pituitary axis involving kisspeptin receptors [47-48]. MD simulations in aqueous solution: Since molecular docking represents only a single snapshot of the interactions, we have performed molecular dynamic simulations to study the protein–ligand interactions in motion contributing for their stable bound conformation and to visualize the effect of ligand binding on protein conformational changes. For better understanding, the effect of kisspeptin on kiss1r, first we have subjected the protein kiss1r (without kisspeptin) to MD simulations, and then after we have performed MD simulations for kiss1r in complex with kisspeptin (figure 7). All the statistically significant results of these MD simulations trajectory analysis have been presented in Supplementary Table 2. To analyze the stability and conformational changes of kiss1r in its native state and in presence of ligand (i.e, kisspeptin) to better understand the effect of kisspeptin binding with kiss1r’s catalytic active site, we have performed the MD simulations of upto 100ns each. 100 nano seconds of simulation time used in this present study is of enough time for the side chain rearrangements in native as well as protein–ligand complexes in order to facilitate the most stable binding conformation [49]. The kiss1r–kisspeptin protein–ligand binding complex with the binding energy of −11.3 kcal/mol obtained using Autodock module was used for carrying out MD simulations; results of the MD simulations are tabulated in Supplementary Table 2. After completion of each simulation, we have studied the graphs for RMSD (Figure 8a), radius of gyration (ROG) (Figure 8b) and intra molecular hydrogen bonds (Figure 8c) contributions as the time dependent function of MD simulations. RMSD of the Kiss1r protein backbone (blue) found to be fluctuating upto 9.3 with a mean of 7.0 Å for the MD simulation. These averaged graphs suggest that the protein structure is unstable during most of the simulated time. When the RMSD for the trajectory of Kiss1rcomplexed with kisspeptin (red) using its initial model as a reference structure was analyzed it has shown much lowered and stabilized backbone RMSD with up to 4.4 with a mean of 3.72 Å, which is an indication of tighter bound protein-ligand bound complex. Compactness of the protein structure can be predicted using ROG [50]. We have analyzed the Kiss1r’s ROG in order to study the protein compactness variations along the 11
simulated time scale. Throughout the simulation time, Kiss1r in its apo state was measured to be having a ROG in the range of 18.5 to 19.3 Å and with an average of 19.0, with this ROG data, we can say that the protein is quite stable without any major note worthy expansion/contraction in the overall protein structure in its apo state. When the ROG graph for the Kiss1r-kisspeptin complex has evidence that the protein Kiss1r has slightly contracted as the simulation progresses in presence of kisspeptin by maintaining an average of 18.9 Å within a range of 18.3 to 19.4 Å statistically (Supplementary table 2). To further understand the underlying forces contributing for the stability of the kiss1r in its apo state as well as in complex with kisspeptin structure, we have analyzed the total number of intra molecular hydrogen bonds present in the kiss1r during the simulated timescale. Figure 8c have evidenced that the kiss1r protein structure contains an averaged 248 intra molecular hydrogen bonds within a range of 217 to 285, whereas, an averaged 257 in a range of 216 to 288 were noted for the kiss1r protein in presence of kisspeptin (Supplementary Table 2). This increase in the total number of intra molecular hydrogen bonds can be correlated with the contraction in the overall protein structure as evident with an above discussed radius of gyration results. To understand the underlying energy landscapes, we have performed a separate 100ns metadynamic simulations with selected variables at the tip of the peptide ligand and the bottom of the protein receptor. As per the metadynamic results, it was revealed that the free energy for the complex has observed to be as upto -10.2 Kcal/mol (Figure 9). From the figure 9, structural differences between the pre and post simulation structures numbered as 1 and 2 can be observed. Notably, opening of the protein channel upon binding of the kisspeptin can be observed, which was noted towards the end of the simulated timescale of 100ns. This metadynamics simulation results coupled with molecular dynamic simulation results, further enhanced our knowledge towards understand the plausible underlying mode of action predicted by Nocillado et al., [11] who demonstrated in in vitro that the full length Kiss2r in fish preferentially signal via PKC with kisspeptin-10 amidated peptide as a ligands and Gutierrez et al.,[13] showed kisspeptin-10 as full kiss1r agonistic activity and exhibit high receptor affinity and biopotency. Conclusion: A present study reveals kiss1r structure as a highly organized hepta helical transmembrane bundle. A set of residues that interacts with the kisspeptin resulting in its functionality has been elucidated. Complete revelation of this important receptor’s sequence is of high value and present computationally modeled structure of this receptor helped in 12
shedding light on its plausible molecular mechanism.
Considering dynamic simulation
results of RMSD, ROG, intra molecular hydrogen bonds and energy landscapes collectively, it can be said that current kiss1r protein structure is quite stable throughout the simulated time, and can be used as a reliable model for further studies. Moreover, analyzing the binding mode and the knowledge of crucial interactions could provide insight for the design of new drugs for treating diseases associated with the reproductive, nervous, endocrine, immune, and cardiovascular systems, among others Author’s contribution: MA, IR and NS have done bench work. PN MB and JKS design primers and analysed the sequence of kiss1r. SHB performed and analyzed the computational modelling and molecular dynamic simulation studies. MA, SHB and RS drafted the manuscript. All authors have read and approved the manuscript for submission.
Conflict of Interest: All authors declared that they have no conflict of interest.
Additional Information: Competing financial interests: The authors declare no competing financial interests.
13
References: 1. Kotani M, Detheux M, et.al. The metastasis suppressor gene KiSS-1 encodes kisspeptins, the natural ligands of the orphan G protein-coupled receptor GPR54. The Journal of Biological Chemistry 276:34631-34636 (2001). 2. Ohtaki T et.al., Metastasis suppressor gene KiSS-1 encodes peptide ligand of a Gprotein-coupled receptor. Journal of Nature; 411:613–617 (2001). 3. Roux de N et.al., Hypogonadotropichypogonadism due to loss of function of the KiSS1derived peptide receptor GPR54. Proceedings of the National Academy of Sciences U S A; 100: 10972–1097 (2003). 4. Topaloglu
AK
et.al.,
Inactivating
KISS1
mutation
and
hypogonadotropichypogonadism. New England Journal of Medicine 366: 629–635 (2012). 5. Roa, J., Tena-Sempere, M., KiSS-1 system and reproduction: comparative aspects and roles in the control of female gonadotropic axis in mammals. General and Comparative Endocrinology. 153, 132–140 (2007). 6. Oakley, A. E., Clifton, D. K. and Steiner, R. A., Kisspeptin signaling in the brain. Endocrinology Review, 30: 713–743 (2009). 7. Roa, J., Navarro, V.M., Tena-Sempere, M., Kisspeptins in reproductive biology: consensus knowledge and recent developments. Biology of Reproduction 85, 650– 660 (2011). 8. Lee, J. H., Miele, M. E., Hicks, D. J., Phillips, K. K., Trent, J. M., Weissman, B. E. and Welch, D. R., 1996. KiSS-1, a novel human malignant melanoma metastasissuppressor gene. Journal of the National Cancer Institute 88 (23): 1731-1737 9. Whitlock, B.K., Daniel, J.A., Amelse, L.L., Tanco, V.M., Chameroy, K.A. and Schrick, F.N., 2015. Kisspeptin receptor agonist (FTM080) increased plasma concentrations of luteinizing hormone in anestrous ewes. PeerJ, 3 (2015): p.e1382. 10. Decourt, C., V. Robert, K. Anger, M. Galibert, J-B. Madinier, X. Liu, H. Dardente et al. "A synthetic kisspeptinanalog that triggers ovulation and advances puberty." Scientific reports 6 (2016): 26908 11. Nocillado, J. N., J. Biran, Y. Y. Lee, B. Levavi-Sivan, Alejandro S. Mechaly, Y. Zohar, and A. Elizur. "The Kiss2 receptor (Kiss2r) gene in Southern Bluefin Tuna, Thunnusmaccoyii and in Yellowtail Kingfish, Seriolalalandi–functional analysis and isolation of transcript variants." Molecular and cellular endocrinology 2, 211-220. (2012): 14
12. Scott, G., Ahmad, I., Howard, K., MacLean, D., Oliva, C., Warrington, S., Wilbraham, D. and Worthington, P., Double‐ blind, randomized, placebo‐ controlled study of safety, tolerability, pharmacokinetics and pharmacodynamics of TAK‐ 683, an investigational metastin analogue in healthy men. British journal of clinical pharmacology, 75(2),.381-391(2013.). 13. Gutiérrez-Pascual, E., Leprince, J., Martínez-Fuentes, A. J., Ségalas-Milazzo, I., Pineda, R., Roa, J., &Pinilla, LIn vivo and in vitro structure-activity relationships and structural conformation of Kisspeptin-10-related peptides. Molecular pharmacology, 76(1), 58-67. (2009). 14. Pinilla, L et.al. Kisspeptin and reproduction: Physiological roles and regulatory mechanism. Physiolology. Review 92: 1235–1316 (2012). 15. Parhar, I. S., Ogawa, S. and Sakuma, Y., Laser-captured single digoxigenin-labeled neurons of gonadotropin-releasing hormone types reveal a novel G protein coupled receptor (Gpr54) during maturation in cichlid fish. Endocrinololgy 145: 3613–3618 (2004). 16. Kim, D.G et.al. Molecular coevolution of neuropeptides gonadotropin-releasing hormone and kisspeptin with their cognate G protein-coupled receptors. Frontier of Neuroscience. 6 (3), 1-6 (2012). 17. Pasquier, Jet.al., Comparative evolutionary histories of kisspeptins and kiss1rs in vertebrates reveal both parallel and divergent features. Frontier of Endocrinololgy.3, 173 (2012). 18. Escobar, S et.al., Expression of kisspeptins in the brain and pituitary of the European sea bass (Dicentrarchus labrax). Journal of Comparative Neurology.521, 933–948 (2013). 19. Zmora, N et.al., Differential and gonad stage-dependent roles of Kisspeptin1 and Kisspeptin2 in reproduction in the modern teleosts, Morone species. Biology of Reproduction 86, 177 (2012). 20. Tena-Sempere, M., Felip, A., Gómez, A., Zanuy, S., Carrillo, M., Comparative insights of the kisspeptin/kiss1r system: lessons from nonmammalian vertebrates. General and. Comparative. Endocrinololgy.175, 234–243 (2012). 21. Biran, J., Ben-Dor, S. and Levavi S.B., Molecular identification and functional characterization of the kisspeptin/kiss1r system in lower Vertebrates. Biology of Reproduction 79: 776–786 (2008).
15
22. Lee, Y.R et.al. Molecular evolution of multiple forms of kisspeptins and GPR54 receptors in vertebrates. Endocrinology 150, 2837–2846 (2009). 23. Fairgrieve Marian R et.al.
Molecular characterization of the gonadal kisspeptin
system: Cloning, tissue distribution, gene expression analysis and localization in sablefish (Anoplopoma fimbria). General and. Comparative. Endocrinololgy.225: 212–223(2016). 24. Mechaly, A.S., Vinas, J., Piferrer, F., Identification of two isoforms of the Kisspeptin1 receptor (kiss1r) generated by alternative splicing in a modern teleost, the Senegalese sole (Solea senegalensis). Biology of Reproduction. 80, 60–69 (2009). 25. Okuzawa, K., Puberty in teleosts. Fish Physiololgy and Biochemtry. 26, 31–41 (2002). 26. Kumari, K., Pathakota, G. B., Annam, P. K., Kumar, S., & Krishna, G. Characterisation and Validation of House Keeping Gene for Expression Analysis in Catlacatla (Hamilton). Proceedings of the National Academy of Sciences, India Section B: Biological Sciences, 85(4), 993-1000(2015). 27. Livak KJ and Schmittgen TD. Analysis of relative gene expression data using realtime quantitative PCR and the 2−ΔΔCT 543 method. Methods; 25: 402-408 (2001). 28. Sela, I., Ashkenazy, H., Katoh, K. and Pupko, T. GUIDANCE2: accurate detection of unreliable alignment regions accounting for the uncertainty of multiple parameters. Nucleic Acids Research, Jul 1; 43 (2015) (Web Server issue): W7-W14.; doi: 10.1093/nar/gkq443 (2015). 29. Drummond, A. J., Suchard, M. A., Xie, D., & Rambaut, A. Bayesian phylogenetics with BEAUti and the BEAST 1.7. Molecular biology and evolution, 29(8), 1969-1973 (2012). 30. Drummond, A. J., Ho, S. Y., Phillips, M. J., &Rambaut, A. Relaxed phylogenetics and dating with confidence. PLoSBiol, 4(5), e88(2006). 31. Gernhard, T. The conditioned reconstructed process. Journal of theoretical biology, 253(4), 769-778 (2008). 32. Shandilya, Ashutosh, et al. "A plausible mechanism for the antimalarial activity of artemisinin: A computational approach." Nature Scientific reports 3 (2513): 1-6 (2013). 33. Reddy, S. V. G., et al. "Molecular docking and dynamic simulation studies evidenced plausible immunotherapeutic anticancer property by Withaferin A targeting 16
indoleamine 2, 3-dioxygenase." Journal of Biomolecular Structure and Dynamics 33(12):2695-709: 1-15 (2015). 34. Goodsell, D. S., Morris, G. M., & Olson, A. J. Automated docking of flexible ligands. Journal of Molecular Recognition, 9, 1–5(1996). 35. Hussain BS &Naresh Kumar, K. Ligand and structure based virtual screening studies to identify potent inhibitors against herpes virus targeting gB–gH–gL complex interface as a novel drug target, Open access Scientific reports1(12), 1–13. (2012).doi:10.4172/ scientific reports.566 36. Morris, G. M et.al. Automated docking using a Lamarckian genetic algorithm and an empirical binding free energy function. Journal of Computational Chemistry, 19, 1639–1662 (1998). 37. Shivakumar, D et.al. Prediction of absolute salvation free energies using molecular dynamics free energy perturbation and the OPLS force field. Journal of Chemical Theory and Computation, 6, 1509–1519 (2010). 38. Jorgensen, W. L et.al. Comparison of simple potential functions for simulating liquid water. The Journal of Chemical Physics, 79, 926–935(1983). 39. Shinoda, W., & Mikami, M. Rigid-body dynamics in the isothermal-isobaric ensemble: A test on the accuracy and computational efficiency. Journal of Computational Chemistry, 24, 920–930 (2003). 40. Nosé, S. A unified formulation of the constant temperature molecular dynamics methods. The Journal of Chemical Physics, 81, 511 (1984). 41. Servili, A et.al. Organization of two independent kisspeptin systems derived from evolutionary-ancient kiss genes in the brain of zebrafish. Endocrinology 152, 1527– 1540 (2011). 42. Li, S., Zhang et.al. Structural and functional multiplicity of the kisspeptin/GPR54 system in goldfish (Carassius auratus). J. Endocrinol. 201, 407–418 (2009). 43. Muir, A. I et.al. AXOR12, a novel human G protein-coupled receptor, activated by the peptide KiSS-1. J. Biol. Chem., 276: 28969–28975 (2001). 44. Tariq, A.R., et.al. Kiss1 and Kiss1 receptor expression in the rhesus monkey testis: a possible local regulator of testicular function. Cent. Eur. J. Biol. 8, 968–974 (2013). 45. Roseweir, A. K., & Millar, R. P. The role of kisspeptin in the control of gonadotrophin secretion. Human reproduction update, 15(2), 203-212(2009). 46. Rather, M. A., I. A. Bhat, P. Gireesh-Babu, A. Chaudhari, J. K. Sundaray, and R. Sharma. "Molecular characterization of kisspeptin gene and effect of nano– 17
encapsulted kisspeptin-10 on reproductive maturation in Catla catla." Domestic animal endocrinology 56, 36-47. (2016) 47. Rather, Mohd Ashraf, Irfan Ahmad Bhat, Pravesh Kumar Rathor, P. Gireesh-Babu, Aparna Chaudhari, Sundaray Jeetendra Kumar, and Rupam Sharma. "In silico analysis and expression studies of kisspeptin gene in C. catla." Journal of Biomolecular Structure and Dynamics 1,1-12 (2016) 48. Rather, Mohd Ashraf, Irfan Ahmad Bhat, Niti Sharma, Rupam Sharma, Aparna Chaudhari, and Jeetendra Kumar Sundaray. "Molecular characterization, tissue distribution of Follicle‐ Stimulating Hormone (FSH) beta subunit and effect of kisspeptin‐ 10 on reproductive hormonal profile of Catla catla (Hamilton, 1822)." Aquaculture Research 47, 89-100 (2014). 49. DuBay, K. H., & Geissler, P. L. Calculation of proteins’ total side-chain torsional entropy and its influence on protein–ligand interactions. Journal of Molecular Biology, 391, 484–497 (2009). 50. Lobanov, M. Y., Bogatyreva, N. S., & Galzitskaya, O. V. Radius of gyration as an indicator of protein structure compactness. Molecular Biology, 42, 623–628 (2008).
18
Figure 1: Nucleotide and deduced amino acid sequences of cDNA encoding kiss1r in C. catla. Kiss1r protein consisted of seven putative transmembrane domains (TMD) showed in different coloured rectangles. The start and stop codons are bolded and highlighted in red colour. Dotted black rectangular square nucleotides are 5' untranslated regions and lower nucleotides underlined indicate 3' untranslated region.
19
Figure 2: Phylogenetic tree of the kisspeptin receptor of different vertebrates. Maximum clade credibility tree obtained through Bayesian phylogenetic analysis using BEAST tool. The numbers on the nodes are the posterior probability values ascertaining the strength of the analysis. The brown colour highlighted indicates kiss1r of C. catla.
20
Figure 3: Homology modelled structure of kiss1rshowing A) Front view B) Top view and C) Bottom view of the receptor displaying the arrangement of the seven transmembrane helices
21
Figure 4: a) Two-dimensional model displaying the overall architecture of the kiss1r. Residues are shown single letter codons in circles. b) Ramachandran plot of modelled kiss1r.
22
Figure 5: Tissue specific distribution of kiss1r mRNAs in the male and female C. catla as determined by Real-time PCR. Results are expressed as means ± SEM (n = 4) for each sex. Asterisks above the bars represent significant differences between sexes (P < 0.05).
23
Figure 6: a, b) Docking snapshot of kisspeptin binding at the extracellular surface opening of the kiss1r top and front view respectively c) Molecular interactions responsible for the binding of kisspeptin (blue) inside the kiss1r (interacting residues were shown as red colored sticks). Hydrogen bonds have been shown in black dotted lines along with their bond lengths, d, e) represents the binding surface of the kiss1r showing the hydrogen donor/acceptor region and hydrophobicity map respectively.
24
Figure 7: Molecular dynamic simulations snapshots of Kiss1r in its apo state (top row) and kiss1r in complex with kisspeptin (bottom row). Kiss1r has been represented in ribbon structure, whereas kisspeptin is represented in space filling model in blue color.
25
Figure 8: a) RMSD plot of kiss1r in its apo state (blue) and kiss1r in complex with kisspeptin (red). b) ROG plot of kiss1r in its apo state (blue) and kiss1r in complex with kisspeptin (red). c) Total number of intra molecular hydrogen bonds plot of kiss1r in its apo state (blue) and kiss1r in complex with kisspeptin (red).
26
Figure 9: Free energy landscape analysis from kiss2receptor and kisspeptin complex metadynamics simulation results. Numerics 1 and 2 represents the pre and post simulation representative snapshots of the complex.
27
Figure 10: Preparation steps of kisspeptin as ligand for docking calculations: a) represents the initial imported kisspeptin peptide with 54 rotatable bonds enabled by default. b) represents the kisspeptin peptide with 43 rotatable bonds enabled via making amide bonds non-rotatable and c) represents the final kisspeptin structure with 11 rotatable bonds used for further studies made by inversing the above step selection.
28