Fusion with Promiscuous Gα16 Subunit Reveals Signaling Bias at Muscarinic Receptors

A complex evaluation of agonist bias at G-protein coupled receptors at the level of G-protein classes and isoforms including non-preferential ones is essential for advanced agonist screening and drug development. Molecular crosstalk in downstream signaling and a lack of sufficiently sensitive and selective methods to study direct coupling with G-protein of interest complicates this analysis. We performed binding and functional analysis of 11 structurally different agonists on prepared fusion proteins of individual subtypes of muscarinic receptors and non-canonical promiscuous α-subunit of G16 protein to study agonist bias. We have demonstrated that fusion of muscarinic receptors with Gα16 limits access of other competitive Gα subunits to the receptor, and thus enables us to study activation of Gα16 mediated pathway more specifically. Our data demonstrated agonist-specific activation of G16 pathway among individual subtypes of muscarinic receptors and revealed signaling bias of oxotremorine towards Gα16 pathway at the M2 receptor and at the same time impaired Gα16 signaling of iperoxo at M5 receptors. Our data have shown that fusion proteins of muscarinic receptors with α-subunit of G-proteins can serve as a suitable tool for studying agonist bias, especially at non-preferential pathways.


Introduction
G-protein-coupled receptors (GPCRs) are the largest family of human membrane proteins that transmit signals into a cell through heterotrimeric G-proteins. GPCRs represent the primary target for drug development with potential application in essentially all clinical fields. They mediate a broad range of physiological processes by driving multiple intracellular effectors through various classes of G-proteins. Individual GPCRs preferentially couple to the particular class of G-proteins but they can also successfully activate others [1][2][3]. This coupling promiscuity was observed in both artificial systems with over-expressed GPCRs and native cells [4,5]. Besides G-proteins, GPCR can couple with β-arrestins which desensitize and scaffold G-protein-driven signaling pathways [6]. The multiplicity of signaling leads to the high complexity of the functional response of GPCRs to agonist stimulation.
Structurally different agonists induce specific changes in the GPCRs leading to stabilization of agonist-specific conformations that can lead to non-uniform agonist-specific modulation of signaling pathways. This preferential orientation of the signaling of a GPCR towards a subset of its signal transducers is termed signaling bias [7]. An agonist biased to a particular G-protein pathway may promote therapeutically desired signaling while simultaneously avoiding side effects mediated by activation of others, especially in conditions with well-understood pathophysiology [8][9][10]. For example, melanocortin receptor 4 (MC4R) agonist melanotan II produces its anorectic effects through coupling to G q/11 and its adverse cardiovascular effects through G s coupling, suggesting potential therapeutic benefit in obesity for G q/11 -biased ligands [11].
The accurate evaluation of agonist bias regarding individual G-protein pathways is crucial for preclinical drug development. However, it is a difficult task given the high complexity of GPCRs signaling. The molecular crosstalk that can occur among downstream effector molecules may bring in further complexity [12]. One of the most challenging tasks is to develop a suitable technique for the analysis of the signaling pathway of interest with high sensitivity and sufficient selectivity that is free from the interference of other signaling pathways. Measurement of second messengers struggles with molecular crosstalk of signaling pathways. The analysis of coupling of GPCRs with individual G-proteins is difficult due to the presence of others that interact concurrently with a given signaling pathway, especially in studies of non-preferential signaling pathways [13] or in studies of signaling pathways mediated by individual isoforms of given G-protein.
Muscarinic signaling is implicated in numerous pathologic events, such as the promotion of carcinoma cell growth, early pathogenesis of neurodegenerative diseases in the central nervous system of Alzheimer's and Parkinson's, schizophrenia, drug addiction, pain, and also in some internal diseases, e.g., asthma or overactive bladder [14,15]. As of now, no affinity-based selective agonists of individual muscarinic receptors have been discovered, due to the high homology of the orthosteric binding site among individual muscarinic subtypes [16][17][18].
Selective targeting on the G i/o versus G q/11 mediated pathway by biased agonist could be a way to achieve selectivity to even or odd muscarinic subtypes [19]. Moreover, agonists biased to individual isoforms of G-proteins could lead to tissue-specific activation of mAChRs, due to the predominant expression of some G-proteins in specific tissues (e.g., G o in the central nervous system or G 16 in hematopoietic cells) [10,20]. We have focused on variations in the G 16 signaling pathway that was not studied so far, is rare and may lead to very specific effects (e.g., tissue-specific activation). We have analyzed variation in the G 16 signaling profile among individual subtypes of muscarinic receptors.
To reveal and properly quantify putative agonist bias to certain G-proteins and their isoforms, especially non-preferential ones, among individual subtypes of muscarinic receptors, a system that is sufficiently sensitive and specific is required. Furthermore, 1:1 Gα-receptor stoichiometry would simplify the analysis and interpretation of found agonist bias. We assume that fusion proteins of muscarinic receptors with Gα subunit of interest could serve as a convenient tool to screen agonist bias towards the particular Gα among individual subtypes of muscarinic receptors. Importantly, we expect that tight fusion of a receptor with a particular Gα prevents the coupling of other competing Gα to the receptor. If so, the signaling of a pathway of interest can be selectively analyzed. Fusion proteins of GPCR and α-subunit of G-protein were used to study individual G-protein mediated pathways in several studies [21][22][23][24][25][26]. We validate our assumptions on an example of fusion proteins of individual muscarinic receptors and non-canonical promiscuous Gα 16 subunit.
Muscarinic signaling via non-canonical G 16 G-protein may play a relevant physiological role. At the protein level, G 16 expression is only detected in highly specific cell types (hematopoietic and epithelial cells) characterized by a high rate of cell turnover [27]. G 16 mediated signaling may play role in immune response [28] and tumor cell growth [29].
Promiscuous Gα 16 efficiently couples to any subtype of muscarinic receptor. That leads to the phospholipase C-activation, resulting in the formation of inositol phosphates. We performed a binding and functional analysis of these constructs using eleven structurally different muscarinic agonists. We demonstrated agonist-specific activation of non-canonical G 16 pathway varying among individual subtypes. Additionally, we compared the signaling of agonists oxotremorine and iperoxo at Gα 16 _fused, Gα 16 co-transfected, and wild types of M 2 and M 5 muscarinic receptors and revealed signaling bias of oxotremorine towards Gα 16 pathway at the M 2 receptor and at the same time impaired Gα 16 signaling of iperoxo at M 5 receptors.

Description of Fusion Proteins
Fusion proteins (denoted M 1 _Gα 16 through M 5 _Gα 16 ) were constructed from individual subtypes of muscarinic receptors M 1 -M 5 and α-subunit of G 16 G-protein. The α-subunit was tightly connected to the C-terminus of the respective receptors as described in the Methods. Palmitoylation sites at helix 8 of receptors, as well as at N-terminus of Gα 16 , were preserved to ensure their anchoring to the membrane. Complete sequences of fusion proteins are shown in Supplementary Materials.

Homology Models of Fusion Proteins
To test whether fusion proteins respect the natural arrangement of the receptor and G-protein α-subunit that allows their successful coupling with no serious occurring disturbance to the structure arrangement, we have built homology models of M 1 _Gα 16 and M 2 _Gα 16 fusion proteins. Homology modeling resulted in a good model free of unusual structural features. Overlays of fusion proteins with cryo-EM of M 1 + Gα 11 (6OIJ) and M 2 + Gα o (6OIK) receptor-G-protein complexes [30] are shown in Figure 1. The stability of structure was verified by running molecular dynamics (MD) of fusion proteins in complex with G-protein βγ-dimer in membrane/water system. Analysis of MD trajectories by Simulation Quality Analysis tools of Maestro confirmed the stability of the structures (Supplementary Materials Figure S3). No structural rearrangements occurred during 120 ns of MD. Insertion of the C-terminus of Gα 16 to G-protein binding site at the receptor located between transmembrane helix 3 and 6 corresponds to the insertion of Gα 11 at M 1 and Gα o at M 2 . At the M 1 receptor, the position of the C-terminus of Gα 11 and Gα 16 are practically identical. However, at M 2 , Gα 16 is inserted under a sharper angle than Gα o . Insertion of the C-terminus of Gα 16 instead of Gα 11 or Gα o to the G-protein binding site did not induce any major change in the receptor conformation.  Table S1).

Lack of Coupling of Gα 16 _Fused Receptors with Endogenous G-Proteins
Muscarinic receptors are able to activate multiple G-proteins. Preferentially, muscarinic receptors M 1 , M 3 , and M 5 couple with Gα q/11 and M 2 and M 4 receptors with Gα i/o G-proteins. All muscarinic subtypes efficiently activate non-canonical promiscuous G-protein (G 16 ) followed by activation of the appropriate signaling pathway (phospholipase C-activation and generation of IPX). Based on molecular models, we expected that fusion of muscarinic receptors with Gα 16 subunit would sterically prevent coupling of endogenously expressed G-proteins. To this end, we analyzed changes in cAMP level, mediated by endogenous G i/o and G s proteins, after activation of wt and Gα 16 _fused M 2 and M 4 receptors by agonist carbachol. Basal level of cAMP was determined in presence of 10 µM adenylate cyclase activator forskolin. Values of basal level determined as % of incorporated radioactivity varied in the range of 2.5-3% and are the same in cells expressing wt and fused receptors. Level of cAMP was calculated as fold over basal level. We demonstrate that tight fusion with Gα 16 prevented the coupling of preferential Gα i/o and non-preferential G s to M 2 and M 4 receptors. While carbachol stimulated accumulation of IP X at all Gα 16 _fused receptors (Supplementary Materials Table S3), at M 2 _Gα 16 and M 4 _Gα 16 , did not change the level of cAMP, whereas at wt M 2 and M 4 , carbachol inhibited cAMP synthesis via preferential G i/o G-proteins at submicromolar concentrations and stimulated it via non-preferential G s G-proteins at micromolar concentrations ( Figure 2). Thus, the fusion proteins pass the signal solely through the fused Gα subunit.

Binding and Functional Analysis of Gα 16 Fused Receptors
Eleven structurally different agonists, varying in the binding mode to muscarinic receptors, potency, and efficacy to activate muscarinic receptors (arecoline, carbachol furmethide, iperoxo, McN-A343, N-desmethylclozapine, oxotremorine, pilocarpine, xanomeline, JR-6, and JR-7), were used for pharmacological evaluation of the fusion proteins. Structures of tested agonists are shown in Supplementary Materials Figure S2 The affinity of tested agonists to fused proteins was assayed in competition experiments with 1nM [ 3 H]NMS, calculated according to Equation (4), and is summarized in Table 1. All agonists completely inhibited [ 3 H]NMS binding to fused proteins. All tested agonists displayed only low-affinity binding, except for iperoxo at M 1 _Gα 16 and JR6 at M 4 _Gα 16 . Affinities of low-affinity binding of tested agonists (carbachol, oxotremorine, pilocarpine, JR6, and JR7) to wt and Gα 16 _fused muscarinic receptors were compared. Data are summarized in Supplementary Materials Table S2. The affinity of carbachol was slightly lower at all Gα 16 _fused receptors than at corresponding wt. The decrease in affinity was observed also for pilocarpine, JR7, and oxotremorine (except oxotremorine at M 1 and JR7 at M 2 ). On the other hand, JR6 had a higher affinity at all Gα 16 _fused receptors, especially at M 2 _Gα 16 affinity of JR6 was 23-times higher than at wt M 2 .
by endogenous Gi/o and Gs proteins, after activation of wt and Gα16_fused M2 and M4 receptors by agonist carbachol. Basal level of cAMP was determined in presence of 10 µM adenylate cyclase activator forskolin. Values of basal level determined as % of incorporated radioactivity varied in the range of 2.5-3% and are the same in cells expressing wt and fused receptors. Level of cAMP was calculated as fold over basal level. We demonstrate that tight fusion with Gα16 prevented the coupling of preferential Gαi/o and nonpreferential Gs to M2 and M4 receptors. While carbachol stimulated accumulation of IPX at all Gα16_fused receptors (SM Table S3), at M2_Gα16 and M4_Gα16, did not change the level of cAMP, whereas at wt M2 and M4, carbachol inhibited cAMP synthesis via preferential Gi/o G-proteins at submicromolar concentrations and stimulated it via non-preferential Gs G-proteins at micromolar concentrations ( Figure 2). Thus, the fusion proteins pass the signal solely through the fused Gα subunit.

Binding and Functional Analysis of Gα16 Fused Receptors
Eleven structurally different agonists, varying in the binding mode to muscarinic receptors, potency, and efficacy to activate muscarinic receptors (arecoline, carbachol furmethide, iperoxo, McN-A343, N-desmethylclozapine, oxotremorine, pilocarpine, xanomeline, JR-6, and JR-7), were used for pharmacological evaluation of the fusion proteins. Structures of tested agonists are shown in Supplementary Materials (SM) Figure S2

Functional Response of Gα 16 _Fused Muscarinic Receptors to Agonists
The fusion with non-canonical promiscuous G-protein Gα 16 couples all subtypes of muscarinic receptors to phospholipase C-activation and generation of IP X . The level of IP X was measured by radio-chromatographic separation. Basal level (in absence of agonist) varied in range of 2-3% of incorporated radioactivity and was the same in cells expressing wt and individual fused receptors. Level of IP X in presence of individual concentrations of tested agonists was calculated as fold over basal level. Parameters of accumulation of IP X as a functional response of fused proteins to stimulation by a tested agonist, EC 50 and E MAX , are summarized in (Supplementary Materials Table S3). To calculate the coefficient of operational efficacy τ of functional response of individual Gα 16 _fused receptors to tested agonists, the system E MAX was determined from functional responses to the agonists carbachol, oxotremorine, and pilocarpine according to the procedure described recently [31]. Then, the τ value was used for the calculation of the equilibrium dissociation constant K A . The values of τ and K A calculated according to Equation (6) are summarized in Table 2.
G i/o -biased muscarinic partial agonists JR6 and JR7 did not stimulate the accumulation of IP X at any fused protein. Although fusion with Gα 16 led to an increase in the affinity of JR6 to all fusion proteins, JR6 and JR7 induced conformation incompatible with activation of the Gα 16 signaling pathway. Except for JR6 and JR7, all tested agonists stimulated accumulation of IPx at all Gα 16 _fused receptors.
Quantification of Agonist Bias towards Individual Gα 16 Fused Receptors.
To compare agonist specific activation of G 16 mediated pathway among individual muscarinic subtypes and to quantify agonist bias towards individual Gα 16 _fused receptors, intrinsic activities relative to carbachol (RA i ) were calculated according to Equation (8) from the E MAX and EC 50 values of accumulation of inositol phosphates (Supplementary Materials Table S3). Values of RA i are summarized in Table 2 and plotted in Figure 3. Interestingly, M 2 super-agonist iperoxo [32,33] Table S4, Figure S1). The variability in bias among agonists eliminates the possibility that protein fusion introduced a bias towards some of the receptors.
We analyzed activation of IP X pathway by agonist carbachol, oxotremorine, and iperoxo at Gα 16 _fused receptors, receptors co-transfected with Gα 16 -subunit and wt M 2 and M 5 receptor ( Figure 4, Table 3).  Figure S1). The variability in bias among agonists eliminates the possibility that protein fusion introduced a bias towards some of the receptors.  subunits to M2 may take place) which indicates bias of oxotremorine towards Gα16 mediated pathway at M2 receptor. Interestingly, agonist iperoxo had higher operational efficacy τ than carbachol at all Gα16_fused receptors, except M5_Gα16 (Table 2). At fused M5_Gα16, τ value of iperoxo was almost 30% lower than τ values of carbachol. At the rest of Gα16_fused receptors, τ values of iperoxo were greater than τ for carbachol, least by 40 % (M4) and most nearly 7-fold (at M2). In contrast to M5_Gα16, iperoxo stimulated accumulation of IPX at M5-wt with higher operational efficacy than carbachol. At co-transfected system M5 + Gα16, iperoxo and carbachol stimulated IPX accumulation with comparable efficacy (Figure 4, Table 3) which indicates impairment of Gα16 signaling of super-agonist iperoxo at fused M5 receptor.   Agonist induced coupling with Gα 16 : Data show better coupling of Gα 16 fused M 2 receptor in comparison to the co-transfected system for carbachol and oxotremorine. At IP X pathway, the value of equilibrium dissociation constant expressed as the negative logarithm (pK A ) for reference agonist carbachol as well as tested agonist oxotremorine, was higher at Gα 16 _fused receptors than at M 2 receptors co-transfected with Gα 16 (Table 3), indicating a better coupling in the case of the fusion protein. The better coupling to G 16 at fused receptors M 1 _Gα 16 through M 5 _Gα 16 than at co-transfected variants is obvious also from comparison of ( Table 2 with our previous data Randakova et al. [28]). The pK A value for reference agonist carbachol and oxotremorine, as well as pilocarpine or xanomeline, was higher at Gα 16 _fused receptors than at corresponding wt receptors co-transfected with Gα 16. The increase in pK A ranged from 3-fold for xanomeline at M 1 to 63-fold for oxotremorine at M 3 . Table 3. Comparison of parameters of functional response of variants of M 2 and M 5 receptor. Parameters of agonist-induced functional response EC 50 and E MAX were obtained by fitting Equation (6) to data from measurement of the accumulation of inositol phosphates. Operational efficacy τ, agonist equilibrium dissociation constant K A and agonist relative intrinsic activity RA i were calculated according to Equations (6)-(8), respectively. EC 50  In contrast, pK A of iperoxo at the fused M 5_ Gα 16 was lower than at co-transfected M 5 +Gα 16 , indicating worse coupling of the fusion protein ( Table 3). The high variability in the observed shift in pK A excludes a possibility of the systemic artifact caused by protein fusion.
Comparison of operational efficacies of selected agonists: In comparison to the cotransfected system M 2 +Gα 16 , the increase in operational efficacy τ of both reference agonist carbachol as well as tested agonist oxotremorine to stimulate the non-canonical accumulation of IP X at M 2 _G α16 (Table 3) indicates the higher sensitivity of measurement of functional response at Gα_fused receptors. Oxotremorine had higher operational efficacy than carbachol at M 2 _Gα 16 and M 5 _Gα 16 (Table 2). At the rest of the Gα 16 _fused receptors, the operational efficacies of oxotremorine and carbachol were the same. In contrast, oxotremorine stimulated accumulation of IP X at M 2 + Gα 16 and M 2 wt with efficacy comparable (Table 3) or lower [19] to carbachol. Operational efficacies τ of functional responses of carbachol and oxotremorine at M 2 _Gα 16 and M 2 + Gα 16 (Figure 4) are summarized in Table 3. In other words, at M 2 _G α16 fusion protein (where M 2 receptor signals only via G α16 ) oxotremorine had higher efficacy than in co-transfected system (where binding of other Gα subunits to M 2 may take place) which indicates bias of oxotremorine towards Gα 16 mediated pathway at M 2 receptor.
Interestingly, agonist iperoxo had higher operational efficacy τ than carbachol at all Gα 16 _fused receptors, except M 5 _Gα 16 (Table 2). At fused M 5 _Gα 16 , τ value of iperoxo was almost 30% lower than τ values of carbachol. At the rest of Gα 16 _fused receptors, τ values of iperoxo were greater than τ for carbachol, least by 40 % (M 4 ) and most nearly 7-fold (at M 2 ). In contrast to M 5 _Gα 16 , iperoxo stimulated accumulation of IP X at M 5 -wt with higher operational efficacy than carbachol. At co-transfected system M 5 + Gα 16 , iperoxo and carbachol stimulated IP X accumulation with comparable efficacy (Figure 4, Table 3) which indicates impairment of Gα 16 signaling of super-agonist iperoxo at fused M 5 receptor.

Discussion
In this study, we show that fusion proteins of receptor and α-subunit of G-protein are a suitable tool for studying agonist bias. We demonstrate it on the example of muscarinic receptors fused with Gα 16 subunit and 11 muscarinic agonists whose signaling profile (bias) varies among receptor subtypes.
Analysis of signaling bias of muscarinic receptors, concerning G-protein mediated signaling, has several pitfalls. Coupling promiscuity of muscarinic receptors leads to molecular crosstalk in downstream signaling. For example, calcium ions released upon activation of G q/11 IP X pathway modulate some adenylate cyclases and thus cAMP signaling. In turn, βγ-dimers released from G i/o G-proteins modulate some calcium channels and thus calcium signaling [36,37]. Moreover, signals of non-preferential pathways are usually weak, thus, highly sensitive methods are needed. The main obstacle, in the study of the non-preferential G-protein pathways, is the competition of different (mainly preferential) Gα-subunits for the binding site at a given receptor. Activation of a non-preferential G-protein pathway may play important roles in processes characterized by fluctuation in an expression of individual G-proteins or GPCRs, e.g., immune cell maturation [28], progression of cancer [38], or Parkinson s disease [39].
Several tools including G-protein-specific pharmacological inhibitors or toxins [40], C-terminus mimicking peptides [41], small interfering RNA [42,43], using artificial systems with limited endogenous G-proteins [44][45][46] or reconstitution of purified receptors and G-proteins in the artificial membrane [47,48] limit the signal mediated by certain G-proteins. Techniques like the immunoprecipitation with specific Gα antibodies [2,49], resonance energy transfer techniques, where bioluminescent (BRET) or fluorescent (FRET) donors and acceptors are fused on the C-terminus of the GPCR and in one of the subunits of the G-protein [50,51] were used to study specific GPCR-G-protein interactions. Although these methods diminish or eliminate signaling crosstalk, they are not aimed at high sensitivity.
Receptor_Gα fusion proteins are well described to study the activation of individual Gprotein mediated signaling pathways at many GPCR [21][22][23][24][25][26]. We demonstrate their use to study agonist bias at non-canonical G 16 pathway among individual subtypes of muscarinic receptors. Gα 16(15) is expressed only in highly specific cell types such as hematopoietic and epithelial cells [27], which are characterized by a high rate of cell turnover. Muscarinic receptors expressed in these cells appear to be involved in the regulation of diverse cellular activities including immune response [28], cell proliferation, or cell differentiation [52,53].
The engineering of receptor-transducer fusion proteins seems to be an effective strategy to target cellular effectors more efficiently and specifically [21]. Fusion proteins enable the study of signaling mediated by G-proteins up to the level of individual G-proteins isoforms. Moreover, receptor-G-protein fusion forces a 1:1 stoichiometry and ensures efficient coupling of the given receptor to an attached Gα subunit. Receptor-G-protein stoichiometry is a relevant aspect of signaling bias and should be taken into account in the screening of biased agonists [54].
We have created fusion proteins of individual muscarinic receptors (M 1 -M 5 ) and non-canonical promiscuous Gα 16 subunit and performed detailed binding and functional analysis of these constructs using 11 structurally different muscarinic agonists to evaluate the suitability of such fusion proteins to study agonist bias. Structurally different agonists vary in interactions in the orthosteric binding site of the muscarinic receptor [55]. The portfolio of used agonists included reference balanced full agonist carbachol, classic muscarinic agonists arecoline, furmethide, pilocarpine, oxotremorine, super-agonist iperoxo [32,33], bitopic agonists xanomeline [56], and McN-A343 [57], and Gi/o-biased agonists JR6 and JR7 [19] (Supplementary Materials Figure S2).
The use of Gα_fused receptors for analysis of signaling bias is conditioned by the full preservation of binding and functional properties of both receptor and Gα subunit. In the preparation of the construct, palmitoylation sites, at the C-terminus of the receptors [58], and the N-terminus of the Gα subunit [59], that mediate interaction with the membrane, were maintained. That is essential for keeping the native conformation of a receptor as well as G-protein. Comparison of homology models of prepared constructs M 1 _Gα 16 and M 2 _Gα 16 with cryo-EM structures of receptor-G-protein complexes M 1 + Gα 11 and M 2 + Gα o [30] confirmed the natural arrangement of the receptor and Gα in fusion proteins ( Figure 1). At the M 1 receptor, insertion of the C-terminus of related Gα 16 and Gα 11 subunits into the G-protein binding site of the receptor is practically identical. On the other hand, at the M 2 receptor, evolutionarily more distant Gα 16 and Gα o differ in the angle at which they are inserted into the G-protein binding site. Gα-specific insertion of C-terminus into the intracellular cavity of cognate GPCR was observed in 3D structures of GPCR-G-protein complexes [30,[60][61][62][63] and demonstrated using molecular dynamics (MD) as well [64]. Furthermore, the fusion of muscarinic receptors with Gα 16 subunit did not affect the binding affinity of the labeled antagonist [ 3 H]N-methylscopolamine at any fusion protein (Supplementary Materials Table S1), indicating that fusion did not markedly influence receptor conformation.
The signaling of interest can be selectively analyzed when the binding of other competing G-proteins to the receptor is excluded. We hypothesized that tight fusion of the receptor with a particular Gα subunit prevents the binding of other G-proteins. The M 2 and M 4 receptors preferentially inhibit cAMP synthesis via Gα i/o G-proteins and can also couple with non-preferential Gα s to activate cAMP synthesis [65,66]. In contrast to the wt M 2 and M 4 receptors (Figure 2), carbachol did not induce changes in cAMP level at fused M 2 _Gα 16 and M 4 _Gα 16 receptors, indicating no coupling to endogenous G i/o or G s G-proteins. It suggests that, unlike some fusion constructs [67], our directly Gα 16 _fused constructs indeed prevent the access of competitive Gα subunits to the receptor.
The binding analysis has shown that in contrast to wild-type (wt) receptors, at Gα 16 _fused constructs, almost all tested agonists displayed only low-affinity binding. Gprotein binding to a receptor might, in turn, allosterically influence ligand binding [68,69]. The absence of high-affinity binding of most agonists to Gα 16 _fused receptors may be either due to lack of pre-coupling of Gα 16 to the receptor or receptor is pre-coupled Gα 16 that binds GDP [43]. Since the decrease in the value of equilibrium dissociation constant (K A ) of agonists at Gα 16_ fused receptors in comparison to wt receptor co-expressed with Gα 16 (Table 3, Table 2 versus our previous data [28]) indicates pre-coupling, the absence of high-affinity binding indicates pre-coupling Gα 16 that binds GDP [43]. G 16 G-protein is efficiently capable to couple all muscarinic subtypes (M 1 -M 5 ) via phospholipase C activation (IP X accumulation). Thus, it may be possible to analyze the activation of all muscarinic subtypes using one assay (measurement of the accumulation of inositol phosphates, IP X ) and demonstrate agonist-specific activation of this pathway. Signaling bias among individual Gα 16 _fused receptors was calculated from relative intrinsic activities RA i to reference agonist carbachol [70] (Table 2). RAi values can be easily calculated for several pathways and many ligands and quickly compared. In principle, for a single signaling pathway and two or more receptors, a ligand that has greater RAi at one receptor than at other(s) is biased to a given pathway at that receptor. Additionally, we analyzed our data also by conventionally used bias factor [35]. Data are summarized in (Supplementary Materials Table S4) and plotted (Supplementary Materials Figure S1). Quantification of agonist bias obtained by both ways was the same, showing that a quick comparison of RAi factors is sufficient and that analysis was conducted correctly. Presented data demonstrate differences in the pattern of the Gα 16 pathway activation at five subtypes of Gα 16 _fused muscarinic receptors after stimulation by structurally different agonists. It points to variations in the compatibility of agonist-specific conformations with Gα 16 coupling and activation. While some agonists have quite balanced Gα 16 pathway activation patterns-such as pilocarpine, oxotremorine, or xanomeline-profound bias towards individual Gα 16 _fused muscarinic receptors was observed for agonists McN-A-343 (towards M 2 _Gα 16 ) and iperoxo (towards M 3 _Gα 16 ). McN-A343 is a bitopic agonist, capable of stimulating the G q pathway while incapable of stimulating G s at M 1 expressed in CHO cells [64]. We have shown that McN-A343 successfully activates G 16 pathway at all muscarinic subtypes with a bias towards M 2 . Interestingly, M 2 super-agonist iperoxo displayed bias towards M 3 _Gα 16 over other Gα 16 _fused receptors. It was demonstrated that iperoxo-based dualsteric compounds exert bias G i/o over G s pathway at M 2 [71] but exert bias to G q/11 over G i/o signaling at the M 1 receptor [72]. On the other hand, M 1 -preferring agonist N-desmethylclozapine displayed bias to M 1 _Gα 16 and M 3 _Gα 16 over the rest of the subtypes. It points to huge variability in signaling depending on the combination of a ligand-receptor-pathway system, promising a chance to find agonists with a bias to the desired pathway at the desired receptor subtype.
Comparison of parameters of functional response of selected agonists at Gα 16 _fused and Gα 16 cotransfected wt receptors suggest better coupling of fused Gα subunit. The better coupling of Gα 16 in fusion proteins was demonstrated by a decrease in the value of equilibrium dissociation constant (K A ) of agonists at Gα 16_ fused receptors (Table 2 vs. our previous data Randakova et al. [19], Table 3), except iperoxo at M 5 _Gα 16 (Table 3, discussed below). The elimination of interaction with other competitive Gα subunits as well as fusion alone could lead to better coupling of fused Gα subunits. The operational equilibrium dissociation constant K A quantifies an affinity of agonist to the conformation that initiates a given signaling pathway. Thus, it can be considered as one of the coupling parameters.
Furthermore, the operational efficacy τ to stimulate the non-canonical accumulation of IP X induced both by reference agonist carbachol and tested agonist oxotremorine is higher at fusion protein M 2 _Gα 16 than at co-transfected system M 2 + G α16 ( Table 3). The better coupling (both the decrease in K A and increase in τ) indicates that the fusion protein strategy is highly sensitive and thus suitable for detection and analysis of lowefficacy pathways.
Despite the high sensitivity, we did not detect accumulation of IPx induced by G i/o biased agonist JR6 and JR7 (Table 2) at any fused protein. These data further support the true G i/o bias of these novel agonists and also support the suitability of these fusion systems in the analysis of signaling bias.
Our data demonstrate that oxotremorine stimulates accumulation of IPx at M 2 _Gα 16 more efficiently in comparison with co-transfected system M 2 +Gα 16 , where the competition of endogenous G i/o and G q/11 occurs and more efficiently than at wt M 2 via endogenous G q/11 (Table 3). In our previous study of Randáková et al. [19], oxotremorine displayed lower RA i to stimulate the accumulation of IP X in the co-expressed system M 2 + Gα 16 than in the presented study. This discrepancy can be explained by different levels of expression of Gα 16 in co-expressed systems and points to the advantage of using fusion proteins with 1:1 stoichiometry for easier spotting of agonist bias. In comparison with our previous data [19], oxotremorine exerts bias towards IPx accumulation (via M 2 _Gα 16 ) over cAMP inhibition via G i/o at wt M 2 . Signaling bias of agonist oxotremorine to G 16 over G i/o and G q/11 pathway at M 2 receptor would be hard to reveal and quantify without fusion proteins due to signaling crosstalk or could be overlooked due to competition with other G-proteins. We show that using fusion proteins for this analysis can be very practical.
Furthermore, we demonstrate impairment of Gα 16 signaling of super-agonist iperoxo at the M 5 receptor. Besides worse coupling (lower pK A ) of iperoxo to fused M 5 _Gα 16 (Table 3), unlike other Gα 16 _fused receptors, super-agonist iperoxo stimulated accumulation of IP X at M 5 _Gα 16 with lower operational efficacy than reference agonist carbachol. In contrast at wt M 5 receptors expressed in CHO cells, iperoxo stimulated accumulation of IP X through cognate Gα q/11 with higher operational efficacy than carbachol. In CHO cells expressing wt M 5 co-transfected with Gα 16 , operational efficacy for carbachol and iperoxo was the same (Figure 4, Table 3), which could be explained by competition of Gα 16 with endogenous preferential G q/11 . Combined data thus indicate incompatibility of active M 5 receptor conformation specific to iperoxo with Gα 16 coupling and activation.

Construct Preparation
Constructs containing sequences of human variants of muscarinic acetylcholine receptors M 1 -M 5 fused with the human variant of Gα 15 subunit (also known as Gα 16 [73]) were prepared, and new stable cell lines of Chinese hamster ovary (CHO) expressing these fusion proteins were generated. Plasmids pcDNA3.1 coding human receptors M 1 -M 5 and Gα 16 subunit were obtained from Missouri S&T cDNA Resource Center (Rolla, MO, USA). Plasmid pCMV6-A-Hygro containing hygromycin as a mammalian selection marker was purchased from Origene (Rockville, MD, USA). The coding sequence for Gα 16 subunit and subsequently sequences for M 1 -M 5 receptors and were subcloned into the pCMV6-A-Hygro vector using restriction endonucleases. To this end, restriction site AflII at the N-terminus of the Gα 16 subunit and AgeI at the C-terminus of receptor sequences were created. Both parts were connected via short GATRARS linker, corresponding to the C-terminal amino acids in the M 2 sequence and N-terminal amino acid of the Gα 16 subunit. Cysteines needed for palmitoylation of receptors (C 435 at M 1 , C 457 at M 2 , C 561 at M 3 , C 470 at M 4 , and C 512 at M 5 ) were preserved. Sequences of all fusion proteins are in the Supplementary Material.

Homology Modeling
Homology models of fusion proteins were constructed as hybrid models using YASARA software, Biosciences (Vienna, Austria) [74]. For M 1 _Gα 16

Molecular Dynamics
The homology model of fusion proteins and structure of M 1 receptor in complex with G 11 G-protein (6OIJ) were aligned on the receptor part using MUSTANG [75]. The βγ-dimer from the 6OIJ structure was added to the homology model. To evaluate the stability of homology models, conventional molecular dynamics (MD) was simulated using Desmond/GPU ver. 6.1, D. E. Shaw Research (New York, NY, USA). The simulated system consisted of a receptor-G-protein complex in 1-palmitoyl-2-oleoyl-sn-glycero-3phosphocholine membrane set to receptor helices in water and 0.15 M NaCl. The system was first relaxed by the standard Desmond protocol for membrane proteins. Then 120 ns of NPγT (Nosé-Hoover chain thermostat at 300 K, Martyna-Tobias-Klein barostat at 1.01325 bar, isotropic coupling, Coulombic cut-off at 0.9 nm) molecular dynamics without restrains was simulated. The quality of molecular dynamics simulation was assessed by Simulation Quality Analysis tools of Maestro. CHO cells expressing individual Gα 16 _fused muscarinic receptors were grown to confluence in 75 cm 2 flasks in Dulbecco's modified EaglE s medium supplemented with 10% fetal bovine serum, Life Technologies (Carlsbad, CA, USA). One million cells were subcultured in 100 mm Petri dishes. Cells were washed with phosphate-buffered saline and harvested by mild trypsinization for functional experiments or manually for binding experiments on day five after subculture. After harvesting cells were centrifuged for 3 min at 250× g.

Cell Culture and Membrane Preparation
Membranes from CHO cells were prepared for binding experiments. The pellets of harvested cells were suspended in the ice-cold homogenization medium (100 mM NaCl, 20 mM Na-HEPES, 10 mM EDTA, pH = 7.4) and homogenized on ice by two 30 sec strokes using a Polytron homogenizer Ultra-Turrax; Janke & Kunkel GmbH & Co. KG, IKA-Labortechnik, (Staufen, Germany) with a 30-sec pause between strokes. Cell homogenates were centrifuged for 5 min at 1000× g. The supernatant was collected and centrifuged for 30 min at 30,000× g. Pellets were suspended in the washing medium (100 mM NaCl, 10 mM MgCl2, 20 mM Na-HEPES, pH = 7.4), left for 30 min at 4 • C, and then centrifuged again for 30 min at 30,000× g. The resulting membrane pellets were kept at −80 • C until assayed.

Radioligand Binding Experiments
All radioligand binding experiments were optimized and carried out as described by El-Fakahany and Jakubik [76]. Briefly, membranes (approximately 10 µg of membrane proteins per sample) were incubated in 96-well plates for 3 h at 25 • C in the incubation medium (100 mM NaCl, 20 mM Na-HEPES,10 mM MgCl 2 , pH = 7.4). In the case of the M 5 receptor, which has very slow kinetics of binding, the incubation time was extended to 5 h. Incubation volume for competition and saturation experiments with [ 3 H]NMS was 400 µL or 800 µL, respectively.
In saturation experiments of binding of [ 3 H]NMS six concentrations of the radioligand (ranging from 63 to 2000 pM) were used. Agonist binding was determined in competition experiments with 1 nM [ 3 H]NMS. Nonspecific binding was determined in the presence of 10 µM unlabeled atropine. Incubation was terminated by filtration through Whatman GF/C glass fiber filters, Whatman (Maidstone, GB using a Brandel harvester, Brandel (Gaithersburg, MD, USA). Filters were dried in a microwave oven (3 min, 800 W), and then solid scintillator Meltilex A was melted on filters (105 • C, 90 s) using a hot plate. The filters were cooled and counted in a Microbeta scintillation counter, PerkinElmer Waltham, MA, USA).

Measurement of Production of cAMP
Agonist-induced changes in the cAMP level were analyzed at Gα 16 _fused M 2 and M 4 receptors and M 2 and M 4 wild types. The level of cAMP was determined in radiochromatographical separation of [ 3 H]-cAMP as described previously [4]. To determine levels of cAMP, cells in suspension were pre-incubated for 1 h with 0.4 µM

Accumulation of Inositol Phosphates
The functional response of Gα 16 _fused muscarinic receptors was measured as an agonist-stimulated accumulation of inositol phosphates (IPX) using radiochemical chromatography as described previously [4]. The assay was performed in cells in suspension. IPX was determined after separation on ion-exchange columns Dowex 1 × 8-200, Sigma (St.Louis, MO, USA). Harvested cells were resuspended in Krebs-HEPES buffer (KHB; 138 mM NaCl; 4mM KCl; 1.3 mM CaCl 2 ; 1mM MgCl 2 ; 1.2 mM NaH 2 PO 4 ; 20 mM HEPES; 10 mM glucose; pH adjusted to 7.4) and centrifuged 250 g for 3 min. Cells were resuspended in KHB supplemented with 500 nM [ 3 H]myo-inositol, ARC (St.Lous, MO) and incubated at 37 • C for 1 h. Then they were washed once with an excess of KHB, resuspended in KHB containing 10 mM LiCl, and incubated for 1 h at 37 • C in the presence of indicated concentrations of agonists. The total reaction volume was 800 µL. Incubation was terminated by the addition of 0.5 mL of stopping solution (chloroform: methanol: 35% HCl; 2: 1: 0.1) and placed in 4 • C for 1 h. An aliquot (0.6 mL) of the upper (aqueous) phase was taken and loaded onto ion-exchange columns. Columns were washed with 10 mL of deionized water and 20 mL of 60 mM ammonium formate/5 mM sodium borate solution. IPX were collectively eluted from columns by 4 mL of 1 M ammonium formate-0.1 M/formic acid buffer. Level of IPx is expressed in dpm (decay per minute). Data are expressed as fold over basal level (after subtraction of blank value), the bottom (basal) is equal to 1.

Data and Analysis
Experiments were independent, using different seedings of CHO cells. Binding experiments were carried out in three experiments with samples in quadruplicates and functional assays were carried out at least in three experiments with samples in triplicate. Experimenters were blind to tested agonists.
After subtraction of non-specific binding (binding experiments) or background/blank values (functional experiments) data were normalized to control values determined in each experiment. IC 50 and EC 50 values and parameters derived from them (Ki and K A ) were treated as logarithms. All data were included in the analysis, no outliers were excluded. In statistical analysis value of p < 0.05 was taken as significant for all data. In multiple comparison tests ANOVA with p < 0.05 was followed by Tukey HSD post-test (p < 0.05). Data were processed in Microsoft office, analyzed, and plotted using the program Grace. The statistic was calculated using R (www.r-project.org, accessed on 13 September 2021).

[ 3 H]NMS Saturation Binding
The equilibrium dissociation constant (K D ) and maximum binding capacity (B MAX ) were determined in the saturation experiments. Non-specific binding in the presence of 10 µM atropine was subtracted to determine specific binding. Free concentration of [ 3 H]NMS was calculated by subtraction of values of specific binding from the final concentration of [ 3 H]NMS calculated from measurements of added radioactivity. Equation (1) was fitted to the data.
where y is specific binding at free concentration x. K D values are expressed as negative logarithms and B MAX values as pmol of binding sites per mg of membrane protein.

Competition Binding
The binding of tested agonists was determined in competition experiments with 1 nM [ 3 H]NMS fitting of Equation (2) for one-site competition or Equation (3) for twosite competition y = 100 − 100 * x x + IC 50 (2) where y is specific radioligand biding at concentration x of competitor expressed as a percent of binding in the absence of a competitor, IC 50 is concentration causing 50% inhibition of radioligand binding, flow is the fraction of low-affinity binding sites expressed in percents. Inhibition constants K I for analyzed agonists were calculated as where IC 50

Functional Response
The potency of analyzed agonists (EC 50 ) to induce maximal response (E MAX ) were obtained by fitting Equation (5) to the data from measurement of the accumulation of inositol phosphates, y =1 + E MAX − 1 * x nH EC nH 50 +x nH (5) where y is a functional response at a concentration of tested compound x, E MAX is the apparent maximal response to the tested compound, EC 50 is concentration causing halfmaximal effect and nH is slope factor (Hill coefficient). EC 50 values are expressed as negative logarithms and E MAX values as folds over basal.

Operational Model of Functional Agonism
The operational efficacy coefficient τ [77] was determined by fitting Equation (6) to data from the functional assay. y = E MAX * τ * x K A + (τ + 1) * x (6) where y is a functional response at a concentration of tested compound x, E MAX is the maximal response of the system, K A is the equilibrium dissociation constant. Equation (6) was fitted to data from functional experiments. Equation (6) was fitted to data by the two-step procedure described earlier [31]. In the first step, system E MAX was determined using carbachol, oxotremorine, and pilocarpine as internal standards by global fit to all data for a given receptor subtype and signaling pathway. In the second step, Equation (6) with E MAX fixed to the value determined in the first step was fitted to individual experimental data sets.

Relative Intrinsic Activity
For comparison of effects of agonists at different receptors fused with alpha Gα 16 to IPX signaling pathways, relative intrinsic activity (R Ai ) was calculated according to Griffin et al. [70].
where τa and K Aa are half-effective concentration and apparent maximal response to the tested compound, respectively. As Hill coefficients were equal to one, R Ai values were calculated according to Equation (8).
where EC 50a and E MAXa are half-effective concentration and apparent maximal response to the tested compound, respectively.

Signaling Bias
For receptors activating two or more signaling pathways, a ligand that has greater R Ai value for one pathway than for other(s) is biased to that pathway. Analogically, for a single signaling pathway and two or more receptors, a ligand that has greater R A i at one receptor than at other(s) is biased to a given pathway at that receptor.

Conclusions
The analysis of agonist bias at individual G-protein mediated pathways including non-preferential ones, plays a relevant role in the agonist screening and the development of drugs with reduced side effects that temper their clinical use. Our data showed that fusion proteins of muscarinic receptors and Gα subunits can serve as a suitable approach to analyze agonist bias and to serve as a convenient screening tool. Fusion proteins provide 1:1 receptor Gα stoichiometry, which makes quantification of agonist bias easier. We demonstrate that fusion of muscarinic receptors with Gα 16 limits access of other competitive Gα subunits to the receptor. That, in turn, makes it easier to quantify signaling via the non-canonical Gα 16 . We demonstrated agonist-specific activation of G 16 mediated pathway among individual subtypes of muscarinic receptors. We have confirmed functional selectivity of novel muscarinic agonists JR6 and JR7 for G i/o signaling pathway [19]. Furthermore, our data revealed signaling bias of oxotremorine towards non-canonical G 16 at M 2 and impairment of iperoxo mediated signaling through G 16 , regarding G i/o and G q/11 for M 2 and G q/11 for M 5 G-proteins.