Patterns of Positive Selection in Seven Ant Genomes

1Ants constitute an extremely successful lineage of animals which has colonized virtually all ecosystems on Earth. The pivotal feature at the basis of this ecological success is their highly social system with a reproductive division of labor, where one or a few queens specialize in reproduction, whereas workers conduct all the colony tasks such as brood care, nest maintenance, and food collection. In this article, we take advantage of the recent availability of seven sequenced ant genomes to perform a genome-wide scan for positive selection on amino-acid changes in protein coding genes during the evolution of the ant lineage. We addressed three main questions.

2First, we compared the amount of positive selection in functional categories of genes. Previous large-scale scans for positive selection in animals indicated that positive selection predominantly affects certain types of genes, such as those involved in evolutionary arms races, sexual selection, or conflicts with pathogens. Such genes experienced positive selection events recurrently on broad evolutionary time scales, and it is likely that they contribute to a fraction of the positive selection events that occurred in the ant lineage. To identify these genes, we reasoned that they likely also were under positive selection in other insect lineages. A systematic comparison of the targets of positive selection from published studies in insects is not straightforward because genome-wide scans for positive selection were often performed with different methods in different lineages. For example, a positive selection scan on 12 Drosophila species (all solitary) used the site test of Codeml, which is aimed at detecting recurrent positive selection events affecting particular sites of a protein, whereas a scan on 10 bee species (including solitary, primitively social, and highly social species) used the branch test, which tends to detect positive selection events affecting a large number of sites of a protein but during a limited period of time. To perform a robust comparison of the genes that were under positive selection in ants and other insects, we conducted similar scans for positive selection in ants and the flies and bees outgroups. An example of genes expected to be repeatedly under positive selection in insects are genes involved in defense and immunity. On the basis of the observed smaller set of immunity genes in the honeybee compared with Drosophila melanogaster, it has been suggested that selective pressure on these genes might have been relaxed in social insects, perhaps because they have social hygienic behaviors. However, the addition of several newly sequenced insect genomes revealed that the important gene complement in fruit fly is a derived character. We used our data sets to test whether there was evidence for weaker positive selection on immunity genes in ants and bees compared with flies.

3Next, we aimed at detecting sets of genes involved in functions likely to reflect ant-specific adaptations. We focused on three main adaptations. The first relates to the wide range of coordinated collective behaviors associated with division of labor in ant societies. Complex cooperative behaviors occur among nestmates for tasks such as communal nest construction and defense, brood rearing, social hygienic behavior, and collective foraging. It has been suggested that the evolution of social interactions may be tracked down to molecular changes affecting nervous system development and function. In particular it may translate into increased rates of positive selection on nervous system-related genes, as documented in primitively social lineages of bees, which evolved social behaviors independently from ants. Complex collective behaviors also require efficient communication systems that are essentially mediated by chemical signaling in social insects. Ants identify nestmates from non-nestmates, as well as ants from other species, through their scent. Individuals also use various types of pheromones as alarm signals and to mark their trails and territories. It has therefore been suggested that genes involved in chemical signaling, notably pheromone production and perception, should experience increased positive selection in ants compared with solitary insects. A manually curated data set of 873 olfactory receptor genes (ORs) allowed us to conduct a test for increased positive selection on these genes in ants.

4The second type of potential molecular adaptation relates to phenotypic plasticity among castes. Although queens and workers usually develop from totipotent eggs, they display dramatic morphological and physiological differences. Queens are often larger, have wings, and have much more highly developed ovaries than workers that often are sterile and lack a sperm storage organ. In most species, the differences between castes result from developmental differences induced by environmental factors rather than genetic differences. We therefore investigated whether there was evidence for increased positive selection in genes and pathways potentially involved in developmental plasticity.

5A third and interesting type of ant-specific adaptation relates to the extremely long lifespan of ant queens, which can live more than 20 years in some species. This corresponds to a 100-fold increase in lifespan compared with solitary insects. The variation in lifespan among castes is also remarkable, with queens living up to 10 times longer than workers and 500 times longer than males. So far, a limited number of molecular candidates have been identified to explain this pattern, mainly inspired from work in Drosophila. We therefore investigated whether there was evidence of positive selection on genes that have previously been associated with aging in model organisms. It is possible that positive selection acted on the same sets of genes in the bee lineage, where queens also live longer than other castes and than solitary insects, but such a signal should not be observed in short-lived species of the Drosophila lineage. To further assess the link between positive selection and aging, we investigated whether genes that experienced positive selection in ants were genes shown in D. melanogaster to be differentially expressed between old and young individuals, and between oxygen-stressed and control individuals.

6Finally, we investigated whether there was a difference in the level of positive selection between genes showing biased expression in queens, workers, and males. The efficiency of natural selection acting on an advantageous mutation—and thus the probability of its long-term fixation—is proportional to its effect on fitness. The fitness effects of mutations in genes that are expressed only in nonreproductive workers are indirect, so everything else being equal, selection should be less efficient at fixing them than mutations on genes expressed in queens and males. This could translate into lower levels of positive selection on genes expressed specifically in workers compared with males and queens. We therefore analyzed previously published microarray data from the red fire ant Solenopsis invicta, and compared the amount of positive selection between groups of genes varying in the level of caste-biased expression.

7To detect positive selection episodes that acted on protein-coding genes during the evolution of the ant lineage, the branch-site test of Codeml was run on 4,261 protein alignments of single-copy orthologs composed of four to seven ant and three to five outgroup species (see Materials and Methods). All branches that led to ant species in each gene family tree (including 2 hymenopteran and 13 ant branches; fig. 1) were successively tested for the presence of episodic positive selection. As many as 1,832 single-copy orthologs families (43%) displayed a signal of positive selection (at 10% false discovery rate [FDR]) in at least one of the branches tested (supplementary table S1, Supplementary Material online). In 91% of the significant alignments, at least one residue targeted by positive selection could be identified with a posterior probability greater than 0.9 (Bayes Empirical Bayes test; fig. 2). There was evidence for positive selection in at least one branch of the ant lineage for 830 (20%) of the genes analyzed. For 74% of them positive selection was specific to ants and not observed in the basal hymenopteran branches #7 and #8 (fig. 1). The 10 gene families with the most significant test values in the ant lineage are given in table 1.

8The proportion of positively selected genes varied significantly across the different branches tested (χ2 test, P < 1e-15; table 2), similarly to previous analyses with experimental and simulated data sets. This likely results at least in part from lower power of the branch-site test in shorter branches. Consistent with this view, there was a significant correlation in our data set between the length of tested branches and the test score (log-likelihood ratio; Spearman correlation ρ = 0.41, P < 1e-15). Additional analyses ruled out the hypothesis that false positives caused by convergence problems of the test, selective constraints acting on synonymous sites, saturation of synonymous substitution rate dS, or sequencing errors could be responsible for this pattern (supplementary text, Supplementary Material online).

9Taken together, these results demonstrate that positive selection was common in the evolution of the ant genes. The proportion of significant genes was similar in magnitude in the outgroup data set of 10 bees analyzed with the same methodology (20%; supplementary table S2, Supplementary Material online), but even higher in the outgroup data set of 12 flies (36%; supplementary table S3, Supplementary Material online). This difference might reflect biological differences between the lineages, such as effective population size NE, but also differences in the topology and branch lengths of the species trees, which influence the power to detect positive selection events in protein alignments (see supplementary text, Supplementary Material online).

10To compare the amount of positive selection experienced by different functional categories of genes, we classified genes based on their Gene Ontology (GO) annotation in D. melanogaster orthologs, and performed a gene set enrichment test using for each gene family a score reflecting the overall occurrence of positive selection in the ant lineage (Materials and Methods; supplementary text, Supplementary Material online). Such an approach of grouping genes enables a more sensitive search for positive selection, while buffering the impact of potential false positives (e.g., from remaining alignment errors or GC-biased gene conversion events which are difficult to distinguish from real positive selection signals; see supplementary text and supplementary tables S20–S24, Supplementary Material online). Twenty-four functional categories of genes were significantly enriched for positively selected genes in the ant lineage (at 20% FDR; table 3). A large number of them (11 out of 24) were related to mitochondria and mitochondrial activity. The other significant categories were related to nervous system development, behavior, immunity, protein translation and degradation, cell membrane, and receptor activity. Thus, positive selection apparently targeted a diverse array of gene functions during the evolution of the ant lineage.

11To identify GO categories that experienced positive selection not only in ants but also in other insects, we reanalyzed the fly and bee data sets with the same methodology used for the ant data set. These analyses revealed 106 GO categories significantly enriched for flies and 38 for bees (tables 4 and 5; supplementary tables S4 and S5, Supplementary Material online). We investigated which categories were enriched for positively selected genes in the three lineages. The first group of genes commonly enriched in ants, flies, and bees was related to proteolysis. This group included 4 of the 24 significantly enriched GO categories in ants (“proteolysis,” “metallopeptidase activity,” “peptidase regulator activity,” and “hydro-lyase activity”), 8 of the 106 GO categories enriched in flies (“serine-type endopeptidase activity,” “endopeptidase activity,” “proteolysis,” “metalloendopeptidase activity,” “peptidase activity,” “peptidase activity, acting on L-amino acid peptides,” “metallopeptidase activity,” and “exopeptidase activity”), and 6 of the 38 GO categories enriched bees (“amine metabolic process,” “metalloendopeptidase activity,” “metallopeptidase activity,” “signal transduction,” “cellular amine metabolic process,” and “cellular amino acid metabolic process”).

12The second group of genes enriched for positive selection signal in ants, flies, and bees was involved in response to stimuli. There was an enrichment of the GO category “receptor activity” in the three lineages as well as the GO categories “transmembrane receptor activity” and “olfactory receptor activity” in flies. This class of genes plays a pivotal role in the interactions between individuals and their environment. In addition, the GO categories “response to biotic stimulus” and “response to other organism” were enriched in flies, and the GO category “response to hormone stimulus” was enriched in bees. In ants, “response to ecdysone” and “response to steroid hormone stimulus” were marginally significant (FDR = 21%; supplementary table S6, Supplementary Material online).

13Some functions were enriched for positively selected genes in only two of the three lineages. These included GO categories related to immunity that were enriched in ants and flies, and some categories related to metabolism which were enriched in flies and bees. Evidence for positive selection on immunity-related functions in ants came from a significant enrichment of the GO categories “immune system development” and “hemopoietic or lymphoid organ development,” the organ that produces during larval development the cells mediating the immune response in insects. Seven GO categories related to immunity were also enriched in flies (“defense response,” “immune system process,” “regulation of antimicrobial humoral response,” “regulation of immune effector process,” “regulation of defense response,” “immune response,” and “antimicrobial humoral response”). The absence of significant enrichment for related categories in bees might reflect a lack of power of the gene set enrichment test, because the set of immunity genes is small in the honeybee and the data set analyzed was further depleted in genes with immunity-related functions (supplementary text and table S7, Supplementary Material online). Consistent with this interpretation, there was a trend in the direction of an enrichment, although nonsignificant, for 5 of the 7 tested GO categories related to immunity in bees (data not shown), suggesting that immunity might be a common target of positive selection in insects.

14The second set of GO categories enriched in two of the three insect lineages included various metabolic processes and their regulators, with metabolism of chitin, aminoglycan, carbohydrate, polysaccharide, glucose, hexose, glycerolipid, and phosphatidylinositol being enriched in flies and metabolism of lipid, amino-acid, nucleotide, and phosphorus being enriched in bees. There was no significant enrichment for GO categories related to metabolism in ants, but some categories were close to significance (e.g., “chitin metabolic process” and “rRNA metabolic process,” with FDR = 21% and 24%, respectively; supplementary table S6, Supplementary Material online). Metabolic functions, such as amino-acid, fatty acid, lipid, or RNA metabolism, were significantly enriched in ants when we used KEGG pathways annotation instead of the GO to perform the gene set enrichment test, as well as when the single-copy orthologs data set was reanalyzed with another multiple alignment method and a different quality filtering method (supplementary tables S25 and S26 and supplementary text, Supplementary Material online). It thus seems that metabolism is a common target of positive selection in insects.

15Two of the GO categories enriched in ants (“neuron recognition” and “adult locomotory behavior”) might potentially be linked to the evolution of neural systems and behavior (table 3). The first category was “neuron recognition.” However, GO categories related to neural systems were also enriched in a nonsocial hymenoptera lineage (“regulation of synaptogenesis,” “mushroom body development,” and “memory” on branch #8, basal to the Hymenoptera; fig. 1 and supplementary table S9, Supplementary Material online), and in the branches leading to primitively social lineages of bees (“synapse,” “synapse organization,” “regulation of synaptic growth at neuromuscular junction”; data not shown) (also reported in ), suggesting that positive selection on neural system genes in ants might not be directly associated with the emergence of social behaviors.

16The second GO category enriched in ants was “adult locomotory behavior.” This category was not enriched in any of the other tested lineages. The three genes contributing most to the positive selection signal in this GO category were DCX-EMAP, turtle, and beethoven. Mutational analyses of these genes in Drosophila suggest that they play an important role in sensory perception functions. Adult flies carrying a piggyBac insertion in DCX-EMAP are uncoordinated and deaf and display loss of mechanosensory transduction and amplification. Turtle plays an essential role in the execution of coordinated motor output in complex behaviors in flies, notably regarding the response to tactile stimulation. Finally, beethoven is involved in male courtship behavior, adult walking behavior, and sensory perception of sound in flies. This suggests that positive selection might have been important for the evolution of sensory perception functions in ants.

17A specific analysis of ORs did not provide support for the evolution of sociality being associated with increased levels of positive selection on ORs. A scan for positive selection across branches of a tree gathering 873 manually annotated ORs from two ants (Pogonomyrmex barbatus and Linepithema humile) and the solitary wasp Nasonia vitripennis (see Materials and Methods) revealed that positive selection was pervasive, with 277 branches (23%) displaying significant signals for positive selection (fig. 3 and supplementary fig. S1, Supplementary Material online). However, positive selection was detected in only 19% of the 929 branches leading to ant species, whereas as many as 40% of the 156 branches leading to wasps were under positive selection (Fisher's exact test P = 7.6e-9).

18None of the GO categories enriched for positively selected genes in ants could be linked to phenotypic plasticity (i.e., caste differences). In particular, there was no evidence of a significant enrichment for GO categories related to morphology or morphogenesis in the ant lineage. Another enrichment test using annotations obtained from mutant phenotypes in D. melanogaster (which are more relevant than GO annotations to analyze genes involved in morphogenesis since genes sets mostly refer to anatomical structures) also provided no clear support for positive selection on genes associated with phenotypic plasticity in ants (supplementary table S12, Supplementary Material online).

19However, among the genes with the highest support for positive selection in ants (table 1), two genes had a role in wing development (Guanine nucleotide exchange factor in mesoderm and Methylthioribose-1-phosphate isomerase) and one in larval development (Dis3), suggesting that even if positive selection did not act consistently on large sets of genes related to morphogenesis, it could have acted specifically on a few individual genes.

20Eleven GO categories enriched for positively selected genes in ants were related to mitochondrial activity (e.g., “mitochondrial electron transport,” “mitochondrial matrix,” “mitochondrial respiratory chain,” “NADH dehydrogenase ubiquinone activity,” and “oxidoreductase activity”; table 3). The mitochondrial processes under positive selection were not restricted to respiration and energy production, but also included translation (“organellar ribosome,” “mitochondrial small/large ribosomal subunit”). GO categories related to mitochondria were also enriched for positively selected genes on many individual branches of the ant lineage analyzed separately (supplementary table S9, Supplementary Material online), and in a larger data set including duplicated genes analyzed with the site test (see Materials and Methods; table 1, supplementary tables S13 and S19, Supplementary Material online). This indicates that recurrent events of positive selection occurred on genes with mitochondrial functions during the evolution of the ant lineage. In contrast, mitochondria-related GO categories did not display any enrichment for positively selected genes in flies and bees (tables 4 and 5), despite a high power to detect it on the respective data sets (supplementary text, Supplementary Material online). Similarly, no mitochondrial function was significantly enriched in the branches #7 and #8, basal to the ant lineage (fig. 1 and supplementary table S9, Supplementary Material online), reinforcing the idea that increased positive selection on mitochondria is restricted to the ant lineage.

21Of note, none of the 13 protein coding genes from the mitochondrial genome was included in our main data set because the mitochondrial genomes of most of the ant species analyzed were not annotated. Our results thus reflect positive selection on nuclear genes encoding proteins that function in the mitochondrion. We annotated mitochondrial genomes in 5 of the 7 ant species analyzed and tested whether positive selection could also be detected on the mitochondrial genomes themselves (Materials and Methods). However, we did not find evidence for positive selection on these alignments, neither with the branch-site test (supplementary table S15, Supplementary Material online) nor with the site test (supplementary table S14, Supplementary Material online).

22There was a significant enrichment for positively selected genes in the ant orthologs of D. melanogaster genes that were downregulated in 61-day-old flies compared with 10-day-old flies, based on a published microarray analysis (P = 0.011; below Bonferroni threshold α = 0.05/4 = 0.0125; supplementary table S16, Supplementary Material online).

23Two other genes known to be involved in aging were among the top-scoring genes for positive selection in our data set. The first was Tequila, which has been shown to be associated with aging in an experimental evolution study in D. melanogaster. The other was mitochondrial trifunctional protein α subunit, whose knock-out also reduces lifespan in D. melanogaster. Although not in the list of top hits, Sod2 (superoxide dismutase [Mn], mitochondrial), a gene known to have antioxidant activity and whose overexpression has been shown to be associated with increased lifespan in some strains of D. melanogaster, underwent positive selection at the base of the Hymenoptera lineage (FDR = 0.0073) and in the Acromyrmex echinatior branch (FDR = 9.6e-8).

24There was a marginally significant enrichment for positively selected genes among genes with biased expression in adult workers in S. invicta (effect size = 1.2, P = 0.025; not significant after Bonferroni correction α = 0.05/6 = 0.0083; supplementary table S17, Supplementary Material online) and a stronger enrichment for genes with queen-biased expression in adults (effect size = 1.8, P = 0.0028). Surprisingly, however, there was a pattern of weaker enrichment for genes with male-biased expression in adults (effect size = 1.04, P = 0.2381). At the pupal stage, we did not detect a significant enrichment for positively selected genes among any group of genes showing caste-biased expression. But similarly to the adult stage, the enrichment effect size was higher for genes with queen-biased expression (effect size = 1.2) than for genes with worker-biased expression (effect size = 1.1), and it was the lowest for genes showing male-biased expression (effect size = 1.06).

25In this article, we report results from a genome-wide scan for positive selection in protein-coding sequences of seven ant genomes, using the rigorous branch-site model of Codeml with stringent data quality control. Positive selection was detected in the ant lineage for 20% of the gene families analyzed. This proportion is similar in magnitude to the values observed in the other two insect lineages that we reanalyzed in this study: 20% in the 10 bee species and 36% in the 12 flies species.

26Our analysis identified similarities in patterns of positive selection between the ants and other insect lineages. Notably, at the broadest phylogenetic scale that our data sets allowed us to study, functional categories related to proteolysis, metabolism, response to stimuli, and immunity, were enriched for positively selected genes in ants, bees, and flies. Interestingly, studies in mammals, fishes, and urchins also provided evidence for positive selection on similar functional categories. Recurrent positive selection on such long evolutionary time scales is typical of genes involved in the interaction with changing environments or in conflict and competition, such as evolutionary arms races between sexes or between different species, which cause the perpetuation of adaptations and counter-adaptations in competing sets of coevolving genes. It is notable that positive selection patterns on these categories of genes do not seem to reflect or be strongly affected by the large life-history differences between lineages analyzed here, for example the evolution of eusociality in the hymenopteran lineages. In particular, our results on immunity-related genes challenge the hypothesis that hygienic behaviors in social insects could have relaxed the selective pressure on immune genes, since this should be reflected in reduced levels of positive selection on these genes.

27Our analysis indicated that genes involved in neurogenesis were under positive selection in ants and the primitively social lineages of bees. It was previously hypothesized that stronger selection on genes related to brain function and development should be observed in eusocial Hymenoptera species due to high cognitive demands associated with social life. However, our results are not consistent with this prediction because we also uncovered signs of positive selection at the base of the Hymenoptera lineage, i.e., before the evolution of sociality. Interestingly, a similar pattern had previously been reported with brain morphological data. A comparative analysis of insects showed that the size of mushroom body started to increase at the base of the Euhymenopteran (Orussioidea + Apocrita) lineage, approximately 90 My before the evolution of sociality in the Aculeata, and that there was no clear correlation between the size of brain components and the levels of sociality or cognition capabilities. To account for this observation, tried to identify factors, other than sociality, that may have placed unique selective pressure on brain evolution in species of the Hymenoptera lineage. Based on the observation that there was less positive selection on neurogenesis genes in highly social bees than in primitively social bees, they proposed that cognitive challenges might be associated with the mode of colony founding in social Hymenoptera. In particular, primitively social bees, which transit from a solitary phase during the process of colony founding to a social phase, could experience higher cognitive needs than highly social bees, which never go through a solitary phase. However, our results are also inconsistent with this model since increased positive selection was observed before the evolution of sociality in Hymenoptera. A comprehensive survey of positive selection on neurogenesis genes in Hymenoptera species, including species basal to the lineage, is required to identify precisely when the selective regime of these genes started to change, and in which hymenopteran sublineages it was maintained.

28Our results also challenge the hypothesis that genes involved in chemical signaling experienced increased positive selection in social insects. The analysis of olfactory receptor repertoires in two ants and a nonsocial wasp indicates that positive selection on amino-acid substitutions was surprisingly less frequent in ant than in wasp branches. Given the limited number of species used in this analysis, future work should concentrate on generating extensive annotation of olfactory receptors from more Hymenoptera as well as outgroup species to identify characters or traits that could be associated with the pattern of positive selection on olfactory receptors.

29Although our analyses did not provide support for previous hypotheses about the expected effect of sociality on gene evolution, we identified several interesting functional categories which were enriched for positively selected genes exclusively in the ant lineage, possibly reflecting ant-specific adaptations. The most consistent and robust result was that genes functioning in the mitochondria were particularly likely to be under positive selection. Mitochondrial activity plays an important role in the process of reproductive isolation and speciation, interactions with endosymbionts such as Wolbachia, diseases, and aging. In that respect it is notable that the evolution of sociality has been accompanied by a nearly 100-fold increase in lifespan of queens compared with their solitary ancestors. Three lines of evidence suggest that increased lifespan of queens might be related to increased positive selection on mitochondrial genes in the ant lineage.

30First, lifespan extension, not only in insects but also in other lineages such as birds and bats, appears to be associated with decreased production of reactive oxidative species (ROS). ROS are a normal by-product of cellular metabolism. In particular, one major contributor to oxidative damage is hydrogen peroxide (H2O2), which is produced from leaks of the respiratory chain in the mitochondria. Positive selection in ants on genes functioning in the mitochondria may thus reflect selection to increase mitochondrial efficiency and reduce ROS production. Interestingly, positive selection on genes with mitochondrial functions was previously documented in the bat lineage, which include species with exceptional longevity. In the bat Myotis lucifugus, ROS production was also shown to be significantly lower than in two similar sized mammal species (a mouse and a shrew) although the metabolic rates, and thus mitochondrial activity, of the former were much higher because of flight demands.

31Second, on the basis of gene expression data obtained in the fire ant S. invicta, our analyses revealed that positive selection was strongest on genes with queen-biased expression, intermediate on genes with worker-biased expression, and weakest on genes with male-biased expression. This association between levels of positive selection and caste-biased differences in gene expression cannot be simply accounted by differences in expression levels of mitochondrial genes (which are enriched for positively selected genes in ants) since in S. invicta mitochondrial genes are significantly less expressed in queens than in workers at the larval stage, and not differentially expressed at the adult stage (supplementary fig. S2, Supplementary Material online). The finding of higher levels of positive selection for genes more highly expressed in the castes with the longer lifespan (queens can live decades in some species, whereas workers have lifespan in the order of months, and males in the order of days) suggests that increased positive selection on queen-specific genes could be related to longer lifespan.

32Third, our analyses showed that the levels of positive selection were higher on orthologs of genes which are down-regulated during aging in flies. These genes include numerous energy metabolism genes, and their downregulation in old flies is thought to reflect a decline of normal and functional mitochondria with age. The finding of increased levels of positive selection on genes whose expression declines at older ages suggests that the function of these genes might be improved in ants, potentially delaying the loss of normal activity in mitochondria with age. It would be interesting to test if parallel mechanisms also evolved in the ant lineage to maintain the expression of these genes and delay the decline of mitochondria activity through lifespan in queens.

33In contrast to ants, there was no evidence of elevated levels of positive selection on mitochondrial functions in bees. As most social species, bees also evolved longer queen lifespans (more than 2 years) compared with males and workers (a few weeks). There are four possible explanations for the difference between ants and bees in the level of positive selection on mitochondrial genes. First, lifespan differences between castes are less pronounced in bees, where queens live up to 2–5 years, than in ants, where queens can live up to 30 years, possibly resulting in lower selective pressure to increase lifespan in bees than in ants. Second, because eusociality evolved independently in ants and bees it is possible that extended queen lifespans evolved by different molecular mechanisms. For example, vitellogenin may play a more central role for aging in bees than ants. Third, the evolution of mitochondria-related genes may have been differently constrained in ants and bees. For example, metabolic rates differ greatly between flying bee workers and non-flying ant workers because flight is an energetically costly behavior requiring highly elevated metabolic rates. Because metabolism and mitochondrial activity are closely connected, lower metabolic rates in ants might have alleviated functional constraints on mitochondria-related genes, allowing selection to act on lifespan extension. Fourth, the GC content in bee genomes was shown to be lower than in ant genomes. Some parts of the bee genomes, in particular their mitochondrial genomes, display extreme bias in nucleotide composition, which leads to significant effect on both the codon usage patterns and amino-acid composition of proteins and may have interfered with the action of positive selection.

34If positive selection acted to optimize the functioning of mitochondria in ants, it could be expected that the mitochondrial genome itself should be targeted by positive selection. However, mitochondrial genes generally exhibit very low dN/dS ratios and there was no clear evidence in our results for positive selection on the 13 genes of the mitochondrial genome itself. This suggests that innovations related to mitochondrial activity could arise more easily on nuclear genes, whereas mitochondrial genes seem more likely to maintain conserved core functionalities.

35In conclusion, this study provides a detailed analysis of the extent of positive selection events on protein-coding genes in seven ant species. Because false positives are a major concern for whole-genome scans for positive selection, we used a conservative methodology. We also reanalyzed data in bees and flies with the same methods to permit an unbiased and robust comparison of positive selection between lineages. The comparison between these three lineages provided interesting perspectives on the evolution of genes implicated in immunity, neurogenesis, and olfaction, and allowed us to pinpoint positive selection events that were specific to the ant lineage. In particular, we found that the evolution of extreme lifespan in ants was associated with positive selection on genes with mitochondrial functions, suggesting that a more efficient functioning of mitochondrial genes might have been an important step toward the extreme lifespan extension that characterizes this lineage. It would be interesting to complement this study by scans for genes under lineage-specific strong or relaxed purifying selection, to get a more global picture of natural selection patterns in ant genomes, and uncover additional genes that could have played a significant role during the evolution of the ant lineage.

36Protein-coding gene sequences of the seven ant genomes were downloaded from the Hymenoptera Genome Database (http://hymenopteragenome.org/ant_genomes/, last accessed April 24, 2014).

37The complete annotated gene sets were OGS_1.0 for Acromyrmex echinatior, OGS_1.2 for Atta cephalotes, OGS_2.2.3 for Solenopsis invicta, OGS_1.2 for Pogonomyrmex barbatus, OGS_3.3 for Camponotus floridanus, OSG_1.2 for Linepithema humile, and OGS_3.3 for Harpegnathos saltator. Coding sequences of five outgroup species were downloaded from the Hymenoptera Genome Database for the honey bee (Apis mellifera Amel_pre_release2) and the jewel wasp (Nasonia Vitripenis OGS_v1.2), from Flybase for the fruit fly (Drosophila melanogaster FB5.29), from BeetleBase for the flour beetle (Tribolium castaneum Tcas_3.0), and from vectorBase for the body louse (Pediculus humanus PhumU1.2).

38Gene families were obtained from a custom run of the OrthoDB pipeline for the Ant Genomic Consortium (http://cegg.unige.ch/orthodbants and http://bioinfo.unil.ch/supdata/Roux_positive_selection_ants/orthoDB_run.zip, last accessed April 24, 2014; pipeline of OrthoDB release 4). Briefly, OrthoDB implements a Best Reciprocal Hit clustering algorithm based on all-against-all Smith–Waterman protein sequence comparisons. The longest alternatively spliced form of genes is used. The orthologous groups are built at different taxonomic levels and it is possible to query for specific phyletic profiles by combining the criteria of absent, present, single-copy, multicopy, or no restriction, for each species within the studied clade.

39Gene families including strictly one ortholog in each of the 12 species were selected (2,756 gene families). Because annotations of the studied genomes are likely to be incomplete, families with a few missing genes—gene losses or unannotated genes—were included, with the restriction that at least four genes out of the seven ant species, and three genes out of the five outgroup species should be present in the gene family. have shown that among the seven ant species, there were generally few lost or missing genes, apart from S. invicta (less than 400 S. invicta genes were missing in single-copy orthologs families) and Ac. echinatior (<150 Ac. echinatior genes were missing in single-copy orthologs families). Our gene family selection criteria allow handling such a moderate amount of missing genes in families. In order to transfer functional annotations from D. melanogaster, only families including a fruit fly ortholog were retained. With these criteria, the number of OrthoDB groups in the data set increased to 4,337. All gene families were assumed to follow the species tree topology (fig. 1). The exclusion of families that experienced gene duplication facilitates the comparison of branches between gene families, and keeps our analysis from biases related to differential duplication among lineages and among genes, and to the consequences of duplication. Finally, results on single-copy orthologs can be easily compared with previously published studies using similar gene family topologies.

40Basic sequence quality features were first controlled as in. CDS (coding sequences) whose length was not a multiple of 3 or did not correspond to the length of the predicted protein, or that contained an internal stop codon, were eliminated; the longest CDS of genes showing multiple isoforms was retained; CDS shorter than 100 nt were eliminated.

41Because misalignment errors can be an important source of false positives in genome-wide scans for positive selection in coding sequences, we took great care at filtering the potentially problematic sites in the alignments. The quality filtering pipeline used here is adapted from the pipeline of the Selectome database release 4 (http://selectome.unil.ch, last accessed April 24, 2014). The multiple alignment of the protein sequences in each gene family was computed by M-Coffee from the T-Coffee package v8.93, which combines the output of different aligners. Similarly to Ensembl Compara (see http://www.ensembl.org/info/docs/compara/homology_method.html [last accessed April 24, 2014] for more details), four different aligners were used for M-Coffee (mafftgins_msa, muscle_msa, kalign_msa, and t_coffee_msa). M-Coffee outputs a consensus of four alignments from the different aligners, and a quality score for each residue based on the concordance of the alignment at each position by different aligners. Scores lie between 0, if a residue was not aligned at the same position by the different aligners, and 9 if it is reliably aligned at the same position in all cases. Reliably aligned residues with a score of 7 or above were retained. We used the heuristic algorithm of MaxAlign v1.1 to detect and remove sequences badly aligned as a whole (gap-rich sequences) in the multiple sequence alignments. When a sequence was removed, the gene family was realigned and refiltered using M-Coffee. Families left with less than four sequences were discarded because of insufficient power to detect positive selection. The protein alignments were reverse-translated to nucleotide alignments using the seq_reformat utility of the T-Coffee package.

42We used a stringent Gblocks filtering (v0.91b; type = codons; minimum length of a block = 4; no gaps allowed) to remove gap-rich regions from the alignments, as these are problematic for positive selection inference. The large memory requirements of M-Coffee for long alignments led us to use only Gblocks without M-Coffee scoring if the length of the alignment was greater than 9,000 nt.

43After filtering, our data set included 4,261 gene families with an average of 10.4 branches per family to test (fig. 1; 44,306 branches to test; median = 10 branches per family). The mean length of filtered alignment was 1,133 nt (median = 885 nt), ranging from a minimum of 54 nt to a maximum of 22,248 nt. Of note, lost or missing genes in families affect the topology of the trees and the possibility to compare equivalent branches of different families. In total, our data set contains 36,681 branches (83%) in 4,256 families which corresponded to the canonical topology defined by the species tree (fig. 1) and could be compared across families (e.g., table 2).

44Our analyses are likely to underestimate the genome-wide number of positive selection events because 1) single-copy orthologs tend to evolve under stronger purifying selection than multicopy gene families, 2) the ant genomes still lack good annotation of gene models and single-copy orthologs gene families could be missed, and 3) we filtered out unreliable parts of sequence alignments including fast evolving residues that are difficult to align. The last point is balanced by the fact that conserved regions might be more prone to positively selected substitutions and that the removal of unreliable regions seems to increase the power to detect positive selection.

45Another data set gathered all gene families from the OrthoDB database that could pass our quality filters, and notably families that experienced duplications. The CDS were filtered as described earlier. Amino-acid sequences were aligned using PAGAN version 0.47. The program GUIDANCE (v1.1) was used to assess alignment confidence and mask unreliably aligned residues. The combination of a phylogeny-aware aligner (PAGAN replaces PRANK and is based on the same principle) and of this filtering algorithm was shown to perform the best in recent benchmark studies on simulated data. Gene family phylogenies were built using RAxML (v7.2.9) from the amino-acid sequences, with the LG matrix and the CAT model. Amino-acid alignments were reverse-translated into the corresponding codon alignments. This resulted in 6,186 families tested, with an average of 11 genes, and an average length of filtered alignment of 3,129 nt (median of 2,385 nt, ranging from a minimum of 192 nt to a maximum of 20,556 nt).

46Contigs corresponding to mitochondrial genomes could be downloaded for five of the seven ant genomes (Ac. echinatior, At. cephalotes, S. invicta, P. barbatus, and L. humile). They were submitted to MITOS, a web server for the annotation of metazoan mitochondrial genomes (http://mitos.bioinf.uni-leipzig.de/index.py, last accessed April 24, 2014). This gave us the predicted coordinates of 13 mitochondrial protein-coding genes in these species. Frameshift errors or incomplete gene predictions were manually corrected. Mitochondrial genes from the outgroup species Ap. mellifera, N. Vitripenis, and T. castaneum were downloaded from GenBank (accession numbers: L06178; EU746609.1, and EU746613.1; AJ312413.2 and NC_003081.2, respectively). Mitochondrial genes from D. melanogaster were downloaded from Flybase at ftp://ftp.flybase.net/genomes/Drosophila_melanogaster/dmel_r5.43_FB2012_01/fasta/dmel-dmel_mitochondrion_genome-CDS-r5.43.fasta.gz (last accessed April 24, 2014). The alignment and filtering steps for the 13 mitochondrial gene families were identical to the data set of single-copy orthologs nuclear gene families (see above). A total of 119 branches were tested in this data set (average of 9.2 and median of 9 branches per family; average length of filtered alignment of 641 nt and median of = 621 nt, ranging from a minimum of 39 nt to a maximum of 1,413 nt).

47Single-copy ortholog gene families from the twelve sequenced Drosophila species were downloaded from ftp://ftp.flybase.net/12_species_analysis/clark_eisen/alignments/ (last accessed April 24, 2014) (files “all_species.guide_tree.longest.cds.tar.gz” and “all_species.guide_tree.longest.translation.tar.gz”). The alignment and filtering steps for these gene families were identical to the data set of single-copy ortholog gene families used for the ant analysis. Out of 6,698 initially downloaded Drosophila gene families, 3,749 (56%) passed our filters and could be tested for positive selection, resulting in 77,495 branches tested (average of 20.7 and median of 21 branches per family; average length of filtered alignment of 876 nt and median of 708 nt, ranging from a minimum of 15 nt to a maximum of 14,535 nt).

48Single-copy ortholog gene families from 10 bee species were downloaded from http://insectsociogenomics.illinois.edu/ (last accessed April 24, 2014). This set of gene families is incomplete as it is derived from the sequencing of expressed sequence tags (using 454 Life Science/Roche GS-FLX platform) from nine bee species, and from gene models of the honey bee Ap. mellifera. The alignment and filtering steps for these gene families were identical to the data set of single-copy ortholog gene families used for the ant analysis. Out of 3,647 initially downloaded gene families, 2,256 (62%) passed our filters and could be tested for positive selection, resulting in 20,169 branches tested (average of 8.9 and median of 9 branches per family; average length of filtered alignment of 611 nt and median of 528 nt, ranging from a minimum of 27 nt to a maximum of 3,945 nt).

49We used the updated branch-site test of Codeml from the package PAML v4.4c to detect Darwinian positive selection experienced by a gene family in a subset of sites in a specific branch of its phylogenetic tree. This test was previously used in genome-wide scans for positive selection in various lineages and is used by the Selectome project (http://selectome.unil.ch, last accessed April 24, 2014). It is acknowledged to be more sensitive for the detection of positive selection than branch tests or site tests, because it does not average the signal over all codons in the alignment (branch test) nor over all branches of the phylogeny (site test). It is also robust to relaxation of purifying selection (ω close to 1) since this scenario is accounted for in the null model. The alternative model is contrasted to the null model using a likelihood-ratio test (LRT), where log-likelihood ratios are compared to a chi-square distribution with 1 degree of freedom. Previous studies have reported the branch-site test to be conservative in this setup. We did not use the ω estimates to infer the strength of positive selection because they were shown to be unreliable.

50In the absence of a specific a priori hypothesis regarding which branches to test for positive selection, our implementation runs the test multiple times on each gene family, successively changing the branch selected as foreground. The branches considered as foreground are highlighted in red in figure 1. This approach was shown to be legitimate if P-values from the successive tests are corrected for multiple testing. We applied a FDR correction over all the P-values treated as one series (number of branches tested × number of gene families tested). In the ant single-copy orthologs nuclear data set, we analyzed a maximum of 15 branches leading to the 7 ant species, summing to 44,306 tests performed. In the ant mitochondrial data set, we analyzed a maximum of 11 branches leading to 5 ant species, summing to 119 tests (branches in red in supplementary fig. S3, Supplementary Material online). In the Drosophila single-copy orthologs data set, we analyzed a maximum of 21 branches, leading to a total of 77,495 tests (supplementary fig. S4, Supplementary Material online). Finally in the bee data set, we analyzed a maximum of 17 branches, leading to a total of 20,169 tests (supplementary fig. S5, Supplementary Material online).

51All computations were performed using Slimcodeml (release May 4, 2011), an optimized version of Codeml, based on the release 4.4c of the PAML package (downloadable at http://selectome.unil.ch/cgi-bin/download.cgi, last accessed April 24, 2014). Slimcodeml was estimated to run the branch-site models about 1.77 times faster than the original Codeml thanks to the use of external standard libraries for linear algebra calculations and specific optimizations for the computer architecture used. We verified on a subset of the gene families that the results given by Slimcodeml were identical with the original Codeml. Examples of Slimcodeml/Codeml control files used are provided in supplementary text, Supplementary Material online. For the ant mitochondrial data set, Codeml was used with the option “icode = 4” to use the Invertebrate mitochondrial genetic code (http://www.ncbi.nlm.nih.gov/Taxonomy/Utils/wprintgc.cgi#SG5, last accessed April 24, 2014).

52The branch-site model is known to display convergence problems in the calculation of likelihoods, leading to negative or artificially large log-likelihood ratios. We thus launched three independent runs for both the alternative and null hypotheses, for each branch of each gene family, and kept the best likelihood value of each run to calculate the log-likelihood ratio. Of note, the likelihood differences observed across the three runs were most of the time very small. Even after reconciliation of three replicate runs, we still observed a number of negative log-likelihood ratios (8% of the tests—most of them very close to 0). In such cases, we manually set the log-likelihood ratios to 0 (meaning nonsignificance). We recorded the largest differences in likelihood values between the three independent runs in both fixed and alternative models (d). The distribution of differences was bimodal, with a first major mode at d = 0, gathering most data, and a second minor mode at d ∼ 1. A cutoff at d = 0.004 clearly separated the two peaks. We used this stringent cutoff (d > 0.004) to eliminate all tests with potential convergence problems in the fixed and alternative models (see supplementary text and table S23, Supplementary Material online).

53Values of dN and dS were calculated with parameters extracted from Codeml results files (.mlc files).

54All calculations were performed on the SIB Vital-IT cluster in Lausanne (http://www.vital-it.ch/, last accessed April 24, 2014). All three runs and the two hypotheses of each test were performed on the same node of the cluster.

55The site test of Codeml from the package PAML v4.4e, allowing the dN/dS ratio (ω) to vary among sites, was run on the extensive data set of 6,186 families (see above). We contrasted the null model M8a (beta and ω with ω = 1) to the alternative model M8 (beta and ω with ω ≥ 1) with 11 site classes. Examples of Codeml control files used are provided in supplementary text, Supplementary Material online. Similar to the branch-site test, we launched three independent runs for both the alternative and null hypotheses for each gene family and kept the best likelihood value of each run for the LRT (supplementary table S19, Supplementary Material online). The likelihood ratios were compared to a chi-square distribution with 1 degree of freedom as recommended in PAML user’s guide (http://abacus.gene.ucl.ac.uk/software/pamlDOC.pdf, last accessed April 24, 2014).

56The program nhPhyml was used to estimate the G + C content at third codon positions at each node of the gene family trees (topology fixed, transition/transversion ratio estimated, alpha parameter estimated with eight categories). Following, we calculated the shift in GC3 content at each branch as the difference between GC3 contents at the nodes delimitating that branch.

57Olfactory receptors are difficult to process in automated pipelines since they are characterized by dynamic patterns of duplications and pseudogenization during evolution. Furthermore, the sequences of ORs are highly variable and notoriously difficult for automatic gene annotation. Accordingly, our main data set of single-copy orthologs was depleted in genes involved in olfaction (supplementary tables S7, S8, S10, and S11, Supplementary Material online) and GO categories related to olfaction could not be tested for enrichment of positively selected genes because they included too few annotated genes. We therefore used a more comprehensive data set of 873 manually annotated protein-coding sequences of OR genes (excluding suspected pseudogenes) provided by Hugh Robertson for P. barbatus (291 genes), L. humile (320 genes), and N. vitripennis (262 genes). Nucleotide sequences were translated and amino-acid sequences were aligned using MAFFT. Unreliably aligned residues were masked using GUIDANCE based on 32 bootstrap samples and a cutoff of 0.2 that was chosen so that the 15% lowest scoring residues are masked. Phylogeny was reconstructed using RAxML with the JTT substitution matrix, the CAT approximation, and 100 bootstrap samples. Because the resulting gene tree was too large for an analysis with the branch-site test of Codeml, we divided it into 16 smaller subtrees, each containing less than 100 leaves. Branches with as high as possible bootstrap support were chosen as splitting points. The 16 subtrees include all ant sequences but only 105 N. vitripennis sequences. The sequences from each subtree were realigned using PRANK version 100701 and reverse-translated into corresponding codon alignments. GUIDANCE was used to mask unreliably aligned codons (0.8 cutoff). Phylogeny was reconstructed using RAxML as above. Out of 1,744 branches in the initial tree, 1,400 branches from the subtrees were tested using the branch-site test of Codeml (see above), and the computation was successful (both null and alternative hypotheses) for 1,184 branches. Significant branches are highlighted in red in figure 3 and in supplementary figure S1, Supplementary Material online. Full results of the branch-site test on all 16 clades are shown in supplementary table S18, Supplementary Material online. A full tree with branch names and bootstrap values is provided as supplementary figure S1, Supplementary Material online. Newick trees of the 16 individual subtrees along with annotation of tested branches are available in supplementary text, Supplementary Material online.

58GO annotations for gene families were taken from the annotation of the D. melanogaster gene member they include (downloaded from http://flybase.org/static_pages/downloads/FB2011_02/go/gene_association.fb.gz, last accessed April 24, 2014). The annotation of children GO categories was propagated to their parent categories following the GO graph structure. GO categories mapped to 10 genes or less were discarded for the enrichment analysis.

59To identify over- and underrepresented functional categories present in the data sets used in this study, the package topGO version 2.4 of Bioconductor was used. A Fisher's exact test was used, combined with the “elim” algorithm of topGO, which decorrelates the graph structure of the GO to reduce nonindependence problems. The reference set was constituted of all OrthoDB families including a D. melanogaster gene with GO annotation. GO categories with an FDR < 20% are reported.

60Regarding the functional enrichment of genes targeted by positive selection, the Fisher's exact test approach has been criticized because it imposes the choice of an arbitrary cutoff to dichotomize genes into “significant” and “nonsignificant” categories. This leads to a loss of information and limits the power and robustness of this method. To test for GO functional categories for enrichment for positively selected genes, we instead used a gene set enrichment approach, which tests whether the distribution of scores of genes from a gene set differs from the whole data set scores distribution, allowing the detection of gene sets that contain many marginally significant genes. Different implementations for this approach have been proposed. The most widely used is the gene set enrichment analysis (GSEA), but it was shown to perform relatively poorly. Here, we used a SUMSTAT test: for a given gene set g including n genes, the SUMSTAT statistic is defined as the sum of scores of the n genes. This statistic was shown to be more sensitive than a panel of other methods, while controlling well for the rate of false positives. To be able to use the distribution of log-likelihood ratios of the positive selection test—which follows a chi-square distribution with 1 degree of freedom and spans several orders of magnitude—as scores in the SUMSTAT test, we applied a fourth root transformation as variance stabilizing method. This transformation conserves the ranks of gene families (see http://udel.edu/∼mcdonald/stattransform.html, last accessed April 24, 2014). According to the Central Limit Theorem, the distribution of SUMSTAT scores from random gene sets approaches a normal distribution whose mean and variance derives from the mean and variance of the scores of the complete list of tested genes G: and We performed bidirectional tests against this distribution to test whether the SUMSTAT statistic for a given gene set is higher or lower than expected by chance, corresponding to respectively enrichment or depletion for positively selected genes in this gene set. We verified the accuracy of this methodology by drawing an empirical null distribution for each gene set of size n found in the real data set, based on scores of 10,000 gene sets of same size n randomly picked from the whole data set. The distribution of SUMSTAT scores of these randomized gene sets approximates closely a normal distribution, even when the set size is small (supplementary fig. S6, Supplementary Material online). This makes the SUMSTAT test less computationally intensive than other gene set enrichment approaches (e.g., GSEA) where the null distribution cannot be inferred mathematically and randomizations have to be performed for each individual test. We verified that a GSEA approach gave broadly similar results (not shown).

61Because different gene sets sometimes share many genes in common, the list of significant gene sets resulting from enrichment tests is usually highly redundant. We implemented the “elim” algorithm from the Bioconductor package topGO, to decorrelate the graph structure of the GO. Briefly, the GO categories are tested recursively starting from the deeper levels of the GO tree, and the genes annotated to these significant categories are removed from all their parent categories. As the tests for different categories are not independent, it is not clear whether classical approaches to assess the FDR (e.g., ) are accurate. Thus, we calculated empirically an FDR at each P-value threshold by performing 100 randomizations where the scores of gene families were permuted and the gene set enrichment test rerun. The FDR is estimated as where at a given P-value threshold N0 represents the mean number of false positives obtained in the randomizations and Nt represents the number of positives obtained with the real data set. The FDR obtained with this approach was in good agreement with the Benjamini–Hochberg FDR. GO categories with an FDR < 20% are reported. Functional categories depleted in positive selection reflect the most conserved sets of functional categories, under the action of purifying selection. These are not discussed in this article.

62The gene set enrichment test ran on each individual branch of the tree with results of the branch-site test yields heterogeneous results, probably resulting from differences in power of the branch-site test on different branches of the phylogeny (supplementary table S9, Supplementary Material online; only branches Sinv, Pbar, Hsal, #3 and #6 show some significant categories at FDR 20%). This test could also be sensitive to false positive results of the branch-site test (e.g., GC-biased gene conversion, discussed in supplementary text, Supplementary Material online). Thus, we designed a test less sensitive to these problems. We considered a unique score per gene family reflecting the evidence for positive selection globally in the ant lineage, the mean of the branch-site test scores on the 13 individual ant branches. This scoring scheme should unveil functional categories of genes that experienced extensive and probably recurrent episodes of positive selection in the ant lineage, but is not strictly equivalent to using the results of a site test on ants branches, since it allows the detection of gene families with positive selection events affecting different sites on different branches. We also checked that in most cases, the enriched categories were not significant only because of a single outlier gene with a strong positive selection score, but displayed a significant shift in the distribution of positive selection scores of numerous genes (supplementary fig. S7, Supplementary Material online).

63Finally, as a sanity check, the gene set enrichment test was also performed using KEGG pathways annotation. KEGG pathways and the mapping to D. melanogaster genes were downloaded with the KEGG REST API (http://www.kegg.jp/kegg/rest/keggapi.html, last accessed April 24, 2014). Because hierarchical relationships among KEGG pathways are limited, we did not use the “elim” decorrelation algorithm. Pathways mapped to more than 10 genes were retained. In total, 51 KEGG pathways were tested.

64Mutant phenotype annotations of D. melanogaster genes were extracted from Flybase. The following ontologies were downloaded from the OBO foundry: the Flybase controlled vocabulary ontology (http://www.berkeleybop.org/ontologies/obo-all/flybase_vocab/flybase_vocab.obo, last accessed April 24, 2014), the Drosophila anatomical ontology (http://www.berkeleybop.org/ontologies/obo-all/fly_anatomy/fly_anatomy.obo, last accessed April 24, 2014), and the Drosophila developmental stages ontology (http://www.berkeleybop.org/ontologies/obo-all/fly_development/fly_development.obo, last accessed April 24, 2014). The relationships between genes and alleles, and between alleles and phenotypes (anatomical and developmental ontology categories) were extracted from Flybase (ftp://ftp.flybase.net/releases/FB2012_01/reporting-xml/FBgn.xml.gz, last accessed April 24, 2014; “derived_pheno_class” and “derived_pheno_manifest” entities). The information on gain or loss-of-function alleles was extracted from the file ftp://ftp.flybase.net/releases/FB2012_01/reporting-xml/FBal.xml.gz (last accessed April 24, 2014) (loss of function: controlled vocabulary term FBcv:0000287 and child terms; gain of function: FBcv:0000290 and child terms). The annotation of child phenotypic categories (anatomy of development) was propagated to their parent categories following the respective ontologies structures.

65To perform an enrichment analysis based on mutant phenotypes in fruit fly, we used the SUMSTAT test. Because the annotation is scarcer than the GO annotation, we used only the categories mapped to more than five genes for the enrichment analysis. The reported results include the annotation for gain and loss-of-function alleles. We observed very similar results when gain-of-function alleles were removed from the annotation (not shown).

66Microarray expression data from S. invicta were provided by the authors upon request. These included expression levels of clones of the spotted microarray used, as well as the list of genes identified to be overexpressed in each of the three castes (workers, queens, and males), both at pupal and adult stages. The mapping of clones to the gene model of S. invicta (OGS_2.2.3) was provided by Y. Wurm, and is similar to the mapping used in. If multiple clones mapped to the same gene, the average signal was used for expression. For differential expression, we used the results of the original study (BAGEL analysis, where a clone was considered to be differentially expressed between conditions if the Bayesian posterior probability was P < 0.001, corresponding to an FDR ∼ 5%). A gene was considered differentially expressed if at least one clone mapped to it was found differentially expressed. Expression data were available for 1,327 genes of the single-copy orthologs data set, including 603 genes overexpressed in at least one condition. We ran a SUMSTAT gene set enrichment test on the sets of genes with caste-specific expression (pupal male, pupal queen, pupal worker, adult male, adult queen, and adult worker). P-values were obtained by comparison to an empirical distribution created with 10,000 randomizations of gene scores.

67Aging and oxidative stress associated genes were obtained from a microarray study in D. melanogaster comparing the expression of genes in 10-day-old flies to 61-day-old flies, and flies exposed to 100% O2 for 7 days to controls. We tested the enrichment for positively selected genes (SUMSTAT test) in four gene sets constituted of up and downregulated genes in both contrasts. P-values were obtained by comparison to an empirical distribution created with 10,000 randomizations of gene scores.

68Genes with mitochondrial function were identified as those mapped to any of the 310 GO categories including “mitochondria*” in their names or synonym names (using the search engine on http://amigo.geneontology.org/, last accessed April 24, 2014). Three hundred and thirteen of the identified genes had available microarray expression data in S. invicta.

69Raw and filtered alignments used in these analyses track files for the alignment editor Jalview, Codeml control files and result files can be downloaded at http://bioinfo.unil.ch/supdata/Roux_positive_selection_ants/Roux_et_al_datasets.tar.gz (last accessed April 24, 2014).

70A simple web interface displaying gene families, GO mapping, Codeml results, and alignments (through a Jalview applet) is available at http://bioinfo.unil.ch/supdata/Roux_positive_selection_ants/families.html (last accessed April 24, 2014). Jalview tracks display the regions used or filtered out in the original protein alignments, as well as the residues found to be under positive selection by Bayes Empirical Bayes in all the branches tested for each of the three replicate runs (fig. 2).

71Supplementary tables S1–S26, figures S1–S7, and supplementary text are available at Molecular Biology and Evolution online (http://www.mbe.oxfordjournals.org/).