Distinct Evolutionary Trajectories of Neuronal and Hair Cell Nicotinic Acetylcholine Receptors

Abstract The expansion and pruning of ion channel families has played a crucial role in the evolution of nervous systems. Nicotinic acetylcholine receptors (nAChRs) are ligand-gated ion channels with distinct roles in synaptic transmission at the neuromuscular junction, the central and peripheral nervous system, and the inner ear. Remarkably, the complement of nAChR subunits has been highly conserved along vertebrate phylogeny. To ask whether the different subtypes of receptors underwent different evolutionary trajectories, we performed a comprehensive analysis of vertebrate nAChRs coding sequences, mouse single-cell expression patterns, and comparative functional properties of receptors from three representative tetrapod species. We found significant differences between hair cell and neuronal receptors that were most likely shaped by the differences in coexpression patterns and coassembly rules of component subunits. Thus, neuronal nAChRs showed high degree of coding sequence conservation, coupled to greater coexpression variance and conservation of functional properties across tetrapod clades. In contrast, hair cell α9α10 nAChRs exhibited greater sequence divergence, narrow coexpression pattern, and great variability of functional properties across species. These results point to differential substrates for random change within the family of gene paralogs that relate to the segregated roles of nAChRs in synaptic transmission.

Functional nAChRs are homomeric, comprising five identical subunits, or heteromeric, formed by at least two different subunits (Karlin 2002;Millar and Gotti 2009;Zoli et al. 2015). The rules that govern the combinatorial assembly of functional nAChRs are complex and for the most part unknown. Some nAChR subunits can combine with other numerous, albeit specific, subunits. In contrast, some subunits can only form functional receptors with a limited subset (Karlin 2002;Millar and Gotti 2009). This segregation has given rise to subgroups of vertebrate nAChRs named for their initially described tissue of origin and main functional location. Thus, neuronal nAChRs are formed by multiple combinations of a2-a7 (and a8 in nonmammals) and b2-b4 subunits, comprising a wide combinatorial range with alternative stoichiometries (Millar and Gotti 2009;Zoli et al. 2015Zoli et al. , 2018. Muscle receptors show tighter coassembly rules. They have a typical a1 2 b1cd (or a1 2 b1ed) stoichiometry and do not coassemble with nonmuscle subunits (Mishina et al. 1986;Millar and Gotti 2009). Finally, the hair cell nAChR has a very strict coassembly constraint, being formed exclusively by a9 and a10 subunits (Elgoyhen et al. 2001;Sgard et al. 2002;Scheffer et al. 2007). Although a9 subunits can form functional homomeric receptors (Elgoyhen et al. 1994;Lipovsek et al. 2014), these are unlikely to play a major role in inner ear hair cells in vivo (Vetter et al. 2007). Although the a9 and a10 subunits were initially classified as a neuronal, the a9a10 receptor does not appear to be functionally present in the brain (Morley et al. 2018). A consequence of the differences in coassembly rules between the three subgroups of nAChRs is that, although muscle cells can toggle between at least two receptor variants and neurons are capable of expressing a great diversity of nAChRs, hair cells express only one type.
The complement of nAChR subunits is highly conserved across vertebrates and even more so in tetrapods (Dent 2006;Pedersen et al. 2019). This suggests a high gene family-wide negative selection pressure for the loss of paralogs and highlights the functional relevance of each individual subunit across the vertebrate clade. However, given the major differences in coexpression patterns and coassembly rules that distinguish the subgroups of nAChRs, in particular the contrast between neuronal and hair cell receptors, distinct evolutionary trajectories most likely took place in different members of the family. On the one hand, if only one (functionally isolated) receptor is expressed by a given cell type, then selection pressure could have acted on stochastic changes of the coding sequence of component subunits that altered receptor function. Such changes would not have resulted in deleterious effects on other nAChRs given the restricted expression pattern of the subunits. The aforesaid process could have dominated the evolutionary history of the hair cell a9a10 receptor. On the other hand, widely expressed subunits that coassemble into functional receptors with multiple others have been most probably subjected to strong negative selection pressure. Changes in the coding sequence that lead to changes in functional properties may have had deleterious effects on alternative receptor combinations expressed in different cell types. In this case, functional diversification could have arisen from stochastic changes in the expression patterns of receptor subunits, resulting in a given cell changing the subtype of receptor it expresses, while preserving individual subunit functionality. Such processes could have dominated the evolutionary history of neuronal nAChRs. These contrasting theoretical scenarios bring forward three predictions for the evolutionary history of nAChRs across the vertebrate clade. First, isolated (hair cell) subunits may show coding sequence divergence, whereas widespread (neuronal) subunits a high degree of coding sequence conservation. Second, isolated subunits may show low coexpression variation, whereas widespread subunits a great variability in coexpression patterns. Finally, isolated receptors may present divergent functional properties across species resulting from changes in coding sequences, whereas widespread receptors may show highly conserved functional properties when formed by the same subunits.
To test these predictions, we studied the molecular evolution and the variability in expression pattern, coupled to coassembly potential, of vertebrate nAChR subunits. Furthermore, we performed a comprehensive comparative analysis of the functional properties of the hair cell (a9a10) and the two main neuronal (a7 and a4b2) nAChRs from three representative species of tetrapods. We present strong evidence supporting the notion that, within the family of paralog genes coding for nAChR subunits, receptors belonging to different subgroups (hair cell or neuronal) underwent different evolutionary trajectories along the tetrapod lineage that were potentially shaped by their different coexpression patterns and coassembly rules. We propose that it is the difference in the most probable substrate for random change and subsequent selection pressure that separate the contrasting evolutionary histories of nAChRs subgroups.

Amino Acid Sequence Divergence and Coexpression Patterns Differentiate Neuronal from Hair Cell nAChR Subunits
The comparative analysis of gene paralogs provides the opportunity to test predictions about the evolutionary history of a gene family. To study the degree of conservation of coding sequences in subgroups of nAChRs, we performed an exhaustive evaluation of sequence divergence of vertebrate nAChR subunits. The analysis included amino acid sequences from 11 species of birds, reptiles, and amphibians, all groups that were notably underrepresented in previous work (Ortells and Lunt 1995;Tsunoyama and Gojobori 1998;Le Novère et al. 2002;Dent 2006;Franchini and Elgoyhen 2006;Lipovsek et al. 2012) and are of particular importance to help resolve differences between clades. Overall, we analyzed 392 sequences from 17 nAChR subunits belonging to 29 vertebrate species (supplementary table S1, Supplementary Material online, and fig. 1). Based on sequence identity, the family of nAChR subunits can be split into four groups: a, non-a, a7-like, and a9-like (figs. 1 and 2, supplementary fig. S1, Supplementary Material online, and Karlin [2002] and Millar and Gotti [2009]). a subunits are characterized by the presence of two consecutive cysteine residues in the ACh binding site. In the pentameric assembly, a subunits provide the principal component of the ACh binding domain and they require the coassembly with non-a subunits to form functional receptors. Non-a subunits lack the consecutive cysteines and contribute the complementary component of the ACh binding domain. a7-like subunits the possess consecutive cysteines and are capable of assembling homopentameric receptors, thus providing both the principal and complementary component of the ACh binding domain (figs. 1 and 2 and Karlin [2002]). Finally, a9-like subunits also have the consecutive cysteines and only form receptors with each other (Elgoyhen et al. 2001). As previously reported with smaller data sets (Franchini and Elgoyhen 2006;Lipovsek et al. 2012), the present extended analysis showed that a10 subunits are Nicotinic Receptors Evolution . doi:10.1093/molbev/msz290 MBE unique in presenting a segregated grouping of orthologs with nonmammalian a10 subunits as a sister group to all a9 subunits, and mammalian a10 subunits an outgroup to the a9/ nonmammalian a10 branch (figs. 1 and 2). This may relate to the overall low %seqID of all vertebrate a10 subunits, coupled to high sequence conservation within individual clades (supplementary table S2, Supplementary Material online), together with the higher rate of nonsynonymous substitutions reported for the mammalian clade (Franchini and Elgoyhen 2006;Lipovsek et al. 2012). Next, we evaluated whether an accumulation of clade-specific amino acid changes was present in any nicotinic subunit that could point to clade-specific changes in protein function. To achieve this, we searched for site-specific evolutionary shifts in amino acid biochemical state between clades using the DIVERGE 3.0 software (Gu et al. 2013). This analysis predicts amino acid sites that may be involved in between-clade functional divergence against the background of neutral evolution (Gu 2006). It first estimates a type II functional divergence coefficient (h II ) that indicates site-specific evolutionary shifts in amino acid biochemical state between clades (i.e., positively charged amino acid on a given site in clade 1 and negatively charged amino acid on the same site in clade 2). Subsequently, it computes the posterior probability that each individual site contributes to the clade-specific functional diversification. We found that both a9 and a10 subunits presented h II values significantly >0 (supplementary table S3 figure 2A, obtained with variation rates among sites modeled by a gamma distribution. Red branches, mammals; yellow branches, sauropsids; green branches, amphibians; blue branches, fish; light blue branches, coelacanth. Shadings denote the different groups of subunits: light greens, a subunits; dark green, non-a subunits; purple, a7-like subunits; orange, a9-like subunits. The tree was built using minimum evolution method and pairwise deletion for missing sites. The optimal tree with a sum of branch length of 47.32565515 is shown. For clarity, the percentage of replicate trees in which the associated taxa clustered together in the bootstrap test (1,000 replicates) are shown only next to the branches that separate different subunits. The tree is drawn to scale, with branch lengths in the same units as those of the evolutionary distances used to infer the tree. Marcovich et al. . doi:10.1093/molbev/msz290 MBE identified, in a9 subunits, the aspartate/asparagine (mammals/sauropsid) substitution in the extracellular domain and the alanine/aspartate (mammals/sauropsid) substitution at the exit of the channel pore ( fig. 2B, asterisks) that we previously reported to be involved in functional differences in calcium permeability between mammalian and avian receptors (Lipovsek et al. 2014). In contrast, neuronal subunits failed to show between-clade functional divergence at the sequence level (fig. 2B and supplementary table S3, Supplementary Material online). Overall, the extended molecular evolution analysis of nAChR subunits showed that while neuronal subunits were highly conserved across all vertebrate clades, a9 and a10 hair cell subunits showed a greater degree of between-clade sequence divergence that differentiates mammalian and sauropsid subunits.
The capability to coassemble into functional receptors and the coexpression of nAChR subunits within a given cell delineate the complement of receptors that shape its nicotinic ACh response. Numerous heterologous expression experiments and subunit-specific pharmacological studies have outlined a comprehensive repertoire of functionally validated pentameric assemblies (supplementary table S6, Supplementary Material online, and references therein). However, no systematic gene expression analysis that explores the potential spectrum of nAChRs in a given cell type has been performed. In order to evaluate coexpression patterns, we performed a meta-analysis of gene expression data from ten recent mouse single-cell transcriptomic studies (supplementary table S4, Supplementary Material online). Starting with the gene expression matrices and the cell types identified by each study, we used a Bayesian approach (Kharchenko et al. 2014) to estimate the likelihood of a gene being expressed at any given average level in a given cell type ( fig. 3A and supplementary table S5 and fig. S3, Supplementary Material online). We next combined these data with the catalog of validated nAChR pentamers (supplementary table S6, Supplementary Material online) and outlined the potential complement of pentameric receptors present in each cell type, by identifying the subunit combinations that are present within a 10-fold, 100-fold, or 1,000-  Taken together the analysis of amino acid sequences and coexpression patterns indicates that, although the a9 and a10 subunits have a restricted expression pattern, they also have the highest clade-specific sequence divergence. On the contrary, neuronal nAChRs present higher sequence conservation together with a widespread expression pattern.

Divergence of Biophysical Properties Differentiates Neuronal from Hair Cell nAChRs
The restricted coexpression pattern of hair cell a9 and a10 nAChR subunits ( fig. 3) and their exclusive coassembly with each other (Elgoyhen et al. 2001;Sgard et al. 2002;Scheffer et al. 2007), together with their higher level of between-clade sequence divergence ( fig. 2  Supplementary Material online) may have resulted in highly conserved functional properties of neuronal receptors assembled from the same subunits. We experimentally tested this prediction by performing a comprehensive side-by-side comparison of the functional properties of the two main neuronal, a4b2 and a7 nAChRs, and the hair cell a9a10 nAChR from three representative species of tetrapod clades. To this end, we injected Xenopus laevis oocytes with the corresponding cRNAs and characterized the biophysical properties of ACh responses.
The conserved desensitization profiles observed for both types of neuronal receptors was in stark contrast with that of a9a10 receptors. Although rat and chicken a9a10 receptors showed similar and somewhat slow desensitization, with 60-65% of current remaining 20 s after the peak response to ACh (P ¼ 0.9999), frog a9a10 nAChRs exhibited significantly higher desensitization when compared with rat a9a10 (P ¼ 0.0051) and chicken a9a10 (P ¼ 0.0042) receptors Neuronal nAChRs are potentiated by extracellular divalent cations such as Ca 2þ via an allosteric and voltageindependent mechanism (Mulle et al. 1992;Vernino et al. 1992). The rat a9a10 nAChR, on the other hand, is both potentiated and blocked by physiological concentrations of extracellular divalent cations (Weisstaub et al. 2002). Blockage occurs in the millimolar range, is voltage dependent, and proposed to occur as a result of calcium permeation (Weisstaub et al. 2002). To perform a comparative analysis of calcium modulation on rat, chicken, and frog a4b2, a7, and a9a10 receptors, responses to near-EC 50 concentrations of ACh were recorded in normal Ringer's solution at a range of Ca 2þ concentrations and normalized to the response at 1.8 mM Ca 2þ . For the neuronal a4b2 and a7 receptors from all three species, a similar potentiation pattern was observed, with increasingly higher responses to ACh at greater Ca 2þ concentrations ( fig. 6, top and middle panels). In contrast, responses of a9a10 nAChRs from rat, chicken, and frog  (10)  Representative responses of a4b2, a7, and a9a10 nAChRs to a 60 s (for a4b2 and a9a10) or 30 s (for a7) application of 100 lM ACh for all a4b2 and amniotes a9a10, and 1 mM ACh for all a7 and frog a9a10 nAChRs. Bottom panels. Percentage of current remaining 20 s (a9a10 and a4b2) or 5 s (a7) after the peak response, relative to the maximum current amplitude elicited by ACh. Bars represent mean 6 S.E.M., open circles represent individual oocytes (n ¼ 4-10). **P < 0.01, one-way ANOVA followed by Dunn's test (a4b2 nAChRs) or Kruskal-Wallis followed by Holm Sidak's test (a7 and a9a10 nAChRs).
FIG. 6. Extracellular Ca 2þ potentiates neuronal nAChRs but differentially modulates a9a10 nAChRs. ACh response amplitude as a function of extracellular Ca 2þ concentration. ACh was applied at near-EC50 concentrations (10 lM ACh for all a4b2, rat and chick a9a10 nAChRs and 100 lM ACh for all a7 and frog a9a10 nAChRs). Current amplitudes recorded at different Ca 2þ concentrations in each oocyte were normalized to the response obtained at 1.8 mM Ca 2þ in the same oocyte. V h : À90 mV. Bars represent mean 6 S.E.M., open circles represent individual oocytes (n ¼ 4-12). *P < 0.05, **P < 0.01, ***P < 0.005, and ****P < 0.0001, paired t-test (rat and frog a4b2 nAChRs and all a7 and a9a10 nAChRs) or Wilcoxon matched pair test (chick a4b2 nAChR)-comparing 0.5 mM Ca 2þ versus 3 mM Ca 2þ . Calcium permeation through nAChRs holds great functional significance for the activation of calcium-dependent conductances and intracellular signaling pathways (Dani and Bertrand 2007). Receptors containing the a7 subunit have high calcium permeability (S egu ela et al. 1993), whereas receptors containing a4b2 subunits have a lower contribution of calcium to the total current (Haghighi and Cooper 2000;Lipovsek et al. 2012). Amniote inner ear hair cell a9a10 nAChRs show differences in the extent of calcium permeability (Lipovsek et al. 2012(Lipovsek et al. , 2014. In order to perform a comparative analysis of the extent of the calcium component of ACh responses, we studied the differential activation of the oocyte's endogenous calcium-activated chloride current (ICl Ca ). In oocytes expressing a recombinant receptor with high calcium permeability, the ICl Ca is strongly activated upon ACh application (Barish 1983). Incubation of oocytes with the membrane-permeant fast Ca 2þ chelator BAPTA-AM subsequently abolishes the Cl À component of the total measured current (Gerzanich et al. 1994). Figure 7 shows representative responses to ACh before and after a 3-h incubation with BAPTA-AM for a4b2, a7, and a9a10 nAChRs from rat, chicken, and frog. Although ACh-evoked currents were only slightly affected in all a4b2 nAChRs denoting no major calcium influx (70-80% of current remaining after BAPTA incubation, fig. 7, left panel), all a7 receptors showed a strong reduction of the ACh response after BAPTA incubation (30-40% of current remaining after BAPTA), indicating significant calcium permeation ( fig. 7, middle panel). Thus, no interspecies differences in the proportion of calcium current for both the low calcium permeant a4b2 and the highly calcium permeant a7 neuronal receptors were observed (table 1). Conversely, and as previously reported (Lipovsek et al. 2012(Lipovsek et al. , 2014, we observed a marked difference in the extent of calcium current between the rat and chicken a9a10 receptors (P ¼ 0.0299- fig. 7, right panels and table 1). Moreover, the percentage of remaining response after BAPTA-AM incubation for the frog a9a10 receptor (fig. 7, right panels) was similar to that of the rat receptor (P ¼ 0.9999-table 1) and significantly different from that of the chicken a9a10 nAChR (P ¼ 0.0013-table 1).
Neuronal nAChRs are characterized by a strong inward rectification, with negligible current at depolarized potentials (S egu ela et al. 1993;Haghighi and Cooper 2000). This is proposed to be a relevant feature for their roles as presynaptic modulators of neuronal transmission (Haghighi and Cooper 2000). On the other hand, amniote a9a10 nAChRs exhibit a peculiar current-voltage (I-V) relationship due to a considerable outward current at positive potentials (Elgoyhen et al. 2001;Lipovsek et al. 2012). In order to perform a comparative analysis of the rectification profiles of neuronal and hair cell nAChRs, we obtained I-V curves and determined the ratio of current elicited at þ40 mV to that at À90 mV (I þ40 /I À90 ). All neuronal nAChRs exhibited I-V curves with marked inward rectification with no significant interspecies differences for either a4b2 or a7 tetrapod receptors, presenting I þ40 /I À90 values below 1 ( fig. 8, left and middle panels and table 1). On the contrary, each of the hair cell a9a10 nAChRs analyzed presented a unique I-V profile. As previously reported, rat a9a10 receptors showed significant outward current at depolarized potentials and greater inward current at hyperpolarized potentials with a I þ40 /I À90 ratio close to 1 ( fig. 8, right panels and Elgoyhen et al. [2001]). The chicken a9a10 receptor showed outward current similar to that of the rat a9a10 FIG. 7. Unlike neuronal nAChRs, a9a10 nAChRs exhibit differential Ca 2þ contribution to the total inward current. Top panels. Representative responses to near-EC 50 concentration of ACh (10 lM ACh for all a4b2 and amniotes a9a10 nAChRs and 100 lM ACh for all a7 and frog a9a10 nAChRs) in oocytes expressing a4b2, a7, and a9a10 nAChRs, before (gray traces) and after (color traces) a 3-h incubation with BAPTA-AM (V h ¼ À70 mV). Bottom panels. Percentage of the initial response remaining after BAPTA incubation. Bars represent mean 6 S.E.M., open circles represent individual oocytes (n ¼ 4-10). ****P < 0.0001, one-way ANOVA followed by Dunn's test (a4b2 and a7 nAChRs) or Kruskal-Wallis followed by Holm Sidak's test (a9a10 nAChRs).
receptor, but the inward current was smaller ( fig. 8, top right panel), resulting in a significantly higher I þ40 /I À90 ratio (P ¼ 0. 0229-table 1). Surprisingly, the frog a9a10 receptor showed an I-V profile similar to that of neuronal nAChRs, with strong inward rectification, almost no outward current at depolarized potentials ( fig. 8, right panels) and a I þ40 /I À90 below 1, significantly different to that obtained for chick (P ¼ 0.0406-table 1) and rat (P < 0.0001-supplementary table S9, Supplementary Material online) a9a10 receptors.

Comparative Functional Analysis of Neuronal and Hair Cell nAChRs Shows Distinct Evolutionary Trajectories
Altogether, the characterization of individual functional properties of tetrapod nAChRs showed a stark contrast between neuronal and hair cell nAChRs. In order to concomitantly analyze the diversification, or conservation of receptor function, we performed principal component analysis (PCA) on all the functional variables measured on a4b2, a7, and a9a10 receptors from the three species (supplementary tables S8 and S9, Supplementary Material online). The first two principal components (PC) accounted for 82% of the variability ( fig. 9). Moreover, the distribution of receptors in PCA space reflected their overall functional differences and similarities. Both neuronal a4b2 and a7 receptors occupied distinct regions, more distant in PC1 than in PC2 denoting that these receptors differ more on ACh apparent affinity, desensitization, and calcium permeability than they do on rectification and calcium modulation ( fig. 9, inset). Also, a4b2 and a7 receptors from the different species were located very close together, reflecting the lack of interspecies differences in functional properties. In contrast, the hair cell a9a10 receptors from the different species were distant from each other in PCA space, denoting their extensive functional divergence ( fig. 9). Interestingly, the frog a9a10 nAChR was closer to the a7 receptors than to its amniote counterparts, highlighting the overall functional similarity between the amphibian a9a10 and a7 nAChRs.
Amino acid sequence phylogenies, coexpression/coassembly patterns, and functional experiments support the hypothesis of differential evolutionary trajectories for neuronal and hair cell receptors. Further insight could be attained by the evaluation of the receptors present in the last common ancestor of rat and chicken (amniote ancestor) and of rat, chicken, and frog (tetrapod ancestor). To tackle this, we inferred the character state for the functional properties of the a4b2, a7, and a9a10 ancestral receptors (see Materials and Methods). We then projected this predicted functional states onto the functional PCA space of the extant receptors ( fig. 9). As expected, the predicted ancestral a4b2 and a7 amniote and tetrapod receptors were located close to their extant counterparts. In contrast, the predicted ancestral states of the a9a10 receptors were localized halfway between the frog a9a10 and the amniote chicken and rat a9a10 extant receptors, reflecting the greater evolutionary distance covered in functional space by the hair cell receptors. In summary, the data portrayed in figure 9 describes a plausible scenario for the functional evolutionary trajectories undertaken by two neuronal and the hair cell nAChRs within the tetrapod lineage.

Discussion
The expansion of ion channel families on the vertebrate stem branch was mostly driven by two rounds of whole genome duplications (Roux et al. 2017). This has played a crucial role in the evolution of nervous systems and provided raw FIG. 8. Hair cell, but not neuronal, nAChRs show differential current-voltage relationships across species. Top panels. Representative I-V curves obtained by the application of voltage ramps (À120 to þ50mV, 2 s) at the plateau response to 3 lM ACh for a4b2 and a9a10 or by the application of 100 lM ACh at different holding potentials for a7 nAChRs. Values were normalized to the maximal agonist response obtained for each receptor. Bottom panels. Ratio of current amplitude at þ40 mV relative to À90 mV for each oocyte. Bars represent mean 6 S.E.M., open circles represent individual oocytes (n ¼ 5-11). *P < 0.05 and ****P < 0.0001, one-way ANOVA followed by Dunn's test (a4b2 and a7 nAChRs) or Kruskal-Wallis followed by Holm Sidak's test (a9a10 nAChRs).
Nicotinic Receptors Evolution . doi:10.1093/molbev/msz290 MBE material that enabled the diversification of cell types (Liebeskind et al. 2015) resulting in the complexity reached by vertebrate brains (e.g., Sugino et al. 2019). Among the different ion channel families, Cys-loop receptors are within those that underwent the greatest expansion (Liebeskind et al. 2015). The entire extant complement of nAChRs, which includes 17 different subunits (a1-a10, b1-b4, d, e, and c), was already present in the last common ancestor of all vertebrates (Pedersen et al. 2019). With the exception of some fish species that acquired a11, b1.2, and b5 subunits, for which no expression or functional data have yet been reported (Pedersen et al. 2019), and the loss of the a7-like a8 subunits in mammals, the complement of original vertebrate nAChR subunits has been remarkably conserved. Moreover, nAChRs are unique in that subgroups of the family have distinct roles in synaptic transmission, either in the nervous system, inner ear hair cells, or the neuromuscular junction (Mishina et al. 1986;Dani and Bertrand 2007;Elgoyhen and Katz 2012). In this work, we performed in-depth analyses of coding sequence molecular evolution, subunit coexpression patterns, and comparative functional properties of neuronal and hair cell receptors to explore the potential impact of these segregation of nAChRs subgroups on their respective evolutionary histories. We found that neuronal subunits showed high degree of coding sequence conservation, coupled to greater coexpression variability and conservation of functional properties across tetrapod clades. In contrast, the hair cell a9a10 nAChR showed greater sequence divergence, a highly restricted coexpression pattern and a great degree of variability of functional properties across species. These results indicate that along the tetrapod lineage, neuronal and hair cell nAChRs underwent alternative evolutionary trajectories.

Functional Conservation of Neuronal nAChRs
The observation that the biophysical functional properties of neuronal a4b2 and a7 nAChRs were conserved in the three tetrapod species analyzed relates to their high degree of amino acid sequence conservation (figs. 1 and 2 and supplementary tables S2 and S3, Supplementary Material online). Cholinergic innervation is pervasive, with almost every area of the brain being influenced by nicotinic signaling (Dani and Bertrand 2007). Moreover, the expression of neuronal nAChR subunits is widespread in the central and peripheral nervous system ( fig. 3), where they assemble in multiple combinations (supplementary table S6, Supplementary Material online). Thus, randomly acquired changes in the coding sequence of a given subunit might be deleterious for receptor function in a multitude of heteropentameric assemblies present in diverse neuronal types. Such changes would therefore be under strong negative selection pressure. This is in agreement with the low degree of divergence observed for neuronal subunits and the absence of signatures of positive selection in the protein-coding sequences of a4, b2, and a7 subunits (Lipovsek et al. 2012) and other brain expressed genes (Haygood et al. 2010). However, neurons could alternatively resort to reshuffling of the complement of nAChR subunits being expressed to attain functional divergence. Therefore, and as for most brain expressed genes (Hoekstra and  (table 1). Square symbols represent rat nAChRs, circles represent chick nAChRs, and triangles represent frog nAChRs, a4b2 nAChRs are shown in shades of green, a7 nAChRs in shades of purple, and a9a10 nAChRs in shades of orange. The projected locations of inferred functional states are shown for amniote (stars) and tetrapod (crosses) ancestral receptors and colored in yellow (a4b2), blue (a7), or pink (a9a10). Inset. Biplot of the relative contribution of the five biophysical properties to PC1 and PC2.  Haygood et al. 2010), random changes in noncoding regions that lead to differential expression patterns across brain areas or species may have played a substantial role in delineating the evolutionary trajectories of neuronal nAChRs. Hitherto undescribed evolutionary processes of a higher complexity level that involve changes in the expression pattern, or function, of chaperon proteins that influence post-translational assembly and surface expression of neuronal nAChR subunits (Koperniak et al. 2013;Gu et al. 2016Gu et al. , 2019Matta et al. 2017;Dawe et al. 2019) may provide an additional substrate for functional divergence.

Differential Coexpression Patterns of Neuronal and Hair Cell nAChRs
Our meta-analysis of expression patterns across the mouse brain highlights numerous instances of potential functional variability and diversification, even between closely related neuronal types ( fig. 3, fig. 3), indicating that they may express differential levels of functionally distinct a4b2 containing (a4b2*) receptors. Inclusion of the a5 subunit can alter a4b2* receptor properties substantially, increasing ACh sensitivity, desensitization kinetics, and Ca 2þ permeability (Ramirez-Latorre et al. 1996;Tapia et al. 2007;Kuryatov et al. 2008). In addition, incorporation of the b3 subunit to the a4b2* receptor also increases ACh sensitivity, without significantly affecting Ca 2þ permeability (Tapia et al. 2007;Kuryatov et al. 2008). Moreover, VTA dopaminergic neurons also showed expression of the a6 and a3 subunits, both of which can coassemble with a4, a5 and/or b2 subunits, greatly enhancing the complexity of individual nAChRs-mediated responses of VTA dopaminergic neurons to modulatory cholinergic input. Another interesting example is provided by layer VI cortical pyramidal neurons, whose activity is modulated by a dense cholinergic innervation from the basal forebrain. Here, ACh elicits robust excitatory responses acting on a4b2* nAChRs, with layer VI being one of the few cortical areas that express the accessory a5 nAChR subunit (Proulx et al. 2014). Cortical neurons that project to both the ventral posteromedial nucleus (VPM) and the posteromedial complex of the thalamus express significantly higher levels of the a5 subunit than neurons projecting to the VPM only (supplementary fig. S3A, Supplementary Material online). This suggests VPM-only projecting neurons could have a lower density of a4a5b2 compared with a4b2 nAChRs, potentially contributing to differential cholinergic modulation of these subtypes of layer VI neurons, that also show differences in excitability (Landisman and Connors 2007).
Hair cells of the inner ear express high levels of a9 and a10 subunits, along with a number of neuronal (b2 and b4) and muscle (a1, b1, and c/e) subunits ( fig. 3 and Roux et al. [2016] and Scheffer et al. [2007]). a9a10 nAChRs mediate fast efferent neurotransmission to cochlear and vestibular hair cells in the inner ear  and are characterized by unique pharmacological and biophysical properties (Elgoyhen et al. 1994(Elgoyhen et al. , 2001Katz et al. 2000;Weisstaub et al. 2002). Most notably, nicotine, the diagnostic agonist of the nAChRs family, does not act as an agonist of a9a10 receptors, but as a competitive antagonist (Elgoyhen et al. 2001). In inner ear hair cells, no responses to nicotine application are detected Scheffer et al. 2007), indicating that functional muscle or neuronal nAChRs are not present at the plasma membrane. The presence of neuronal and muscle subunits mRNA may result from redundant or residual transcriptional regulation mechanisms. Moreover, similar "leaky" expression of muscle subunits was detected in a number of neuronal types ( fig. 3).
Expression of a9 and a10 subunits is restricted to inner ear hair cells, with a few interesting exceptions. In the inner ear, spiral ganglion neurons (SGNs) provide afferent innervation to cochlear hair cells and express a range of neuronal nAChR subunits ( fig. 3 and Hiel et al. [1996] and Morley et al. [1998]). Interestingly, low levels of a9 and a10 subunits are present in SGNs ( fig. 3 and Shrestha et al. [2018]) with similarly low coexpression detected by two independent single-cell RNA sequencing studies (Petitpre et al. 2018;Sun et al. 2018). If this low level of a9 and a10 mRNA proves to be more than "transcriptional noise," then SGNs may be unique among neurons in expressing the hair cell a9a10 receptor. This might be related to the shared developmental origin of SGNs and hair cells at the otic placode (Fritzsch and Beisel 2004;Arendt et al. 2016) and could open the possibility that in addition to neuronal nAChRs, which are thought to partly mediate the nicotinic effect of lateral olivocochlear terminals on afferent dendrites (Reijntjes and Pyott 2016), a9a10 nAChRs may also play a role. Finally, in dorsal root ganglia neurons, a9 expression was not detected, whereas a10 is present at very low levels in only a few subtypes ( fig. 3 and Usoskin et al. [2015]). These observations support those reported by quantitative PCR (qPCR) and functional assays (Hone et al. 2012) and provide further evidence that the participation of a9* nAChRs in pain processes is not due to its presence in dorsal root ganglia neurons (Rau et al. 2005;Hone et al. 2012).
Of note, a9 and a10, together with other nAChR subunits, are expressed in other peripheral, nonneuronal, tissues (McIntosh et al. 2009;Zoli et al. 2018). A plausible autocrine/paracrine effect of ACh in these cells can be served by a multiple and redundant battery of nAChRs that might play a signaling function in these peripheral tissues (Zakrzewicz et al. 2017;Criado 2018). Due to the redundancy in pathways for ACh signaling, it is unlikely that the function of a9a10 nAChRs in these peripheral tissues provided the selection forces that shaped the accumulation of nonsynonymous changes on this receptor.

The a9a10 nAChR and the Evolution of the Efferent Olivocochlear System
The observations that a9 and a10 genes are only coordinately transcribed in inner ear hair cells, together with their ability to only form functional heteromeric receptors when coassembled with each other but not with other nAChR subunits Nicotinic Receptors Evolution . doi:10.1093/molbev/msz290 MBE (Elgoyhen et al. 1994(Elgoyhen et al. , 2001Scheffer et al. 2007), support our hypothesis that evolutionary changes in the hair cell receptors may have been focused at the coding sequence. Accordingly, vertebrate a9 and a10 subunits exhibit significant sequence divergence (supplementary table S3, Supplementary Material online, and fig. 2), with mammalian a10 subunits showing a higher than expected accumulation of nonsynonymous substitutions that were positively selected (Franchini and Elgoyhen 2006;Lipovsek et al. 2012). In addition, both a9 and a10 subunits show a high number of clade-specific (mammalian vs. sauropsid) functionally relevant amino acid changes ( fig. 2, supplementary table S3, Supplementary Material online, and Lipovsek et al. [2014]). Consequently, the biophysical properties of a9a10 receptors drastically changed across tetrapod species ( fig. 9 and Lipovsek et al. [2012Lipovsek et al. [ , 2014). Since the primary function of the a9a10 receptor is at the postsynaptic side of the olivocochlear synapse, it can be hypothesized that clade-specific differences in efferent modulation of hair cell activity could have shaped the functional properties of a9a10 receptors. Upon the transition to land, the hearing organs of tetrapods underwent parallel evolutionary processes, mainly due to the independent emergence of the tympanic middle ear, at least five times, in separate groups of amniotes (Manley 2017). This was followed by the independent elongation of the auditory sensory epithelia that extended the hearing range to higher frequencies and the elaboration of passive and active sound amplification mechanisms that lead to the fine tuning of sound detection (Dallos 2008;Manley 2017). More importantly, mammals and sauropsids underwent a parallel diversification of hair cell types, segregating, partially in birds but completely in mammals, the phonoreception and sound amplification functions (Koppl 2011). Efferent innervation to hair cells, mediated by a9a10 nAChRs, is an ancestral feature common to all vertebrate species (Sienknecht et al. 2014). In the auditory epithelia, it modulates sound amplification and followed the hair cell diversification pattern: In birds, it mainly targets short hair cells, whereas in mammals, it targets outer hair cells (Koppl 2011). The latter developed a clade-specific sound amplification mechanism driven by the motor protein prestin and termed somatic electromotility (Dallos 2008). Prestin, together with bV giant spectrin, a major component of the outer hair cells' cortical cytoskeleton which is necessary for electromotility, shows signatures of positive selection in the mammalian clade that may relate to the acquisition of somatic electromotility (Franchini and Elgoyhen 2006;Cortese et al. 2017). Thus, the mammalian clade-specific evolutionary processes observed in both the a9 and a10 subunits (Franchini and Elgoyhen 2006;Lipovsek et al. 2012Lipovsek et al. , 2014 may be related to overall changes in the efferent olivocochlear systems of this clade that is tasked with the modulation of prestin-driven somatic electromotility. A recent high throughput evolutionary analysis identified 167 inner ear expressed genes with signatures of positive selection in the mammalian lineage (Pisciottano et al. 2019). These inner ear genes, including those encoding the a9 and a10 nAChR subunits, can be considered as hotspots for evolutionary innovation in the auditory system across species.
This scenario provides a context for evaluating the relationship between evolutionary trajectories and the functional role of a9a10 receptors. In mammals, the high calcium influx through a9a10 receptors activates large conductance, voltage, and low-calcium-sensitive BK potassium channels mediating hyperpolarization of outer hair cells in higher frequency regions of the cochlea (Wersinger et al. 2010). In contrast, in short hair cells from the chicken basilar papillae, hyperpolarization is served by the ACh-dependent activation of high calcium sensitive SK potassium channels (Fuchs and Murrow 1992;Samaranayake et al. 2004). Moreover, in contrast to adult mammalian hair cells where efferent fibers directly contact outer hair cells, but not the inner hair cells that release glutamate to activate afferent auditory fibers, efferent innervation in birds and amphibians coexists with afferent innervation in the same hair cells. Calcium influx in these hair cells could therefore result in efferent-triggered, ACh-mediated release of glutamate to auditory afferents due to calcium spill over, bypassing sound mechanotransduction. Thus, limiting the extent of calcium influx through a9a10 nAChRs may be paramount to avoid confounding sensory inputs. In this hypothetical scenario, the low calcium permeability of the avian a9a10 nAChR or the very high desensitization kinetics of the amphibian a9a10 nAChR that restrict calcium load could be related to the aforementioned selection pressure.

Subgroups of nAChRs and Differential Sources of Functional Divergence
Our observations on expression pattern, coding sequence, and functional divergence provide further evidence in support of the proposition that a9 and a10 are not a subtype of brain nicotinic subunit (for review see Morley et al. [2018]), but form a group of their own, characterized by unique expression profile, pharmacological, and biophysical properties (Elgoyhen et al. 1994(Elgoyhen et al. , 2001Katz et al. 2000;Weisstaub et al. 2002) and evolutionary history.
The contrasting evolutionary trajectories of neuronal and hair cell receptors, with functional variability stemming from combinatorial coexpression for the former and changes in coding sequence for the latter, support the notion of differential substrates for random change and ensuing functional divergence. For neuronal subunits, the source of random variability may have been rooted on changes in regulatory sequences. In contrast, for the hair cell receptor, random changes in the coding sequence were fixed throughout the evolutionary history of the tetrapod lineage. Interestingly, muscle subunits showed relatively low levels of coding sequence conservation ( fig. 2) and, via combinatorial coassembly, muscle cells can toggle between at least two receptor variants (Mishina et al. 1986). This places muscle receptors in between the two extremes of hair cell (isolated) and neuronal (widespread) receptors. A comparative functional study of muscle receptors would further test our hypothesis, with the prediction that a modest level of functional divergence may be encountered, but outweighed by the functional differences between muscle receptor variants. Marcovich et al. . doi:10.1093/molbev/msz290 MBE In summary, the present work provides evidence in support of different evolutionary trajectories for neuronal and hair cell nAChRs. These may have resulted from the differential substrates for random change that dominated evolutionary processes in each receptor subgroup: diversity of coexpression/coassembly patterns for neuronal subunits, changes in coding sequence for hair cell subunits. Finally, the simultaneous analysis of coding sequences, expression patterns, and protein functional properties generated new insights into the evolutionary history of gene paralogs, thus providing further context for the role of nAChRs in neuronal and hair cell synaptic transmission.

Materials and Methods
All experimental protocols were carried out in accordance with the guides for the care and use of laboratory animals of the National Institutes of Health and the Institutional Animal Care and Use Committee of the Instituto de Investigaciones en Ingenier ıa Gen etica y Biolog ıa Molecular, "Dr. H ector N. Torres."

Phylogenetic Analysis of Vertebrate nAChRs Subunits
All sequences were downloaded from GenBank (www.ncbi. nlm.nih.gov/genbank; last accessed December 2018), UCSC (http://genome.ucsc.edu/; last accessed December 2018), and Ensembl (www.ensembl.org; last accessed December 2018) databases. The signal peptides of all sequences were excluded from the analysis since they are not present in the mature functional protein. Accession numbers are listed in supplementary table S1, Supplementary Material online. All sequences were visually inspected, and missing and/or incorrect exons were obtained from the NCBI Genome Project traces database (http://blast.ncbi.nlm.nih.gov/Blast.cgi; last accessed December 2018). Sequence alignment was performed using ClustalW on the MEGA7 software (www.megasoftware.net; last accessed December 2018; Kumar et al. 2016), with the following parameters: for pairwise alignments, gap opening penalty: 10, gap extension penalty: 0.1; for multiple alignments, gap opening penalty: 10, gap extension penalty: 0.2. Protein weight matrix: Gonnet. Residue specific and hydrophilic penalties: ON. Gap separation distance: 4. End gap separation: OFF and no negative matrix was used. The delay divergent cutoff was 30%. The full alignment for the nicotinic subunits from representative vertebrate species is available in supplementary file 1, Supplementary Material online, in fasta format.
The phylogenetic tree of all nAChR subunits was built using the MEGA7 software. Positions containing alignment gaps and missing data were eliminated only in pairwise sequence comparisons. The final data set contained a total of 773 positions. The evolutionary distances (i.e., number of amino acid substitutions per site) were computed using the JTT matrix-based method (Jones et al. 1992). The Neighbor-Joining algorithm (Saitou and Nei 1987) was used to generate the initial tree and branch support was obtained by bootstrap test (1,000 replicates) (Felsenstein 1985). The evolutionary history was inferred using the minimum evolution method (Rzhetsky and Nei 1992). A first tree was generated with variation rate among sites modeled by a gamma distribution ( fig. 1) and a second tree assuming uniform variation rates among sites (supplementary fig. S1, Supplementary Material online). Overall, tree topology was similar between both methods.
Average percentage sequence identity was calculated for each subunit using the percentage of sequence identity between each pair of sequences (supplementary file 2, Supplementary Material online) from the same category for all sequences and for within or between mammalian and/or sauropsid sequences. For a10 subunits, the average percentage of sequence identity was also calculated for the nonmammalian paraphyletic group. Values obtained are summarized in supplementary table S2, Supplementary Material online.

Functional Divergence Analysis
The DIVERGE 3.0 software was used to test for functional diversification of nAChR subunits between the mammalian and sauropsid clades. DIVERGE predicts amino acid sites that may be involved in between-clade functional divergence against the background of neutral evolution (Gu 2006). Briefly, it estimates the type II functional divergence coefficient (h II ) that indicates site-specific evolutionary shifts in amino acid biochemical state between clades and then uses a Bayesian approach to compute the posterior probability that each individual site contributes to the clade-specific functional diversification. Type II sites represent amino acids that are highly conserved within each clade, but in a different biochemically state between clades (i.e., positively charged in clade 1 and negatively charged in clade 2).
To implement the DIVERGE analysis, multiple alignments of protein sequences for each individual subunit were generated using the MEGA7 software as described above. The highly variable amino-terminal signal peptides and intracellular domains were excluded from the analysis. Phylogenetic trees, with a topology corresponding to the species tree, were constructed for each nicotinic subunit by maximum likelihood and the JTT matrix-based method (Jones et al. 1992). The rates among sites were modeled as a gamma distribution. All positions with <95% site coverage were eliminated. a8, b1, and e nAChR subunits were not included in the analysis due to lack of mammalian and/or sauropsid sequences to perform suitable comparisons. Subsequently, protein sequence alignments and their corresponding phylogenetic trees (supplementary file 3, Supplementary Material online) were used as input data for DIVERGE 3.0 type II functional divergence analysis (Gu et al. 2013), with default parameters. h II and site-specific posterior probabilities were calculated for each subunit. A h II value significantly >0 (P < 0.05) indicates that residues conserved within each group have undergone radical changes in amino acid identity between groups (Gu 2006). z-Scores were used to test the significant difference of h II coefficients against the null hypothesis (h II ¼ 0) that implies no sites are present in the protein that reflect between-clade functional divergence (supplementary table  S3

Analysis of nAChR Subunit Expression in Single-Cell RNAseq Data Sets
A meta-analysis of single-cell gene expression data from ten studies was performed to describe the expression patterns of nAChR subunits across cell types. Processed gene expression data tables were obtained from ten single-cell RNAseq studies that evaluated gene expression in retina (Shekhar et al. 2016), inner ear sensory epithelium (Burns et al. 2015;McInturff et al. 2018), spiral ganglion (Shrestha et al. 2018), ventral midbrain (La Manno et al. 2016, hippocampus (Cembrowski et al. 2018), cortex (Chev ee et al. 2018, hypothalamus (Romanov et al. 2017), visceral motor neurons (Furlan et al. 2016), and dorsal root ganglia (Usoskin et al. 2015). Accession numbers, cell types inferred, and number of cells analyzed are summarized in supplementary table S4, Supplementary Material online. For all data sets, we used the published gene expression quantification and cell type labels. Each data set was analyzed separately. For the retina data set, we used the Smart-Seq2 sequencing data from Vsx2-GFP positive cells (Shekhar et al. 2016). For the gene expression quantification, we only analyzed four cell types that had a minimum number of cells in the data set that allowed reliable fitting to the error models: BC1A, BC5A, BC6, and RBC. From the layer VI somatosensory cortex data set (Cembrowski et al. 2018), we used a subset of the expression matrix that corresponds to day 0 (i.e., control, undisturbed neurons) of their experimental manipulation. For the hypothalamic neurons data set (Romanov et al. 2017), we used a subset that contained only neurons from untreated (control) mice and only quantified gene expression on the ten broad cell types identified. From the ventral midbrain dopaminergic neurons data set ), we used a subset comprising DAT-Cre/tdTomato positive neurons from P28 mice. For the SGNs data set, we used a subset comprising type I neurons from wild type mice (Shrestha et al. 2018). For the utricle hair cell data sets, we used the normalized expression data of McInturff et al. (2018). For the cochlear hair cell data, we used the normalized expression data from Burns et al. (2015) and continued our analysis with only the ten cochlear hair cells identified. For the visceral motor neurons data set (Furlan et al. 2016), we excluded the neurons that were "unclassified" from further analysis. For the dorsal root ganglia data set (Usoskin et al. 2015), we used a subset containing only successfully classified neurons that were collected at room temperature. Inspection of all data sets for batch effects was performed using the scater package (version 1.10.1) (McCarthy et al. 2017). All data analysis was implemented in R (version 3.5.1) and Bioconductor (version 3.8) (http://www.bioconductor.org/; last accessed December 2018), running on RStudio (version 1.1.456) (http://www.rstudio.com/; last accessed December 2018).
The publicly available expression matrices for a number of the data sets contained raw counts (retina, hippocampus, hypothalamus, midbrain, and visceral motor neurons). For each of these data set separately, we performed a normalization step using the scran package (version 1.10.2) (Lun et al. 2016) that computes pool-based size factors that are subsequently deconvolved to obtain cell-based size factors.
The normalized expression matrices and cell type information were used as input to quantify cell-type-specific gene expression. Analysis was performed using the scde package (version 1.99.1) (Kharchenko et al. 2014). We modeled the gene expression measurements from each individual cell as a weighted mixture of negative binomial and low-magnitude Poisson distributions. The former accounts for the correlated component in which a transcript is detected and quantified, whereas the latter accounts for drop-out events in which a transcript fails to amplify when present. The weighting of the two distributions is determined by the level of gene expression in the cell population (Kharchenko et al. 2014). We then used these error models to estimate the likelihood (joint posterior probability) of a gene being expressed at any given average level in a given cell type (Kharchenko et al. 2014). Probability distributions for all nAChR subunit genes detected in all the cell types analyzed are shown in supplementary figure S2, Supplementary Material online. This whole transcriptome analysis provides accurate estimations of gene expression levels, thus allowing for the comparison of individual genes within a given cell type (i.e., the complement of nAChR subunits) or the analysis of expression level differences between cell types of the same data set (i.e., change in expression level of nAChR subunits between neuronal subtypes). Inferred mean expression values are summarized in supplementary table S5, Supplementary Material online, and normalized mean expression values are depicted in figure 3.
Subsequently, we combined the information about the complement of nAChR subunits for each cell type with a comprehensive catalog of experimentally validated subunit combinations (supplementary table S6, Supplementary Material online, and references therein). We identified the subunit combinations that were present, in each cell type, within a 10-fold, 10-to 100-fold, or 100-to 1,000-fold range of expression level or absent all together ( fig. 3B). Admittedly, this analysis approach overlooks the complexities of posttranslational modifications, receptor assembly, role of chaperone proteins, and transport to the plasma membrane. However, it provides conservative estimates of the maximum potential of combinatorial assembly of subunits and thus a maximum for the repertoire of nAChR assemblies that could be present at the cell membrane.

Expression of Recombinant Receptors in X. laevis Oocytes and Electrophysiological Recordings
Rat and chick a9 and a10 cDNAs subcloned into pSGEM, a modified pGEMHE vector suitable for X. laevis oocyte expression studies, were described previously (Elgoyhen et al. 1994(Elgoyhen et al. , 2001Lipovsek et al. 2012). Rat a4, b2, and a7 subunit cDNAs subcloned into pBS SK(À) (Agilent Technologies, Santa Clara, CA) were kindly provided by Dr Jim Boulter (University of California, Los Angeles, CA). Chicken a4 and b2 subunit cDNAs subcloned into pCI (Promega, Madison, WI) were Marcovich et al. . doi:10.1093/molbev/msz290 MBE kindly provided by Dr Isabel Bermudez-D ıaz (Oxford Brookes University, Oxford, UK). Chicken a7 subunit cDNA cloned into pMXT was kindly provided by Dr Jon Lindstrom (University of Pennsylvania) and was then subcloned into pSGEM between HindIII and SalI sites. Frog a4, b2, a7, a9, and a10 subunits were cloned from whole brain Xenopus tropicalis cDNA. Total RNA was prepared from whole brains using the RNAqueous -Micro kit AM1931 (Ambion, Thermo Fisher Scientific, Boston, MA). First strand cDNA synthesis was performed using an oligodT and the ProtoScript Taq RT-PCR kit (New England Biolabs, Ipswich, MA). Second strand synthesis was performed with the NEBNext mRNA Second Strand Synthesis Module kit -E6111S (New England Biolabs). Full length cDNAs for each subunit were PCR amplified (MultiGene 60 OptiMaxTM Thermal Cycler -Labnet International Inc. Edison, NJ) using specific primers (supplementary table S7, Supplementary Material online). PCR products were sequenced and subcloned into pSGEM between EcoRI and XhoI sites for a9, a7, and b2 nAChRs subunits, between HindIII and XhoI sites for the a10 subunit and between EcoRI and HindIII sites for the a4 subunit. All expression plasmids are readily available upon request.
Capped cRNAs were in vitro transcribed from linearized plasmid DNA templates using the RiboMAX Large Scale RNA Production System-T7 (Promega). The maintenance of X. laevis, as well as the preparation and cRNA injection of stage V and VI oocytes, has been described in detail elsewhere (Katz et al., 2000). Briefly, oocytes were injected with 50 nl of RNasefree water containing 0.01-1.0 ng of cRNAs (at a 1:1 molar ratio for a9a10 and a4b2 receptors) and maintained in Barth's solution (in mM): NaCl 88, Ca(NO 3 ) 2 0.33, CaCl 2 0.41, KCl 1, MgSO 4 0.82, NaHCO 3 2.4, HEPES 10, at 18 C.
Electrophysiological recordings were performed 2-6 days after cRNA injection under two-electrode voltage clamp with an Oocyte Clamp OC-725B or C amplifier (Warner Instruments Corp., Hamden, CT). Recordings were filtered at a corner frequency of 10 Hz using a 900BT Tunable Active Filter (Frequency Devices Inc., Ottawa, IL). Data acquisition was performed using a Patch Panel PP-50 LAB/1 interface (Warner Instruments Corp.) at a rate of ten points per second. Both voltage and current electrodes were filled with 3 M KCl and had resistances of $1 MX. Data were analyzed using Clampfit from the pClamp 6.1 software (Molecular Devices, Sunnyvale, CA). During electrophysiological recordings, oocytes were continuously superfused (10 ml min À1 ) with normal frog saline composed of (in mM): 115 NaCl, 2.5 KCl, 1.8 CaCl 2 , and 10 HEPES buffer, pH 7.2. In order to minimize the activation of the oocyte's native Ca 2þ -sensitive chloride current (ICl Ca ) by inward Ca 2þ current through the nAChRs, all experiments, unless otherwise stated, were carried out in oocytes incubated with the membrane-permeant Ca 2þ chelator 1,2-bis (2-aminophenoxy)ethane-N,N,N 0 ,N 0 -tetraacetic acid-acetoxymethyl ester (BAPTA-AM; 100 lM) for 3 h prior to electrophysiological recordings. This treatment was previously shown to effectively chelate intracellular Ca 2þ ions and, therefore, to impair the activation of the ICl Ca (Gerzanich et al. 1994). All recordings were performed at À70 mV holding potential, unless otherwise stated.

Biophysical Properties of nAChRs
Concentration-response curves were obtained by measuring responses to increasing concentrations of ACh. Current amplitudes were normalized to the maximal agonist response in each oocyte. The mean and S.E.M. values of the responses are represented. Agonist concentration-response curves were iteratively fitted, using Prism 6 software (GraphPad Software Inc.), with the equation: I/I max ¼ AnH/(AnH þ EC 50 nH), where I is the peak inward current evoked by the agonist at concentration AnH, I max is current evoked by the concentration of agonist eliciting a maximal response, EC 50 is the concentration of agonist inducing half-maximal current response, and nH is the Hill coefficient.
Desensitization of Ach-evoked currents was evaluated via prolonged agonist applications. The percentage of current remaining 5 s (for a7 nAChRs) or 20 s (for a4b2 and a9a10 nAChRs) after the peak of the response was determined for each oocyte.
The effects of extracellular Ca 2þ on the ionic currents through nAChRs were studied by measuring the amplitudes of the responses to a near-EC 50 concentration of ACh (10 lM for all a4b2 and amniote a9a10 nAChRs, and 100 lM for all a7 and frog a9a10 nAChRs) on extracellular Ca 2þ ranging from nominally 0 to 3 mM at a holding potential of À90 mV (Weisstaub et al. 2002). These experiments were carried out in oocytes injected with 7.5 ng of an oligonucleotide (5 0 -GCTTTAGTAATTCCCATCCTGCCATGTTTC-3 0 ) antisense to connexinC38 mRNA (Arellano et al. 1995;Ebihara 1996) to minimize the activation of the oocyte's nonselective inward current through a hemigap junction channel that results from the reduction of external divalent cation concentration. Current amplitudes at each Ca 2þ concentration were normalized to that obtained in the same oocyte at a 1.8 mM Ca 2þ .
I-V relationships were obtained by applying 2 s voltage ramps from À120 to þ50 mV from a holding potential of À70 mV, at the plateau response to 3 lM ACh. Leakage correction was performed by digital subtraction of the I-V curve obtained by the same voltage ramp protocol prior to the application of ACh. Generation of voltage protocols and data acquisition were performed using a Patch Panel PP-50 LAB/1 interface (Warner Instruments Corp.) at a rate of ten points per second and the pClamp 7.0 software (Axon Instruments Corp., Union City, CA). Current values were normalized to the maximum amplitude value obtained for each oocyte. The fast desensitizing a7 receptors had negligible plateau currents. For these receptors, responses to 100 lM ACh were obtained at different holding potentials and normalized to the amplitude response at À120 mV in the same oocyte. Table 1 summarizes the biophysical properties and statistical comparisons from rat, chicken, and frog a4b2, a7, and a9a10 receptors.

Statistical Analysis
Shapiro-Wilks normality test was conducted using custom routines written in R v3.4.1, through RStudio software v1.0.153. Statistical significance was determined using either parametric paired t-test or one-way ANOVA followed by Nicotinic Receptors Evolution . doi:10.1093/molbev/msz290 MBE Holm Sidak's test, or nonparametric Wilcoxon or Kruskal-Wallis tests followed by Dunn's tests conducted using Prism 6 software (GraphPad Software Inc.). A P < 0.05 was considered significant.

Reagents
All drugs were obtained from Sigma-Aldrich (Buenos Aires, Argentina). ACh chloride was dissolved in distilled water as 100 mM stocks and stored aliquoted at À20 C. BAPTA-AM was stored at À20 C as aliquots of a 100 mM stock solution in dimethylsulfoxide, thawed and diluted into Barth's solution shortly before incubation of the oocytes. ACh solutions in Ringer's saline were freshly prepared immediately before application.

Inference of Character State of Functional Properties of Ancestral Receptors
The pipeline followed to infer the ancestral character state of biophysical properties of nAChRs is schematized on supplementary figure S4, Supplementary Material online. Briefly, we first reconstructed the ancestral tetrapods and amniote DNA sequences of the a4, a7, a9, a10, and b2 nAChRs subunits (supplementary file 4, Supplementary Material online). For that purpose, we used the same ortholog sequences that were used to construct the phylogenetic tree of figure 1, together with a species tree with no branch lengths obtained from Ensembl (https://www.ensembl.org/info/about/speciestree.html; last accessed December 2018). Inferred ancestral DNA sequences for the amniote and tetrapod nodes (supplementary file 4, Supplementary Material online) were obtained, for all three codon positions, by the maximum likelihood method (Nei and Kumar 2000) under the Tamura-Nei model (Tamura and Nei 1993), on the MEGA7 software (Kumar et al. 2016). The initial tree corresponds to the provided Species Tree with very strong restriction to branch swap. The rates among sites were treated as a gamma distribution. All positions with <95% site coverage were eliminated.
Subsequently, multiple alignments including extant rat, chick, and frog and ancestral amniote and tetrapod amino acid sequences were performed using MEGA7 for each studied nAChR (supplementary file 5, Supplementary Material online). The sequence identity was used to assign branch length values to a7, a4b2, and a9a10 nAChRs trees, corresponding to 1-SeqID (supplementary table S10, Supplementary Material online). Theoretical concatemeric constructions were built for the heteromeric a9a10 nAChRs considering the descripted prevalent (a9) 2 (a10) 3 stoichiometry (Plazas et al. 2005). For the a4b2 nAChRs, the high sensitivity (a4) 2 (b2) 3 stoichiometry was used to generate the theoretical concatemeric receptor (supplementary file 5, Supplementary Material online). The resulting trees (supplementary fig. S5, Supplementary Material online) were used, together with the biophysical properties experimentally determined for the extant receptors (table 1) as input data for ancestral character inference. ACh sensitivity (EC 50 values), desensitization rates (% of remaining I 20 s after ACh peak), Ca 2þ modulation (I elicited by ACh at Ca 0.5 mM/ Ca 3 mM), Ca 2þ permeability (% of remaining I after BAPTA treatment), and rectification profile (I þ40 mV /I À90 mV ) for the ancestral amniote and tetrapod receptors were inferred using the ace function from the APE package v5.2 (Paradis et al. 2004) implemented in R v3.4.1 and RStudio v1.0.153. We used the Brownian motion model (Schluter et al. 1997), where characters evolve following a random walk fitted by maximum likelihood (Felsenstein 1973)

Supplementary Material
Supplementary data are available at Molecular Biology and Evolution online.