All gene sequences used in this study were manually curated by us and published in Chipman et al. (2014). In brief, we ran a number of similarity searches on the S. maritima genome to identify candidate S. maritima GR and IR genes. We then superimposed the results of these searches, available EST information, and the automated gene predictions on the genome to manually annotate each chemosensory gene and obtain the dataset we used in the present study (see details in Chipman et al., 2014; S. maritima IR and GR CDS sequences are included in the Supplementary Material).

The methodological approach we used herein to detect signs of positive selection depends on sequence variation and its accuracy is reduced if it is applied to very divergent paralogs. To circumvent this problem, we first identified recent clades based on synonymous divergence in both the GR and the IR gene families within which separate analyses would be carried out (Figure 1). For that, we first obtained a multiple alignment of amino acid sequences with MAFFT v. 7 (Katoh and Standley, 2013) using default parameters and used it in a maximum likelihood search to obtain a gene tree using the program RAxML v.8 (Stamatakis, 2006) with the PROTGAMMAWAG model. We were aware that the S. maritima GRs and IRs most likely do not represent species-specific chemosensory expansions since we have low phylogenetic coverage in this part of the arthropod tree (the most recent common ancestor of centipede and insects from which we have IR sequences was ~700 mya); chemoreceptors from other myriapod species are expected to spread out across the S. maritima clades. Nevertheless, some of the considered sub-clades (Figure 1) are relatively recent and they could in fact represent recent paralogous expansions. Both gene families showed expanded clades with more than five recently diverged paralogs (Figure 1). In the GR family, six main clades were identified, while in the IR family there were four clades. One IR clade (Clade 1, Figure 1A) was further subdivided into two subclades to maximize the number of analyzed positions and the power of the ML analyses. Notice that these clades only include the paralogous genes belonging to the centipede specific expansion. We then obtained separate multiple alignments of the nucleotide sequences of the genes included in each clade with the program MAFFT (Katoh and Standley, 2013) for the statistical analyses. Pseudogenes, i.e., fragments with early stop codons and/or frameshift mutations, were excluded from the analysis since they have inflated mutation rates that are not adaptive.

We used the codon-based substitution models (Nielsen and Yang, 1998; Yang et al., 2000) implemented in codeml program (PAML 4.4; Yang, 2007) to estimate non-synonymous to synonymous substitution rate (d_N/d_S) ratios (ω) across paralog sequences. To determine whether there is significant evidence of evolution by positive selection in recently diverged paralogs, we used codeml to estimate the goodness of fit of different models to the observed data using a maximum likelihood approach. First, we fit to the data a model with one single average ω (Model M0) in order to obtain the branch lengths to be used as initial values for more complex models. Then we estimated model parameters and the log likelihood (L) of the data under two alternative models: M7, a model with ten classes of sites with beta distributed ω-values in the interval (0–1) (no positively selected sites allowed) and M8, where an extra class of sites (with proportion p_s) constrained to have ω_s≥1 is added to the beta model. We compared the L of these two models by using the likelihood ratio test (LRT; α = 0.05 after controlling for the false discovery rate—FDR; Benjamini and Hochberg, 1995). As pointed out by Swanson et al. (2003) the M7 vs. M8 comparison may result (in some particular cases) in a high proportion of false positives (significant tests in the absence of positive selection). In fact, the test only indicates if there is a class of sites with a ω> 1 but not if this ratio is significantly greater than 1. Hence, we also used the M8a model, which was proposed as an alternative null hypothesis in which the extra class in M8 is fixed to ω = 1, making a more refined (and conservative) test for the existence of positively selected sites (when compared to model M8). We assumed that the asymptotic null distribution of the LRT statistic is a c² (df = 2) for the M8 vs. M7 comparison and a 50:50 mixture of point mass 0 and c² (df = 1) for M8 vs. M8a comparison (Self and Liang, 1987; Swanson et al., 2003). All analyses were repeated using different ω starting values to avoid local optima in the maximum likelihood calculation. The Bayes Empirical Bayes (BEB; Yang et al., 2005) analysis under model M8 was used to identify codons under positive selection. Briefly, this method calculates, for each site in the alignment, the posterior probability (PP) of belonging to each of the different omega site classes defined in the model (11 classes in M8). Sites with PP > 0.5 of belonging to the class of sites with ω> 1 are candidates to have been under positive selection.

To evaluate whether sites with a probability of being under positive selection have a random or patterned distribution across the different protein domains, we mapped the amino acids corresponding to these codons in the predicted functional domains of the IR and GR receptors. We used different pieces of information to delimit these domains. For IRs, we first predicted the transmembrane segments in the multiple sequence alignment of each clade using TMHMM Server v. 2.0 (http://www.cbs.dtu.dk/services/TMHMM/) with default settings. Then we delimited the pore loop, S1 and S2 ligand binding domains of S. maritima IRs using the alignments of the amino acid sequences of Drosophila and human iGluRs and Drosophila IRs in Benton et al. (2009) as a guide (based on the sequence similarity with some highly conserved residues).

For GRs, the prediction of the transmembrane domains is less clear-cut because not all of the domains are as prominently hydrophobic as most transmembrane proteins such as GPCRs, channels, and transporters, and various transmembrane prediction programs regularly under-predict them and sometimes over-predict them; they consequently also often indicate the incorrect orientation in the membrane (e.g., Benton et al., 2006; Smart et al., 2008; Zhang et al., 2011; Hull et al., 2012). The orientation in the membrane is known for the conserved Odorant Receptor Co-Receptor or Orco protein and other ORs (Benton et al., 2006; Lundin et al., 2007; Smart et al., 2008; Tsitoura et al., 2010; Hull et al., 2012), which are clearly related to the GR family in the insect chemoreceptor superfamily (Robertson et al., 2003), as well as for one GR (Zhang et al., 2011). It is the opposite of the GPCR orientation, with the N-terminus internal to the cell. The transmembrane domains were therefore defined on the basis of a combination of hydrophobic regions highlighted in CLUSTALX (Larkin et al., 2007) alignments and Kyte-Doolittle hydropathy plots (Kyte and Doolittle, 1982).

We tested whether amino acid sites with PP > 0.5 of being one of the positive selected sites in the BEB analysis were homogeneously distributed across functional domains. For that, we built a contingency table (2 × n table, where n is the number of domains defined in each receptor) to compare the relative number of sites with PP > 0.5 with the number of sites with PP < 0.5 across domains. In this way, the total number of sites of each domain is implicitly taken into account, removing the effect of domain length. We used the Fisher's exact test and FDR (False Discovery Rate—when doing multiple comparisons across clades) to obtain P-values.

The number of genes per clade varied from 5 to 20 in the IRs and the total tree length in number of substitutions per codon varied from 2.59 to 13.24 (Table 1). These medium-to-high levels of nucleotide divergence have been found to maximize power and accuracy of the LRT (Anisimova et al., 2001). The ω-values across alignment (codeml model M0) never reached 1, being largest (ω = 0.579) for the IR Clade 1A (Table 1). The lowest value (ω = 0.221) was observed in Clade 2, which was also the oldest clade as inferred from the longer branch lengths separating its members (Figure 1A). Among the GRs, the number of genes per analyzed clade varied from 7 to 14 and the maximum tree length was 7.99 in GR Clade 6 (Table 1). The ω-values were in a similar range (0.337–0.646) to that observed in the IR family, although a little higher, on average (Table 1). Interestingly, this is in accordance with the GR family having on average the lowest selective constraint among chemosensory families observed in Drosophila (Sánchez-Gracia et al., 2009).

All the estimated ω-values are relatively high as compared to values found in ortholog comparisons across several gene families in Drosophila (Clark et al., 2007). Chemosensory genes, however, may have extraordinarily low functional constraints, with high ω-values even in among-ortholog comparisons as compared to the ω-values obtained for other Drosophila gene families (Sánchez-Gracia et al., 2009, 2011). In some Drosophila chemosensory ortholog comparisons, ω-values were as high as the ones observed here among paralogs; these were outliers, however, as gene family averages of ω estimates from ortholog comparisons were never higher than 0.25 (Sánchez-Gracia et al., 2011). The relatively high ω-values we observed among S. maritima IR and GR paralogs is in accordance with a relaxation of purifying selection as expected in duplicated genes (Innan and Kondrashov, 2010) and observed in chemosensory gene families in other species (Smadja et al., 2009; Almeida et al., 2014). On the other hand, ω-values below one suggest that although selective constraints are relaxed and probably some codons might be under positive selection, most codons are evolving under purifying selection in these genes.

The presence of positively selected sites accounting for the increased ω-values observed in the S. maritima paralogs comparisons was statistically supported in IR Clade 1 and IR Clade 4, and all six GR clades (Table 2). Among the IR clades, the strongest signal of positive selection was in the larger and more recently expanded clade of centipede IRs (Clade 1, LR P-values in Table 2), both considering the clade as whole and separately for each of the two main sub-clades (Clade 1A and 1B, Figure 1A). The range of ω-values of the positively selected sites in IRs (1.2–2.7), estimated under the M8 model, is slightly smaller but similar to those estimated for GRs (1.4–5.5). The proportion of sites estimated to have an ω> 1 ranged from p₁ = 4.6–26.7% in IRs and p₁ = 2.1–20.8% in GRs, which indicates high levels of positive selection acting during the diversification of these receptors. Evidence of positive selection in recently duplicated chemosensory genes has also been found in Drosophila using similar site-specific methods (Croset et al., 2010) or using a branch-model approach (Almeida et al., 2014). The application of the latter methods demonstrated that the evidence of positive selection is stronger in the divergence among recent paralogs than in orthologous diversification. In the present work we applied for the first time a model accounting for heterogeneous selective pressure across sites to study the diversification of chemoreceptor paralogs in a non-insect species.

Given the highly significant LRT attesting for the presence of positively selected sites in the analyzed paralogs, we used the BEB analysis to predict the putative location of such sites (Supplementary Tables S1–S10). The BEB analysis, however, had a low performance in pinpointing the specific codon sites affected (i.e., positions with strong posterior probability (PP) to belong to the ω > 1 site class). The number of sites with PP > 0.95 in IRs was very low (e.g., 3 out of 21 sites with probability to be under positive selection in IR proteins of Clade 1), with no site with PP > 0.99 (the sites with PP > 0.95 and PP > 0.7 belonging to the positively selected class were 9 and 22, respectively). In the GR family, the number of sites with high probability is somewhat higher; the GR Clade 4 had the highest number of sites with PP > 0.95 (8), followed by GR Clade 6 (7).

To investigate how selective pressure is distributed across functional elements, we mapped the location of all amino acid sites under positive selection (PP > 0.5) in the predicted functional domains of the IR and GR proteins. The structure of the IRs is characterized by an N-terminal domain (N-term), which is highly variable in size, a bipartite ligand-binding domain (LDB), three transmembrane domains (M1–M3), an ion channel domain composed of a pore loop (P-loop) located between M1 and M2 and a short C-terminal domain (C-term) (Croset et al., 2010; Traynelis et al., 2010). The LDB recognizes specific ligands that trigger the opening of the ion channel pore (Armstrong et al., 1998). The structure of the GRs consists of seven transmembrane domains (M1-7), three extracellular short loops (ECL1-3), three intracellular long loops (ICL1-3), an extracellular C-terminus (C-term), and an intracellular N-terminus (N-term) (Clyne et al., 2000). The GRs are evolutionarily related to the ORs, sharing the same basic structure. Although little is known about the specific function of the GRs' domains, some information has been gathered on the ORs' domains that could be extrapolated. For instance, several studies have shown that the outer edge of different transmembrane domains affect ligand specificity (Nichols and Luetje, 2010; Pellegrino et al., 2011; Leary et al., 2012; Hughes et al., 2014).

In a lineage specific diversification of chemosensory genes, such as the observed expansions of IRs and GRs in centipedes, it would be expected that the main domains under positive selection would be the ones affecting ligand recognition/specificity. Positive selection for diversification in these domains would lead to an increase in the number of chemicals able to be sensed by the system. In agreement with this expectation, using a branch-site model, Croset et al. (2010) identified a positively selected site (PP > 0.95) in the LBD of a recently duplicated IR of Drosophila mojavensis. Here, we found that the vast majority (87%) of the IR residues predicted to be under positive selection (including all sites predicted with PP > 0.95) are in the extracellular loops (Figures 2A–D, Supplementary Tables S1–S4). In the proteins of IR Clade 1, most of these residues were located in the S1 and S2 domains, while those of the IR Clade 4 preferentially accumulated in the large N-terminal domain. We also inferred a number of positively selected sites in the P-loop domain, suggesting that such structure might also have an important role in the functional diversification of IR paralogs. This pattern is specific to this kind of ion channel receptor. Importantly, the accumulation of positively selected sites in specific domains (in this case the LBD and P-loop domains) is significant and not explained by the relative length (in number of amino acids) of these domains within the protein (P = 1.27×10⁻¹⁵ or P = 6.28×10⁻¹² taking together the sites analyzed for proteins of IR Clade 1 and IR Clade 4, or IR Clade 1A, IR Clade 1B and IR Clade 4, respectively), according to the results of the Fisher's exact test. Considering each clade separately, we found a significant departure from a homogeneous distribution of candidate sites across domains in all cases (P<0.006, after controlling for FDR), except for IR Clade1A (P = 0.424). The LBD and the P-loop, therefore, have made a major contribution to the functional diversification of IR paralogs in centipedes.

Similarly, in the GR family the majority of the candidate sites are located in outer parts of the proteins (Figures 3, 4). Overall, the ECL1-3 and the outer sections of the transmembrane domains accumulate more candidate sites than expected given their relative length within the protein (P = 0.021). This is very evident for GR Clades 1, 3, 5, and 6 (P<0.02 after controlling for FDR; Figures 3B,D–F, Supplementary Tables S5–S10). In these proteins the sites with PP > 0.95 also grouped in the ECL1-3 and transmembrane domains (P = 1.61×10⁻¹²; Figure 4), indicating an important functional role of these amino acids. Although there is no functional study focused on GRs, several authors have identified specific amino acids that are involved in ligand specificity of particular ORs. For instance, Nichols and Luetje (2010) found that the outer edge of TM3 of DmOr85b affects ligand specificity; Pellegrino et al. (2011) found that a Val91Ala polymorphism at the outer edge of TM2 of DmOr59b affects ligand-specificity; Leary et al. (2012) found that an Ala148Thr substitution at the outer edge of TM3 in a moth pheromone receptor mediated the ability to detect a new pheromone component; and Hughes et al. (2014) found that mutation of Ala195 on the outer edge of TM4 in Anopheles gambiae Or15 greatly affects ligand specificity. Hence, we expect that the positively selected sites we found in GR proteins at positions equivalent to those in the ORs have similar ligand-binding functions.

The present results are remarkable in view of the notable differences in the molecular structure and transmembrane configuration of these two receptor families (GRs have seven transmembrane receptors with an extracellular C-terminal domain, while IRs have an inverted topology with three transmembrane helices and an intracellular C-terminus). Similar results were obtained with the aphid GRs, in which more than 70% of positively selected sites in paralog comparison were identified in the extracellular region that includes the putative ligand-binding domains (Smadja et al., 2009). Such uneven distribution of positively selected sites across these proteins is reinforcing evidence of selective processes acting on key functional domains. Hence, we can hypothesize that functional diversification of chemosensory receptors would have been largely driven by adaptive changes affecting the ligand-binding specificity and/or sensitivity. Importantly, our results point to specific codons that can, in the future, be targeted in functional experiments to determine ligand specificity and channel characteristics in the GR and IR genes of S. maritima.

Very little is known about the molecular apparatus of the chemosensory system of non-insect arthropods. This is a first contribution that explores the potential of molecular evolutionary analyses to contribute to the understanding of functional diversification in chemoreceptors. We found that positive selection has had a role in the evolution of the GR and IR repertoires in S. maritima and that adaptive diversification has happened mostly in the ligand-binding interacting domains of these receptors. The next steps in the study of these gene families are the characterization of their expression patterns and of their ligand specificity. The identification of chemosensory genes in other non-insect arthropods would also greatly contribute to our understanding of chemosensation in these less studied organisms.

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.