Glutamate receptors are divided in two unrelated families: ionotropic (iGluR), driving synaptic transmission, and metabotropic (mGluR), which modulate synaptic strength. The present classification of GluRs is based on vertebrate proteins and has remained unchanged for over two decades. Here we report an exhaustive phylogenetic study of GluRs in metazoans. Importantly, we demonstrate that GluRs have followed different evolutionary histories in separated animal lineages. Our analysis reveals that the present organization of iGluRs into six classes does not capture the full complexity of their evolution. Instead, we propose an organization into four subfamilies and ten classes, four of which have never been previously described. Furthermore, we report a sister class to mGluR classes I-III, class IV. We show that many unreported proteins are expressed in the nervous system, and that new Epsilon receptors form functional ligand-gated ion channels. We propose an updated classification of glutamate receptors that includes our findings.
Nerve cells or neurons communicate with each other by releasing specific molecules in the gap between them, the synapses. The sending neuron passes on messages through packets of chemicals called neurotransmitters, which are picked up by the receiving cell with the help of receptors on its surface. Neurons use different neurotransmitters to send different messages, but one of the most common ones is glutamate.
There are two families of glutamate receptors: ionotropic receptors, which can open or close ion channels in response to neurotransmitters and control the transmission of a signal, and metabotropic receptors, which are linked to a specific protein and control the strength of signal.
Our understanding of these two receptor families comes from animals with backbones, known as vertebrates. But the receptors themselves are ancient. We can trace the first family back as far as bacteria and the second back to single-celled organisms like amoebas. Vertebrates have six classes of ionotropic and three classes of metabotropic glutamate receptor. But other multi-celled animals also have these receptors, so this picture may not be complete.
Here, Ramos-Vicente et al. mapped all major lineages of animals to reveal the evolutionary history of these receptors to find out if the receptor families became more complicated as brain power increased. The results showed that the glutamate receptors found in vertebrates are only a fraction of all the types that exist. In fact, before present-day animal groups emerged, the part of the genome that holds the ionotropic receptor genes duplicated three times. This formed four receptor subfamilies, and our ancestors had all of them. Across the animal kingdom, there are ten, not six, classes of ionotropic receptors and there is an extra class of metabotropic receptors. But only two subfamilies of ionotropic and three out of four metabotropic receptor classes are still present in vertebrates today.
The current classification of glutamate receptors centers around vertebrates, ignoring other animals. But this new data could change that. A better knowledge of these new receptors could aid neuroscientists in better understanding the nervous system. And, using this technique to study other families of proteins could reveal more missing links in evolution.
Glutamate is the principal excitatory neurotransmitter in the central nervous system of animals (Fonnum, 1984; Danbolt, 2001; Pascual-Anaya and D'Aniello, 2006). It acts on two families of structurally unrelated receptors: ionotropic glutamate receptors (iGluRs), which are ligand-gated ion channels and G-protein coupled receptors (GPCRs), known as metabotropic glutamate receptors (mGluRs) (Sobolevsky et al., 2009; Conn and Pin, 1997). While fast excitatory neurotransmission is mediated by iGluRs, metabotropic receptors modulate synaptic transmission strength. iGluRs are formed by four subunits, which can be traced back to bacteria (Tikhonov and Magazanik, 2009). The current classification of iGluR subunits includes six classes: α-amino-3-hydroxy-5-methyl-4-isoxazolepropionic acid (AMPA) receptors, Kainate receptors, N-methyl-D-aspartate (NMDA) receptors (actually comprising three classes: NMDA1-3) and Delta receptors (Traynelis et al., 2010). iGluR subunits of the same class assemble into homo- or heterotetramers (Karakas and Furukawa, 2014; Kumar et al., 2011) and their ligand selectivity is dictated by a small number of residues located in the ligand-binding domain (Traynelis et al., 2010). Accordingly, NMDA subunits GluN1 and GluN3 as well as the Delta subunit GluD2 bind glycine and D-serine, while all subunits from the AMPA and Kainate classes bind glutamate (Traynelis et al., 2010; Kristensen et al., 2016). Metabotropic glutamate receptors are class C GPCRs and as such are formed by a single polypeptide. mGluRs also appeared before the emergence of metazoans, being present in unicellular organisms such as the amoeba Dictyostellium discoideum (Taniura et al., 2006). mGluRs are presently organized into three classes (I, II and III) and all their members respond to glutamate (Conn and Pin, 1997; Pin et al., 2003).
While the phylogeny of the two families of GluRs is well characterized in vertebrates, that of the entire animal kingdom is only poorly understood. The few studies on iGluR evolution outside vertebrates concentrate on a few phyla, leaving many proteins unclassified (Greer et al., 2017; Brockie et al., 2001; Janovjak et al., 2011; Kenny and Dearden, 2013). Similarly, the vast majority of mGluRs described so far fall into the three classes described in vertebrates (Krishnan et al., 2013; Kucharski et al., 2007; Dillon et al., 2006). Although, the existence of three insect mGluRs that cluster apart from classes I-III led to propose the existence of a fourth class (Mitri et al., 2004). Here we present what to our knowledge is the most comprehensive phylogenetic study of ionotropic and metabotropic GluRs along the animal kingdom. We have favored the use of more slow-evolving species for the construction of phylogenetic trees. These species are particularly amenable to phylogenetics (Simakov et al., 2013; Simakov et al., 2015; Putnam et al., 2007) as they arguably present lower rates of molecular evolution than other organisms. Our work shows that metazoan evolution of GluRs is much more complex than previously thought. iGluRs present an overall organization into four subfamilies that were already present in the last ancestor of all metazoans. Vertebrate species only retain members of two of these subfamilies. Furthermore, we identify many lineage-specific gains, losses or expansions of GluR phylogenetic groups. Finally, we present experimental evidence showing that unreported GluRs found in the basally divergent chordate Branchiostoma lanceolatum (amphioxus) are highly expressed in the nervous system and that members of the unreported Epsilon subfamily, the most phylogenetically spread among unreported groups, can form functional ligand-gated ion channels.
We have performed a systematic phylogenetic study of iGluR evolution across the animal kingdom. To increase the confidence on iGluRs evolutionary history phylogenetic trees have been generated using two independent methods (Bayesian inference and Maximum-likelihood (ML), Figure 1 and Figure 1—figure supplement 1). Our analysis indicates that the family of iGluRs experienced key duplication events that define its present organization into four previously unreported subfamilies, of which two contain the extensively studied vertebrate classes. Assuming ctenophores as the sister group to all other animals (Moroz et al., 2014; Ryan et al., 2013), our data suggest that the three major duplication events leading to this four subfamilies occurred before the divergence of current animal phyla (see Figure 2 for a summary scheme of iGluRs evolution). The first of these duplications produced the separation of the Lambda subfamily, the second lead to divergence of the NMDA subfamily and the third to the split between Epsilon and AKDF subfamilies.
Bayesian phylogeny of metazoan ionotropic glutamate receptors.
Ionotropic glutamate receptor subfamilies are indicated in colored boxes at the right. Sequences belonging to the same class are highlighted together by dashed lines and the class name is also shown. Green circles highlight the three duplications occurred before the divergence of the ctenophore lineage that lead to these four subfamilies. Posterior probabilities are shown at tree nodes and protein names at the end of each branch. Tree branches are colored based on phylum, as indicated in the legend. For unreported phylogenetic groups, names of proteins predicted to bind glycine or glutamate are highlighted in yellow or orange, respectively. Protein names from non-vertebrate species are composed of four parts: (i) 'GluR#', where # is a code denoting class or subfamily (A, AMPA; K, Kainate; F, Phi; D, Delta; Akdf, AKDF; E, Epsilon; N, NMDA and L, Lambda); (ii) a number, or range of numbers, denoting orthologous vertebrate protein(s), if any; (iii) a Greek letter to identify non-vertebrate paralogs, if any and (iv) a three-letter species code. iGluRs from A. thaliana were used as an outgroup. All information on species and proteins used is given in Figure 1—source data 2. Phylogenetic reconstruction was performed using Bayesian inference. The amino acid substitution model used was Vt + G + F, number of generations: 14269000, final standard deviation: 0.007016 and potential scale reduction factor (PSRF): 1.000. Scale bar denotes number of amino acid substitutions per site. Although the GluAkdf2_Tad protein localizes to the Delta class in this tree, we do not consider this molecule as a confident member of this class. This is because the statistical support provided by the Bayesian analysis is low and because the Maximum-likelihood analysis (see Figure 1—figure supplement 1) does not position this protein in the Delta branch.
Schematic representation of iGluRs metazoan evolution.
(a) Summary tree showing the evolution of iGluR subfamilies and classes in the metazoan lineages investigated. Each branch corresponds with one lineage. Phylogenetic subfamilies are represented by yellow boxes and classes by blue boxes. The four subfamilies present in the ancestor of all current metazoan lineages are shown at the base of the tree. Duplications of subfamilies in ancestors of current lineages are indicated. When a class or subfamily is lost in a lineage or in an ancestor, the corresponding box is crossed out with a red cross. Mollusca and annelida are lophotrochozoans and priapulida and arthropoda ecdysozoans. In priapulida NMDA2s and NMDA3s were not investigated. (b) Table indicating the presence or absence of iGluR subfamilies and classes in the metazoan lineages investigated. When a phylogenetic group is present in a lineage it is indicated by a green tick and if it is absent by a red cross. The last column shows the total number of groups found in each phylum. The last row shows the number of phyla where each phylogenetic group is present.
Image credit: Placozoa, author Oliver Voigt, licensed under CC BY-SA 3.0 Germany license; source https://commons.wikimedia.org/wiki/File:Trichoplax_mic.jpg; P caudatus, author Shunkina Ksenia, licensed under CC BY 3.0 source https://commons.wikimedia.org/wiki/File:Priapulus_caudatus.jpg; Hemichordata, released under GNU Free Documentation License, source https://commons.wikimedia.org/wiki/File:Eichelwurm.jpg; Cephalochordata, author Hans Hillewaert, licensed CC BY-SA 4.0 International license, source https://commons.wikimedia.org/wiki/File:Branchiostoma_lanceolatum.jpg.
The Lambda subfamily is the most phylogenetically restricted, as we could only identify it in porifers. Thus, Lambda would have been lost in two occasions, in the lineage of ctenophores and in a common ancestor of placozoans, cnidarians and bilaterals. On the other hand, the Epsilon subfamily is the best represented among non-bilaterians, being present in all non-bilaterian phyla investigated. Including in porifers, although we could only identify one Epsilon in sponges, GluE_Ifa from the demosponge Ircinia fasciculata. Our data also indicate that this subfamily has been lost in multiple occasions along metazoan evolution, as we could not find it in the protostome, echinoderm or vertebrate species investigated. Interestingly, all ctenophore iGluRs identified, which have been previously reported (Alberstein et al., 2015), belong to the Epsilon subfamily. Thus, this phylum would have lost NMDA, Lambda and AKDF proteins. Contrarily, ctenophores would have experienced an important expansion of Epsilon iGluRs, as we report 17 and 10 of these proteins in the two species with genomic information available, M. leidyi and P. bachei, respectively.
Although we have not identified NMDA receptors in ctenophores, porifers and placozoans our analysis indicates that this subfamily was already present in the last common ancestor of metazoans. This is because the topology of the tree shows that NMDAs appear in the phylogeny at the same level as the Epsilon subfamily, which has representatives in all non-bilateral phyla. According to our data, NMDA1s on the one hand and NMDA2s and NMDA3s on the other contain members of the cnidarian phylum. Although we have only been able to identify one member more closely related to NMDA2 and NMDA3 than NMDA1 (GluN2/3_Nve), its position in the phylogeny is very well supported by both analyses performed. This indicates that a specific duplication occurred in the ancestor of bilaterians originating NMDA2s and NMDA3s. Moreover, we have also identified a cnidarian-specific NMDA class, that we have termed NMDA-Cnidaria, this class presents representative proteins in 3 of the four species investigated. Among bilaterals we have observed conservation of all NMDA classes with the exception of NMDA2s in echinoderms, which are absent from the two species examined. Interestingly, studied cnidarian species substantially expanded their NMDA subfamily repertoire, with at least six members in Nematostella vectensis.
In bilaterians the AKDF subfamily diversified into the known AMPA, Kainate and Delta classes, but also into a fourth new class that we have termed Phi. The phylogenetic spread of these classes is quite variable, as AMPA and Kainate are in all bilateral phyla investigated but Delta and Phi are more restricted. Deltas are almost completely absent from ecdysozoan species, as we could only find a single member of this class in priapulids (P. caudatus) and none in arthropods or nematodes. Similarly, Deltas are poorly represented in mollusks and, with the available data, absent in annelids. Finally, we could only identify Phi proteins in cephalochordates, hemichordates and echinoderms, indicating that this class might be lost in the lineages of protostomes and olfactores (i.e. vertebrates and urochordates). The AKDF subfamily also includes proteins from the non-bilateral phyla of porifera, placozoa and cnidarian. The exact organization of these proteins into classes is not as straightforward as for bilateral proteins. The Bayesian and ML analysis only agree in the position of 12 iGluRs from the sponge O. carmela, these would constitute the only clear class in non-bilaterals, which we have termed AKDF-Oca.
Another example of a multiple lineage-specific event that occurred during animal evolution of iGluRs can be observed in the evolution of AMPA and Kainate proteins among protostomes. The general iGluRs phylogeny (Figure 1) suggests that ecdysozoan species have expanded their repertoire of Kainate subunits when compared with lophotrochozoans (e.g. mollusks, annelids), since C. teleta and L. gigantea only presents one and two genes coding for Kainate receptors, respectively. Contrarily, we found more AMPA subunits in lophotrochozoans than in ecdysozoan species. To investigate whether the two protostome lineages have alternatively expanded genes coding for AMPA or Kainate subunits we conducted a phylogenetic analysis of these two classes using eight species of ecdysozoans and seven of lophotrochozoans with well-characterized genomes (Figure 3 and Figure 3—figure supplement 1). Nematodes were left out of the analysis as they lack Kainate receptors (Brockie et al., 2001). This analysis retrieved 40 lophotrochozoan genes coding for AMPA subunits but only 15 coding for Kainates. The opposite scenario was observed in the genomes of ecdysozoan species, with 10 AMAP and 40 Kainate proteins,. Yet, among ecdysozoans the priapulid P. caudatus has two AMPA and two Kainate subunits, indicating that the expansion of Kainate receptors might be exclusive to arthropods. Overall the AMPA:Kainate ratio resulted to be around 1:4 in ecdysozoans and 4:1 in lophotrochozoans.
Bayesian phylogeny of AMPA and Kainate classes in protostomes.
Ionotropic glutamate receptors classes are indicated at the right. Posterior probabilities are shown at tree nodes and protein names at the end of each branch. Tree branches are colored based on phylum, as indicated in the legend. Protein names from non-vertebrate species are composed of four parts: (i) 'GluR#', where # is a one letter code denoting class (A for AMPA and K for Kainate); (ii) a number, or range of numbers, denoting orthologous vertebrate protein(s), if any; (iii) a Greek letter to identify non-vertebrate paralogues, if any and (iv) a three-letter species code. GluN1s from chordates were used as an outgroup. All information on species and proteins used in this phylogeny is given in Figure 3—source data 2. Phylogenetic reconstruction was performed using Bayesian inference. The amino acid substitution model used was Vt + I + G, number of generations: 8868000, final standard deviation: 0.0072 and potential scale reduction factor (PSRF): 1.001. Scale bar denotes number of amino acid substitutions per site.
All proteins from unreported groups (i.e. subfamilies and classes) present well-conserved sequences in iGluR domains, including transmembrane domains or residues involved in receptor tetramerization (Figure 1—figure supplement 2 and Figure 1—source data 1). Three-dimensional (3D) models of two Epsilon subunits from amphioxus (GluE1 and GluE7) indicate that their general fold is well preserved (Figure 1—figure supplement 3a). The only noticeable distinction in proteins from these groups is an insertion in the intracellular loop between the first and second transmembrane domains in Epsilon proteins. This insertion is particularly distinct in ctenophore iGluRs, having been termed as the cysteine-rich loop (Alberstein et al., 2015) (Figure 1—figure supplement 4). We have also identified a sequence difference among Epsilon proteins. Ctenophore iGluRs have two cysteines that form a disulfide bond at loop 1 of the ligand binding domain (Alberstein et al., 2015), which are also present in NMDA proteins. Nevertheless, this element is absent from the remaining members of the Epsilon subfamily.
The 'SYTANLAAF' motif, essential for channel gating (Traynelis et al., 2010), is also well conserved in most sequences, in particular the second, fourth and fifth residues (Figure 1—figure supplement 2). Nevertheless, all members of the Lambda subfamily and some proteins of the Phi class present lower levels of conservation in this sequence. Whether these changes have a functional impact is something that will require further investigation. The Q/R site (Q586, residue numbering according to mature rat GluA2) and the acidic residue located four positions downstream D/E590 (Figure 1—figure supplement 4) are involved in calcium permeability and polyamine block of AMPA and Kainate receptors (Bowie and Mayer, 1995; Koh et al., 1995; Kamboj et al., 1995). Of these two positions the latter is much better conserved, especially outside ctenophores and the Lambda subfamily. We have identified an acidic residue at position 590 in 84 out of 122 iGluRs from unreported groups, including cnidarian NMDAs. Yet, only 1/3 of these proteins present a glutamine (Q) at position 586. This includes most AKDFs and Epsilon proteins from non-ctenophores, contrarily, none of the Phi subunits presents a Q586.
The key ligand binding residues involved in fixing the amino acid backbone (α−amino and α−carboxyl) are Arg485 and an acidic residue at position 705 (Naur et al., 2007; Armstrong and Gouaux, 2000; Mayer, 2005; Furukawa et al., 2005; Yao et al., 2008). These two positions are well conserved in 94 of the 122 proteins from unreported groups, suggesting that their endogenous ligand is an amino acid (see Figure 1—figure supplement 3b for a 3D representation of ligand binding by GluE1 and Figure 1—figure supplement 5 for an alignment of iGluR residues involved in ligand binding). The residue changes found in the remaining 28 proteins would render them unable to bind an amino acid (Figure 1—figure supplement 5). This is are particularly common among class Phi proteins from amphioxus and in NMDA-Cnidaria.
Residues involved in ligand selectivity show higher variability. These are located at positions 653 and 655, and are occupied by glycine and threonine in glutamate-binding proteins and by serine and a non-polar residue in glycine-binding iGluRs. However, a recent study of ctenophore receptors has found that position 653 can be occupied by serine or threonine in glutamate-binding iGluRs, and by an arginine in glycine-binding subunits (Alberstein et al., 2015). Based on this previous knowledge we have predicted the ligand specificities of most previously unreported receptors. The preferred ligand could be confidently predicted for 72 out of the 94 proteins with well-conserved residues involved in fixing the amino acid backbone.
Interestingly, all unreported groups comprise glycine- and glutamate-specific iGluRs. Gly-specific receptors slightly outnumber those predicted to respond to glutamate (overall ratio about 3:2). The Lambda subfamily would include three proteins specific for glutamate and one for glycine, while seven remain with an unknown selectivity. Of note, the protein predicted to bind glycine (GluL5_Oca) displays an arginine at position 653, a feature which had only been reported in ctenophores (Alberstein et al., 2015). This residue would form a salt bridge with Glu423, which is key for glycine selectivity in ctenophores (Alberstein et al., 2015). Most Epsilon and AKDF proteins would preferably bind glycine, although ctenophores present a similar number of Epsilon receptors predicted to respond to glycine or glutamate (Figure 1) (Alberstein et al., 2015). In the Phi class we also found a similar number of receptors binding glycine and glutamate. Finally, we could only predict binding specificity for two of the 9 NMDA-Cnidaria proteins, as they present many changes in the residues involved in either amino acid backbone binding or side chain recognition.
Interestingly, the 22 proteins for which we could not confidently predict their ligand selectivity (Figure 1—figure supplement 5), present a limited number of residues occupying position 653 and 655, suggesting constrained evolution. Of these: (i) nine present residues with negative polarity at both positions, being candidates to bind glutamate, (ii) six present a Gly653 and a non-polar residue at position 655, and thus are candidates to bind glycine, (iii) five proteins, all from the Branchiostoma genus, present a tyrosine at position 653. A structural model of one of these receptors, GluE7 (Figure 1—figure supplement 3c), shows that a Tyr653 aromatic side chain would occupy the ligand-binding pocket, strongly suggesting that amino acid binding would be blocked. Finally, (iv) two proteins present a phenylalanine in either of the two positions and remain unclassified.
We used quantitative PCR (qPCR) to investigate gene expression levels of all iGluR subunits identified in B. lanceolatum, including those from the Epsilon and Phi groups. All 24 B. lanceolatum iGluR subunits identified in silico were found expressed in amphioxus, with the exception of Grie5 (Figure 4a). Furthermore, they all showed a significantly higher expression in the nerve cord as compared to the whole body, suggesting tissue-enriched expression. While we observed low expression levels for Epsilon genes coding for subunits with a tyrosine at position 653 (Grie5-8), which according to the 3D model would block the ligand-binding pocket, the expression of Grif1-2, also presenting the same tyrosine, reach much higher levels, comparable to those of subunits from the Kainate, Delta or NMDA classes. Thus, the presence of a tyrosine at position 653 does not appear to be directly correlated with low expression levels.
Expression and functional analysis of amphioxus iGluRs.
(a) iGluRs mRNA expression (mean and standard deviation) in Branchiostoma lanceolatum. Bars show average relative expression of B. lanceolatum (amphioxus) iGluR genes as determined by qPCR. Filled bars represent whole body and open bars nerve cord expression levels. Note that all genes show significantly enriched expression in the nerve chord relative to the whole body, with the exception of Grie5 (Student's t-test, n = 3). Expression level in the nerve chord is compared across genes of the same class. Statistics: pair comparisons were done by Student's t-test, n = 3, multiple comparisons were done by one-way ANOVA followed by Tukey's Post-Hoc test, n = 3. Significance levels: ***p < 0.001, **p < 0.01 and *p < 0.05; ns, not significant. (i) AMPA class. (ii) Kainate class. (iii) Phi class. (iv) Delta class. (v) Epsilon subfamily. (vi) NMDA classes. (b) Multiple sequence alignment of iGluRs transmembrane region M4 containing residues involved in tetramerization, these are indicated by a black frame. Higher amino acid conservation is represented by increasing intensity of blue background and by a bar chart at the bottom. Sequences included are GluE1 and GluE7 from amphioxus and representatives of human iGluRs. (c) Immunoblot of chimeric GluE1 and GluE7, containing the signal peptide from rat GluA2, expressed in HEK293T cells. Proteins were detected using the immuno-tags (c-Myc and HA, respectively) located after the rat signal peptide. Protein extracts from non-transfected cells were loaded as negative controls. (d) Immunofluorescence of HEK293T cells expressing rat GluA2 (top), cMyc-tagged GluE1 (middle) or HA-tagged GluE7 (bottom). Both non-permeabilized and permeabilized conditions are shown. (e) Immunoblot of tetrameric rat GluA2, GluE1 and GluE7 expressed in HEK293T cells. Amphioxus proteins were detected using the immuno-tags (c-Myc and HA, respectively) located at the N-terminus of each sequence. Protein extracts from non-transfected cells were loaded as negative controls.
Amphioxus genes coding for GluE1 and GluE7 were synthesized in vitro and transiently expressed in HEK293T cells for functional studies. Wild-type GluE1 and GluE7, which are not predicted to have a canonical signal peptide by SignalP 4.1 (Nielsen, 2017), expressed well but were not trafficked to the plasma membrane (Figure 4—figure supplement 1a–d), even though residues involved in tetramerization (Salussolia et al., 2013) are well conserved (Figure 4b). We thus synthesized new variants of these genes with the signal peptide from rat GluA2 (Figure 1—figure supplement 1cd). These constructs also expressed well (Figure 4c) and now were efficiently trafficked to the plasma membrane, as indicated by the staining observed in non-permeabilized cells (Figure 4d). Furthermore, analysis of receptor oligomerization, performed using non-denaturing gel electrophoresis and immunoblot, clearly indicates that both proteins form homotetramers in vitro (Figure 4e).
We next investigated the gating properties of two Epsilon proteins from amphioxus, GluE1 and GluE7. The presence of a serine and a tryptophan at positions 653 and 704, respectively, suggested that GluE1 would bind glycine. Indeed, neither glutamate nor aspartate elicited a response in our experimental settings. Instead, glycine application was able to elicit an inward whole-cell current at a membrane potential of −60 mV (Figure 5a). Interestingly, the chemically related amino acids alanine and D-serine only generated very low responses, indicating a high selectivity of the GluE1 homotetramer for glycine.
Glycine activates an amphioxus homomeric Epsilon receptor.
(a) Representative homomeric GluE1 (from B. lanceolatum) whole-cell currents evoked by a rapid pulse (500 ms) of different amino acids (10 mM) in HEK293T cells. Left and right glycine-mediated currents denote agonist application before and after alanine, D-serine and glutamate applications respectively for ruling out run-down of the currents. (b) Representative GluE1 responses to 10 mM glycine at different membrane voltages (from −80 to +80 mV in 20 mV steps). Note that a strong inward rectification can be observed even in the absence of added polyamines in the intracellular solution. Inset: currents at negative membrane voltages are shown. (c) Current-voltage relationship for peak currents evoked by glycine (500 ms, 10 mM) applied to whole HEK293T cells containing homomeric GluE1 subunits normalized for the current at −80 mV (n = 3) fitted to a 5th order polynomial function. Error bars represent SEM. (d) Homomeric GluE1 glycine-mediated currents recorded at different time intervals by using a paired pulse protocol. (e) Rate of recovery of desensitization fitted to a single exponential of time constant 10.8 s (n = 3–5). Plot shows the average ratio values (P2/P1) and SEM (error bars).
The Epsilon receptor displayed a strong inward rectification, even in the absence of added polyamines in the intracellular solution (Figure 5b,c). This behavior is characteristic of unedited AMPA and Kainate receptors displaying a glutamine (Q) and an acidic residue at positions 586 and 590, respectively (Bowie and Mayer, 1995; Koh et al., 1995; Kamboj et al., 1995) and GluE1 presents a glutamine and an aspartic acid at these positions (Figure 1—figure supplement 4). Glycine-mediated currents showed a slow rate of recovery from desensitization when compared with AMPA or Kainate mammalian receptors, requiring 20–25 seconds until a complete recovery was achieved and a full response of the same magnitude could be recorded (Figure 5d,e). Similar observations have been made with ctenophore receptors activated by glycine in which the recovery from desensitization has an unusually long time constant of 81 seconds (Alberstein et al., 2015).
Finally, functional studies on receptors formed by GluE7 did not retrieve any positive results. None of the following amino acids: glutamate, aspartate, asparagine, glycine, alanine or D-serine elicited a response in our experimental system. We hypothesize that, as predicted by the 3D model, the presence of a tyrosine at position 653 renders a homomeric form of this receptor unable to function as an amino acid-gated ion channel.
We next performed a phylogenetic study of metabotropic glutamate receptors (Figure 6 and Figure 6—figure supplement 1). This analysis has revealed that the three historical mGluR classes (I to III) have a sister group. Following the current nomenclature we have named this as class IV. The existence of this class had already been proposed on the bases of three insect proteins (Mitri et al., 2004). Yet, here we show that this class is actually present in all bilateral phyla, excluding vertebrates. Furthermore, we also show that class IV appeared together with classes I-III before radiation of bilateral lineages. We have identified clear orthologues to class I-IV in porifers, placozoans and cnidarians but not in ctenophores. These are organized into four classes, two from cnidarians, and one from placozoans and porifers (Figure 6). We have also identified non-bilaterian mGluRs that fall outside the above-mentioned classes. Unfortunately, the Bayesian and ML phylogenies do not agree on the exact organization of these early divergent mGluRs, except for the fact that they diverge prior to bilaterian classes. For this reason we have left these sequences unclassified. Whether these sequences belong to one, or even multiple classes that would have been lost in bilateral organisms is something that will require further investigation.
Bayesian phylogeny of metazoan metabotropic glutamate receptors.
Identified metabotropic glutamate receptor classes from bilateral and non-bilateral organisms are indicated by colored boxes at the right. Dashed boxes further highlight individual classes from bilateral organism. Posterior probabilities are shown at tree nodes and protein names at the end of each branch. Tree branches are colored based on phylum, as indicated in the legend. Protein names from non-vertebrate species are composed of four parts: (i) 'mGluR', followed by a number, or range of numbers, denoting orthologous vertebrate protein(s), if any (for Class IV and group I-II-III-IV proteins, the name is followed by the name of the class/group); (ii) a Greek letter to identify non-vertebrate paralogs, if any and (iv) a three-letter species code. GABA-B receptors from vertebrates were used as an outgroup. All information on species and proteins used in this phylogeny is given in Figure 6—source data 2. Phylogenetic reconstruction was performed using Bayesian inference. The amino acid substitution model used was WAG + I + G + F, number of generations: 5327000, final standard deviation: 0.004788 and potential scale reduction factor (PSRF): 1.001. Scale bar denotes number of amino acid substitutions per site.
Although all class IV proteins show well conserved sequences overall (Figure 6—figure supplement 2a, Figure 6—figure supplement 3 and Figure 6—source data 1), two residues critical for glutamate binding, Arg78 and Lys409, are non-conservatively replaced by non-polar or acidic residues in all class IV proteins identified (Figure 6—figure supplement 2a, residue numbering corresponds to human mGluR1). These changes are predicted to hamper glutamate binding and, indeed, functional studies of a class IV receptor from fruit fly indicated that it does not respond to this amino acid (Mitri et al., 2004). All class IV proteins would share this feature. On the other hand, residues involved in contacts with the amino acid backbone are well conserved (Figure 6—figure supplement 2a), suggesting that these proteins might bind an amino acid other than glutamate. Similarly, mGluR residues from most non-bilaterian sequences involved in binding the amino acid backbone are highly conserved. Among non-bilaterian proteins the residues involved in glutamate binding are only conserved in approximately half of the proteins from classes orthologous to I-II-III-IV. Finally, we investigated mGluRs expression in amphioxus following the same procedure described for iGluRs. All five amphioxus mGluRs showed an enriched expression in the nerve cord, including the two class IV genes. Noticeably, these two genes showed significantly higher expression levels than orthologues of vertebrate classes (Figure 6—figure supplement 2b).
We have performed what to our knowledge is the most comprehensive phylogenetic study of metazoan glutamate receptors. This has revealed that their evolutionary history is much more complex than what is currently acknowledged, especially for the family of iGluRs. Our study has also revealed the existence of unreported phylogenetic groups in both ionotropic and metabotropic glutamate receptors. Importantly, our data indicate that the evolution of glutamate receptors has not occurred in an unequivocal incremental manner only in those clades with more elaborated neural systems, but it has rather followed an scattered lineage-specific evolutionary history. This means that certain lineages have experienced the gain, loss, expansion or reduction of specific phylogenetic groups.
Our phylogenetic analysis indicates that the family of iGluRs is actually divided into four unreported subfamilies that we have termed Lambda, Epsilon, NMDA and AKDF. Interestingly, this general organization was already present in the last common ancestor of all metazoans and later duplications within NMDA and AKDF subfamilies resulted in the formation of well-known iGluR classes. The other two subfamilies are absent from the majority of model species used in neuroscience research. The NMDA subfamily diversified into classes NMDA1-3 but also into the NMDA2/3 and NMDA-Cnidaria. Similarly, the AKDF subfamily diversified into the AMPA, Kainate and Delta classes, but also into the previously unreported Phi class. We have also identified and AKDF class exclusive to porifers, represented by sequences form O. carmela. Most well-studied iGluR classes are the result of duplications in ancestors of current bilateral species, >650 million years ago (mya) (Kumar et al., 2017), only class NMDA1 originated earlier, as cnidarians present members within this class. The Epsilon subfamily, which includes all iGluRs from ctenophores, is the only subfamily present in all non-bilateral phyla investigated, including sponges. It is thus the subfamily presenting a larger phylogenetic spread, as it is also present in hemichordates and in non-vertebrate chordates. On the other hand, the unreported Phi class shows a more restricted phylogenetic spread, as it is present only in three deuterostome phyla. Moreover, Lambda proteins seem restricted to Porifers, which constitutes an interesting evolutionary case due to maintenance of a glutamate receptor family in a phylum without nervous system.
The phylogenetic analysis of metabotropic glutamate receptors has allowed us to unambiguously establish the existence of a sister group to the well-known classes I, II and III. Following the present nomenclature we have named this as class IV. This class had been previously proposed based on the identification of three insect mGluRs that did not cluster with members of known classes (Mitri et al., 2004). Here we show that class IV is not restricted to insects, but is actually present in all bilaterian phyla investigated, with the exception of vertebrates where this class has been lost. Interestingly, as it occurs for most well-known iGluR classes, mGluR classes I-IV appeared simultaneously in the ancestor of bilaterals. Our phylogenetic analysis also indicates that the non-bilateral phyla of cnidarians, placozoans and porifers present clear orthologues to classes I-IV, which are organized into four classes, while we failed to find any in the early-branching ctenophores. Finally, we were unable to confidently classify many non-bilateral mGluRs, which might constitute one or more classes.
We have identified many examples of lineage-specific evolutionary events. These would antagonize with a model in which species with less elaborated nervous systems would present GluR families with lower complexity. The most noticeable examples are: (i) the absence of all subfamilies but Epsilon in analyzed ctenophores, (ii) the loss of Delta receptors from arthropods, nematodes and annelid species investigated, (iii) the loss of the Epsilon subfamily in vertebrates, echinoderms and protostomes, (iv) the loss of the Phi class in vertebrates and studied protostomes, (v) the specific expansion of Kainate receptors in arthropods, which contrasts with the expansion of AMPA receptors in its sister lineages of mollusks and annelids, (vi) the large expansion of the Epsilon subfamily in ctenophores, placozoans and cephalochordates and, finally (vii) the loss of mGluR class IV in vertebrates.
Along the same line, it is interesting to note that amphioxus (B. belcheri and B. lanceolatum), with a simple nervous system, have over 20 genes encoding iGluRs, while mammals have 18. Other non-vertebrate species also present large numbers of iGluRs, including the 19 iGluRs identified in the sponge O. carmela or the 17 present in the ctenophore M. leidyi, to mention a few. Similarly, the cnidarian A. digitifera and the ctenophore M. leidyi have seven mGluRs each, while the placozoan T. adhaerens presents eleven, three more than the eight mGluRs found in the human genome. The large number of GluRs found in many non-vertebrate animals suggests that there has been an evolutionary trend to increase their number in many metazoan lineages.
Our experimental results suggest that unreported receptors would play a role in the nervous system, as Epsilon, Phi and mGluR class IV genes are highly expressed in the nerve cord of amphioxus. Nevertheless, whether all these proteins are expressed at the synapse and act as neurotransmitter receptors is an issue that will require further investigation. Their presence in other tissues, such as sensory organs, cannot be ruled out. Those receptors showing more divergent sequences, particularly in residues involved in ligand binding, might respond to other molecules. For instance, they could behave as chemoreceptors, as it is the case of antennal receptors found in insects (Croset et al., 2010; Benton et al., 2009).
Proteins from all unreported groups generally present a good conservation of residues involved in binding the amino acid backbone, indicating that their ligand would be an amino acid or a closely related molecule. Interestingly, we could identify proteins predicted to bind either glycine or glutamate in all unreported iGluR subfamilies and classes. If our functional predictions are correct, the ability to recognize one or the other amino acid would have emerged repeatedly in all unreported iGluR phylogenetic groups. Unexpectedly, the nature of the residues conferring amino acid specificity indicates that only a minority of proteins from unreported GluR groups would respond to glutamate. Sequence analysis and structural considerations strongly suggest that class IV mGluRs will not bind glutamate and that among non-bilateral mGluRs only a minority, belonging to classes orthologous to I-II-III-IV, are predicted to bind to this neurotransmitter. Similarly, among unreported iGluR groups, the number of proteins binding glycine outnumbers those binding glutamate. Interestingly, we report a glycine-binding poriferan protein (GluL5_Oca) with a structural feature that had only been reported in ctenophores (Alberstein et al., 2015). This is an Arg653 that through establishing a salt bridge with Glu423 confers glycine specificity (Alberstein et al., 2015). We thus report that this structural element is not exclusive to ctenophores. We have also identified iGluR subunits with important changes in critical ligand binding residues, indicating that they might have evolved new biological functions, for example, response to other, as yet unidentified small molecules.
The activation of Epsilon receptors by glycine has been experimentally corroborated by electrophysiological analysis of homotetrameric receptors composed by GluE7 from M. leidy (Alberstein et al., 2015) and GluE1 from amphioxus (this study). In our hands the amphioxus receptor showed a very high selectivity for glycine, since ion currents could not be elicited by chemically related amino acids such as serine or alanine. Glycine-binding Epsilon subunits from phyla other than ctenophores present structural features similar to those from glycine-binding iGluRs in vertebrates. The greater number of glycine receptors found in non-vertebrate species could be related to the higher abundance of this amino acid in their nerve cord as compared with the mammalian brain (Pascual-Anaya and D'Aniello, 2006).
Altogether, our phylogenetic analysis and experimental findings have uncovered the complex evolution of glutamate receptors within the metazoan kingdom. Our data indicate that the classification of iGluRs is not restricted to the six classes currently recognized. Instead, iGluRs are organized into four subfamilies: Lambda, Epsilon, NMDA and AKDF and ten classes with varying phylogenetic spread. With the data available, the NMDA subfamily is organized into classes NMDA1, NMDA 2, NMDA3, NMDA-Cnidaria and NMDA2/3, while subfamily AKDF contains classes AMPA, Kainate, Delta, Phi and AKDF-Oca. Both NMDA2/3 and AKDF-Oca are represented by sequences from only one species, further sequencing of non-bilateral species will be required to fully demonstrate their existence. Furthermore, the evolution of mGluRs has generated a sister group to classes I, II and III, class IV. We have also identified classes of non-bilaterian mGluRs orthologous to I-II-III-IV. We propose that the classification of these two families of GluRs, key to the physiology of the nervous system, has to be updated to include our findings.
Phylogenetic analysis were performed with sequences from at least two species from each of the following metazoan phyla: Porifera, Ctenophora, Placozoa, Cnidaria, Lophotrochozoa, Ecdysozoa, Hemichordata, Chordata and Vertebrata, with the exception of placozoans for which only one species is available. When possible, we chose slowly evolving species. The complete lists of species used for iGluR phylogenies are given in Figure 1—source data 2. Species used in the phylogeny of metabotropic glutamate receptors are listed in Figure 6—source data 2. Sponge sequences were taken from (Riesgo et al., 2014), B. lanceolatum sequences were retrieved from unpublished genomic and transcriptomic databases (access was kindly provided by the Mediterranean Amphioxus Genome Consortium), A. digitifera and P. flava sequences were obtained from the Marine Genomics Unit (Simakov et al., 2015; Shinzato et al., 2011) and P. bachei sequences from NeuroBase (Moroz et al., 2014).
GluR sequences were identified using homology-based searches in a two-tier approach. Mouse glutamate receptors were used as search queries (iGluRs: Gria1-4; Grik1-5; Grid1-2, Grin1, Grin2A-D and Grin-3A-B; mGluRs: mGluR1-8). In a first search GluR homologs were identified using the BLASTP tool (Altschul et al., 1990) with default parameters. Subject sequences with an E-value below 0.05 were selected as candidate homologs. These were re-blasted against the NCBI database of 'non-redundant protein sequences' using the same BLAST tool. If the first hit obtained in the reciprocal BLAST was a glutamate receptor the sequence was included in the phylogenetic analysis. In a second stage the same mouse sequences were used to perform TBLASTN searches against genomic and, when available, transcriptomic databases. Subject sequences not identified in the first tear and having an E-value below 0.05 were selected as candidate homologs. These were re-blasted using BLASTX against the NCBI 'non-redundant protein sequences' database. Finally, if the first hit of this search was a glutamate receptor the sequence was also included in the phylogenetic analysis. Identified iGluR sequences in which less than four residues of the SYTANLAAF motif (Traynelis et al., 2010) were conserved were not considered for the final phylogenetic analysis. mGluR sequences lacking two or more of the seven transmembrane regions were also discarded. The complete reference lists of all iGluRs used in the final phylogeny are given in files Figure 1—source data 2. The reference list of metabotropic glutamate receptors is presented in Figure 6—source data 2. The alignments used for the phylogenetic analysis of iGluRs, mGluRs and AMPAs and Kainates from protostomes are provided in Figure 1—source data 3, Figure 3—source data 1 and Figure 6—source data 3.
The iGluR tree was constructed with 224 sequences identified in 26 non-vertebrate species (Figure 1—source data 2). The tree also included 18 iGluR sequences from vertebrates and two iGluR proteins from A. thaliana, used as an outgroup (Chiu et al., 2002). The phylogenetic analysis of AMPA and Kainate classes in protostomes was inferred using 110 sequences from 15 protostome species (Figure 3—source data 2) and 37 sequences from deuterostomes, of which 4 GluN1 proteins were used as an outgroup. The mGluR tree was constructed with 149 proteins from 29 non-vertebrate species, 38 mGluRs from vertebrate species and 10 sequences from vertebrate metabotropic GABA receptors, used as an outgroup (Figure 6—source data 2).
Protein sequences were aligned with the MUSCLE algorithm (Edgar, 2004), included in the software package MEGA6 (Tamura et al., 2013) with default parameters. ProtTest v3.4.2 was used to establish the best evolutionary model (Darriba et al., 2011). Trees were constructed using MrBayes v3.2.6 (Ronquist et al., 2012) for Bayesian inference and IQ-TREE (Nguyen et al., 2015) for Maximum-likelihood analysis. For Bayesian inference phylogenies were stopped when standard deviation was below 0.01 and its value was fluctuating but not decreasing. Markov chain Monte Carlo (MCMC) was used to approximate the posterior probability of the Bayesian trees. Bayesian analyses included two independent MCMC runs, each using four parallel chains composed of three heated and one cold chain. Twenty-five % of initial trees were discarded as burn-in. Convergence was assessed when potential scale reduction factor (PSRF) value was between 1.002 and 1.000. In Maximum-likelihood analysis the starting tree was estimated using a neighbor-joining method and branch support was obtained after 1000 iterations of ultrafast bootstrapping (Hoang et al., 2018). Gene/protein names were given based on their position in the tree. Phylogenetic trees were rendered using FigTree (http://tree.bio.ed.ac.uk/software/figtree/). Phylogenetic calculations were performed at the IBB - UAB heterogeneous computer cluster 'Celler' and at the CIPRES science gateway (RRID: SCR_008439) (Miller et al., 2010).
Branchiostoma lanceolatum adults were collected in the bay of Argelès-sur-Mer, France (latitude 42° 32' 53' N and longitude 3° 03' 27' E) with a specific permission delivered by the Prefect of Region Provence Alpes Côte d'Azur. B. lanceolatum is not a protected species. Animals were kept in tanks with seawater at 17°C under natural photoperiod.
Adult amphioxus (B. lanceolatum) were anesthetized in 0.1% diethyl pyrocarbonate (DEPC; Sigma, D5758) PBS buffer. Animals were sacrificed by cutting the most anterior part of the body. The nerve chord was surgically extracted from the animal while submerged in DEPC-PBS using a magnifying glass. Individual nerve chords were snap frozen in liquid nitrogen and stored at −80°C until use. RNA was extracted from whole animals or from dissected nerve chords. Ten nerve chords were used for each RNA extraction, so that biological variability between individuals could be normalized. The tissue was homogenized in 1 mL of TRI Reagent (Sigma, T9424) using a Polytron homogenizer. Homogenates were transferred into an Eppendorf tube and incubated 5 min at room temperature (RT) before adding 100 µL of 1-bromo-3-cloropropane. Tubes were vigorously mixed by vortexing for 10–15 s, incubated 15 min at RT and centrifuged at 13000 rpm for 15 min at 4°C. RNA was precipitated from the aqueous phase with 500 µL of isopropanol and 20 µg of glycogen. Tubes were frozen for 1 hr at −80°C and then thawed, incubated at RT for 10 min and centrifuged at 13000 rpm for 10 min at 4°C. The RNA pellet was washed twice with 500 µL of 75% ethanol and air-dried. cDNA was synthesized from 0.5 µg of total RNA. One µL of Oligo(dT)15 (Promega), 1 µL of 10 mM dNTP mix (Biotools), RNA and DEPC distilled water were mixed in a PCR tube to a final volume of 14 µL. This mix was incubated at 65°C for 5 min in a T100 Thermal Cycler (BioRad). After cooling tubes on ice for 1 min, we added 4 µL of First Strand 5x buffer, 1 µL of 0.1 M DTT and 1 µL of SuperScript III (Invitrogen). Tubes were placed in a T100 Thermal Cycler (BioRad) with the following program: 60 min at 50°C, 15 min at 70°C. RNA expression levels were determined using qPCR and the GAPDH gene used as a reference. Primers used for qPCR analysis of iGluRs are in Figure 4—figure supplement 2 and those used for mGluR qPCR in Figure 6—figure supplement 4. qPCR data for iGluRs and mGluRs are given in Figure 4—source data 1 and Figure 6—source data 4, respectively.
cDNA from nerve chord and whole body samples was diluted 1:10 for the glutamate receptor gene reactions, and 1:100 for the reference gene reaction. For each gene 2.5 µL of diluted cDNA were added to 5 μL of iTaq Universal SYBR Green Supermix (Bio-Rad), along with 0.5 µL of each primer and 1.5 µL of RNase free water. qPCR was run in a C1000 Touch thermocycler combined with the optic module CFX96. Three technical replicates were performed for all genes analyzed. Primer pairs were designed to detect the expression levels of each glutamate receptor (Figure 4—figure supplement 2 and Figure 6—figure supplement 4). B. belcheri glutamate receptor sequences were aligned with the genomic sequence of B. lanceolatum, and high identity fragments were used to design primers. All primers were 20–25 base pair long, had GC content over 40–45% and a Tm between 60–65°C. Primers were designed to obtain amplicons between 140–270 base pairs. Values of normalized expression were statistically analyzed using GraphPad Prism5. No outliers were identified and no data points were excluded. Comparisons between whole body and nerve chord expression levels were done with Student's T-Test for unpaired samples or the Welch variant of the Student's T-Test for samples with different variance. For multiple comparisons between the expression levels of genes belonging to the same class one-way ANOVA analysis was performed using Tukey's Post-Hoc test.
Grie1 and Grie7 genes were selected for transient expression in the mammalian cell line HEK293T. We prepared two constructs for each gene. We first introduced an immuno-tag in the N-terminus before the first element of secondary structure. For Grie1 we used the c-Myc tag, which was placed after residue 39, and for Grie7 we used the hemagglutinin (HA) tag introduced after residue 10 of the wild-type sequence. The second set of constructs prepared substituted the wild type N-terminal sequence for the signal peptide from rat GluA2 while maintaining the immuno-tags (Figure 4—figure supplement 1). Codon-optimized genes for expression in human cells were synthesized and cloned into pICherryNeo (Addgene, 52119) and pIRES2_EGFP (Addgene 6029–1) by the Invitrogen GeneArt Gene Synthesis service.
All expression experiments were done with a mycoplasma-free HEK293T cell line kindly provided by Prof. F. Ciruela (Universitat de Barcelona) and purchased from the American Type Culture Collection (ATCC, CRL-3216, RRID: CVCL_0063). The ATCC has confirmed the identity of HEK293T by STR profiling (STR Profile; CSF1PO: 11,12; D13S317: 12,14; D16S539: 9,13; D5S818: 8,9; D7S820: 11; TH01: 7, 9.3; TPOX: 11; vWA: 16,19; Amelogenin: X). After the purchase of the cell line, mycoplasma tests are performed in the laboratory on every new defrosted aliquot. The kit used for mycoplasma detection is PlasmoTest (Invivogen, code: rep-pt1).
HEK293T cells were maintained in Dulbecco's Modified Eagle Medium (DMEM) supplemented with 10% FBS and 1% Antibotic-Antimycotic (Gibco) in a humidified incubator at 5% CO2 air and 37°C. The day before transfection, cells were plated onto poly-D-lysine coated coverslips in 6-well plates, to reach 60–80% confluence. HEK293T cells were transiently transfected with the following plasmids: empty pIRES2-EGFP, pIRES2-EGFP containing the Grie7_Bbe gene, empty pICherryNeo and pICherryNeo containing Grie1_Bla. Cells were transfected using 3 μg of polyethylenimine and 1 μg of plasmid DNA for each ml of non-supplemented DMEM. Cells were incubated 4–5 hr with transfection medium without supplementation, which was then removed and replaced by supplemented medium. Twenty-four hours after transfection the medium was removed and cells were washed 3 times with PBS. For surface receptor staining, cells were blocked in 2% BSA in PBS for 10 min at 37°C, and incubated for 25 min at 37°C with primary antibodies against HA (Covance, MMS-101P, RRID: AB_291259), c-Myc (Cell Signalling, 2272S, RRID: AB_10692100) or GluA2 (Millipore, MAB397, RRID: AB_2113875). HA and GluA2 antibodies were diluted 1:200 and c-Myc 1:100 in DMEM without supplementation. Cells were washed 3 times with PBS, fixed in 4% paraformaldehyde (PFA) for 15 min at RT, rinsed in PBS and incubated 1 hr at 37°C with secondary antibodies Alexa Fluor 555 donkey anti-mouse IgG (H + L) (A-31570, Invitrogen, RRID: AB_2536180) and Alexa Fluor 647 goat anti-rabbit IgG (H + L) highly cross-adsorbed (Life Technologies, A-21245, RRID: AB_2535813), diluted 1:1000 and 1:500 in PBS, respectively. Finally, coverslips were washed and mounted onto slides with Fluoroshield with DAPI (Sigma-Aldrich, F6057). For intracellular labeling cells were first fixed in 4% PFA for 15 min at RT, permeabilized with 0.2% Triton X-100 in PBS for 10 min, and finally blocked with PBS containing 2% BSA and 0.2% Triton X-100 for 20 min. Primary antibodies against HA (Covance, MMS-101P, RRID: AB_291259) and GluA2 (Millipore, MAB397, RRID: AB_2113875) were diluted 1:1000 and c-Myc (Cell Signalling, 2272S, RRID: AB_10692100) antibody was prepared at 1:100 in PBS. Incubation lasted 25 min at 37°C. Secondary antibody incubations and coverslip mounting were done in the same way as for non-permeabilized cells. Cells were examined using a confocal laser-scanning microscope (Zeiss LSM 700) with a 63x oil objective.
HEK293T cells were grown in 6-well plates as described previously and transfected with plasmids expressing amphioxus GluE1, GluE7 or GluA2. Twenty-four hours after transfection cells were rinsed with PBS and the content of 4 wells was resuspended in solubilization buffer (PBS containing 2% N-dodecyl-α-maltopyranoside (DDM; D310HA, Anatrace) and the protease inhibitors mix cOmplete EDTA-free Protease Inhibitor Cocktail, Roche). Cell lysates were homogenized in a Dounce homogenizer in ice with 20 strokes and kept under orbital agitation for 1 hr at 4°C. Lysates were centrifuged at 89000xg in a Beckman TLA120.2 rotor for 40 min at 4°C. The supernatant containing solubilized membrane proteins was recovered in a new tube and stored at −20°C until used.
For native gel electrophoresis proteins were resolved in a Mini-PROTEAN TGX Gel 4–20% (Bio-Rad). Samples were mixed with Native Sample Buffer (Bio-Rad) and run along with HiMark Pre-Stained Protein Standard (Life Technologies). Electrophoresis was performed in ice at a constant voltage of 100 V for 180 min. Gels were transferred at constant current (35 mA) to polyvinylidene fluoride (PVDF) membranes overnight (16–18 hr) at 4°C. After transfer, membranes were blocked for 1 hr with Odyssey Blocking Buffer (Li-cor) in TBS, and incubated overnight at 4°C with primary antibodies anti-HA (Covance, MMS-101P, RRID: AB_291259), anti-c-Myc (Cell Signaling, 2272S, RRID: AB_10692100) or anti-GluA2 (Millipore, MAB397, RRID: AB_2113875) diluted 1:1000 in TTBS (TBS containing 0.05% Tween-20). After three 15 min washes in TTBS, membranes were incubated with donkey anti-mouse (Li-cor, 926–32212, RRID: AB_621847) and donkey anti-rabbit (Li-cor, 926–68073, RRID: AB_10954442) diluted 1:7500 in TTBS for 1 hr. Blots were analyzed in an Odyssey scanner (Li-cor).
For denaturing gel electrophoresis (SDS-PAGE) protein lysates were denatured by adding loading sample buffer 10x (500 mM Tris-HCl pH 7.4, 20% SDS, 10% β-mercaptoethanol, 10% glycerol and 0.04% bromophenol blue), and incubated for 5 min at 95°C. Protein lysates were loaded in a 10% SDS- polyacrylamide gel and separated at a constant current (25 mA). Gels were transferred at a constant voltage of 100 V for 90 min in ice. Membranes were blocked for 1 hr with Odyssey Blocking Buffer in TBS, and incubated overnight at 4°C with the same primary antibodies at the same dilution as for native gels in TBS containing 0.1% Tween 20. After three 15 min washes in TTBS, membranes were incubated with secondary antibodies as above. Blots were analyzed in an Odyssey scanner.
Models for full-length GluE1 and GluE7 were generated with RaptorX (Källberg et al., 2012) based on deposited three-dimensional crystal structures of the full-length AMPA-subtype ionotropic glutamate receptor from Rattus norvegicus, GluA2, bound to competitive antagonists (PDB codes 4U4G (Yelshanskaya et al., 2014) and 3KG2 (Sobolevsky et al., 2009), respectively). Models of their respective ligand binding domains were generated with SWISS-MODEL (Biasini et al., 2014) using the atomic-resolution crystal structure of the rat GluA2 LBD bound to glutamate as template (PDB code 4YU0). Model quality was assessed with MolProbity (http://molprobity.biochem.duke.edu/, RRID: SCR_014226). MolProbity scores for all models are given in Figure 1—source data 4. Models were inspected with MIFit (Smith, 2010) and figures were prepared with PyMOL (www.pymol.org).
Cells were visualized with an inverted epifluorescence microscope (AxioVert A.1, Zeiss) and were constantly perfused at 22–25°C with an extracellular solution containing (in mM): 145 NaCl, 2.5 KCl, 2 CaCl2, 1 MgCl2, 10 HEPES and 10 glucose (pH = 7.42 with NaOH; 305 mOsm/Kg). Microelectrodes were filled with an intracellular solution containing (in mM): 145 CsCl, 2.5 NaCl, 1 Cs-EGTA, 4 MgATP, 10 HEPES (pH = 7.2 with CsOH; 295 mOsm/Kg). Electrodes were fabricated from borosilicate glass (1.5 mm o.d., 1.16 i.d., Harvard Apparatus) pulled with a P-97 horizontal puller (Sutter Instruments) and polished with a forge (MF-830, Narishige) to a final resistance of 2–4 MΩ. Currents were recorded with an Axopatch 200B amplifier filtered at 1 KHz and digitized at 5 KHz using Digidata 1440A interface with pClamp 10 software (Molecular Devices Corporation).
Whole-cell macroscopic currents were recorded from isolated or coupled pairs of mCherry or EGFP positive HEK293T cells. Rapid application (<1 ms exchange) of agonists (500 ms pulses) at a membrane potential of −60 mV was achieved by means of a theta-barrel tool (1.5 mm o.d.; Sutter Instruments) coupled to a piezoelectric translator (P-601.30; Physik Instrumente). One barrel contained extracellular solution diluted to 96% with H2O and the other barrel contained 10 mM of the amino acid solution. For measuring current-voltage relationships, 500 ms agonist jumps were applied at different membrane voltages (−80 mV to +80 mV in 20 mV steps) and peak currents were fitted to a 5th order polynomial function. To study recovery from desensitization, a two-pulse protocol (500 ms each) was used in which a first pulse was applied followed by a second pulse at different time intervals (from 2.5 s to 25 s). The paired pulses were separated 30–60 s to allow full recovery from desensitization. To estimate the percentage of recovery, the magnitude of peak current at the second pulse (P2) was compared with the first one (P1). Electrophysiological recordings were analyzed using IGOR Pro (Wavemetrics Inc.) with NeuroMatic (Jason Rothman, UCL, RRID: SCR_004186).
The Journal of Neuroscience 21:1510–1522.
Molecular Biology and Evolution 19:1066–1082.
Molecular Biology and Evolution 35:518–522.
The Journal of Physiology 486 ( Pt 2:297–303.
The Journal of Physiology 486 ( Pt 2:305–312.
Molecular Biology and Evolution 34:1812–1819.
Proceedings of the Gateway Computing Environments Workshop 14:1–8.
The Journal of Biological Chemistry 279:9313–9320.
59–73, Predicting secretory proteins with SignalP, 1611, Methods in Molecular Biology, New York, Springer.
International Journal of Biological Sciences 2:87–92.
Molecular Biology and Evolution 31:1102–1120.
Journal of Biological Chemistry 281:12336–12343.
Neuroscience and Behavioral Physiology 39:763–773.
Thank you for submitting your article "Evolution of Glutamate Receptors Reveals Three Unreported Classes and Divergent Phyla-Specific Adaptations" for consideration by eLife. Your article has been reviewed by three peer reviewers, including Leon D. Islas as the Reviewing Editor and Reviewer #1, and the evaluation has been overseen by Richard Aldrich as the Senior Editor. The following individuals involved in review of your submission have agreed to reveal their identity: Mark L Mayer (Reviewer #3).
This manuscript presents an extensive phylogenetic analysis of ionotropic and metabotropic glutamate receptors. The authors find two new families of iGluRs a and one family of mGluRs. Based on these new families, they also propose a new evolutionary tree for both classes of receptors. This work also provides functional data on one of the newly found families, AKD1 and show that they can work as excitatory glutamate receptors and are expressed in the nervous system of amphioxus, an early vertebrate.
Three reviewers have seen your manuscript. Although there is enthusiasm for the new findings reported, the reviewers found that there are major deficiencies with the analysis that led to the proposal of the new families of ionotropic glutamate receptors. In particular, the fact that the proposed AK1 family only shows up in one of the genealogies seems troubling and might be an indication that the methods used to get the phylogenies might not be optimal. The consensus among reviewers is that the phylogenetic analysis should be repeated using more robust methods (see the reviews), in order to discard the possibility that the AK1 group might be an artifact. If this is the case, it is possible that the AKD1 group and the new metabotropic receptor family will still be present, which is a significant finding. The proposed changes should be carried out in order for the manuscript to be reconsidered.
1) If amino acid sequences were used for phylogeny reconstruction (after alignment of nucleotide sequences), then this needs to be clarified in the methods. However, if DNA sequences were used during phylogeny reconstruction, then both phylogenetic analyses presented in Figure 1 need to be repeated to make sure the proposed new classes of receptor are recovered.
2) For the expression experiments of AKD1 receptors, authors need to clarify if maybe these receptors do no not express functionally, since they seem to lack signal peptides.
3) The analysis done to claim loss of certain receptors in some phyla (e.g. Cnidarians) need to be redone and the interpretation needs to be in line with the results of the analysis, since as it stands, seems overreaching.
The paper by Ramos-Vicente and collaborators presents bioinformatics and functional data pertaining the discovery of new classes of both ionotropic and metabotropic glutamate receptors. The authors have incorporated these new families into a revised evolutionary tree of glutamate receptors. These findings are exciting and expand our knowledge of the evolution of these important membrane proteins.
In the Introduction the following phrase needs clarification: Thesespecies present a genomic structure, gene sequence and gene contentresembling that of their ancestors…. We do not know the genome of their ancestors and therefore this is an inference.
Structural models of two proteins from the new AKD1 class from amphioxus where constructed using known iGluR structures. Although this is at first glance reasonable, a less biased approach should be taken. Models should also be produced letting the software choose the appropriate structural template. Also, the MolProbity scores should be given as well as a discussion of the probability that the sequences adjust to the templates.
In the Discussion section, there is no correlation between the number of iGluR genes and the nervous system complexity of an organism. Is it possible that this might reflect the fact that in lower organisms glutamate receptors also have non-synaptic functions, as chemoreceptors?
Please include more information on the iGluRs found on A. thaliana and their relation to the newly proposed phylogeny of iGluRs.
In the figures that plot mean values, it is not specified what measure of variability is being used (i.e. standard error, variance, etc.).
In the manuscript "Evolution of Glutamate Receptors Reveals Three Unreported Classes and Divergent Phyla-Specific Adaptations", the authors conduct a phylogenetic analysis of animal glutamate receptors (both ionotropic and metabotropic) and claim to have found 3 new phylogenetically-distinct classes (2 ionotropic and 1 metabotropic). They propose that the classification mechanism for glutamate receptors must be revised to account for these new gene families. Further, they provide data from structural modeling and functional expression to argue that ionotropic receptors in one of these classes, AKD1, are gated by glycine and not glutamate.
While the manuscript addresses and interesting subject worthy of a high profile publication, there are issues with the phylogenetic analysis that must be addressed or clarified to understand how robustly the claims for new classes of ionotropic glutamate receptors are supported.
1) There appears to be a major flaw in the methods for constructing the phylogeny, at least based on my interpretation of the Materials and methods section. The authors state that they aligned coding DNA sequences in MEGA6 with the "codon option" activated. The codon option will avoid introducing gaps into codons, but it does not correct for systematic biases that can get introduced into protein phylogenies that cover a wide range of divergent species with distinct codon usage biases. For this reason, it is generally accepted that such phylogenies need to be constructed with direct use of the amino acid sequences, as stated in the documentation for MEGA6. If amino acid sequences were used for phylogeny reconstruction (after alignment of nucleotide sequences), then this needs to be clarified in the methods. However, if DNA sequences were used during phylogeny reconstruction, then both phylogenetic analyses presented in Figure 1 need to be repeated to make sure the proposed new classes of receptor are recovered.
2) There are major issues with congruency between the phylogenies presented in Figure 1 with and without the ctenophore glutamate receptors, and this directly impacts the proposed new groups that the authors identify. Specifically, the AK and AKD2 groups in the main figure are split do not hold up in the analysis with the ctenophores included. While AKD2 is not a firmly-proposed new class, AK is, and this is concerning. It certainly suggests that proposing the AK class is an over-reach of the data.
3) The two phylogenies do not use the same names for the majority of the same genes. This needs to be corrected so that the position of a gene within the two phylogenies can be readily discerned. For instance, in the proposed AKD1 group, two genes in the phylogeny with the ctenophores simply contain XP numbers.
4) Posterior probabilities calculated in Bayesian phylogeny reconstructions are not equivalent to Bootstrap support as the authors suggest – this should be corrected. Furthermore, posterior probabilities tend to be much higher than one would get for bootstrap values if the same alignment was used to reconstruct a phylogeny under alternate methods such as Maximum Likelihood. In my experience, while branches with bootstrap support of >60% in Maximum likelihood tend to have relatively stable position in a phylogeny as sequences are added or subtracted, you don't get a similar stability of branch position until a Bayesian posterior probability is over 90%. This is borne out here by the fact that the AK group and the base of the greater AK/Delta/AMPA/Kainate branch, which contain several posterior probabilities below 90%, are not stable between the two phylogenies. The Bayesian approach is a valid one, but here it seems unable to resolve one of the key groups that the authors propose. It isn't a valid approach to simply remove the ctenophore channels to get the desired result. Confidence in the AK group could be increased if it could be recovered by an alternate method of analysis such as Maximum Likelihood, and/or was recovered with high support in phylogenies based on amino acid sequences (point 1).
5) The methods of receptor identification and species selection are reasonable for identifying most of the glutamate receptors in a species or phylogenetic group, but they are by no means exhaustive. Receptors will be missed if automated gene predictions miss key exons or miss genes entirely. While this is not a major issue for the phylogenetic classification of receptors that have been found, it is a major problem for claims of receptor loss or absence in various phylogenetic groups. For example, the lack of NR2 subunits in Nematostella is used to claim a loss of NR2 in Cnidarians. While this might end up being true, the claim requires use of all data available from Nematostella (genome and transcriptomes in addition to gene predictions) and at least several of the other phylogenetically-diverse cnidarians for which data is available. If still nothing is found, then it may be safe to propose absence across this phylogenetic group. In this particular case, it isn't even clear if the absence of NR2 subunits would represent a loss, as the phylogeny could be interpreted to say that NR2 subunits simply evolved after the split between cnidarians and bilaterians. Throughout the paper, claims of loss in various groups should be reviewed to determine if data has been analyzed sufficiently to make the claim. If not, the claim in question could be removed or additional searches could be conducted to support the claim.
This paper presents a comprehensive phylogenetic analysis of metazoan glutamate receptor evolution, resulting in discovery of new families of glutamate receptor ion channels and GPCRs in phyla excluded from previous studies. Using amphioxus as a model system, mRNA expression in native tissue is established for some of these novel receptors via quantitative PCR, and using heterologous expression systems, cell surface expression and functional ion channel activity is established for one of the novel ion channels.
The work is of interest because it gives deeper insight into the evolution of this important class of neurotransmitter receptor than prior studies, which had largely focused on individual phyla, e.g. recent papers on ctenophores, molluscs or insects. Many of these studies had previously reported either receptors with novel properties (e.g. glycine activated homomeric ctenophore iGluRs), and/or the expansion of individual gene families in selected phyla, but did not do so in the context of a comprehensive analysis of metazoan phyla as presented here, and had failed to identify the AK and AKD1 families reported here.
As detailed below, there are issues in the current manuscript which in combination preclude publication until addressed. These include the need to provide a better presentation and interpretation of sequence alignment data; a more complete description of the sequences used for analysis; in addition, there are numerous errors which in places make the data presented very challenging to interpret. All of this can be addressed by careful rewriting.
1) It is unclear why the authors essentially exclude ctenophore iGluRs from their classification. In their revised scheme they identify eight iGluR families expressed in animals, but for a truly comprehensive study of metazoan iGluRs this needs revision to include at least two additional ctenophore subtypes, which are activated by glutamate and glycine respectively. Figure 1 should be replaced by the current Figure 1—figure supplement 4 using appropriate shading.
2) In multiple cases the sequence alignments are very difficult to interpret because they are presented without any shading to indicate identity/conservation and in addition lack coloring to indicate biochemically important residues (e.g. Cys residues known to form disulfide bonds in iGluRs of known structure etc). This should be addressed for ALL figures containing sequence alignments as commonly done by other groups. Notably, Figure 1—figure supplement 6 shows very poor conservation of the highly conserved M3 region highlighted in yellow for the four GluAk_Spu subunits, without any comment noting this.
3) Information should be given in a Table detailing the following for each of the novel AK, AKD1 and AKD2 predicted amino acid sequences: length (documenting unusual insertions or deletions separately for the ATD, LBD and TMD domains) and whether a predicted signal peptide is present and the predicted cut site. A good example can be found in Table S1 of Alberstein et al., 2015. For example, GluAk-beta_Nve appears to lack alpha helix D, a key structural element in the LBD.
4) I am concerned that for both of the novel iGluRs selected for heterologous expression, it was necessary to introduce a vertebrate signal peptide. Are these receptors actually expressed in native tissue? For GluAkd1-eta native expression is likely impossible because the signal peptide is a non consensus sequence that is much too short to form a signal peptide. For GluAkd1-alpha the sequence is much longer than normal, but the Sig3P server does predict cleavage.
5) Subsection "Phylogenetic identification of two unreported iGluR classes". If it is most probable, as stated, that AKD2 is a long-branch attraction artifact does this warrant (i) inclusion of AKD2 as a separate family in Figure 1 and (ii) speculation that AKD2 may represent the ancestral genes of AMPA, kainate and Delta receptors? Revise text to make it clear that the earlier hypotheses that start this section are actually the least probable.
Thank you for resubmitting your work entitled "Evolution of Glutamate Receptors Reveals Three Unreported Classes and Divergent Phyla-Specific Adaptations" for further consideration at eLife. Your revised article has been evaluated by Richard Aldrich (Senior Editor), a Reviewing Editor, and two reviewers.
1) The discussion and classification of what constitutes a new gene family here is somewhat arbitrary and confusing to follow, even though the data is quite clear. In the iGluRs, they have identified the Epsilon and Phi families, a cnidarian-specific clade of NMDA-like receptors and a cnidarian-specific sister group to the AMPA/Kainate/Delta/Phi group of bilaterians. These cnidarian-wide clades are mentioned at certain points in the manuscript, but the authors only headline the Epsilon and Phi families, even though Phi is more restricted in terms of phylogenetic spread than the cnidarian-specific groups. It might be better to say they have found one new type conserved across the metazoa (Epsilon), and 3 phylogenetically-restricted new types in cnidarians and deuterostomes. The authors should not arbitrarily state that their study has increased the number of iGluR classes from 6 to 8 when the data shows they have changed it from 6 to 10. Similarly, for the mGluRs their data shows that they have found: (1) that an mGluR type previously identified as insect-specific is widespread in bilaterians (they have not found a new family as stated, but interesting new information about an existing family), (2) three new cnidarian-specific classes, and (3) a class restricted to basal metazoan lineages including ctenophores, sponges, placozoans and perhaps cnidarians. Support for this latter class is just as strong as for any of the other classes, but it is puzzlingly not presented as a new class in the manuscript. The authors actually state that no classes of mGluRs span multiple phyla, but this ignored class does. The data seems to indicate that the number of mGluR types has moved from 3 to 8, not 3-4.
2) Throughout the manuscript, discussion of when the various types of receptors arose and what phyla/phylogenetic groups lost various receptor types is not interpreted correctly in light of recent advances in our understanding of the metazoan phylogeny. The current view is that ctenophores are the earliest-diverging metazoan phyla, followed by Sponge. Placozoans and Cnidarians group separately from ctenophores and sponges with the bilateria. Within this Placozoan/Cnidarian/Bilaterian group, the most accepted (though not universally accepted) order of divergence is placozoans > cnidarians > bilaterians. With this view of the animal phylogeny in mind: a) There is no evidence to support the claim that NMDA receptors were lost in ctenophores, sponges and placozoans. The simplest interpretation of the data is that NMDA receptors evolved after the divergence of these clades in an ancestor of cnidarians and bilaterians. This needs to be corrected. b) The Epsilon family, to which all ctenophore iGluRs belong, is likely to be the oldest and thus possibly the ancestral metazoan iGluR lineage. Epsilon must have first appeared in an ancestral metazoan prior to the divergence of ctenophores. It would thus be found in an ancestor of all extant metazoans, not just an ancestor of ctenophores/placozoans/cnidarians. Given the animal phylogeny above, that ancestor would also be the ancestor sponges and bilaterians.c) I can find no mention of sponges (porifera) in the results or discussion related to iGluRs. Are iGluRs not found in sponges? Or were they left out of the analysis? If the latter, that is a major oversight that needs to be corrected (and the only experiment I would see needing to be redone). If the former, the authors need to state that none were found, and then they could suggest that the Epsilon family might have been lost from sponges. There is no evidence presented that other iGluR classes were present prior to the divergence of sponges, so the absence of other classes in sponge would not represent gene loss. d) For the mGluRs, the ignored classes with sponge/ctenophore/placozoan/cnidarian sequences could be interpreted as the ancestral mGluR class, given that it is the only class present in the earliest diverging metazoan lineages.
3) The discussion of NMDA receptor evolution is not clear in all places. The data say that an NR2/NR3 ancestor evolved in the cnidarian/bilaterian ancestor. NR2 and NR3 then separated in a bilaterian ancestor after the divergence of bilaterians and cnidarians. In addition, cnidarians have a separate phyla-specific class of NMDA-like subunits. In Figure one, the cnidarian NR2/3 sequence and phyla-specific class are incorrectly lumped together as NMDA-like subunits. This needs to be corrected.
4) Alberstein (2015) reported that all ctenophore iGluRs contain a conserved disulfide bond in loop 1 of the LBD that is found also in NMDA receptors but not in AMPA receptors (their Figure 1, Figure 1—figure supplement 3 and Figure 1—figure supplement 4). This is an unusual structural feature, which combined with activation of some ctenophore iGluRs by glycine lead Alberstein to speculate that these were NMDA receptor precursors. Is this feature found in all of the Epsilon class? If not, is it conserved in sub branches of this clade?
In the first place we want to sincerely thank the reviewers for their constructive positive and comments on our work. For the past four months we have worked very hard to address all of them. We present a revised article that includes extended versions of our previous analysis and a new phylogenetic study. The main conclusions of the article remain essentially unchanged.
1) Phylogenies are now presented using protein sequences (Comment 2.1).
3) Protein sequences from ctenophores have been included in the final trees (Comments 2.1, 2.2 and 3.1).
4) We have extended our sequence homology searches in order to identify sequences that might have been missed previously (see the 'Materials and methods' section in the revised version of the manuscript for a description of the new search strategy). This has resulted in the identification of 9 additional iGluRs and 6 extra mGluRs from the species we presented in the previous version of the manuscript. Two of the newly found iGluRs have allowed us to identify the Delta class in Lophotrochozoans. The rest were from the cnidarian N. vectensis and have been key to identify two cnidarian-specific groups of iGluRs. (Comment 2.5).
5) We have reviewed all claims regarding i/mGluR class loss in a particular phylum. This has been done by incorporating more species from that phylum in the analysis and by performing taxon-specific BLAST searches in the NCBI database. We now have at least two species in each phylum, with the exception of placozoans, for which only T. adhaerens is available. As a result, the new version of the iGluR phylogeny is based on the comparison of sequences from 25 species, as opposed to the 15 used previously (see Author response image 1). In a similar manner, for the mGluR phylogeny we have incorporated 11 additional species. (Author response image 2).
Species used in iGluR phylogeny.
In bold new species used in the present version of the manuscript.
Species used in the mGluR phylogeny.
In bold new species used in the present version of the mansucript.
6) The identification of missed sequences and the use of new species have resulted in an iGluR phylogeny with 233 sequences as opposed to the previous one that presented 171. Similarly, the mGluR phylogeny now includes 123 sequences, while previously it contained 103.
7) We have renamed the AKD1 and AK iGluR classes. AKD1 is now referred to as 'Epsilon' and AK as 'Phi'. We consider that these names are more appropriate since the fourth iGluR class discovered was named after the fourth Greek letter (Delta). We thus though that unreported classes could receive the names of subsequent Greek letters. Gene and protein names from the Epsilon class start with Grie/GluE and those of the Phi class by Grif/GluF.
8) Finally, we have incorporated a new phylogenetic analysis to infer the phylogenetic history of AMPA and Kainate classes in protostomes (Figure 1—figure supplement 2 and Figure 1—figure supplement 3; see Author response image 3 for a list of species used in this phylogeny).
Species used in the phylogeny of AMPA and Kainate classes in protosotomes.
9) Taking into account the three phylogenies we present sequences from 45 different species. In this new version of the manuscript we include, as 'Source Data', the protein alignments used for each phylogeny. These are: Figure 1—source data 4, Figure 1—source data 5 and Figure 4—source data 3.
1) All subunits of ionotropic glutamate receptors from ctenophores belong to the Epsilon class (former AKD1 class).
2) The Phi class (former AK) now incorporates sequences from hemichordates but no longer includes sequences from the cnidarian N. vectensis. Phi proteins are only present in echinoderms, hemichordates and non-vertebrate chordates. For this reason, we hypothesize that this class was formed in an ancestor of deuterostomes. The addition of more cnidarian species to the phylogeny has revealed the existence of a cnidarian-specific group of non-NMDA iGluRs that contains the former AKs from N. vectensis.
3) In addition to this group there is a second cnidarian-specific group of NMDA-like iGluRs.
4) In the previous phylogeny a N. vectensis gene grouped with NMDA3 proteins (GRIN3_A/B_Nve). Yet in the present phylogeny this protein (current name GluN2-3_Nve) appears as the pre-bilaterian orthologue of both NMDA2 and NMDA3 classes. As we have not been able to identify members of the NMDA2 or NMDA3 classes in any of the cnidarian species investigated, the current phylogeny indicates that neither class NMDA2 (as proposed in the previous version of the manuscript) nor class NMDA3 would exist in pre-bilaterians.
5) Proteins from the former AKD2 class no longer cluster together. Therefore, we no longer refer to them as potentially constituting another class.
6) The existence of cnidarian- and placozoan-specific groups of non-NMDA iGluRs indicates that the duplication leading to the AMPA/Kainate/Phi/Delta branch occurred in pre-bilaterians.
7) Changes in the Delta class: we have been able to identify members of the Delta class in a few (three) lophotrochozoan species. Among ecdysozoans, priapulids might have retained members of this class as well. However, our phylogeny indicates that arthropods and nematodes, which account for the vast majority of ecdysozoan species, have lost this class.
8) The current version of S. purpuratus (sea urchin) genome does not include genes coding for members of any of the three NMDA classes. The inclusion of Acanthaster planci (starfish) in the revised phylogeny has allowed to identify members of classes NMDA1 and NMDA3 in the Echinodermata phylum. Class NMDA2 might have been lost in echinoderms, genomic data from other echinoderm species could help clarify this issue.
9) Author response image 4 and Author response image 5 summarize the distribution of iGluR classes among metazoan phyla in the previous version of the manuscript and inthe current version, respectively.
Previous distribution of iGluR classes among metazoan phyla.
Current distribution of iGluR classes among metazoan phyla.
# We could only find 3 species with members of the Delta class in Lophotrochozoans.
* Among Ecdysosoans members of the Delta class are only found in Priapulids. We have not found genes coding for subunits of Delta receptors in arthropod or nematode species.
1) We have identified three cnidarian-specific groups of mGluRs.
2) We now report Class IV proteins in echinoderms.
3) We have found Class I members in a few ecdysozoan species, including priapulids. Previously, we could only identify Class I mGluRs in lophotrochozoans amongst protostomes.
4) Author response image 6 summarizes the present distribution of mGluR classes among metazoan phyla.
Current distribution of mGluR classes among metazoan phyla.
Our detailed responses to all issues raised by the three reviewers follow below.
1.1) In the Introduction the following phrase needs clarification: Thesespecies present a genomic structure, gene sequence and gene contentresembling that of their ancestors…. We do not know the genome of their ancestors and therefore this is an inference.
1.2) Structural models of two proteins from the new AKD1 class from amphioxus where constructed using known iGluR structures. Although this is at first glance reasonable, a less biased approach should be taken. Models should also be produced letting the software choose the appropriate structural template. Also, the MolProbity scores should be given as well as a discussion of the probability that the sequences adjust to the templates.
To address this point models for the full-length proteins from the Epsilon (formerly AKD1) class were generated again with RaptorX (http://raptorx.uchicago.edu/), which automatically selects the appropriate template for modeling. In our hands, RaptorX outperforms all other 3D prediction software when dealing with distant homologs, as indicated also by the results of the latest CASP contests (http://raptorx.uchicago.edu/about/). For modeling the well-conserved ligand-binding domain of GluE1 we used SWISS-MODEL (https://swissmodel.expasy.org/). Here, from a large number of equally appropriate templates, based on sequence similarity alone, we have chosen the atomic-resolution structure (1.26 Å) of the topologically equivalent rat GluA2 domain, which had been refined to an Rfactor of 12.9% (Rfree: 15.7%).
In response to the reviewer's request, we include two pairs of Supplementary Tables (Figure 1—figure supplement 8) containing all relevant MolProbity information for models of full-length GluE1 (formerly GluAkd1_alpha) and GluE7 (GluAkd1_eta), before and after refinement with ModRefiner (https://zhanglab.ccmb.med.umich.edu/ModRefiner/). Another Supplementary Table (Figure 1—figure supplement 9) presents the MolProbity results for the 3D model of the Gly-bound ligand-binding domain (LBD) from GluE1. As can be seen, MolProbity values compare well with those of structures solved at the resolution of the respective templates. The quality of the generated models is underscored by their low P-values (P-value is the likelihood of a predicted model of being worse than the best of a set of randomly-generated models for the studied protein (or domain). Thus, the smaller the P-value, the higher quality the model).
1.3) In the Discussion section, there is no correlation between the number of iGluR genes and the nervous system complexity of an organism. Is it possible that this might reflect the fact that in lower organisms glutamate receptors also have non-synaptic functions, as chemoreceptors?
"Whether all these proteins are expressed at the synapse and act as neurotransmitter receptors is an issue that will require further investigation. While, at least in amphioxus, all iGluRs are expressed in the central nervous system, their presence in other tissues, such as sensory organs, cannot be ruled out. Those receptors showing more divergent sequences, particularly in residues involved in ligand binding, might respond to molecules that are not neurotransmitters. For instance, they could behave as chemoreceptors in sensory organs, as it is the case of antennal receptors found in insects.
1.4) Please include more information on the iGluRs found on A. thaliana and their relation to the newly proposed phylogeny of iGluRs.
We have not identified the Arabidopsis thaliana iGluRs used in the phylogeny. These had been previously reported (please see Chiu et al., 2002 for a phylogenetic characterization of A. thaliana iGluRs). In the revised version of the manuscript we include this reference. We apologize for not having properly referenced these proteins in our previous version of the manuscript. Furthermore, ionotropic glutamate receptors from A. thaliana had been previously used as an outgroup for phylogenetic analysis of C. elegans iGluRs (please see Brockie et al., 2001, in our manuscript). It is for this reason that we choose to use them for the same purpose.
1.5) In the figures that plot mean values, it is not specified what measure of variability is being used (i.e. standard error, variance, etc.).
We have included this information in the legends to all figures presenting mean values (Figure 2, Figure 3 and Figure 4—figure supplement 2B).
2.1) There appears to be a major flaw in the methods for constructing the phylogeny, at least based on my interpretation of the Materials and methods section. The authors state that they aligned coding DNA sequences in MEGA6 with the "codon option" activated. The codon option will avoid introducing gaps into codons, but it does not correct for systematic biases that can get introduced into protein phylogenies that cover a wide range of divergent species with distinct codon usage biases. For this reason, it is generally accepted that such phylogenies need to be constructed with direct use of the amino acid sequences, as stated in the documentation for MEGA6. If amino acid sequences were used for phylogeny reconstruction (after alignment of nucleotide sequences), then this needs to be clarified in the methods. However, if DNA sequences were used during phylogeny reconstruction, then both phylogenetic analyses presented in Figure 1 need to be repeated to make sure the proposed new classes of receptor are recovered.
In the previous version of the manuscript we presented trees constructed with DNA sequences. This has been addressed and now we only present phylogenies constructed with protein sequences. The 'Materials and methods' section has been re-written accordingly in the revised version of the manuscript (please see subsection "Phylogenetic analyses").
2.2) There are major issues with congruency between the phylogenies presented in Figure 1 with and without the ctenophore glutamate receptors, and this directly impacts the proposed new groups that the authors identify. Specifically, the AK and AKD2 groups in the main figure are split do not hold up in the analysis with the ctenophores included. While AKD2 is not a firmly-proposed new class, AK is, and this is concerning. It certainly suggests that proposing the AK class is an over-reach of the data.
We now include sequences from ctenophores in the final phylogenies of iGluRs and mGluRs. Furthermore, for each phylogenetic analysis we present trees constructed with a Bayesian and a maximum likelihood method (previously we only used the Bayesian approach). We have not found any relevant congruency issues between the two methods used in any of the three phylogenies presented (please see Figure 1 and Figure 1—figure supplement 1; Figure 1—figure supplement 2 and Figure 1—figure supplement 3; and Figure 4 and Figure 4—figure supplement 1).
As just mentioned, the revised iGluR phylogeny now includes sequences from the ctenophore M. leidyi (used previously). However we also include: (i) another member of the same phylum (P. bachei), (ii) three other cnidarians (O. faveolata, A. digitifera and H. magnipapillata), (iii) a second hemichordate (P. flava), (iv) a second echinoderm (A. planci), (v) a third cephalochordate (B. floridae), (vi) a member of the ecdysozoan lineage of priapulids (P. caudatus) and (vii) two other lophotrochozoans: A. californica and C. gigas.
The inclusion of ctenophores, as well as all these other species, in the phylogeny of iGluRs has resulted in an unexpected finding, namely, that all ctenophore iGluRs would in fact belong to the Epsilon class. In this revised version of the phylogeny we have also identified two cnidarian-specific groups of iGluRs.
Finally, we no longer present the AKD2 group as a potential new class.
2.3) The two phylogenies do not use the same names for the majority of the same genes. This needs to be corrected so that the position of a gene within the two phylogenies can be readily discerned. For instance, in the proposed AKD1 group, two genes in the phylogeny with the ctenophores simply contain XP numbers.
This issue has been addressed. Bayesian and maximum likelihood phylogenies present the same protein names in the revised version of the manuscript.
2.4) Posterior probabilities calculated in Bayesian phylogeny reconstructions are not equivalent to Bootstrap support as the authors suggest – this should be corrected. Furthermore, posterior probabilities tend to be much higher than one would get for bootstrap values if the same alignment was used to reconstruct a phylogeny under alternate methods such as Maximum Likelihood. In my experience, while branches with bootstrap support of >60% in Maximum likelihood tend to have relatively stable position in a phylogeny as sequences are added or subtracted, you don't get a similar stability of branch position until a Bayesian posterior probability is over 90%. This is borne out here by the fact that the AK group and the base of the greater AK/Delta/AMPA/Kainate branch, which contain several posterior probabilities below 90%, are not stable between the two phylogenies. The Bayesian approach is a valid one, but here it seems unable to resolve one of the key groups that the authors propose. It isn't a valid approach to simply remove the ctenophore channels to get the desired result. Confidence in the AK group could be increased if it could be recovered by an alternate method of analysis such as Maximum Likelihood, and/or was recovered with high support in phylogenies based on amino acid sequences (point 1).
In the previous version of our manuscript we inadvertently referred to posterior probabilities as bootstraps. This has now been corrected (Please see figure legends of phylogenetic trees).
For all phylogenetic analysis we now present two trees, one using a Bayesian approach and a second based on a Maximum likelihood method. In the two iGluR trees construed with these two different methods previously unreported classes are recovered and they have the same position in both trees. Including class AK, currently refer to as class Phi.
2.5) The methods of receptor identification and species selection are reasonable for identifying most of the glutamate receptors in a species or phylogenetic group, but they are by no means exhaustive. Receptors will be missed if automated gene predictions miss key exons or miss genes entirely. While this is not a major issue for the phylogenetic classification of receptors that have been found, it is a major problem for claims of receptor loss or absence in various phylogenetic groups. For example, the lack of NR2 subunits in Nematostella is used to claim a loss of NR2 in Cnidarians. While this might end up being true, the claim requires use of all data available from Nematostella (genome and transcriptomes in addition to gene predictions) and at least several of the other phylogenetically-diverse cnidarians for which data is available. If still nothing is found, then it may be safe to propose absence across this phylogenetic group. In this particular case, it isn't even clear if the absence of NR2 subunits would represent a loss, as the phylogeny could be interpreted to say that NR2 subunits simply evolved after the split between cnidarians and bilaterians. Throughout the paper, claims of loss in various groups should be reviewed to determine if data has been analyzed sufficiently to make the claim. If not, the claim in question could be removed or additional searches could be conducted to support the claim.
We have also included more cnidarian species in our analyses to account for possible bias caused by using only Nematostella vectensis (sea anemone). These are: Acropora digitifera (stony coral), Hydra magnipapillata (hydrozoa), Orbicella faveolata (stony coral), Exaiptasia pallida (sea anemone) and Stylophora pistillata (stony coral).
We have reviewed all claims regarding class loss in a particular phylum. This has been done by including more species from that phylum in the analysis and by performing taxon-specific Blast searches in the NCBI database. In the previous version we used 15 species for iGluR phylogenetic reconstruction and 18 for the mGluR tree. Now we use 25 and 29 species, respectively (please see Figure 1—source data 2 and Figure 4—source data 2, for a list of all species used in i/mGluR phylogenies).
3.1) It is unclear why the authors essentially exclude ctenophore iGluRs from their classification. In their revised scheme they identify eight iGluR families expressed in animals, but for a truly comprehensive study of metazoan iGluRs this needs revision to include at least two additional ctenophore subtypes, which are activated by glutamate and glycine respectively. Figure 1 should be replaced by the current Figure 1—figure supplement 4 using appropriate shading.
We now include two ctenophora species in the final phylogenies of iGluRs and mGluRs. To better account for this phylum we include sequences from Pleurobrachia bachei, as before we only used sequences from M. leidyi. Three of the previously reported proteins from P. bachei (Alberstein et al., 2015) did not present a sufficiently conserved SYTANLAAF motif in our alignments (please see Materials and methods section). For this reason, they were not included in the final tree, these are: PbiGluR5, PbiGluR11 and PbiGluR12. As we now incorporate most iGluRs from these two ctenophore species into the phylogeny we do include members of two subtypes of ctenophore iGluRs (activated by glutamate and glycine) as requested by the reviewer.
3.2) In multiple cases the sequence alignments are very difficult to interpret because they are presented without any shading to indicate identity/conservation and in addition lack coloring to indicate biochemically important residues (e.g. Cys residues known to form disulfide bonds in iGluRs of known structure etc). This should be addressed for ALL figures containing sequence alignments as commonly done by other groups. Notably, Figure 1—figure supplement 6 shows very poor conservation of the highly conserved M3 region highlighted in yellow for the four GluAk_Spu subunits, without any comment noting this.
We thank the reviewer for this comment. We have amended all figures with sequence alignments following the reviewer indications.
3.3) Information should be given in a Table detailing the following for each of the novel AK, AKD1 and AKD2 predicted amino acid sequences: length (documenting unusual insertions or deletions separately for the ATD, LBD and TMD domains) and whether a predicted signal peptide is present and the predicted cut site. A good example can be found in Table S1 of Alberstein et al., 2015. For example, GluAk-beta_Nve appears to lack alpha helix D, a key structural element in the LBD.
We now include Tables detailing all this information, as suggested by the reviewer (please see Figure 1—source data 1, Figure 4—source data 1). When there is a difference in domain length between mammalian proteins and non-vertebrate i/mGluRs this is, in most cases, because non-vertebrate proteins are shorter. We think that these 'potential deletions' are in fact due to poor genome sequence data, as they are usually very large, especially in the ATD domain. For instance, the sequence of GluAKbeta_Nve (now GluCni2_Nve), that the reviewer refers to, only presents 140 LBD residues (while mammalian ones present >250).
3.4) I am concerned that for both of the novel iGluRs selected for heterologous expression, it was necessary to introduce a vertebrate signal peptide. Are these receptors actually expressed in native tissue? For GluAkd1-eta native expression is likely impossible because the signal peptide is a non consensus sequence that is much too short to form a signal peptide. For GluAkd1-alpha the sequence is much longer than normal, but the Sig3P server does predict cleavage.
We provide experimental evidence for the expression of the genes Grie1 and Grie7 in the cephalochordate B. lanceolatum (Figure 2A); these genes code for GluE1 (formerly GluAkd1_α) and GluE7 (GluAkd1_η), respectively. Importantly, we show that these two genes present an enriched expression in the nervous cord of amphioxus, indicating that their main biological function is performed in this tissue. Furthermore, the expression level of Grie1 is similar to that observed for genes coding for subunits of AMPA, Kainate or NMDA receptors from amphioxus (Figure 2A).
Secondly, as indicated by the reviewer, a signal peptide (SP) can be identified in GluE1, although this is longer than canonical SPs found in vertebrates. Nevertheless, the existence of long SPs, of over 50 residues, is well accepted (Nielsen, 2017) and therefore SignalP includes an option to predict them. A simple explanation for the lack of a SP in GluE7, as predicted by SignalP4.1, might be that the 5' end of this gene has not been completely sequenced, missing its SP.
Nevertheless, we have also considered the possibility that SignalP4.1 might be less efficient at identifying SPs from non-vertebrates. We have consulted Dr. Henrik Nielsen, the key scientist responsible for developing this software, who has communicated us that invertebrate sequences are less represented in the training sets used to develop SignalP. He has also explained us that SignalP 4.1 is less efficient at identifying long SPs than SPs with canonical lengths.
This analysis suggests that SignalP4.1 is not as efficient at identifying SPs in non-vertebrate iGluRs as it is for vertebrate ones. The increase in SPs predicted when long SPs are allowed suggests that invertebrate species would have a higher proportion of iGluRs with long SPs. Nevertheless, differences in genomic data quality between vertebrate and invertebrate species might explain part of the differences observed.
In our mRNA expression analysis (Figure 2A) all iGluRs from amphioxus were found expressed in the nervous cord, yet SignalP4.1 fails to identify a SP in approximately half of them. Importantly, SignalP4.1 not only fails to identify SPs in members of the new classes reported in this article, it also misses SPs in proteins from known classes, such as in Kainate or NMDA2 (see Author response image 8).
Altogether these observations suggest that SignalP4.1 is less efficient at predicting SPs in invertebrate species, which in turn might present more iGluRs with long SPs. This is in accordance with the lower number of sequences from non-vertebrate species available for the training sets used to develop SignalP, as reported by Dr, Nielsen.
Therefore, the fact that SignalP4.1 does not predict a SP in an invertebrate iGluR should not be taken as a very strong indication for the absence of a functional SP. Unfortunately, to the best of our knowledge there are no comprehensive studies on the evolution/conservation of SPs along the metazoan kingdom, and the evolution of the machinery involved in SP processing has neither been systematically investigated. Thus, we do not know if all SPs from non-vertebrates would be recognized by mammalian cells (in our particular case human HEK293 cells). It is plausible that the mammalian machinery involved in SP identification and protein trafficking to the plasma membrane is not able to properly recognize SPs from cephalochordate proteins, which diverged from the main metazoan lineage over 650 million years ago (TimeTree: http://www.timetree.org).
Finally, we also want to indicate that GluE1 and GluE7 both expressed well in HEK293 cells without the SP from rat GluA2, although they were not efficiently trafficked to the plasma membrane (please see Author response image 9, and Figure 2—figure supplement 1C,D for immunoblots of whole HEK293 cell lysates 48 hours post-transfection of wild-type GluE1 and GluE7). Note that protein degradation is not observed, which generally is a good indicator of protein stability. In our opinion, the fact that these two proteins, despite their considerable degree of divergence from mammalian homologues (the average protein identity between members of the Epsilon class from amphioxus and mouse iGluRs is around 25%) can be efficiently expressed in HEK293 cells suggests that they are functional, that is, that they are expressed in native tissue. In this regard, we wish to note that non-expressed genes rapidly accumulate deleterious mutations becoming pseudogenes that are usually not identified by standard genomics projects.
Expression of GluE1 and GluE7 in HEK293 48 hours after transfection.
We therefore consider that the fact that wild-type GluE1 and GluE7 could not be trafficked to the plasma membrane of HEK293 cells nor their lack of canonical SPs (as predicted by SignalP4.1) warrant sufficient evidence to discard their expression in native tissue. On the contrary, we think that the following evidence strongly supports their functional expression in native tissue: (i) we do observe expression of both genes at the mRNA level, (ii) GluE1 has a SP, although a long one, (iii) GluE1 presents expression levels similar to those of GluAs, GluKs or GluNs from amphioxus, (iv) both proteins express well in HEK293 cells and, finally, (v) GluE1 forms functional ligand-gated ion channels.
3.5) Subsection "Phylogenetic identification of two unreported iGluR classes". If it is most probable, as stated, that AKD2 is a long-branch attraction artifact does this warrant (i) inclusion of AKD2 as a separate family in Figure 1 and (ii) speculation that AKD2 may represent the ancestral genes of AMPA, kainate and Delta receptors? Revise text to make it clear that the earlier hypotheses that start this section are actually the least probable.
In the revised iGluR phylogeny we no longer consider the AKD2 class.
We regret that our classification was confusing to follow. In order to resolve this problem we now present a different classification system that we think is clearer as well as more accurate. As requested by the reviewers we have included in our lists of 'Classes' phylogenetic groups found in only one phylum, which in the previous version of the manuscript were classified as 'Groups'. This modification together with the inclusion of porifers in the phylogenies has resulted in the identification of 10 iGluR classes and 8 mGluR classes.
In iGluRs the classification system presents 3 levels: family, subfamily and class. In mGluRs only 2, family and class. Although there is not a total agreement in the field of neurosciences, we believe that the term 'family' is generally used to refer to all iGluRs or all mGluRs. Thus, we reserve the use of 'family' for this porpoise; there are therefore only 2 families of GluRs (iGluRs and mGluRs).
We now give more strength to the idea that metazoan iGluRs are organized into 4 monophyletic groups. Three of these groups were already present with the same topology in our previous version of this manuscript, although we did not discuss them as such. The incorporation of porifers has resulted in the addition of a fourth main group to the previous three. We have given the name of 'subfamily' to each of these groups, to distinguish them from 'classes'. Classes are contained within subfamilies. Due to the topology of both iGluR trees (Bayes and ML) we propose that these 4 subfamilies arose from 3 gene duplication events that occurred prior to the radiation of all present-day metazoan phyla (please see Figure 1).
Two of these subfamilies include the well-known vertebrate classes. The first one, subfamily 'NMDA', includes classes NMDA1, NMDA2, NMDA3, NMDA2/3 and NMDA-Cnidaria. The other subfamily, which we have termed 'AKDF', has 5 classes: AMPA, Kainate, Delta, Phi and AKDF-Oca. The later exclusively contains sequences from one sponge (O. carmela), but it is present in the Bayes and ML phylogenies with very good statistical support. The same is true for class NMDA2/3, for which we have only identified 1 sequence from N. vectensis. Following reviewers' advice, we now refer to these groups as 'classes'.
All these classes (but AKDF-Oca) were already present in the previous version of our manuscript containing the same proteins and presenting the same topology. Beyond AKDF-Oca we have found several proteins from porifer, placozoan and cnidarian species that clearly belong to the AKDF subfamily. There are two important considerations about these proteins. In the first place one of them (from T. adhaerens, in green) is in the branch of Deltas in the Bayes tree. Yet the ML tree does not place this sequence in this class. Therefore, we do not consider it as being a Delta member. In the second place, the branches involving these sequences present different topologies in the Bayes and ML trees. For this reason, we have left them unclassified.
The former Epsilon class is now referred to as Epsilon subfamily, as it is at the same phylogenetic level as the other subfamilies. Finally, as mentioned earlier, the addition of sequences from sponges has allowed identifying a fourth subfamily, being located at the very base of the tree, that we have termed 'Lambda'.
Finally, regarding the phylogeny of mGluRs, the addition of new poriferan species has introduced some changes regarding the organization of non-bilaterian classes. Please see below.
2) Throughout the manuscript, discussion of when the various types of receptors arose and what phyla/phylogenetic groups lost various receptor types is not interpreted correctly in light of recent advances in our understanding of the metazoan phylogeny. The current view is that Ctenophores are the earliest-diverging metazoan phyla, followed by Sponge. Placozoans and Cnidarians group separately from ctenophores and sponges with the bilateria. Within this Placozoan/Cnidarian/Bilaterian group, the most accepted (though not universally accepted) order of divergence is placozoans > cnidarians > bilaterians.
We are aware of the current debate in the field about which metazoan phylum is the earliest-diverging one. Several important publications (Moroz et al., 2014; Ryan et al., 2013 or Dunn et al., 2008, among others) advocate for ctenophores as the earliest metazoan phylum; yet, other authors support a Porifera-First model (please see, for instance, Simon et al., 2017 and Feuda et al., 2017). We believe that on the matter of which is the earliest-diverging metazoan phylum, the jury is still out, and that unanimous consensus has not yet been reached. Nevertheless, in order to account for the recommendations of the reviewers, we have interpreted our phylogenies considering a ctenophores-first model.
With this view of the animal phylogeny in mind: a) There is no evidence to support the claim that NMDA receptors were lost in ctenophores, sponges and placozoans. The simplest interpretation of the data is that NMDA receptors evolved after the divergence of these clades in an ancestor of cnidarians and bilaterians. This needs to be corrected.
We disagree with this remark. We believe that since the branch of the NMDA subfamily appears at the same level that the branch with ctenophore iGluRs (See Figure 1), NMDA protein/s must have appeared before the divergence of Ctenophores and lost in this phylum.
b) The Epsilon family, to which all ctenophore iGluRs belong, is likely to be the oldest and thus possibly the ancestral metazoan iGluR lineage. Epsilon must have first appeared in an ancestral metazoan prior to the divergence of ctenophores. It would thus be found in an ancestor of all extant metazoans, not just an ancestor of ctenophores/placozoans/cnidarians. Given the animal phylogeny above, that ancestor would also be the ancestor sponges and bilaterians.
While we realize that the possibility presented by the reviewers is very parsimonious and would agree well with the ctenophores-first model of animal evolution, the present iGluR phylogeny indicates that there were 4 subfamilies of iGluRs in the last common ancestor to present-day metazoans: Epsilon, NMDA, AKDF and Lambda.
c) I can find no mention of sponges (porifera) in the results or discussion related to iGluRs. Are iGluRs not found in sponges? Or were they left out of the analysis? If the latter, that is a major oversight that needs to be corrected (and the only experiment I would see needing to be redone). If the former, the authors need to state that none were found, and then they could suggest that the Epsilon family might have been lost from sponges. There is no evidence presented that other iGluR classes were present prior to the divergence of sponges, so the absence of other classes in sponge would not represent gene loss.
Thanks to the identification of the Compagen Database (http://www.compagen.org), which we previously were not aware of, we could incorporate porifers into the iGluR phylogeny. We have also included poriferan species from Compagen into the mGluRs phylogeny.
d) For the mGluRs, the ignored classes with sponge/ctenophore/placozoan/cnidarian sequences could be interpreted as the ancestral mGluR class, given that it is the only class present in the earliest diverging metazoan lineages.
We agree that in our previous manuscript we did not discuss in detail the large branch close to the base of the tree including sponge, ctenophore, placozoan and cnidarian sequences. Although, as stated by the reviewers, this branch was consistently found in the Bayes and ML trees with good statistical support. The reason why we did not discuss it in detail is that we were puzzled by the presence of a more basal branch that only included cnidarian sequences (Cni3). We found difficult to integrate these two groups into a parsimonious explanation of the early evolution of mGluRs.
In any case, the incorporation of new porifers in the mGluR phylogeny has changed this part of the tree. Now the former Cni3 group has merged with another former cnidarian-specific group (Cni2) and is located next to the former Cni1 group. These two branches of cnidarian mGluRs are clearly orthologue to classes I-IV in the present phylogeny. Similarly, we found clear orthologues of these 4 classes in placozoans and sponges. Altogether these non-bilaterian orthologues of classes I-IV are organized in 4 branches, with identical topology in bates and ML trees (Figure 6 and Figure 6—figure supplement 1). Following reviewers' advice, we now refer to these branches as classes. Thus, in total, we report 8 classes of mGluRs.
As previously mention we have re-formulated the overall classification of iGluRs. We hope that the explanation regarding NMDA classes evolution is now clear. We have also amended the error in Figure 1.
This is certainly a relevant point as it could have important evolutionary implications. It is something we looked into but did not include in the manuscript, as these two cysteines are clearly absent from all other Epsilon proteins. Furthermore, these two residues are also absent from all non-bilaterian AKDF proteins and most Lambda members. In our alignment (Figure 1—source data 3) the second (most C-term) cysteine forming this bond is very well aligned in all NMDAs and Epsilons from ctenophores. The only sequence among unreported groups presenting a good alignment of this cysteine is GluL1_Oca (a Lambda member from O. Carmela). This sequence also presents a cysteine 9 residues upstream, which could participate in a S-S bond. Nevertheless, experimental data will probably be needed to confirm if this bond is formed.
We thank I Gich (Biomedical Research Institute Sant Pau) for support in biostatistics and Dr. Xavier Daura and Oscar Conchillo (Computational Biology Group and Data Center from the Institut de Biotecnologia i Biomedicina (UAB). DRV, GG, RRV, JL, PFP, and AB work at a research institute supported by the CERCA Programme/Generalitat de Catalunya.
© 2018, Ramos-Vicente et al.