Chromosome-level genome assembly reveals homologous chromosomes and recombination in asexual rotifer Adineta vaga.

Bdelloid rotifers are notorious as a speciose ancient clade comprising only asexual lineages. Thanks to their ability to repair highly fragmented DNA, most bdelloid species also withstand complete desiccation and ionizing radiation. Producing a well-assembled reference genome is a critical step to developing an understanding of the effects of long-term asexuality and DNA breakage on genome evolution. To this end, we present the first high-quality chromosome-level genome assemblies for the bdelloid Adineta vaga, composed of six pairs of homologous (diploid) chromosomes with a footprint of paleotetraploidy. The observed large-scale losses of heterozygosity are signatures of recombination between homologous chromosomes, either during mitotic DNA double-strand break repair or when resolving programmed DNA breaks during a modified meiosis. Dynamic subtelomeric regions harbor more structural diversity (e.g., chromosome rearrangements, transposable elements, and haplotypic divergence). Our results trigger the reappraisal of potential meiotic processes in bdelloid rotifers and help unravel the factors underlying their long-term asexual evolutionary success.

INTRODUCTION

Sexual reproduction and recombination are prevalent throughout eukaryotes, despite the substantial evolutionary costs such as the twofold cost of males or the cost of recombination that breaks up coadapted gene combinations (, ). Several eukaryotic species appear to have evolved adaptations that reduce these costs of males, for example, by producing males only facultatively as in cyclical parthenogens [e.g., Brachionus plicatilis ()] or by retaining a modified meiosis rescuing diploidy without fertilization by males [e.g., Diploscapter pachys ()]. Very few, however, appear to have renounced sex and recombination completely by abolishing males and meiosis. Theory predicts that, in the absence of recombination during meiosis, physical linkage among loci reduces the effectiveness of selection upon individual loci, resulting in a decreased rate of adaptation and the accumulation of mildly deleterious mutations (). Obligate asexuals are therefore suitable model systems to gain general insights into the long-term consequences of the lack of recombination and sexual reproduction.

Bdelloid rotifers are notorious ancient asexual animals. The longevity [>60 million years (Ma)] of the bdelloid rotifer clade and their diversity (>400 morphospecies) challenge the expectation that obligatory asexual animal lineages, in which recombination and outcrossing are absent, are evolutionary dead-ends. Historical observations (or lack thereof) had yielded a consensus that bdelloid rotifers do not produce male or hermaphrodite individuals () and that they are strictly parthenogenetic without any meiosis (, ). Moreover, the initial description of the structure of Adineta vaga genome, lacking colinear homologous scaffolds, was irreconcilable with meiosis (). A draft genome assembly of the closely related bdelloid species Adineta ricciae found colinearity between homologous regions but could not verify it at chromosome scale (), which was also the case for previous studies based on a handful of genomic regions (–). The presence or the absence of an ameiotic genome structure in bdelloids therefore remained unresolved, and a chromosome-scale assembly appeared critical.

Besides its asexual evolution, the bdelloid rotifer A. vaga also became a model species for its extreme resistance to desiccation, freezing, and ionizing radiation, with implications for space research (, ). Both prolonged desiccation, encountered in their ephemeral limno-terrestrial habitats, and ionizing radiation induce oxidative stress and massive genome breakage that A. vaga seems able to cope with, retaining high survival and fecundity rates while efficiently repairing DNA damage (–). Maintaining such long-term survival and genome stability following DNA fragmentation likely requires the use of homologous recombination (HR), at least in the germ cells. Given the supposed absence of homologous chromosomes in A. vaga (), the exact nature of their double-strand break (DSB) repair mechanism remains elusive.

Recent studies have provided evidence for recombination in bdelloid rotifers. These include a drop of linkage disequilibrium (LD) with increasing distance between genomic loci in A. vaga (), signatures of gene conversion (, ), heterozygosity levels within the range of those reported for sexual metazoans (, , ), and reports of allele sharing between bdelloid individuals from the wild (, –). While recombination likely takes place in bdelloid rotifers, its underlying mechanisms remain unknown. Recombination might theoretically occur in a mitotic or meiotic cellular context, involve short genomic regions or canonical chromosome pairing, and take place between homologous and nonhomologous (i.e., ectopic) loci. The interpretation of these recombination events has yet to be reconciled with the long-standing apparent absence of males and meiosis in bdelloid rotifers () and to account for their ubiquity in semiterrestrial habitats where frequent desiccation occurs, inducing DNA DSBs (, , ).

Here, we present a high-quality chromosome-level genome assembly of A. vaga. This is pivotal to tackle the contradictions between its putative ameiotic structure and the footprints of recombination, possibly associated to DSB repair and desiccation. We combined the use of short reads (Illumina), long reads (ONT and PacBio), and chromosome conformation capture data (Hi-C) with three assembly methods to successfully assemble the A. vaga genome. We provide the first telomere-to-telomere haploid and phased assemblies of a parthenogenetic lineage, paving the way to study genome evolution in an asexual clade. Using a newly developed and publicly available tool to detect horizontal gene transfers (HGTs), Alienomics, we annotated HGT candidates (HGTc) and confirmed that A. vaga has the highest number of HGTs across all animals. A. vaga has a diploid genome made of six pairs of homologous chromosomes, refuting the ameiotic structure previously described for this genome () and challenging the complete absence of meiosis in one of the most notable asexual animal clade. In addition, by observing large tracks of losses of heterozygosity (LOHs), we show that large-scale recombination between homologous chromosomes occurs in A. vaga. The possibility of chromosome pairing in A. vaga, during either a meiotic-like or mitotic process, allows a reinterpretation of the signatures of LD decay and allele sharing. Until now, the lack of chromosome-scale assemblies of parthenogenetic genomes hampered the investigation of the impact of meiosis, recombination, outcrossing, or their absence on entire genomes. Moreover, characterizing homologous chromosomes as potential templates for DNA repair through HR in A. vaga is an important landmark in understanding of bdelloid extreme resistance. This high-quality genome assembly of A. vaga (AV20) is also timely for comparative biology within rotifers and protostomians, extending the list of chromosome-level genomes in overlooked phyla.

RESULTS AND DISCUSSION

A diploid genome with a tetraploid past

Distinct independent genome assembly procedures, relying on different assumptions regarding ploidy levels [Bwise (), NextDenovo (), and Falcon ()], were first used on a combination of short and long sequencing reads. These assemblies were then scaffolded using Hi-C data and instaGRAAL (), revealing similar chromosome-level assemblies and genome size estimations consistent with flow cytometry measurements (Fig. 1A and figs. S1 and S2). All pairwise alignments of the three independent assemblies (referred to as “phased” without ploidy assumption, “haploid,” and phased “diploid”; see Fig. 1A) confirmed chromosome-level synteny and converged toward an identical genome structure, with the six longest scaffolds from the haploid assembly (hereafter named “AV20”) being each colinear to exactly two long scaffolds from the phased assembly (Fig. 1B; see also figs. S3 to S5). To validate these assemblies, we performed fluorescence in situ hybridization (FISH) analyses with three pairs of fluorescent probe libraries complementary to separate parts of chromosomes 2, 5, and 6 from the AV20 assembly (Fig. 1B, right). For each pair of probes (one green and one red), two individual chromosomes were labeled with little or no overlap between both signals (Fig. 1C). Chromosome painting on the karyotype of 12 chromosomes of A. vaga () was consistent with our chromosome-scale assemblies, showing that the A. vaga genome is diploid, being composed of six pairs of colinear homologous chromosomes.

Fig. 1.

The genome structure of A. vaga is diploid.

(A) Outline of the three genome assembly approaches underlined by different assumptions on genome ploidy, with median read coverage for all sequencing technologies indicated on the left and estimated with respect to the AV20 haploid genome assembly. The haploid genome size estimate of A. vaga obtained by flow cytometry (under the assumption that the genome is diploid), as well as the summary statistics of the genome assemblies, is given. The number of chromosomes corresponds to the number of scaffolds longer than 10 Mbp. Ploidy levels of assemblies are indicated by the KAT plots of k-mers distribution (first and second peaks correspond to heterozygous and homozygous k-mers, respectively; red and purple indicate haploidy and diploidy, respectively). (B) Circos plot of the pairwise colinearity between the haploid AV20 and the phased Bwise genome assemblies, depicted by colored links and obtained using nucmer. Synteny blocks within AV20 genome (between homoeologous copies) are depicted as gray links and were obtained using MCScanX. Coverage along scaffolds of both AV20 and the phased assembly is depicted as gray histograms and was computed on the basis of Illumina reads from sample GC047403. Thin black bars on the scaffold ideograms correspond to ribosomal RNA (rRNA) genes. Schematic position of the FISH probe libraries on chromosome pairs 2, 5, and 6 is indicated on the corresponding AV20 chromosomes. (C) Fluorescence microscopy images of the 12 chromosomes of A. vaga [4′,6-diamidino-2-phenylindole (DAPI) staining] with chromosome pairs 2, 5, and 6 highlighted by oligo painting using the FISH probe libraries depicted in (B).

We compared our new AV20 assembly to the previously published draft genome assembly [hereafter named “AV13” ()]. None of the previously described colinearity break points and palindromes were retrieved in the new AV20 genome, indicating that these were likely assembly artifacts resulting from erroneous scaffolding of phased haplotypes (figs. S6 and S7). Chromosome-level colinearity, albeit weaker than between homologous chromosomes, was also observed between pairs of homoeologous (or ohnologous) chromosomes in the AV20 genome, a signature confirming the previously reported paleotetraploidy of A. vaga (, , ) (gray links in Fig. 1B). The three chromosome pairs 1, 2, and 3 are homoeologous to the three pairs 4, 5, and 6, respectively. A. vaga is thus a diploid, paleotetraploid species in which the level of synteny between homoeologous chromosomes is high. Notably, 30.8% of the genes have homoeologous copies within conserved synteny blocks (see Materials and Methods) and with an average nucleotide divergence of about 13% (fig. S8).

Recombination between homologous chromosomes causes loss of heterozygosity

The discovery of homologous chromosomes in the oldest known asexual animal clade represents a major shift for studies of ancient asexuals and leads us to reconsider the possibility of HR in A. vaga. One potential genetic consequence of recombination between homologous chromosomes is loss of heterozygosity (LOH). We measured and compared heterozygosity along the chromosomes of three A. vaga samples cultured from the same ancestral laboratory strain that never underwent stresses causing recombinogenic DSBs and that were sequenced at three distinct time points (2009, 2015, and 2017; Fig. 2A and table S1). Mean single-nucleotide polymorphism (SNP) heterozygosity (i.e., divergence between homologous chromosomes) was around 1.7% [horizontal line in Fig. 2A; similar to previous reports (, )]. We observed large regions (from 100 kb to 4.5 Mb) that were fully homozygous, except for a few SNPs, in specific isolates, while heterozygous in others (numbered tokens in Fig. 2A). Note that a few homozygous tracks were associated with coverage variation and could have been caused by a hemizygous deletion (when coverage drops by approximately 50%, e.g., event 5 in Fig. 2A) or by the high density of repeated sequences (e.g., event 13 in Fig. 2A). Given the genealogy of these laboratory lines, we argue that the large homozygous tracks that are associated with homogeneous median coverage are signatures of allelic recombination events causing LOH (Fig. 2B).

Fig. 2.

Heterozygosity dynamics in A. vaga.

(A) Heterozygosity and coverage distributions of three independent A. vaga samples from the same laboratory strain along the six chromosomes. Samples are labeled by the date of the extraction of their DNA (i.e., 2009, 2015, and 2017). Data from 2009 were used to assemble the previous version of A. vaga genome (). Lines indicate short read coverage (normalized), and filled areas indicate the percentage of heterozygosity (y axis). Chromosome lengths (x axis) are in Mb. Mean SNP heterozygosity (1.7%) is depicted by the horizontal black line. (B) Schematic reconstruction of heterozygosity evolution among three samples from the same initial A. vaga laboratory strain. Note that each sample had its own independent evolution, and the exact sequence and timing of loss of heterozygosity (LOH) events is unknown. LOH events noted with a small asterisk might correspond to deletions given the drop of coverage associated with the absence of heterozygosity.

These LOHs appear to accumulate through time as some are shared by two samples (e.g., event 12 in Fig. 2), while others appeared in only one of these two samples (e.g., event 2 in Fig. 2). Note that no ancestral LOH was found that would be shared by all of the strains. This is likely because large LOH events increase the chance of exposing recessive deleterious mutations and are thus likely to be selectively eliminated in nature, maintaining the relatively homogeneous heterozygosity level in the ancestral laboratory strain (Fig. 2A). The presence of LOH tracks on all six chromosome pairs in the three laboratory samples over a relatively short period of time (i.e., several years; Fig. 2B) might be due to the culturing conditions, allowing possible bottlenecks and relaxed selection. Recombination occurring along the entire chromosomes, instead of being restricted to the telomeres only (), invalidates the hypothesis that an Oenothera-like meiosis underlies their reproductive mode [in agreement with a recent study ()]. Overall, these LOH tracks combined with the recently reported LD decay () represent a clear footprint that molecular processes involving recombination between homologous chromosomes occur in the germ line of A. vaga.

Recombination could be accidental or programmed

Theoretically, recombination between homologous chromosomes resulting in inheritable LOH could occur in the germ line during mitotic repair of accidental DSBs or when handling programmed DSBs during meiosis [potentially induced by the Spo11 protein (, )]. DSBs can be repaired by different recombination pathways, but LOH of large chromosome regions without coverage reduction (e.g., events 1, 6 to 9, 11, and 12; Fig. 2) can primarily arise from two processes, break-induced replication (BIR) and the formation of crossing-over (CO). BIR is a mechanism of mitotic recombination characterized by replication fork progression over hundreds of kilobases on the repair template (). When involving allelic loci, it causes LOH in the segment, extending from the break point site until the end of the chromosome. If a double BIR occurs, switching templates from the homologous chromosome back to the sister or the original chromatid, an LOH tract, possibly long, that does not encompass the telomere is produced (). Such LOH could also be generated by the recombinational repair of one or two DSBs, respectively, leading to CO (i.e., a reciprocal genetic exchange between chromosomes). Compared to BIR, CO is, however, a minor pathway in mitotically cycling cells () that preferentially takes place between sister chromatids and therefore remains genetically silent ().

Alternatively, programmed DSBs during meiosis can produce large LOH tracks by favoring CO formation between homologous chromosomes (). LOH signatures in the A. vaga genome could therefore be acquired through meiotically induced recombination instead of during mitosis. Several mechanisms of meiotic parthenogenesis, globally referred to as automixis, have been described in various species such as in Daphnia pulex (), Artemia parthenogenetica (), or Apis mellifera capensis (). If automixis occurs in A. vaga, the heterozygosity patterns observed here (Fig. 2), in which maternal heterozygosity is conserved along chromosomes due to the nonsegregation of homologous chromosomes, while large LOH tracks (likely counterselected in nature) could result from their CO recombination, are genetically equivalent to what is referred to as central fusion automixis (). Nevertheless, no cytological evidence of any meiotic process has been described so far in bdelloid rotifers. Whether recombination is a key feature of the reproductive mode of bdelloids (through programmed DSBs during a modified meiosis) or whether it is mainly driven by desiccation resistance mechanisms (through accidental DSB repair in the germ line), or both, remains an open question. Whichever mechanism is involved, recombination likely plays a major role in the long-term evolution of the A. vaga genome.

Dynamic subtelomeric regions

We found a low amount of transposable elements (TEs) in A. vaga (Fig. 3 and fig. S9). By combining two approaches to annotate both TE-like elements, including repeated sequences, as well as canonical TEs (i.e., the EDTA and REPET pipelines), we detected 6.6% of TE-like elements and 1.9% of canonical TEs, predominantly located in subtelomeric regions (Fig. 3). In addition, rotifer-specific telomeric repeats [i.e., (TGTGGG)n ()] were detected at the extremities of every scaffold of the AV20 assembly, indicating that they correspond to telomeric and subtelomeric regions and that AV20 reached a chromosome-level assembly (fig. S10). Most consensus TE sequences were found at low copy numbers (i.e., 96% of canonical TE consensus sequences are present in less than six copies in AV20; see fig. S9). Notably, terminal inverted repeat (TIR) DNA transposons (i.e., Class-II TIRs) were quantitatively dominant (48% of all TEs) among the low amount of TEs in A. vaga genome (figs. S9 and S11). These results are in line with previous studies of TEs in bdelloids (, , , ). Using sequence similarity between a TE copy and its consensus as a proxy for how recent this copy is, we found that Class-II TIRs and Class-I LINEs (long interspersed nuclear elements) and LTRs (long terminal repeats) had high average similarity to their consensus sequences, suggesting that they have been at least recently active in A. vaga genome (fig. 12). Investigating putative endogenous viral elements (EVEs) in A. vaga revealed very few viral-like sequences (i.e., 94 loci scattered along the six chromosomes; Fig. 3), with potential donor candidates belonging to the group of large double-stranded DNA viruses. However, none of these EVE candidates had definitive viral origins, as their similarity was not restricted to viral sequences and there was no conservation of viral gene synteny.

Fig. 3.

DNA content of the haploid AV20 genome assembly.

Conservation of synteny of HGTc is depicted by colored links between the six chromosome pairs. Violet links correspond to the synteny block of HGTc between nonhomoeologous chromosomes. Localization of alternative haplotigs, removed before genome scaffolding, is depicted by gray links. Distribution of repeated elements, genes, HGTc, and viral-like sequences are depicted in green, blue, red, and black bars, respectively. Ancient HGTc hotspots are indicated by pink stars.

Syntenic HGTc among nonhomoeologous chromosomes are visible, mainly in subtelomeric regions (violet links on Fig. 3) and may suggest chromosomal rearrangements. Subtelomeric regions are also the regions on which almost all divergent haplotypes (i.e., haplotigs corresponding to uncollapsed haplotypes during genome assembly process) were located (gray links on Fig. 3). Overall, these subtelomeric regions in A. vaga not only are enriched in canonical TEs, TE-like elements, HGTc, and viral-like sequences but also retain a higher haplotypic divergence (i.e., uncollapsed haplotigs) and most chromosomal rearrangements. When accounting for coding sequences only, no distinct increase or decrease of heterozygosity could, however, be observed in subtelomeric regions (fig. S13). At this stage, it is therefore unclear whether HR rate covaries with distance from telomeres in A. vaga. Nevertheless, our results suggest that subtelomeric regions seem more prone to chromosomal rearrangements, incorporation of foreign DNA (TEs and HGTs), and structural variations such as putative allelic deletions (see LOH events 3 and 10 in Fig. 2), evolving faster than the rest of the genome.

HGTs in A. vaga genome

The acquisition of foreign DNA has been hypothesized to play an important role in bdelloid evolution (). HGTs could circumvent some deleterious effects of the lack of sexual outcrossing, and the occasional integration of foreign DNA could trigger adaptation (, , , , ). No automated tool existed to detect HGTs; therefore, we developed Alienomics, an innovative pipeline to detect both HGTs and contaminants in a genome assembly. Alienomics combines several genomic parameters such as gene taxonomy, GC (guanine-cytosine) content, and sequencing depth, also taking into account gene integration into the genome using synteny and expression data, to detect HGTs from nonmetazoan species. In contrast with the overall low amount of TEs, many HGTc (2679, about 8.3% of all genes) were detected in the A. vaga AV20 genome assembly, confirming previous reports of the highest HGT content among metazoans (, , , ). HGTc were enriched in subtelomeric regions as previously reported (), although many HGTc were distributed along the chromosomes and two visible local hotspots were detected outside of the subtelomeric regions (pink stars in Fig. 3). One HGT hotspot was associated with a slight increase of interstitial telomeric repeats (fig. S10) that could represent a signature of an ancient chromosome fusion. Overall, the heterogeneity in HGTc density between subtelomeric regions and the rest of the genome could be explained either by varying rates of HGT incorporation along the chromosomes or by varying successful integration of HGTs within the genome through selection.

Using both MCScanX and Alienomics outputs, we determined that 257 foreign genes (9.6% of all HGTc) had conserved their synteny across homoeologous chromosomes, including the HGTc hotspots notably visible on homoeologous chromosomes 1 and 4 (stars on Fig. 3). These horizontal transfer events therefore occurred before the ancestral tetraploidization of modern bdelloids. This amount of ancient HGTc is, however, likely underestimated as any loss or translocation of an ancient HGTc copy would break the ancestral synteny. When looking specifically at HGTc that occurred before the tetraploidization, we observed an enrichment of genes involved in DNA recombination and DNA ligation being part of the DNA DSB repair pathways, among other enriched functional categories (see table S2). These HGTc might have allowed bdelloids to resist and overcome massive DNA breakage. When analyzing all HGTc, we found that they are enriched in genes involved in oxidation reduction and carbohydrate metabolic processes () as well as in the response to nitrosative stress (see table S2). Acquisition of HGTs might therefore be central in their resistance to extreme desiccation and toward more efficient homeostasis. Overall, these results are in line with previous studies, suggesting that HGTs have been continuously acquired within bdelloid rotifers, even before their tetraploidization (, ). However, if bdelloids have the same low rate of HGT acquisition from other individuals of the same species than from nonmetazoans (12.8 HGT/Ma), HGT is possibly insufficient to compensate for the plausible lack of outcrossing in bdelloid rotifers (). Nevertheless, a high rate of HGT acquisition from distinct species might be deleterious for A. vaga, because recombination between homologous chromosomes appears to be central to the maintenance of heterozygosity and/or genome structure in this species.

Reasoning on bdelloid rotifers’ reproductive mode

Bdelloid rotifer species are both supposedly devoid of males and prone to integrate foreign DNA (through HGTs) into their genome. In this context, several reports of allele sharing between bdelloid individuals sampled from the wild triggered a debate whether they could exchange genetic content at all and whether this might be done through HGT or through sexual reproduction (, –, –). At a first glance, showing that homologous chromosomes exist and recombine in A. vaga could be viewed as a support to the hypothesis that bdelloids might undergo meiotic sexual reproduction (). However, this hypothesis has yet to be reconciled with the absence of both males and canonical meiosis in bdelloid rotifers, and here, we speculate on the mechanisms of HR in A. vaga. The three A. vaga lineages analyzed here (Fig. 2) were kept in hydrated conditions, leaving few opportunities for desiccation-induced, accidental DNA DSBs. Moreover, a much lower heterozygosity than for A. vaga has been observed in two obligate aquatic bdelloid rotifer species (i.e., in the genus Rotaria in which the upper limit of homologous divergence ranged between 0.033 and 0.075), also described as asexual and never experiencing desiccation. Both these observations favor the hypothesis that HR in bdelloids could be caused by programmed DNA DSBs during a meiotic-like process. Frequent and programmed recombination would cause LOH in hydrated A. vaga (Fig. 2) and would have lowered heterozygosity in the obligately aquatic Rotaria species.

Whatever the underlying mechanism, the observed recombination signatures in bdelloid rotifers are compatible with the three hypotheses proposed to explain the previous reports of allele-sharing patterns in bdelloid rotifers: (i) allele sharing may be an artifact due to undetected contamination between cultures, either during colony culture itself or during sample preparation for sequencing (–); (ii) allele sharing is the result of horizontal genetic transfers between bdelloid individuals through unknown molecular mechanisms, possibly associated with desiccation (but not for the nondesiccating species) and potentially linked to the high propensity of bdelloids to retain nonmetazoan genes into their genomes (, , , , , ); and (iii) allele sharing is caused by cryptic sexual reproduction (), with sex events being rare enough that males, sperm, fertilization, and meiosis were never observed, but sufficiently frequent to leave a distinctive footprint in every population sample studied so far (, , , ). The mechanism behind the observed signatures of genetic exchanges between bdelloid individuals remains puzzling, and therefore, the significance of outcrossing in this asexual lineage is still unclear. We anticipate that the chromosome-level genome assembly of A. vaga presented here will stimulate future population genomic studies that will help to determine the cause of these allele-sharing patterns.

Long-term asexual evolution

This high-quality telomere-to-telomere assembly firmly establishes A. vaga as a model system to study long-term asexual evolution. Homologous chromosomes are present in the bdelloid species A. vaga and might well occur in all bdelloid rotifers, as colinear pairs of sequenced fosmids were found in the two distinct bdelloid species A. vaga and Philodina roseola, with each colinear pair in one species resembling a colinear pair in the other species (). The observed long LOH tracks indicate the existence of long-range HR (Fig. 4). Whether this occurs during a meiotic-like parthenogenetic mode of reproduction or in a mitotic context during frequent repair of accidental DNA DSBs remains speculative. Recombination (mitotic or meiotic) could increase the rate of gene conversion in asexual lineages, a signature previously observed in A. vaga (). Gene conversion, particularly when slightly biased, could correct deleterious mutations and reduce the rate of clonal deterioration () or even speed up the fixation of beneficial mutations (). However, besides signatures of LOH via allelic recombination, we also observed LOH via deletions in the genome of A. vaga. The random accumulation of LOH events could expose deleterious recessive mutations in asexuals through loss of complementation (). Our new chromosome-level genome assembly of asexual A. vaga therefore provides a critical tool to be able to evaluate the relative benefits of these recombination events on their long-term evolution and paves the way for studies on genome dynamics in A. vaga.

Fig. 4.

Schematic representation of the karyotype of A. vaga.

The 12 chromosomes correspond to six pairs of homologous chromosomes (i.e., a diploid genome structure) sharing the same color. Ancestral genome hybridization (or whole-genome duplication) led to the existence of pairs of homoeologs (or ohnologs), represented in different but similar colors. Gray blocks linking homologous chromosomes and their corresponding numbered tokens depict LOH events produced by HR (see also Fig. 2). The length of chromosomes (in Mb) is indicated by the scale on the right.

In general, asexual populations suffer from the absence of gene shuffling with other individuals, and the long-term evolutionary success of bdelloid rotifers in the absence of outcrossing therefore remains puzzling. It is important to try to discriminate between the consequences of the two aspects underlying sex: recombination and outcrossing. Theoretical work on population genetics showed that selection could be at least as efficient in automictic lineages than in sexuals under certain circumstances [e.g., effective population size and recombination rates ()]. It is therefore conceivable that the combination of potentially large populations, a relatively high level of heterozygosity (or mutation rate), and specific recombination rates might explain how A. vaga maintains a delicate balance between losing and accumulating heterozygosity and how it adapts and persists in the long term. Unfortunately, critical knowledge about bdelloid biology (e.g., quantitative estimates of mutation and recombination rates) is still lacking to determine whether they might circumvent the lack of outcrossing through a genetic equivalent of automixis. Outcrossing through sexual reproduction might speed up adaptation by allowing the combination of independently evolved alleles within the same individual, but might not be essential for bdelloid rotifers, especially if a high frequency of HGT is also taking place. Despite presenting the highest amount of HGTs among animals, our results also suggest that bdelloid rotifers might have to balance the acquisition rate of HGTs, a source of functional novelties, with the maintenance of faithful homology between chromosomes for HR. Overall, our work reinforces the hypothesis that recombination is critical for lineage longevity. Ancient asexual animals without a minimal rate of recombination, programmed through meiotic processes and/or accidental through their lifestyle, might not exist at all.

MATERIALS AND METHODS

Complete description of Materials and Methods can be found in the Supplementary Materials.

Genome sequencing and assembly

Cultures of the A. vaga AD008 laboratory strain were processed to obtain the following sequencing data: about 350× coverage of WGS 250–base pair (bp) paired-end Illumina reads, 200× coverage of PacBio RSII long reads, 125× coverage of ONT long reads, and 75-bp paired-end Illumina reads of Hi-C libraries. Three independent genomes were assembled using Bwise () (on Illumina short reads), NextDenovo () (on ONT long reads), and Falcon () (on PacBio long reads). Uncollapsed haplotypes in the ONT-based assembly were detected and discarded using purge_dups. This NextDenovo assembly was then polished on the basis of Illumina short reads and PacBio long reads using HyPo and referred to as AV20. All assemblies were scaffolded with instaGRAAL () using Hi-C data. Ploidy level and genome size were confirmed through interpretation of k-mers spectra using KAT () as well as by synteny analyses between the three genome assemblies presented here and the previously published AV13 genome assembly using MCScanX () (synteny blocks detection), nucmer () (genome alignment), and D-GENIES () (dot plot visualization). Flow cytometry measurements further confirmed this genome size, and FISH was used to confirm the ploidy level. For this, we designed three pairs of oligos on three different chromosomes using the OligoMiner pipeline (), and embryos fixed at one-cell stage in methanol were permeabilized (i.e., mechanical squash, Triton, and saponin). Oligo probes were added to embryos followed by DNA denaturation (at 92°C) before hybridization of probes and DNA (at 37°C). Karyotype images were obtained using confocal imaging.

Genome annotation

TE-like elements and canonical TE consensus were predicted from the AV20 genome assembly using EDTA () and TEdenovo pipeline (, ) (part of the REPET pipeline). TE-like elements were then annotated using TEannot. Genes were annotated using funannotate () as described below. First, repeated elements were masked using BEDTools (). Second, part of the available RNA sequencing (RNA-seq) reads were mapped onto the genome, while the other part of RNA-seq reads were used to produce a de novo assembled transcriptome that was subsequently aligned onto the genome as part of the PASA pipeline (). Third, PASA annotations, the de novo assembled transcriptome, the metazoan BUSCO database, and the protein UniProt database were used as inputs within funannotate predict function. This function first produced ab initio gene annotations using GeneMark-ES, which were then used along with transcripts and protein data to train Augustus to generate a second set of gene annotations. Fourth, funannotate used EVidenceModeler as a weighted approach to combine gene annotations from PASA, GeneMark, and high-quality annotations from Augustus into a single, integrated, gene annotation set. Last, InterProScan 5 () was used to produce functional annotations to these genes. These functional annotations, along with the BUSCO metazoan database, were then used as input for the funannotate annotate function with default parameters to produce the final functional annotation for the predicted genes. Alienomics, a newly designed pipeline, was used to detect HGTs. This approach combines GC content, coverage, sequence similarity, taxonomic information, expression level, and synteny information to detect both HGTc (i.e., alien genes integrated into host scaffolds) and potential contaminants (i.e., alien genes present on alien scaffolds). Note that our approach can only detect transfers from alien source outside of a given clade (here, metazoans). Viral-like genes were detected by performing a DIAMOND BLASTX search on AV20 scaffolds using all viral proteins extracted from the nr database of the National Center for Biotechnology Information (NCBI) (February 2020) to the exception of Retroviridae and Hepadnaviridae.

Genome analyses

Coverage along AV20 scaffolds was computed using read mapping with BWA-MEM. These mappings were used for measuring heterozygosity in three samples (i.e., GC047403, BXQF, and ERR321927) by genotyping them using GATK (HaplotypeCaller function with -ERC GVCF option). The resulting gvcf files were combined (CombineGVCFs function) and were then jointly genotyped (GenotypeGVCFs function). The variants were then filtered to only retain SNPs using custom bash and perl scripts. Divergence between homoeologous chromosomes was assessed by producing a self-alignment of AV20 genome using nucmer, which was then filtered to only retain genomic alignments between homoeologous regions ranging from 500 to 10,000 bp. MCScanX was used to detect synteny among HGTc, and custom scripts were used to detect strictly homoeologous HGTc synteny blocks stemming from the paleotetraploidization of bdelloids (i.e., ancient HGTc). Gene Ontology (GO) terms from the functional annotation were extracted for the 2422 recent HGTc; the 257 ancient HGTc were then compared to the entire gene set of the AV20 genome containing 32,378 proteins. Enrichment analyses were performed using the topGO package with a Fisher test and the “elim” algorithm ().

Reappraisal of the AV13 genome assembly

ONT reads were trimmed using porechop () and were mapped onto the previous AV13 genome assembly using NGMLR (). The AV20 and AV13 genome assemblies were aligned together using Sibelia (). This genome alignment was used to localize putative break points, which were then inspected using Tablet () to evaluate whether ONT reads confirmed previously reported break points in the AV13 genome assembly. Synteny blocks from the Sibelia alignment between AV20 and AV13 were analyzed using a custom perl script (available at https://github.com/jnarayan81/huntPalindrome) to localize putative palindromes, and circos plots of these locations were inspected to evaluate the existence of palindromes in the AV20 genome.