Causes and consequences of a complex recombinational landscape in the ant Cardiocondyla obscurior

Initial linkage mapping and genome assembly refinement and annotation

To improve the current genome assembly (Cobs2.1) to chromosome resolution, we produced an initial ordering of 57,557 markers (out of 342,508 raw biallelic SNPs) into 29 LGs (Supplemental Fig. S1) using whole-genome data from 100 F₂ offspring of a single F₁ mother obtained from a cross between a male and a queen belonging to two Brazilian populations of C. obscurior. These populations were specifically chosen to maximize the number of heterozygous markers and minimize runs of homozygosity in the experimental cross. Although some small LGs may be physically connected to each other or to larger LGs, they were not combined owing to large genetic distances (∼50 cM), indicating that they segregate independently from other LGs.

The new genome (Cobs3.1) spans 193.05 Mb and is an improvement over the previous version (Cobs2.1) with ∼96% (184.4 Mb) of the physical assembly placed into pseudomolecules and a scaffold N50 of 7.38 Mb compared with 6.29 Mb of Cobs2.1 (Supplemental Table S1). Using the Liftoff tool (v.1.6.3) (Shumate and Salzberg 2021), we converted the coordinates of 20,944 genes (36,599 transcripts) annotated in Cobs2.1 to Cobs3.1.

Transposable elements (TEs) exhibit a bimodal distribution in the genome of C. obscurior with regions of high TE content (“TE islands”) in a genomic background of low TE content (low-density regions [LDRs]) (Supplemental Fig. S2A), as described previously (Schrader et al. 2014; Errbii et al. 2021). A total of 30 TE islands are located at the center (e.g., LGs 1, 2, 3, and 5) and the ends (e.g., LGs 4, 8, and 9) of LGs. We also observed secondary TE accumulations on several LGs (e.g., 1, 5, and 10).

Tandem repeat (TR) sequences are commonly associated with centromeres across various species, including ants (Melters et al. 2013; Huang et al. 2016). Upon annotating TRs in C. obscurior, we observed an overlap of these sequences with TE islands (Supplemental Fig. S2B). Candidate centromeric repeats from five different ant species (Melters et al. 2013; Huang et al. 2016) were not enriched in TE islands, suggesting a distinct repertoire of repeats in centromeres of C. obscurior, consistent with prior findings of rapid divergence among centromeric TRs in ants (Huang et al. 2016). Together, these results suggests that TE islands may serve as potential centromeres for the chromosomes of C. obscurior.

The recombination landscape of C. obscurior

To estimate recombination rates, we used the updated assembly (Cobs3.1) for a second round of linkage mapping with 57,309 high-quality markers (out of 342,616 raw biallelic SNPs) covering 84% (162.9 Mb) of the physical assembly and revealing 1061 effectively recombining positions. It is important to note that the number of recombining positions is independent of marker density and is instead determined by the occurrences of crossover events in the experiment. The final map (1315.4 cM in 29 LGs) (Supplemental Fig. S3) showed a positive correlation between genetic and physical positions of markers, indicating a high-quality map (Supplemental Fig. S4, correlation coefficients are given in the figure). Based on the proportion of the physical assembly covered by markers, the average recombination rate in C. obscurior is 8.1 cM/Mb. Although higher than observed in other nonsocial insects and animal species (Supplemental Table S2), this estimate falls within the range of other ants (Sirviö et al. 2006, 2011b) and is lower than those of the honey bee (Beye et al. 2006; Shi et al. 2013; Liu et al. 2015; Rueppell et al. 2016).

Recombination rate varied considerably among the 29 LGs (Table 1) with LGs <5 Mb having a higher recombination rate (median = 14.2 cM/Mb, 95% CI [8.52, 18.06]) compared with LGs >10 Mb (median = 7.5 cM/Mb, 95% CI [3.65, 10.29]) and intermediate LGs (≥5 Mb and ≤10 Mb; median = 7.5 cM/Mb, 95% CI [4.60, 9.55]). Physical and genetic distances correlated positively (Spearman's ρ = 0.57, P = 0.0014) (Supplemental Fig. S5A), whereas average recombination rate per LG showed a significant negative correlation (Spearman's ρ = −0.39, P = 0.039) with physical lengths of LGs (Supplemental Fig. S5B).

Eight LGs showed very little or no recombination because of a reduced number of crossovers (LGs 4, 6, 12, 17, 18, 26, 28, 29) (Table 1). Although the number of markers varied between LGs (Supplemental Fig. S6), there was no correlation between average recombination rates and marker count (Spearman's ρ = −0.087, P = 0.65) per LG (Supplemental Fig. S5C). In the following, we investigated potential factors contributing to the variation in the number of crossovers and, consequently, recombination rate in C. obscurior.

Recombination rate distribution and genome architecture in C. obscurior

To investigate the genomic determinants of recombination in C. obscurior, we estimated the local recombination rate in 250 kb nonoverlapping windows and analyzed these with respect to various sequence parameters (Fig. 1A; Supplemental Figs. S7, S8). Because crossovers can only be detected between two markers, recombination rate estimates were restricted to regions between the first and last markers on each chromosome, excluding the chromosome ends.

Figure 1.

Recombination rate variation and genome architecture in C. obscurior. (A) Circos plot showing the following inward: (SV) the location of identified structural variants (dark green), (RR) recombination rates across the genome (cM/Mb) with introgressions depicted with purple bars, (TE) TE content (“LDRs”, blue; “TE islands”, orange) with the location of TE islands shown as black bars, (Ex) exon content, and (π) the log₁₀-transformed nucleotide diversity in two natural populations for C. obscurior from Itabuna (gold) and Una (maroon). C. obscurior queen with brood (photo: L. Schrader). (B) Recombination rates between TE-poor regions (“LDRs” in blue) and TE-rich regions (“TE islands” in orange). (C) Recombination rates in introgressed regions (“LDRs introgr.” In purple) compared with the remainder of LDRs (“LDRs other” in blue). (D) Levels of nucleotide diversity (π) in nonrecombining (“LDRs nonrec.”) compared with recombining regions LDRs (“LDRs rec.”) in two populations of C. obscurior from Itabuna and Una. (****) P < 0.0001, (**) P < 0.01.

Open in new tabDownload PowerPointLink to figure

The recombination rate in C. obscurior ranged from 0 to 82.2 cM/Mb, with ∼38% (∼70 Mb) of the analyzed sequence showing minimal to no recombination (e.g., on LGs 6, 12, 17, and 18) (Fig. 1A; Supplemental Fig. S7A). Note that the recombination rate estimates were not affected by marker density. Although a correlation was observed between marker density and recombination rate in 250 kb windows (Spearman's ρ = 0.075, P = 0.049) (Supplemental Fig. S9A), this correlation disappeared when considering only recombining windows (Spearman's ρ = 0.045, P = 0.37) (Supplemental Fig. S9B).

Markers were unevenly distributed across the LGs, with 23 out of the 29 LGs having median intermarker distances below 0.64 kb (range: 0.12 to 256.97 kb) (Supplemental Table S3). Three LGs 6, 12, and 18 showed minimal recombination and displayed particularly large intermarker distances (medians = 210.27, 37.84, and 256.97 kb, respectively), 90% of which were <750 kb (Supplemental Table S3). Because of chiasma interference (Sturtevant 1915; Muller 1916), in which the occurrence of one crossover disrupts the formation of another nearby, these intervals should still be sufficient to capture crossovers. Unless exceptionally rare double crossovers occur within these short intervals, we expect an average of one crossover every 2.7 Mb (Liu et al. 2015) to 7.14 Mb (Sirviö et al. 2011b), according to studies in the honey bee and harvester ants, respectively.

Recombination rates are significantly reduced in TE islands (median = 0 cM/Mb), which are gene-poor, compared with LDRs (median = 3.11 cM/Mb) in which gene density is higher (Wilcoxon test: W = 43644, P = 5.25 × 10⁻¹³) (Fig. 1B). Note that the decrease in recombination rate within TE islands is not because of a lack of markers (Supplemental Fig. S6; Supplemental Table S4). LDRs on their own show no significant correlation between recombination rate and TE content (Spearman's ρ = 0.016, P = 0.74) or exon content (Spearman's ρ = −4 × 10⁻³, P = 0.93) (Supplemental Fig. S8).

Beside TE islands, several other regions in LDRs also showed little or no recombination (Fig. 1A). Some of these nonrecombining regions are characterized by high sequence divergence between grandparental haplotypes (gp-d_XY) (Supplemental Fig. S7B). In these regions, one of the two grandparental haplotypes showed high sequence similarity to Asian and European populations of C. obscurior (Supplemental Fig. S7C), indicating introgression from the Old-World lineage (Errbii et al. 2021). Based on a cutoff of threefold mean gp-d_XY (∼0.001), these introgressions occupy a large fraction (>40%) of LGs 9, 10, 17, and 21 and smaller regions on LGs 1 and 2. After accounting for differences in the number of introgressed (53) and nonintrogressed (330) windows, we found that regions with high gp-d_XY recombine significantly less (median = 1.8 cM/Mb) than the rest of the genome (median = 8.52 cM/Mb) (Wilcoxon test: W = 1043.024, P = 0.046) (Fig. 1C). Several inversions overlap with these introgressions (on LGs 10, 17, and 21) (Supplemental Table S5). Together, this indicates that strong sequence divergence among homologous chromosomes impairs recombination in C. obscurior, at least in part caused by polymorphic inversions.

Other nonrecombining regions on LGs 1, 4, 6, 8, and 12 showed no signs of introgression and contained no large-scale inversions that could explain local suppression of recombination (Fig. 1A; Supplemental Fig. S7A,B; Supplemental Table S5). Considering that meiotic recombination disrupts the linkage between segregating loci, we thus hypothesized that the absence of recombination in these regions could result from epistatic selection acting to reduce recombination in regions harboring coadapted loci (Turner 1967; Charlesworth 2016). To explore this possibility, we screened these regions for signs of reduced genetic variation and Tajima's D, measures commonly used in population genetics to infer signatures of selection (Nielsen 2005). The nonrecombining regions on LGs 1, 4, 6, 8, and 12 showed significantly reduced levels of nucleotide diversity in natural populations (Itabuna_median = 2.18 × 10⁻⁵; Una_median = 4.04 × 10⁻⁵) compared with the remainder of the LDRs (excluding introgressions; Itabuna_median = 1.79 × 10⁻⁴; Una_median = 2.03 × 10⁻⁴; Wilcoxon test; Itabuna: W = 5559.5, P < 2.2 × 10⁻¹⁶; Una: W = 6289.5, P < 2.2 × 10⁻¹⁶) (Fig. 1D). Similarly, Tajima's D was also significantly reduced in these nonrecombining regions (Itabuna_median = −1.86; Una_median = −2.1) compared with the remainder of the LDRs (excluding introgressions; Itabuna_median = −0.99; Una_median = −1.37; Wilcoxon test; Itabuna: W = 9547, P = 1.32 × 10⁻¹⁵; Una: W = 9860, P = 1.07 × 10⁻¹⁴) (Supplemental Figs. S7D, S10).

We further applied a generalized linear model with recombination as the response variable and found that TE content, sequence divergence, and epistasis were all significant explanatory variables (Supplemental Table S6).

Finally, there was no correlation between recombination rate and GC content (overall: Spearman's ρ = 0.044, P = 0.35; intergenic: Spearman's ρ = 0.025, P = 0.6) (Supplemental Figs. S7E, S8) in C. obscurior.

Meiotic recombination contributes to genetic variation in natural populations and promotes intraspecific divergence

In support of a substantial effect of recombination on patterns of genetic variation, we found a positive correlation of recombination rates with measures of nucleotide diversity (Itabuna: Spearman's ρ = 0.28, P = 2.1 × 10⁻⁹; Una: Spearman's ρ = 0.25, P = 9.9 × 10⁻⁸) and Tajima's D (Itabuna: Spearman's ρ = 0.21, P = 7.1 × 10⁻⁶; Una: Spearman's ρ = 0.3, P = 5.6 × 10⁻¹¹) calculated in 250 kb windows in LDRs (Fig. 2).

Figure 2.

Relationship between recombination rates and measures of diversity in populations of C. obscurior. (Left) Correlation between recombination rates and nucleotide diversity and Tajima's D. (Middle) Correlation between recombination rates and nucleotide diversity at synonymous (π_S) and nonsynonymous (π_NS) sites. (Right) Correlation between recombination rates and levels of genetic differentiation (F_ST) and absolute divergence (d_XY) between lineages of C. obscurior. The x-axis represents binned recombination rates, with the upper limit of each bin indicated and “TEI” referring to TE islands (orange boxplot). (ρ) Spearman's rank correlation coefficient, with significance shown by asterisks: (****) P < 0.0001, (***) P < 0.001, (**) P < 0.01. All correlations were performed using the raw data. Lines and shaded areas are linear regressions and 95% confidence intervals.

Open in new tabDownload PowerPointLink to figure

Recombination rates were also positively correlated with levels of diversity at both nonsynonymous (π_NS; Itabuna: Spearman's ρ = 0.21, P = 1 × 10⁻⁵; Una: Spearman's ρ = 0.22, P = 3.7 × 10⁻⁶) and synonymous sites (π_S; Itabuna: Spearman's ρ = 0.26, P = 3.9 × 10⁻⁸; Una: Spearman's ρ = 0.22, P = 2.3 × 10⁻⁶) (Fig. 2), but not with their ratio (π_NS/π_S; Itabuna: Spearman's ρ = 0.0062, P = 0.91; Una: Spearman's ρ = 0.0018, P = 0.97) (Supplemental Fig. S11), indicating that meiotic recombination does not correlate with positive or purifying selection.

Consistent with models of linked selection in which diversity is reduced in regions of low recombination (Cutter and Payseur 2013; Burri 2017), there was a negative correlation between F_ST and recombination in two different population pairs (F_ST [Itabuna-Tenerife]: Spearman's ρ = −0.18, P = 1.7 × 10⁻⁴; F_ST [Una-Tenerife]: Spearman's ρ = −0.12, P = 9.9 × 10⁻³) (Fig. 2). However, there was also a significant positive association between recombination and d_XY in two additional population pairs (d_XY [Itabuna-Leiden]: Spearman's ρ = 0.41, P < 2.2 × 10⁻¹⁶; d_XY [Una-Taiwan]: Spearman's ρ = 0.37, P = 8.6 × 10⁻¹⁶) (Fig. 2). Together, this indicates that recombination contributes to within-population polymorphism and to the divergence of C. obscurior lineages. Future studies should explore whether the same patterns persist when using recombination rates estimated for a population of the Old World lineage.

Recombination is associated with small and large genome rearrangements

Ectopic recombination (i.e., crossovers between nonallelic homologous loci), particularly in repeat-rich regions, has been shown to produce structural variants (SVs) (Liu et al. 2012; Savocco and Piazza 2021). To investigate to what extent the high recombination rates in C. obscurior are correlated with SVs, we used whole-genome sequencing data generated previously for 16 individuals from five populations (Errbii et al. 2021).

Analysis of 481 SVs (129 deletions, 54 duplications, and 298 inversions) showed a positive correlation between recombination rates and SV frequencies and their lengths in LDRs (SVs length: Spearman's ρ = 0.14; P = 0.0025; SVs frequency: Spearman's ρ = 0.15, P = 8.1 × 10⁻⁴) (Fig. 3A,B). When analyzed separately, none of the individual SV types correlate significantly with recombination (Supplemental Fig. S12). Additionally, genetic diversity at 95,841 short insertions and deletions (indels) showed a strong and significant positive association with recombination (Spearman's ρ = 0.54, P < 2.2 × 10⁻¹⁶) (Fig. 3C). A similar pattern was observed when genic and intergenic indels were analyzed separately (Supplemental Fig. S13).

Figure 3.

Association between recombination rates and genome dynamics in C. obscurior. Correlations between recombination rates and length (A) and frequency (B) of SVs and indel polymorphism (C) in 16 individual workers. The x-axis represents binned recombination rates, with the upper limit of each bin indicated. (ρ) Spearman's rank correlation coefficient is indicated for all structural variants (SVs) and separately for deletions (DEL), duplications (DUP), and inversions (INV), with significance shown by asterisks: (****) P < 0.0001, (***) P < 0.001, (**) P < 0.01. All correlations were performed using the raw data. Red lines and shaded areas are linear regressions and 95% confidence intervals.

Open in new tabDownload PowerPointLink to figure

No association between recombination and evolutionary rates

Meiotic recombination is expected to enhance the efficacy of selection by uncoupling linkage between targets of selection and mitigating the impact of Hill–Robertson interference (HRI) (Hill and Robertson 1966; Betancourt et al. 2009; Campos et al. 2014). To examine this relationship, we used omega (ω) as a proxy, that is, the ratio of nonsynonymous (d_N) to synonymous (d_S) substitutions (Webster and Hurst 2012). In regions of low recombination, increased genetic interference should compromise both positive selection for advantageous variants and purifying selection against deleterious ones (Betancourt et al. 2009). Given that the majority of nonsynonymous mutations are deleterious, a negative correlation between recombination rate and ω is expected (Haddrill et al. 2007; Betancourt et al. 2009; Bolívar et al. 2016). However, we found no significant correlation between recombination and d_N (Spearman's ρ = 0.008, P = 0.61), d_S (Spearman's ρ = −0.015, P = 0.35), or ω (Spearman's ρ = 0.019, P = 0.24) for 3916 single-copy ortholog genes in LDRs (Supplemental Fig. S14).

No association between recombination and gene expression bias

In honey bees, genes with caste-biased expression in brains were enriched in regions of high recombination, suggesting that caste-specific transcriptional plasticity is facilitated by the effects of recombination (Kent et al. 2012; Liu et al. 2015; Wallberg et al. 2015; Jones et al. 2019). In C. obscurior genes showing caste-biased brain expression are as yet unknown. However, plastic gene expression has been studied in this species in the context of caste differentiation (Schrader et al. 2015, 2017). To assess a potential association between recombination and developmental gene expression plasticity, we reanalyzed this RNA-seq data set of third instar larvae of different castes and sexes (Supplemental Fig. S15). Based on 6807 expressed genes, we did not find a significant association between recombination and gene expression plasticity over multiple contrasts (Spearman's ρ = 0.014, P = 0.26) (Supplemental Fig. S16).

Functional analysis of genes in regions of high recombination

To explore whether regions of high recombination rates in C. obscurior show signatures of functional specialization, we conducted Gene Ontology (GO) enrichment analyses on genes contained in the 10% most frequently recombining regions in the genome. These regions showing recombination rates >22.3 cM/Mb and containing 1541 genes were significantly enriched for fatty acid elongase and reductase genes (Table 2; Supplemental Table S7). Although this implies a potential association between recombination and cuticular hydrocarbon (CHC) biosynthesis, it is important to acknowledge the limitations of these findings, particularly the failure to determine the specific genes directly affected by recombination.

Linkage mapping and recombination rate for C. obscurior

The high recombination rates reported for the social Hymenoptera are not easily explained and have not been thoroughly investigated in ants. To start filling this gap, we produced a high-resolution genetic map. At least 26 of the 29 LGs represent full chromosomes, whereas the remaining small LGs are most likely unlinked fragments of larger chromosomes or extra ones in the mapping population. The updated chromosome-level genome of C. obscurior (Cobs3.1) confirms the architecture of the previous assembly-based versions and that it is the smallest ant genome described so far (Schrader et al. 2014; Errbii et al. 2021).

Average recombination rate in C. obscurior is 8.1 cM/Mb, which is within rates reported for two ants (Sirviö et al. 2006, 2011b) and is lower compared to several honey bee species (Hunt and Page 1995; Beye et al. 2006; Liu et al. 2015; Rueppell et al. 2016). This value is likely an underestimate because of two reasons. First, we excluded unlikely double-crossover events as false positives. Although such double crossovers are highly improbable owing to chiasma interference (Sturtevant 1915; Muller 1916), they may occur at low frequencies. Second, recombination in our experiment was suppressed in a large fraction of the genome (∼38%), partly owing to introgressed haplotypes in the grandparental cross. These regions are likely to recombine freely in the absence of introgression, which should increase average recombination rates quite significantly.

Two factors can affect the resolution of the map. First is the number of markers and their distribution across the genome. Despite the clustering of nonrecombining markers, which reduced the effective number of markers available for recombination rate estimation (from 57,309 markers to 1061 clusters), the map maintained its density. The average intermarker distance in this study stands at 1.24 cM and is relatively lower than the other two available ant maps (A. echinatior: 145 markers or an intermarker distance of 14.32 cM [Sirviö et al. 2006] and P. rugosus: 215 markers or an intermarker distance of 16.55 cM [Sirviö et al. 2011b]). A second factor that affects map resolution is the size of the mapping population. Here, we used data from 100 F₂ individuals from a single queen, comparable to previous maps in ants (Sirviö et al. 2006, 2011b) and bees (Liu et al. 2015, 2017; Rueppell et al. 2016) and generally considered sufficient for recombination rate estimation and the construction of a well-resolved linkage map.

Recombination rates varied among the 29 LGs, showing a negative correlation with chromosome size, a feature shared among many species (Hunt and Page 1995; Jensen-Seaman et al. 2004; Pessia et al. 2012; Haenel et al. 2018). Several LGs in our map were ∼50 cM long, aligning with the expectation under the one obligate crossover rule. This raises the question whether the high recombination rates observed in eusocial species can be attributed to an increase in the number of chromosomes with obligate crossovers. However, chromosome numbers of eusocial and solitary species of Hymenoptera do not differ significantly (Ross et al. 2015). Furthermore, these high rates can also not be linked to genome size reductions in social Hymenoptera, as smaller genomes do not generally exhibit higher recombination rates in animals (Stapley et al. 2017). Thus, neither chromosome number nor genome size is sufficient to explain the high recombination rates of social Hymenoptera.

The factors affecting recombination rates in C. obscurior

Recombination rates, calculated in 250 kb windows, varied considerably (0 to 82.2 cM/Mb), resulting in a complex recombinational landscape. Such regions of extremely high recombination rates have also been observed in the honey bee (Beye et al. 2006; Liu et al. 2015). Consistent with previous observations indicating higher recombination rates in gene-dense regions and open chromatin (Tiley and Burleigh 2015; Marand et al. 2017), gene-rich LDRs displayed higher recombination rates compared with gene-poor TE islands (Schrader et al. 2014; Errbii et al. 2021). This supports the interpretation that TE islands are (peri-)centromeric regions, where recombination is usually suppressed (Beadle 1932; Lambie and Roeder 1986; Talbert and Henikoff 2010). TE insertions can accumulate in regions with reduced recombination because ectopic recombination as a mechanism to remove TEs is not available (Dolgin and Charlesworth 2008; Kent et al. 2017). In accordance, centromere structure and location in heterochromatin can locally suppress recombination (Talbert and Henikoff 2010), allowing accumulations of TEs in these regions to evolve secondarily. Future cytological examinations are necessary to confirm that these regions are indeed centromeres.

Sequence divergence resulting from introgression between distant populations of C. obscurior (Errbii et al. 2021) also inhibit recombination, in line with studies in other taxa (Liharska et al. 1996; Opperman et al. 2004; Li et al. 2006; Delame et al. 2019). Particularly, inversions that occur in some of these introgressions can significantly impact recombination (Noor et al. 2001; Rieseberg 2001; Stevison et al. 2011) by disrupting the pairing of homologous chromosomes and synapsis formation during meiosis (Chakraborty et al. 2019).

Recombination was also suppressed in regions that showed low TE content and no signs of introgression, as well as significantly reduced levels of genetic diversity in two populations. Such patterns can emerge under epistatic selection (Otto and Lenormand 2002; Hadany and Comeron 2008; Takahasi 2009; Fuller et al. 2020; Navarro-Domínguez et al. 2022), which acts to reinforce linkage between coadapted and/or essential loci through recruitment of recombination modifiers (Otto and Feldman 1997). These modifiers become associated with high-fitness allele combinations favored by selection, rendering any alternative allelic combination deleterious and making recombinant offspring less viable (Altenberg and Feldman 1987; Barton 1995; Singhal et al. 2019). Future studies examining the functional role of these regions are necessary to resolve whether they indeed evolved because of epistatic selection of coevolved loci.

Associations between GC content and recombination were described in the honey bee (Wallberg et al. 2015; Jones et al. 2019) and other organisms (Meunier and Duret 2004; Pessia et al. 2012), most likely as a result of GC-biased gene conversion (Webster and Hurst 2012). However, our investigation in C. obscurior, consistent with recent findings in butterflies (Torres et al. 2023), revealed no significant correlation between recombination rates and GC content. This mirrors patterns seen in selfing plants like Arabidopsis thaliana, in which it has been suggested to result from increased inbreeding and a low N_e (Marais et al. 2004).

No support for hypotheses regarding the role of recombination in social evolution

It has been suggested that recombination contributes to diversity around genes involved with social life. By increasing the efficiency of natural selection (Felsenstein 1974; Comeron et al. 2008), high recombination rates may promote the evolution of caste-specific genes (Kent et al. 2012; Kent and Zayed 2013). Studies in Drosophila have suggested a link between recombination and the efficacy of selection (Betancourt and Presgraves 2002; Haddrill et al. 2007; Campos et al. 2014), whereas other reports found no compelling evidence for such an association (Webster and Hurst 2012), in line with data from humans (Bullaughey et al. 2008). Consistent with a recent study in honey bees (Jones et al. 2019), in C. obscurior we found no association between recombination rate and rates of d_N/d_S, used as a proxy for selection efficiency. Hence, our finding does not support the hypothesis that high recombination rates promoted the evolution of sociality by reducing HRI effects (Hill and Robertson 1966) and allowing the evolution of genes important for social life.

Alternative models suggested that elevated rates of recombination in social Hymenoptera may increase resistance to pathogens and facilitate behavioral and morphological differentiation by increasing genetic diversity at the colony level (Shykoff and Schmid-Hempel 1991; Hughes and Boomsma 2004; Wilfert et al. 2007; Liu et al. 2015; Wallberg et al. 2015; Jones et al. 2019). We found no enrichment of genes involved in immunity or behavioral plasticity in regions of high recombination. Additionally, there was no correlation between recombination rates and genes showing plastic expression in larvae of C. obscurior, indicating that high recombination rates are also not related to caste differentiation in this species.

GO term enrichment analyses revealed a significant enrichment of elongase and reductase genes in regions of high recombination rates. These genes are primarily known for their roles in regulating fatty acids' elongation and their conversion into fatty alcohols during the biosynthesis of CHCs (Holze et al. 2021), which play an important role in nestmate recognition in ants (Hefetz 2007; Sprenger and Menzel 2020), including C. obscurior (Drakula et al. 2023). Gene family expansions of elongases and reductases have been suggested to contribute to the evolution of CHC complexity (Finck et al. 2016; Hartke et al. 2019) and pheromone communication (Tupec et al. 2019) in ants. To what extent high recombination rates could contribute to the function or evolution of this key component of ant communication has so far not been explored. Although the direct impact of recombination on these genes or their surrounding regulatory sequences remains unclear, it is plausible that frequent recombination at these loci could facilitate the production of unique colony odors by increasing genotypic variation between colonies in genes involved in the CHC biosynthesis pathway (Holze et al. 2021). Future studies should investigate whether this pattern can be found in other ants and whether ectopic recombination might have facilitated the expansion of these gene families.

Recombination promotes genetic diversity of haplodiploid social insects with low N_e

In haplodiploids, N_e is expected to be three-fourths that of diploid species, theoretically reducing genetic variation and constraining adaptation to environmental variation (Wright 1931; Kimura and Crow 1964). In eusocial haplodiploids, N_e and, consequently, genetic variation are reduced further (Metcalf et al. 1975; Pamilo et al. 1978; Pamilo and Crozier 1997; Romiguier et al. 2014), as only a fraction of the individuals of a population reproduces (but for other potential determinants of N_e in Hymenoptera, see Weyna and Romiguier 2021). Studies in taxa beyond Hymenoptera, such as social shrimp (Chak et al. 2022) and termites (Roux et al. 2024), corroborate the association between eusociality and reduced N_e. So, how is genetic variation as the substrate for natural selection and species evolution maintained in eusocial haplodiploids?

Our study reveals a positive correlation between genetic diversity in populations of C. obscurior and recombination rates, consistent with previous findings (Begun and Aquadro 1992; Sirviö et al. 2006; Spencer et al. 2006; Roesti et al. 2013; Jones et al. 2019). Recombination rates also correlated positively with SV frequency and indel polymorphisms, likely owing to ectopic recombination—a significant mechanism driving SV emergence in various organisms, including honey bees (Rueppell et al. 2016), Drosophila (Montgomery et al. 1991; Delprat et al. 2009), and humans (Liu et al. 2012; Robberecht et al. 2013; Startek et al. 2015). Thus, high recombination rates may function as a compensatory mechanism for the reduced genetic diversity observed in haplodiploid social Hymenoptera with low N_e. Although there is some reasoning that increased recombination may evolve in small populations (Otto and Barton 2001), further investigations are necessary to elucidate the relationship between N_e and recombination rate evolution in social and nonsocial haplodiploids, particularly using empirical approaches.

Recombination is associated with genome dynamics and divergence in C. obscurior

The correlation between genetic diversity and recombination in C. obscurior has important implications for its genome dynamics and evolutionary processes. Through hitchhiking and background selection, linked selection in regions with low recombination leads to lower diversity and a skew toward rare alleles (Charlesworth et al. 1993; Comeron et al. 2008; Sella et al. 2009; Smukowski and Noor 2011; Ellegren and Galtier 2016). Conversely, the divergence of populations by genetic drift is expected to be stronger in regions where recombination rates and, thus, genetic diversity are high (Hellmann et al. 2003; Ellegren and Galtier 2016). Accordingly, in C. obscurior, recombination rates are positively correlated with genetic diversity and levels of d_XY. However, there was no correlation with ω, consistent with a role of recombination in driving within-species divergence but not between-species divergence. Moreover, C. obscurior is a highly inbred tramp species and its populations experience genetic bottlenecks, reducing diversity in this species (Heinze 2017). Along with genetic hybridization and an increased TE activity (Errbii et al. 2021), introduced populations of this species may benefit from the effects of recombination on genetic diversity, providing yet another possible resolution to the “genetic paradox of invasive species” (Frankham 2005).

Mapping population

C. obscurior originates from Southeast Asia, and colonies across the globe belong to two (or more) genetically distinct lineages (Errbii et al. 2021). In this study, we performed a cross involving a male and a queen collected from two New World populations separated by ∼50 km (Itabuna and Una, respectively, in Brazil). These samples were collected in 2018 under permit 63371-1 issued by the Brazilian Ministério do Meio Ambiente. In haplodiploid species, haploid males develop from unfertilized eggs that represent the direct product of meiosis in females—the only recombining sex—making them ideal for linkage mapping. However, C. obscurior queens are not very productive (Jaimes-Nino et al. 2022) and cannot produce sufficient males for linkage mapping. To compensate, the F₁ queen was mated with her brother (herein F₁ male) to sample a total of 100 F₂ individuals (12 males and 88 workers) at the pupal stage.

DNA extraction and sequencing

DNA was extracted from the two grandparents, the F₁ parents, and 100 F₂ individuals using a modified CTAB (0.75 M NaCl, 50 mM Tris/HCl at pH 8.0, 10 mM EDTA, 1% hexadecyltrimethylammonium bromide) protocol (Doyle and Doyle 1987). Libraries were prepared using the Illumina DNA prep kit with five to six PCR cycles, and 150 bp paired-end read sequencing was performed on an Illumina NovaSeq platform at the Cologne Center for Genomics. The grandmother and the F₁ queen were sequenced to an average coverage of ∼60×, the males to an average coverage of ∼27×, and the workers to an average coverage of ∼42×.

Read mapping and variant calling

The quality of the raw reads was assessed using FastQC (v.0.11.7) (https://www.bioinformatics.babraham.ac.uk/projects/fastqc/), and Trimmomatic (v.0.38) (Bolger et al. 2014) was used to remove short and low-quality reads as well as adapter sequences (options: ILLUMINACLIP:NexteraPE-PE.fa:2:30:10 SLIDINGWINDOW:4:20 MINLEN:40). The resulting paired reads were mapped to the Cobs2.1 assembly (Errbii et al. 2021) using BWA-MEM (v.0.7.17) (Li and Durbin 2009) with default parameters. The mapping quality was checked using Qualimap (v.2.2.2-dev) (Supplemental Table S8; Okonechnikov et al. 2016), and duplicate reads were marked using Picard's MarkDuplicates (v.2.20.0). Joint-variant calling was performed using GATK's HaplotypeCaller (v.4.1.2.0) (McKenna et al. 2010) with default parameters.

Variant filtering and phasing

After removing indels, the remaining 342,508 biallelic single-nucleotide polymorphisms (SNPs) were screened for high-quality markers using the following criteria: First, markers must be genotyped in both F₁ parents and the grandparents. Second, in the F₁ queen, markers must be heterozygous, whereas in the grandparents, they must be homozygous for different alleles or heterozygous in the grandmother. Because males are haploid, each SNP is expected to be hemizygous in males. Third, consequently, heterozygous markers in at least one male, likely owing to copy number variations, were discarded. Fourth, in BCFtools (v.1.14) (Danecek et al. 2021), genotypes with Mendelian errors were removed. Markers were further filtered using VCFtools (v.0.1.16) (Danecek et al. 2011) to only include those with a minor allele frequency (MAF) > 0.2, genotyped across >80% of samples, a mean coverage between 19× and 50×, and at least five reads per genotype. In total, we identified 57,558 high-quality markers that were heterozygous in the F₁ queen and were used as the basis for phasing (Supplemental Fig. S17).

The mapping population included 88 diploid F₂ workers, each carrying the F₁ father's haplotype and one of the two recombining haplotypes from the F₁ queen (Supplemental Fig. S17). Using a custom script (Supplemental Code), we excluded the F₁ male contribution from each of the F₂ workers, retaining only the maternal contribution.

Scaffold anchoring into chromosomes

We improved the current version of the genome (Cobs2.1) by correcting misassembled (i.e., chimeric) and misoriented scaffolds. We used MSTmap (v.1.0) (Wu et al. 2008), with parameters specified in Supplemental Table S9, to order all 57,558 markers into LGs according to the method of Gadau (2009). We then generated a BED file with the corrected genomic coordinates and orientation of scaffolds (Supplemental Table S10) based on the initial genetic map (Supplemental Fig. S1). Using BEDTools’ getfasta (v.2.30.0) (Quinlan and Hall 2010) and Cobs2.1, we extracted and concatenated genomic sequences belonging to the same LG.

To accurately split chimeric scaffolds, we first checked for the existence of gaps separating the parts belonging to different LGs. When no gap was found, scaffolds were split after the last marker of the first part. For misoriented scaffolds, we inspected the correlation between the physical and genetic positions of markers. When the correlation was positive, the original orientation was kept, and when negative, the reverse complementary sequence was used. Scaffolds lacking markers and those that could not be mapped to any of the LGs were kept unplaced. The final chromosome-level assembly contained 29 LGs and 84 unplaced sequences.

Gene and (tandem) repeat annotation

For gene annotation, Liftoff (v.1.6.3) (Shumate and Salzberg 2021) was used to transfer the Cobs2.1 annotation (Errbii et al. 2021) to Cobs3.1 with the default parameters and the -polish option. For repeat annotation, we used RepeatMasker (v.4.0.7) (https://www.repeatmasker.org) and a previous TE library (Errbii et al. 2021). TE islands were defined as previously described (Errbii et al. 2021).

The tandem repeats finder tool (v.4.09.1) (Benson 1999) was used with the recommended parameters (2 5 7 80 10 50 2000) to identify TRs in the genome. Additionally, TRs from five ant species (Melters et al. 2013; Huang et al. 2016) were used with RepeatMasker to explore their enrichment in the genome of C. obscurior. Further details on annotation procedures are provided in the Supplemental Methods.

Estimation of recombination rate

To calculate local recombination rates, Cobs3.1 was used as a reference following the pipeline described above for read mapping, variant calling, filtering, and phasing.

We visually inspected all markers and discarded erroneous genotypes and those that indicated unlikely double-crossover events. These events were primarily identified as single-locus markers as well as blocks of two to three adjacent markers, tens to hundreds of base pairs long (Supplemental Figs. S18–S25), showing genotypic differences from immediately adjacent markers. Our reasoning to exclude these markers, was based on the expectations that double crossovers should be exceedingly rare due the suppression via chiasma interference (Sturtevant 1915; Muller 1916), in which the occurrence of one crossover disrupts the formation of another on the same chromosome. The effect of this interference remains robust over long distances, ranging from hundreds of kilobases to tens of megabases (Hillers 2004; Fowler et al. 2018). We also excluded regions where we could not identify crossovers unambiguously owing to low mapping accuracy and/or copy number variations (e.g., Supplemental Fig. S26).

Because all 84 unplaced scaffolds had insufficient markers (two or fewer), we only used markers on the 29 refined LGs. The final map was based on 57,532 high-quality markers that were clustered into 29 LGs (Supplemental Fig. S3) using MSTmap (Supplemental Table S9). Crossovers were defined as genotype changes between the maternal haplotypes, with most crossover tracks spanning >100 kb (Supplemental Figs. S22, S24, S25).

For local recombination rates, markers on each LG were ordered according to their physical distance, and we removed 223 markers showing discrepancies between their physical and genetic positions. Subsequently, nonrecombining markers were clustered according to their shared genetic positions, resulting in 1061 effectively recombining marker clusters. Local recombination rates were then inferred by calculating the ratio of the genetic (in cM) to the physical (Mb) distances between adjacent clusters. Average distance between crossover events was 149 kb (95% CI [122, 177]). To ensure that genomic regions had a proportional influence on estimates of recombination rate and to account for the nonuniform distribution of crossovers along the LGs, we calculated a weighted mean of recombination rate in 250 kb nonoverlapping windows. To this end, each window was assigned the recombination rate of the corresponding intercluster interval. When a window matched multiple intervals, the mean recombination rate of those intervals, weighted by their physical length, was assigned (Supplemental Code).

Estimating genetic diversity and differentiation

Nucleotide diversity (π) and Tajima's D were calculated in 250 kb nonoverlapping windows using two Pool-seq data sets of C. obscurior, a pool of 30 workers from Itabuna, Brazil (Errbii et al. 2021), and a pool of 50 workers from a ∼50 km distant location in Una, Brazil. Following the method of Errbii et al. (2021), we filtered paired reads, mapped them to Cobs3.1, and generated an indel-free mpileup file for each population using SAMtools (v.1.15.1) (Li et al. 2009). We used the PoPoolation (Kofler et al. 2011a) Variance-sliding.pl script to calculate π and Tajima's D. For π at synonymous (π_S) and nonsynonymous (π_NS) sites, we used the PoPoolation Syn-nonsyn-sliding.pl script, the mpileup files, and the codon and nonsynonymous codon length tables from PoPoolation.

For genetic differentiation (F_ST) between lineages of C. obscurior, we combined the alignment files of the two pools with a third alignment file of a pool of 16 workers from Tenerife, Spain (Errbii et al. 2021), to create a single mpileup file using SAMtools. We then converted this file into a synchronized file and used the fst-sliding.pl script under PoPoolation2 (Kofler et al. 2011b) to calculate F_ST in 250 kb nonoverlapping windows.

Calculation of absolute divergence

To explore genetic divergence between lineages of C. obscurior (d_XY) and the grandparental haplotypes (gp-d_XY), we calculated absolute divergence in 250 kb windows using scripts available at GitHub (https://github.com/simonhmartin/genomics_general).

Published paired-end short reads data from single workers of the Old World (Taiwan and Leiden) and New World (Itabuna and Una) lineages of C. obscurior (Errbii et al. 2021) were filtered and mapped to Cobs3.1 as described above. Note that reads from both lineages can be aligned to the reference genome with high confidence (Supplemental Fig. S27). Further details on SNP calling and filtering are given in the Supplemental Methods.

SV identification

To identify SVs, we adopted an approach by Mérot et al. (2023) using published paired-end short read data from 16 workers from five populations (Errbii et al. 2021). Reads were filtered and mapped to Cobs3.1, and read duplicates were removed using Picard's MarkDuplicates. SVs were identified using three different tools: Manta (v.1.6.0) (Chen et al. 2016); smoove (v.0.2.8) (https://github.com/brentp/smoove), which is based on LUMPY (Layer et al. 2014); and DELLY (v.0.8.1) (Rausch et al. 2012), all run with default parameters. To identify genomic inversions potentially affecting recombination in the F₁ queen, we used the three tools and followed the same approach. Subsequently, both sets of variants were filtered following protocols detailed in the Supplemental Methods.

For indels, we used the raw variant calls generated by GATK's HaplotypeCaller for the 16 workers. Indels were extracted and subjected to GATK's recommended hard filtering criteria (QD < 2.0; QUAL < 30.0; FS > 200, and ReadPosRankSum < −20.0). For further details on the filtering process and diversity calculation, see the Supplemental Methods.

Genome parameters estimation

BEDTools and BEDOPS (v.2.4.37) (Neph et al. 2012), along with the gene annotation, were used to calculate gene content. Overall and intergenic GC contents were calculated using the BEDTools nuc function in 250 kb nonoverlapping windows.

Estimating substitution rates

To explore the correlation between recombination and the efficiency of selection, we calculated omega (ω) for single-copy ortholog genes between C. obscurior and three ant species (Monomorium pharaonis [GCA013373865v2], Solenopsis invicta [UNIL_Sinv_3.0], and Ooceraea biroi [Obir_v5.4]) and the parasitoid wasp Nasonia vitripennis (Nvit_psr_1.1). Genomic and protein sequences, and gene annotations, were downloaded from Ensembl (release-54). To estimate the branch-wide rate ω and rates of synonymous (d_N) and synonymous mutations (d_S), we applied HyPhy's aBSREL (v.2.5.48; adaptive Branch-Site Random Effects Likelihood) (Smith et al. 2015) and the FitMG94 model, implemented in HyPhy. Further details on single-copy ortholog retrieval and ω calculation are given in the Supplemental Methods.

Gene expression bias

To investigate the relationship between local recombination rate and developmental gene expression bias, we calculated the plasticity index (Schrader et al. 2017), capturing gene expression bias across multiple caste/morph contrasts. High values indicate caste-specific gene expression, whereas low values suggest uniform expression across castes. RNA-seq data of worker, queen, wingless male, and winged male third instar larvae (n = 7 individuals each) were retrieved from the NCBI BioProject database (https://www.ncbi.nlm.nih.gov/bioproject/) under accession number PRJNA237579 (Schrader et al. 2015, 2017). For detailed information on read filtering and mapping, differential expression analysis, and computation of plasticity indice, see Supplemental Methods.

GO term enrichment analysis

To test whether genes found in regions of high recombination rates are enriched for specific functions, we ran a GO enrichment analysis. We defined these regions as windows with recombination rates falling within the top 10% (recombination rate > 22.3 cM/Mb). We used BEDTools’ intersect to extract 1541 genes located in these regions and ran topGO (v.2.46.0) (https://bioconductor.org/packages/release/bioc/html/topGO.html) to identify biological processes enriched in these regions. TE islands and regions affected by introgression were excluded from the analyses.

Correlation analyses and plotting

All correlations using the raw data and plots were produced in R (v.4.1.3) (R Core Team 2020). To explore the associations of TE content, sequence divergence, and epistasis with recombination rates per 250 kb, a generalized linear model was fitted using the glm function in R. We used Circos (Krzywinski et al. 2009) for graphical visualization in R. LGs 9, 10, 11, 17, and 21 were excluded from all correlation analyses, as introgression is evident in a significant fraction of these LGs. Further, LGs 26 and 29 were excluded because of low marker coverage (Supplemental Fig. S6).