Identification of Key Regions Mediating Human Melatonin Type 1 Receptor Functional Selectivity Revealed by Natural Variants.

Melatonin is a hormone mainly produced by the pineal gland and MT1 is one of the two G protein-coupled receptors (GPCRs) mediating its action. Despite an increasing number of available GPCR crystal structures, the molecular mechanism of activation of a large number of receptors, including MT1, remains poorly understood. The purpose of this study is to elucidate the structural elements involved in the process of MT1's activation using naturally occurring variants affecting its function. Thirty-six nonsynonymous variants, including 34 rare ones, were identified in MTNR1A (encoding MT1) from a cohort of 8687 individuals and their signaling profiles were characterized using Bioluminescence Resonance Energy Transfer-based sensors probing 11 different signaling pathways. Computational analysis of the experimental data allowed us to group the variants in clusters according to their signaling profiles and to analyze the position of each variant in the context of the three-dimensional structure of MT1 to link functional selectivity to structure. MT1 variant signaling profiles revealed three clusters characterized by (1) wild-type-like variants, (2) variants with selective defect of βarrestin-2 recruitment, and (3) severely defective variants on all pathways. Our structural analysis allows us to identify important regions for βarrestin-2 recruitment as well as for Gα12 and Gα15 activation. In addition to identifying MT1 domains differentially controlling the activation of the various signaling effectors, this study illustrates how natural variants can be used as tools to study the molecular mechanisms of receptor activation.

Next-generation sequencing (NGS) is rapidly increasing our understanding of disease causes and susceptibilities in humans as well as of the adaptation and evolution of organisms. Human DNA contains 10,000 to 11,000 nonsynonymous genetic variants on average,[1] although only few of them are expected to lead to disease, while the rest are considered neutral. Information on such variants in databases has been often limited to basic annotations like “neutral” or “deleterious”. Deleteriousness can be approximated with available information on sequence conservation and molecular functionality, but such data are rarely combined. Furthermore, even a strong association between a mutation and a clinical trait does not provide specific information on the corresponding changes in protein function and/or structure.[2,3] Since just a single genetic variant can change the local structure of a protein or its interactions with other proteins, it can actually have various functional effects.[4] Thus, genetic variants could be instrumental in deciphering molecular determinants and pathways.

Many proteins display several functions by interacting with different protein partners; however, most protein-interaction maps ignore structural details and the mechanisms of protein–protein interaction. This is the case for many G protein-coupled receptors (GPCRs), which can interact with different G proteins, GPCR-regulated kinases (GRKs), and βarrestins, triggering different signaling pathways. The molecular determinants controlling these different functions are still generally not well understood. It is now known that different ligands can preferentially activate specific GPCR signaling pathways leading to distinct physiological effects. Variants at specific positions in the receptor structure can also modify the activation mechanism and thus the signaling profile of a receptor[5] by inhibiting one or several pathways or, conversely, by promoting others. This concept can be referred to as genetically defined functional selectivity or biased signaling.[6]

We took advantage of the abundance of genetic variants in MTNR1A (encoding the melatonin receptor 1 [MT1]) to explore the link between structure and functional selectivity. Melatonin (MLT) is mainly released by the pineal gland in circadian and seasonal manners.[7] The MT1 is one out of two melatonin receptors involved in the physiological action of this neuro-hormone.[8] MT1 expression has been reported in several tissues and organs[9] including the retina,[10,11] brain,[12] pancreas, and immune system.[13,14] MT1 is involved in various physiological functions including circadian entrainment,[15] sleep,[16] dopamine level regulation,[17] and glycemic homeostasis.[18]

In the present study, we identified 36 MTNR1A variants and investigated their capacity to activate seven G protein subtypes and to recruit βarrestin-2. We found that these variants can be clustered in three distinct functional groups providing insights into specific structural features of the MT1 underlying functional selectivity.

Results

Signaling Profile of Wild-Type MT1

To obtain a comprehensive profiling of the signaling repertoire of MT1, we used bioluminescence resonance energy transfer (BRET)-based biosensors to assess G protein activation measuring Gα protein dissociation from the Gβγ complex[19] (Figure S1A). We measured βarrestin-2 recruitment to the activated receptor measuring its proximity with the plasma membrane using enhanced bystander BRET (ebBRET)[20] (Figure S1B). Melatonin stimulation in cells expressing the wild-type (WT) human MT1 leads to the activation of Gαi1, Gαi2, Gαi3, GαoA, GαoB, Gα12, and Gα15 (Figure ). This is consistent with previous findings suggesting a predominant coupling to Gαi for this receptor. Gα15, an atypical member of the Gαq/11 family, has been shown to couple in a promiscuous manner to many distinct GPCRs.[21] Because the human embryonic kidney 293 (HEK293) cells used do not express Gα15, its potential activation by MT1 was assessed by probing the downstream activation of protein kinase C (PKC) using a BRET-based PKC activity sensor[22] (Figure S1C) upon Gα15 expression. The data show that no MLT-promoted PKC activation could be detected in the absence of Gα15, but the addition of Gα15 resulted in a robust PKC activation illustrating the ability of MT1 to functionally engage this protein (Figure ).

Figure 1

MT1-WT activates Gαi/o, Gα12, Gα15, and recruits βarrestin-2. (A) MLT concentration curves for MT1-WT. Data were plotted using nonlinear regression curves with a fixed Hill slope equal to 1. (B) MT1 response for Gαi1, Gαi2, Gαi3, GαoA, GαoB, Gα12, Gαq, Gα11, Gα15, and Gαs activation and βarrestin-2 recruitment at saturated MLT concentration (100 nM) of the corresponding curves in A. Statistical analysis was performed using one-way ANOVA followed by a Dunnett’s posthoc test. Data represents the mean ± SEM of 5 to 9 experiments. Data were fitted in GraphPad Prism 9. * p ≤ 0.05, ** p ≤ 0.01, **** p ≤ 0.0001. βarr2: βarrestin-2.

Identification of MT1 Variants and Cell Surface Expression

To identify natural variants, we sequenced the two exons of the MTNR1A gene in 8687 unrelated European individuals. We accurately identified 36 genetic variants, including 34 rare variants with a minor allele frequency (MAF) of <1% (Table S1). The two frequent variants are nonsynonymous variants. Among the 34 rare variants, we found 2 nonsense variants introducing STOP codons (p.W93*, p.Y170*), one variant on a donor splice site (c.184 + 1G > T), one frameshift variant (c.999_1014del/p.W333*), one loss of STOP codon variant (c.1053A > T/p.*351Yext* 18 (YKSTTFRVRWTRCARPRS)) and 29 nonsynonymous variants (Table S1). As we previously found that rare deleterious variants in MTNR1B (encoding the MT2 receptor) were significantly associated with increased risk of type 2 diabetes,[23] we assessed the association between the coding variants of MTNR1A and the risk of type 2 diabetes or obesity. Neither the frequent MTNR1A variants (i.e., p.G166E-rs28383653 and p.A266 V-rs28383652) nor the 34 rare MTNR1A variants were associated with the risk of type 2 diabetes or obesity in both children and adults (Table S2).

Then, we determined the functional properties of all 32 nonsynonymous MT1 variants in comparison to the MT1-WT in terms of cell surface expression, G protein activation, βarrestin-2 recruitment, and binding of the radioligand 2(125I)-iodomelatonin (125I-MLT).

Cell surface expression of MT1 variants was assessed by enzyme-linked immunosorbent assay (ELISA), detecting the Flag tag located at the extracellular N-terminus of the receptor. Figure S2 shows that similar expression levels were observed for all variants compared to MT1-WT, indicating that the results of the functional assays for the different MT1 variants would not be severely affected by differences in the number of receptors at the cell surface.

Signaling Profiles of MT1 Variants

Melatonin concentration–response curves were generated for each MT1 variant and MT1-WT, and both the EC50 (potency) and the maximal agonist-mediated response (maximal efficacy; Emax) were determined for βarrestin-2 recruitment and for those G proteins for which activation by MT1-WT was observed: Gαi1, Gαi2, Gαi3, GαoA, GαoB, Gα12, and Gα15 (Table , Figure S3, and Table S3). Each signaling pathway was studied using a different biosensor, and consequently, comparison of BRET intensities can only be done between receptors for the same biosensor but not between different pathways for a same variant. Representative examples of the signaling profiles for the MT1-WT, and three melatonin-responsive variants are shown in Figure and clearly indicate the differential impact of variants including selective effects on the potency and/or efficacy toward the different pathways.

Table 1

Comparison of Efficacies and Potencies of the Different MT1 Receptors to Activate Gα Proteins and Recruit βarrestin-2a

Functional profiling of Gα protein activations and βarrestin-2 recruitment by MT1-WT and MT1 variants. Data represents the mean ± SEM of 3-16 independent experiments with repeats in quadruplicate. NR (No response) denotes that the experimental parameter could not be determined due to lack of a concentration-response curve. The data were analyzed by comparing independent curves.

Figure 2

Melatonin concentration–response curves for G protein activations and βarrestin-2 recruitment of the MT1-WT and representative receptor variants with distinct signaling profiles. MLT concentration curves for (A) MT1-WT, (B) V52A, (C) R125C, and (D) L138P receptors. Data were plotted using nonlinear regression curves with a fixed Hill slope equal to 1. Data points represent means ± SEM of 14 to 16 experiments (A) or 3 to 5 experiments (B, C, D). βarr2: βarrestin-2, WT: wild-type. See also Figure S1. Data were fitted in GraphPad Prism 9.

To provide a visual representation of the relative signaling profiles of each variant for the different signaling pathways, we used a radial graph representation including the Emax and the EC50 for each pathway (Figure S4). We used a non-negative matrix factorization (NMF) and K-means clustering method to divide the variants into groups according to their functional impact similarities and differences. The different signaling profiles could be grouped in three distinct clusters. Figure presents the superposition of the individual signaling profiles of each member of a given cluster providing a global view of the cluster signaling characteristics.

Figure 3

Radial graph representing the different variant profiles compiled in 3 clusters. Compilation of the different MT1 variant radial graph profiles in three clusters. On each radius, maximal agonist-induced efficacy (Max) and potency (EC50) are indicated. MT1-WT profile is represented in red and mutant profiles are in blue. The values were normalized to fit a −1 to +1 scale using the following formula: (variant – WT)/(variant + WT) as indicated in the methods. A loss of potency or efficacy for a specific protein results in the decrease of the blue area. Data were fitted in Microsoft Excel 2016.

Cluster 1 comprises 21 variants distributed over transmembrane (TM) domains 2, 3, 4, 5, and 6, helix 8, the N-terminus (Nterm), extracellular loops (ECL) 2 and 3, and the C-terminus (Cterm) (p.G18RNterm, p.I88V2.65, p.N91Y2.68, p.G96D3.21, p.A157V4.55, p.G166E4.64, p.A180TECL2, p.I212T5.63, p.V221M5.72, p.K228R6.25, p.A266TECL3, p.A266VECL3, p.S267GECL3, p.R307S8.55, p.V321MCterm, p.D326NCterm, p.V331FCterm, p.K334NCterm, p.T340ICterm, p.N342SCterm, and p.V345ICterm) (superscripts represent Ballesteros–Weinstein nomenclature[24]) that did not significantly impact any of the signaling pathways considered in response to MLT. Cluster 2 comprises seven variants (p.V52A1.57, p.R54W1.59, p.S87L2.64, p.H131R3.56, p.I257F6.54, p.I309T8.57, and p.C314RCterm) that cause a total loss of βarrestin-2 recruitment and a decrease of Emax and/or EC50 for some Gα proteins, the most impacted being Gα12; however, for each of these variants, the activation of at least one G protein is not negatively impacted by the mutation. Cluster 3 comprises three variants (p.I112N3.37, p.R125C3.50, p.L138PICL2) (intracellular loop (ICL)) with a total loss of βarrestin-2 recruitment and a global reduction in the activation of all G protein pathways.

To exclude the possibility that differences in the signaling properties, in particular EC50 values, could result from reduced binding affinity to MLT, we determined the dissociation constant (Kd) of the melatonin receptor-specific radioligand 125I-MLT for the 20 variants with decreased EC50 values and for the melatonin-unresponsive p.P80L2.57. Kd values were not statistically different between the variants and the MT1-WT, except for p.P80L2.57, which showed a complete loss of binding capacity (Table ). Our cell surface ELISA data show that p.P80L2.57 expresses properly at the cell surface excluding a defect in trafficking (Figure S2). The p.P802.57 amino acid is part of a YPYP motif conserved in melatonin receptors. Our observation is consistent with previous studies suggesting that this motif located in transmembrane helix 2 (TM2) is crucial for both stability and activation of MT1.[25]

Table 2

MT1 Variant Receptors Bind to Melatonin as well as the MT1-WT Except P80La

receptors	Kd (pM) ± SEM
Wild Type	101 ± 9.3
V52A	99 ± 6.9
R54W	178 ± 25
P80L	No binding
S87L	136 ± 19
I88V	79 ± 4.5
G96D	81 ± 8.7
I112N	244 ± 75
R125C	227 ± 61
H131R	62 ± 11
L138P	171 ± 61
G166E	63 ± 1.2
I212T	101 ± 5.0
I257F	168 ± 7.2
A266V	224 ± 43
S267G	186 ± 54
R307S	40 ± 2.5
I309T	75 ± 5.6
C314R	48 ± 7.8
K334N	395 ± 68
N342S	140 ± 43
V345I	69 ± 1.8

Melatonin dissociation constant measured by 2-(125I)Iodomelatonin binding experiment on MT1 variants with at least one signaling defect. Data represents means ± SEM of at least three experiments.

Evolutionary Action Analysis of MT1 Variants

To assess whether the extent of the functional impact of the naturally occurring variants could be estimated by the evolutionary impact of the individual substitution, we scored the Evolutionary Action (EA) of the nonsynonymous variants. EA is a bioinformatics algorithm designed to predict the impact of mutational variants[26] in which the EA score varies from 0, for variations with the most benign impact, to 100, for those with the most adverse effects. As expected, cluster 1 is comprised of variants of low EA scores, cluster 2 of medium scores, and cluster 3 of high scores. Figure shows the statistically significant correlation (R2 = 0.580, P < 0.0001) between EA scores and the experimentally determined signaling phenotypic impact for each of the variants, indicating that the severity of the functional outcome of SNPs can be reasonably well predicted using trace evolution-based methods such as EA.

Figure 4

Evolutionary action of MT1 variants and correlation with signaling impairment. (A) An evolutionary action (EA) score was calculated for each variant. EA ranges from 0 to 100 with a score of 0 predicted to be benign and a score of 100 predicted to be highly impactful or detrimental to protein function. Scores are colored based on EA (EA = 0 (blue) to EA = 100 (red)). (B) Functional defects of every MT1 variant were correlated with the EA score. The 3 clusters were represented with the cluster 1 in blue, the cluster 2 in yellow, and the cluster 3 in red. Data were fitted in GraphPad Prism 9 by linear regression analysis, P and R2 values obtained for the overall correlation.

Structural Determinants of the Functional Impact of MT1 Variations

To investigate the structural basis of the functional impact of the MT1 variants, we analyzed the position and topology of the individual mutated residues of Clusters 1–3 on the three-dimensional models of the active form of MT1 coupled to Gαi or βarrestin-2 (see Methods).

Variants in Cluster 1 (with no impact on the receptor signaling profile) are spread throughout the structure (Figure A,B), appearing mostly at the surface of the extra- and intracellular domains and generally far from the ligand-binding pocket[25,27] or from any of the known GPCR activation microswitches.

Figure 5

Location of the mutated residues and topology of Clusters 1–3. (A) Distribution of the mutated residues depicted in a 2D snake plot of the MT1 receptor. Residues in Cluster 1 are colored gray, in Cluster 2 blue, and in Cluster 3 red; these colors are preserved in all panels. (B) Location of Cluster 1 mutations in the three-dimensional structure of inactive MT1. The receptor is colored in a rainbow spectrum (TM1 blue - TM7/helix8 red). (C) Location of Cluster 2 mutations in active MT1 bound to βarrestin-2. (D) Location of Cluster 3 mutations in active MT1 bound to βarrestin-2 and Gαi. (E) Mutation I112N3.37 is located near the key functional PIF (or triad) motif. (F) Mutation P80L2.57 is located near the binding site (shown as gray spheres) of melatonin (shown as sticks).

The seven variants in Cluster 2 (p.V52A1.57, p.R54W1.59, p.S87L2.64, p.H131R3.56, p.I257F6.54, p.I309T8.57, and p.C314RCterm) are associated with a total loss of βarrestin-2 recruitment and a decrease of different Gα activations, especially for Gα12 and/or Gα15 Emax (Table and Figure S3), but at least one Gα protein is not impacted. Interestingly, four of these variants are located at the interface between the cytoplasmic side of TM1 (p.V52A1.57 and p.R54W1.59) and helix 8 (H8) (p.I309T8.57 and p.C314RCterm) (Figure C). The short amphipathic H8 lies adjacent to the inner leaflet of the plasma membrane, anchored through a palmitoylated cysteine (p.C314 in MT1). In other GPCRs, H8 has been shown to play a crucial role in the interaction of the receptor with βarrestin-2—leading to receptor internalization[28] and its palmitoylation has been shown to enhance the recruitment of βarrestin-2.[29] The structural models show that both H8 and ICL1 of MT1 can interact with the finger loop of βarrestin-2 (Figure C). Thus, we hypothesize that p.V52A1.57, p.R54W1.59, p.I309T8.57, and p.C314RCterm variants affect βarrestin-2 recruitment by altering the local environment of H8 in MT1. Both p.I309T8.57 and p.C314RCterm affect directly H8 by weakening hydrophobic interactions with nearby phospholipids or by removing the palmitoylation site, respectively. On the other hand, p.V52A1.57 and p.R54W1.59—located at the interface between TM1 and ICL1—would affect H8 indirectly. These results are consistent with findings on other GPCRs; for instance, hydrophobic amino acids in H8 of the mouse odorant receptor mOR-S6 have been suggested to form a hydrophobic core with TM1 responsible for the proper positioning of H8, and mutations of these residues increase the flexibility of H8 and destabilize its structure.[30] Not surprisingly, the loss of βarrestin-2 recruitment by these variants leads to a loss of internalization of the receptor (Figure S5) as measured by ebBRET between the RlucII-fused receptors and rGFP located in early endosomes[20] In addition, in Cluster 2, the p.H131R3.56 variant is located at the interface between the cytoplasmic end of TM3 and ICL2, at a location where it could interact directly with βarrestin-2 (Figure C). This variation also places a positively charged residue next to the highly conserved p.Y1263.51, which might alter its function by establishing a spurious cation−π interaction. Interestingly, this region has been proposed to be a secondary site for βarrestin interaction in other GPCRs.[31,32] The p.I257F6.54 variant is located at the extracellular side of TM6 (Figure C), a site already shown to be important for the active melatonin-bound form of MT1, favoring the conformational changes necessary for signal transmission.[33] Clement et al. proposed that, similarly to rhodopsin,[34] ECL2 of MT1 stabilizes TM5, TM6, and TM7 upon activation, particularly through the interaction between ECL2 and p.I2576.54. Consequently, we suggest that the variation p.I257F6.54 alters the interactions with ECL2 decreasing the stability of the MT1 active form. Furthermore, p.I2576.54 is located one turn above the highly conserved proline-kink at p.P2536.50, a key residue in shaping TM6. The extracellular sides of TM6 and TM7—near position p.I2576.54—form important interactions to maintain the fold of the receptor,[35] and TM7 is known to be important for the interaction between GPCRs and βarrestin.[36] Thus, the p.I257F6.54 variant might also destabilize the TM6/7 interface and prevent the appropriate location of TM7 impairing interactions between MT1 and βarrestin-2. The last variant (p.S87L2.64) of Cluster 2 is located away from known potential functional sites (Figure C), which makes it difficult to rationalize its effect from a structural point of view.

We propose that residues from Cluster 2 contribute to an allosteric network between the binding pocket and the domains involved in the engagement of βarrestin-2. The fact that these variations do not affect the efficacy or potency of MT1 to activate Gαi/o family members rules out the possibility that the functional defects result from a more general misfolding of the receptor. Regarding this specific loss of βarrestin-2 recruitment and Gα12 and/or Gα15 activation, we cannot exclude the possibility that a modest nonspecific structural impairment has a more important impact on the pathways that are less efficiently coupled to the receptor. Indeed, the lower potency of MT1 to activate Gα12 and Gα15 vs the Gαi/o family members suggests that these pathways are indeed less efficiently coupled to the receptor.

The burden of the seven variants of Cluster 2 were not found to be associated with an increased risk of type 2 diabetes or obesity (Table S1 and Table S2).

The three variants of Cluster 3 (p.I112N3.37, p.R125C3.50, and p.L138PICL2) (Figure D) have a more severe effect on the signaling profile of MT1. On one hand, as in Cluster 2, they completely abolish βarrestin-2 recruitment (Figure S1), but moreover, they result in loss of potency (EC50)—but not in efficacy (Emax) for p.I112N3.37 and p.L138PICL2—in all signaling pathways compared to MT1-WT. The possibility that this loss of potency could be due to a loss of binding affinity for MLT could be ruled out, since we observe no difference in binding affinity for these three variants compared to MT1-WT (Table ). We suggest instead, that these variants lead to a destabilized active state and, thus, to a less efficient activation.

The variant p.L138PICL2 of cluster 3 is located in ICL2, which forms a short α-helix in several GPCRs when activated that can then interact with G proteins[37] and βarrestins (Figure D). The exchange of a leucine for a proline is expected to either disrupt the α-helix or at least prevent an optimal position, which would lead to a less efficient interaction with G proteins and βarrestin-2 and the observed decrease in potency.

The variant p.I112N3.37 is located in the middle of TM3 (Figure E) near residues p.I1153.40 and p.P1995.50, which are part of a key switch (the “triad” or PIF motif) for GPCR activation.[38] The p.I112N3.37 variant would alter the packing and physicochemical properties of this region, leading to a less efficient activation of the PIF motif and an impaired transduction of the conformational changes through TM6 and TM5 that lead to the formation of the active state. As in p.L138PICL2, the p.I112N3.37 variant can still activate G proteins (its Emax is unaffected), but these events are less efficient most likely due to an impaired conformational change in TM5 and TM6 leading to the observed decrease in potency.

The p.R125C3.50 variant (Figure D) changes a key residue of the D/NRY motif, known to be essential for GPCR function and which is the only natural variants directly involved in the binding of βarrestins and G proteins. The p.R125C3.50 variant shows a total loss of signaling capacity, except for GαoB which remains well activated by this variant, an observation that remains without explanation at the current state.

The burden of the three variants of Cluster 3 were not associated with an increased risk of type 2 diabetes or obesity (Table S2). The p.P80L2.57 variant, with a completely impaired MLT binding (Figure F), is located in the YPYP motif (residues p.Y792.56 to p.P822.59) that shapes the ligand binding pocket by modulating a kink in TM2.[25,39] Thus, the proline to leucine substitution most likely modifies the local structure of TM2 and disrupts the binding pocket, leading to its inability to bind MLT. This variant was very rare in the Genome Aggregation Database (gnomAD; v2.1.1) with a MAF below 0.00005 and was found in a normal-glucose, obese individual from the present study.

Discussion

Here, we took advantage of the large diversity of naturally selected variants in the coding region of the human MTNR1A gene encoding the MT1 receptor to obtain mechanistic insights into the function of this receptor. Notably, all impactful mutants show a loss-of-function phenotype, and no gain of function for the pathways engaged by the WT receptor was observed. The possible gain of function for the pathways not engaged by the WT receptor was not assessed. One mutant, p.P80L2.57, was a complete loss-of-function mutant as it could not bind MLT. Based on the functional profiles of the MLT-responsive mutants, their location in the MT1 structure, and the general knowledge about the role of specific regions in receptor function, we were able to classify these impactful natural variants into two categories: those with severely impaired overall signaling locating near regions identified as key regions of general GPCR activation and those with major defects on more specific signaling functions, i.e., βarrestin-2 recruitment, that are located in distinct regions of ICLs and H8.

Our study reports the first exhaustive characterization of the G protein activation profile of the MT1 receptor. Upon MLT stimulation, this receptor can activate a broad set of G proteins (Gαi1, Gαi2, Gαi3, GαoA, GαoB, Gα12, and Gα15) in addition to the recruitment of βarrestin-2 in transfected HEK293 cells. These results are consistent with the notion that MT1 inhibits cyclic adenosine monophosphate (cAMP) production through pertussis-toxin-sensitive G proteins.[40] Here, we provide direct evidence for the ability of MT1 to engage the entire family of Gαi/o proteins. We also show that in the Gαq/11 family only Gα15 (but not Gαq or Gα11) can be directly activated. Notably, Gα15 is an atypical member of the Gαq/11 family that has been found to be activated by a large number of GPCRs.[41] The lack of observed coupling to Gαq and Gα11 is not due to a lack of sensitivity of the biosensors or the HEK293 cells used, since coupling to more conventional Gαq/11 family members (Gαq and Gα11) can be measured in HEK293 cells for other GPCRs.[42,43] Coupling of MT1 to Gα15 is likely to be of physiological relevance as MT1 expression has been reported in Gα15-expressing lymphocytes in which MLT promotes diacylglycerol production.[44] MLT-promoted activation of Gα15 appears to be specific to MT1, as it was not observed for MT2 in our previous study using the same biosensors.[5] However, given the high sequence similarity between MT1 and MT2 for cluster 2 residues, no conclusion on the determinants of Gα12 and Gα15 coupling could be drawn. Our study also reveals for the first time the coupling of MT1 to Gα12, a G protein linked to the activation of small G proteins of the Rho family and cytoskeleton modulation. Coupling to Gα12 appears to be specific to MT1, as it was also not observed for MT2.[5] Coupling of MT1 to Gα12 is consistent with previous studies reporting cytoskeleton changes upon stimulation of MT1.[45] In conclusion, our study shows that MT1 couples not only to Gαi/o proteins, like MT2, but shows a more expended profile by coupling in addition to Gα12 and Gα15 proteins.

Previous reports showed that frequent and rare variants in the MTNR1B gene encoding the MT2 melatonin receptor are associated with type 2 diabetes (T2D) risk.[23] The most extensively studied frequent rs10830963 variant, located in the unique intron, affects insulin secretion, most likely during the development of prediabetic fasting hyperglycemia.[46] Carriers of a rare variant, located in the coding region, display a number of adverse sleep, circadian, and caloric intake phenotypes that might participate in their T2D risk.[47] We therefore evaluated whether the genetic MT1 variants identified in our cohort of 8687 Europeans could be associated with increased risks of T2D or obesity. We were unable to detect any association either in children or in adults. The absence of association with obesity is consistent with previous observations that did not find strong support for a major role of MLT in body weight regulation and lipid metabolism in adult humans.[46] They are contrasting with rodent data showing that MLT has a beneficial effect on body weight regulation by limiting body weight gain and adipose tissue mass[48] as well as the phenotype of MT1 knockout mice that develop selective leptin resistance in the arcuate nucleus, an important hypothalamic center for body weight control.[49] Studies in rodents indicate also a general beneficial role of melatonin on glucose homeostasis[49] and MT1 knockout mice show systemic insulin resistance, most likely because of decreased insulin sensitivity and impaired hepatic glucose production.[50,51] There exists some information in humans supporting the beneficial role of melatonin on glucose regulation, but the role of MT1 remains unknown (see ref (46) for review). Our data on nonsynonymous MTNR1A gene variants or on genetically loss-of-function variants (i.e., p.W93*, p.Y170*, c.184 + 1G > T, c.999_1014del/p.W333*; Table S2) do not support a major role of MT1 in T2D development. Further studies, including copy number variations, will be necessary to fully address this point, as the MTNR1A gene shows an exceptionally high number of gene duplications compared to other GPCR genes.[2] In conclusion, these observations indicate important species-specific differences concerning the role of MT1 in T2D and obesity development between humans and rodents. Existing data are compatible with nonredundant functions of both MLT receptors in humans with a predominant role of MT2 in the regulation in glucose homeostasis. Whether the variants significantly impacting G protein signaling and/or βarrestin recruitment could be associated with specific phenotypes in humans remains to be determined.

The signaling profiles obtained for 31 naturally occurring MT1 variants reveals different signaling signatures that allowed us to classify them in three different clusters: Cluster 1 with no effect in the signaling profiles; Cluster 2 showing a major reduction in βarrestin-2 recruitment and different Gα protein activations, in particular Gα12 and/or Gα15; and Cluster 3 with severe impairing of overall signaling. Interestingly, EA scores matched the experimental results, which cross-validates the measured signaling profiles and the predictive power of evolutionary-based approaches.

The analysis of the location of these different clusters in three-dimensional structural models of active MT1 allowed us to propose specific mechanisms for the effect of the variations on receptor activation and signaling. The effects in Cluster 2 can be explained in most of the cases by alterations in the structure/dynamics of H8 that impair receptor-βarrestin-2 interactions. On the other hand, each variant in Cluster 3 is located near key regions for GPCR activation (the helical ICL2, the D/NRY motif, and the PIF motif); these variations may thus alter the signaling profile by altering the environment of these activation switches and impairing their function.

Interestingly, we found specific mutations that have a greater impact on βarrestin-2 recruitment, exemplifying how distinct receptor domains and allosteric interaction networks are involved in functional selectivity, i.e., the selective engagement of distinct signaling pathways. Remarkably, we found natural MT1 variants with signaling defects in apparently healthy people. However, the lack of detectable disease association may be due to the rare frequency of the variants that prevent robust statistical analysis. Alternatively, as MLT can act through two receptors—MT1 and MT2—one could think that MT2 could compensate for any loss of MT1 activity. However, it has been recently shown that MT2 can only activate Gαi/o family members and to recruit βarrestin-2.[5] Thus, such compensation from MT2 would only lead to partial signaling in response to MLT production. Most striking defects would be expected in individuals carrying MT1 defective variants in those cells/tissues in which MT2 is unable to compensate the MT1 defect, either because of the absence of MT2 expression or because of a MT1-specific function. Prominent examples are lymphocytes and MCF-7 human breast cancer cells, both expressing exclusively MT1 receptors. In addition, differential roles of MT1 and MT2 in sleep and circadian rhythm regulation have been suggested.[16] Based on our observation that several variant carriers show a complete loss of βarrestin-2 recruitment to MT1, the study of such variants in these biological system should be informative on the role of MT1-mediated βarrestin-2 signaling.

In conclusion, we have shown that sequencing and the functional characterization of rare variants combined to their computational analysis in relation with structural modeling represent valuable tools to identify key functional regions or residues in a receptor and allow prediction and rationalization of the effect of mutations on functional selectivity.

Methods

Study Participants

We investigated 8687 blood DNA samples from several population studies: (1) 3443 participants from the D.E.S.I.R. 9-year prospective study including middle-aged men and women from western France;[52] (2) 2642 participants who were recruited and followed up either by the CNRS UMR8199 (Lille, France), by the Department of Nutrition of Hotel-Dieu Hospital (Paris, France), or by the Centre d’Etude du Polymorphisme Humain (CEPH, Saint-Louis Hospital, Paris, France);[53,54] (3) 1567 participants who were recruited from the Department of Endocrinology of the Corbeil-Essonnes Hospital (Corbeil-Essonnes, France);[53] (4) 631 participants from the French Haguenau regional cohort study;[55] (5) 404 participants from the French Fleurbaix-Laventie Ville Santé study.[56] The study protocols were approved by local ethics committees. All participants older than 18 years signed an informed consent form. Oral assent from children or adolescents was obtained and parents (or legal guardian) signed an informed consent form.

In participants older than 18 years, obesity was defined as body mass index (BMI) ≥ 30 kg.m–2, while normal weight was defined as BMI < 25 kg.m–2. In children and adolescents younger than 18 years old, obesity was defined as BMI-for-age ≥95th percentile, while normal weight was defined as BMI-for-age <85th percentile according to the Centers for Disease Control and Prevention (CDC) growth charts. Type 2 diabetes (T2D) was defined as fasting plasma glucose ≥7.0 mmol L–1, and/or use of hyperglycemia treatment; while normal glucose was defined as fasting plasma glucose <6.1 mmol L–1 without hyperglycemia treatment.

Next-Generation Sequencing, Detection of MTNR1A Coding Variants, and Their Association Analyses

NGS targeting the coding exons of MTNR1A (NM_005958.4) was done as previously described.[57] In brief, target enrichment was performed according to the manufacturer’s protocol (NimbleGen SeqCap EZ) for Illumina sequencing. Sequence reads were mapped to the human genome (hg19/GRCh37) using Burrows-Wheeler Aligner (version 0.7.13). The variant calling was performed using Genome Analysis ToolKit (GATK; version 3.3). The annotation of variants was performed using the Ensembl Perl Application Program Interfaces (version 75). In the participants, all MTNR1A coding variants had a coverage higher than 8 reads and a QUAL score higher than 50. No variant had more than 5% of missing genotype (i.e., with a depth of coverage below 8 reads or a QUAL score below 50) across the participants, and no participant had more than 5% of missing genotypes (i.e., with a depth of coverage below 8 reads or a QUAL score below 50) across MTNR1A. The location of each variant was determined relative to the translation initiation codon using Human Genome Variation Society (HGVS) nomenclature for the description of variants. Among adults, we performed a case-control study for obesity adjusted for age, sex, the presence of T2D and genetic ancestry (controls were normal-weight adults), and a case-control study for T2D adjusted for age, sex, BMI, and genetic ancestry (controls were normal glucose adults with age above 40). Among children and adolescents, we performed a case-control study for obesity adjusted for age and sex (controls were normal-weight children/adolescents). Ancestry assessment was done as previously described.[57] The rare variants (with a MAF below 1%) were analyzed as a single cluster using MiST. MiST provides a score statistic S(π) for the mean effect (π) of the cluster, and a score statistic S(τ) for the heterogeneous effect (τ) of the cluster. The overall p-value is the combined p-values Pπ and Pτ, from the Fisher’s procedure. The frequent variants (with a MAF above 1%) were separately analyzed using a standard logistic regression adjusted for covariates defined above.

Cell Culture and Transfection

HEK293T cells were grown in complete medium (Dulbecco’s modified Eagle’s medium supplemented with 10% fetal bovine serum, 4.5 g·L–1 glucose, 100 U mL–1 penicillin, 0.1 mg mL–1 streptomycin, and 1 mM glutamine) (Invitrogen, CA). Transient transfections were performed using PEI (Sigma-Aldrich), according to manufacturer’s instructions.

Analysis of Cell Surface Receptor Expression

ELISAs were performed using the mouse antiflag antibody (Sigma-Aldrich, cat# F3165, RRID: AB_259529) and a horseradish peroxidase-conjugated mouse immunoglobulin G (IgG) whole-sheep secondary antibody (Sigma-Aldrich, cat# A6782, RRID: AB_258315).

Bioluminescence Resonance Energy Transfer (BRET) Measurement

Gαi1/i2/i3/oA/oB activations was assessed by measuring change of BRET between RlucII-Gαi1/i2/i3/oA/oB and GFP10-Gγ2 upon treatment of HEK293 cells cotransfected to express receptor and BRET sensors with MLT (Sigma-Aldrich).[5,19] For RlucII-Gαs, RlucII-Gαz, RlucII-Gαq, and RlucII-Gα12, we used the appropriate RlucII-tagged Gα subunits,[5] GFP10-Gγ1 was used given that optimal responsive window is generally obtained using these BRET sensor pairs for these Gα proteins. In all cases, the BRET pairs were used in the presence of untagged Gβ1 protein. Gα15 activation by receptors was assessed using a unimolecular biosensor. Since MT1 stimulation by melatonin does not activate PKC in our HEK293T cells in the absence of Gα15 but does so in its presence, we used the PKC activation as a surrogate assay to measure Gα15 activation by MT1. PKC activation was assessed by measuring BRET with this unimolecular biosensor. The BRET2-based PKC sensor consists of GFP10, 2 phospho-sensing domains (FHA1 and FHA2), a long and flexible linker, 2 specific phospho-substrate sequences, a RlucII, and the diacyl glycerol (DAG)-binding “C1b” domain of PKCδ added to increase signal specificity. Once phosphorylated, the phospho-substrates interact with their respective phospho-sensing domain, which stabilizes a conformation characterized by a higher BRET signal. βarrestin-2 recruitment was assessed by measuring BRET between RlucII-βarrestin-2 and rGFP-CAAX (prenylation CAAX box of KRas) demonstrating the recruitment of the βarrestin-2 to the plasma membrane.[20] Transfected cells were plated in 96-well white microplates (100 μL per well; Greiner). Forty-eight hours after transfection, cells were washed with Dulbecco’s phosphate-buffered saline, and then Tyrode’s buffer was added. After the addition of MLT, the cells were incubated for 10 min at 37 °C, and coelenterazine 400a (2.5 μM; NanoLight Technology) was added 5 min before BRET reading in a Synergy Neo Microplate reader (BioTek) (acceptor filter: 515 ± 15 nm; donor filter: 410 ± 40 nm).

To monitor receptor internalization, ebBRET was measured between MT1-WT or variants (R54W, I309T, or C314R) fused to RlucII and the endosomal compartment marker rGFP-FYVE.[20] Prolume purple (2.5 μM; NanoLight Technology) is added for 6 min followed by MLT stimulation and readings are recorded on a Mithras (Berthold scientific) microplate reader (acceptor filter: 515 ± 10 nm; donor filter: 400 ± 35 nm).

Net BRET represents the BRET signal minus the signal in the absence of a GFP10-tagged sensor. ΔBRET refers to the difference in net BRET recorded from cells treated with agonist and cells treated with vehicle. For the agonist concentration–response curves, the responses for the variants were expressed as a % of the response observed for the WT receptor.

Evolutionary Action

The evolutionary action (EA) of a mutation on fitness is given by[26]where ϕ is the fitness phenotype, γ is the genotype sequence, and f is the evolutionary function between them: f(γ) = ϕ. Although the function f describes the fitness landscape and is unknown, its derivative is formally given by eq , the terms of which can be approximated from evolutionary information. Specifically, dy is the magnitude of mutation, which can be calculated in terms of log odds of amino acid substitutions in the case of a coding mutation, and ∇f is the gradient of the fitness landscape, which can be calculated from the functional sensitivity to mutations at each sequence positions, which can be calculated by the coupling between sequence variations and phylogenetic divergences.[58] In practice, the EA score is normalized by percentile rank from 0 to 100, where a variant is relatively more deleterious as the score close to 100. This approach tends to perform well against other predictors of variant effect in blind challenges.[59] To predict the effect of each variant within the MT1 receptor, EA score is calculated using MTNR1A multiple sequence alignment retrieved from the nonredundant protein sequence (nr) database (https://blast.ncbi.nlm.nih.gov/Blast.cgi).

Statistical Analysis of Functional Analyses

Data are presented as means ± SEM of indicated “n” independent experiments, performed at least in triplicate to ensure the reliability of single values. The level of probability used in this article is defined as follows: NS: nonsignificant, *p < 0.05, **p < 0.01, ***p < 0.001 relative to control.

LogEC50 and agonist-induced Emax values were obtained following nonlinear regression (curve fit) with four parameters of data from a minimum of eight different concentrations per experiment, repeated at least three times independently. All assays performed in this study were previously validated and demonstrated the robustness and variability of the procedure using this number of independent experiments for concentration–response curves. Data and statistical analysis were performed using GraphPad Prism software version 6.0.

For each assay, the agonist-induced Emax value for every MT1 variant was normalized as a percentage of the maximal MLT-stimulated response of the WT receptor (set at 100) monitored in parallel with the receptor variant. The dose–response curves were fitted to an operational model of agonism designed by Kenakin and Christopoulos[60,61] to obtain log(τ/KA) values for the WT receptor and its variants. Normally, an agonist is set as a reference agonist, against which within-pathway comparisons for the same receptor to other agonists can be made and expressed as Δlog(τ/KA). Here, within-pathway comparisons were made between MT1 variants and the WT receptor. Normalized difference was calculated on values corresponding to agonist-induced Emax, and logEC50 to fit a–1 to +1 scale using the following formula: (variant – WT)/(variant + WT). In the case of logEC50, before normalization, the antilogs were first calculated and then were fitted to the aforementioned formula. Positive and negative values represent mutations with better or worse responses, respectively, than those of the WT receptor. Subsequently, logEC50, agonist-induced Emax, and Δlog(τ/KA) values were expressed as means ± SEM of the indicated number of experiments (n). Statistical analysis for logEC50 and agonist-induced Emax was performed by comparing independent fits with a global fit that shared the selected parameter.

Computational Analysis

To quantify similarities between signaling perturbations resulting from each mutation, we adapted the NMF and K-means clustering process originally described in ref (6). We performed this analysis on an input matrix of 31 mutations and 24 signaling parameters spanning 8 signaling pathways (Gαi1/Gαi2/Gαi3/GαoA/GαoB/Gα12/Gα15/βarrestin-2) in response to melatonin. The mean and STD for each phenotypic measure were used to propagate experimental error through the clustering process. Explicitly, 1000 input matrices were generated by independently sampling from the experimental error (Mean + STD) for each phenotypic data point. Each of these sampled input matrices were then processed separately using the NMF/K-means method. The values for each phenotype were scaled to avoid bias introduced by differences in phenotype scale usingfor every mutation i and every parameter j, where maximumj and minimumj values were specific to each parameter column (e.g., EC50, Emax, τ/Ka). In this way, every parameter was normalized between 0 and 1. The dimensionality of these 1000 input matrices were then independently reduced using NMF to identify the key signaling patterns. K-means clustering was used, as described, to cluster the [mutation × feature] (31 × K) basis vector. As previously described, clustering results for each sampled matrix were averaged across K = 2 to K = 6, where K is the number of features for NMF and the number of clusters for K-means (KNMF = KK-means). Final cluster assignment, quantifying mutation similarity, was the average over all 1000 sampled matrices. Clustering frequency was converted into distance matrix using Pearson’s correlation and visualized as a heatmap and dendrogram.

Radioligand Binding Experiments

Radioligand saturation binding experiments were performed in transient transfected HEK293T cell lines. Saturation binding assays were carried out on crude membrane preparations with increasing concentrations of 2-[125I]MLT (PerkinElmer Life Sciences) to determine dissociation constants (Kd) and total melatonin receptor expression (Bmax) of wild-type and mutant MT1 receptors using PRISM software (GraphPad). Nonspecific binding was determined using 10 μM melatonin.

Three-Dimensional Modeling

The location of the mutant variants in the inactive state of MT1 was mapped onto the crystal structure of MT1 bound to ramelteon (PDBid: 6ME2).[25] The location of the mutant variants in the active state of MT1 was mapped onto a three-dimensional model of the receptor built by homology modeling using as templates (i) the cryo-EM structure of the human mu-opioid receptor (33% sequence similarity) in complex with Gαi (PDBid: 6DDE)[62] for the transmembrane bundle (TM1–TM7; residues 21–290), and (ii) the crystal structure of the human adenosine A2A receptor (26% sequence similarity) (PDBid: 2YDO)[63] for the last 18 residues of the model comprising ICL4 and helix 8 (residues 291–308). We chose A2A receptor as a template because, although it couples to Gαs, it has a high sequence similarity with MT1 in the transmembrane domain (see above) and, importantly, is phylogenetically close to it (the melatonin receptor cluster and the MECA receptor cluster are neighbors in the α group of class A GPCRs[64]). A structure-based sequence alignment between the templates and MT1 was obtained from the GPCRdb.[65] This initial alignment was manually refined using Chimera[66] to adjust the gaps in the loop regions. Using this alignment and the structural templates, 3D models of MT1 in an active conformations were built using Modeler ver. 9.16.[67] The cysteine bridge between p.C1003.25 in TM3 and p.C177ecl2.50 in the second ECL was explicitly defined during model building. All models were subjected to 300 iterations of variable target function method (VTFM) optimization and thorough (md_level = refine.slow) molecular dynamics and simulated annealing optimization within Modeler, and scored using the discrete optimized protein energy potential.[67] The 20 best-scoring models were analyzed visually, and a suitable model (in term of low score and structure of the loops) was selected.

To visualize the arrangement of the MT1/Gαi and MT1/arrestin complexes, the structures of Gαi (PDBid: 6DDE)[62] and visual βarrestin-1 (PDBid: 5W0P)[68] were structurally aligned onto the model of active MT1 (using the receptor transmembrane domains for the fit) with molecular graphics software (The PyMOL Molecular Graphics System, Version 2.2.3, Schrödinger, LLC.).