Structural Consequences of the Polymorphism Q223R in the Human Leptin Receptor: A Molecular Dynamics Study

Leptin Receptor (LEPR) is a component of a signaling pathway related to appetite and energy expenditure. Single Nucleotide Polymorphisms (SNP) of Leptin receptor gene (lepr) have been proposed as possible modulator of adipose tissue and body weight. The main phenomenological consequence reported of these SNPs is the modulation of the LEP-LEPR interaction promoting the weight gain. Particularly, Q223R polymorphism has been associated with human obesity in some populations. In this work, we analyze the structural effects of Q223R substitution in a model of the extracellular region of LEPR comparing the stability between LEPR-Q and its Q223R variant (rs1137101) by Molecular Dynamics (MD) simulations. These results showed different behavior between both molecules after one nanosecond (ns) of simulation and significant differences in the secondary structure content were evidenced.


INTRODUCTION
Leptin Receptor (LEPR) is a membrane protein of the cytokine I class family that is a regulator of food intake and energy expenditure and is expressed in the nucleus arcuatus in the hypothalamus, where it stimulates several orexigenic neuropeptides (Woods and Stock, 1996;Loos and Bouchard, 2003;Dalamaga et al., 2013;Carrillo-Vazquez et al., 2013). Genetic mutations in LEPR have been associated to obesity in some populations. (Matsuoka et al., 1997;Silver et al., 1997;Wauters et al., 2001;Heo et al., 2002;Quinton et al., 2001;Guizar-Mendoza et al., 2005). Different theories about the origin of obesity have been proposed, one of these theories is the phenomenon called Leptin resistance and it is related to the pathologically increase of Leptin (LEP) levels due to saturation of the transport of this hormone in obese persons (Banks et al., 1996;Bjorbaek et al., 1999;Bahrenberg et al., 2002;Scarpace et al., 2005;Myers et al., 2010;Fuentes et al., 2010). The interaction of Leptin with its receptor has been investigated through amino acid replacements and deletion experiments demonstrating that some regions affect the LEP-LEPR interaction (Bahrenberg et al., 2002;Devos et al., 1997;Fong et al., 1998;Iserentant et al., 2005;Peelman et al., 2004;. Therefore, Leptin resistance could be related to mutations and SNPs of Leptin receptor gene (lepr) and recently associated to increase risk of insulin resistance and metabolic syndrome (Bjorbaek et al., 1999;Bjorbaek et al., 1997;Enriori et al., 2006;Phillips et al., 2010). In particular, Q223R (rs1137101) polymorphism is one of the most studied SNPs in LEPR that has been reported with these properties Science Publications (Gotoda et al., 1997;Thompson et al., 1997;White et al., 1997;Yiannakouris et al., 2001;Chagnon et al., 2000;Duarte et al., 2006;Fairbrother et al., 2007;Guizar-Mendoza et al., 2005).

AJABS
An in silico study published by Hiroike and coworkers explained LEP-LEPR interaction, they constructed a model of LEPR using the Granulocyte Colony-Stimulating Factor Receptor (GSF-R-PDB 1BGC) as template. By substituting amino acids and a subsequent minimization of energy, the human LEPR protein was modeled. As a result, the model was consistent with the experimental results of Fong and coworkers with a stoichiometry 2:2 of LEP-LEPR complex (Hiroike et al., 2000). Another report describes the LEP-LEPR interaction by several three-dimensional models using different computational strategies implemented in MOE, MODELLER, T-COFFE and other softwares. The template for this modeling was CHR domain (PDB ID: 1I1R) and Ig-like domain (1P9M PDB file) of the IL6 receptor, to obtain a hexameric complex model. Thus, by a 2:4 relation, a platform for mutagenesis and residue conservation studies was offered. Consequently, an effort for obtaining a better atomic description of the LEPR stability in its polymorphic variants is not pointless, particularly when such studies are made in variants associated with obesity (Bahrenberg et al., 2002;Devos et al., 1997;Peelman et al., 2004;Zabeau et al., 2003).
Computational simulations of the stability of proteins, in particular Molecular Dynamics (MD) technique and its variations permit to formulate mechanistic hypotheses on the movements of a protein in atomic terms. Although the solution of Newtonian equations of thousands of particles requires numerical methods, the growing efficiency of software and hardware makes more affordable these calculations (Huet and Derreumaux, 2006;Perryman et al., 2004;Della-Longa and Arcovito, 2013;Ilizaliturri-Flores et al., 2013).
The aim of the current study was to investigate the structural consequences in LEPR when a single amino acid glutamine (Q) is substituted by arginine (R) in its position 223. Two models were constructed representing both genetic variants of the human Leptin Receptor, namely the Wild Type (LEPRWT) and with the polymorphism Q223R (LEPRR). Once the MD simulation was reached for both models, the geometric parameters that provide details on the stability of LEPR were calculated. An increased instability of the protein due to 223R is evidenced and functional consequences for LEPR are proposed such as a variation in its binding with LEP and its homodimerization. Therefore, obesity phenotype could be associated to the alteration of these processes via Leptin-resistance in this pathological condition.

MATERIALS AND METHODS
The sequence of human LEPR was obtained from Uniprot database with ID P48357. To find the template for homology model we use Blast program in order to construct the homology model of LEPRWT and LEPRR. As for LEPR, two previously reported structures presented the necessary homology domains for the modeling, in this case IL-6 Receptor (1IL-R) served as the template. Using several techniques of structure modeling (ab initio, threading and homology modeling software and servers) we obtained the models of the LEPRWT and the mutant LEPRR. All models were analyzed stereochemically by Ramachandran Plot tool in Wincoot software (Emsley and Cowtan, 2004;Emsley et al., 2010;Joosten et al., 2011). PyMOL software (Delano, 2002;Lill and Danielson, 2011) and Magic Fit tool of Swiss-PDB Viewer (Guex and Peitsch, 1997;Herdy et al., 2004) provided the images of the LEPR models and generated the structural alignments and superpositions, respectively. Molecular Dynamics (MD) simulations were performed by GROMACS 4.0 software. The system was constituted by the LEPR protein and water; afterwards we performed an energy minimization to obtain the initial coordinates for the phase of simulations of production. MD simulations were performed using a NVT ensemble at 300 K during 1 ns (Hess et al., 2008;Rosas-Trigueros et al., 2011).

RESULTS
First, to obtain the most suitable models of LEPRWT and LEPRR, some conditions were considered. The model must contain the region where Q223R polymorphism is located and all the neighboring amino acids of the polymorphism. Also, the template used must be similar in structure or show a high homology to the Leptin receptor. In both models, the stereochemistry was analyzed in terms of the best composition of psi and phi angles in a Ramachandran plot. Initially, the ab initio modeler Rosetta server (Bonneau et al., 2001) was used to predict the structure of the protein and its variant, but the length of 615 residues of the sequence of interest was intractable for this server. When the threading modelers Phyre, Loop and I-Tasser were used, 10, 209 and 4 models were produced, respectively; however, the Science Publications AJABS resulting structures showed a not negligible difference with other reported structures (Bennett-Lovsey et al., 2008;Meller and Elber, 2001;Teodorescu et al., 2004;Tobi and Elber, 2000;Wu et al., 2007;Zhang, 2007;2008). Furthermore, these models generated by these softwares, which take into account the alignment and folding, presented a lot of gaps in the sequence and in consequence a destabilization of the region of interest. Also, various models showed inconsistencies in the secondary structure when they were compared with the templates previously selected. In the third approach, the homology modeling programs 3DJigsaw, EasyPred and Genod3D predicted two, three and three models, respectively (Bates et al., 2001;Bates and Sternberg, 1999;Combet et al., 2002;Guex et al., 1999;Lambert et al., 2002). The most suitable 3-D models of LEPRWT and LEPRR were predicted by EasyPred3D Server (http://www.fundp.ac.be/sciences/biologie/urbm/bioinfo/ esypred/). Both models used the IL6-R as template and include the region of the polymorphism Q223R. It is worthy to mention that the template presented the elements to generate a model that contains the motifs that are common in the cytokine family. Three motifs of this family are mainly conserved: fibronectin, Ig-like and cytokine domains. The results of the BLAST tool confirmed the correct election of the work of Hiroike and Peelman, where the GCSF-R and the IL-6-R played the role of templates. The sequences and atomic coordinates were obtained from Protein Data Bank, (PDB code 1I1R for IL6R and 2D9Q for GCSFR). An alignment of the sequences evidenced the conserved regions that are necessary for the models. The result of the LANLVIEW analysis showed the percentages of homology between LEPRWT and GCSF-R were 24.5%, whereas 24.1% when the IL6-R sequence was used. Regardless the higher homology percentage of the GCSF-R, we decided to employ the IL6-R model as a template for LEPRWT and LEPRR due to the previously reported work (Bazan, 1990;Hiroike et al., 2000;Peelman et al., 2004;. Stereochemical properties of both models were analyzed and optimized by WinCoot program by the Ramachandran Plot tool, obtaining acceptable results in terms of allowed angles (Data not shown). Additionally, LEPRWT and LEPRR models included the extracellular domain of the receptor and the residue 223, whose change offers the two polymorphic forms studied in this work. LEPRWT was exclusively modeled from residue 23 to 322, atomic coordinates of the remaining regions (from 1 to 22 and from 159 to 174) were infeasible to obtain. Using this sequence with the substitution of Q by R in the position 223 (LEPRR), both sequences of LEPRWT and LEPRR were submitted to the servers employed. As for LEPRR, we obtained a model of the extracellular region that spans from residue 32 to 322, i.e., eight residues more than in the LEPRWT model, although there is a gap from residue 156 to 172. Both models were compared by a structural alignment with IL6-R by the MAGIC FIT Tool of SPDBViewer software in order to visualize and understand which regions were used to generate the model. LEPRWT and LEPRR models were submitted to MD simulations during 1 ns and the last snapshot was stereochemically analyzed. To gain insight in the global behavior of the LEPR when the polymorphism Q223R is present, the Root Square Mean Deviation (RSMD) and the fluctuation of the alpha carbons were calculated (Bruschweiler, 2003;Rosas-Trigueros et al., 2011). Firstly, the RMSD showed that the equilibration was reached more rapidly in the WT (150 ps) model than in the mutant model (200 ps). Even after 800 ps, the LEPRR structure could hardly be considered as an equilibrated system, suggesting that Q confers more stability to LEPRWT in contrast to R in the LEPRR ( Fig.  1a and b). RMSD fluctuation analysis of alpha carbons showed that the two structures have notable differences in the region constituted from residue 143 to 233 (Fig.  1c, Blue Square). Therefore, the presence of the amino acid R could produce instability of the secondary structure. Specifically, four peaks (red vertical arrows) around Q223R (purple vertical dotted line) represent the corresponding residues with more fluctuation as a consequence of the arginine. Also, we observed that the entire region modified the localization of the carbons nearby. A superposition of the initial structures of LEPRWT and LEPRR, in terms of the MD simulations at the time equal zero (t = 0) was analyzed. In the LEPRWT, we observed that Q is part of a beta sheet; however, in the LEPRR with the substitution by R the EASYPred program was unable to model this secondary structure even using the same template (IL6R) (Fig. 2a).
To complete the analysis of the alpha carbons fluctuations, a superposition of the LEPRWT and LEPRR models at two different times (t = 0 and 1 ns) was performed. Comparing the movements of the regions nearby residue 223, LEPRWT showed subtle differences between the initial and final time of the MD simulation: the position of the Q223 is almost invariable and the secondary structure is still conserved (Fig. 2b). Similarly, the LEPRR structure was compared with itself at 0 and 1 ns. Structural alignment evidenced a divergence between both models at distinct times. Furthermore, a loss of secondary structure is accompanied by a large displacement of the structure. Snapshot at 1 ns showed a significant decreasing of secondary structure in the region near the residue 233, suggesting a destabilization of this domain . Besides, the beta-sheet nearest to residue 223 is no longer structured at 1 ns (Fig. 3a). Also, superposition of MD final structures of LEPRWT and LEPRR showed a higher structural divergence in the latter than wild type (Fig.   3b), whereby this genetic variant could induce a slight loss of secondary structure with potential physiological consequences (Bjelkmar et al., 2009).

DISCUSSION
Gathering the information produced by the modeling and MD simulation, it is suggested that structural changes due to a substitution of Q by R in the Leptin receptor could alter its functionality. The first structural differences were evidenced by the modeling programs although LEPRWT and LEPRR models shared the Science Publications AJABS template. The lack of secondary structure in the LEPRR model suggested that the change of amino acid renders the modeling algorithm produces alternative structures even when it utilizes two almost identical sequences. Nevertheless some evidence shows that the accuracy of the models has been enough to investigate the structural function (Gront et at., 2012). Also, some reports have evidenced that the polymorphisms could affect the protein's functionality such as its interaction with ligands, the dimerization, or signalling. As for LEPR, fluctuations of the alpha carbons support the hypothesis that the change of a single amino acid could modify the LEP-LEPR interaction or LEPR function (Basu et al., 2011;Mahmoud et al., 2011;Marie et al., 2012;Yates and Sternberg, 2013).
The analysis of the RMSD of the LEPRWT shows that equilibrium could be associated to the secondary structure due to its conservation along the simulation. Some experimental and in silico studies have proposed that beta sheets is a hallmark of the stability of the protein structure. In this work, the LEPRWT model has more beta sheets than the LEPRR model; possibly the difference observed in the RMSD is partially due the secondary structure (Varley et al., 1993;Fogolari et al., 2011;Jager et al., 2009). LEPRR model presents less secondary structure, producing an unexpected flexibility and instability in its structure, this was confirmed by the RMSD and the RMSF values. As we mentioned, the beta sheets produce a faster stabilization of the structure and in the neighborhood of the residue 223. Also, differences in the atomic fluctuations induced by the SNP can be observed in the increase of movement of alpha carbons, suggesting that the instability of this region could affect the molecular contacts with the Leptin receptor.
These findings indicate that the polymorphism Q223R could produce molecular changes in the Leptin receptor because the conserved substitution of glutamine (neutral amino acid) by arginine (charged amino acid) involves also a modification in the physicochemical environment of the extracellular region of the LEPR. Extensive evidence of mutant proteins confirms this hypothesis and new methods of prediction of damaging missense mutation will offer new evidence of the structural consequences of the polymorphisms, particularly in the membrane receptors (Guerois et al., 2002;Adzhubei et al., 2010). This work sheds light to understand the role of Q223R in modulating the activity of the Leptin receptor and in consequence of susceptibility to gain weight. A full-length modeling of the receptor and a longer time of the simulation would be required to confirm this result Klepeis et al., 2009). Additionally, we should note that other polymorphisms of the LEPR could modulate the dynamical behavior of the extracellular region of the receptor; thus, the combination of particular mutations in the protein might promote obesity. Even more, as it is known, other nongenetic factors could contribute to Leptin-resistance and finally to favor the obese phenotype.

CONCLUSION
The polymorphism Q223R in the LEPR could modify its functionality, decreasing its stability due to the loss of secondary structure and increasing the mobility of the regions near the substitution.

ACKNOWLEDGEMENT
This research was funded by SIP-IPN (Project 20130466) and COFAA-IPN and SNI-CONACYT.