\documentclass[10pt]{article}

\usepackage{fullpage}
\usepackage{setspace}
\usepackage{parskip}
\usepackage{titlesec}
\usepackage[section]{placeins}
\usepackage{xcolor}
\usepackage{breakcites}
\usepackage{lineno}
\usepackage{hyphenat}





\PassOptionsToPackage{hyphens}{url}
\usepackage[colorlinks = true,
            linkcolor = blue,
            urlcolor  = blue,
            citecolor = blue,
            anchorcolor = blue]{hyperref}
\usepackage{etoolbox}
\makeatletter
\patchcmd\@combinedblfloats{\box\@outputbox}{\unvbox\@outputbox}{}{%
  \errmessage{\noexpand\@combinedblfloats could not be patched}%
}%
\makeatother


\usepackage{natbib}




\renewenvironment{abstract}
  {{\bfseries\noindent{\abstractname}\par\nobreak}\footnotesize}
  {\bigskip}

\titlespacing{\section}{0pt}{*3}{*1}
\titlespacing{\subsection}{0pt}{*2}{*0.5}
\titlespacing{\subsubsection}{0pt}{*1.5}{0pt}


\usepackage{authblk}


\usepackage{graphicx}
\usepackage[space]{grffile}
\usepackage{latexsym}
\usepackage{textcomp}
\usepackage{longtable}
\usepackage{tabulary}
\usepackage{booktabs,array,multirow}
\usepackage{amsfonts,amsmath,amssymb}
\providecommand\citet{\cite}
\providecommand\citep{\cite}
\providecommand\citealt{\cite}
% You can conditionalize code for latexml or normal latex using this.
\newif\iflatexml\latexmlfalse
\providecommand{\tightlist}{\setlength{\itemsep}{0pt}\setlength{\parskip}{0pt}}%

\AtBeginDocument{\DeclareGraphicsExtensions{.pdf,.PDF,.eps,.EPS,.png,.PNG,.tif,.TIF,.jpg,.JPG,.jpeg,.JPEG}}

\usepackage[utf8]{inputenc}
\usepackage[T2A]{fontenc}
\usepackage[polish,ngerman,english]{babel}



\usepackage{float}








\begin{document}

\title{No evidence for accumulation of deleterious mutations and fitness
degradation in clonal fish hybrid: Abandoning sex without regrets}



\author[1]{Jan Kočí}%
\author[2]{Jan Röslein}%
\author[3]{Jan Pačes}%
\author[4]{Jan Kotusz}%
\author[5]{Karel Halacka}%
\author[6]{Ján Koščo}%
\author[6]{Jakub Fedorčák}%
\author[1]{Nataliia Iakovenko}%
\author[2]{Karel Janko}%
\affil[1]{University of Ostrava}%
\affil[2]{Affiliation not available}%
\affil[3]{Institute of Molecular Genetics of the ASCR, v. v. i.}%
\affil[4]{University of Wroclaw}%
\affil[5]{Institute of Vertebrate Biology, v.v.i.}%
\affil[6]{University of Prešov}%


\vspace{-1em}



  
  \date{June 03, 2024}


\begingroup
\let\center\flushleft
\let\endcenter\endflushleft
\maketitle
\endgroup





\selectlanguage{english}
\begin{abstract}
Despite its inherent costs, sexual reproduction is ubiquitous in nature,
and the mechanisms to protect it from a competitive displacement by
asexuality remain unclear. Popular mutation-based explanations, like the
Muller's ratchet and the Kondrashov's hatchet, assume that purifying
selection may not halt the accumulation of deleterious mutations in the
non-recombining genomes, ultimately leading to their degeneration.
However, empirical evidence is scarce and it remains particularly
unclear whether mutational degradation proceeds fast enough to ensure
the decay of clonal organisms and to prevent them from outcompeting
their sexual counterparts. To test this hypothesis, we jointly analyzed
the exome sequences and the fitness-related phenotypic traits of the
sexually-reproducing fish species and their clonal hybrids, whose
evolutionary ages ranged from F1 generations to 300 ky. As expected,
mutations tended to accumulate in the clonal genomes in a time-dependent
manner. However, contrary to the predictions, we found no trend towards
increased non-synonymity of mutations acquired by clones, nor higher
radicality of their amino-acid substitutions. Moreover, there was no
evidence for fitness degeneration in the old clones compared to that in
the younger ones. In summary, although a purifying selection may still
be relaxed in the asexual genomes, our data indicate that its efficiency
is not drastically decreased. Even the oldest investigated clone was
found to be too young to suffer fitness consequences from a mutation
accumulation. This suggests that mechanisms other than mutation
accumulation may be needed to explain the competitive advantage of sex
in the short term.%
\end{abstract}%



\sloppy




\par\null

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\selectlanguage{ngerman}é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 `\ldots{} 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 \emph{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\selectlanguage{polish}ń, \selectlanguage{english}\& 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 \emph{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 \emph{C. elongatoides} with either \emph{C. taenia}
or\emph{C. 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 \emph{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 \textasciitilde{}0.3 Mya (Janko,
Culling, et al., 2005, Majt\selectlanguage{ngerman}ánová et al., 2016). This lineage, referred
to as the `hybrid clade I', originated in the Ponto-Caspian region as a
diploid \emph{C.} \emph{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 \selectlanguage{english}- sex transitions, the emergence of
asexuality in the \emph{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\selectlanguage{polish}š \selectlanguage{english}et al. (2019)
documented niche segregation between the \emph{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 \emph{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 by\emph{d}
\textsubscript{N}/\emph{d} \textsubscript{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\selectlanguage{polish}š\selectlanguage{english}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 \emph{C. paludica} as an outgroup (Tab. 1, Fig. 1).
Samples of parental species represent all major
phylogroups,\emph{sensu}Janko, 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
\emph{Cobitis}reference 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 \textless{} 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 \emph{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
published\emph{C. 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 in\emph{d} \textsubscript{N}/\emph{d}
\textsubscript{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
\emph{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 \emph{C. tanaitica} is fixed for
introgressed elongatoides-derived mitochondrion and appears paraphyletic
to \emph{C. elongatoides} , potentially indicating repeated
introgression events in the past. Therefore, we also left one internal
branch within \emph{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-specific\emph{d} \textsubscript{N}/\emph{d}
\textsubscript{S}ratios, 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
that\emph{d} \textsubscript{N}/\emph{d} \textsubscript{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 \emph{d} \textsubscript{N}/\emph{d}
\textsubscript{S}ratios 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 \textasciitilde{}0.001 of positions
contained more than two variants, similar to observations by
Ament-Vel\selectlanguage{ngerman}á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.

\emph{d} \textsubscript{N}/\emph{d} \textsubscript{S}ratio

The\emph{d} \textsubscript{N}/\emph{d} \textsubscript{S}ratio, i.e. the
ratio of the number of nonsynonymous substitutions per nonsynonymous
site (\emph{d} \textsubscript{N}) to the number of synonymous
substitutions per synonymous site (\emph{d} \textsubscript{S}), is a
common technique used to detect selection in pairs of protein-coding
sequences. Comparisons of\emph{d} \textsubscript{N}/\emph{d}
\textsubscript{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 \emph{d}
\textsubscript{N}/\emph{d} \textsubscript{S} values due to relaxed
purifying selection (e.g. Pellino et al., 2013). However, due to the
hybrid origin of asexuals, their\emph{d} \textsubscript{N}/\emph{d}
\textsubscript{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 their\emph{d} \textsubscript{N}/\emph{d}
\textsubscript{S} are lower than those of intraspecific polymorphisms
segregating within parental species. Thus, we also calculated\emph{d}
\textsubscript{N}/\emph{d} \textsubscript{S} ratios for all possible
\emph{elongatoides-taenia} and
\emph{elongatoides-tanaitica}interspecific pairs, which served as null
distributions to be compared with real asexual hybrids. This approach is
based on the rationale that\emph{d} \textsubscript{N}/\emph{d}
\textsubscript{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 of\emph{d} \textsubscript{N}/\emph{d}
\textsubscript{S}distributions 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 calculated\emph{d} \textsubscript{N}/\emph{d}
\textsubscript{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
(\emph{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
of\emph{d} \textsubscript{N}/\emph{d} \textsubscript{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 \emph{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/SL\textsuperscript{3}; 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 (\selectlanguage{polish}Š\selectlanguage{english}imkov\selectlanguage{ngerman}á
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\selectlanguage{polish}č\selectlanguage{ngerman}ák, Ko\selectlanguage{polish}šč\selectlanguage{english}o,
Hala\selectlanguage{polish}č\selectlanguage{english}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\selectlanguage{polish}č\selectlanguage{english}ka, Luskov\selectlanguage{ngerman}á, and Lusk (2000), and Juchno and Boron (2010).
Differences in size-structure of oocytes among females were estimated
according to Hala\selectlanguage{polish}č\selectlanguage{english}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 \emph{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, and\emph{d} \textsubscript{N}/\emph{d}
\textsubscript{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 \emph{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 \emph{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 by\emph{d}
\textsubscript{N}/\emph{d} \textsubscript{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, \emph{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's\emph{r} = 0.94, \emph{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 (\textasciitilde{}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 (\textasciitilde{}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. \emph{C. elongatoides} possessed the highest average
heterozygosity of all parental species, and \emph{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 (\textasciitilde{}10\%)
in\emph{elongatoides-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: \emph{t} -test p =
0.006; EN vs EEN: \emph{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 \textless{} 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 from\emph{d}
\textsubscript{N}/\emph{d} \textsubscript{S} ratios and mutation
frequency spectra

A more formal test of selection was based on the per-gene distributions
of\emph{d} \textsubscript{N}/\emph{d} \textsubscript{S}ratios 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
\emph{d} \textsubscript{N}/\emph{d} \textsubscript{S} ratio of asexuals
with that estimated from \emph{elongatoides-taenia}
or\emph{elongatoides-tanaitica} sequence pairs, as well as estimates
from laboratory F1 hybrid offspring. The usefulness of this approach was
corroborated by the fact that \emph{elongatoides-taenia}
pairwise\emph{d} \textsubscript{N}/\emph{d} \textsubscript{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 \textless{}
1e-8). However, the\emph{d} \textsubscript{N}/\emph{d}
\textsubscript{S}ratios 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, F\textsubscript{1, 56} = 115.13, p = 3.346e-15) and also a
significant effect of season, with spring samples having higher SSI
(t\textsubscript{95} = 8.34, p \textless{}\textless{} 0.000), HSI
(t\textsubscript{97} = 3.75, p = 0.000), and GSI (t\textsubscript{97} =
4.740, p \textless{}\textless{} 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\selectlanguage{polish}ł\selectlanguage{english}otnica River site had smaller relative SL/Wt ratios
(F\textsubscript{1, 56} = 33.56, p = 3.291e-07) and slower growth rates
(F\textsubscript{1, 251} = 103.3101, p \textless{} 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,
F\textsubscript{1, 56} = 0.0018, p = 0.9667), CC (t\textsubscript{98} =
-0.838, p = 0.404; Fig. 6B), HI (t\textsubscript{90} = 0.265, p =
0.791), HSI (t\textsubscript{98} = 1.08, p = 0.282), VI
(t\textsubscript{58} = 1.713, p = 0.92), or in growth rates
(F\textsubscript{1, 251} = 2.11, p = 0.147), except for a significantly
higher SSI in old clones (t\textsubscript{95} = -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: t\textsubscript{17} = 2.617, p = 0.018, in fall:
t\textsubscript{78} = 5.968, p \textless{} 0.000). Egg counts followed
this pattern, showing statistical significance in Spring samples:
absolute fecundity: t\textsubscript{17} = 3.069, p = 0.007; relative
fecundity: t\textsubscript{17} = 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 \textless{}\textless{}
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
\textasciitilde{}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 \textgreater{}3000 studied nuclear
and mitochondrial loci, as evident from higher proportions of
nonsynonymous and more radical mutations, from higher\emph{d}
\textsubscript{N}/\emph{d} \textsubscript{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 compared\emph{d}
\textsubscript{N}/\emph{d} \textsubscript{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-expected\emph{d}
\textsubscript{N}/\emph{d} \textsubscript{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
oldest\emph{Cobitis} clonal lineage is further indicated by its
occupation of the largest distribution range of all \emph{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 (\emph{elongatoides-tanaitica}
versus\emph{elongatoides-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 \emph{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 \emph{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
\emph{Cobitis}distribution 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 \emph{Cobitis}
genomotypes, including the polyploids (Hala\selectlanguage{polish}č\selectlanguage{english}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 \textasciitilde{}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\selectlanguage{ngerman}ö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 \emph{Cobitis} .

Refuting cryptic sexuality in predominantly clonal organisms is
notoriously difficult (e.g. Birky, 2010), but available evidence argues
against its role in the \emph{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 the\emph{Cobitis}
.

First, even the oldest \emph{Cobitis} clones accumulated LOH events in
only \textasciitilde{}10\% of the studied genes over
\textasciitilde{}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 \textless{}
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 \textasciitilde{}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 \emph{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\selectlanguage{polish}š \selectlanguage{english}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 \emph{Cobitis} . Instead, we propose that the
apparent absence of signs of clonal decay, as in the \emph{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 \emph{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 \selectlanguage{polish}Š\selectlanguage{english}ipo\selectlanguage{polish}š \selectlanguage{english}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\selectlanguage{polish}ř\selectlanguage{english}F/2016 to JKo\selectlanguage{polish}č \selectlanguage{english}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. \emph{PLoS ONE} , \emph{4} (3), e4969.
doi:10.1371/journal.pone.0004969

Ament-Vel\selectlanguage{ngerman}ásquez, S. L., Figuet, E., Ballenghien, M., Zattara, E. E.,
Norenburg, J. L., Fernández-Álvarez, F. A., \ldots{} Galtier, N. (2016).
Population genomics of sexual and asexual lineages in fissiparous ribbon
worms (Lineus, Nemertea): hybridization, polyploidy and the Meselson
effect. \emph{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. \emph{Proceedings of the National
Academy of Sciences} , \emph{93} (2), 906--907.

Bates, D., Maechler, M., Bolker, B., \& Walker, S. (2015). Fitting
Linear Mixed-Effects Models Using lme4. \emph{Journal of Statistical
Software, 67} (1), 1-48. doi:10.18637/jss.v067.i01.

Birky, C. W. (2010). Positively Negative Evidence for
Asexuality.\emph{Journal of Heredity} , \emph{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.), \emph{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,\emph{Cobitis taenia} : a model-based approach. \emph{Folia
Biologica} ,\emph{51 Suppl} , 55--60.

Barto\selectlanguage{polish}š, \selectlanguage{english}O., R\selectlanguage{ngerman}öslein, J., Kotusz, J., Paces, J., Pekárik, L., Petrtýl,
M., \ldots{} \& Juchno, D. (2019). The legacy of sexual ancestors in
phenotypic variability, gene expression and homoeolog regulation of
asexual hybrids and polyploids. \emph{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. \emph{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. \emph{Nature Communications} , 8(1), 873. doi:
10.1038/s41467-017-01002-8

Burgarella, C., Gayral, P., Ballenghien, M., Bernard, A., David, P.,
Jarne, P., \ldots{} Glémin, S. (2015). Molecular Evolution of Freshwater
Snails with Contrasting Mating Systems. \emph{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.\emph{Nature} , \emph{348} (6300), 454--455.
doi:10.1038/348454a0

Choleva, L., Janko, K., De Gelas, K., Bohlen, J., \selectlanguage{polish}Š\selectlanguage{english}lechtov\selectlanguage{ngerman}á, V., Rábová,
M., \& Ráb, P. (2012). Synthesis of Clonality and Polyploidy in
Vertebrate Animals by Hybridization Between Two Sexual
Species.\emph{Evolution} , \emph{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 (\emph{Cobitis} ; Teleostei), despite
Clonal Reproduction of Hybrids. \emph{PLoS ONE} , \emph{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.\emph{Bioinformatics, 21} (18), 3674--3676.
https://doi.org/10.1093/bioinformatics/bti610

Dedukh, D., Majtánová, Z., Marta, A., P\selectlanguage{polish}š\selectlanguage{english}eni\selectlanguage{polish}č\selectlanguage{english}ka, M., Kotusz, J., Kl\selectlanguage{ngerman}íma,
J., \ldots{} Janko, K. (2019). Parthenogenesis as a solution to hybrid
sterility: the mechanistic basis of meiotic distortions in clonal and
sterile hybrids. \emph{BioRxiv} , 663112. https://doi.org/10.1101/663112

DePristo, M. A., Banks, E., Poplin, R., Garimella, K. V., Maguire, J.
R., Hartl, C., \ldots{} Daly, M. J. (2011). A framework for variation
discovery and genotyping using next-generation DNA sequencing
data.\emph{Nature Genetics} , \emph{43} (5), 491--498.
doi:10.1038/ng.806

Doncaster, C. P., Pound, G. E., \& Cox, S. J. (2000). The ecological
cost of sex. \emph{Nature} , \emph{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. \emph{Molecular Biology and
Evolution} , \emph{26} (9), 2097--2108. doi:10.1093/molbev/msp119

Fedor\selectlanguage{polish}č\selectlanguage{ngerman}ák, J., Ko\selectlanguage{polish}šč\selectlanguage{english}o, J., Hala\selectlanguage{polish}č\selectlanguage{english}ka, K., \& Manko, P. (2017). Growth
differences in different biotypes of the hybrid complex of \emph{Cobitis
elongatoides} \selectlanguage{ngerman}× \emph{Cobitis tanaitica} (Actinopterygii: Cypriniformes:
Cobitidae) in the Okna River (Danube River basin), Slovakia. \emph{Acta
Ichthyologica Et Piscatoria} , \emph{47} , 125--132.
doi:10.3750/AIEP/02059

Flot, J.-F., Hespeels, B., Li, X., Noel, B., Arkhipova, I., Danchin, E.
G. J., \ldots{} Van Doninck, K. (2013). Genomic evidence for ameiotic
evolution in the bdelloid rotifer \emph{Adineta vaga} . \emph{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. \emph{Evolutionary Biology} , \emph{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., \ldots{} Galtier, N. (2013). Reference-Free Population
Genomics from Next-Generation Transcriptome Data and the
Vertebrate--Invertebrate Gap. \emph{PLoS Genetics} , \emph{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.\emph{Nature} , \emph{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.\emph{Evolution} , \emph{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.\emph{Systematic Biology} , \emph{59} (3), 307--321.
doi:10.1093/sysbio/syq010

Hala\selectlanguage{polish}č\selectlanguage{english}ka, K., Luskov\selectlanguage{ngerman}á, V., \& Lusk, S. (2000). Fecundity of \emph{Cobitis
elongatoides} in the Nová \selectlanguage{polish}Ř\selectlanguage{ngerman}í\selectlanguage{polish}š\selectlanguage{english}e Reservoir. \emph{Folia Zoologica} , (49),
141--150.

Halligan, D. L., \& Keightley, P. D. (2003). How many lethal
alleles?\emph{Trends in Genetics: TIG} , \emph{19} (2), 57--59.
doi:10.1016/S0168-9525(02)00045-8

Hamilton, W. D. (1980). Sex versus Non-Sex versus Parasite.\emph{Oikos}
, \emph{35} (2), 282--290. doi:10.2307/3544435

Hartfield, M., \& Keightley, P. D. (2012). Current hypotheses for the
evolution of sex and recombination. \emph{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 \emph{Timema} Stick Insects. \emph{Molecular
Biology and Evolution} , \emph{29} (1), 401--408.
doi:10.1093/molbev/msr237

Hollister, J. D., Greiner, S., Wang, W., Wang, J., Zhang, Y., Wong, G.
K.-S., \ldots{} Johnson, M. T. J. (2015). Recurrent loss of sex is
associated with accumulation of deleterious mutations in\emph{Oenothera}
. \emph{Molecular Biology and Evolution} , \emph{32} (4), 896--905.
doi:10.1093/molbev/msu345

Howe, D., \& Denver, D. (2008). Muller's Ratchet and compensatory
mutation in \emph{Caenorhabditis briggsae} mitochondrial genome
evolution. \emph{BMC Evolutionary Biology} , \emph{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).\emph{Journal of Experimental Zoology Part A: Comparative
Experimental Biology} , \emph{305A} (6), 513--523. doi:10.1002/jez.a.283

Janko, K., Bohlen, J., Lamatsch, D., Flaj\selectlanguage{polish}š\selectlanguage{english}hans, M., Epplen, J. T., R\selectlanguage{ngerman}áb,
P., \ldots{} \selectlanguage{polish}Š\selectlanguage{english}lechtov\selectlanguage{ngerman}á, V. (2007). The gynogenetic reproduction of
diploid and triploid hybrid spined loaches (\emph{Cobitis} : Teleostei),
and their ability to establish successful clonal lineages---on the
evolution of polyploidy in asexual vertebrates. \emph{Genetica}
,\emph{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. \emph{Biology Direct} , \emph{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. \emph{Evolution}
,\emph{62} (5), 1264--1270.

Janko, K., Flaj\selectlanguage{polish}š\selectlanguage{english}hans, M., Choleva, L., Bohlen, J., \selectlanguage{polish}Š\selectlanguage{english}lechtov\selectlanguage{ngerman}á, V.,
Rábová, M., \ldots{} others. (2007). Diversity of European spined
loaches (genus \emph{Cobitis} L.): an update of the geographic
distribution of the \emph{Cobitis taenia} hybrid complex with a
description of new molecular tools for species and hybrid
determination.\emph{Journal of Fish Biology} , \emph{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. \emph{Evolution} , \emph{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 (\emph{Cobitis} ; Teleostei) using the
analysis of mitochondrial DNA variation. \emph{Molecular Ecology}
,\emph{14} (10), 2991--3004.

Janko, K., Kotusz, J., De Gelas, K., Slechtová, V., Opoldusová, Z.,
Drozd, P., \ldots{} Balá\selectlanguage{polish}ž, \selectlanguage{english}M. (2012). Dynamic formation of asexual
diploid and polyploid lineages: multilocus analysis of
\emph{Cobitis}reveals the mechanisms maintaining the diversity of
clones. \emph{PloS One} , \emph{7} (9), e45384.
doi:10.1371/journal.pone.0045384

Janko, K., Pa\selectlanguage{polish}č\selectlanguage{english}es, J., Wilkinson-Herbots, H., Costa, R. J., Roslein, J.,
Drozd, P., \ldots{} Choleva, L. (2018). Hybrid asexuality as a primary
postzygotic barrier between nascent species: On the interconnection
between asexuality, hybridization and speciation. \emph{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 (\emph{Cobitis} , Cobitidae, Pisces) from
the Black Sea basin that are morphologically similar to \emph{C. taenia}
, with the description of a new species. \emph{Folia Zoologica}
,\emph{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.\emph{Evolution} , \emph{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. \emph{Evolution} , \emph{61} (11),
2728--2735. doi:10.1111/j.1558-5646.2007.00233.x

Juchno, D., Arai, K., Boro\selectlanguage{polish}ń, \selectlanguage{english}A., \& Kujawa, R. (2016). Meiotic
chromosome configurations in oocytes of \emph{Cobitis taenia} and its
polyploid hybrids. \emph{Ichthyological Research} , 1--4.
doi:10.1007/s10228-016-0556-1

Juchno, D., \& Boron, A. (2010). Fecundity of the spined
loach,\emph{Cobitis taenia} (Pisces, Cobitidae) and natural
allopolyploids of\emph{Cobitis} from a diploid-polyploid population.
\emph{Folia Zoologica} . Retrieved from
http://www.highbeam.com/doc/1P3-2013942551.html

Juchno, D., Jab\selectlanguage{polish}ł\selectlanguage{english}o\selectlanguage{polish}ń\selectlanguage{english}ska, O., Boro\selectlanguage{polish}ń, \selectlanguage{english}A., Kujawa, R., Leska, A., Grabowska,
A., \ldots{} Lao, M. (2014). Ploidy-dependent survival of progeny
arising from crosses between natural allotriploid \emph{Cobitis} females
and diploid \emph{C. taenia} males (Pisces, Cobitidae). \emph{Genetica}
,\emph{142} (4), 351--359. doi:10.1007/s10709-014-9779-0

Kondrashov, A. S. (1993). Classification of Hypotheses on the Advantage
of Amphimixis. \emph{Journal of Heredity} , \emph{84} (5), 372--387.

Kotusz, J. (2008). Morphological Relationships Between Polyploid Hybrid
Spined Loaches of the Genus \emph{Cobitis} (Teleostei: Cobitidae) and
their Parental Species. \emph{Annales Zoologici} , \emph{58} , 891--905.
doi:10.3161/000345408X396800

Kotusz, J., Popio\selectlanguage{polish}ł\selectlanguage{english}ek, M., Drozd, P., De Gelas, K., \selectlanguage{polish}Š\selectlanguage{english}lechtov\selectlanguage{ngerman}á, 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 \emph{Cobitis taenia} hybrid
complex.\emph{Biological Journal of the Linnean Society} , \emph{113}
(1), 220--235. doi:10.1111/bij.12329

Lefort, V., Longueville, J.-E., \& Gascuel, O. (2017). SMS: Smart Model
Selection in PhyML. \emph{Molecular Biology and Evolution} ,\emph{34}
(9), 2422--2424. doi:10.1093/molbev/msx149

Li, H., \& Durbin, R. (2009). Fast and accurate short read alignment
with Burrows-Wheeler transform. \emph{Bioinformatics (Oxford, England)}
,\emph{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. \emph{Journal of Botany} , \emph{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
(\emph{Poecilia formosa} ). \emph{BMC Evolutionary Biology} , 8(1), 88.
https://doi.org/10.1186/1471-2148-8-88

Maciak, S., Janko, K., Kotusz, J., Choleva, L., Boro\selectlanguage{polish}ń, \selectlanguage{english}A., Juchno, D.,
\ldots{} Konarzewski, M. (2011). Standard Metabolic Rate (SMR) is
inversely related to erythrocyte and genome size in allopolyploid fish
of the \emph{Cobitis taenia} hybrid complex. \emph{Functional Ecology}
,\emph{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. \emph{Trends in Parasitology} , 32(11), 855--862.
https://doi.org/10.1016/j.pt.2016.08.006

Majt\selectlanguage{ngerman}á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 (\emph{Cobitis} , Cypriniformes,
Teleostei).\emph{PloS One} , \emph{11} (1), e0146872.
doi:10.1371/journal.pone.0146872

McKenna, A., Hanna, M., Banks, E., Sivachenko, A., Cibulskis, K.,
Kernytsky, A., \ldots{} DePristo, M. A. (2010). The Genome Analysis
Toolkit: a MapReduce framework for analyzing next-generation DNA
sequencing data. \emph{Genome Research} , \emph{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. \emph{MBio} , \emph{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. \emph{Molecular Biology and Evolution} , \emph{27}
(4), 954--963. doi:10.1093/molbev/msp300

Otto, S. P. (2007). The evolutionary consequences of
polyploidy.\emph{Cell} , \emph{131} (3), 452--462.
doi:10.1016/j.cell.2007.10.022

Otto, S. P., \& Whitton, J. (2000). Polyploid Incidence and
Evolution.\emph{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 \emph{Daphnia pulex} . \emph{Evolution}
,\emph{59} (4), 800--813.

Paland, S., \& Lynch, M. (2006). Transitions to Asexuality Result in
Excess Amino Acid Substitutions. \emph{Science} , \emph{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.\emph{Nature}
, \emph{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 \emph{Ranunculus auricomus} complex: examining the effects of
hybridization and mutation accumulation. \emph{Molecular Ecology}
,\emph{22} (23), 5908--5921. doi:10.1111/mec.12533

Pertea, G. (2015). fqtrim: v0.9.4 release. \emph{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. \emph{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. \emph{Bioinformatics (Oxford,
England)} , \emph{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.\emph{Journal of Evolutionary Biology} , \emph{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 \emph{Poeciliopsis} (Atheriniformes:
Poeciliidae). \emph{Proceedings of the National Academy of Sciences} ,
89(1), 348--352. https://doi.org/10.1073/PNAS.89.1.348

Ricker, W. E. (2010). \emph{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. \emph{Genetics} , \emph{198} (3), 1059--1069.
doi:10.1534/genetics.114.170472

Schön, I., \& Martens, K. (2003). No slave to sex. \emph{Proceedings of
the Royal Society of London. Series B: Biological Sciences} ,\emph{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.\emph{Molecular Ecology} , \emph{18} (1), 28--42.
doi:10.1111/j.1365-294X.2008.03992.x

Sémon, M., \& Wolfe, K. H. (2007). Consequences of genome
duplication.\emph{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. \emph{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. \emph{Biology Direct} ,\emph{5} , 14.
doi:10.1186/1745-6150-5-14

Stajich, J. E., Block, D., Boulez, K., Brenner, S. E., Chervitz, S. A.,
Dagdigian, C., \ldots{} Birney, E. (2002). The Bioperl toolkit: Perl
modules for the life sciences. \emph{Genome Research} , \emph{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. \emph{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., \ldots{} Holliday, J. A. (2016). Exome capture from
the spruce and pine giga-genomes. \emph{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 \emph{Daphnia pulex} . \emph{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., \ldots{} DePristo, M. A. (2013). From
FastQ data to high confidence variant calls: the Genome Analysis Toolkit
best practices pipeline. \emph{Current Protocols in Bioinformatics}
,\emph{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., \ldots{} Schartl, M. (2018). Clonal polymorphism and high
heterozygosity in the celibate genome of the Amazon molly. \emph{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.), \emph{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. \emph{Journal of Evolutionary Biology} ,\emph{12}
(6), 1003--1012. doi:10.1046/j.1420-9101.1999.00119.x

{[}dataset{]} Ko\selectlanguage{polish}č\selectlanguage{ngerman}í J., Röslein J., Pa\selectlanguage{polish}č\selectlanguage{english}es J., Kotusz J., Hala\selectlanguage{polish}č\selectlanguage{english}ka K.,
Ko\selectlanguage{polish}šč\selectlanguage{english}o J., Fedor\selectlanguage{polish}č\selectlanguage{ngerman}á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\selectlanguage{polish}š, \selectlanguage{english}NI, performed the
experiments. JKo\selectlanguage{polish}č, \selectlanguage{english}JR, JP, JKot, KJ analysed data. JKo\selectlanguage{polish}č, \selectlanguage{english}JKot, JR, and
KJ contributed to writing the MS.

Figures\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/image1/image1}

\end{center}
\end{figure}

\textbf{Fig. 1:} Sampling sites of Cobitis study materials. E, T, N -
genotypes of \emph{C. elongatoides} , \emph{C. taenia} , \emph{C.
tanaitica} , respectively; number of specimens in parentheses. Red area
= range of \emph{C. taenia} ; yellow area = range of \emph{C.
elongatoides} ; blue area = range of \emph{C. tanaitica} . Dashed line
delimitates range of the oldest asexual lineage, the `hybrid clade I'.\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/image2/image2}

\end{center}
\end{figure}

\textbf{Fig. 2:} Maximum Likelihood-based reconstruction of phylogenetic
relationships of mtDNA haplotypes as input for RELAX analysis. Asexual
branches were tested (blue) against intraspecific sexual branches, which
served as reference branches (red) in the RELAX model. Blue box
represents lineages of the oldest asexual lineage, the `Hybrid clade I'.
Letter codes after each sample indicate its genomotype.\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/image3/image3}

\end{center}
\end{figure}

\textbf{Fig. 3:} A) Counts of synonymous (Ns) and nonsynonymous (Nn)
mutations. The left panel shows comparisons of polymorphisms segregating
within (EE, NN, TT) and fixed between (ExN, ExT, TxN) sexual species.
The right panel shows averages (and s.d.) for private asexual SNPs per
each hybrid genomotype. Letters correspond to haploid genomes as
follows: E = \emph{C. elongatoides} , N = \emph{C. tanaitica} , T =
\emph{C. taenia} ; `x' stands for interspecific comparison, so that e.g.
ExT indicates comparison performed on \emph{C. elongatoides} and
\emph{C. taenia} . B\textbf{)} Boxplots of observed heterozygosities
within sexual and asexual individuals.\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/image4/image4}

\end{center}
\end{figure}

\textbf{Fig. 4:} Correlation between accumulation of private mutations
and loss of heterozygosity (LOH) events. Green and blue lines represent
second-order polynomial models fitted to diploid data: g(x) corresponds
to the best fit, while f(x) corresponds to the model with the quadratic
coefficient forced to positive values. Note, that the second model fits
convex parabole, consistent with the assumption that accumulating LOH
events erase increasing proportion of private SNPs. Letters correspond
to haploid genomes as follows: E = \emph{C. elongatoides} , N = \emph{C.
tanaitica} , T = \emph{C. taenia} .\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/image5/image5}

\end{center}
\end{figure}

\textbf{Fig. 5:} Radicality of amino acid substitutions and\emph{d}
\textsubscript{N}/\emph{d} \textsubscript{S} ratios. A) boxplots of
PAM100 score distributions of intraspecific, interspecific and private
asexual SNPs (the latter split for ancient and recent clones). Lower
score means more radical change. B) Boxplots of\emph{d}
\textsubscript{N}/\emph{d} \textsubscript{S} ratio within sexual
species, interspecific pairs, artificial F1 hybrids and hybrid
genomotypes split into ancient and recent clones. The horizontal bars
indicate the distribution medians, while diamonds indicate the means. C)
Histograms on linear and logarithmic (insert) scales of absolute
frequencies of private asexual SNPs within identified clonal lineages. X
axis indicates the number of specimens of the same clone that share
given SNP.\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/image6/6}

\end{center}
\end{figure}

\textbf{Fig. 6:} Comparison of selected fitness related traits between
seasons and ancient vs. recent clones. A) shows the oocyte size
distributions between ancient (red) and recent (green) clones during the
spring (squares) and autumn (diamonds) seasons; B-D) depict Clark's
condition factor, gonadosomatic index (GSI), and absolute fecundities.

Tables

\textbf{Tab. 1: Overview of samples used in genomic part of the study.}\selectlanguage{english}
\begin{longtable}[]{@{}llll@{}}
\toprule
\textbf{genomotype} & \textbf{\# of individuals} & \textbf{category} &
\textbf{Clone ID}\tabularnewline
EE & 8 & sexual & NA\tabularnewline
NN & 4 & sexual & NA\tabularnewline
TT & 10 & sexual & NA\tabularnewline
LL & 1 & outgroup & NA\tabularnewline
EEN & 10 & old clone & 9 (10)\tabularnewline
EET & 5 & young clone & 1 (3), 16 (2)\tabularnewline
ET & 7 & F1, young clone & F1 (2), 2d (1), 10d (2), 12d (1), 67
(1)\tabularnewline
EN & 2 & old clone & NA\tabularnewline
ETT & 4 & young clone & 53 (1), 61 (1), 65 (1)\tabularnewline
\bottomrule
\end{longtable}

Upper case letters represent a haploid genome of these species: E
=\emph{C. elongatoides} ; L = \emph{C. paludica} ; N = \emph{C.
tanaitica} ; T = \emph{C. taenia} ; for example EET stands for triploid
with two genomes of \emph{C. elongatoides} and one of \emph{C. taenia} .
Categories were used to divide datasets for use in linear models
(LME/GLM/GLME), where sexuals were also used for interspecific
differences between pairs of species. F1 denotes F1 laboratory hybrids.
Clone ID denote the membership of given individual to particular
multilocus lineage inferred from SNP data, which represent independent
clonal lineage, with number of individuals per lineage in parentheses.

\textbf{Tab. 2: Locality information and number of
\emph{Cobitis}specimens used in the study concerning fitness proxy
traits.}\selectlanguage{english}
\begin{longtable}[]{@{}llllllll@{}}
\toprule
\textbf{River} & \textbf{Site} & \textbf{Site} & \textbf{Site} &
\textbf{No. spec. old clone} & \textbf{No. spec. old clone} &
\textbf{No. spec. young clones} & \textbf{collection
dates}\tabularnewline
& \textbf{Latitude} & \textbf{Longitude} & \textbf{(spring/ fall)} &
\textbf{(spring/ fall)} & \textbf{(spring/ fall)} & \textbf{(spring/
fall)} &\tabularnewline
Skora & 51\selectlanguage{ngerman}° 17' & 16\selectlanguage{ngerman}° 05' & (3/5) & (3/5) & (4/20) & (4/20) &
21.04.2005, 12.09.2013\tabularnewline
Z\selectlanguage{polish}ł\selectlanguage{english}otnica & 51\selectlanguage{ngerman}° 33' & 17\selectlanguage{ngerman}° 41' & 0/12 & 0/12 & 0/12 & 0/12 &
11.09.2013\tabularnewline
Polska Woda & 51\selectlanguage{ngerman}° 31' & 17\selectlanguage{ngerman}° 30' & 0/23 & 0/23 & 0/4 & 0/4 &
12.09.2013\tabularnewline
Sw\selectlanguage{polish}ę\selectlanguage{english}drnia+ & 51\selectlanguage{ngerman}° 49' & 18\selectlanguage{ngerman}° 15' & 0 & 0 & 0/4 & 0/4 &
11.09.2013\tabularnewline
Budkowiczanka & 50\selectlanguage{ngerman}° 52 ' & 18\selectlanguage{ngerman}° 02' & (1/0) & (1/0) & (2/0) & (2/0) &
22.04.2005\tabularnewline
Wierzbiak & 51\selectlanguage{ngerman}° 14' & 16\selectlanguage{ngerman}° 14' & (2/0) & (2/0) & (2/0) & (2/0) &
21.04.2005\tabularnewline
Sumina & 50\selectlanguage{ngerman}° 09' & 18\selectlanguage{ngerman}° 24' & (6/0) & (6/0) & 0 & 0 &
22.04.2005\tabularnewline
total & & & 12/40 & 12/40 & 10/40 & 10/40 &\tabularnewline
\bottomrule
\end{longtable}

+the sample where the parental species was \emph{C. taenia} , on the
remaining sites parental was \emph{C. elongatoides}

\selectlanguage{english}
\FloatBarrier
\end{document}

