Genomic Features of Parthenogenetic Animals

Abstract Evolution without sex is predicted to impact genomes in numerous ways. Case studies of individual parthenogenetic animals have reported peculiar genomic features that were suggested to be caused by their mode of reproduction, including high heterozygosity, a high abundance of horizontally acquired genes, a low transposable element load, or the presence of palindromes. We systematically characterized these genomic features in published genomes of 26 parthenogenetic animals representing at least 18 independent transitions to asexuality. Surprisingly, not a single feature was systematically replicated across a majority of these transitions, suggesting that previously reported patterns were lineage-specific rather than illustrating the general consequences of parthenogenesis. We found that only parthenogens of hybrid origin were characterized by high heterozygosity levels. Parthenogens that were not of hybrid origin appeared to be largely homozygous, independent of the cellular mechanism underlying parthenogenesis. Overall, despite the importance of recombination rate variation for the evolution of sexual animal genomes, the genome-wide absence of recombination does not appear to have had the dramatic effects which are expected from classical theoretical models. The reasons for this are probably a combination of lineage-specific patterns, the impact of the origin of parthenogenesis, and a survivorship bias of parthenogenetic lineages.

The switch from sexual reproduction to obligate, female-producing parthenogenesis (thelytoky) has occurred repeatedly among animals and is phylogenetically widespread, with several thousand parthenogenetic animal species described (Bell 1982;van der Kooi et al. 2017;Liegeois et al. 2020). Parthenogenesis is predicted to have many consequences for genome evolution since gamete production via meiosis is heavily modified and the restoration of somatic ploidy levels via fertilization no longer takes place. Predicted consequences of parthenogenesis include the accumulation of deleterious mutations (Muller 1964;Felsenstein 1974;Keightley and Otto 2006), as well as changes in intragenomic heterozygosity levels (Birky 1996;Balloux et al. 2003) and transposable element (TE) dynamics (Hickey 1982). In the present study, we evaluate whether parthenogenesis indeed generates these predicted genomic signatures by reanalyzing and comparing the published genomes of 26 parthenogenetic animal species ( Figure 1). Previous genome studies were unable to address the question of how parthenogenesis affects genome evolution because they focused on individual lineages. Yet because parthenogenesis is a lineage-level trait, disentangling causes of parthenogenesis from lineage-level characteristics require replication across independently evolved instances of parthenogenesis. Our study includes species from at least 18 independently evolved parthenogenetic lineages from 5 different animal phyla, providing us with the unique opportunity to detect general consequences of parthenogenesis that are not solely consequences of lineage-specific evolution. Furthermore, we study the same features in all genomes, whereas many of the original studies focused on different genomic features (Figure 1), which thus far precluded broad comparisons across different parthenogenetic groups.
Because the predicted consequences of parthenogenesis are strongly affected by how asexuality evolved from the sexual ancestor (Box 1) as well as by the cellular mechanisms underlying parthenogenesis (Box 2), we include biological differences among species in our comparisons. For example, some parthenogenetic species have evolved via hybridization (Box 1), which generates incipient parthenogens with high intragenomic heterozygosity and can result in increased activity of TEs (Arkhipova and Rodriguez 2013;Neiman et al. 2014;Rodriguez and Arkhipova 2018). In such instances, it can be difficult to disentangle the consequences of hybridization from those of parthenogenesis. Similarly, some cellular mechanisms underlying parthenogenesis involve meiotic divisions, with a secondary restoration of somatic ploidy levels, while others do not. In the former case, heterozygosity is expected to decay rapidly, while in the latter case, it could be maintained or even increase over time (Engelstädter 2017). Finally, because the genome studies differed in their focus and in the methods used, we reanalyzed the published genomes with standardized approaches. Whenever possible, we conducted quantitative comparisons between groups of parthenogenetic species. However, for interpretation, it is important to consider that the available genomes are neither a random nor a fully representative sample of parthenogenetic animals.
Using the 26 genomes, we studied 9 genomic features that have been proposed to be affected by parthenogenesis. Four of them represent classical theoretical predictions for consequences of asexuality on genome evolution, namely, consequences for mutation accumulation, positive selection, TE dynamics, and intragenomic heterozygosity (see below). The 5 remaining ones are unusual genomic features that were observed in individual parthenogenetic species and suggested to be linked to their mode of reproduction: horizontally acquired genes, palindromes, gene conversion, gene family expansions, and gene losses, see below.
We first reviewed the literature for information on all 9 genomic features in the 26 parthenogenetic species (Figure 1). However, the methods used to evaluate a given genomic feature varied among studies of parthenogenetic species (Supplementary Figure 3), which precluded quantitative comparisons. Hence we developed a standardized pipeline to quantify heterozygosity, TE load, the frequency of horizontally acquired genes, and palindromes. Heterozygosity and TE loads were quantified using raw sequencing reads to avoid biases  Figure 1. Genome features studied in parthenogenetic animal species. The phylogeny displays the taxonomic relationships of the 26 sequenced parthenogenetic animal species considered here, representing at least 18 independent transitions to parthenogenesis from 5 different animal phyla. Species that might derive from the same original transition are grouped in triangles. The color of the circle indicates the cellular mechanism of parthenogenesis and the number inside the circle the ploidy of the species (see Supplementary Table 1 for details). We note M. floridensis as triploid, as shown by our analyses, even though it is reported as diploid in the original paper; see Supplementary Material S1 for details (Ranallo-Benavidez et al. 2020). Each original genome paper explored a given set of genome features: the green cells represent cases where the genomic feature was quantified (values are indicated); the gray cells represent studies where the genomic features were addressed with respect to parthenogenesis, but the results were not quantitatively comparable to other studies. We reanalyzed heterozygosity, palindromes, TEs, and horizontal gene transfer in this study; the discussion of the remaining features is based on the analyses reported in the individual genome studies (Abad et al. 2008;Xu et al. 2011;Flot et al. 2013;Tucker et al. 2013;Oxley et al. 2014;Kraaijeveld et al. 2016;Bast et al. 2016;Szitenberg et al. 2016;Koutsovoulos et al. 2016;Hashimoto et al. 2016;Arkhipova et al. 2017;Faddeeva-Vakhrusheva et al. 2017;Jiang et al. 2017;Hiraki et al. 2017;Fradin et al. 2017;Blanc-Mathieu et al. 2017;Szitenberg et al. 2017;Warren et al. 2018;Nowell et al. 2018;Lindsey et al. 2018;Gutekunst et al. 2018;Schiffer et al. 2018;Smith et al. 2019;Schiffer et al. 2019). Findings for mutation accumulation and adaptive evolution refer to comparisons between sexual and asexual species and are reported with respect to theoretical predictions (yes: as predicted, no: opposite to predictions, inconclusive: no difference). e/(g*nt): event per generation per nucleotide; l/s: number of lost genes among the studied genes related to sexual reproduction.

Box 1: Transitions to Parthenogenesis
Meiotic sex and recombination evolved once in the common ancestor of eukaryotes (Cavalier-Smith 2002). Parthenogenetic animals therefore derive from a sexual ancestor, but how transitions from sexual to parthenogenetic reproduction occur can vary and have different expected consequences for the genome (Neiman et al. 2014).

Hybrid Origin:
Hybridization between sexual species can generate hybrid females that reproduce parthenogenetically (Schultz 1973;Neiman et al. 2014). Parthenogenesis caused by hybridization can generate a highly heterozygous genome, depending on the divergence between the parental sexual species prior to hybridization. Hybridization can also result in a burst of TE activity (Arkhipova and Rodriguez 2013).

Intraspecific Origins:
Endosymbiont Infection. Infection with intracellular endosymbionts (such as Wolbachia, Cardinium, or Rickettsia) can cause parthenogenesis, a pattern that is frequent in species with haplodiploid sex determination (van der Kooi et al. 2017). This type of transition often (but not always) results in fully homozygous lineages because the induction of parthenogenesis frequently occurs via gamete duplication (Box 2).

Spontaneous Mutations/Contagious Parthenogenesis.
Spontaneous mutations can also underlie transitions from sexual to parthenogenetic reproduction. In addition, parthenogenetic females of some species produce males that mate with females of sexual lineages, and thereby generate new parthenogenetic strains (contagious parthenogenesis). In both cases, the genomes of incipient parthenogenetic lineages are expected to be very similar to those of their sexual relatives and subsequent changes should be largely driven by the cellular mechanism underlying parthenogenesis (Box 2).

Box 2: Cellular Mechanisms of Parthenogenesis
In sexual species, offspring are generated through the fusion of male and female gametes. In parthenogens, females generate diploid (or polyploid) offspring from unfertilized oocytes via different cellular mechanisms. The mechanism is predicted to affect genome evolution and especially heterozygosity levels. For details see (Suomalainen et al. 1987;Engelstädter 2017).
Functionally Mitotic Parthenogenesis (Apomixis). Under functionally mitotic asexuality, no ploidy reduction occurs and offspring are clones of their mother. Note that because gametogenesis and egg production are tightly linked to meiotic divisions in sexual species, functionally mitotic parthenogenesis derives from meiosis and does therefore generally not correspond to mitotic divisions in terms of molecular mechanisms. Many functionally mitotic parthenogens notably still feature vestiges of the ancestral meiotic divisions such as homologous chromosome pairing or chromosome condensation typical for meiosis.
Meiotic Parthenogenesis (Automixis). Under meiotic parthenogenesis, meiotic divisions occur partially or completely, but somatic ploidy levels are maintained via different mechanisms. Depending on how recombination occurs, some of these mechanisms are equivalent to functionally mitotic parthenogenesis, even though meiosis is fully maintained.
Endoduplication. A duplication of the entire chromosome set occurs before normal meiosis, during which ploidy is reduced again. If recombination occurs between identical chromosome copies rather than between chromosome homologs, endoduplication is functionally mitotic, that is, produces offspring that are clones of their mother.
Central Fusion and Terminal Fusion. Under these 2 mechanisms, somatic ploidy levels are restored through the introduced by differences in genome assembly quality, while analyses of horizontally acquired genes and palindromes were based on the published genome assemblies because of methodological requirements. We relied on published information (without reanalysis) for evaluating mutation accumulation, positive selection, gene family expansions, gene losses, and gene conversion because these genome features cannot be studied using individual genomes of parthenogenetic species. Indeed, analyzing these features requires genomes of sexual relatives and/or population genomic data which is not available for the vast majority of the 26 parthenogenetic species.

Overview of Species and Genomes Studied
We reanalyzed the published genomes of 26 asexual animal species with the aim of identifying general genomic signatures of asexualilty. The 26 species correspond to at least 18 independent transitions to asexuality and cover a broad taxonomic range, including chordates, rotifers, arthropods, nematodes, and tardigrades. In addition to covering this taxonomic range, these asexual species vary in the cellular mechanisms underlying parthenogenesis, in the mechanisms that caused the transition to asexuality, ploidy, as well as in other biological aspects (Figure 1, Supplementary  Tables 1 and 2). This variation allows us to assess whether parthenogenesis generates universal genomic signatures independently of species-specific traits.
The cellular mechanisms underlying parthenogenesis have been reported in 22 of the 26 species. Eleven of them reproduce via functionally mitotic parthenogenesis, while the 11 remaining species have different types of meiotic parthenogenesis (Figure 1). Information on how asexuality evolved is available for 16 of the 26 sequenced species (Supplementary Table 1). A hybrid origin has been suggested for 10 of these, based on the identification of parental species and/ or genomic analyses. Endosymbionts are the most likely cause of asexuality in 4 species (the springtail, both wasps, and the thrips), and spontaneous mutation in 2 (the ant and the cape honey bee).
Most if not all predicted consequences of asexuality are expected to accumulate over time, meaning that their effect size, as well as the power to detect them, is increased in old asexual lineages. However, estimating the age of asexual lineages is difficult and always associated with large uncertainties (Schurko et al. 2009;Neiman et al. 2009). Therefore we did not include quantitative comparisons among asexuals with respect to their age. However, because our set of species comprises some asexuals believed to be "ancient" (i.e., several million years old, see Supplementary Table 1), we discuss, where appropriate, potential age effects in a qualitative manner.

Recomputed Genomic Features
We combined different methods into a complete pipeline that collects published assemblies, sequencing reads, and genome annotation data from online databases, and automatically computes the focal genome features. The pipeline is available at https://github. com/KamilSJaron/genomic-features-of-parthenogenetic-animals. We provide a summary of each method below, with technical details (e.g., specific parameters used for each method) indicated in the Supplementary Methods. For some species, additional genomes to the ones we used were available, but we did not include them because of low data quality and/or unavailable Illumina reads (this was the case for one sample of Meloidogyne incognita, Meloidogyne floridensis, and multiple samples of Daphnia pulex (Abad et al. 2008;Tucker et al. 2013;Lunt et al. 2014

Heterozygosity
To compare intragenomic heterozygosity among species with different ploidy levels, we estimated heterozygosity as the proportion of sites with more than one allele present among all homologous genome regions (consistent with Lokki 1976). In diploid species, this genomewide heterozygosity can correspond to the divergence between alleles (homolog heterozygosity), or, if the species has a history of hybridization, to the divergence of gene copies derived from different species (hereafter homoeologs, following the terminology of Glover et al. 2016). In polyploid species, heterozygosity can be a combination of homolog and homoeolog divergence. We distinguished homolog and homoeolog heterozygosity whenever possible, or inferred a "composite heterozygosity" (the sum of the 2) otherwise. To avoid biases fusion of 2 of the 4 meiotic products (products separated during the first meiotic division merge under central fusion, products separated during the second division merge under terminal fusion). In the absence of recombination, central fusion generates offspring that are clones of their mother, while terminal fusion generates fully homozygous offspring. The consequences for heterozygosity are opposite under inverted meiosis, where chromatids are separated during meiosis I and chromosomes during meiosis II. For example, terminal fusion with an inverted sequence of meiosis and no recombination (not shown here) generates offspring that are clones of their mother (see Lenormand et al. 2016 for a recent review).
Gamete Duplication. After a full meiosis, a haploid meiotic product undergoes duplication. This results in a diploid, but fully homozygous offspring. stemming from variable genome assembly qualities, we estimated heterozygosity directly from sequencing reads using kmer spectra analysis (Ranallo-Benavidez et al. 2020), except for 4 bdelloid rotifers where the heterozygosity levels between homoeologs exceeded the range quantifiable by this method (Supplementary Material S2).
Kmer spectra analysis is based on the decomposition of sequencing reads into all possible subsequences (kmers) and generating a coverage histogram of the kmers (kmer spectrum). We performed the kmer spectra analysis using GenomeScope 2.0 (Ranallo-Benavidez et al. 2020). Genomescope fits a mixture model of evenly spaced negative binomial distributions to the kmer spectrum. Fits are then used to estimate heterozygosity. With the exception of bdelloid rotifers, we are not able to directly compare the divergence of individual haplotypes (because such a comparison requires phased genomes). However, we are able to measure the haplotype structure on a per-locus basis. These approaches notably allow us to distinguish biallelic from triallelic loci in triploid organisms.

Transposable Elements
We quantified TEs using DnaPipeTE (Goubert et al. 2015). This method uses the haploid genome size to subsample sequencing reads to a low coverage of 0.5×. These subsampled reads are then assembled using an assembler (Trinity) that can deal with uneven coverages and is able to split assembled regions with few differences (including different TE families). The assembled sequences largely correspond to repetitions as non-repetitive genome regions present in the subsampled reads drop out at this stage because the coverage of such regions is too low for assembling. The assembled sequences are annotated by homology using a database of known TEs. The output of the method is the number of sampled nucleotides assembled and annotated as different types of repeats, and fractions are calculated as the numbers divided by the total number of sampled nucleotides. Our reported values of TE loads include only repeats that were annotated as TEs, that is, we did not include "unknown" repeats which consist of tandem repeats (satellite repeats), duplications, or very divergent/unknown TEs. Note that in addition to the 26 parthenogenetic genomes (Figure 1), we also analyzed one sexual species Procambarus fallax (SAMN06115719) for comparison to its asexual sister species Procambarus virginalis.

Palindromes
Palindromes are formed of 2 homologous reverse complementary sequences on the same chromosome (Supplementary Figure 2). Palindromes can facilitate gene conversion and therefore help to escape mutational meltdown via Muller's ratchet (Marais et al. 2010;Trombetta and Cruciani 2017). To test if they play such a role in asexual organisms we identified palindromes using colinearity analysis implemented in the program MCScanX . Detected collinear blocks were filtered to contain only reverse complementary collinear blocks on the same chromosome, since only such structures have the capacity to form a hairpin (Supplementary Figure 2). This method for palindrome identification depends on genome assemblies. Palindromes are less likely to be detected in highly fragmented assemblies, and artificial palindromes can be generated by erroneous scaffolding (see also Nowell et al. 2018). Our analyses assume that there are no systematic scaffolding errors in the published assemblies, meaning that our list of palindromes includes false positives that are generated by misassemblies in the published reference genomes. Palindrome identification methods rely on genome annotations, which are available for 23 of the 26 asexual species (all except D. pulex, Apis mellifera capensis, and Aptinothrips rufus). We screened these 23 genomes for the presence of palindromic arrangements (see Supplementary Methods for details).

Horizontal Gene Transfer
We systematically estimated the percentage of non-metazoan HGT candidates (HGT C ) in the 23 of the 26 asexual species with available gene annotations using a sequence comparison based approach, following (Nowell et al. 2018). For each species, we compared the set of annotated genes to the UniRef (Lespinet et al. 2002) and UniProtKB/Swiss-Prot protein databases to identify genes of likely non-metazoan origin (Suzek et al. 2015). We considered non-metazoan genes as HGT candidates only if they were on a scaffold that also encoded at least one gene of unambiguous metazoan origin, to control for potential contamination in the genome assemblies (see Supplementary Methods for details).

Results and Discussion
There has been an accumulation of genome studies of individual parthenogenetic animal species ( Figure 1) and most of these studies have suggested that one or a few specific genome features uncovered in their focal species might be linked to parthenogenesis. To test whether any of these genomic features were in fact general consequences of parthenogenetic reproduction rather than lineagespecific traits, we evaluated evidence from published genomes of 26 parthenogenetic animals in the light of theoretical expectations for genome evolution under parthenogenesis. We quantified heterozygosity, TE loads, the frequency of horizontally acquired genes and palindromic sequence arrangements using standardized methods, and combined these analyses with a review of published information on mutation accumulation, positive selection, gene family expansions, and gene losses.

Mutation Accumulation and Positive Selection
One of the classical predictions linked to parthenogenesis is that it reduces the efficacy of selection (Fisher 1930;Muller 1932;Muller 1964;Hill and Robertson 1966;Felsenstein 1974;Keightley and Otto 2006). This reduction occurs because linkage among loci in asexual species prevents selection from acting individually on each locus, resulting in different forms of selective interference (Otto 2020). This selective interference can result in a faster accumulation of deleterious mutations and a slower rate of adaptation. While there is accumulating evidence for these processes in experimental evolution studies (e.g., Kaltz and Bell 2002;Goddard et al. 2005;McDonald et al. 2016), their impact for natural populations remains unclear (Neiman et al. 2018). Genomic data can provide the basis for studying mutation accumulation and adaptation (positive selection) in natural populations.
In summary, results from genome-wide studies addressing the prediction of deleterious mutation accumulation in natural populations of parthenogenetic species are equivocal. More studies are therefore needed. A major constraint for studying deleterious mutation accumulation, and the reason why it was not studied in most genome studies of parthenogenetic species (Figure 1), is that it requires homologous gene sets from sexual outgroups for comparison. These species are either unknown or not included in most published genome studies of parthenogens.
The same constraints likely explain why no study has thus far directly addressed adaptive evolution in the genome of a parthenogenetic species. The question of adaptive evolution was addressed indirectly in the Amazon molly, by studying the amount of segregating variation at immune genes (where variation is known to be beneficial). The authors found very high diversities at immune genes (Warren et al. 2018). However, these were difficult to interpret because standing variation was not compared to that in sexual relatives and because the Amazon molly is a hybrid species. Hence the high diversity could be a consequence of the hybrid origin rather than of parthenogenesis. Furthermore, sexual and parthenogenetic sister species differ in other aspects besides their reproductive mode, which complicates interpretations of such individual comparisons.

Heterozygosity
Intragenomic (individual-level) heterozygosity is the nucleotidic divergence between the haploid genome copies of an individual. In a panmictic sexual population, the level of intragenomic heterozygosity corresponds to the level of genetic diversity in a population (the amount of variation observed between DNA sequences from different individuals). This is however not the case in asexual populations, which are, by definition, not panmictic (Balloux et al. 2003).
Intragenomic heterozygosity in asexual organisms is expected to depend on 3 major factors: (1) the mechanism of transition to parthenogenesis (which determines the initial level of heterozygosity, expected to be high for hybrid origins; Box 1), (2) the cellular mechanism underlying parthenogenesis (which determines whether heterozygosity will increase or decrease over time; Box 2), and (3) how long a species has been reproducing asexually (because the effect of parthenogenesis accumulates over time).
As expected, all of the species with a known hybrid origin of parthenogenesis display high heterozygosity levels (1.73-8.5%, Figure 2). By contrast, species with an intraspecific origin of parthenogenesis show low heterozygosity levels (0.03-0.83%, Figure 2). In addition to hybrid origins, polyploidy may also contribute to high heterozygosity in parthenogens as heterozygosity is higher in polyploid (1.84-33.21%) than diploid species (0.03-5.26%). It is impossible to disentangle the effects of hybrid origin from polyploidy on heterozygosity in our dataset as across the 26 species, hybrid origins are correlated with polyploidy. Six of the 11 polyploids in our sample are of hybrid origin, while for the 5 others a hybrid origin is supported by our results (see below), even though it was not suggested previously. While the correlation between hybrid origin and polyploidy in our dataset is striking, it is important to note that this correlation would most likely be weaker in a random sample of parthenogenetic animals. Indeed, many polyploid parthenogenetic animals are not of hybrid origin, including several well-studied species such as the New Zealand mudsnail Potamopyrgus antipodarum, the bush cricket Saga pedo, or the bagworm moth Dahlica triquetrella. None of these has a published genome yet, which precludes their inclusion in our study.
The heterozygosity levels present at the inception of parthenogenesis should decay over time for most forms of meiotic parthenogenesis (Schön et al. 2009;Engelstädter 2017) (see also Box 2). In functionally mitotic parthenogens, heterozygosity is expected to increase over time as haplotypes can accumulate mutations independently of each other (generating the so-called "Meselson effect") (Birky 1996). However, gene conversion can strongly reduce haplotype divergence and, if high enough, can even result in a net loss of heterozygosity over time, even under functionally mitotic parthenogenesis (Birky 1996;Flot et al. 2013). In spite of the prediction that the cellular mechanism of parthenogenesis should affect heterozygosity, it appears to have no detectable effect on heterozygosity levels once we control for the effect of hybrid origins (Figure 2). Indeed, heterozygosity levels in the 11 functionally mitotic parthenogens are high, but all these species are of hybrid origin. Furthermore, 9 of the 11 species are polyploid (the diploid species are the Amazon molly and the nematode Diploscapter pachys). Conversely, all the species with meiotic parthenogenesis are diploid. This is expected given that polyploidy can generate problems during meiosis (reviewed in Cifuentes et al. 2010), but complicates the interpretation of heterozygosity levels among species with different cellular mechanisms of parthenogenesis. Nevertheless, for the species that are not of hybrid origin, it is interesting to note that different forms of meiotic parthenogenesis (including gamete duplication, terminal fusion, and central fusion) are associated with similarly low heterozygosity levels. This suggests that although the rate of heterozygosity loss is expected to vary according to mechanisms of parthenogenesis, this variation is only relevant very recently after transitions to parthenogenesis, and no longer affects heterozygosity among established parthenogenetic species. Consistent with this view, the 2 non-hybrid parthenogens with relatively higher heterozygosity are the youngest ones in our dataset: 40 years for the cape bee (Smith et al. 2019) and 100 years for the raider ant (Oxley et al. 2014). Alternatively, variation in heterozygosity caused by different forms of meiotic parthenogenesis may be too small to be detected with our methods.

Heterozygosity Structure in Polyploids
In polyploids, the estimated genome-wide heterozygosity can be generated by a single-genome copy that is highly divergent while 33.2 28.9 27.9 4.9 4.3 1.8 5 2.6 5.2 5.7 2.4 5.5

P . v ir g in a li s P . d a v id i M . fl o ri d e n s is M . e n te ro lo b ii M . in c o g n it a M . a re n a ri a M . ja v a n ic a A . v a g a A . ri c c ia e R . m a c ru ra R . m a g n a c a lc a ra ta
triploids tetraploids others are similar, or by homogeneous divergence across all copies present, or a combination of these. We, therefore, characterized the most likely origin of heterozygosity for the polyploid species in our dataset. The heterozygosity of 2 of the 5 triploid species in our dataset (the crayfish P. virginalis and nematode M. floridensis) is generated mostly by biallelic loci (Figure 3, Supplementary Materials S3). The very low proportion of triallelic loci in these 2 genomes suggests an AAB structure of the 2 genomes, where 2 of the haploid genome copies (A) are nearly identical and the last copy (B) is the carrier of the observed heterozygosity. This AAB model is in agreement with the previous genomic analysis of the parthenogenetic crayfish P. virginalis (Gutekunst et al. 2018). This previous analysis suggested that P. virginalis has a non-hybrid origin, which would imply that the divergence between the A and B haploid genomes corresponds to the heterozygosity present in the sexual ancestor P. fallax from which P. virginalis derived only 30 years ago. However, our estimation of the divergence of the third genome copy B (1.8%) exceeds by far our estimation of the heterozygosity in P. fallax (0.76%; Supplementary Materials S4). We, therefore, suggest a hybrid origin of P. virginalis, similar to the well established hybrid origins of the other polyploid asexuals in our dataset. The second triploid species in our dataset with an AAB genome structure is the root-knot nematode M. floridensis, which was previously mistaken for diploid (Supplementary Materials S1). This genome structure distinguishes M. floridensis from the other triploid Meloidogyne species (whose genomes are comprised of a larger portion with an ABC structure, Figure 3), which might suggest that the origin of triploidy in M. floridensis is independent of the origin of triploidy in the other Meloidogyne species.
In the tetraploid species, the biallelic loci can be sorted into those where one genome copy carries the alternative allele (yellow portions in Figure 3), and those where 2 genomic copies carry the alternative allele (pink portions). The genomes of the 2 tetraploid Meloidogyne species contain high proportions of all heterozygosity structures ( Figure 3) suggesting a complex genomic structure, such as AABC or ABCD. Alternatively, this signal could be also caused by partial aneuploidy, which is common among Meloidogyne species (Triantaphyllou 1981).
Haplotype divergence can provide information on the origin of parthenogenesis in polyploid species: in asexual polyploids of hybrid origin we expect, and indeed observe, highly heterogeneous divergences among haplotypes, while in those of intraspecific origin we expect homogeneous divergences. Notably, divergence between the haplotypes of the 4 bdelloid rotifers is highly asymmetric (Figure 3). When tetraploidy was first discovered in bdelloids, it was proposed that it stemmed from either a whole-genome duplication or a hybridization event in their ancestor (Mark Welch et al. 2008). Studies of bdelloid rotifers traditionally refer to the divergent haplotypes as "ohnologs" (e.g., Flot et al. 2013;Nowell et al. 2018), which, following the unified vocabulary of Glover et al. (Glover et al. 2016) would imply that the diverged haplotypes are products of a wholegenome duplication. However, the most parsimonious explanation for the highly asymmetric divergence of the different bdelloid haplotypes is a hybrid origin. Referring to the diverged haplotypes as homoeologs, therefore, appears more appropriate. Our analyses also indicate that the allelic heterozygosity varies extensively among bdelloid rotifer genera. Divergence is very low in Rotaria (0.49% in Rotaria magnacalcarata and 0.125% Rotaria macrura) but relatively high in Adineta vaga (2.4%) and in Adineta ricciae (5.5%). There is currently no good explanation for the higher allelic heterozygosity in A. ricciae compared to A. vaga, but our analyses of the A. ricciae sequencing reads suggest that different ploidy levels in the 2 species could play a role (Supplementary Materials S2). Additionally, variable gene conversion and mutation rates could also contribute to the different allelic heterozygosities. In A. vaga, it has been suggested that gene conversion reduces divergence between homologs in some genome regions (Flot et al. 2013) but there are currently no direct estimates for gene conversion rates in rotifers. Independently of the mechanisms causing the differences between bdelloids, it is important to note that with such low levels of divergence between homologs, there can be no strong genome-wide "Meselson effect" in bdelloid rotifers (see also Flot et al. 2013). It remains possible that the subset of genomic regions with divergence between homologs in Adineta features allele phylogenies as expected under the "Meselson effect." This is the case in the asexual unicellular eukaryote Trypanosoma brucei gambiense: some genome regions feature high heterozygosity and allele phylogenies as expected under the "Meselson effect," while others are largely homozygous (Weir et al. 2016). Again, it remains unknown why there is such extensive heterogeneity in divergence across the genome in this species. A possible explanation is that the heterozygous genome regions are the consequence of ancient introgression, and that gene conversion rates are low in such regions because of their very high heterozygosity (see Conclusions).

Palindromes and Gene Conversion
Palindromes are duplicated regions on a single chromosome in reverse orientation. Because of their orientation, palindromes can align and form hairpins, which allows for gene conversion within duplicated regions (Supplementary Figure 2). Palindrome-mediated gene conversion was shown to play a major role in limiting the accumulation of deleterious mutations for non-recombining human and chimpanzee Y chromosomes (Rozen et al. 2003;Marais et al. 2010;Trombetta and Cruciani 2017). Indeed, approximately onethird of coding genes on these Y chromosomes occur in palindromes, and the highly concerted evolution of palindromic regions indicates that the rates of gene conversion are at least 2 orders of magnitude Only species with at least one palindrome detected are listed in the table. Rows in bold highlight species with more than 100 genes detected in palindromes. a The detected number of palindromes in these species exceeds the number reported in the corresponding genome articles (17 in A. vaga and 11 in F. candida). This is because we included individual genes in palindromic arrangements, whereas the original genome studies only included genes if they were in palindromic synteny blocks of at least 5 genes. See also Supplementary Material S5. higher in the palindromes than between homologous chromosomes (Rozen et al. 2003). The reports of palindromes in the genomes of the bdelloid rotifer A. vaga (Flot et al. 2013) and of the springtail Folsomia candida (Faddeeva-Vakhrusheva et al. 2017) led to the hypothesis that palindromes could play a similar role in asexual organisms by reducing deleterious mutation accumulation in the absence of recombination. However, the potential benefit of palindromemediated gene conversion depends on the portion of genes in palindromic regions (Marais et al. 2010). In addition to identifying palindromes, it is therefore important to also quantify the number of genes affected by palindrome-mediated gene conversion.
We identified 19 palindromes in A. vaga, 16 in F. candida, and 1 to 4 palindromes in 8 additional genomes (Table 1). Not a single palindrome was detected in the remaining 13 species where the available data allowed us to identify palindromes (see Methods and Supplementary Materials S2 for details). The frequency of palindromes had no phylogenetic signal; for example, although we found 19 palindromes in A. vaga, we found no palindromes in the 3 other bdelloid rotifers (in agreement with Nowell et al. 2018). However, it is important to note our analyses assume that there are no systematic scaffolding errors in the published assemblies, meaning that our list of palindromes includes false positives that are generated by mis-assemblies in the published reference genomes. There is also no indication for major rearrangements being present solely in very old asexuals; among the very old asexuals, the non-A. vaga rotifers along with the Diploscapter nematodes have either zero or only a single palindrome.
Adineta vaga and F. candida are the only 2 species with more than 100 genes potentially affected by palindrome-mediated gene conversion, but even for these 2 species, the overall fraction of genes in palindromes is very small (1.23% and 0.53% respectively). The fraction of genes in the other 7 species ranges between 0.01% and 0.16%, suggesting that palindromes do not play a major role in the genome evolution of any of the parthenogenetic lineages analyzed. Our findings substantiate the conclusion of a previous study (Nowell et al. 2018) that major genomic rearrangements and the breaking of gene syntenies do not occur at high rates in asexual organisms. They appear to occur at rates similar to those known in recombining genome portions of sexual species (Fan and Meyer 2014;Chen et al. 2017).
Gene conversion can also occur outside of palindromic regions between alleles, for example when double-stranded DNA breaks are repaired using the homologous chromosome as a template (Lee et al. 2009;Hum and Jinks-Robertson 2017). This can, in theory, contribute to the loss of heterozygosity under all forms of parthenogenesis, but allelic gene conversion rates have only rarely been studied in parthenogenetic species-or sexual ones for that matter. Allelic gene conversion rates are estimated differently in different studies and are therefore difficult to compare: in the water flea D. pulex, they were estimated to amount to approximately 10 −6 locus −1 generation −1 (Xu et al. 2011;Tucker et al. 2013;Keith et al. 2016), and in the Amazon molly Poecilia formosa to 10 −8 (Warren et al. 2018). Up to 11% of the genome of the nematode D. pachys (Fradin et al. 2017) is suggested to be homozygous as a consequence of gene conversion, and studies have also argued for an important role of gene conversion for genome evolution in root-knot nematodes (Szitenberg et al. 2017) and rotifers (Flot et al. 2013;Nowell et al. 2018), although no quantitative estimates are available for these groups.

Transposable Elements
TEs are DNA sequences that can autonomously change positions in a genome via various "cut-and-paste" and "copy-and-paste" mechanisms (Burt and Trivers 2006;Wicker et al. 2007). TEs can invade genomes even though they generally provide no adaptive Figure 4. Percentage of nucleotides annotated as TEs in parthenogenetic genomes. Both the TE load and frequency of TE classes vary substantially between individual parthenogenetic lineages. The TE classes are: class I "cut-and-paste" DNA transposons (DNA), and class II "copy-and-paste" long interspersed nuclear elements or autonomous non-LTR elements (LINEs), short interspersed nuclear elements or non-autonomous non-LTR elements (SINEs), long terminal repeat elements (LTR), and rolling-circle elements (Helitron).
advantage to the individual carrying them (Doolittle and Sapienza 1980;Le Rouzic et al. 2007; Hua- Van et al. 2011). To the contrary, new TE insertions in coding or regulatory sequences disrupt gene functions and cause deleterious effects in the host; only very rarely can specific insertions be co-opted to acquire novel, adaptive functions for the host (Hua- Van et al. 2011). In sexual organisms, TEs can spread through panmictic populations because of their ability to rapidly colonize new genomes (Hickey 1982;Zeyl et al. 1996). At the same time, sexual reproduction facilitates the purging of deleterious TE insertions, because recombination, segregation, and genetic exchange among individuals improve the efficacy of selection (Wright and Finnegan 2001;Nuzhdin and Petrov 2003). In the absence of sex, TEs could therefore accumulate indefinitely, which led to the prediction that TEs could frequently drive the extinction of parthenogenetic lineages. Only parthenogenetic lineages without active TEs, or with efficient TE suppression mechanisms, would be able to persist over evolutionary timescales (Wright and Finnegan 2001;Dolgin and Charlesworth 2006). Consistent with this view, initial investigations in bdelloid rotifers reported extremely low TE loads (Arkhipova and Meselson 2000). This prompted the authors to suggest that bdelloid rotifers could have been able to persist in the absence of canonical sex for over 40 million years thanks to their largely TE-free genomes.
Our analysis of parthenogenetic animal genomes does not support the view that bdelloid rotifers have unusually low TE contents, even though our methods tend to underestimate TE content in bdelloids relative to other parthenogenetic species (see Supplementary Methods). The estimated TE content of bdelloid rotifers (0.7% to 9.1%) is comparable to other parthenogenetic animal taxa (median 6.9%, Figure 4), all of which are considerably younger than the bdelloids. Across the 26 genomes, there was a large variation in total TE content, overall ranging from 0.7% to 17.9%, with one species, the marbled crayfish, reaching 34.7%. The extreme repeat content in the latter is largely inherited from its sexual ancestor, as the TE loads we estimated for the close sexual relative P. fallax are nearly identical (32.5%). Nevertheless, the abundance of TEs in parthenogenetic animal genomes appears to be generally lower than in sexual species, which range typically from 8.5% to 37.6% (median: 24.3%) (Canapa et al. 2015). These low abundances could stem from the evolution of reduced TE activity in parthenogenetic species (Charlesworth and Langley 1986;Bast et al. 2019), and/or appear if successful parthenogenetic taxa generally derive from sexual ancestors with largely inactive TEs and low TE contents. However, whether the apparently lower TE content in the 26 genomes is indeed linked to parthenogenesis remains an open question as TE loads are known to be highly lineage-specific (Kofler et al. 2012;Ågren et al. 2015;Bast et al. 2016;Szitenberg et al. 2016).
In addition to other lineage-specific characteristics, the cellular mechanisms underlying parthenogenesis could also affect TE loads. For example, many forms of meiotic parthenogenesis can allow for the purging of heterozygous TE insertions, given the loss of heterozygosity between generations (Box 2). However, in the genomes analyzed here, we did not find any effect of cellular mechanisms on TE loads (Supplementary Figure 8), likely because the expected effect of the cellular mechanisms is small relative to lineage-specific mechanisms. Moreover, host TE suppression mechanisms can contribute to the inactivation and subsequent degeneration of TE copies over time, independently of the cellular mechanism of parthenogenesis (Aravin et al. 2007;Hua-Van et al. 2011). Similarly, we did not find any difference in TE content according to hybrid versus intraspecific origin of asexuals (Supplementary Figure 8), even though mismatches between species-specific TEs and silencing machineries can result in increased TE activity in hybrids (Fontdevila 2005;Ågren and Wright 2011;Arkhipova and Rodriguez 2013;Romero-Soriano et al. 2017;Rodriguez and Arkhipova 2018).

Horizontal Gene Transfer
Parthenogenetic species could harbor many genes acquired via horizontal gene transfer (HGT) as a consequence of relaxed selection on pairing of homologous chromosomes (see also Gene Family Expansions). It has also been proposed that HGTs might represent an adaptive benefit that allows for the long-term maintenance of parthenogenesis (Danchin et al. 2010).
We found no evidence for parthenogenesis favoring the retention of HGTs. The majority of the species for which the available data allowed for HGT inferences (16 out of 23) showed a low proportion of candidate HGT genes, with ~1% candidate HGT (Supplementary Table 4). In agreement with previous findings, we identified elevated levels of candidate HGT in the 4 bdelloid rotifer species A. ricciae (10%), A. vaga (10.6%), R. macrura (8.4%), and R. magnacalcarata (7.2%) and the springtail F. candida (6.3%). However, we also detected unexpectedly high levels of candidate HGT in the ant O. biroi and the wasp T. pretiosum (Supplementary Table 4). To evaluate a potential link between elevated candidate HGT levels and parthenogenesis in hexapods we additionally quantified candidate HGT in a number of published genomes of sexual species from the same order or superfamily as the parthenogenetic species. This revealed a similarly high proportion of candidate HGT in these sexual relatives, suggesting that whatever the cause of high candidate HGT in these taxa is a general characteristic of hexapods and not linked to the switch to parthenogenetic reproduction. The precise nature of these putative foreign genes in hexapod genomes remains unknown.

Gene Family Expansions
Most genome papers scan for expansions of specific gene families. Such expansions are then discussed in the light of the focal species' biology. The expansion of specific gene families per se is thus generally a species-specific trait (Lespinet et al. 2002) that is not related to parthenogenesis. To our knowledge, the only example of a gene family expansion that could be directly associated with parthenogenesis is the diversification of the RNA silencing machinery of TEs in bdelloid rotifers (Flot et al. 2013). TEs are expected to evolve reduced activity rates in parthenogenetic hosts (see Transposable Elements), and an improved RNA silencing machinery could be the mechanism underlying such reduced activity rates.
Functionally mitotic parthenogenesis might facilitate variation in gene copy numbers between homologous chromosomes as a consequence of relaxed constraints on chromosome pairing. Gene family expansions (and contractions) could therefore be more extensive and be retained more frequently in parthenogenetic than sexual species. To test this hypothesis, an overall comparison of gene family expansions in sexual and parthenogenetic sister species is needed (see Supplementary Materials S5). Four studies have surveyed gene family expansions in parthenogenetic species as well as in (sometimes distantly related) sexual counterparts, but these studies found no differences between reproductive modes (Faddeeva-Vakhrusheva et al. 2017;Warren et al. 2018;Lindsey et al. 2018;Schiffer et al. 2018). However, only 2 of the 4 studies are based on parthenogens with functionally mitotic parthenogenesis (i.e., where chromosome pairing is not required), and additional studies are therefore needed to address the question of whether parthenogenesis affects gene family expansions.

Gene Loss
Parthenogenetic animals are predicted to lose genes underlying sexual reproduction traits, including male-specific traits and functions (e.g., male-specific organs, spermatogenesis), as well as female traits involved in sexual reproduction (e.g., pheromone production, sperm storage organs) (van der Kooi and Schwander 2014). In the absence of pleiotropic effects, gene loss is expected due to mutation accumulation in the absence of purifying selection maintaining sexual traits, as well as to directional selection to reduce costly sexual traits (Schwander et al. 2013). Some gene loss consistent with these predictions is documented. For example, the sex determination genes xol-1 and tra-2 are missing in the nematode D. coronatus (Hiraki et al. 2017). Furthermore, genes believed to be involved in male functions harbor an excess of deleterious mutations in the wasp Leptopilina clavipes , which could represent the first step towards the loss of these genes. However, a similar excess of deleterious mutations in genes with (presumed) male-specific functions was not detected in the Amazon molly P. formosa (Warren et al. 2018).
Parthenogenetic species are further predicted to lose genes specific to some meiotic processes that no longer take the place during egg production (Schurko and Logsdon 2008). The genes involved in meiosis have been studied in 7 parthenogens (the Amazon molly, 4 bdelloid rotifers, and 2 Diploscapter nematodes). There was no apparent loss of meiosis genes in the Amazon molly P. formosa (Warren et al. 2018). The majority of meiosis genes were recovered in bdelloid rotifers (Flot et al. 2013;Nowell et al. 2018), and furthermore, nearly all of the meiosis genes missing in bdelloid rotifers were also absent in sexual monogonont rotifers, indicating that meiosis gene loss is not associated with the evolution of obligate parthenogenesis (Hanson et al. 2013). Both Diploscapter nematodes also lack certain meiosis genes, but it is unknown whether these genes are also missing in sexual relatives (Hiraki et al. 2017;Fradin et al. 2017). As much as the idea is appealing, there does not seem to be any support for the predicted loss of meiotic genes in functionally mitotic parthenogens. We note that the lack of our understanding of meiosis on the molecular level outside of a few model organisms (particularly yeast and Caenorhabditis elegans) makes the interpretation of gene loss (or absence thereof) difficult. This is best illustrated by the fact that losses of presumed "core" meiosis genes have also been reported in different sexual species, where meiosis is clearly fully functional (Tvedte et al. 2017).
In summary, some gene loss consistent with the loss of different sexual functions has been reported in several parthenogenetic species. However, there is no striking pattern relative to sexuals, and a clear interpretation of gene loss in parthenogenetic species is problematic because the function of the vast majority of genes is unknown in these non-model organisms.

Conclusions
We reanalyzed 26 published genomes of parthenogenetic animals to identify genomic features that are characteristic of parthenogenetic animals in general. Many of the original genome studies highlighted one or a few specific features in their focal species and suggested that they might be linked to parthenogenesis. However, our analyses combined with reviewing published results show that none of these genome features appear to be a general consequence of parthenogenesis, given that none of them was replicated across even a majority of analyzed species.
The variation among genomes of parthenogenetic species is at least in part due to species-or lineage-specific traits. But variation among the features detected in the published single-genome studies is also generated by differences in the methods used. Such differences are often less obvious, yet they can be critical in our assessment of genome diversity among animals. In this work, we thus re-analyzed several key genome features with consistent methods. To minimize the potentially confounding effects of differences in assembly quality, we have utilized methods that analyze directly sequencing reads. For example, re-estimating heterozygosity levels directly from reads of each species allowed to show a strong effect of hybrid origin, but not of cellular mechanism of parthenogenesis ( Figure 2). Another advantage of using the same methods for each species is that it diminishes the "researcher degrees of freedom" (Simmons et al. 2011;Gelman and Loken 2013;Wallach et al. 2018). For example, the analysis of polyploid genomes requires choosing methods to call heterozygosity and ploidy. By providing a common framework among species, we have shown that homoeolog divergence is very diverse among polyploid asexuals.
We have identified hybrid origin as the major factor affecting heterozygosity levels across all parthenogenetic animal species with available genomic data. This is consistent with the conclusions of 2 studies that focused on individual asexual lineages: hybridization between diverse strains explains heterozygosity in Meloidogyne root-knot nematodes and in Lineus ribbon worms (Ament-Velásquez et al. 2016;Szitenberg et al. 2017). This rule applies more generally to all the species analyzed with known transitions to parthenogenesis, but it is important to highlight that all the non-hybrid species in our dataset are hexapods. Thus in principle, the low heterozygosity could be a hexapod specific pattern, for example, due to high mitotic gene conversion rates in hexapods. The taxonomic range of the sequenced species is wide but we are missing several clades rich in parthenogenetic species, such as mites or annelids (Norton and Palmer 1991;Veresoglou et al. 2015). These clades would be useful foci for future genomic studies of parthenogenetic species. Independently of the findings of such future studies, our results suggest that heterozygosity loss via meiosis and/or gene conversion plays a significant and highly underappreciated role in the evolution of parthenogenetic species of intraspecific origin and support the theoretical argument that one of the main benefits of sex could be the masking of recessive deleterious mutations (referred to as "complementation") (Archetti 2004;Archetti 2010). Conversely, high rates of heterozygosity loss could also allow for the purging of deleterious mutations, as in highly selfing species (e.g., Glémin and Ronfort 2013;Szövényi et al. 2014). Such purging could help explain why most of the genome-scale studies did not find support for the theoretical expectation that parthenogenetic reproduction should result in increased rates of deleterious mutation accumulation (see Mutation Accumulation and Positive Selection). More generally, given the major differences in genome evolution for parthenogens of intraspecific versus hybrid origin, our study calls for future theoretical approaches on the maintenance of sex that explicitly consider the loss versus the maintenance of heterozygosity in asexuals.
In our evaluation of the general consequences of parthenogenesis, we were not able to take 2 key aspects into account: survivorship bias of parthenogenetic lineages, and characteristics of sexual ancestors. How often new parthenogenetic lineages emerge from sexual ancestors is completely unknown, but it has been speculated that in some taxa parthenogenetic lineages might emerge frequently, and then go extinct rapidly because of negative consequences of parthenogenesis. In other words, parthenogens that would exhibit the strongest consequences of parthenogenesis, as predicted by theoretical models, are expected to go extinct the fastest. Such transient parthenogens remain undetected in natural populations, because research focuses on parthenogenetic species or populations, and not on rare parthenogenetic females in sexual populations. Indeed, most of the species included in our study have persisted as parthenogens for hundreds of thousands to millions of years. They might thus be mostly representative of the subset of lineages that suffer weaker consequences of parthenogenesis. Finally, the key constraint for identifying consequences of parthenogenesis is that almost none of the published genome studies of parthenogenetic animals included comparisons to close sexual relatives. This prevents the detection of specific effects of parthenogenesis, controlling for the variation among sexual species-which is extensive for all of the genome features we analyzed and discussed in our study. Overall, despite the importance of recombination rate variation for understanding the evolution of sexual animal genomes (e.g., Stapley et al. 2017;Bachtrog 2013), the genome-wide reduction 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 lineagespecific patterns, differences according to the origin of parthenogenesis, and survivorship bias of parthenogenetic lineages.