The genomic basis of army ant chemosensory adaptations.

The evolution of mass raiding has allowed army ants to become dominant arthropod predators in the tropics. Although a century of research has led to many discoveries about behavioural, morphological and physiological adaptations in army ants, almost nothing is known about the molecular basis of army ant biology. Here we report the genome of the iconic New World army ant Eciton burchellii, and show that it is unusually compact, with a reduced gene complement relative to other ants. In contrast to this overall reduction, a particular gene subfamily (9-exon ORs) expressed predominantly in female antennae is expanded. This subfamily has previously been linked to the recognition of hydrocarbons, key olfactory cues used in insect communication and prey discrimination. Confocal microscopy of the brain showed a corresponding expansion in a putative hydrocarbon response centre within the antennal lobe, while scanning electron microscopy of the antenna revealed a particularly high density of hydrocarbon-sensitive sensory hairs. E. burchellii shares these features with its predatory and more cryptic relative, the clonal raider ant. By integrating genomic, transcriptomic and anatomical analyses in a comparative context, our work thus provides evidence that army ants and their relatives possess a suite of modifications in the chemosensory system that may be involved in behavioural coordination and prey selection during social predation. It also lays the groundwork for future studies of army ant biology at the molecular level.

INTRODUCTION

Army ants are dominant arthropod predators in tropical ecosystems worldwide (Brady, 2003; Kronauer, 2009). It has been estimated that a given square metre of rain forest is subject to an army ant raid more than once a day on average (Kaspari & O’Donnell, 2003), with some species depleting ~25% of invertebrate biomass in such events (Kaspari et al., 2011). The key to the ecological success of army ants is mass raiding, a form of obligate collective foraging in which colonies attack and overwhelm live prey (Gotwald, 1995; Kronauer, 2020; Rettenmeyer, 1963; Schneirla, 1971). Army ant mass raids contain up to hundreds of thousands of individuals collectively searching for and exploiting food resources (Brady, 2003; Gotwald, 1995; Kronauer, 2009; Schneirla, 1971). This behaviour has evolved at least twice within the ant subfamily Dorylinae, giving rise to separate Old World and New World army ant lineages that display striking convergence in behaviour and morphology (Borowiec, 2017). A suite of other traits have evolved in parallel with mass raiding, including nomadism; morphologically specialized, wingless queens; and reproduction via colony fission (Brady, 2003; Gotwald, 1995; Kronauer, 2009; Schneirla, 1971). Together, these traits are referred to as the “army ant adaptive syndrome” (Gotwald, 1982). Both army ant lineages display elevated rates of diversification (Borowiec, 2017), demonstrating the evolutionary success of this lifestyle. Army ants have proven excellent models for studying social predation, and several behavioural and physiological adaptations for this behaviour have been described (Brady, 2003; Bulova et al., 2016; Gotwald, 1995; Kronauer, 2009; Rettenmeyer, 1963; Schneirla, 1971). However, molecular adaptations facilitating social predation—a premise for the ecological dominance of army ants in tropical rainforests—remain poorly understood, and no genomic studies of army ants have been conducted.

Like most other ants (Wilson, 1965), army ants have evolved a sophisticated pheromone communication system to coordinate social behaviour, including navigation, recruitment of nestmates to prey sources, alerting nestmates of danger, and coordinating defence of the queen and care for the brood (Brown, 1959; Chadab & Rettenmeyer, 1975; Schneirla, 1971; Torgerson & Akre, 1970). The chemosensory system is responsible for recognizing these pheromones (Dumpert, 1972; Ozaki et al., 2005; Pask et al., 2017; Sharma et al., 2015; Slone et al., 2017; Trible et al., 2017; Yan et al., 2017). Cuticular hydrocarbons (CHCs), a diverse class of chemicals that are displayed on the surface of arthropods, play a particularly prominent role in ant communication by conveying species, colony and caste identity (Liebig, 2010; van Zweden & D’Ettorre, 2010). Chemosensation is also important for ant foraging by facilitating prey detection and discrimination (Lang & Menzel, 2011; Manubay & Powell, 2020; Mirenda et al., 1980; Witte et al., 2010), and here again, CHCs play an important role (Lang & Menzel, 2011; Manubay & Powell, 2020; Witte et al., 2010). While CHCs and other army ant pheromones have been studied in some detail (e.g., Bagneres et al., 1991; von Beeren et al., 2018; Brückner et al., 2018; Keegans et al., 1993), nothing is known about the army ant chemosensory system and how it compares with that of other ants. Genomic and neuroanatomical studies in the nonarmy ant doryline Ooceraea biroi, which possesses a subset of army ant traits (Kronauer, 2009), have revealed an especially elaborate chemosensory system, with several features that have so far not been observed in other ants (McKenzie et al., 2016; McKenzie & Kronauer, 2018; Oxley et al., 2014; Trible et al., 2017). Chemosensory receptor genes in the 9‐exon odorant receptor (OR) gene family, which detect CHCs, are particularly expanded in O. biroi, as are associated neural structures in the brain (T6 glomeruli in the antennal lobe; McKenzie et al., 2016). Furthermore, workers of this species possess a unique arrangement of CHC‐detecting sense hairs (basiconic sensilla; McKenzie et al., 2016). Specifically, in O. biroi, basiconic sensilla occur in an unusually high density along the ventral (with respect to the body axis while antennae are held naturally) surface of the female antennal club, but nowhere else on the antenna. Given the paucity of molecular or neurobiological data within the Dorylinae, it is unknown whether any of these features are conserved across dorylines and whether they might be linked to the evolution of army ant‐like lifestyles.

The conspicuous Eciton burchellii (Figure 1a)—an abundant and ecologically impactful species in Neotropical forests (Kaspari et al., 2011; O’Donnell et al., 2007)—is perhaps the best studied army ant (Gotwald, 1995; Rettenmeyer, 1963; Schneirla, 1971). Like other army ants, E. burchellii is nomadic, reproduces via colony fission, possesses morphologically specialized queens and forages via mass raiding (Gotwald, 1995; Rettenmeyer, 1963; Schneirla, 1971). E. burchellii stands out among New World army ants in being an epigaeic swarm raider, meaning that foragers fan out across the surface of the forest floor into a large “raiding carpet” (Schneirla, 1971). This species consumes more nonant prey than most other New World army ant species: ~50% by mass according to Franks (1982), 7% of prey items according to Hoenle et al. (2019). However, the majority of E. burchellii prey items belong to a few species of Camponotus carpenter ants, indicating that, like most other army ants, E. burchellii is still a specialized ant predator (Hoenle et al., 2019).

FIGURE 1

The genome of Eciton burchellii. (a) An E. burchellii major (photo by D.J.C.K.). (b) Phylogeny and genome assembly sizes of representative ant species. Phylogeny after Moreau and Bell (2013). Genus abbreviations are as follows: D. = Dinoponera, H. = Harpegnathos, O. = Ooceraea, E. = Eciton, Li. = Linepithema, Ps. = Pseudomyrmex, La. = Lasius, Cam. = Camponotus, Po. = Pogonomyrmex, S. = Solenopsis, Car. = Cardiocondyla, At. = Atta, Ac. = Acromyrmex. (c) Normalized number of N. vitripennis gene matches in the genomes of representative ant species, as calculated by the Blat algorithm with a 50% alignment coverage cutoff, and relativized to genome completeness as estimated by busco. Data from E. burchellii are in orange

To investigate the genomic and chemosensory basis of predation in E. burchellii, we sequenced the species’ genome and worker and male antennal and whole‐body transcriptomes, and investigated worker and male antennal and antennal lobe anatomy. Compared with other ants, we found that E. burchellii has a compact genome, with a particular paucity of interspersed repeats. E. burchellii also shows a general reduction in gene number, but no large, significant contractions in any specific gene family. Correspondingly, E. burchellii possesses relatively few chemosensory genes, with reduced repertoires of most OR gene subfamilies. However, as in O. biroi, the hydrocarbon‐sensitive 9‐exon OR subfamily is expanded. Gene evolution analysis revealed that these expansions happened independently in both species, indicating parallel selection pressures on OR genes. In addition, we detected an unusually high number of macroglomeruli in the antennal lobes, which may play a role in the recognition of nonhydrocarbon pheromones. Finally, we observed striking similarities in the peripheral and central chemosensory systems of E. burchellii and O. biroi. Integrating all data, our study suggests that E. burchellii has evolved specific adaptations in the olfactory system, potentially for processing semiochemicals in the contexts of predation and social information transfer. This suggests that the chemosensory system has played an important role in the evolution of army ants and other dorylines, and lays a foundation for studying the molecular basis for social predation.

METHODS

Genome sequencing and assembly

From a single colony, Eciton burchellii (subspecies E. burchellii foreli) workers and males were collected at La Selva Biological Station in Heredia Province, Costa Rica, frozen alive at −80°C, and split into two packages containing dry ice (see RNA sequencing below). No queens were collected due to the difficulty of sampling army ant queens and the fact that queen removal destroys the whole colony. One package contained only males and was transported to the University of Wisconsin‐Madison, Madison, WI, USA. Genomic DNA was extracted from one male using a modified version of the Genomic‐tip DNA extraction protocol for mosquitoes and other insects (Qiagen) as previously described (Suen et al., 2011). Three libraries for Illumina sequencing were constructed, including a 3‐ and an 8‐kbp insertion mate‐pair library, and sequenced at the University of Wisconsin‐Madison Biotechnology Center on two independent 2×300‐bp paired‐end Illumina MiSeq runs, yielding 10 and 16 million reads, respectively. A third shotgun library was sequenced on an Illumina HiSeq 2500 with 2×150‐bp paired‐end reads to yield 157 million reads. Reads were trimmed and prepared with nextclip (Leggett et al., 2014), and assembled using allpaths‐lg (Gnerre et al., 2011). The final assembly contained 190,515,515 bp on 2,492 contigs, scaffolded onto 449 scaffolds for a final length of 192,351,860 bp including gaps. The sequencing data and genome assemblies are available at the NCBI SRA and genome databases (Bioproject PRJNA632625).

RNA sequencing

The second package contained workers and males and was transported to the Rockefeller University. The antennae of 10 workers, spanning a range of worker sizes including representatives of each size class (minims, media, submajors and majors), were separated from bodies at 4°C. All antennae and all remaining bodies were then transferred separately into two vials containing Trizol reagent (Thermo Fisher Scientific), and subsequently the two samples were homogenized using a TissueLyser II (Qiagen). Five males were processed in the same way, except that only one of the antennae‐less bodies was used due to the exceptionally large size of E. burchellii males. Chloroform was then added to the Trizol, and the aqueous phase was retrieved using Phase Lock columns (5Prime) and cleaned using Qiagen RNeasy columns following the manufacturer's protocol. Illumina cDNA libraries were prepared at the Rockefeller University Genome Resource Center and sequenced on an Illumina HiSeq 2500. Reads were aligned to the assembled genome with star (Dobin et al., 2013), and gene expression was calculated using cufflinks (Trapnell et al., 2010).

Manual annotation of chemosensory genes

Odorant receptor genes, ionotropic receptor genes, gustatory receptor genes, odorant binding protein genes and chemosensory protein genes were annotated manually, largely following Oxley et al. (2014). Briefly, gene loci were located using tblastn (Altschul et al., 1990) and exonerate (Slater & Birney, 2005), with queries from Ooceraea biroi, Solenopsis invicta and Pogonomyrmex barbatus. Alignments from tblastn (Altschul et al., 1990) and exonerate (Slater & Birney, 2005) were then pieced together manually in the Apollo genome annotation editor (Lee et al., 2009), using RNA‐sequencing (RNA‐seq) data and muscle (Edgar, 2004) alignment to query proteins to generate complete annotations in cases where blast and exonerate alignments were insufficient. Gene models encoding over 200 amino acids with premature stop codons, frame shifts, missing splice sites, or truncated or missing exons (which could not be explained by an assembly gap) were designated as pseudogenes, while gene models encoding over 200 amino acids with no such issues were designated as putatively functional.

Automated gene annotation

The ab‐initio gene predictors augustus (Stanke et al., 2006) and genemark (Lomsadze et al., 2005) were trained within the BRAKER1 pipeline (Hoff et al., 2016) using all RNA‐seq reads aligned to the genome with tophat. trinity (Haas et al., 2013) was used to assemble RNA‐seq reads into putative transcripts, using the “jaccard clip” option to avoid gene fusions from overlapping untranslated regions (UTRs). The ab‐initio gene predictor snap (Korf, 2004) was trained with genes predicted by the CEGMA pipeline (Parra et al., 2007) as well as with trinity‐assembled transcripts which met all of the following criteria: (i) complete open reading frames, (ii) a blastp alignment to a Drosophila melanogaster gene (FlyBase release 6) which was at least 95% of the length of both the target and query proteins, and (iii) exonerate alignment to the genome with at least 95% of the possible alignment score. augustus, genemark and snap were then run in the MAKER2 pipeline (Holt & Yandell, 2011). Protein evidence for MAKER2 consisted of official gene sets for D. melanogaster and Apis mellifera, and NCBI refseq predictions for Acromyrmex echinatior and Ooceraea biroi. blastx evidence was excluded as it led to numerous fusions in tandem arrays. “EST” evidence consisted of trinity‐assembled RNA‐seq reads as well as tophat junctions, and the “avoid EST fusion” option was used because many UTRs overlapped. maker repeat masking (repeat masker and repeat runner programs, www.repeatmasker.org) was used with custom‐trained repeat libraries constructed by repeat modeler (www.repeatmasker.org). Final gene models from MAKER2, as well as all ab‐initio predictions from both MAKER2 and BRAKER1 and alignment‐based models from CEGMA, were then input into evidence modeler (Haas et al., 2008), along with exonerate‐aligned proteins and trinity‐assembled RNA‐seq evidence (as used in MAKER2) and PASA‐assembled transcript evidence from unguided and genome‐guided trinity RNA‐seq assemblies. Gene models with no significant blast hits in D. melanogaster or O. biroi, no expression in our RNA‐seq data, and no PFAM‐A (Finn et al., 2010) nonretrotransposon domain (as determined by interpro scan; Zdobnov & Apweiler, 2001) were then excluded, manual annotations of chemosensory genes were added, and exons from the evidence modeler gene set that overlapped manual annotations were deleted.

Comparative genomics

Genomes from 12 additional ant species (Bonasio et al., 2010; Konorov et al., 2017; Nygaard et al., 2011; Oxley et al., 2014; Patalano et al., 2015; Rubin & Moreau, 2016; Schrader et al., 2014; C.D. Smith et al., 2011a; Smith et al., 2011b; Suen et al., 2011; Wurm et al., 2011) were downloaded from the Hymenoptera Genomes Database (Elsik et al., 2016) or the NCBI Assembly database. Genome assembly statistics were calculated using quast (Gurevich et al., 2013). Genome repeat content, including transposable element content, was calculated for each genome using the repeatmodeler2 pipeline (Flynn et al., 2020). To obtain an approximation of gene complement independent of annotation effort/quality, we used blat (Kent, 2002) to align the NCBI RefSeq gene set (version 101) from Nasonia vitripennis to each ant genome, using a 50% alignment coverage cutoff to ensure that genes were not counted twice. This gene set was chosen because N. vitripennis is an outgroup to all ants, and therefore equally related to all ants, and because the N. vitripennis genome is exceptionally well annotated. Apis mellifera also matches these criteria and is more closely related to ants. However, it possesses a slightly reduced gene complement (Elsik et al., 2014). The busco program (version 3.0.2) (Simão et al., 2015) was used to determine genome completeness for all ant species using the Hymenopteran ODB9 gene set, and the number of N. vitripennis gene blat hits was relativized to the percentage of complete genes recovered by busco to control for genome completeness effects on estimated gene content. For analyses of protein family and pathway evolution, protein sequences for all species were functionally annotated with GO terms, PFAM domains and KEGG pathways using eggnog‐mapper (Huerta‐Cepas et al., 2017). ℓ1ou (Khabbazian, Kriebel, Rohe, & An′e, 2016) was used to detect significant shifts in gene repertoire evolution for total gene number, GO Slim terms, PFAM domains, KEGG pathways and OR subfamilies (see below) in an explicitly phylogenetic framework. ℓ1ou estimates the number and position of shifts in phenotypic optima of quantitative traits (e.g., gene number) along phylogenetic trees using Ornstein–Uhlenbeck (OU) models. Statistical support for proposed shifts was calculated using likelihood ratio tests against models where the phenotypic optimum was not allowed to shift along a specified branch (using the “candid.edges” option of the “estimate_shift_configuration” function).

Odorant receptor evolutionary analyses

Amino acid sequences for ORs from eight additional ant species and two outgroup hymenopterans were obtained from other studies (McKenzie et al., 2016; Oxley et al., 2014; Robertson et al., 2010; C.D. Smith et al., 2011a; Smith et al., 2011b; Zhou et al., 2012) and aligned using mafft with the linsi parameter set (Katoh et al., 2005). raxml (Stamatakis, 2014) was used to reconstruct a maximum‐likelihood phylogeny for these genes under the LG +gamma model of protein evolution (Le & Gascuel, 2008), and the raxml rapid bootstrap algorithm (Stamatakis et al., 2008) was used to assess node support. The LG +gamma model was chosen as it was consistently recovered as the best or second best (after JTT +gamma) model in our previous OR phylogenetic analyses (McKenzie et al., 2016; McKenzie & Kronauer, 2018). For gene birth and death analysis, this tree was imported into notung (Chen et al., 2000) and reconciled with a species tree following the topology of Moreau et al. (2006), Brady et al. (2006) and Borowiec et al. (2017), and the divergence times from Borowiec et al. (2017), with an edge weight threshold of 70% bootstrap support (i.e., notung was allowed to rearrange nodes with lower than 70% bootstrap support when it led to a more parsimonious gene‐tree/species‐tree reconciliation). notung was then used to count gene duplications and losses along each branch of the species tree. We calculated birth and death rates on a per‐gene basis along each branch in the species tree as Log2((ancestral repertoire +duplication events or loss events)/ancestral repertoire)/branch length, which represents the number of times a gene family doubled or halved per unit time. Net expansion was calculated similarly as Log2((ancestral repertoire +duplication ‐ loss events)/ancestral repertoire)/branch length. Positive selection analysis was performed with the codeml algorithm of the paml package (Yang, 2007), with coding sequences aligned based on their translated peptide sequences using the mafft aligner with the ginsi parameter set (Katoh et al., 2005). codeml site models 1a and 2a, 7 and 8, and 8a and 8 were compared using likelihood ratio tests.

Neuroanatomical analyses

Additional E. burchellii workers and males were collected at La Selva Biological Station and preserved in 70% ethanol. Brains were later dissected from one male and one medium‐sized worker in 70% ethanol and further dehydrated in 80%, 90%, 95% and 3 × 100% ethanol. Brains were then cleared and mounted in methyl salicylate. Imaging was performed on a Zeiss LSM 880 laser scanning confocal microscope with 488‐nm laser stimulation used to image tissue autofluorescence. We found that ethanol fixation and preservation was sufficient to preserve antennal lobe morphology, and in fact, this method performed better than immersing whole ants in 4% PFA for autofluorescence imaging (data not shown). Glomeruli were reconstructed manually using the segmentation editor plugin in the Fiji distribution of imagej (Schindelin et al., 2012), and assigned to input tracts following Zube et al. (2008). Potentially due to the general decrease in glomerulus number outside of the T6 cluster, glomeruli innervated by tracts T1, T2 and T4 could not be unambiguously distinguished. These input tracts were therefore not separated in our analysis. Glomerulus volume was calculated using the 3D analysis plugin in the 3D imagej suite (Ollion et al., 2013). Glomerulus radius was calculated following Kelber et al. (2009) as the radius of a sphere with equal volume to the glomerulus in question. We found that this transformation (essentially a cube‐root transformation) yielded more normal distributions compared to raw volume measurements. Putative macroglomeruli were identified using two criteria. First, following Kelber et al. (2009), we identified glomeruli with radii at least five standard deviations larger than the median glomerulus radius. Because this cutoff is arbitrary, we also tested whether larger glomeruli belonged to a size distribution separate from the remaining glomeruli by fitting statistical distributions to our measures of glomerular radii using finite Gaussian mixture modelling with an expectation maximization algorithm implemented in the R package mclust (Scrucca et al., 2016). The Bayesian information criterion (BIC) was used to determine the best fitting model parameters, including the number of distribution components/clusters. A Shapiro–Wilk's test was used to assess whether the clusters identified by mclust matched the Gaussian distribution assumed by the model.

Scanning electron microscopy

Eciton burchellii workers were fixed immediately after collection in Karnovsky solution (2.5% glutaraldehyde and 2% paraformaldehyde, 0.1 m phosphate buffer, pH 7.4) and post‐fixed in 1% osmium tetraoxide for 1 hr. Samples were dehydrated in a gradient of ethanol (30, 50, 70, 80, 90, 95, 100%), critical point dried with CO2 (Leica EM CDP300, critical point dryer), and coated with gold (Eiko IB5 ion coater). Male specimens stored at −80°C were washed in 50% DMSO overnight, and then washed in distilled water at least 10 times prior to fixation. All samples were imaged with an Hitachi S‐3700N scanning electron microscope.

Statistics

All other statistics were calculated in R or with the stats package from SciPy (R Core Team, 2013; Virtanen et al., 2019). For Pearson's correlation analysis, normality of residuals was assessed using the Shapiro–Wilk's test for normality.

RESULTS

Assembly of 157 million paired‐end 150‐bp Illumina HiSeq shotgun reads and 26 million 300‐bp Illumina MiSeq mate‐pair reads (with 3‐ and 8‐kbp inserts) yielded a 192‐Mbp draft genome. Contig N50 was 187 kbp, while scaffold N50 was 2,433 kbp. This genome is relatively small compared to other ants, but within the range of variation (Figure 1b; Table S1). busco (Simão et al., 2015) analysis of genome completeness showed that 97.1% of all widely conserved hymenopteran genes were present in this draft assembly, making this one of the most complete ant genomes reported to date (Table S1).

Our annotation procedure yielded a final set of 12,670 genes, including 499 manually annotated chemosensory genes. Of these genes, 6,941 had a blastp alignment to a Drosophila melanogaster protein (FlyBase release 5) that spanned >50% of the length of the D. melanogaster sequence. The total number of identified genes is small relative to the other ant genomes we examined with similar annotation strategies (Lasius niger and Cardiocondyla obscurior), and the number of D. melanogaster blastp hits is likewise small compared with most of the other ant genomes we examined (Table S1). Eciton burchellii also had the fewest Nasonia vitripennis blat matches against the genome itself relative to busco completeness (Figure 1c; Table S1), indicating that the reduced gene number we observed is due to a true reduction in gene content, rather than an annotation artefact. However, phylogenetic analysis of trait shifts using OU models implemented in ℓ1ou (Khabbazian et al., 2016) did not detect a statistically significant shift in total gene repertoire size along the branch leading to E. burchellii (Figure S1a). Gene count as approximated by our methods showed only a weak correlation with genome size (Pearson's correlation, R 2 =.32, p =.045), even though E. burchellii is on the low end for both. To elucidate the evolutionary forces leading to gene repertoire reduction in E. burchellii, we used ℓ1ou to test for significant repertoire size shifts for gene GO terms, PFAM domains and KEGG pathways along the phylogenetic tree branch leading to E. burchellii. Two GO terms showed significant expansion in the E. burchellii lineage: GO:0030705 (cytoskeleton‐dependent intracellular transport; false discovery rate [FDR]‐adjusted p =.042, Figure S1b) and GO:0003729 (mRNA binding; FDR‐adjusted p =.042, Figure S1c). However, both of these GO terms also show expansion in the L. niger and C. obscurior lineages (Figure S1b,c). All three of these species were annotated with in‐house pipelines using similar ab‐initio predictors (Table S1), while all other annotations used were from the NCBI RefSeq pipeline. We therefore suspect that these expansions are technical artefacts rather than biologically meaningful. One PFAM domain (PF00412; LIM domain) also showed trait deviation in the three maker‐annotated lineages (expanded in E. burchellii [FDR‐adjusted p =.019] and in C. obscurior, contracted in L. niger; Figure S1d), while another (PF07524; Bromodomain associated) showed significant contraction in E. burchellii (FDR‐adjusted p =.0011) potentially unrelated to annotation method (Figure S1e). The contraction observed in PF07524 was rather small, however, from a median of 11 genes in most ants to seven in E. burchellii. No significant shifts in the E. burchellii lineage were found for any KEGG pathway. Together, our analyses failed to recover any reduction in gene functional categories which might be driving the overall trend of reduced gene numbers in E. burchellii.

Analysis of repeat content showed fairly typical amounts of LINE and LTR element transposons and all types of clustered repeats in E. burchellii, but notable reductions in SINE and DNA element transposons as well as unclassified interspersed repeats (Table S1). As a consequence, E. burchellii has the second lowest repeat content of any ant examined in this study, second only to the only ant with a smaller genome, C. obscurior (Table S1). In fact, we found that total repeat content was highly correlated with genome size across all species (R 2 =.95, p <.00001), suggesting that this is a major driver of genome size evolution in ants. This correlation held for every subtype of repeat except for simple repeats and low‐complexity sequences (FDR‐adjusted p <.05 for SINEs, LINEs, LTR elements, DNA elements, unclassified interspersed repeats and satellites). This contrasts with the findings of Schrader et al. (2014), who found the opposite pattern in their analysis of eight species using a custom TE annotation pipeline.

Our manual annotation of chemosensory genes in E. burchellii yielded 318 intact, putatively functional OR gene models, as well as 12 CSPs, 13 OBPs, 20 GRs and 23 IRs (Dataset S1). We also recovered 109, two and two putative pseudogenes in the OR, GR and IR families, respectively. The numbers of CSPs, OBPs and IRs all fall into the typical ranges for ants (Table S1). The number of GRs is lower than in many ants, but it is similar to the repertoire size seen in Harpegnathos saltator and O. biroi, which along with E. burchellii belong to phylogenetic lineages (the poneroids [H. saltator] and dorylines [O. biroi and E. burchellii]) that branched off before the other ants with high‐quality chemosensory gene annotations (Borowiec et al., 2017).

Strikingly, E. burchellii possesses one of the smallest putatively functional OR repertoires known for an ant (Figure 2a,b). This repertoire reduction is primarily due to high rates of gene death in the ancestral lineage leading to the two sequenced dorylines, E. burchellii and O. biroi, followed by relatively high rates of gene death and lower rates of gene birth in the E. burchellii lineage compared to the O. biroi lineage (Figure 2a; Figure S2a,b, Datasets S2–S4). In contrast to this overall trend, the 9‐exon subfamily of ORs is expanded in E. burchellii compared to most other ants (Figure 2b). Phylostratigraphy showed elevated rates of gene duplication in 9‐exon ORs in the ancestral lineage leading to E. burchellii and O. biroi, as well as along the branches leading to each species (Figure 2c; Figure S2c). The high rate of 9‐exon OR gene death in the ancestral lineage indicates that most extant gene expansions are specific to either E. burchellii or O. biroi (Figure S2d, Dataset S5). The overall net reduction in OR repertoire in the common ancestor of E. burchellii and O. biroi is driven by a stark net reduction in non‐9‐exon ORs, due to high rates of gene death (Figure 2d; Figure S2e,f). To determine whether OR repertoire size evolution in E. burchellii and O. biroi relative to other ants represented a significant evolutionary shift vs. random trait drift, we tested for trait optimality shifts along the tree using ℓ1ou for all ORs and for the five largest ant OR subfamilies (9‐exon, L, V, U and E). We found a significant expansion in total OR number along the branch leading to the ants, but none within the ants (Figure S3a). Along the branch leading to E. burchellii and O. biroi, we observed a significant expansion in 9‐exon ORs and a significant contraction in V ORs (p =.02 for both subfamilies; Figure S3b,c), with no significant changes in the other large OR subfamilies.

FIGURE 2

Evolution of odorant receptor gene repertoires in Hymenoptera. (a) Phylogeny of 11 species of Hymenoptera with branches coloured according to the rate of net change in OR gene repertoire size. Species abbreviations are as follows: Nvit: Nasonia vitripennis, Amel: Apis mellifera, Hsal: Harpegnathos saltator, Obir: Ooceraea biroi, Ebur: Eciton burchellii, Lhum: Linepithema humile, Cflo: Camponotus floridanus, Pbar: Pogonomyrmex barbatus, Sinv: Solenopsis invicta, Acep: Atta cephalotes, Aech: Acromyrmex echinatior. (b) Number of intact ORs in 11 species of Hymenoptera. OR numbers are broken down by subfamily for the five largest odorant receptor subfamilies. (c) Hymenoptera phylogeny with branches coloured according to the rate of net change in 9‐exon OR gene repertoire size. (d) Hymenoptera phylogeny with branches coloured according to the rate of net change in OR gene repertoire size excluding 9‐exon ORs

To determine whether misannotation of functional genes as pseudogenes due to sequencing errors or stop codon readthrough (e.g., Prieto‐Godino et al., 2016) could be influencing our results, we characterized the expression, OR subfamily and number of putative inactivating features (frame shifts, splice site mutations or premature stop codons) of all putative pseudogenes (Table S2). Fifteen pseudogene models had a single putative inactivating feature, five pseudogene models had two putative inactivating features and two pseudogene models had three putative inactivating features. The remaining 87 pseudogene models were classified as “extensively pseudogenized” based on having four or more putative inactivating features or else having missing exons that could not be explained by an assembly gap. Of the 22 genes with three or fewer putative inactivating features, only 11 showed expression at >2 FPKM in any sex/tissue. Ten of these 11 potential false‐pseudogene models were 9‐exon subfamily ORs.

To investigate evolutionary forces acting on the recent E. burchellii expansions, we examined site‐specific positive selection in the six species‐specific 9‐exon OR expansions in E. burchellii which contained more than four genes, following Engsontia et al. (2015). Site model tests revealed that one of these expansions showed statistically significant signs of positive selection at some sites under the most robust model test (M1a vs. M2a; LRT, p =.009), and one additional expansion was significant under the more powerful tests (M7 vs. M8 and M8a vs. M8; LRT, p =.007 and 0.02, respectively) (Table S3).

Most E. burchellii 9‐exon ORs belong to a single clade (9e‐α ORs sensu McKenzie et al., 2016), and genes in this clade comprise over half of all intact ORs in the E. burchellii genome (Figure 3a). RNA‐seq of pooled antennae from workers and males, pooled bodies (without antennae) from workers and a body (without antennae) from a single male showed that all but one of the 9e‐α ORs are expressed specifically in worker antennae (Figure 3b; Table S4). Most of the remaining ORs showed a mixture of weak male‐ and worker‐biased expression, although 21 additional ORs (mostly in the E and U subfamilies) were expressed specifically in worker antennae, and five ORs were expressed specifically in male antennae. Eleven ORs were expressed in the bodies of workers and/or the male (four and two in workers and the male, respectively, and five in both). Of the ORs classified as expressed in bodies (at a cutoff of 2.0 FPKM), the median expression within those bodies was 10.0 FPKM, while the median expression in the antennae of all ORs classified as expressed in the antennae was 11.6 FPKM (Mann–Whitney U test, p =.26). These findings indicate genuine (though not necessarily functional) expression of a small number of ORs in nonantennal tissues.

FIGURE 3

Tissue‐ and sex‐specific expression of Eciton burchellii odorant receptors. (a) Phylogeny of E. burchellii ORs, with branches belonging to the five largest ant OR subfamilies coloured by subfamily. 9e‐α ORs are designated by a grey box. (b) Heat maps of gene expression in antennae and antennae‐less bodies of E. burchellii workers and males. The first four columns represent absolute gene expression in each tissue, while the last represents fold enrichment in male vs. worker antennae. Expression data should be interpreted cautiously given the lack of biological replicates

To test whether the expansion of 9‐exon ORs was mirrored by an expansion of T6 glomeruli as seen in O. biroi (McKenzie et al., 2016), we next examined the olfactory neuroanatomy of E. burchellii workers and males. Three‐dimensional reconstruction of antennal lobe glomeruli from confocal micrograph stacks of E. burchellii worker and male brains confirmed that E. burchellii workers do indeed have an expanded T6 cluster of glomeruli compared with nondoryline ants with published antennal lobe reconstructions (Camponotus floridanus, C. japonicus and Atta vollenweideri [Kelber et al., 2010; Nakanishi et al., 2010; Zube et al., 2008]) (Figure 4; Table S5). Consistent with the lack of expression of 9‐exon ORs in male antennae, males lacked an expansion of T6 glomeruli (Figure 4; Table S5). E. burchellii also exhibits a corresponding reduction in non‐T6 glomeruli, possessing fewer non‐T6 glomeruli than any other ant yet examined (Table S5).

FIGURE 4

Reconstructed antennal lobes of an Eciton burchellii worker and male, with glomeruli coloured according to input tract. Clusters T1, T2 and T4 could not be reliably distinguished and are all coloured gold. T3 is coloured purple, T5 green, T6 blue and T7 yellow. The antennal nerve is shown in brown. P = posterior (n‐anterior), D = dorsal (n‐ventral), M = medial, L = lateral

Apart from the general reduction in glomerulus number, we observed many exceptionally large glomeruli in clusters T1–4 in workers (Figure 5), which resembled pheromone‐responsive “macroglomeruli” observed in other insect species (Boeckh & Selsam, 1984; Kuebler et al., 2010; Matsumoto & Hildebrand, 1981). Two E. burchellii glomeruli passed the threshold proposed by Kelber et al. (2009) for defining macroglomeruli, and Gaussian mixture modelling showed that 16 glomeruli were significantly larger than all other glomeruli (Figure 5). None of the larger glomeruli belonged to the T6 cluster of glomeruli. We also observed many relatively large glomeruli in E. burchellii males (Figure S4). However, none qualified as macroglomeruli according to the definition of Kelber et al. (2009). Although mclust also found two distinct populations of glomeruli in E. burchellii males (Figure S4), the distribution of glomerular radii in the larger cluster deviated significantly from a Gaussian distribution (Shapiro–Wilk test, p =.004), suggesting that the models tested by mclust did not fit our male data well.

FIGURE 5

Optical slices from a confocal micrograph stack of the Eciton burchellii worker antennal lobe (top and lower left) and a plot of the size distribution of worker glomeruli (lower right). The two largest glomeruli, which pass Kelber et al.’s (2009) threshold for “macroglomeruli,” are indicated by green and magenta arrows. Bars in the size distribution plot are coloured according to their corresponding mclust size distribution cluster classification

Scanning electron microscopy of E. burchellii worker antennae (Figure 6) revealed that basiconic sensilla were present at high density along the ventral surface of the six terminal antennal segments (segments 7–12; directionality is relative to body axis with antennae held in the resting position). This pattern is reminiscent of that observed in O. biroi (McKenzie et al., 2016). A small number of additional basiconic sensilla were detected on the ventral surface of antennal segment 6, while segments 1–5 had no basiconic sensilla. As in other ants, males lacked basiconic sensilla entirely (Figure S5).

FIGURE 6

Scanning electron micrographs of Eciton burchellii worker antennae showing basiconic sensilla restricted to the ventral surface (as defined relative to the body axis while antennae are held naturally) of the seven terminal antennal segments. (a) Whole antennal flagellum with segments numbered (1–12). (b) Magnified terminal antennal segment. A representative basiconic sensillum is circled in both images

DISCUSSION

The New World army ant Eciton burchellii is an extensively studied, prevalent and ecologically important arthropod predator in Neotropical ecosystems. By providing genomic resources for this species, our study will facilitate future investigations into the molecular basis of social predation. We found that E. burchellii possesses a relatively compact genome, with a general reduction in gene complement compared to other ant species. Examination of evolutionary dynamics in functional subgroups of genes did not show accelerated loss of any specific functional subgroup, suggesting relatively even gene loss across the E. burchellii genome. Analysis of repeat content across ants shows that it is tightly correlated with genome size, indicating that repeat content expansion and contraction probably plays a major role in shaping genome size in ants—as previously observed for many different types of organisms (Lynch, 2007). A paucity of interspersed repeats in particular is probably largely responsible for the compact genome size of E. burchellii.

Corresponding to the general reduction in gene complement, we found that E. burchellii possesses one of the smallest intact odorant receptor gene repertoires among the ants we examined in this study. However, the 9‐exon subfamily of ORs is expanded in E. burchellii, and in fact this species possesses the second largest 9‐exon OR repertoire of any ant with high‐quality OR annotations examined in this study. We also found that E. burchellii has an unusually large number of T6 glomeruli in the antennal lobe, and an unusually dense region of basiconic sensilla along the ventral surface of the antenna. Almost all 9‐exon ORs are expressed specifically in worker antennae in E. burchellii. While our lack of replicates prevents us from drawing statistical conclusions about any given OR, this pattern is consistent with observations from Ooceraea biroi, Camponotus floridanus and Harpegnathos saltator (McKenzie et al., 2016; Zhou et al., 2012). All of these data are consistent with previous work which indicated that 9‐exon ORs are expressed in OSNs that are housed in basiconic sensilla and innervate T6 glomeruli (McKenzie et al., 2016). Intriguingly, the clonal raider ant O. biroi, which together with E. burchellii belongs to the ant subfamily Dorylinae, also possesses an expanded 9‐exon OR repertoire and T6 glomerulus cluster, as well as an unusually dense region of basiconic sensilla on the ventral surface of the antennal club (McKenzie et al., 2016). We also observed signatures of positive selection in two of six E. burchellii‐specific 9‐exon OR expansions. This mirrors findings of Engsontia et al. (2015) in O. biroi, where four of six species‐specific 9‐exon OR expansions also showed signatures of positive selection. Engsontia et al. (2015) found only one other positively selected species‐specific clade, a 9‐exon OR expansion in C. floridanus, suggesting that doryline 9‐exon ORs are under unusually high selective pressure.

Both 9‐exon ORs and basiconic sensilla have been shown to respond to CHCs, an important class of ant pheromones used in nest‐mate recognition, fertility signalling and division of labour (Holman et al., 2010; Liebig, 2010; Ozaki et al., 2005; Pask et al., 2017; Sharma et al., 2015; Sturgis & Gordon, 2012; van Zweden & D’Ettorre, 2010). Hydrocarbons have also been shown to function in species recognition and prey selection in ants (Kidokoro‐Kobayashi et al., 2012; Lang & Menzel, 2011; Manubay & Powell, 2020), and basiconic sensilla are specifically implicated in species recognition in the ant Formica yessensis and prey discrimination in the wasp Liris niger (Anton & Gnatzy, 1998; Kidokoro‐Kobayashi et al., 2012). Because both O. biroi and E. burchellii are believed to be specialist ant predators (Franks, 1982; Hoenle et al., 2019; Rettenmeyer, 1963; Tsuji & Yamauchi, 1995), the expansion of the hydrocarbon detection system in O. biroi and E. burchellii may represent adaptations for prey discrimination, particularly for discriminating between prey and nonprey ant species. Unfortunately, O. biroi and E. burchellii are the only doryline ants with sequenced genomes. Additionally, these species share further life history traits not seen in the other ant species in our study, namely nomadism and cyclical reproduction (Oxley et al., 2014; Rettenmeyer, 1963; Schneirla, 1944; Tsuji & Yamauchi, 1995). Future comparative genomics studies with increased phylogenetic and ecological resolution will therefore be necessary to conclusively determine the evolutionary significance, if any, of 9‐exon OR gene expansions.

Recently, two studies have been published examining the genomic consequences of ecological specialization in socially parasitic inquiline (Schrader et al., 2021) and kidnapping (Jongepier et al., 2021) ant species. While these were published too recently to be included in our in‐depth comparative analyses, it is worth noting some striking similarities and differences in these ant groups relative to E. burchellii. Specifically, both studies found large contractions in the OR family in social parasites compared with their free‐living relatives (Schrader et al., 2021, Jongepier et al., 2021). This suggests that olfactory system simplification is a common theme of ecological specialization in ants. However, both studies also showed the strongest contractions in the 9‐exon ORs, while we find a major expansion in this gene subfamily in E. burchellii. This suggests that there may be fundamental differences in the olfactory processes related to finding host colonies to parasitize vs. distinguishing between prey species during raids. Of note, most inquilines and kidnapping ants parasitize closely related ant species (Schrader et al., 2021, Jongepier et al., 2021), and thus they may be ancestrally able to detect volatile host pheromones.

Our neuroanatomical investigation of the E. burchellii AL also revealed many exceptionally large glomeruli, all located outside of the T6 cluster and thus probably not involved in hydrocarbon perception. Enlarged glomeruli, termed macroglomeruli, are often involved in the detection of ecologically particularly important odours, such as sex pheromones in male moths, cockroaches and bees (Boeckh & Selsam, 1984; Matsumoto & Hildebrand, 1981; Sandoz, 2006), and food odours in specialist flies (Ibba et al., 2010). Notably, in the leaf‐cutting ant Atta vollenweideri, large workers (majors) possess a single macroglomerulus that responds to trail pheromones (Kuebler et al., 2010). Workers in C. floridanus also possess a single enlarged glomerulus, although its function is unknown (Zube et al., 2008). Most insects with macroglomeruli possess one or two such structures. However, our analysis of glomerulus size distribution suggested that 16 glomeruli in E. burchellii workers belong to a distinctly larger size class. As they are located outside of the T6 cluster, we hypothesize that these glomeruli are involved in volatile pheromone detection, such as alarm, trail and recruitment pheromones. Functional studies will be necessary to confirm pheromone responses in these glomeruli. While our study was restricted to males and medium‐sized workers, future studies on neuroanatomical differences between worker subcastes may also prove informative, given the above‐mentioned major‐restricted macroglomeruli observed in Atta. Intriguingly, previous work on AL size in Eciton has shown a reduction in AL size in Eciton majors relative to the brain as a whole (O’Donnell et al., 2018). Future work should also examine differences in neuroanatomy and genomics of different E. burchellii populations, as recent work has found striking genetic differentiation between E. burchellii subspecies (Winston et al., 2017).

Overall, our findings suggest that doryline ants possess a suite of chemosensory system traits that have so far not been observed in other ants. These include an expanded repertoire of putative hydrocarbon sensors, and a dense clustering of hydrocarbon‐sensitive sensilla along the ventral surface of the antennae. These traits may represent adaptations for specialized predation and myrmecophagy. Moving forward, it will be important to investigate the chemosensory systems of other army ant‐like species with separate evolutionary origins, such as the ponerine genera Simopelta and Leptogenys, the amblyoponine genus Onychomyrmex, and the army ant‐like members of the myrmicine genus Carebara (previously Pheidologeton). The E. burchellii genome will also provide an invaluable resource for future studies of the genomics of social predation and army ant biology in general.

Data Accessibility

Sequencing data and genomic assemblies and annotations have been deposited in the NCBI SRA and genome databases under the BioProject accession no. PRJNA632625. All other data are available as supplemental material accompanying this paper.

CONFLICT OF INTEREST

N/A.

AUTHOR CONTRIBUTION

N/A.

Fig S1‐5

Click here for additional data file.

Table S1

Click here for additional data file.

Table S2

Click here for additional data file.

Table S3

Click here for additional data file.

Table S4

Click here for additional data file.

Table S5

Click here for additional data file.