Introduction:
The evolutionary success of sex has been one of the greatest puzzles evolutionary biologists have pondered for over a century (Weismann, 1889). Although, theoretically, sexually-reproducing organisms should be outcompeted by the asexual organisms due to inherent costs of sex, they remain predominant in the world of eukaryotes. Numerous hypotheses have been proposed to explain this paradox (see Shcherbakov, 2010 for one recent hypothesis, reviewed in e.g., Hartfield & Keightley, 2012). Some are linked to inherent ecological properties of reproductive modes, such as the penalization of asexuals for their inability to rapidly generate novel variants, thus increasing their intraspecific competition (Doncaster, Pound, & Cox, 2000; Peck, Yearsley, & Waxman, 1998) or vulnerability to pathogens (Hamilton, 1980). Another popular class of theories emphasizes an increased extinction risk of the clones (ultimately decreasing their diversification), or the gradual deterioration of their genomes (hereafter referred to as mechanisms of ‘clonal decay’). These are based on the hypothesis that every asexual lineage will acquire harmful mutations, and selection will not be efficient enough to maintain the fittest genome (Kondrashov, 1982) because individual mutations may not escape their genomic background due to a whole-genome linkage, leading to accumulation of mutations in a stochastic (Muller’s ratchet) or a deterministic (Kondrashov’s hatchet) manner (Kondrashov, 1993). Consistent with this hypothesis, most asexual lineages appear to be evolutionarily short-lived, and occupy terminal positions on the tree of life (reviewed in Janko, Drozd, Flegr, & Pannell, 2008; Schwander & Crespi, 2009). Several genomic studies indeed have found evidence of higher rates of accumulation of nonsynonymous mutations in asexuals (e.g. Henry, Schwander, & Crespi, 2012; Hollister et al., 2015; Howe & Denver, 2008; S. G. Johnson & Howard, 2007; Neiman, Hehman, Miller, Logsdon, & Taylor, 2010; Paland & Lynch, 2006).
Despite broad consensus that accumulation of mutation should ultimately lead to the demise of asexual lineages, it remains unclear whether the process proceeds fast enough to exterminate the asexuals before they become competitive threats to the sexual populations. Some argue that higher rates of mutation accumulation can considerably harm the clonal populations even over a short period of time (e.g. Neiman et al., 2010), especially in synergy with other mechanisms (Pound, Cox, & Doncaster, 2004; West, Lively, & Read, 1999). However, a non-negligible number of studies have failed to confirm higher loads of mutations in asexual organisms (e.g. Allen, Light, Perotti, Braig, & Reed, 2009; Brandt et al., 2017; Naito & Pawlowska, 2016; Pellino et al., 2013; Warren et al., 2018), and several studies have also demonstrated comparable or higher net diversification rates in asexual clades than in the sexual ones, suggesting these lineages may not inherently suffer from a decreased diversification, an increased extinction, or both (Fontaneto, Tang, Obertegger, Leasi, & Barraclough, 2012; M. T. J. Johnson, FitzJohn, Smith, Rausher, & Otto, 2011; Liu et al., 2012).
Such contradictions do not necessarily invalidate detrimental effects of mutation accumulation as a factor, but they indicate that the full explanation is more complex. For example, asexuals may have some mechanisms that temporarily delay the process of clonal decay, such as ‘minimal sex’, ameiotic recombination, gene conversions, increasing ploidy with consequent genome refreshment, and masking mutations with heterozygous states (e.g. Flot et al., 2013; Halligan & Keightley, 2003; Loewe & Lamatsch, 2008; Sémon & Wolfe, 2007). It is also possible that relaxation of selection stemming from asexuality is not drastic enough to universally cause the expected trends in mutation accumulation patterns (Allen et al., 2009). It thus remains a major challenge in evolutionary biology to identify the short-term-acting mechanisms that protect the sexuals from competitive displacement by the asexuals.
So called sexual-asexual complexes, i.e., taxa in which related sexual and asexual forms coexist, play prominent roles in studies of the persistence of sex since both types of organisms enter the same evolutionary arena for mutation, selection, and drift (Birky & Barraclough, 2009; Janko, Drozd, & Eisner, 2011). In some complexes, asexual lineages may emerge from the sexual populations very dynamically (Janko et al., 2012), which poses an additional question regarding the role of mutation accumulation process in saving sex: even if individual clonal lineages become debilitated in a relatively short term, whole asexual populations, composed of multiple frequently emerging clones, may escape extinction and pose a serious short-term or long-term threat to the sexual counterparts (Paland, Colbourne, & Lynch, 2005). Interestingly, Janko et al. (2008) showed that continuous influx of new clones into an asexual population causes inevitable replacement of the old clones by drift without invoking an age-dependent decay of the clonal lineages. This implies that ages attained by the individual clonal lineages may not correspond to their theoretical maxima determined by an age-dependent decay, but rather by the rate of influx of the new clones into an asexual community of a finite size. This drift-like clonal turnover process proved sufficient to explain some phylogenetic data (Janko, 2014), suggesting that in some cases clones may vanish from a population before mutation accumulation or other such processes compromise their fitness.
Even less is known about the direct impact of mutation accumulation on clone fitness. Indeed, if increased mutation accumulation in asexuals is relevant to the maintenance of sex, it should have phenotypic consequences when relatively old asexual lineages ‘… suffer in comparison to younger asexual lineages (and sexuals) with regard to important traits such as mitochondrial function and the rate of offspring production and population growth’ (Neiman et al. 2010). Unfortunately, studies investigating fitness-related traits are scarce, and their results are contradictory; age-dependent fitness deterioration has been demonstrated in some lineages of viruses, bacteria, and yeasts (Andersson & Hughes, 1996; Chao, 1990; Goddard, Godfray, & Burt, 2005), but not in water frogs with genomes transmitted clonally for over 25 ky (Guex, Hotz, & Semlitsch, 2002).
In this study, we focused on the Cobitis taenia hybrid complex of the European bottom dwelling loaches, and examined whether the genomes of sexual, young asexual, and old asexual lineages differed in accumulation of nonsynonymous mutations, and whether the old and the young clones differed in traits that impacted their fitness.
Our model taxon consists of several parapatrically distributed sexual species that inhabit Europe and coexist over much of the continent with their hybrids that regularly establish successful, persistent clonal lineages with a gynogenetic reproductive mode (i.e. females lay unreduced eggs that require activation by sperm from sexual males to trigger development). Production of clonal eggs is achieved via “pre-meiotic endomitosis” when oogonial chromosomes are duplicated, and subsequent meiotic divisions do not yield any variability since crossovers occur among identical copies of chromosomes (Juchno, Arai, Boroń, & Kujawa, 2016, Dedukh et al. 2019). Phylogenetic and crossing experiments (Choleva et al., 2012; Janko, Culling, Rab, & Kotlik, 2005; Janko et al., 2012; Janko, Vasil’ev, et al., 2005; Janko, Bohlen, et al. 2007) showed that range shifts related to the Pleistocene climatic oscillations have repeatedly placed the Cobitis parental species into a secondary reproductive contact, provoking the dynamic formation of hybrid clones. Their diversity is created in a two-step process (Fig. 1 insert). First, diploid clonal F1 females form through crosses of the parental species C. elongatoides with either C. taenia orC. tanaitica (Choleva et al. 2012; Janko et al. 2018); for simplicity, the haploid genomes of those species have been hereafter denoted as E, T, and N, respectively, and hybrid forms with particular combinations of the parental genome (so-called ‘genomotypes’) have been denoted by these letters, e.g. ETT indicates a triploid hybrid with one elongatoides and two taenia genomes. In the second step, mating with sexual males may occasionally result in the incorporation of sperm genomes leading to the formation of secondary triploid clones. Tetraploids are occasionally produced in an analogous way, but generally die before maturity (Janko, Bohlen, et al., 2007; Juchno et al., 2014). Contemporary Cobitis populations are composed of a mixture of clones of quite different evolutionary ages. Some clones are very recent. Other clones originated during the early-Holocene with establishment of the secondary contact zones between the parental species. Molecular dating using mtDNA data indicates that the origin of the oldest known lineage dates to ~0.3 Mya (Janko, Culling, et al., 2005, Majtánová et al., 2016). This lineage, referred to as the ‘hybrid clade I’, originated in the Ponto-Caspian region as a diploid C. elongatoides-tanaitica (EN) hybrid, and during its subsequent expansion it has incorporated additional genomes on multiple occasions, resulting in a monophyletic clonal assemblage of the EN, the EEN, the ENN, and the ETN clonal genomotypes. Owing to a lack of evidence on the observation of asex → sex transitions, the emergence of asexuality in the Cobitis appears to be a unidirectional process that has been more or less continuously occurring for at least hundreds of thousands of years, which should be enough to observe the putative effects of clonal decay (Loewe & Lamatsch, 2008). Mechanisms ensuring such long-term coexistence with sexual forms have not been properly understood, although Kotusz et al. (2014) and Bartoš et al. (2019) documented niche segregation between the Cobitis forms, suggesting some role of ecological factors.
Dynamic emergence of asexuality, the presence of diploid and polyploid clones, and a relatively well-resolved evolutionary history and ecology all make the Cobitis an attractive model for investigation of genomic and phenotypic consequences of asexual reproduction, allowing, among other objectives, the disentanglement of polyploidy effects from asexuality and/or hybridization. Here, we simultaneously analyzed the genetic variability in nuclear and mitochondrial genes, together with fitness data, to test predictions of mutation accumulation hypotheses. Specifically, we investigated whether asexual genomes carried traces of relaxed purifying selection which may be manifested byd N/d S ratios closer to 1 along the ‘asexual’ branches of the phylogenetic trees (Wertheim et al., 2014), by generally higher rates of accumulation of the nonsynonymous mutations, or by more radical amino acid substitutions in clonal genomes (Pellino et al., 2013; Sharbrough, Luse, Boore, Logsdon, & Neiman 2018). We further compared mutation frequency spectra to test whether the spread of novel nonsynonymous variants differed between the sexual and the asexual lineages (Hollister et al., 2015). Finally, we extended our analyses to the phenotypic data from additional individuals and tested whether traits related to the body-condition and fecundity tended to decay in the aging clones.
Material & Methods
Origins and numbers of analyzed specimens are listed in Figure 1 and Tables 1 & 2. The taxonomical identities of all fish used in this study were determined with a set of previously published and validated species diagnostic markers including allozyme and microsatellite loci and flow-cytometric determination of ploidy (Janko, Flajšhans, et al., 2007; Janko et al., 2012).
Genomic data
Exome data were acquired using the exome capture method described in Janko et al. (2018) and briefly below. Using probes for targeted enrichment of gDNA loci we collected sequence data from nuclear exomic and mtDNA segments of 22 females belonging to three hybridizing sexual species, 26 hybrids of different ploidy, and the genome composition and one specimen of C. paludica as an outgroup (Tab. 1, Fig. 1). Samples of parental species represent all major phylogroups,sensuJanko, Culling, et al. (2005), to capture maximum coverage of intraspecific variability in our data. Hybrid specimens include two or more individuals from several independent clonal lineages characterized previously by Janko et al. (2012), to compare genetic variability within and between clonal lineages.
Sequencing and SNP calling
Isolated gDNA was sheared with Bioruptor, tagged by indices, pooled, and hybridized to custom-designed sequence-capture probes. Captured fragments were sequenced on an Illumina NextSeq in 75 bp paired-end (PE) mode. Fastq files were trimmed by the fqtrim tool (Pertea, 2015) with minimum read lengths of 20 bp, and 3’ end trimming when average base quality within a sliding window drops below 15, which is a commonly recommended threshold (e.g. Suren et al., 2016).
To further process the reads, we used a published Cobitisreference transcriptome (GGJF00000000; Janko et al., 2018), from which Janko et al. (2018) discarded potentially paralogous contigs, as described in Gayral et al. (2013). In short, we considered as potentially paralogous contigs those that possessed spurious heterozygosity patterns when the same heterozygotic positions occurred across distantly related species (including the outgroup). Such patterns are unlikely to originate from meaningful biological processes, and probably result from mapping of paralogous reads onto one contig present in the reference. Loci with such properties were removed from the reference. For the purposes of this study, we annotated the reference transcriptome with Blast2GO software (Conesa et al., 2005). Contigs were aligned against the nr database (18.11.2015) by Blastx 2.2.31. Hits with e-values < 0.0001 were accepted as significant. Otherwise, default settings were used. Subsequently, we identified the longest open reading frame (ORF) within each annotated contig using the getorf tool from Embassy package (v6.6.0.0).
We aligned reads to this reference with BWA MEM (Li & Durbin, 2009), followed by Picard tools v1.140 to mark duplicates (Broad Institute, http://broadinstitute.github.io/picard). Individuals’ variants were called with the GATK v3.4 HaplotypeCaller tool, and all individuals were jointly genotyped using the GenotypeGVCFs tool (DePristo et al., 2011; McKenna et al., 2010; Van der Auwera et al., 2013). Variant recalibration was based on a previously used database of species-diagnostic positions (Janko et al., 2018) representing a learning set for the variant quality score recalibration tool VariantRecalibrator, and all variants were then filtered with the ApplyRecalibration tool. All resulting high confidence SNPs with coverage ≥ 5 were transferred to a database using custom SQL/python scripts.
Four samples represent an experimental family including two C. elongatoides-taenia F1 hybrid individuals and both their parents, which allowed us to experimentally verify the efficiency of sequencing and SNP calling by comparing F1 hybrids to their parents. We also compared our NGS data with Sanger sequences of several loci obtained by Choleva et al. (2014) to assess the accuracy of our method.
Tests of Relaxed Selection in Mitochondrial genomes
Because haploid mtDNA data are not biased by the problem of phasing, we first analyzed mitochondrial genomes by mapping reads to publishedC. elongatoides mitochondrion (accession no. NC_023947.1) and applying the RELAX software (Wertheim et al. 2014) to individual haplotypes to test the hypothesis that selection is relaxed along asexual branches, resulting ind N/d S ratios closer to 1 than along sexual branches. Since the software requires phased sequences representing single ORF with no heterozygous sites or stop codons, we discarded all noncoding sequences (D-loop and tRNAs), and the reverse-oriented ND6 gene from our mtDNA reference. We also masked several observed heterozygous positions and stop codons with ’N’ using the HyPhy package (Pond, Frost, & Muse, 2005), since these may represent sequencing errors or mutated premature termination codons. In case of overlapping ORFs in several mitochondrial genes, we split the sequence at the start codon of the second ORF, and padded the 3’ end of the first ORF with ’N’ to keep the correct reading frame of the whole alignment.
The final alignment of coding mtDNA was uploaded to the PhyML server (v3.0 build 20120412; Guindon et al., 2010), and phylogenetic trees were computed with default settings. A substitution model was selected based on AIC using Smart Model Selection (SMS 1.8.1.; Lefort, Longueville, & Gascuel, 2017; Fig. 2, Table S1). We then used both the alignment and Newick tree file in HyPhy 3.8 to perform relaxed selection analysis in RELAX with the assumption of strict sex → asex transitions, which is justified based on known pathways of asexual hybridization in Cobitis (Choleva et al., 2012). The software requires two types of branches to be specified in each tree. As test branches, we selected those leading exclusively to asexual individuals (i.e. asexual branches). As reference branches, we considered the intraspecific tree branches within all three sexual species (i.e. sexual branches; Fig. 2). We left the interspecific branches unclassified, since sequence evolution at the interspecific level is likely to bear traces of stronger negative selection, as many mildly deleterious mutations segregating within species would not become fixed between species. Choleva et al. (2014) showed that C. tanaitica is fixed for introgressed elongatoides-derived mitochondrion and appears paraphyletic to C. elongatoides , potentially indicating repeated introgression events in the past. Therefore, we also left one internal branch within C. tanaitica unclassified, as we could not determine whether it represents an intra- or interspecific branch (see Fig. 2). We also left unclassified the branches leading to laboratory progeny used in F1 crosses, as they represent non-natural experimental strains with unknown reproductive modes in the F1 generation.
We compared the fit to data of the two models allowing site-specificd N/d Sratios, i.e. the null model assuming the same ratios for asexual (test) and sexual (reference) branches of the ML phylogenetic tree (see Fig. 2), and the alternative model assuming thatd N/d S ratios differ between types of branches.
Tests of Relaxed Selection in nuclear DNA
While mitochondrial data allowed the application of sophisticated methods (i.e. RELAX), we were also interested in trends in nuclear DNA, which is highly heterozygous in hybrids and thus had to be analyzed differently. Specifically, we inspected whether ratios of synonymous and nonsynonymous mutations differ between sexual and asexual populations by comparing the frequency spectra in each population, and testing whether radicality of nonsynonymous substitutions differs among genomotypes. Finally, as a more formal test of selection, we compared per-gene distributions of d N/d Sratios among datasets. Details are provided below.
Types of SNPs, N/S ratios, site frequency spectra, and radicality of amino acid substitutions
All nDNA analyses are based on ORFs that represent annotated protein-coding sequences with assigned GO terms. For methodological simplicity (see e.g. Burgarella et al., 2015), we constrained all downstream analyses to bi-allelic SNP positions, i.e. those with a maximum of two alternative alleles observed across all investigated specimens (note that less than ~0.001 of positions contained more than two variants, similar to observations by Ament-Velásquez et al., 2016). We further masked with ’N’ all triplets containing more than one polymorphic site since we could not be sure about their amino acid translation without knowing their phase. Also, to enter the SNP database, each site had to be successfully sequenced in at least 80% of individuals. After filtering, each SNP was categorized into one of the three following categories. First, we identified sites segregating within any of the three parental species (i.e. intraspecific sexual polymorphisms). Second, we scored SNP variants that were monomorphic within each species but differed between some pairs of species (i.e. fixed interspecific sexual polymorphisms). Finally, we noted the so-called ‘private asexual SNPs’, where all parental sexual species appeared invariant, but asexual hybrids carried apomorphic alleles. We assumed that such substitutions generally accumulated after the origin of hybrid clonal lineages (see below).
To determine synonymity of any mutation, we translated nucleotide sequences into amino acids (including premature stop codons) using R package seqinr and compared numbers of synonymous and nonsynonymous mutations among all SNP types (intraspecific and interspecific sexual, and private asexual). Although some tools exist for direct comparisons of selection efficiency and distributions of mutation fitness effects between populations (Eyre-Walker & Keightley, 2009) and have even been applied to sexual-asexual comparisons (e.g Hollister et al., 2015), some of their assumptions (segregation and recombination) make their usage for asexuals problematic, hence we did not apply them. Instead, we also estimated the frequency of each SNP in sexual populations and within each clonal lineage, and tested whether patterns differ between synonymous and nonsynonymous substitutions. For this analysis, we assumed that the ancestral state of each intraspecific or private asexual SNP corresponded to the variant present in the sister species, while the other allele was assumed to be mutated.
We also tested whether radicality of acquired nonsynonymous mutations could be higher in asexuals because lack of segregation in asexuals should relax selection against mutations with strong effects, especially when the clone possesses another functional copy of a given gene (Hollister et al. 2015; Sharbrough et al. 2018). We performed two tests of this hypothesis. We first calculated for all categories of SNPs the proportions of polymorphic sites carrying premature stop codons, which likely represent strong-effect mutations. Second, we scored the ”radicality” of nonsynonymous mutations with BLOSUM90 and PAM100 substitution matrices, and compared score distributions among sexual and asexual datasets. We have chosen these matrices due to low numbers of mutation per gene, which may leave some amino acid substitutions unobserved.
d N/d Sratio
Thed N/d Sratio, i.e. the ratio of the number of nonsynonymous substitutions per nonsynonymous site (d N) to the number of synonymous substitutions per synonymous site (d S), is a common technique used to detect selection in pairs of protein-coding sequences. Comparisons ofd N/d S ratios between pairs of sexual and asexual lineages is a suitable tool to address the efficiency of selection with the prediction that asexuals should possess higher d N/d S values due to relaxed purifying selection (e.g. Pellino et al., 2013). However, due to the hybrid origin of asexuals, theird N/d S values are not directly comparable to those from sexual species. Many of asexuals’ polymorphisms are inherited from substitutions fixed between parental species brought together by hybridization, and we know theird N/d S are lower than those of intraspecific polymorphisms segregating within parental species. Thus, we also calculatedd N/d S ratios for all possible elongatoides-taenia and elongatoides-tanaiticainterspecific pairs, which served as null distributions to be compared with real asexual hybrids. This approach is based on the rationale thatd N/d S ratios of simulated hybrids provide variability that would be expected in hybrids without the effects of any process operating on asexual genomes.
For comparisons ofd N/d Sdistributions within and between species and genomotypes, we selected ORFs with at least 50 intact codons fully resolved in all individuals, to maintain representative sequence lengths. We calculatedd N/d S ratios for all pairs of sequences using the Nei-Gojobori evolutionary pathway method and the Jukes-Cantor model as implemented in the Bioperl module Bio::Align::DNAStatistics (http://www.bioperl.org; Stajich et al., 2002). Since some sequence pairs differed in synonymous or nonsynonymous mutations, but not both, which would introduce zeros into numerators or denominators ofd N/d S ratios, we added a small constant to each estimate (0.01; Pellino et al., 2013).
Estimating fitness effects in old vs. recent clones
Inspection of SNPs, albeit very instructive, may only provide an indirect proxy for inferences about fitness decay in clones. The next logical step, albeit rarely undertaken, is to test for age-dependent fitness deterioration of clones. To make such a comparison, we examined additional samples from 7 sites in the Central European hybrid zone, where asexuals belonging to recent and ancient lineages coexist. All specimens used for phenotypic analyses were genetically assessed by standard markers to determine their genomotype (Janko et al. 2007).
Because diploid and triploid Cobitis hybrids differ in many traits (Bobyrev, Burmensky, Vasil’ev, Kriksunov, & Lebedeva, 2003; Juchno & Boron, 2010; Kotusz, 2008; Maciak et al., 2011), we restricted analyses to triploid females, and investigated traits related to body condition, growth, and fecundity, predicting more prominent fitness decay in older clones compared to recent ones (Neiman et al., 2010). For this part of the study we captured and analyzed 48 EEN females belonging to the ancient (hybrid clade I) and 52 females belonging to recent clonal lineages (EET and ETT genomotypes; Fig. 1, Tab. 2). To address seasonal variation, we performed sampling before and after the reproductive season (April and September, respectively). We caught the fish by electrofishing, sacrificed them by overdose of 2-di-phenoxyethanol, and fixed them in 4% formaldehyde after removing a piece of fin tissue for genotyping.
Each specimen was measured (SL; 1 mm accuracy), weighed (Wt; 1g) and dissected for internal organs: heart, gonad, liver, and spleen, which were weighed. The bodies were cooked, then re-measured and re-weighed after dissection of the vertebral column and entrails (We). Weight and length data were used to estimate standard body-condition coefficients: relative SL/Wt relationship (LWR) and Clark’s condition factor (CC = We/SL3; Ricker, 2010), followed by anatomical traits, such as the heart-somatic index (HI), reflecting heart functional capacity, spleen-somatic index (SSI) as a measure of immunocompetence, the hepato-somatic index (HSI) — a measure of energy reserves (Šimková et al., 2015), and vertebrae index (VI), which assesses skeleton ossification. Size of vertebrae was measured along the horizontal axis using digital calipers (accurate to 0.1 mm) of the three sequential vertebrae starting at the plane of the first dorsal pterygiophore. Age-growth dynamics were estimated according to Fedorčák, Koščo, Halačka, and Manko (2017) on the basis of the total number of annuli on the vertebral body (distances from the vertebra’s centre to each annual ring, and to the vertebra edge were measured to the nearest 0.005 mm). Growth dynamics were interpreted as an indication of inter-clonal competition ability for resources within the shared niche. Fecundity, either absolute (i.e. total number of pre-matured oocytes in a female) or relative (i.e. absolute fecundity per weight) was estimated as the number of exogenous vitellogenic oocytes in the gonad according to Halačka, Lusková, and Lusk (2000), and Juchno and Boron (2010). Differences in size-structure of oocytes among females were estimated according to Halačka et al. (2000) by measuring the diameters of 250 randomly chosen oocytes. We note that these latter measures were performed on a subset of fish of similar size and age to avoid biases due to age differences (Tab. 2).
To further test for fecundity differences, we allowed four and three hybrid females belonging to ancient (EEN) and recent (ETT) clones, respectively, to mate at will with C. elongatoides males under seminatural conditions. We kept each triploid hybrid female in an identical tank, together with a single male, with spawning substratum following the design of Bohlen (1999). After each spawning event, we counted eggs from each clutch and let them develop under standardized conditions until they hatched and reached the third fry stage, when we evaluated their survival rate.
Statistical Analyses of Data
We used R packages nlme (Pinheiro et al., 2016) and lme4 (Bates, Maechler, Bolker, & Walker, 2015) to analyze differences among fish types in clutch sizes, survival rates, SNP data, site frequency spectra, amino acid radicality, andd N/d S ratios. To address specifics of each dataset, we used several linear models including Linear Mixed Effects models (LME), Generalized Linear Models (GLM), and Generalized Linear Mixed Effects models (GLME). We obtained p-values for the effect in question by Likelihood ratio tests (LRT) against the null model. Details about these statistical methods are provided in Supplementary online materials.
Results
Molecular data
On average, we mapped 17,684,155 reads per sample (all samtools bitwise flags except 0x4) onto a reference Cobitis transcriptome (Janko et al., 2018). Within each of the 20,601 contigs, we identified open reading frames (ORF) altogether providing us with 3,945 nuclear ORFs (3,179,571 bp total length) with annotation and sufficient read coverage to identify SNP variants. Mapping the reads onto a published C. elongatoides mitogenome (NC_023947.1) sufficiently covered on average 82.1% of the reference, with most reads falling into coding regions (10,902 sites).
Test of Relaxed Selection on mtDNA
Smart model selection (SMS), used by the PhyML server, selected GTR + G as the best model of sequence evolution (Table S1); the resulting tree is depicted in Fig. 2. Comparison of the two models in RELAX by likelihood ratio test (LRT) did not indicate a preference for the alternative model over the null model (∆AICc = 0.3; K = 2.98, p = 0.129, LR = 2.30), which assumed no differences between sexual and asexual branches, and suggested that 86.72% of sites are explained byd N/d S values of 0.02, 7.08% by values of 0.44, and the remaining 6.2% by values of 1.56.
Nuclear DNA
SNP categories, relative ages of clones, and reliability of SNP calling
Across all nuclear contigs, we discovered 161,139 bi-allelic SNPs, of which 149,248 were variable within the ingroup. Focusing specifically on annotated ORFs, we categorized 28,878 sites as intraspecific polymorphisms segregating within any of the sexual species. These SNPs were characterized by the prevalence of nonsynonymous substitutions over synonymous ones. We further scored 10,532 fixed interspecific SNP variants, which, in contrast to the intraspecific SNPs, showed the opposite trend, with prevailing synonymous substitutions (Fig. 3A). This difference between SNP types was highly significant (GLME with binomial family of error distribution, LRT against null model, p -value = 7.874699e-05). Moreover, the site-frequency spectra within each species showed that mutations inducing synonymous changes were significantly more frequent than mutations changing the amino-acid (LME, LRT p = 1.61937e-107; Fig. S1).
We further noted 8,840 of the so-called ‘private asexual SNPs’. The hypothesis that such sites most likely accumulated after the hybrid origin of clonal lineages is corroborated by the fact that their proportion in hybrids’ genomes was tightly correlated with the mtDNA distance of any clone from its nearest sexual haplotype (Pearson’sr = 0.94, p -value = 7.135e-14). Specifically, the ‘hybrid clade I’, i.e. the oldest asexual lineage, possessed the most divergent mtDNA haplotypes, and the highest proportion of private SNPs (~1.5 % of all SNPs). Hybrids from Central Europe with a putative Holocene origin (Janko et al., 2012) possessed haplotypes identical or closely related to those segregating in current sexual populations, and simultaneously had a moderate proportion of private SNPs (~0.2-0.4 %), while the two experimental F1 hybrids possessed the least divergent mtDNA haplotypes, and had a very low proportion of private SNPs, which probably represent false-positives due to sequencing errors. The rarity of private mutations in F1s also corroborates the quality of performed SNP calling, which is further strengthened by the fact that our SNP calling recovered 100% matches to sequences of several mitochondrial and nuclear genes previously published by Choleva et al. (2014).
Patterns of heterozygosity
Studied specimens systematically differed in their level of heterozygosity. C. elongatoides possessed the highest average heterozygosity of all parental species, and C. taenia the lowest (Fig. 3B), consistent with its small effective population size (Janko et al., 2018). Asexual individuals had significantly higher genome-wide heterozygosity than sexuals (Fig. 3B; t-test p = 8.963e-14) with 98.5 – 99.8% of ‘private asexual SNPs’ occurring in heterozygous states. The vast majority of fixed interspecific SNPs were in heterozygous states in hybrids’, consistent with expectations. However, every hybrid possessed a small proportion of sites with only one parental allele; i.e. sites with so-called loss of heterozygosity (LOH). The proportion of LOH sites in each individual was rather low, and correlated significantly with putative ages of clones, i.e. either with K2P mtDNA distances from the nearest sexual individual (Pearson’s r = 0.974, 95% c.i. = 0.940-0.988, p-value = 3.016e-16) or with proportion of private asexual mutations (Pearson’s r = 0.955, 95% c.i. = 0.902-0.980, p-value = 3.286e-14). Hence, the proportion of LOH sites was near zero in experimental F1 hybrids, and highest (~10%) inelongatoides-tanaitica (Fig. 4).
Mutation accumulation and radicality in clones and polyploids
Private asexual SNPs had a significantly higher probability of being nonsynonymous than SNPs fixed among parental species (GLME with binomial family of error distribution, LRT p = 0.0008). Both intraspecific sexual polymorphisms and private asexual SNPs had similar proportions of nonsynonymous SNPs (binomial GLME, LRT p = 0.369; Fig. 3A). Nevertheless, we noticed that old clones had a significantly lower proportion of nonsynonymous SNPs than young clones (binomial GLME, LRT p = 0.00055) or than sexual species when measured on intraspecific sexual polymorphisms (binomial GLME, LRT p = 0.12). Altogether, these data indicate no tendency toward a higher nonsynonymous mutation rate in clones, and even suggests the possibility of lower rates in the old clonal lineage.
To evaluate the impact of polyploidization on mutation accumulation, we compared diploid and triploid asexuals within both major types of hybrids (i.e. ET vs EET and ETT in elongatoides-taenia and EN vs. EEN in elongatoides-tanaitica hybrid types). We noticed that triploids possessed significantly higher proportions of private asexual SNPs compared to diploids (Fig. 3A; ET vs EET & ETT: t -test p = 0.006; EN vs EEN: t -test p = 0.001). We also noticed that variance in proportions of private asexual SNPs was generally higher among triploids than among diploids, but these differences were not significant due to small sample sizes (ET vs EET & ETT: F test p = 0.06; EN vs EEN: F test p = 0.4). To compare rates of nonsynonymous mutation accumulation between di- and triploids we applied the binomial GLME separately to elongatoides-taenia and elongatoides-tanaitica hybrids, treating individuals within groups as random factors. Private asexual SNPs in diploid ET hybrids accumulated a significantly higher proportion of nonsynonymous mutations compared to EET or ETT triploids (binomial GLME LRT p = 0.047) but the differences between EN diploids and EEN triploids were not significant (binomial GLME LRT p = 0.42).
When testing the accumulation of large-effect mutations, we found as expected, that interspecific SNPs are significantly less likely to carry premature stop codon mutations than both intraspecific and private asexual SNPs (contingency table, p < 1e-15). Asexuals, when considered as one group, had significantly higher incidence of premature stop codons than intraspecific sexual datasets (contingency table, p = 0.016). In our analyses of old and recent clones, we found that the incidence of stop codons did not differ between intraspecific sexual SNPs and private asexual SNPs in old clones (contingency table, p = 0.99), but their incidence was significantly elevated in young clones (contingency table, p = 2.4e-6).
We also scored the ”radicality” of nonsynonymous mutations using the PAM100 substitution matrix, and compared score distributions among sexual and asexual datasets (Fig. 5A). Again, fixed interspecific nonsynonymous substitutions were significantly less radical than intraspecific segregating mutations (LME, LRT p = 0.0009) or private asexual SNPs (LME, LRT p = 0.004). However, the radicality of mutations segregating within sexual species (intraspecific polymorphisms) did not differ from that of nonsynonymous private asexual SNPs (LME, LRT p = 0.339). We observed no significant differences between old and young clones (LME, LRT p = 0.084). Analogous results were obtained with the BLOSUM90 matrix (data not shown).
No evidence for relaxed selection fromd N/d S ratios and mutation frequency spectra
A more formal test of selection was based on the per-gene distributions ofd N/d Sratios in all genomotypes (Fig. 5B). A subset of 1,993 ORFs passed a threshold of 50 or more codons fully resolved in all investigated individuals. To address the hybrid origin of asexuals’ variability, we compared the d N/d S ratio of asexuals with that estimated from elongatoides-taenia orelongatoides-tanaitica sequence pairs, as well as estimates from laboratory F1 hybrid offspring. The usefulness of this approach was corroborated by the fact that elongatoides-taenia pairwised N/d S values did not differ from those in F1-generation ET hybrids (LME, LRT p = 0.43).
Most genes in all comparisons had values below 1, consistent with the influence of negative selection. Within-species pairwise comparisons yielded significantly higher values than interspecific comparisons, and values observed within asexual genomotypes (LME, LRT p < 1e-8). However, thed N/d Sratios calculated from interspecific pairs of sequences were significantly higher than values from old clones (LME, LRT p = 0.01), young clones (LME, LRT p = 0.02) or asexuals in general (LME, LRT p = 0.004), which was at odds with predictions based on mutation accumulation hypotheses. Moreover, we found no differences between old and young clones (LME, LRT p = 0.35).
Looking at the distribution of private asexual SNPs we found that nonsynonymous substitutions often appeared as singletons, while synonymous substitutions were significantly more likely to be shared by two or more members of the same clone, (Fig. 5C; binomial GLME, LRT p = 5.7e-6). Since the ancient clonal group (hybrid clade I) had the highest sample size, we also looked at its mutation frequency spectra and noticed that nonsynonymous substitutions are again significantly less frequent than synonymous ones (Poisson GLM, LRT p = 1.1e-6).
No evidence for fitness deterioration in old clones
We found a significant correlation between fish length and weight (ANCOVA, F1, 56 = 115.13, p = 3.346e-15) and also a significant effect of season, with spring samples having higher SSI (t95 = 8.34, p << 0.000), HSI (t97 = 3.75, p = 0.000), and GSI (t97 = 4.740, p << 0.000) values, and gonads with higher counts of generally smaller oocytes compared to autumn samples (Fig. 6A). Some effect of location was also present as fishes from the Złotnica River site had smaller relative SL/Wt ratios (F1, 56 = 33.56, p = 3.291e-07) and slower growth rates (F1, 251 = 103.3101, p < 2.2e-16) than those from the remaining sites.
Taking the aforementioned effects into account, we found no differences between old and young clones in condition-related measures: LWR (ANCOVA, F1, 56 = 0.0018, p = 0.9667), CC (t98 = -0.838, p = 0.404; Fig. 6B), HI (t90 = 0.265, p = 0.791), HSI (t98 = 1.08, p = 0.282), VI (t58 = 1.713, p = 0.92), or in growth rates (F1, 251 = 2.11, p = 0.147), except for a significantly higher SSI in old clones (t95 = -2.05, p = 0.043). Looking at traits related to fecundity, the old clonal lineage possessed significantly higher GSI values than the younger clonal lineages (in spring: t17 = 2.617, p = 0.018, in fall: t78 = 5.968, p < 0.000). Egg counts followed this pattern, showing statistical significance in Spring samples: absolute fecundity: t17 = 3.069, p = 0.007; relative fecundity: t17 = 2.723, p = 0.014 (Fig. 6C,D). We also noted differences between lineages in egg size structure. The old clonal lineage possessed larger counts of smaller eggs than the younger clones in the Spring period, while the opposite trend was observed in the autumn (permutation Kolmogorov-Smirnov test, p << 0.0001 in both seasons; Fig. 6A).
A reproductive experiment corroborated the aforementioned differences, since females of the oldest clone had significantly larger clutches compared to those from the younger clones (average of 394 and 176 eggs per clutch, respectively; GLME with female ID as a random factor, Z = 4.364 p = 1.28e-05). We noted generally low survival rates (on average ~16% of eggs per clutch hatched and survived till third fry stage), consistent with previous experiments (Janko, Bohlen, et al., 2007; Juchno et al., 2014; Bohlen pers.comm.), but we found no significant differences between old and recent clones in survival (zero-inflated GLME, Z = -0.001, p = 0.999).
Discussion
Various mechanisms have been proposed to explain why sexual reproduction evolved and became a dominant force in nature but no clear answer has been obtained. We tested whether the mutational deterioration of clones, one of the most cited candidate processes, may serve as an efficient short-term mechanism explaining the stable coexistence of sexual species with clonal lineages, some reaching ages of up to 300 ky.
No evidence for accumulation of deleterious mutations and fitness decay in asexuals:
Negative selection prevailed among >3000 studied nuclear and mitochondrial loci, as evident from higher proportions of nonsynonymous and more radical mutations, from higherd N/d S values observed within the parental species rather than between the parental species, as well as from generally lower intrapopulation frequencies of nonsynonymous mutations compared to the synonymous ones. This is in line with the hypothesis that nonsynonymous or radical intraspecies polymorphisms may be often mildly deleterious, and may be consequently removed by negative selection before reaching fixation.
We further found that clonal genomes accumulated novel mutations in a time-dependent manner, so that the older clonal lineages had higher incidence of private asexual SNPs compared to the younger clones or the experimental F1 hybrids, consistent with mutation accumulation models. Nevertheless, the way in which new mutations accumulated contrasts with expectations of a reduced efficacy of selection in asexuals. Namely, despite hundreds of thousands of years of clonality in some lineages, any potential relaxation of selection was too weak to provide detectable traces in the asexuals’ mitochondria, and there was no evidence for increased proportions of nonsynonymous or radical mutations among private asexual SNPs in the ancient clone. We also comparedd N/d S ratios of real asexuals to values of hybrids simulated from the sequences of the sexual individuals. This approach, which accounts for the effect of fixed interspecific polymorphisms, showed lower-than-expectedd N/d S values in natural clones, with no difference between the old and the recent clones, again contrasting with predictions of mutation accumulation hypotheses.
The fitness-related data are also at odds with the prediction of clonal decay since the ancient ‘hybrid clade I’ outperformed the young clones in some important traits, such as GSI and fecundity. Given that proportions of asexuals in mixed sexual-asexual populations often correlate with their competitive success (Hellriegel & Reyer, 2000; Vrijenhoek & Parker, 2009), the evolutionary success of the oldestCobitis clonal lineage is further indicated by its occupation of the largest distribution range of all Cobitis clones (Janko, Culling, et al., 2005) and its ability to achieve higher densities than co-occurring younger clones at many sites (Janko et al., 2012). We have to note that old and recent clones used in our study have different genomic compositions (elongatoides-tanaitica versuselongatoides-taenia genome combinations, respectively), suggesting their dissimilarity may also reflect inherited differences between parental species rather than asexuality-specific processes. Nevertheless, the absence of a detectable effect on fitness continues to indicate that even the oldest Cobitis clones successfully coexist and compete for resources with recent clones and sexual species.
How do asexuals escape mutational deterioration?
Altogether, the acquired data provide no evidence for reduced efficacy of the selection in asexuals, suggesting that mutation accumulation may not notably hamper the demographic and the competitive performances of the Cobitis clones even after 300 ky of asexuality. What are the reasons behind these observations?
First, the observed patterns were unlikely to be methodological artifacts, since our study covered most of the Cobitisdistribution range and provided a relatively high number of SNPs that were validated in the experimental F1 hybrids and previously published sequences, and estimated values of the phenotypic traits did not deviate from previous reports of the Cobitis genomotypes, including the polyploids (Halačka et al., 2000; Juchno & Boron, 2010). The inherent biological nature of the studied organisms may have affected the data since comparisons between asexuals and their sexual counterparts (or simulated hybrids) conflate the variability frozen by hybridization in the distant past with that segregating within the contemporary sexual populations. Additionally, comparisons between the old and the recent clones potentially suffer owing to the presence of different species in their ancestries. Despite these pitfalls, we clearly documented that mutations did accumulate in time, as predicted by the clonal decay model, but ~300 ky of asexuality did not leave detectable traces of the deleterious mutation accumulation or decreased clonal fitness. This study adds to the growing list of publications that have failed to detect traces of clonal decay (e.g. Pellino et al., 2013; Warren et al., 2018).
What does the growing amount of negative evidence tell us about the validity of mutation-accumulation theories for evolution of sex and clonality? Commonly, negative results have been attributed to special mechanisms enabling investigated asexuals to avoid mutation accumulation, such as deviations from a strict clonality including minimal sex, ameiotic recombination, gene conversions, beneficial effects of polyploidization, or an increased efficiency of DNA repair (e.g. Maciver, 2016; Roach & Heitman, 2014; Schön & Martens, 2003; Warren et al., 2018). However, given the available data, it is unlikely that any such mechanism may have efficiently slowed down the clicks of the ratchet in the Cobitis .
Refuting cryptic sexuality in predominantly clonal organisms is notoriously difficult (e.g. Birky, 2010), but available evidence argues against its role in the Cobitis . Extensive crossing experiments have never documented allelic segregation (Choleva et al., 2012; Janko, Bohlen, et al., 2007; Janko et al., 2018), and Majtánová et al. (2016) reported remarkable stability of the hybrids’ karyotypes with no large-scale chromosomal restructuring or recombination. The current study also shows that the vast majority of the private hybrid SNPs occur in heterozygous states (between 98.5%–99.8%), which corroborates the hypothesized clonal reproduction, with new mutations occurring on one chromosome with little possibility of recombination between the homologues or coalescence between the alleles. Nevertheless, the existence of LOH at fixed interspecific positions indicates that asexual organisms’ alleles may occasionally recombine, convert, or vanish due to hemizygous deletions. However, three lines of evidence argue against a major impact of LOH on the speed of ratchet clicks in theCobitis .
First, even the oldest Cobitis clones accumulated LOH events in only ~10% of the studied genes over ~300 ky, which strikingly contrasts with other asexual organisms, where such processes have been hypothesized to interfere with the accumulation of deleterious mutations. For example, Tucker, Ackerman, Eads, Xu, and Lynch (2013) and Sunnucks, England, Taylor, and Hales (1996) showed orders of magnitude of higher rates of conversions & deletions of individual genes, and even entire chromosomal arms. Second, the efficiency of LOH events in erasing potentially deleterious mutations is expected to increase with the clonal age. This is because, rare LOH events may likely happen on different genes from new mutations in the recent clones, while in the older clones any emerging LOH event has a higher chance of occurring in the genes with existing mutations. Therefore, if LOH events were to attenuate the ratchet by overwriting deleterious mutations, we would expect an exponential correlation between the proportions of LOH and the private asexual SNPs because LOH would erase a disproportionately higher proportion of acquired mutations in old clones compared to the younger ones. As evidenced in Fig. 4, our data offer no support for such a deviation from linearity (note that we focused on diploid clones here since the presence of three homologs in triploids complicates the detection of LOH). Finally, only < 1.5 % of the private SNPs were in homozygous state, probably resulting from LOH events. If we assume that LOH events are more or less symmetrical with respect to preserving or losing the new or ancestral allele, it suggests that only ~3% of the private asexual mutations could have been converted. Of course, if new mutations tend to be deleterious, the selection might have preferred one direction of LOH but altogether, our data indicate that genome restructuring likely has no major effect on the clicks of the ratchet.
Polyploidization is another mechanism that may refresh a clonal genome by temporarily masking deleterious mutations, although it also frees gene copies to mutate, ultimately increasing the per-genome mutation rate (Otto & Whitton, 2000). Consistent with this expectation, we found higher numbers of private asexual SNPs in triploid loaches than in diploids but the differences were relatively small. This may suggest that the contemporary triploids may have spent considerable parts of their history as diploid clones before polyploidization, leaving relatively little time for accumulation of extra mutations. Importantly, we found no evidence for the increased nonsynonymous mutation rates in triploids. Hence, although polyploidization likely plays a very important role in clonal evolution, at least in the Cobitis , its major benefits do not seem to stem from a deleterious mutation masking; after all, the diploid clones have persisted at the Balkan together with their triploid derivatives at least since the last interglacial period (Janko, Culling, Rab, & Kotlik, 2005). Major benefits of triploidy may instead relate to other effects, such as metabolic changes invoked by the modified cell architecture, the altered dosage of parental genomes, or the modified crosstalk between the allopolyploids’ genomes due to different stoichiometric relations between transcription factors and their binding sites (Bartoš et al., 2019; Beukeboom & Vrijenhoek, 1998; Maciak et al, 2011).
The null model of clone persistence assumes delayed effects of accumulated mutations
In summary, although we may not rule out the existence of some mechanisms counteracting the effects of clonality, the aforementioned candidate mechanisms do not appear strong enough to stop or considerably slow the ratchet in the Cobitis . Instead, we propose that the apparent absence of signs of clonal decay, as in the Cobitis and other similar cases, may have a much simpler reason, which we propose as a null hypothesis for any tests of clonal decay. Namely, the simplest plausible explanation is that the relaxation of a purifying selection is not strong enough to allow sufficient mutation accumulation that will compromise the performance of clones over several hundreds of thousands of generations.
Consistent with this explanation, we found similar frequency spectra of non-/synonymous mutations in the sexual and the asexual populations, indicating that selection efficiently restricts the spread of putatively deleterious mutations even in the clones. The counter-intuitive observations of lower nonsynonymous mutation loads and the incidence of premature stop codons in old clones may thus be explained by the classical population genetic theory, since a selection-based removal of the deleterious mutations requires varying amounts of time depending on the selection coefficient, the population size, and the genetic background. It follows that recent clones, despite having acquired lower absolute numbers of mutations, will possess a higher fraction of nonsynonymous private asexual mutations, or more radical SNPs.
Different mutation loads in the old and the young clones may thus reflect a time-lag necessary to remove the slightly deleterious mutations (S. G. Johnson & Howard 2007) rather than some special mechanisms that alleviate the mutation accumulation. Paradoxically, the delay in fitness decay may be caused directly by clonal reproduction because the maintenance of heterozygosity may mask the negative effects of acquired deleterious recessives (Halligan & Keightley, 2003) until new mutations accumulate in high numbers (Otto, 2007) or become exposed to selection in the homozygous states (Guex et al., 2002; Leslie & Vrijenhoek, 1978; Leslie & Vrijenhoek, 1980). This was in line with recent reanalysis of published sequences of asexual animals (Janko et al., 2011), which showed that significant deviations from neutrality, indicative of selection against aging clones, could only be detected in asexual complexes when clones achieved substantial ages beyond 1 million generations.
Hence, even if a clonal lineage were to accumulate deleterious mutations since its birth, it may take considerable time before they start to negatively affect the evolutionary life span of that clone. In younger asexual complexes, the replacement of clonal lineages may be dominated by a more neutral drift-like process of the clonal turnover, which does not depend on the effects of accumulated mutations (Janko et al., 2008). By analogy to a genetic drift, the clonal turnover predicts that a major portion of clonal diversity may be formed by the young clonal lineages, which did not have enough time to spread, while ancient lineages will be less diverse but more geographically widespread. Interestingly, both predictions have been met in the Cobitis and in some other taxa (Janko et al., 2012; Janko, 2014; Quattro et al., 1992).
Conclusions
To date, the existence of successful clones in nature has usually been explained by their special abilities to avoid accumulation of the deleterious mutations. Our study shows that even if clones clearly accumulate mutations as they age, the relaxation of selection due to asexuality may not be sufficient to cause their fitness to deteriorate, even after hundreds of thousands of generations of asexuality. Thus, although mutation accumulation may ultimately drive any clone to extinction, individual clones may be replaced by more recent clones in a drift-like process of the clonal turnover if the influx rates of new clones are high (Janko, 2008). Clones in such complexes may simply not be given an opportunity to become old enough for accumulated mutations to start to manifest their effects. Successful and widespread clones may therefore often represent merely temporary winners in the ongoing turnover race.
This may seem trivial, but it has potentially serious implications for our understanding of the persistence of sex: Mutation accumulation may be relevant to the survival of clones after they reach some substantial age, but before they do so, the persistence of sexual species in mixed sexual-asexual complexes must rely on different mechanisms that offer short-term advantages to sex.
Acknowledgements
We would like to thank Jan Šipoš and Petr Pyszko for valuable advice on linear mixed effect models, Jacek Stefaniak for preparing maps of samples, and Jan Eisner for help with the fitting of second-order polynomial models to our data. This research was supported by The Czech Science Foundation project no. 13-12580S to KJ, University of Ostrava project no. SGS29/PřF/2016 to JKoč and by the National Science Centre in Poland project no. DEC-2011/03/B/NZ8/02095 to JKot. Computational resources were provided by the ELIXIR-CZ project (LM2015047), part of the international ELIXIR infrastructure.
References
Allen, J. M., Light, J. E., Perotti, M. A., Braig, H. R., & Reed, D. L. (2009). Mutational Meltdown in Primary Endosymbionts: Selection Limits Muller’s Ratchet. PLoS ONE , 4 (3), e4969. doi:10.1371/journal.pone.0004969
Ament-Velásquez, S. L., Figuet, E., Ballenghien, M., Zattara, E. E., Norenburg, J. L., Fernández-Álvarez, F. A., … Galtier, N. (2016). Population genomics of sexual and asexual lineages in fissiparous ribbon worms (Lineus, Nemertea): hybridization, polyploidy and the Meselson effect. Molecular Ecology , 25(14), 3356–3369. https://doi.org/10.1111/mec.13717
Andersson, D. I., & Hughes, D. (1996). Muller’s ratchet decreases fitness of a DNA-based microbe. Proceedings of the National Academy of Sciences , 93 (2), 906–907.
Bates, D., Maechler, M., Bolker, B., & Walker, S. (2015). Fitting Linear Mixed-Effects Models Using lme4. Journal of Statistical Software, 67 (1), 1-48. doi:10.18637/jss.v067.i01.
Birky, C. W. (2010). Positively Negative Evidence for Asexuality.Journal of Heredity , 101 (Supplement 1), S42–S45. doi:10.1093/jhered/esq014
Birky, C. W., & Barraclough, T. G. (2009). Asexual Speciation. In I. Schön, K. Martens, & P. van Dijk (Eds.), Lost Sex: The Evolutionary Biology of Parthenogenesis (pp. 201–216). Dordrecht, Heidelberg, London, New York: Springer.
Bobyrev, A., Burmensky, V., Vasil’ev, V., Kriksunov, E., & Lebedeva, E. (2003). Coexistence of triploid and diploid forms of spined loach,Cobitis taenia : a model-based approach. Folia Biologica ,51 Suppl , 55–60.
Bartoš, O., Röslein, J., Kotusz, J., Paces, J., Pekárik, L., Petrtýl, M., … & Juchno, D. (2019). The legacy of sexual ancestors in phenotypic variability, gene expression and homoeolog regulation of asexual hybrids and polyploids. Molecular Biology and Evolution , 36(9), 1902–1920. doi: 10.1093/molbev/msz114
Beukeboom, L. W., & Vrijenhoek, R. C. (1998). Evolutionary genetics and ecology of sperm dependent parthenogenesis. Journal of Evolutionary Biology , 11, 755–782. https://doi.org/10.1007/s000360050117
Brandt, A., Schaefer, I., Glanz, J., Schwander, T., Maraun, M., Scheu, S., & Bast, J. (2017). Effective purifying selection in ancient asexual oribatid mites. Nature Communications , 8(1), 873. doi: 10.1038/s41467-017-01002-8
Burgarella, C., Gayral, P., Ballenghien, M., Bernard, A., David, P., Jarne, P., … Glémin, S. (2015). Molecular Evolution of Freshwater Snails with Contrasting Mating Systems. Molecular Biology and Evolution , 32(9), 2403–2416. https://doi.org/10.1093/molbev/msv121
Chao, L. (1990). Fitness of RNA virus decreased by Muller’s ratchet.Nature , 348 (6300), 454–455. doi:10.1038/348454a0
Choleva, L., Janko, K., De Gelas, K., Bohlen, J., Šlechtová, V., Rábová, M., & Ráb, P. (2012). Synthesis of Clonality and Polyploidy in Vertebrate Animals by Hybridization Between Two Sexual Species.Evolution , 66 (7), 2191–2203. doi:10.1111/j.1558-5646.2012.01589.x
Choleva, L., Musilova, Z., Kohoutova-Sediva, A., Paces, J., Rab, P., & Janko, K. (2014). Distinguishing between Incomplete Lineage Sorting and Genomic Introgressions: Complete Fixation of Allospecific Mitochondrial DNA in a Sexually Reproducing Fish (Cobitis ; Teleostei), despite Clonal Reproduction of Hybrids. PLoS ONE , 9 (6), e80641. doi:10.1371/journal.pone.0080641
Conesa, A., Gotz, S., Garcia-Gomez, J. M., Terol, J., Talon, M., & Robles, M. (2005). Blast2GO: a universal tool for annotation, visualization and analysis in functional genomics research.Bioinformatics, 21 (18), 3674–3676. https://doi.org/10.1093/bioinformatics/bti610
Dedukh, D., Majtánová, Z., Marta, A., Pšenička, M., Kotusz, J., Klíma, J., … Janko, K. (2019). Parthenogenesis as a solution to hybrid sterility: the mechanistic basis of meiotic distortions in clonal and sterile hybrids. BioRxiv , 663112. https://doi.org/10.1101/663112
DePristo, M. A., Banks, E., Poplin, R., Garimella, K. V., Maguire, J. R., Hartl, C., … Daly, M. J. (2011). A framework for variation discovery and genotyping using next-generation DNA sequencing data.Nature Genetics , 43 (5), 491–498. doi:10.1038/ng.806
Doncaster, C. P., Pound, G. E., & Cox, S. J. (2000). The ecological cost of sex. Nature , 404 (6775), 281–285. doi:10.1038/35005078
Eyre-Walker, A., & Keightley, P. D. (2009). Estimating the Rate of Adaptive Molecular Evolution in the Presence of Slightly Deleterious Mutations and Population Size Change. Molecular Biology and Evolution , 26 (9), 2097–2108. doi:10.1093/molbev/msp119
Fedorčák, J., Koščo, J., Halačka, K., & Manko, P. (2017). Growth differences in different biotypes of the hybrid complex of Cobitis elongatoides × Cobitis tanaitica (Actinopterygii: Cypriniformes: Cobitidae) in the Okna River (Danube River basin), Slovakia. Acta Ichthyologica Et Piscatoria , 47 , 125–132. doi:10.3750/AIEP/02059
Flot, J.-F., Hespeels, B., Li, X., Noel, B., Arkhipova, I., Danchin, E. G. J., … Van Doninck, K. (2013). Genomic evidence for ameiotic evolution in the bdelloid rotifer Adineta vaga . Nature , 500(7463), 453–457. https://doi.org/10.1038/nature12326
Fontaneto, D., Tang, C. Q., Obertegger, U., Leasi, F., & Barraclough, T. G. (2012). Different Diversification Rates Between Sexual and Asexual Organisms. Evolutionary Biology , 39 (2), 262–270. doi:10.1007/s11692-012-9161-z
Gayral, P., Melo-Ferreira, J., Glémin, S., Bierne, N., Carneiro, M., Nabholz, B., … Galtier, N. (2013). Reference-Free Population Genomics from Next-Generation Transcriptome Data and the Vertebrate–Invertebrate Gap. PLoS Genetics , 9 (4). doi:10.1371/journal.pgen.1003457
Goddard, M. R., Godfray, H. C. J., & Burt, A. (2005). Sex increases the efficacy of natural selection in experimental yeast populations.Nature , 434 (7033), 636–640. doi:10.1038/nature03405
Guex, G.-D., Hotz, H., & Semlitsch, R. D. (2002). Deleterious alleles and differential viability in progeny of natural hemiclonal frogs.Evolution , 56 (5), 1036–1044.
Guindon, S., Dufayard, J.-F., Lefort, V., Anisimova, M., Hordijk, W., & Gascuel, O. (2010). New algorithms and methods to estimate maximum-likelihood phylogenies: assessing the performance of PhyML 3.0.Systematic Biology , 59 (3), 307–321. doi:10.1093/sysbio/syq010
Halačka, K., Lusková, V., & Lusk, S. (2000). Fecundity of Cobitis elongatoides in the Nová Říše Reservoir. Folia Zoologica , (49), 141–150.
Halligan, D. L., & Keightley, P. D. (2003). How many lethal alleles?Trends in Genetics: TIG , 19 (2), 57–59. doi:10.1016/S0168-9525(02)00045-8
Hamilton, W. D. (1980). Sex versus Non-Sex versus Parasite.Oikos , 35 (2), 282–290. doi:10.2307/3544435
Hartfield, M., & Keightley, P. D. (2012). Current hypotheses for the evolution of sex and recombination. Integrative Zoology , 7(2), 192–209. doi: 10.1111/j.1749-4877.2012.00284.x
Henry, L., Schwander, T., & Crespi, B. J. (2012). Deleterious Mutation Accumulation in Asexual Timema Stick Insects. Molecular Biology and Evolution , 29 (1), 401–408. doi:10.1093/molbev/msr237
Hollister, J. D., Greiner, S., Wang, W., Wang, J., Zhang, Y., Wong, G. K.-S., … Johnson, M. T. J. (2015). Recurrent loss of sex is associated with accumulation of deleterious mutations inOenothera . Molecular Biology and Evolution , 32 (4), 896–905. doi:10.1093/molbev/msu345
Howe, D., & Denver, D. (2008). Muller’s Ratchet and compensatory mutation in Caenorhabditis briggsae mitochondrial genome evolution. BMC Evolutionary Biology , 8 (1), 62. doi:10.1186/1471-2148-8-62
Itono, M., Morishima, K., Fujimoto, T., Bando, E., Yamaha, E., & Arai, K. (2006). Premeiotic endomitosis produces diploid eggs in the natural clone loach, Misgurnus anguillicaudatus (Teleostei: Cobitidae).Journal of Experimental Zoology Part A: Comparative Experimental Biology , 305A (6), 513–523. doi:10.1002/jez.a.283
Janko, K., Bohlen, J., Lamatsch, D., Flajšhans, M., Epplen, J. T., Ráb, P., … Šlechtová, V. (2007). The gynogenetic reproduction of diploid and triploid hybrid spined loaches (Cobitis : Teleostei), and their ability to establish successful clonal lineages—on the evolution of polyploidy in asexual vertebrates. Genetica ,131 (2), 185–194.
Janko, K., Drozd, P., & Eisner, J. (2011). Do clones degenerate over time? Explaining the genetic variability of asexuals through population genetic models. Biology Direct , 6 , 17.
Janko, K., Drozd, P., Flegr, J., & Pannell, J. R. (2008). Clonal turnover versus clonal decay: a null model for observed patterns of asexual longevity, diversity and distribution. Evolution ,62 (5), 1264–1270.
Janko, K., Flajšhans, M., Choleva, L., Bohlen, J., Šlechtová, V., Rábová, M., … others. (2007). Diversity of European spined loaches (genus Cobitis L.): an update of the geographic distribution of the Cobitis taenia hybrid complex with a description of new molecular tools for species and hybrid determination.Journal of Fish Biology , 71 , 387–408.
Janko, K. (2014). Let us not be unfair to asexuals: their ephemerality may be explained by neutral models without invoking any evolutionary constraints of asexuality. Evolution , 68 (2), 569–576. doi:10.1111/evo.12293
Janko, K., Culling, M. A., Rab, P., & Kotlik, P. (2005). Ice age cloning-comparison of the Quaternary evolutionary histories of sexual and clonal forms of spiny loaches (Cobitis ; Teleostei) using the analysis of mitochondrial DNA variation. Molecular Ecology ,14 (10), 2991–3004.
Janko, K., Kotusz, J., De Gelas, K., Slechtová, V., Opoldusová, Z., Drozd, P., … Baláž, M. (2012). Dynamic formation of asexual diploid and polyploid lineages: multilocus analysis of Cobitisreveals the mechanisms maintaining the diversity of clones. PloS One , 7 (9), e45384. doi:10.1371/journal.pone.0045384
Janko, K., Pačes, J., Wilkinson-Herbots, H., Costa, R. J., Roslein, J., Drozd, P., … Choleva, L. (2018). Hybrid asexuality as a primary postzygotic barrier between nascent species: On the interconnection between asexuality, hybridization and speciation. Molecular Ecology , (27), 249–263. doi:10.1111/mec.14377
Janko, K., Vasil’ev, V. P., Rab, P., Rabova, M., Slechtova, V., & Vasil’eva, E. D. (2005). Genetic and morphological analyses of 50-chromosome spined loaches (Cobitis , Cobitidae, Pisces) from the Black Sea basin that are morphologically similar to C. taenia , with the description of a new species. Folia Zoologica ,54 (4), 405–420.
Johnson, M. T. J., FitzJohn, R. G., Smith, S. D., Rausher, M. D., & Otto, S. P. (2011). Loss of Sexual Recombination and Segregation Is Associated with Increased Diversification in Evening Primroses.Evolution , 65 (11), 3230–3240. doi:10.1111/j.1558-5646.2011.01378.x
Johnson, S. G., & Howard, R. S. (2007). Contrasting patterns of synonymous and nonsynonymous sequence evolution in asexual and sexual freshwater snail lineages. Evolution , 61 (11), 2728–2735. doi:10.1111/j.1558-5646.2007.00233.x
Juchno, D., Arai, K., Boroń, A., & Kujawa, R. (2016). Meiotic chromosome configurations in oocytes of Cobitis taenia and its polyploid hybrids. Ichthyological Research , 1–4. doi:10.1007/s10228-016-0556-1
Juchno, D., & Boron, A. (2010). Fecundity of the spined loach,Cobitis taenia (Pisces, Cobitidae) and natural allopolyploids ofCobitis from a diploid-polyploid population. Folia Zoologica . Retrieved from http://www.highbeam.com/doc/1P3-2013942551.html
Juchno, D., Jabłońska, O., Boroń, A., Kujawa, R., Leska, A., Grabowska, A., … Lao, M. (2014). Ploidy-dependent survival of progeny arising from crosses between natural allotriploid Cobitis females and diploid C. taenia males (Pisces, Cobitidae). Genetica ,142 (4), 351–359. doi:10.1007/s10709-014-9779-0
Kondrashov, A. S. (1993). Classification of Hypotheses on the Advantage of Amphimixis. Journal of Heredity , 84 (5), 372–387.
Kotusz, J. (2008). Morphological Relationships Between Polyploid Hybrid Spined Loaches of the Genus Cobitis (Teleostei: Cobitidae) and their Parental Species. Annales Zoologici , 58 , 891–905. doi:10.3161/000345408X396800
Kotusz, J., Popiołek, M., Drozd, P., De Gelas, K., Šlechtová, V., & Janko, K. (2014). Role of parasite load and differential habitat preferences in maintaining the coexistence of sexual and asexual competitors in fish of the Cobitis taenia hybrid complex.Biological Journal of the Linnean Society , 113 (1), 220–235. doi:10.1111/bij.12329
Lefort, V., Longueville, J.-E., & Gascuel, O. (2017). SMS: Smart Model Selection in PhyML. Molecular Biology and Evolution ,34 (9), 2422–2424. doi:10.1093/molbev/msx149
Li, H., & Durbin, R. (2009). Fast and accurate short read alignment with Burrows-Wheeler transform. Bioinformatics (Oxford, England) ,25 (14), 1754–1760. doi:10.1093/bioinformatics/btp324
Liu, H.-M., Dyer, R. J., Guo, Z.-Y., Meng, Z., Li, J.-H., & Schneider, H. (2012). The Evolutionary Dynamics of Apomixis in Ferns: A Case Study from Polystichoid Ferns. Journal of Botany , 2012 . doi:10.1155/2012/510478
Loewe, L., & Lamatsch, D. K. (2008). Quantifying the threat of extinction from Muller’s ratchet in the diploid Amazon molly (Poecilia formosa ). BMC Evolutionary Biology , 8(1), 88. https://doi.org/10.1186/1471-2148-8-88
Maciak, S., Janko, K., Kotusz, J., Choleva, L., Boroń, A., Juchno, D., … Konarzewski, M. (2011). Standard Metabolic Rate (SMR) is inversely related to erythrocyte and genome size in allopolyploid fish of the Cobitis taenia hybrid complex. Functional Ecology ,25 (5), 1072–1078. doi:10.1111/j.1365-2435.2011.01870.x
Maciver, S. K. (2016). Asexual Amoebae Escape Muller’s Ratchet through Polyploidy. Trends in Parasitology , 32(11), 855–862. https://doi.org/10.1016/j.pt.2016.08.006
Majtánová, Z., Choleva, L., Symonová, R., Ráb, P., Kotusz, J., Pekárik, L., & Janko, K. (2016). Asexual Reproduction Does Not Apparently Increase the Rate of Chromosomal Evolution: Karyotype Stability in Diploid and Triploid Clonal Hybrid Fish (Cobitis , Cypriniformes, Teleostei).PloS One , 11 (1), e0146872. doi:10.1371/journal.pone.0146872
McKenna, A., Hanna, M., Banks, E., Sivachenko, A., Cibulskis, K., Kernytsky, A., … DePristo, M. A. (2010). The Genome Analysis Toolkit: a MapReduce framework for analyzing next-generation DNA sequencing data. Genome Research , 20 (9), 1297–1303. doi:10.1101/gr.107524.110
Naito, M., & Pawlowska, T. E. (2016). Defying Muller’s Ratchet: Ancient Heritable Endobacteria Escape Extinction through Retention of Recombination and Genome Plasticity. MBio , 7 (3). doi:10.1128/mBio.02057-15
Neiman, M., Hehman, G., Miller, J. T., Logsdon, J. M., & Taylor, D. R. (2010). Accelerated Mutation Accumulation in Asexual Lineages of a Freshwater Snail. Molecular Biology and Evolution , 27 (4), 954–963. doi:10.1093/molbev/msp300
Otto, S. P. (2007). The evolutionary consequences of polyploidy.Cell , 131 (3), 452–462. doi:10.1016/j.cell.2007.10.022
Otto, S. P., & Whitton, J. (2000). Polyploid Incidence and Evolution.Annual Review of Genetics . https://doi.org/10.1146/annurev.genet.34.1.401
Paland, S., Colbourne, J. K., & Lynch, M. (2005). Evolutionary history of contagious asexuality in Daphnia pulex . Evolution ,59 (4), 800–813.
Paland, S., & Lynch, M. (2006). Transitions to Asexuality Result in Excess Amino Acid Substitutions. Science , 311 (5763), 990–992. doi:10.1126/science.1118152
Peck, J. R., Yearsley, J. M., & Waxman, D. (1998). Explaining the geographic distributions of sexual and asexual populations.Nature , 391 (6670), 889–892. doi:10.1038/36099
Pellino, M., Hojsgaard, D., Schmutzer, T., Scholz, U., Hörandl, E., Vogel, H., & Sharbel, T. F. (2013). Asexual genome evolution in the apomictic Ranunculus auricomus complex: examining the effects of hybridization and mutation accumulation. Molecular Ecology ,22 (23), 5908–5921. doi:10.1111/mec.12533
Pertea, G. (2015). fqtrim: v0.9.4 release. Zenodo . doi:10.5281/zenodo.20552
Pinheiro, J., Bates, D., DebRoy, S., Sarkar, D., & R Core Team. (2016). Nlme: Linear and Nonlinear Mixed Effects Models. R package version 3.1. http://CRAN.R-project.org/package=nlme.
Pond, S. L. K., Frost, S. D. W., & Muse, S. V. (2005). HyPhy: hypothesis testing using phylogenies. Bioinformatics (Oxford, England) , 21 (5), 676–679. doi:10.1093/bioinformatics/bti079
Pound, G. E., Cox, S. J., & Doncaster, C. P. (2004). The accumulation of deleterious mutations within the frozen niche variation hypothesis.Journal of Evolutionary Biology , 17 (3), 651–662. doi:10.1111/j.1420-9101.2003.00690.x
Quattro, J. M., Avise, J. C., & Vrijenhoek, R. C. (1992). An ancient clonal lineage in the fish genus Poeciliopsis (Atheriniformes: Poeciliidae). Proceedings of the National Academy of Sciences , 89(1), 348–352. https://doi.org/10.1073/PNAS.89.1.348
Ricker, W. E. (2010). Computation and Interpretation of Biological Statistics of Fish Populations . Caldwell, New Jersey: The Blackburn Press.
Roach, K. C., & Heitman, J. (2014). Unisexual reproduction reverses Muller’s ratchet. Genetics , 198 (3), 1059–1069. doi:10.1534/genetics.114.170472
Schön, I., & Martens, K. (2003). No slave to sex. Proceedings of the Royal Society of London. Series B: Biological Sciences ,270 (1517), 827–833. doi:10.1098/rspb.2002.2314
Schwander, T., & Crespi, B. J. (2009). Twigs on the tree of life? Neutral and selective models for integrating macroevolutionary patterns with microevolutionary processes in the analysis of asexuality.Molecular Ecology , 18 (1), 28–42. doi:10.1111/j.1365-294X.2008.03992.x
Sémon, M., & Wolfe, K. H. (2007). Consequences of genome duplication.Current Opinion in Genetics & Development , 17(6), 505–512. https://doi.org/10.1016/j.gde.2007.09.007
Sharbrough, J., Luse, M., Boore, J. L., Logsdon Jr, J. M., & Neiman, M. (2018). Radical amino acid mutations persist longer in the absence of sex. Evolution , 72(4), 808-824.
Shcherbakov, V. P. (2010). Biological species is the only possible form of existence for higher organisms: the evolutionary meaning of sexual reproduction. Biology Direct ,5 , 14. doi:10.1186/1745-6150-5-14
Stajich, J. E., Block, D., Boulez, K., Brenner, S. E., Chervitz, S. A., Dagdigian, C., … Birney, E. (2002). The Bioperl toolkit: Perl modules for the life sciences. Genome Research , 12 (10), 1611–1618. doi:10.1101/gr.361602
Sunnucks, P., England, P. R., Taylor, A. C., & Hales, D. F. (1996). Microsatellite and chromosome evolution of parthenogenetic sitobion aphids in Australia. Genetics , 144(2), 747–756. Retrieved from http://www.ncbi.nlm.nih.gov/pubmed/8889535
Suren, H., Hodgins, K. A., Yeaman, S., Nurkowski, K. A., Smets, P., Rieseberg, L. H., … Holliday, J. A. (2016). Exome capture from the spruce and pine giga-genomes. Molecular Ecology Resources . doi: 10.1111/1755-0998.12570
Tucker, A. E., Ackerman, M. S., Eads, B. D., Xu, S., & Lynch, M. (2013). Population-genomic insights into the evolutionary origin and fate of obligately asexual Daphnia pulex . Proceedings of the National Academy of Sciences , 110(39), 15740–15745. https://doi.org/10.1073/pnas.1313388110
Van der Auwera, G. A., Carneiro, M. O., Hartl, C., Poplin, R., Del Angel, G., Levy-Moonshine, A., … DePristo, M. A. (2013). From FastQ data to high confidence variant calls: the Genome Analysis Toolkit best practices pipeline. Current Protocols in Bioinformatics ,43 , 11.10.1-33. doi:10.1002/0471250953.bi1110s43
Warren, W. C., García-Pérez, R., Xu, S., Lampert, K. P., Chalopin, D., Stöck, M., … Schartl, M. (2018). Clonal polymorphism and high heterozygosity in the celibate genome of the Amazon molly. Nature Ecology and Evolution, 2 (4), 669–679. https://doi.org/10.1038/s41559-018-0473-y
Weismann, A. (1889). The significance of sexual reproduction in the theory of natural selection. In E. . Poulton, S. Schönland, & A. . Shipley (Eds.), Essays upon heredity and kindred biological problems. Vol. 1. (pp. 255–332). Oxford: Clarendon Press.
West, S. ., Lively, C. ., & Read, A. . (1999). A pluralist approach to sex and recombination. Journal of Evolutionary Biology ,12 (6), 1003–1012. doi:10.1046/j.1420-9101.1999.00119.x
[dataset] Kočí J., Röslein J., Pačes J., Kotusz J., Halačka K., Koščo J., Fedorčák J., Iakovenko N., Janko K.; 2019; No evidence for accumulation of deleterious mutations and fitness degradation in clonal fish hybrids: Abandoning sex without regrets; Institutional server (IMG CAS, link pending)
Data accessibility
Raw sequencing data will be available from a server at IMG CAS. Input files are provided as supplementary files.
Author contributions
KJ designed the experiments. KH, JF, JKoš, NI, performed the experiments. JKoč, JR, JP, JKot, KJ analysed data. JKoč, JKot, JR, and KJ contributed to writing the MS.
Figures