Evolutionary Genomics of Genes Involved in Olfactory Behavior in the Drosophila melanogaster Species Group.

Previous comparative genomic studies of genes involved in olfactory behavior in Drosophila focused only on particular gene families such as odorant receptor and/or odorant binding proteins. However, olfactory behavior has a complex genetic architecture that is orchestrated by many interacting genes. In this paper, we present a comparative genomic study of olfactory behavior in Drosophila including an extended set of genes known to affect olfactory behavior. We took advantage of the recent burst of whole genome sequences and the development of powerful statistical tools to analyze genomic data and test evolutionary and functional hypotheses of olfactory genes in the six species of the Drosophila melanogaster species group for which whole genome sequences are available. Our study reveals widespread purifying selection and limited incidence of positive selection on olfactory genes. We show that the pace of evolution of olfactory genes is mostly independent of the life cycle stage, and of the number of life cycle stages, in which they participate in olfaction. However, we detected a relationship between evolutionary rates and the position that the gene products occupy in the olfactory system, genes occupying central positions tend to be more constrained than peripheral genes. Finally, we demonstrate that specialization to one host does not seem to be associated with bursts of adaptive evolution in olfactory genes in D. sechellia and D. erecta, the two specialists species analyzed, but rather different lineages have idiosyncratic evolutionary histories in which both historical and ecological factors have been involved.

Introduction

Animals rely on their senses to perceive the world; in this sense olfactory behavior is a crucial character that allows insects to perform tasks of ecological and evolutionary importance such as: finding food, mate recognition, oviposition site selection, avoidance of predators, etc. In the last 15 years a large amount of information about genetic, physiological and molecular aspects of olfaction has been accumulated in invertebrates1,2 and vertebrates.3–5

In the fruit-fly Drosophila melanogaster (Diptera: Drosophilidae), odorants bind to olfactory receptors (Or) that are expressed in olfactory receptor neurons (ORNs) membranes.2 At this point, odorant binding proteins (Obp) expressed in the aqueous lymph that fills the olfactory sensilla on the antenna6 also play important roles in this interaction, for example: helping in solubilisation, transport and presentation of odorants to Or’s and removal and rapid inactivation of odors after reaction with Or’s in ORNs membranes. 7 This first contact of flies’ olfactory molecules with odorants (that take place in the olfactory organs, the third antennal segment and the maxilar palp) constitutes the peripheral level of the Drosophila olfactory system.2 Nevertheless, other gene products, not expressed in the peripheral level, are also necessary for the proper expression of olfactory behavior.8–14 These genes, which encode proteins with quite different functions, occupy central positions in the Drosophila olfactory system. Therefore, to understand the genetic architecture of a phenotypic character identification of all genes, gene-to-gene interactions, gene-to- environment interactions that contribute to its expression is needed, as well as a comprehension of its variational properties.15,16 Several studies have shown that olfactory behavior has a complex genetic architecture since it is orchestrated by many interacting genes that also interact with the environment.8,9,14,17–21

Comparative genomic studies on the evolution of olfactory behavior in Drosophila have focused on particular gene families, such as Or22–26 and Obp gene families.27 Here, we present a comparative genomic study of olfactory behavior in Drosophila including an extended set of genes, besides Or’s and Obp’s gene families, that are known to be involved in the genetic architecture of olfactory behavior. We took advantage of the recent burst of whole genomes sequences and the development of powerful statistical and bioinformatic tools to analyze large sets of genomic data28,29 and tested evolutionary hypotheses regarding the genetic architecture of olfactory behavior in the six species of the D. melanogaster species group for which whole genome sequences are available, namely Drosophila melanogaster, Drosophila simulans, Drosophila sechellia, Drosophila erecta, Drosophila yakuba, Drosophila ananassae.30

Even though in our analysis we did not group genes in gene families, the hypotheses we aim to test are not about individual genes but about groups of genes that share a given characteristic that could help to explain the evolution of olfactory behavior. As pointed out earlier, olfactory genes participate in olfaction at different stages of the life cycle and/or at different positions in the chain of events that builds the olfactory response.2,31,32 In this context, we aim to test the following hypotheses: (i) olfactory genes that participate in the expression of both larval and adult olfactory behavior evolve at a slower rate than genes that contribute either to larval or adult olfactory behavior, since genes with ubiquitous expression are expected to be under stronger functional constraints; (ii) olfactory genes that participate in a peripheral position of the olfaction transduction pathway evolve faster than genes expressed in central positions, since central position genes tend to participate in other signal cascades, besides olfactory. This hypothesis follows from the general notion that more pleiotropic genes may be under stronger purifying selection than less pleiotropic genes.33

Comparative genomic studies between closely related species also offer the unique opportunity to compare the evolutionary rates in lineages with particular ecological and evolutionary histories. Sequenced species in the D. melanogaster species group can clearly be grouped in specialist and generalist species. D. sechellia breeds and feeds only on the decaying fruits of the shrub Morinda citrifolia34–36 and D. erecta is a seasonal specialist of Pandanus candelabrum. 37 Among the remaining species, we found two known quintessential generalists: D. melanogaster and D. simulans; and a less known generalist D. yakuba.38,39 Previous work, based on studies in Or and Obp gene families, proposed that olfactory genes in specialist lineages of the D. melanogaster species group should exhibit faster evolutionary rates than generalists due to directional selective pressures from a one-host environment.23,24,27 To verify these hypotheses we used an extended dataset of olfactory genes, a maximum likelihood model to estimate evolutionary rates at coding sequences and suitable tests aimed to differentiate true cases of positive selection (PS) from relaxation of selective constraints (RSC) or weak positive selection (WS). Our results show that evolutionary rates of olfactory genes are independent of the life cycle stage, and the number of stages in which they participate in olfaction. However, we discovered a relationship between the position that the gene product occupies in the olfactory signal transduction pathway and its evolutionary rate. Also, we found that in the D. melanogaster species group specialization to one host does not seem to trigger bursts of rapid evolution in olfactory genes, but reveals instead idiosyncratic evolutionary patterns in the two specialists species analyzed.

Materials and Methods

We performed a genomic comparative analysis of genes known to be involved in olfactory behavior recovered from the whole genome sequences of the six species of the D. melanogaster species group available so far. The analysis was done with coding sequences (CDS) of 115 D. melanogaster genes, known to be involved in olfaction, and the CDS of the orthologs in the five remaining species; adding to a total of 664 CDS (26 genes under a process of pseudogenization or/and misaligned were excluded from the analysis, so the total number dos not add to 690). CDS were downloaded from FlyBase (sequences were uploaded to FlyBase by40). CDSs were aligned with MUSCLE41,42 using predicted amino acid sequences as templates. Alignment columns containing gaps were removed.

Olfactory genes were categorized in two different functional classes, life cycle stage (larval and/or adult) in which a gene is involved in olfaction and the position (peripheral or central) in which participates in the olfactory system. Regarding the categorization of genes in terms of involvement in larval and/or adult olfactory behavior, we only included those that were tested for olfactory behavior in both life cycle stages. If a gene was reported to participate in the expression of adult olfactory behavior but was never tested in larvae it was not included in the study. The categorization in terms of position in the olfactory system is based on a conception of olfactory behavior as a complex trait orchestrated by many genes with a wide range of molecular functions. ORs and OBPs, which encode for proteins involved in the first contact of the fly with odorant stimuli, were considered as peripheral in the olfactory system. We consider as central a subset of genes with different molecular functions that are known to participate in olfaction but not in the first contact with odorants. This classification of olfactory genes is commonly used in studies of olfaction in Drosophila, as can be seen, for example, in a recent review of molecular aspects.2 The involvement of these genes in olfaction was identified by means of different genetic approaches as loss of function by mutagenesis, rescue experiments and expression techniques like RT-PCR, in situ hybridization and immunohistochemistry. The number of genes included in the analyses along with information on functional categories and the corresponding sources are given in Table 1 (a more detailed description of each gene categorization can be found in Supplementary Table S1).

Table 1

Olfactory genes included in the study.

Categorization	Nº of genes	Reference
Position
Central	53	9, 10, 11, 12, 13, 14, 64, 65, 66, 67
Peripheral	62	2 (and references therein), 68, 69, 70, 71
Stage
Larva	27	2 (and references therein), 14, 64, 68, 69
Adult	34	2 (and references therein), 13, 14, 64, 68, 69, 70, 71
Larva and adult	13	2 (and references therein), 14, 64, 68, 69, 70, 71

Notes: Place of participation of genes in the olfactory system (peripheral: genes that participate in the first contact with odor molecules in the periphery of olfactory response; central: genes known to be involved in the proper expression of olfactory behavior in down stream positions of the olfactory system). Stage: life cycle stage of participation in olfaction (L: larvae; A: adult; LA: larvae and adult). Reference: bibliographical reference from where the information to make the categorization was obtained from: 2 (and references therein), 9, 10, 11, 12, 13, 14, 64, 65, 66, 67, 68, 69, 70, 71. Categorization was based on results obtained in D. melanogaster. A more detailed list of genes used in this paper can be found in Supplementary Table S1.

Although we acknowledge that function, timing and tissue of expression of a given gene may change across species. We assume, given the close phylogenetic relationships among species of the D. melanogaster species group, that the functionality of olfactory genes (as described in D. melanogaster) did not change in any of the branches of the phylogeny of the D. melanogaster species group.

Maximum likelihood tests

Maximum Likelihood (ML) estimates of the ratio ω= dN/dS (dN: nonsynonymous substitution rate; dS: synonymous substitution rate) for each olfactory gene were employed to estimate evolutionary rates and to infer the evolutionary forces acting upon them. Synonymous mutations do not change amino acid sequence, hence dS is assumed neutral with respect to selective pressures on the protein product of the gene; and nonsynonymous mutations change the amino acid sequence, so dN depends on the selective pressures acting upon the protein sequences. The ratio ω= dN/dS, then becomes a measure of selective pressures. On one hand, if nonsynonymous mutations are deleterious, purifying selection will slow their fixation rate and ω values are expected to be lesser than 0.3; when 0.3 < ω< 1 it is interpreted as weak purifying selection. On the other hand, if nonsynonymous mutations are advantageous, they will fix at a higher rate than synonymous mutations, and thus, ω will be greater than 1 and indicative of PS. A ω= 1 stands for neutral evolution. Estimates of ω were obtained applying a free-ratio ML model using the CodeML program of the PAML 4 package.28 Genes with dS = 0 were excluded from the analysis in which we compare median values of ω between groups of categorized genes, since in a comparative framework, the resulting ω values are not biologically informative.

Positive selection was evaluated using two different branch-site models (A and A1) and two different likelihood ratio tests (I and II),43 implemented in the CodeML program of the PAML 4 package28 (from here on called the Zhang-Nielsen-Yang 2005 tests). Branches in the phylogeny were defined a priori as foreground and background lineages. Under these models only the foreground lineage may contain events of PS. Each species, except D. ananassae that was used as out-group, in the D. melanogaster species group and the respective ancestral lineage were tested independently as the foreground lineage. Ancestral sequences were reconstructed by ML from actual species sequences in the phylogenetic tree.

In contrast to the statistical behavior of previous branch-site tests,44 the Zhang-Nielsen-Yang 2005 tests are improved methods of branch-site tests using an ML approach which has proved to be more successful in differentiating PS from RSC.43 Test I compares model M1a against model A. M1a assumes two classes of sites, one with 0 < ωo < 1 and another with ω1 = 1, fixed in all the lineages of the phylogenetic tree. Model A considers four classes of sites: site class 0 which includes conserved codons throughout the tree (0 < ω0 < 1 in all lineages), site class 1 that includes codons evolving neutrally throughout the tree (ω1 = 1 in all lineages), site classes 2a and 2b which include codons conserved or evolving neutrally on the background branches but that become positively selected in the foreground branches (2a: 0 < ω0 < 1 in background lineages and ω2 > 1 in foreground lineages and 2b: ω1 = 1 in background lineages and ω2 > 1 in foreground lineages). The proportion pi of each site class (p0, p1, p2a, p2b) and the mean value of ω2 (ω2a and ω2b) can be estimated from the data using ML methods. Test II compares the null model A1 against model A. Parameters in A1 are equal to those of A with the exception that site classes 2a and 2b are fixed in the foreground with ω2 = 1. As was demonstrated using a simulation, Test I cannot properly distinguish cases of RSC from true events of PS. On the other hand, Test II, by allowing selectively constrained sites in the background to become relaxed under the proportion of site classes with ω2 = 1 set in the foreground of A1, can properly distinguish between these two alternatives, with an acceptable false discovery rate.43 Therefore, we evaluate the results obtained with both tests to distinguish between events of PS from RSC. Since the compared models are nested, likelihood ratio tests were performed and likelihood ratio tests statistics (2Δ l = − 2 [ln[likelihood for null model] – ln[likelihood for alternative model]]) were posteriorly transformed into exact P-values using the pchisq function of the R statistical package.45 Likelihood ratio tests were performed using a χ 2 distribution with d.f. = 2 for Tests I and d.f. = 1 for Test II, which have been shown to be conservative under conditions of PS.43 P-values derived from PS analyses were false discovery rate-adjusted using the method of Benjamini and Hochberg.46

We excluded from the analysis the other 6 species that not belong to D. melanogaster species group for which whole genome information is available, since the inclusion of highly diverged taxa may lead to saturation in synonymous positions, which may produce unreliable estimates of ω.

Phylogenetic relationships between species in the D. melanogaster species group are well resolved, with the exception of controversial results for some members of the melanogaster sub-group, specifically the placement of D. erecta and D. yakuba relative to the D. melanogaster lineage. We only show the results obtained using the best supported topology (D. melanogaster, (D. erecta, D. yakuba))40 (Fig. 1).

Figure 1

Maximum likelihood estimates of median ω of all olfactory genes in each lineage of the D. melanogaster species group analyzed.

For details of parameters and input data used in maximum likelihood models in CodeML see Appendix in Supplementary data.

Multiple regression analysis

In order to rule out the possibility that differences in ω values in our set of genes are related to other variables besides gene categorizations made by us we tested whether physical location of genes in the genome and/or the amount of expression are related to ω by performing a multiple linear regression analysis of ω (estimated by a free ratio model) with physical location (cytological band) and expression level (mRNA abundance in whole adult flies). mRNA abundance data in whole adult flies were obtained from FlyAtlas. 47 The multiple linear regression model we used is:

where Yi is the i-th observation of the dependent variable ω. X1i and X2i are i-th value of the independent variables physical location and expression level, respectively. β 1 and β 2 are unknown parameters representing the rates of change in Y per unit change in X1 and X2, respectively. β 0 is another unknown parameter that represents the intercept of the line.

All statistical analyses and tests were done using R statistical package.45

Results

Maximum likelihood estimates of ω using a free ratio-model revealed that most genes included in the analysis are under purifying selection (ie, ω< 0.3) in all branches of the phylogeny. Median ω values for the set of genes analyzed in each species are: ω = 0.0525, ω = 0.0477, ω = 0.0805, ω = 0.0879, ω = 0.0641, ω = 0.0449 (Fig. 1; Supplementary Table S1). However, four genes exhibited values of ω greater than 1, suggesting the occurrence of PS in their coding regions.

Inferring the forces affecting the evolution of olfactory genes in the D. melanogaster species group

Since estimates of ω from a free-ratio model are not appropriate test for detecting PS, we performed the Zhang-Nielsen-Yang 2005 tests43 for sites under PS in specific branches of the tree. In this case, none of the four olfactory genes that presented a ω> 1 in the free-ratio model reached significance for PS or RSC when tested using Zhang-Nielsen-Yang 2005 tests. The ω values, estimated by a free ratio model, of these genes are only slightly greater than 1, hence Tests I and II both failed to accept either PS or RSC. Instead, two olfactory genes, CG10777 in D. sechellia and Xrp1 in D. erecta, presented signals of PS; and seven showed signs of RSC or weak selection: CG11883, CG32556 and Moesin in D. simulans; discs lost and Or42a in D. erecta; CG16708 and Or82a in D. yakuba (Table 2). In general, the analysis performed on the subset of five closely related species showed that only 9 out of the 565 (1.6%) olfactory genes analyzed exhibited signs of either PS or RSC. The remaining set of genes (98.4%) showed values of ω that were not significantly different from neutral expectations or were compatible with purifying selection.

Table 2

Positive selection and relaxation of selective constraints on olfactory genes.

Gen	Test I		Test II		Biological interpretation: PS or RSC-WS
	d.f.	P-value	d.f.	P-value
CG10777D. sec	2	1	1	0.0037	Gen under PS
Xrp1D. ere	2	2× 10− 5	1	0.0007	Gen under PS
CG11883D. sim	2	0.044	1	1	Gen with signs of RSC-WS
CG32556D. sim	2	0.0054	1	0.1981	Gen with signs of RSC-WS
MoesinD. sim	2	0.0026	1	1	Gen with signs of RSC-WS
discs lostD. ere	2	0.0181	1	1	Gen with signs of RSC-WS
Or42aD. ere	2	0.0136	1	1	Gen with signs of RSC-WS
CG16708D. yak	2	0.0329	1	1	Gen with signs of RSC-WS
Or82aD. yak	2	0.0045	1	0.5590	Gen with signs of RSC-WS

Notes: Significant results for Zhang-Nielsen-Yang 2005 tests are shown. Gene names and species for witch the ortholog gene belongs (in sub-index) are shown.

Abbreviations: d.f., Degrees of freedom; PS, positive selection; RSC-WS, relaxation of selective constraints or weak positive selection.

Evolution of specific groups of olfactory genes in the D. melanogaster species group

The vast amount of genetic information that has accumulated for D. melanogaster in past decades allowed us to assign genes to a priori categories on the basis of: (i) the stage of the life cycle in which a particular gene participates in olfaction -genes can take part in either larval or adult olfactory behavior or in both-, and (ii) the position at which a gene product participates in the transduction pathway that characterizes the olfactory response -genes products can play a role either in the first contact with odor molecules in the periphery of the olfactory response (Or and Obp gene families) or in more central positions of the olfactory transduction pathway. Supplementary Table S1 lists all olfactory genes analyzed and their classification in the different categories.

In the following sections we describe the results of the comparative analyses of evolutionary rates among a priori categorized groups of genes to test the following hypotheses: (i) olfactory genes that participate in the expression of both larval and adult olfactory behavior evolve at a slower rate than genes that contribute to olfactory behavior in only one life cycle stage and (ii) genes that participate in peripheral positions of the olfactory system evolve at a faster rate than genes involved in central positions.

Rates of evolution of olfactory genes and their role in the life cycle

In D. melanogaster the median value of ω for genes involved in larval (L) olfaction (ω= 0.0742) was twice larger than the median ω estimated for genes involved in larval and adult (LA) (ω= 0.0389) and only adult (A) olfaction (ω= 0.0359) (Fig. 2). However, differences among medians were not significant (Fig. 2; Wilcoxon rank sum test, A vs. L: w = 614. P = 0.0243; A vs. LA: w = 184.5. P = 0.3889; L vs. LA: w = 210, P = 0.3249), suggesting that the rate of evolution is independent of the life cycle stage of expression and also of whether the gene participates in olfaction in one (either larval or adult) or in both stages of the life cycle.

Figure 2

Comparisons of gene group categories median ω in D. melanogaster. ML estimates of median ω for gene group categories (life cycle stage of participation in olfaction and position of participation in the olfactory system) for D. melanogaster (horizontal marks).

Notes: Gray circles are free ratio-model ω estimates for individual genes within a given category. To easier the lecture, the plot was truncated for ω values higher than 0.5. Significant differences (highlighted with an S) between median ω of peripheral position genes and central position genes were found using a Wilcoxon rank sum test. A Bonferroni multiple-comparison correction was applied.

We performed the same analysis in the other species of the D. melanogaster species group and observed significant differences between categories in D. simulans, D. erecta and D. ananassae, but not in the other lineages analyzed. In the first two species median ω was significantly closer to 1 for larval genes than for adult genes (Fig. 3; Wilcoxon rank sum test: A vs. L in D. simulans, w = 539, P = 0.0057; in D. erecta, w = 609, P = 0.0154), whereas the median value of ω for L genes was significantly closer to 1 than for LA genes in D. ananassae (Fig. 3; Wilcoxon rank sum test: L vs. LA, w = 258, P = 0.0165). These tests revealed that olfactory genes known to be expressed only in the larval stage evolve faster relative to genes expressed in adults or in both stages, in three of the six species analyzed (Fig. 3). Such difference may be a reflection of different selective pressures (or functional constraints) among species on olfactory genes which help flies to perceive the environment.

Figure 3

Comparisons of the ω for olfactory genes expressed in different life cycle stages for the six species of the D. melanogaster species group. ML estimates of median ω for defined gene groups according to life cycle stage of participation in olfaction for the six species of the D. melanogaster species group.

Notes: Error bars are SEM values (Standard Deviation/N1/2, where N is the number of genes). A Bonferroni multiple-comparison correction was applied. Significant P-values are marked with an S.

Olfactory genes and position in the olfactory transduction pathway

Genes expressed in a central position of the olfactory transduction pathway exhibited a median ω (ω= 0.0326) that was significantly lower than for genes expressed in a peripheral position (ω= 0.0785) in D. melanogaster (Fig. 2; Wilcoxon rank sum test, w = 987, P = 2.276.10−4). The same is true in the other species of the D. melanogaster species group (Fig. 4; Wilcoxon rank sum test: in D. sechellia, w = 1696.5, P = 0.0111; in D. yakuba, w = 2407, P = 1.82 × 10−5; in D. erecta, w = 430, P = 1.54 × 10−8; in D. ananassae, w = 947, P = 9.52 × 10−5) with the only exception of D. simulans. These results suggest stronger selective constraints acting on central than on peripheral genes, a pattern that seems to be independent of the particular evolutionary history and ecological context that characterize the species that belong to the D. melanogaster group.

Figure 4

Comparisons of the ω for olfactory genes that participates in different positions in the olfactory system for the six species of the D. melanogaster species group. ML estimates of median ω for defined gene groups according to position of participation in the olfactory transduction pathway for the six species of the D. melanogaster species group.

Notes: Error bars are SEM values (Standard Deviation/N1/2, where N is the number of genes). A Bonferroni multiple-comparison correction was applied. Significant P-values are marked with an S.

However, these results may be a consequence of the physical location of the genes in the genome and/or the amount of expression and not to the position they occupy in the olfactory system. To rule out this possibility we performed a multiple regression test to evaluate in D. melanogaster whether there is a relationship between ω values and cytological position and/or amount of expression and found non-significant results (F2,112 = 0.58, P = 0.564).

Evolutionary rates in generalist vs. specialist lineages

D. sechellia breeds and feeds on decaying fruits of the shrub Morinda citrifolia34–36 and D. erecta is a seasonal specialist on Pandanus candelabrum.37 Previous works reported that olfactory genes evolve at a faster rate in specialists than in generalists as a consequence of directional selective pressures from a one-host environment23,24,27. In order to investigate this hypothesis we conducted two different tests. First, we compared ω values between generalist and specialist species for a priori defined categories of olfactory genes as described above (peripheral and central position) and used the rest of protein coding genome as control. This organization of the dataset allows us to distinguish between genome wide effects due to demographic events from local effects of PS. Specialists exhibited a significantly faster evolutionary rate than generalists when all olfactory genes were considered (Table 3; Wilcoxon Rank Sum test, Specialists vs. Generalists: w = 3.182, P = 0.00146). However, when we refined the analysis by dividing olfactory genes into peripheral and central, we found that differences between specialists and generalists are only accounted for by peripheral position genes (Table 3; Wilcoxon Rank Sum test, Specialists vs. Generalists: w = 3.793, P < 0.001) and that evolutionary rates of central position genes did not differ between specialists and generalists (Table 3; Wilcoxon Rank Sum test, Specialists vs. Generalists: w = 0.841, P = 0.4). Interestingly, specialists also presented a genome wide acceleration of evolutionary rates (Table 3; Wilcoxon Rank Sum test, Specialists vs. Generalists: w = 25.277, P < 0.0001) relative to generalists. Moreover, comparisons between each specialist against the group of generalists revealed that D. erecta exhibited such acceleration (Table 3; Wilcoxon Rank Sum test, all olfactory genes: w = 3.122, P = 0.0018; peripheral position genes: w = 4.307, P < 0.0001; central position genes: w = 0.141, P = 0.888; rest of protein coding genome: w = 24.797, P < 0.0001) but not D. sechellia (Table 3; Wilcoxon Rank Sum test, all olfactory genes: w = 1.864, P = 0.0623; peripheral position genes: w = 1.609, P = 0.108; central position genes: w = 1.231, P = 0.2183; rest of protein coding genome: w = 14.922, P < 0.0001).

Table 3

Evolutionary analysis of olfactory genes in generalists vs. specialist’s lineages in the D. melanogaster species group.

Species	Peripheral position genes			Central position genes			All olfactory genes			Rest of protein coding genome
	Genes per group	Deviation	P	Genes per group	Deviation	P	Genes per group	Deviation	P	Genes per group	Deviation	P
Specialists vs. Generalists	117–779	0.0426	<0.001	100–052	0.0047	0.4	217–731	0.0315	0.00146	15968–23959	0.0174	<0.0001
D. erecta vs. Generalists	61–179	0.0496	<0.0001	53–152	−0.0052	0.888	114–431	0.0358	0.0018	8288–23959	0.0184	<0.0001
D. sechellia vs. Generalists	56–679	0.024	0.108	47–752	0.0268	0.2183	103–331	0.0266	0.0623	7680–23959	0.016	<0.0001

Notes: Specialists species are D. erecta and D. sechellia; generalist species are D. melanogaster, D. simulans and D. yakuba (D. ananassae was not include in this between-species comparison because, as a results of being used as out-group, the ω estimation for this specie tend to be overestimated compared to the rest). P-values are from Wilcoxon Rank Sum test. A Bonferroni multiple-comparison correction was applied. Significant P-values are in bold.

Secondly, we specifically evaluated if there has been an acceleration of non-synonymous substitution rates in the D. erecta and D. sechellia lineages. To this end, we took advantage of the genome sequences of D. yakuba, a generalist close relative of D. erecta; and D. simulans, a generalist close relative of D. sechellia, to infer the sequences of the corresponding common ancestors. These analyses revealed a significant acceleration of the non-synonymous substitution rate in the branch leading to D. erecta for olfactory genes occupying a peripheral position in the olfactory transduction pathway. However, it should be noted that a similar trend was detected for the rest of the protein coding genome (Fig. 5; Wilcoxon Matched Pairs test, peripheral position genes: N = 53, Z = 5.83, P < 0.0001; rest of protein coding genome: N = 8369, Z = 28.526, P < 0.0001). Another interesting result was that a similar acceleration of the non-synonymous substitution rate was found in the branch leading to the generalist D. yakuba (Fig. 5; Wilcoxon Matched Pairs test, peripheral genes: N = 62, Z = 3.033, P = 0.0024; central genes: N = 53, Z = 2.518, P = 0.0118; rest of protein coding genome: N = 8369, Z = 28.632, P < 0.0001). Therefore, there are footprints of PS for olfactory genes in D. erecta, however, our study also points out that the acceleration of the non- synonymous substitution rate acceleration may not be related to specialization on one host. In contrast, the replacement substitution rate did not appear to be accelerated neither in the branch leading to the specialist D. sechellia (Fig. 5; Wilcoxon Matched Pairs test, peripheral genes: N = 62, Z = 1.606, P = 0.108; central genes: N = 53, Z = 1.171, P = 0.241), nor in the branch of the generalist D. simulans (Fig. 5; Wilcoxon Matched Pairs test, peripheral genes: N = 62, Z = 0.382, P = 0.702; central genes: N = 53, Z = 0.08, P = 0.936). Similarly, we did not detect an acceleration of the evolutionary rates of olfactory genes in the branch leading to D. sechellia/D. simulans common ancestor (Fig. 5; Wilcoxon Matched Pairs test, peripheral genes: N = 62, Z = 0.022, P = 0.982; central genes: N = 53, Z = 0.781, P = 0.435).

Figure 5

Non synonymous substitution rate acceleration analyzes in clades including specialist and generalist species. The mean normalized nonsynonymous substitution rate difference (MNdND) was calculated for specialist lineages D. erecta and D. sechellia, neighbor generalist lineages and for the reconstructed sequence of the common ancestor of D. simulans/D. sechellia and D. erecta/D. yakuba. MNdND = (Mean dN actual linage/Standard Deviation dN actual linage) – (Mean dN ancestor/Standard Deviation dN ancestor). MNdND > 0 indicate non synonymous substitution rate acceleration in derived lineages, MNdND ≤ 0 indicate the absence of non synonymous substitution rate acceleration in derived specialist lineages.

Notes: In red are highlighted the comparisons showing a MNdND significantly different from zero using a Wilcoxon Matched Pairs test. A Bonferroni multiple-comparison correction was applied. Branch lengths are artificial.

Discussion

With regard to general evolutionary patterns we found that purifying selection is the main evolutionary process shaping olfactory genes evolution in the six species of the D. melanogaster species group investigated (Fig. 1; Supplementary Table S1), in agreement with previous studies that only surveyed Or and Obp gene families.22,24–27 However, applying a ML conservative test that effectively distinguishes true cases of PS from RSC, we detected signs of PS in only two olfactory genes, one in D. sechellia (CG10777) and the other in D. erecta (Xrp1) (Table 2). The remaining genes involved in the olfactory system did not show departures from neutral expectations or showed patterns compatible with purifying selection. Regarding studies reporting the occurrence of PS in olfactory genes in D. melanogaster there are conflicting results. Two studies did not find evidence of PS acting upon Obp’s and Or’s26,27 and another two claimed that ten Or genes (Or9a, Or10a, Or19a, Or43a, Or56a, OrN1, OrN2, Or33c, Or42a, Or85e) exhibit the hallmark of PS.22,25 In the present study we included six of the Or genes claimed to have evolved under PS (Or9a, Or10a, Or33c, Or42a, Or43a, Or56a) and found no trace of PS. Although ML codon-models were applied in all studies, subtle differences among studies (genes showing different results among studies are only a minor subset of the Or gene family) may exist because different approaches aimed at detection of positive selection may lead to slight different results even when the same, or very similar, genomic data and phylogeny are used. The aforementioned studies showed that the number of genes found to be under PS is a quite small fraction of the total number of genes analyzed. In our analysis we surveyed a large number of genes involved in olfactory behavior (not just Or and/or Obp genes) and found that only two showed traces of PS and seven exhibited signs of RSC out of a total of 565 genes in five closely related species. More than 98% of our set of genes is evolving neutrally under purifying (negative) rather than positive selection. Therefore, we concur with a recent review of genomic studies of chemosensory gene families stating that “although appealing, adaptationist interpretation might not be justified”.48 Thus, the proposition of a scenario where olfactory behavior is evolving mainly by PS despite weak empirical support should be avoided.

We also tested hypotheses concerning the evolution of genes involved in olfactory behavior. A recent genome-wide study in D. melanogaster revealed a gradient of divergence rates over ontogeny, the earlier the expression of genes the lower the evolutionary rate (ie, embryonic < larval/pupal < adult).49 Moreover, the lavish literature on olfactory behavior is consistent with the idea that larval and adult olfactory behavior are partially decoupled in terms of genetic architecture,20,31,32,50 physiology,51 external anatomy,52,53 and molecular neuroanatomy.2 These observations suggest that evolutionary rates may differ among genes that participate in olfaction in different stages of the life cycle as a result of differential functional constraints and different environmental pressures. Nevertheless, our results show that evolutionary rates of olfactory genes appear to be largely independent of both the life cycle stage and the number of stages in which they are expressed in D. melanogaster (Fig. 2). A similar finding was reported for Or genes24 and confirmed for a wide set of olfactory genes in the present study. Thus, other outcomes of evolutionary processes, rather than changes in coding sequences, should be invoked to account for the partial decoupling of larval and adult olfactory behavior. Then, gene gains and losses22,27 and/or changes in expression patterns may be the key mechanisms underlying these differences.54 We also evaluated the same hypothesis mentioned above in a larger dataset including all species of the D. melanogaster species group for which complete genomes are available. These tests revealed a significant increase in the pace of the evolution of olfactory genes expressed only in larvae (relative to genes expressed in adult and in both stages) in D. simulans, D. erecta and D. ananassae but not in D. melanogaster, D. sechellia and D. yakuba (Fig. 3). Such among-species difference may result from different evolutionary histories reflecting differential selective pressures (or functional constraints) on olfactory genes which help flies to perceive the environment. Particularly important in fruit flies may have been the evolutionary history of host plant utilization. 23,38,39 Alternatively, but not mutually exclusive, differences in population structure across species, either related or unrelated to host usage evolution, may affect the relative roles of natural selection, gene flow and drift55,56 and account for differences in substitution rates heterogeneity across species between genes that participate in different stages of the life cycle.

Concerning comparisons between genes that participate in different positions of the olfaction transduction pathway, stronger purifying selection was detected on central than on peripheral position genes in all species except D. simulans (Fig. 4), a pattern that emerges as reliable and independent of the particular historical and ecological context characterizing the species studied. The most plausible explanation for this observation is that peripheral genes (Or and Obp) only participate in the expression of olfactory behavior, while central genes tend to be more pleiotropic, since they participate not just in olfaction, but also in signal cascades or gene networks that determine other phenotypic traits. Then, more pleiotropic genes tend to be more constrained and evolve at a slower pace. It is well known that pathway-central genes are under strong functional constraints due to extensive pleiotropy.33 In this sense, our present results along with empirical evidence in Drosophila57 and fungi,58,59 and also theoretical work60,61 give strong support to the general hypothesis that more pleiotropic genes evolve at a slower pace.

Finally, we tested the hypothesis that evolutionary rates of olfactory genes may be accelerated in specialists (D. sechellia and D. erecta) as a response to a one host environment as compared to generalist species (D. simulans, D. melanogaster, D. yakuba and D. ananassae). Our comparisons between specialists and generalists species showed that peripheral olfactory genes (but not central) evolved faster in specialists. However, further analysis revealed that this is a genome wide trend and not a particular feature of olfactory genes. We further refine our analysis by comparing D. sechellia and D. erecta individually to the group of generalist species and found that the acceleration in peripheral (but not central) genes can only be detected in D. erecta. In contrast, the genome wide acceleration as well as the lack of acceleration in central genes emerged as features common to both specialists (Table 3). These shared patterns may be a consequence of stronger functional constraints of central genes since they are highly pleiotropic. Moreover, our results concur with current views of network analysis in biology that propose that certain nodes in a network (in this case olfactory central position genes) are ‘essential’ and thus less prone to change.62 D. erecta and D. sechellia represent two independent events of specialization38 in which the evolutionary histories of olfaction genes exhibit certain subtleties that should be mentioned. On one hand, the analysis of non-synonymous substitutions rates in the clade D. simulans/D. sechellia revealed a pattern consistent with the idea of a genome wide acceleration restricted to the D. sechellia lineage, but not particularly for olfactory genes (Fig. 5). Thus, our results suggest that a scenario of directional selective pressures from a one-host-environment driving of the evolution of olfactory genes is not entirely sustainable in D. sechellia. Therefore, in agreement with Gardiner et al. (2008)26 we propose that historical demographic events like small effective population size after speciation from a mainland population55 and/or historical low levels of genetic variation (ie, resulting for historical events not from selective pressures)56,63 cannot be ruled out as determinants of the evolution of the olfactory behavior genetic basis in the specialist D. sechellia. In this vein, Lachaise and Silvain argued, on the basis of historical evidence, that M. citrifolia may not be the ancestral host for this species, which probably evolved on Pandanus candelabrum, and invaded M. citrifolia, which is an introduced species in the Seychelles in recent times.38 On the other hand, the D. yakuba/D. erecta clade depicted a different scenario. Comparisons of each species against the common ancestor revealed a genome wide acceleration of non-synonymous substitutions rates, as well as peripheral genes in both lineages. As in D. sechellia, our results argue against the one-host directional selection hypothesis. Moreover, our phylogenetic comparative analysis reveals that olfactory central genes are an exception to the genome wide acceleration in D. erecta, suggesting that these genes are under stronger constraints in the lineage leading to this specialist. Although there appears to be a relationship between ecological specialization and the evolution of olfactory behavior genetic basis, this relation does not seem to be specifically related to an acceleration of evolutionary rates of olfactory genes. However, the generalized acceleration of evolutionary rates detected in D. erecta may provide new insights to investigate ecological specialization in a wide genomic context. Having complete genomes of other drosophilids with a habitat-specialist or narrow niche-width could help to elucidate the relationship between specialization and genome evolution.

Conclusion

In conclusion, our study shows that evolutionary rates of olfactory genes are mostly independent of the life cycle stage and also of the number of life cycle stages in which they participate in olfaction and, more importantly, that evolutionary rates depend on the position (peripheral or central) that genes occupy in the olfactory system. Our approach is based on a conception that considers olfaction as a complex trait whose genetic bases include Or, Obp and a heterogeneous set of genes (identified as olfaction genes by means of an ample suite of methodological approaches) that do not participate in the periphery of the olfactory system. Actually, our study includes a larger set of genes identified as part of the genetic architecture of olfactory behavior than previous reports22,24–27 allowing us to demonstrate that genes categorized as central in the olfactory system evolve at a slower rate than peripheral genes. In addition, we show that, at variance with expectations of the one host adaptation hypothesis, the analysis of evolutionary rates of olfactory genes does not reveal common patterns in D. erecta and D. sechellia, the two specialists investigated. Our results suggest that specialization to one host does not seem to trigger bursts of rapid evolution in olfactory genes. The relevance of our present report is twofold. Firstly, we confirm the general biological tenet that more pleiotropic genes tend to evolve at a slower pace. Secondly, on a more specific feature concerning the evolution of olfaction in Drosophila we found that positive selection has had a limited incidence and that acceleration of substitution rates in lineages that underwent a process of ecological specialization is a feature extended to the entire protein coding genome not exclusive of olfactory genes.