Genomic features of asexual animals

Evolution under asexuality is predicted to impact genomes in numerous ways, but empirical evidence remains unclear. Case studies of individual asexual animals have reported peculiar genomic features which have been linked to asexuality, including high heterozygosity, a high abundance of horizontally acquired genes, a low transposable element load, and the presence of palindromes. However, it is unclear whether these features are lineage-specific or general consequences of asexuality. We reanalyzed published genomes of 24 asexual animals and found that not a single genome feature is systematically replicated across a majority of these species, suggesting that there is no genomic feature characteristic of asexuality. We found that only asexuals of hybrid origin were characterized by high heterozygosity levels. Asexuals that were not of hybrid origin appeared to be largely homozygous, independently of the cellular mechanism underlying asexuality. Overall, despite the importance of recombination rate variation for understanding the evolution of sexual animal genomes, the genome-wide absence of recombination does not appear to have the dramatic effects which are expected from classical theoretical models. The reasons for this are probably a combination of lineage-specific patterns, impact of the origin of asexuality, and a survivor bias of asexual lineages.


Introduction 34
Sex: What is it good for? The reasons for why most eukaryotes take a complicated detour to 35 reproduction, when straightforward options are available, remains a central and largely 36 unanswered question in evolutionary biology 1,2 . The species that use asexual reproduction 37 as their sole form of replication typically occur at the tips of phylogenies and only few of them 38 have succeeded like their sexually reproducing counterparts 3 . In other words, most asexual 39 lineages may eventually be destined for extinction. These incipient evolutionary failures are 40 however invaluable because by understanding the evolutionary fate of asexual species, 41 something may be learned about the adaptive value of sex. 42 43 An accumulating number of studies have sequenced the genomes of asexually reproducing 44 animals, often with the aim of identifying features that would distinguish them from sexual 45 species (Figure 1). In asexual animals, females produce daughters from unfertilized eggs 46 via so-called thelytokous parthenogenesis (hereafter asexuality) 4 . Asexuality is predicted to 47 have many consequences for genome evolution, since gamete production via meiosis and 48 the restoration of somatic ploidy levels via fertilization no longer take place. Predicted 49 consequences include for example the accumulation of deleterious mutations 5-7 , as well as 50 changes in heterozygosity levels 8,9 and transposable element (TE) dynamics 10 . Some 51 predictions have been tested without genomic data, using a handful of housekeeping genes. 52 However, conclusions based on such small and non-random subsets of genomes can lead 53 to erroneous conclusions 11 . With the advent of high-throughput sequencing it is possible to 54 evaluate classical predictions of asexuality at the genome scale, and furthermore to test new 55 predictions, such as the accumulation of palindromes (see below), which could not be 56 studied with single gene approaches. In the present study, we compare the published 57 genomes of 24 asexual animal species (Figure 1) to assess whether we can identify any key 58 features characteristic of asexual animals. The 24 species comprise four species of bdelloid 59 rotifers, a group that likely persisted and diversified in the absence of canonical sex for over 60 40 million years 12 . Bdelloids have thus far overcome the predicted dead-end fate of 61 asexuality, which raises the question of what mechanisms protect them from extinction, and 62 whether these mechanisms are visible in specific characteristics of their genomes. 63 64 Because the predicted consequences of asexuality are strongly affected by how asexuality 65 evolved from the sexual ancestor (Box 1) as well as by the cellular mechanisms underlying 66 asexuality (Box 2), we include biological differences among asexual species in our 67 comparisons. For example, some asexual species have evolved via hybridization (Box 1), 68 which generates high heterozygosity and can result in increased activity of transposable elements [13][14][15] . In such instances, it can be difficult to disentangle consequences of 70 hybridization from consequences of asexuality. Similarly, some cellular mechanisms 71 underlying asexuality involve meiotic divisions, with a secondary restoration of somatic 72 ploidy levels, while others do not. In the former case, heterozygosity in the asexual species 73 is expected to decay rapidly, while in the latter case, it could be maintained or even increase 74 over time (Box 2). Finally, because genome studies differed in their focus and in the 75 methods used, we reanalyzed published genomes with standardized approaches. Whenever 76 possible, we conducted quantitative comparisons between groups of asexual species. 77 However, for interpretation, it is important to consider that the available genomes are not a 78 random nor representative sample of asexual animals, and that not all of these genomes 79 reflect evolutionarily independent events. 80 81 We uncovered a number of unusual features in the genomes of asexual animals, including 82 extreme loads of transposable elements and highly asymmetric divergence among 83 haplotypes in polyploid species of hybrid origin. However, none of these were systematically 84 replicated across even a majority of analyzed species, let alone all of them, suggesting that 85 there is no universal genomic feature specific to asexual species. We found that a hybrid 86 origin of asexuality was the most important factor affecting heterozygosity, with potential 87 effects of asexuality being masked by effects of hybrid ancestry. Unexpectedly, asexuals 88 that are not of hybrid origin are largely homozygous, independently of the cellular 89 mechanism underlying asexuality. color of the circle indicates the cellular mechanism of asexuality and the number inside the In sexual species offspring is generated through the fusion of male and female gametes. In 131 asexuals, females generate diploid (or polyploid) offspring from unfertilized oocytes via 132 different cellular mechanisms. The cellular mechanism used is predicted to affect genome 133 evolution and especially heterozygosity levels. For details see 4,42 . 134 Mitotic asexuality (Apomixis). Under mitotic asexuality, no ploidy reduction occurs and 135 offspring are clones of their mother. 136 Meiotic asexuality (Automixis). Under meiotic asexuality, meiotic divisions occur partially or 137 completely, but somatic ploidy levels are maintained via different mechanisms. Some of 138 these mechanisms have similar genomic consequences as mitotic asexuality, even though 139 meiosis is involved (for example, endoduplication in hybrid asexuals results in offspring that 140 are clones of their mother. Such mechanisms are often referred to as "functionally mitotic" 141 (or functionally apomictic), especially when the cellular mechanisms are not known in detail 142 but genotyping data suggest that offspring are clones of their mother. 143

Endoduplication.
A duplication of the entire chromosome set occurs before normal 144 meiosis, during which ploidy is reduced again. If recombination occurs between 145 identical chromosome copies rather than between chromosome homologs, 146 endoduplication produces offspring that are clones of their mother. 147 Overview of species and genomes studied 163 We reanalyzed the published genomes of 24 asexual animal species with the aim of 164 identifying general genomic signatures of asexualilty. The 24 species correspond to at least 165 16 independent transitions to asexuality and cover a broad taxonomic range, including 166 chordates, rotifers, arthropods, nematodes and tardigrades. In addition to covering this 167 taxonomic range, the asexual species vary in the cellular mechanisms underlying asexuality, 168 in the mechanisms that caused the transition to asexuality, as well as in other biological 169 aspects (Figure 1, Supplementary Tables 1 & 2). This variation allows us to assess 170 whether asexuality generates universal genomic signatures independently of species-171 specific traits. 172

173
The cellular mechanisms underlying asexuality were studied in 20 of the 24 species. Eight of 174 them use mitotic asexuality, while the 12 remaining species use different types of meiotic 175 asexuality (Figure 1). All but one of the eight species with mitotic asexuality are polyploid, 176 the amazon molly P. formosa being the only diploid studied. Conversely, all but one species 177 with meiotic asexuality are diploid. This is expected given that polyploidy can generate 178 problems during meiosis (reviewed in 45 ). Nevertheless, the nematode Panagrolaimus sp. is 179 characterized by both meiotic asexuality and triploidy 30 (see Supplementary Table 1 for 180 details). 181 182 Information on how asexuality evolved is available for 15 of the 24 sequenced species 183 ( Figure 1). A hybrid origin has been suggested for ten of these. Endosymbionts are the most 184 likely cause of asexuality in four species (the springtail, both wasps and the thrips), and 185 spontaneous mutation in one (the ant). Across the 24 species, hybrid origin is correlated with 186 polyploidy. Six of the 11 polyploids in our sample are of hybrid origin, while for the five others 187 a hybrid origin has thus far not been suggested, but is supported by our results (see below). 188 It is important to note however that there are many polyploid asexual animals that are not of 189 hybrid origin, including several well studied asexual species such as the New Zealand meaning that their effect size as well as the power to detect them increases in old asexual 196 lineages. However, estimating the age of asexual lineages is difficult and always associated 197 with large uncertainties 46,47 . We therefore did not include quantitative comparisons among 198 asexuals with respect to their age. However, because our set of species comprises asexuals 199 believed to be 'ancient' (i.e., several million years old, see Supplementary Table 1) The prediction that deleterious mutations accumulate more rapidly in asexual than sexual 210 lineages has been tested in over twenty groups of different asexual species (reviewed in 51 211 and with three additional studies published since 16,21,52 ), with results generally supporting 212 the prediction. However, in only eight studies were the tests conducted genome wide, while 213 tests in the remaining studies were based on one or a few genes only. Note that four 11,52-54 214 of these studies were based on transcriptomes and are therefore not included in our 215 systematic reanalysis. Among the genome wide tests, results are much more mixed than 216 among the 'single or few genes' studies, raising the question whether the latter are 217 representative of the genome as a whole. Specifically, only two of the eight genome-wide 218 studies found support for deleterious mutation accumulation in asexuals 52,53 . However, two 219 studies found that sexual taxa experienced more deleterious mutation accumulation than 220 asexual taxa 11,19 while the four remaining ones found no differences between sexual and 221 asexual taxa 16,21,24,54 . In the case of the water flea D. pulex, the study specifically reported 222 that earlier inferences of deleterious mutation accumulation under asexuality were incorrect, 223 as the detected deleterious mutations in asexual strains were inherited from the sexual 224 ancestor and did not accumulate after the transition to asexuality 24 . 225 226 In summary, results from genome-wide studies addressing the prediction of deleterious 227 mutation accumulation in asexual species are equivocal. More studies are therefore needed. 228 A major constraint for studying deleterious mutation accumulation, and the reason why it 229 was not studied in most genome studies of asexuals species (Figure 1), is that it requires 230 sexual references for comparison. Such references are either unknown or not included in 231 most published genome studies of asexuals. 232

233
The same constraints likely explain why no study has thus far addressed adaptive evolution 234 in the genome of an asexual species. The question of adaptive evolution was addressed 235 indirectly in the amazon molly, by studying the amount of segregating variation at immune 236 genes (where variation is known to be beneficial). The authors found very high diversities at 237 immune genes 16 . However, these were difficult to interpret because standing variation was 238 not compared to sexual relatives, and because the amazon molly is a hybrid species. Hence 239 the high diversity could be a consequence of the hybrid origin rather than of asexuality. 240 Heterozygosity 241 Expected heterozygosity levels in asexual organisms are influenced by three major factors: 242 (1) the mechanism of transition to asexuality (which determines the initial level of 243 heterozygosity) (Box 1), (2) the cellular mechanism underlying asexuality (which determines 244 whether heterozygosity should increase or decrease over time) (Box 2), and (3) for how long 245 a species has been reproducing asexually (because effects of asexuality accumulate over 246 time). 247

248
As expected, all of the asexual species with a known hybrid origin display high 249 heterozygosity levels (1.73% -8.5%, Figure 2), while the species with an intraspecific origin 250 of asexuality show low heterozygosity levels (0.03% -0.53%, Figure 2). However, it is 251 important to note that hybrid origin is correlated with polyploidy in our dataset, and that haplotype divergence and, if high enough, can even result in a net loss of heterozygosity 262 over time, even under mitotic asexuality 8,17 . In spite of the prediction that the cellular 263 mechanism of asexuality should affect heterozygosity, the cellular mechanism of asexuality 264 appears to have little or no effect on heterozygosity levels once we control for the effect of 265 hybrid origins (Figure 2). However, we have very little power to detect such effects, 266 especially because our dataset does not include any asexual species that uses mitotic 267 asexuality but is not of hybrid origin. Nevertheless, it is interesting to note that species with 268 different forms of meiotic asexuality (including gamete duplication and central fusion) feature 269 similarly low heterozygosity levels. This suggests that although the rate of heterozygosity 270 loss is expected to vary according to mechanisms of asexuality, this variation is only relevant 271 very recently after transitions to asexuality, and no longer affects heterozygosity among 272 established asexual species. Alternatively, variation in heterozygosity caused by different 273 forms of meiotic asexuality may be too small to be picked up with our methods. Heterozygosity in polyploids is estimated via the proportion of sites that differ in at least one 302 of the homologous regions (see Box 3). This means that the estimated genome-wide 303 heterozygosity could be generated by a single haplotype that is highly divergent while others 304 are similar, or by homogeneous divergence across all copies present, or a combination of 305 these. We therefore decomposed genome-wide heterozygosity for each polyploid genome 306 into portions with different divergence structures (Figure 3). 307 308 In polyploid species of hybrid origin (sexual or asexual), heterozygosity is generally driven by 309 divergence between haplotypes originating from different species (hereafter homoeologs, 310 following the terminology of Glover et al 57 ). In our dataset, the polyploid species with 311 confirmed hybrid origins are nematodes in the genera Meloidogyne (five species) and 312 Panagrolaimus (one). As expected, heterozygosity in these species is largely dominated by 313 divergence between homoeologs. In the triploid species, divergence is between a single 314 homoeolog and two similar homologs (yellow portions in Given that in asexual polyploids of hybrid origin we expect and observe highly 320 heterogeneous divergences among haplotypes, while polyploidy of intra-specific origin is 321 predicted to generate homogeneous divergences, haplotype divergences can be used to 322 infer the origin of asexuality in polyploid species. Notably, the highly asymmetric divergence 323 levels between haplotypes in the four bdelloid rotifers (Figure 3) are best explained by a 324 hybrid origin of bdelloids. When tetraploidy was first discovered in bdelloids, it was proposed 325 that tetraploidy stemmed from either a whole genome duplication or a hybridization event in 326 their ancestor 58 . However, studies of bdelloid rotifers traditionally refer to the divergent 327 haplotypes as "ohnologs" (e.g., 17,18 ), which, following the unified vocabulary of Glover et al 57 328 would imply that the diverged haplotypes are products of a whole genome duplication. 329 However, given their likely hybrid origin, referring to them as homoeologs appears more 330 appropriate. 331 332 Whether the crayfish is also of hybrid origin remains as an open question. Its genome 333 features two nearly identical haplotypes and one that is substantially divergent (1.8%, Figure  334 3), which is suggestive of a hybrid origin. However, the divergence of the latter is within the range of heterozygosity commonly observed in sexual species, and therefore we cannot 336 clearly distinguish between an intra-specific or a hybrid origin. 337 338 Our analyses also reveal that the divergence of homologs varies extensively among bdelloid 339 rotifer genera. Divergence is very low in Rotaria (0% in R. magnacalcarata and 0.25% R. 340 macrura) and low in A. ricciae (0.5%) but relatively high in A. vaga (2.4%). The mechanisms 341 causing these differences remain unknown. In A. vaga it has been suggested that gene 342 conversion reduces divergence between homologs in some genome regions 17 . It is possible 343 that rates of gene conversion are higher in Rotaria than Adineta, for unknown reasons.  (Table 1). Not a single palindrome was detected in the remaining 13 species. The 390 frequency of palindromes had no phylogenetic signal; for example, although we found 19 391 palindromes in A. vaga, we found no palindromes in the three other bdelloid rotifers (in 392 agreement with 18 ). There is also no indication for major rearrangements being present 393 solely in very old asexuals; among the very old asexuals, the non-A. vaga rotifers along with 394 the Diploscapter nematodes have either no or only a single palindrome. 395 396 Adineta vaga and F. candida are the only two species with more than 100 genes potentially 397 affected by palindrome-mediated gene conversion, but even for these two species, the 398 overall fraction of genes in palindromes is very small (1.23% and 0.53% respectively). The 399 fraction of genes in the other seven species ranges between 0.01% and 0.16%, suggesting 400 that palindromes do not play a major role in the genome evolution of any of the asexual 401 lineages analyzed. Our findings substantiate the conclusion of a previous study 18   comparable to other asexual animal taxa (Figure 4), all of which are considerably younger 449 than the bdelloids. Across the 24 genomes, there was large variation in total TE content, 450 overall ranging from 6.6% to 17.9%, but with one species, the marbled crayfish, reaching 451 34.7%. Nevertheless, the abundance of TEs in asexual animal genomes appears to be 452 generally lower than in sexual species, which range typically from 8.5-37.6% (median: 453 24.3%) 77 . Whether this difference is indeed driven by asexuality remains an open question 454 as TE loads are known to be highly lineage-specific 20,78 . Furthermore, we annotated TEs in 455 each genome via homology searches in general databases (see methods). This can result in 456 an underestimation of TE loads relative to annotations based on species-specific TE 457 libraries. However, this is unlikely to be of major concern in our study since the methods we 458 used allowed us to identify more TEs than most of the individual genome studies 459 In addition to other lineage-specific characteristics, the cellular mechanisms underlying 464 asexuality could also affect TE loads. For example, most forms of meiotic asexuality can 465 allow for the purging of heterozygous TE-insertions, given the loss of heterozygosity 466 between generations (Box 2). Barring potential gene conversion events, this form of purging 467 cannot occur under mitotic asexuality. However, in the genomes analyzed here, we did not 468 find any effect of cellular mechanisms on TE loads (Supplementary Figure 3), likely 469 because the expected effect of the cellular mechanisms is very small relative to lineage-470 specific mechanisms. Moreover, host TE suppression mechanisms can contribute to the 471 inactivation and subsequent degeneration of TE copies over time, independently of the 472 cellular mechanism of asexuality 72,79 . 473 474 Two asexual animals clearly stand out (Figure 4), one for very low TE contents (the rotifer A. 475 ricciae; <1% of the genome) and one for very high contents (the marbled crayfish P. 476 virginalis >34%). There is currently no known mechanism that could help explain why A. Asexual animals are predicted to lose genes underlying sexual reproduction traits, including 500 male-specific traits and functions (e.g. male-specific organs, spermatogenesis), as well as 501 female traits involved in sexual reproduction (e.g., pheromone production, sperm storage 502 organs) 83 . In the absence of pleiotropic effects, gene loss is expected due to mutation 503 accumulation in the absence of purifying selection maintaining sexual traits, as well as to 504 directional selection to reduce costly sexual traits 84 . Some gene loss consistent with these 505 predictions is documented. For example, the sex determination genes xol-1 and tra-2 are 506 missing in the nematode D. coronatus 28 . Furthermore, genes believed to be involved in 507 male functions harbour an excess of deleterious mutations in the wasp Leptopilina clavipes 508 19 , which could represent the first step towards the loss of these genes. However, a similar 509 excess of deleterious mutations in genes with (presumed) male-specific functions was not 510 detected in the amazon molly P. formosa 16 . 511 512 Species reproducing via mitotic asexuality are further predicted to lose genes specific to 513 meiotic processes 85 . The genes involved in meiosis have been studied in three of eight 514 mitotic parthenogens, as well as in Rotaria rotifers and Diploscapter nematodes, whose 515 cellular mechanisms of asexuality are unknown. Most meiotic genes have been found in the 516 four bdelloid rotifers 17,18 and in both species of Diploscapter nematodes 28,29 . There was also 517 no apparent loss of meiosis genes in the amazon molly P. formosa 16 . As much as the idea 518 is appealing, there does not seem to be any support for the predicted loss of meiotic genes 519 in mitotic asexuals. We note that the lack of our understanding of meiosis on the molecular 520 level outside of few model organisms (particularly yeast and C. elegans) makes the 521 interpretation of gene loss (or absence thereof) difficult. This is best illustrated by the fact 522 that losses of meiosis genes have also been reported in different sexual species, where 523 meiosis is clearly functional 86 . 524 525 In summary, some gene loss consistent with the loss of different sexual functions has been 526 reported in several asexual species. However, a clear interpretation of gene loss in asexual 527 species is problematic because the function of the vast majority of genes is unknown in 528 these non-model organisms. 529 Horizontal gene transfer 530 Asexual species could harbour many genes acquired via horizontal gene transfer (HGT) as a 531 consequence of relaxed selection on pairing of homologous chromosomes. It has also been 532 proposed that HGTs represented an adaptive benefit which allows for the long term 533 maintenance of asexuality 87 . Indeed, bdelloid rotifers have been reported to carry an 534 unusually large amount (6.2% -9.1%) of horizontally acquired genes compared to sexual 535 lophotrochozoan genomes (0.08% -0.7%) 18,88 . Many of these have contributed to adaptive 536 divergence between bdelloid rotifer species 89 . However, there are no other ancient asexuals 537 sequenced and evaluating the role of HGTs in the long-term persistence of asexuality is 538 therefore not possible. In more recent asexuals, levels of HGT appear mostly low, e.g. in 539 Panagolaimus (0.63% -0.66%) and in two tardigrade species (0.8% -0.97%) 18,30,37 . The only 540 genome with a high reported fraction of HGT (2.8%) outside of the rotifers is the springtail F. 541 candida 23 . This is a meiotic asexual, hence a relaxed constraint on chromosome pairing did 542 not contribute to the high retention of horizontally acquired genes. Nevertheless, the 543 presence of a gene for lignocellulose degradation in the springtail and in the root-knot 544 nematode M. incognita, which was likely acquired via HGT in both species, supports an 545 adaptive role of HGT in these asexuals 23,33 . However, such isolated events of adaptive 546 HGTs are not specifically linked to asexuality, since they are reported in sexual species as 547 well 90 . The potential relation of HGT and asexuality will remain unclear until we are able to 548 reliably identify HGTs in more genomes of asexual as well as sexual species. The variation among genomes of asexual species is at least in part due to species-or 586 lineage-specific traits. But variation among the features detected in the published single-587 genome studies is also generated by differences in the methods used. Such differences are 588 often less obvious, and maybe less interesting to discuss, yet they can be critical in our 589 assessment of genome diversity among animals. In this work we thus re-analyzed several 590 key genome features with consistent methods. To minimize the effect of differences in 591 genome quality, we have used in priority robust methods, e.g. based on sequencing reads 592 rather than from assemblies. For example, re-estimating heterozygosity levels directly from 593 reads of each species allowed to show a strong effect of hybrid origin, but not of cellular 594 mechanism of asexuality (Figure 2). Another advantage of using the same methods over all 595 species is that it diminishes the "researcher degrees of freedom" 92-94 . For example, the 596 analysis of polyploid genomes requires choosing methods to call heterozygosity and ploidy. 597 By providing a common framework among species, we have shown that homoeolog 598 divergence is very diverse among polyploid asexuals. 599 600 We have identified hybrid origin as the major factor affecting heterozygosity levels across all 601 asexual animal species with available genomic data. This is consistent with the conclusions 602 of two studies that focussed on individual asexual lineages: hybridization between diverse 603 strains explains heterozygosity in Meloidogyne root knot nematodes and Lineus ribbon 604 worms 34,54 . This rule applies more generally to all the species analysed with known 605 transitions to asexuality, but it is important to highlight that all the non-hybrid species in our 606 dataset are hexapods. Thus in principle the low heterozygosity could be a hexapod specific 607 pattern, for example due to high gene conversion rates in hexapods. The taxonomic range of 608 the sequenced species is wide but we are missing several clades rich in asexual species, 609 such as mites or annelids 95,96 . These clades would be useful foci for future genomic studies 610 of asexual species. Independently of the findings of such future studies, our results suggest 611 that mitotic gene conversion (that acts independently of palindromes) plays a significant and 612 highly underappreciated role in the evolution of asexual species of intraspecific origin. For 613 example, it has been argued that one of the main benefits of sex could be the masking of 614 recessive deleterious mutations (referred to as "complementation") which would be exposed 615 under many forms of meiotic asexuality 97,98 . If gene conversion is indeed pervasive, these 616 arguments would extend to functionally mitotic forms of asexuality. Conversely, high rates of 617 gene conversion could also allow for the purging of deleterious mutations while in the 618 heterozygous state, as in highly selfing species (eg. 99,100 ). Such purging could help explain 619 why most of the genome scale studies did not find support for the theoretical expectation 620 that asexual reproduction should result in increased rates of deleterious mutation 621 accumulation (see section Mutation accumulation and positive selection). More 622 generally, given the major differences in genome evolution for asexuals of intra-specific vs. 623 hybrid origin, our study calls for future theoretical approaches on the maintenance of sex that 624 explicitly consider the loss vs. the maintenance of heterozygosity in asexuals. 625 626 In our evaluation of the general consequences of asexuality, we were not able to take two 627 key aspects into account: survivor bias of asexual lineages, and characteristics of sexual 628 ancestors. How often new asexual lineages emerge from sexual ancestors is completely 629 unknown, but it has been speculated that in some taxa asexual lineages might emerge 630 frequently, and then go extinct rapidly because of negative consequences of asexuality. In 631 other words, asexuals that would exhibit the strongest consequences of asexuality, as 632 predicted by theoretical models, are expected to go extinct the fastest. Such transient 633 asexuals remain undetected in natural populations, because research focuses on asexual 634 species or populations, and not on rare asexual females in sexual populations. Indeed, most 635 of the species included in our study have persisted as asexuals for hundreds of thousands to 636 millions of years. They might thus be mostly representative of the subset of lineages that 637 suffer weaker consequences of asexuality. Finally, the key constraint for identifying 638 consequences of asexuality is that almost none of the published genome studies of asexual 639 animals included comparisons to close sexual relatives. This prevents the detection of 640 specific effects of asexuality, controlling for the variation among sexual species -which is 641 extensive for all of the genome features we analyzed and discussed in our study. Overall, 642 despite the importance of recombination rate variation for understanding the evolution of 643 sexual animal genomes (e.g., 101,102 ), the genome-wide absence of recombination does not 644 appear to have the dramatic effects which are expected from classical theoretical models. 645 The reasons for this are probably a combination of lineage-specific patterns, differences 646 according to the origin of asexuality, and survivor bias of asexual lineages. 647

648
We combined different methods into a complete pipeline that collects published assemblies, 649 sequencing reads, and genome annotation data from online databases, and automatically 650 computes the genome features discussed here. The methods for the different steps in the 651 pipeline are detailed below. The pipeline is available at 652 https://github.com/KamilSJaron/genomic-features-of-asexual-animals. We used this pipeline 653 to gather and analyze the data for 29 sequenced individuals from 24 asexual species. For 654 some species, additional genomes to the ones we used were available, but we did not 655 include them because of low data quality and/or unavailable illumina reads (this was the 656 case for one sample of M. incognita, M. floridensis and multiple samples of D. pulex 24,33,36 ). 657 Overall, the genome features computed were: ploidy, genome size, heterozygosity, 658 heterozygosity, haplotype divergence structure, transposable elements/ repeat content, 659 conserved gene content (see Supplementary Text S3), and palindrome abundance. 660 661 Core genome features (ploidy, haploid genome size, heterozygosity, repetitive fraction of the 662 genome, and characterisation of TE content) were estimated directly from sequencing reads 663 to avoid potential assembly biases in reference genome-based approaches. The raw reads 664 were publicly available for 27 samples and for three more samples shared by authors on 665 request. We cleaned the raw reads by removing adaptors and low quality bases using 666 Skewer (parameters "-z -m pe -n -q 26 -l 21") 103 . 667

668
We used smudgeplot v0.1.3 (available at https://github.com/tbenavi1/smudgeplot) to 669 estimate ploidy levels. This method extracts from the read set unique kmer pairs that differ 670 by one SNP from each other. These kmer pairs are inferred to derive from heterozygous genome regions. The sum of coverages of the kmer pairs is then compared against their 672 coverage ratio. This comparison separates different haplotype structures (Supplementary 673 Figure 1b). The most prevalent structure is then indicative of the overall ploidy of the 674 genome. We used this ploidy estimate in all species, except A. vaga. The most prevalent 675 structure suggested that this species is diploid. A. vaga is well characterized as tetraploid 58 , 676 but we were unable to detect tetraploidy because homoeologs are too diverged to be 677 identified as such by the kmer-based smudgeplot method. 678 679 Using the inferred ploidy levels, we then estimated genome size and heterozygosity using an 680 extended version of GenomeScope 104 . GenomeScope estimates genome wide 681 heterozygosity via kmer spectra analysis, that is, by fitting a mixture model of evenly spaced 682 negative binomial distributions, where the number of fitted distributions is decided given the 683 input ploidy. Estimated distributions correspond to kmers that occur once, twice, etc., in the 684 genome. Fits are then used to estimate heterozygosity, the fraction of repeats in the 685 genome, as well as the 1n sequencing coverage. The latter is subsequently used for 686 estimation of genome size. The definition of heterozygosity for polyploids is not well 687 established (see Box 3), but GenomeScope distinguishes different types of heterozygous 688 loci in polyploids (as shown in Figure 3). 689 690 Kmer spectra analysis is affected by the choice of kmer length. Longer kmers require higher 691 sequencing coverage, but lead to more informative kmer spectra. We have chosen the 692 default kmer size 21 nt for all species except the marbled crayfish, where we chose kmer 693 length 17 nt due to low sequencing coverage. 694

695
We quantified transposable elements using DnaPipeTE v1.2 105 . The method uses haploid 696 genome size (parameter -genome_size) to subsample sequencing reads to 0.5x coverage 697 (parameter -genome_coverage). Subsampled reads are then assembled using an assembler 698 that can deal with uneven coverages, and annotated using the database of known TEs. This 699 process is repeated three times (parameter -sample_number), and the union of results 700 represents the repeat library. Additionally, repeats are annotated as TEs if their sequence 701 matches known TEs by homology (for details see 105 ). Our reported values of TE loads 702 include only repeats that were annotated as TEs, i.e., we did not include 'unknown' repeats 703 which consist of tandem repeats (satellite repeats), duplications or very divergent/unknown 704 We performed collinearity analysis using MCScanX (untagged version released 28.3.2013) 709 106 , allowing even a single gene to form a collinear bloc (parameter -s) if there were fewer 710 than 100 genes in between (parameter -m). The output was then filtered to contain only 711 blocs on the same scaffold in a reverse order. Furthermore we filtered all the homologous 712 gene pairs that have appeared on the same strand. All the remaining blocks are 713 palindromes, blocs built of reverse complementary genes on the same scaffold. See 714 Supplementary Text S2 for more details. 715