Hybrid asexuality as a primary reproductive barrier: on the interconnection between asexuality and speciation

Speciation usually proceeds in a continuum from intensively hybridizing populations until the formation of irreversibly isolated species. Restriction of interspecific gene flow may often be achieved by gradual accumulation of intrinsic postzygotic incompatibilities with hybrid infertility typically evolving more rapidly than inviability. A reconstructed history of speciation in European loaches (Cobitis) reveals that accumulation of postzygotic reproductive incompatibilities may take an alternative, in the literature largely neglected, pathway through initiation of hybrids’ asexuality rather than through a decrease in hybrids’ fitness. Combined evidence shows that contemporary Cobitis species readily hybridize in hybrid zones, but their gene pools are isolated as hybridization produces infertile males and fertile but clonally reproducing females that cannot mediate introgressions. Nevertheless, coalescent analyses indicated intensive historical gene flow during earlier stages of Cobitis diversification, suggesting that non-clonal hybrids must have existed in the past. The revealed patterns imply that during the initial stages of speciation, hybridization between little diverged species produced recombinant hybrids mediating gene flow, but growing divergence among species caused disrupted meiosis in hybrids resulting in their clonality, which acts as a barrier to gene flow. Comparative analysis of published data on other fish hybrids corroborated the generality of our findings; the species pairs producing asexual hybrids were more genetically diverged than those pairs producing fertile sexual hybrids but less diverged than species pairs producing infertile hybrids. Hybrid asexuality therefore appears to evolve at lower divergence than other types of postzygotic barriers and might thus represent a primary reproductive barrier in many taxa.


Introduction
Speciation may be accomplished by various reproductive isolating mechanisms (RIMs) that restrict interspecific gene flow. These RIMs can be initiated either by divergent ecological or sexual selection that creates prezygotic and extrinsic postzygotic reproductive isolation or by gradual accumulation of genetic incompatibilities that cause intrinsic postzygotic isolation, i.e. hybrid sterility and inviability (rev. in [1]). The study of intrinsic postzygotic isolation has a strong tradition in evolutionary biology [2,3] and comparative studies showed that hybridization capability decreases as one moves from closely to distantly related pairs of taxa [4][5][6][7]. Although the rate at which intrinsic postzygotic RIMs accumulate is probably nonlinear [8][9][10] and varies among taxa [11], hybrid infertility generally evolves at lower genetic distances than hybrid inviability [5,12] and according to Haldane's rule, both hybrid sterility and inviability evolve faster in the heterogametic sex [13].
Here we suggest that speciation driven by intrinsic postzygotic barriers may take an alternative, in literature largely neglected, pathway: the accumulation of genetic incompatibilities may lead to distortion of hybrids' reproductive mode towards asexuality already at early stages of the species diversification when hybrids' fitness has not yet considerably decreased. Given that obligate asexuals are generally incapable to introgress their genomes into related sexual species (e.g. [14,15]), such modification of hybrids' reproductive mode may effectively block interspecific gene exchange and contribute to speciation even if other forms of RIMs are absent. This hypothesis is consistent with the observation that asexual organisms often appear as hybrids between sexually reproducing species [16]. Various mutually not exclusive theories attempted to explain the association between hybridization and asexuality (rev. in [17]). Among those, it has been proposed early in 20 th century that the likelihood of initiation of hybrid asexuality could depend on the genetic distance between hybridizing species [18], possibly because the coexistence of distinct genomes within progenitor egg cells may deregulate the pathways underlying meiosis [19][20][21].
Moritz et al. [22] formulated the "balance hypothesis" as a candidate general model of the hybridization-asexuality association. It postulates that asexuality can arise only when the genomes of parental species accumulated enough mutational incompatibilities to alter meiosis or gametogenesis in hybrids yet not enough to seriously compromise hybrid viability or fertility.
From this perspective, one could regard hybrid asexuality as a special case of the accumulation of Dobzhansky-Muller incompatibilities that disrupt critical processes -sexual reproduction in this case -and thus directly contributes to speciation. Unfortunately, there is only scarce empirical support for the link between accumulating divergences among hybridizing genomes and the initiation of asexuality [23] and hence, it remains unclear how or even whether hybridization could initiate the production of clonal gametes [17,24,25].
This study explores the role of hybrid asexuality in speciation by simultaneous reconstruction of the rates of interspecific gene flow and evolution of asexuality during the diversification history of European spined loaches (Cobitis). These bottom-dwelling fish colonized Europe in several phylogenetic lineages during the Tertiary period. While isolated rivers of Mediterranean Europe are inhabited by diverse Cobitis species clustering in deeply divergent phylogenetic lineages (Cobitis Lineage I -V sensu [26]), vast areas of non-Mediterranean Europe were colonized by a single lineage (Cobitis Lineage V) comprising several morphologically and ecologically very similar species. Among these, four closely related species (C. tanaitica, C. taenia, C. taurica and C. pontica) have Ponto Caspian distribution (C. taenia further colonised Northern and Western Europe) while their distant relative, C. elongatoides is distributed throughout the Danubian basin. Their ranges have been fluctuating over the Quaternary [27] and currently display a parapatric distribution with overlapping zones located in Central Europe, the lower Danube, and Southern Ukraine [28] (Fig 1). Hybridization has been documented to take place in these zones and reproductive contact between C. elongatoides and either C. taenia, C. tanaitica or C. pontica has been documented to result in clonally reproducing all-females hybrid lineages [24,28,29]. These hybrid lineages achieved remarkable evolutionary success and colonised most of European continent, some of them having achieved considerable ages (i.e. the so called Hybrid Clades I and II are as ancient as 0.35 and 0.25 Mya, respectively [27,30]). Clonal reproduction of these hybrid lineages theoretically prevents any interspecific gene flow but we have recently documented evidence for intensive gene flow between C. elongatoides and C. tanaitica leading to the fixation of an introgressed C. elongatoides-like mitochondrion over the entire C. tanaitica distribution range [31]. Locality numbers correspond with Table 2. To understand the speciation history in spined loaches and more generally the role of hybrid asexuality in fish speciation, we employed multiple complementary approaches examining a) whether the establishment of reproductive barriers can be primarily accomplished by the formation of hybrid asexuality, b) whether the formation of asexual hybrids does depend on genetic distance between parental species. First, we analyzed the reproductive modes in hybrids between several differently related Cobitis species to reveal their asexuality and the extent of the currently observed reproductive isolation. Second, we performed a detailed population genetic analysis of hybrid zones to test if there is any ongoing introgressive hybridization. Third, we employed coalescent analyses of transcriptomic data to estimate levels and timing of historical gene flow among species.
The combined evidence implies that the diversification of Cobitis species has been accompanied by decreasing intensity of introgressive hybridization. However such a restriction in gene flow has not been accomplished by classical RIMs but rather by the production of asexual hybrids suggesting that asexual hybrids may constitute a primary reproductive barrier between nascent species.

Results
To test whether the formation of asexual hybrids does depend on genetic distance between parental species and how the formation of hybrid asexuality affects the levels of reproductive isolation we performed the following four analyses.

Analysis of reproductive modes of interspecific Cobitis hybrids
Firstly, we tested whether Cobitis hybrids can mediate interspecific gene flow by analysing reproductive modes of natural and artificial hybrids between several differently related Cobitis species.
1. Reproductive modes of natural C. elongatoides-tanaitica hybrids: We successfully analysed five backcrossed families stemming from four natural triploid EEN and one diploid EN hybrid females (letters E and N stand for haploid elongatoides and tanaitica genomes, respectively). Backcross progeny consistently expressed all maternal alleles, suggesting the production of unreduced gametes and lack of segregation. In one allele, single progeny differed from the maternal allele by a single repeat, indicating a mutation event. A part of the progeny also contained a haploid set of paternal alleles indicating that the sperm's genome is sometimes incorporated, leading to a ploidy increase (Table 1 and S1). We also compared allozyme profiles of eggs and somatic tissues of six EEN females, in order to test for hemiclonal reproduction (i.e., hybridogenesis). Hybridogenesis was supposed to lead to reduced, albeit non-segregating gametes [32][33][34], which consistently express allozyme products of one parental species only. Nevertheless, we found that allozyme profiles of all eggs were identical to the somatic tissues of their maternal individuals, further suggesting that natural hybrids do not exclude any genome (Table S2). progeny: 15 individuals contained the complete set of maternal alleles and 9 of them also possessed the haploid set of paternal alleles (Table 1). Such patterns indicate that C. taenia-pontica hybrid females produced not only recombinant sexual gametes but partly also unreduced gametes.
The unreduced eggs either developed into a clonal progeny via gynogenesis or into triploid progeny after a true fertilization with a haploid sperm.

Analysis of C. elongatoides -C. tanaitica hybrid zone and test of ongoing introgressive hybridization
Second, we tested for ongoing gene flow between species by performing population genetic analysis of C. elongatoides -C. tanaitica hybrid zone similar to C. elongatoides -C. taenia analysed in [30]. 818 Cobitis specimens were captured and identified on 26 localities all over the lower Danubian River Basin. Due to morphological similarity, the taxonomical diagnosis was performed by allozyme analysis of diagnostic loci [28], which revealed the presence of C. elongatoides at 11 sites and C. tanaitica at seven sites. 10 localities were inhabited by distantly related species, C.
strumicae, belonging to the subgenus Bicanestrinia (lineage III sensu [26]) (Fig 1 and Table 2). Apart from sexual species, we found various hybrid biotypes, which were mostly polyploid but we also encountered 32 diploid hybrids heterozygotic for species-specific alleles at all diagnostic allozyme loci ( Table 2).  (Table S3), but 28 diploids clustered in eight groups of identical MultiLocus Genotypes (MLG A -MLG H). Such groups of identical individuals were considered as clone mates due to negligible probability of identical genotypes arising from independent sexual events (p < 10 -5 ). Furthermore, we noticed that two clusters of individuals or MLGs were separated by genetic distances, which were significantly lower than would be expected from independent sexual events (p < 0.01). As in [30], both such groups were assigned as MultiLocus Lineages (MLL) and individuals clustering within such lineages were considered as members of the same clone.
Altogether, 29 diploid individuals could be assumed to form eight clonal lineages. Allozymes indicated a hybrid state of all such clonal individuals. No mtDNA haplotypes were shared by C.
elongatoides and C. tanaitica (Fig 2). One clone (MLG B) possessed a haplotype E1, which was shared with C. elongatoides, while the remaining individuals assigned as hybrids possessed haplotypes clustering in the old Hybrid clade I defined in [27].   Fig 3A). Data analysis using the NewHybrids software was slightly sensitive to the type of prior applied to Ɵ but not to π. However, consistent with the Structure, all specimens presumed to be C. elongatoides and C. tanaitica were always assigned as the pure parental species (p > 0.95). Under Jeffrey's prior, the NewHybrids assigned all but one of the abovementioned hybrids into the F1 class with the probability exceeding 95 %; however, it could not decide in the case of MLL E between F1 (p = 0.53), F2 (p = 0.10), and B1 (p = 0.35) states ( Fig 3B). Under the uniform prior, most hybrids were assigned into F1 class, but NewHybrids could not decide between F1 and B1 states of MLL H, MLL E and individual 09BG19K22 ( Fig 3B and Table S3).

Estimation of levels and timing of historical gene flow among the species using transcriptome data
Thirdly, we estimated levels and timing of historical gene flow among the species by analysed single nucleotide polymorphisms (SNP) variability of C. elongatoides, C. taenia, C. tanaitica and C. pontica transcriptomes. As an outgroup we used C. strumicae -a species belonging to subgenus Bicanestrinia whose divergence from the ingroup was set to 17.7 Mya according to [35], which was sampled in isolated Black Sea tributaries outside of the distribution ranges of ingroup species in order to avoid any possible reproductive contact. The assembly of mRNA from five C. taenia specimens comprised 20,385 contigs (potential mRNAs) and was used as a reference for mapping of reads from all species (two individuals of each ingroup species and one outgroup). In total, we identified 187,205 SNPs (Appendix S4), which were then analysed by fitting eight coalescent models assuming different scenarios of species divergence and connectivity [36][37][38][39][40]. The models are characterized in Fig 4 and the results of analyses are described in Table 3. Because the available models only allow fitting pairs of taxa, we separately analysed six datasets corresponding to all pairwise combinations of the four analysed species. Most models were successfully fitted to all datasets but two models did not converge in the case of the C. taenia -C. tanaitica dataset.    This model fitted all pairwise species datasets significantly better than both other models of Wilkinson-Herbots [36,37], namely the isolation model I4 assuming no gene flow since the initial split of the species, and the isolation-with-migration model IM4 assuming ongoing hybridization from the initial split until the present (Likelihood Ratio Test; LRT p << 10 -10 for all comparisons involving C. elongatoides, and p < 0.01 for all other species comparisons). To test if the better fit of IIM7 is due to the correct assessment that interspecific gene flow ceased recently, or merely due to allowing an additional size change at time t 1 , we employed five more models that allow an additional population size change and also relax the assumption of equal population sizes during the migration stage (Fig 4; [ 39,40]).
Pairwise datasets of the three closely related C. taenia, C. tanaitica and C. pontica were better fitted by the isolation-with-initial-migration models (IIM7, IIM8) compared to the I4 and I7 isolation models (LRT p < 0.01 for all comparisons). The generalized isolation-with-migration GIM9 model gave estimates of 0 for the current migration rate and hence reduced to the IIM8 model (where these models could be fitted). However, these models consistently suggested that the The datasets of species pairs including C. elongatoides and any one of C. taenia, C. tanaitica or C. pontica were fitted better by the isolation-with-initial-migration model IIM8 than by isolation models I7 and I4 (LRT p < 0.0002 in all cases) and both isolation-with-migration models, IM4 and IM5 (LRT p << 10 -10 in all cases). The IIM8 model also fitted the data better that the IIM7 model (which assumes equal population sizes during the migration stage) in the case of C. elongatoides -C. taenia and C. elongatoides -C.tanaitica species pairs (lower AIC scores and LRT p < 0.023 and 0.001, respectively) but not for C. elongatoides -C. pontica species pairs where IIM7 had a better AIC score.
The divergence time estimates were consistent across species comparisons and place the initial split of C. elongatoides from the other species at roughly 9 Mya. It is estimated that C. elongatoides exchanged genes with the other species at average rates of less than 0.1 migrants per generation until time t 1 , for which ML estimates vary between 1.65 and 0.49 Mya depending on the data set and model but have generally quite large confidence intervals ( Table 3). Given that the speciation times of the other three species are much more recent (see above) than their split from C. elongatoides, our results suggest that the detected gene flow occurred predominantly between C. elongatoides and the common ancestor of the other three species. The most complex GIM9 model did not improve the fit to the data compared to the IIM8 model (as indicated by the AIC scores of the models and by LRT; p > 0.59 for all cases). Moreover, the GIM9 model estimated that recent migration rates were much lower than historical ones, thus virtually converging to the IIM8 model and confirming that C. elongatoides has been historically exchanging genes with the other species but became isolated in recent times.

Comparative analysis of genetic divergence between hybridizing fish species and types of reproductive isolation including hybrid asexuality
Finally, we tested the generality of our findings, implying that the hybrids' asexuality may represent an intermediate stage of species diversification process. To do so, we investigated the general correlation between genetic divergence of hybridizing pairs of fish species and the dysfunction of their F1 fish hybrids. Dysfunction was measured in terms of infertility, inviability or asexual reproduction of hybrids (see Methods section for details how we amended the dataset from Russell's [5] comparative study to incorporate asexuality). We found that the test, p < 0.01) and always significantly higher than that of Russell's [5] hybrid class 0 (t-test, p < 0.05).
Interestingly, the range of divergences among asexual hybrids is notably similar to the divergences of hybrids where the functionality of one sex is lower than that of the other (hybrid classes 0.5 -1.5; mean = 0.118; s.e. = 0.040).

Discussion
Natural asexual organisms are considered as treasured evolutionary models because they could provide clues about the paradox of sex. Our study offers a conceptually novel view on this phenomenon and suggests that hybrid asexuality may form an inherent stage of the speciation process and constitute an effective reproductive barrier, which arises earlier in speciation than hybrid sterility and inviability. Our conclusion is based on the following findings: (1) We demonstrated that both the intensity of interspecific gene flow as well as the likelihood of formation of asexual hybrids change with increasing genetic distance between parental species. (2) We confirmed that hybrid asexuality creates an effective reproductive barrier to interspecific gene flow.
(3) We demonstrated the general applicability of our hypothesis by finding that hybrid asexuality in fish arises at lower genetic distances than complete hybrid sterility and hybrid inviability and might thus represent a primary reproductive barrier in many taxa.

Accomplishment of speciation in spite of ongoing hybridization
The Cobitis species considered readily hybridize both in natural and experimental conditions [24,28,30], suggesting the absence of strong prezygotic RIMs. However, two independent lines of evidence suggest that contemporary interspecific gene flow between the species with overlapping ranges, i.e., C. elongatoides on one side and C. taenia and C. tanaitica on the other side, is unlikely and that speciation has been completed. First, according to experimental crossings of C.
elongatoides-taenia [24,29] and C. elongatoides-tanaitica (present study), hybrid females do not produce reduced gametes through the "standard" sexual process and do not produce recombinant progeny. Allozyme analyses of oocytes further argue against the production of reduced gametes through hybridogenesis; contrary to evidence from hybridogenetic vertebrates, which exclude one parental species' genome [32][33][34], we found that oocytes of studied Cobitis hybrid females expressed allozyme alleles of both parental species. Unreduced gametes either develop clonally or occasionally incorporate the sperm's genome leading to polyploidy, but are unlikely to enable interspecific gene flow. Hybrid males also do not appear to mediate gene flow since C. elongatoides-taenia hybrid males are infertile [24] and C. elongatoides-tanaitica hybrid males have not been observed in nature.
Second, the completion of speciation is also supported by analyses of C. elongatoides -C.
taenia and C. elongatoides -C. tanaitica hybrid zones. The contact zone between C. elongatoides and C. tanaitica, covering most of the Danube watershed below the Iron Gates Gorge, is larger than the Central European C. elongatoides -taenia hybrid zone [30] and also involves the presence of a third, phylogenetically distant, species C. strumicae, which invaded Danubian stretches during postglacial expansion from the Mediterranean Basin [26]. Nevertheless, both zones have similarities with respect to major evolutionary aspects. Polyploid asexuals dominate in both zones (in the Danube drainage, they also invaded rivers inhabited by C. strumicae [42] and diploids within both zones could consistently be assigned into two types. Mya. Contemporary hybrids of C. elongatoides thus do not appear able to produce recombinant progeny mediating interspecific gene flow. As Cobitis species do not notably differ in ecological requirements, thereby excluding the possibility of strong extrinsic postzygotic isolation, we suggest that the asexuality of hybrid females together with the sterility of hybrid males is the primary reproductive barrier separating these species. Although we found no evidence of contemporary gene flow between C. elongatoides on one side and C. taenia and C. tanaitica on the other side, coalescent analysis of RNAseq data suggested that hybrids that were able to mediate gene flow must have existed in the past. We clearly rejected isolation models, suggesting that gene flow existed after the initial divergence of C. elongatoides from C. taenia, C. pontica and C. tanaitica. Interestingly, datasets including C. elongatoides were significantly better fitted by isolation-with-initial-migration models than by the simpler isolationwith-migration models, which implies that gene flow involving C. elongatoides occurred in the historical period after its initial divergence around 9 Mya but did not occur after t 1 , suggesting the isolation of contemporary populations (Fig 4 and Table 3). As the applied coalescent models assumed gene flow since the initial split of C. elongatoides but compared only pairs of extant species, it is reasonable to conclude that the inferred gene flow occurred between C. elongatoides and the common C. taenia -C. tanaitica -C. pontica ancestor that diversified only recently around ca 0.5 Mya. Unfortunately, the large confidence intervals of the t 1 estimates prevent an unambiguous conclusion as to whether nuclear gene flow continued after the C. taenia -C. tanaitica -C. pontica speciation. Conversely, we found no evidence of gene flow between the closely related C. taenia, C. tanaitica and C. pontica, where coalescent models suggested their recent speciation, probably during the last 0.5 Mya The present results are generally consistent with our previous study based on nine nuclear and one mtDNA loci but two differences were noted. First, Choleva et al. [31] found significant traces of mitochondrial but not nuclear C. elongatoides -C. tanaitica gene flow, which led to the somewhat paradoxical impression that the nucleus has not been affected by hybridization while C. tanaitica's mitochondrion has suffered one of the most massive introgressions among animals (a complete mtDNA replacement). The current analysis of the extended dataset indicated gene flow in the nuclear compartment too, which is more plausible biologically. Second, the previous model-based inference indicated C. taenia -C. tanaitica gene flow in nucleus, which was at odds with the available field data because no C. taenia -tanaitica hybrids have ever been observed in nature and the ranges of both species are closely adjacent but sympatry has not been documented [28]. The present analysis is therefore more in line with current knowledge about Cobitis and the discrepancy with the previous analysis of only few loci may potentially reflect the tendency of Bayesian IM algorithms to inflate estimates of gene flow when the number of loci is low and splitting times are recent [44,45].
Of course, any model-based inferences should be treated with caution. For example, having analysed only pairs of species, our models may be affected by intractable interactions with other species. IM models are robust to this type of violation, for low or moderate levels of gene flow [46], and the inferred lack of hybridization between C. taenia -C. tanaitica -C. pontica implies that the pairwise comparison of any one of these species with C. elongatoides is unlikely to be affected by the existence of the other species. However, other problems may be important, such as uncertainty regarding the relative mutation rates of the different loci (which have been estimated by comparing them with the outgroup), additional population size changes, geographical structure or varying rather than constant intensity of gene flow between t 0 and t 1 . Furthermore, coalescent models assume evolution under neutrality or negative selection, which effectively decreases the locusspecific mutation rate [47], but other types of selection might have affected our data. It is therefore a great advantage that model-based inferences are supported by independent types of data. The existence of gene flow between C. elongatoides and other species is supported by fixation of an C.
elongatoides-like mitochondrial lineage in C. tanaitica [31]. On the other hand, small but non-zero divergences between contemporary C. elongatoides and C. tanaitica mtDNA haplotypes, clonal reproduction of the hybrids studied, and the apparent lack of introgressive hybridization in hybrid zones imply that such gene flow ceased and does not occur at present.

Simultaneous evolution of asexuality and RIMs
Several complementary approaches showed that the initial divergence of C. elongatoides was coupled with the production of hybrids whose reproductive mode enabled more or less intensive gene flow with other three species or common ancestor these species. However, as C. elongatoides diverged from these species, the introgressions became restricted since the major type of hybrids became asexual, which may not effectively mediate the gene flow. Such patterns conform to the "balance hypothesis" [22], which predicts that hybridization between gradually diverging species would initially produce mostly sexual hybrids while successful asexuals would arise at intermediate stages when a hybrid's meiosis is disrupted but fertility is not yet significantly reduced.
Simultaneously, this scenario is consistent with the gradual decline in the species' capability of introgressive hybridization that is expected to evolve along the speciation continuum from weakly separated entities towards strongly isolated species (e.g. [1]). The establishment of asexual hybrids and the gradual accumulation of intrinsic postzygotic RIMs thus appear to be interconnected processes that correlate with divergence of hybridizing species.
Crossing experiments support this idea since distantly related C. elongatoides and C. taenia produced only clonally reproducing F1 females and infertile F1 males [24], while mating of the closely related C. taenia and C. pontica species (this study) resulted in fertile F1 progeny of both sexes that mostly produced recombinant progeny. Such data also conform to the empirical observation that one sex often acquires infertility earlier than the other one [5,9]. Unfortunately, lack of data on Cobitis sex chromosomes prevent speculations as to whether infertility of hybrid males could be related to Haldane's rule or other processes such as the increased sensitivity of spermatogenesis compared to oogenesis to perturbations in hybrids [48]. Interestingly, we noted that some C. taeniapontica female hybrids reproduced clonally. This observation extends the list of previously published asexual lineages and suggests that the emergence of asexual hybrids is a taxonomically widespread phenomenon within the entire family Cobitidae, involving several independent species pairs. These include C. taenia-pontica (laboratory hybrids from this study), C. elongatoides-taenia and C. elongatoides-tanaitica [28] as well as two cases from other related genera of the family [49,50]. Such a widespread and independent emergence of asexuality is not consistent with the phylogenetic constraint hypothesis [51,52], which assumes that the initiation of hybrid asexuals may result from predisposition of only certain parental species or populations to produce such asexuals. Instead, our data are consistent with the Balance Hypothesis and further indicate that hybrid asexuality constitutes an intermediate stage of speciation continuum that may be viewed as a special case of Dobzhansky-Muller incompatibilities. In this respect, it is noteworthy that hybrid asexuality is commonly observed only in one direction of cross [53,54], which is similar to classical patterns observed in hybrid infertility and inviability.
Performed comparative analysis further suggests that the scenario revealed in spined loaches can be extrapolated to other taxa. The inclusion of known asexual hybrids into Russell's [5] study of fish hybrids showed that asexual hybrids appear at intermediate levels of parental divergences between species pairs producing fertile and viable hybrids and those producing infertile hybrids of both sexes ( Fig 5, Appendix S5). In fact, the divergence of asexual-producing species pairs is similar to those producing hybrids with lowered fertility of hybrids of one sex (Russell's hybrid classes 0.5 -1.5), which is in line with the observation that populations of asexual fish often consist of females only.
While Russell [5] corroborated that the severity of a hybrid's dysfunction correlates with the divergence of the hybridizing fish, our amendment suggest that hybrids with a skewed reproductive mode towards asexuality tend to appear at an intermediate stage of the diversification process. The relatively wide range of divergence for the asexual hybrids also agrees with the expectation that accumulation of RIMs follows a variable-rate clock and that different types of incompatibilities accumulate in a noisy manner [11]. Similarly in reptiles, parental species known to produce asexual hybrids typically appear to be somewhat genetically distant rather than being each other's closest relatives [55][56][57], and, divergences among hybridizing species that produce asexual hybrids generally appeared significantly higher than among those producing bisexual hybrids [58].
Interestingly, there also appears to be a continuum in the ability of species pairs to produce hybrid asexuals. At one end, there are dynamically hybridizing species pairs that produce diverse assemblages of asexual hybrids, e.g. Cobitis [24], Pelophylax [51], Poeciliopsis [19]; whereas at the other end there are monoclonal hybrid lineages stemming from one or a few ancient events while attempts to cross their contemporary sexual relatives do not produce clones, e.g., Poecilia [59].
Phoxinus may represent an intermediate case, where historical hybridizations produced a highly diverse asexual assemblage but new clones can no longer be produced [60]. Similarly, hybridization of some species pairs has been documented to produce only asexual hybrids (e.g. Cobitis elongatoides-taenia) while clonal and sexual hybrid females co-occur and fertile hybrid males may exist in other hybridizing pairs (e.g. Fundulus, Rutilus rutilus x Abramis bramma). We suggest that such a variety of patterns may reflect that different species pairs proceeded different distances along the same continuum.
To conclude, we found that the speciation continuum may contain an inherent, previously unnoticed, stage of asexuality. The production of asexual rather than sexual hybrids helps to establish an effective barrier to gene flow even in the absence of other typical forms of RIMs.
Asexuality is phylogenetically widespread but rare phenomenon in animals. However, given that asexual lineages typically have a short life span [61] and that the phase when diverging species could produce asexual hybrids is transient, it is possible that many currently reproductively incompatible species might have historically produced asexual hybrids that have gone extinct. The stage of hybrid asexuality thus certainly requires close attention in research since it might have represented an important mechanism in the speciation of many groups capable of creating asexuals, including arthropods, vertebrates or plants.

Materials and Methods
All experimental procedures involving fish were approved by the Institutional Animal Care and Use  Table 2 indicate the origin of investigated fish specimens used in this study. A priori classification of captured specimens into taxonomic units (species and hybrid types) was based on previously verified diagnostic allozyme loci and PCR-RFLP method or sequencing of the first intron of the nuclear S7 gene following [28]. Ploidy was routinely estimated from gene dose effect using allozymes as described in [28] and also by flow cytometry in cases where we disposed sufficient amount of fixed tissue.
818 individuals were genotyped and those identified as diploids were subsequently subject to analysis of nine microsatellite following the protocols of [24,62]; multiplex 1: loci Cota_006, Cota_010, Cota_027, Cota_068, Cota_093, Cota_111; multiplex 2: loci Cota_032, Cota_033, Cota_041. We also amplified and sequenced 1190 bp fragment of the cytochrome b gene according to [28] in a subset of diploid individuals. The GenAlEx 6.5 software [63] was used to identify clusters of individuals sharing the same microsatellite multilocus genotype (MLG) and to calculate the probability, conditional on observed variability studied loci, that individuals sharing the same MLG arose by independent sexual events (see [30] for details). While some MLG may represent distinct clones, others may belong to the same clonal lineage (the so called 'multilocus lineage' (MLL) and differ only by scoring errors or post formational mutations. We used the approach described in [30] to identify groups of MLGs forming such mutilocus lineages. Phylogenetic relationships among newly discovered and previously published mtDNA haplotypes [30,42] were reconstructed with the medianjoining network [64], and drawn using NETWORK software (http://www.fluxusengineering.com/netwinfo.htm).
2.2 Detection of interspecific gene flow: microsatellite and allozyme data were used to detect admixed individuals from the C. elongatoides -C. tanaitica hybrid zone by two Bayesian clustering methods that are based on different algorithms but have similar power for hybrid detection [65]. We first used admixture model implemented in Structure 2.3.3 [66] to compute the parameter q, i.e. the proportion of an individual's genome originating in one of the two inferred clusters, corresponding to the parental species C. elongatoides and C. tanaitica. The analysis was based on runs with 10 6 iterations, following a burn-in period of 5*10 4 iterations. Ten independent runs for number of populations varying from K= 1 to K= 10. The best value of K was chosen following the method proposed by Evanno et al. [67], with Structure Harvester [68], which takes into account the rate of change in the log likelihood between successive K values.
Second, the Bayesian clustering method implemented in NewHybrids 1.1 [69] was used to compute the posterior probability that an individual in a sample belongs to one of the defined hybrid classes or parental species. The eight genotypic classes investigated were: C. elongatoides, C.
The two backcross types included those having 75 % of their genome originated from the backcrossing species (B1 generation) and those having 95 % of their genome originated from the backcrossing species (further-generation backcrosses). Posterior distributions were evaluated by running five independent analyses to confirm convergence. We started with different random seeds, performed 10 4 burn-in iterations and followed by 500,000 Monte Carlo Markov Chain iterations without using prior allele frequency information. Analyses were run for four combinations of prior distributions (uniform or Jeffreys for Ɵ and π parameters) to explore the robustness of the results [69].
To minimize the effect of clonal propagation, we used only one randomly chosen representative of each unique MLL and repeated the analysis several times to check its robustness against the particular choice of MLL representatives [30]. We also repeated those analyses with either Cota_006 or Cota_041 locus removed due to their possible linkage [62]. Locus Cota_027 was removed since it does not amplify in C. elongatoides.
3. Experimental analysis of reproductive modes of hybrids: We performed crossing experiments under conditions described in [24] and genotyped the parents and their progeny using the microsatellite multiplex 1 to verify their mode or heredity. Two types of crossing were performed: 3.1 First, we analyzed the reproductive mode of naturally occurring diploid and triploid hybrids by individually crossing wild caught EN and EEN females (see Table 1) with males of the sexual species.
Several authors [32][33][34] also reported more complex reproductive modes when one parental genome is excluded prior to meiosis while the other is clonally transmitted. This type of so-called hybridogenetic reproduction leads to reduced, albeit non-segregating gametes and may be detected by comparison of allozyme profiles of somatic and germinal tissues. In order to test for premeiotic genome exclusion in Cobitis hybrids, we have therefore compared allozyme profiles of somatic tissues of hybrid females with those from their ova using the standard allozyme methods adapted for Cobitis [28].  Table 1 and S1).
4. Detection of historical admixture: In order to test for periods of historical or contemporary gene flow on a genome-wide scale, we applied coalescent models to SNP data from mRNA. Initial 1886536 reads (648620753 bp) from C. taenia were filtered to remove low quality part of reads, adaptors and primer sequences using Trimmomatic software [70]. Technical PCR multiplicates were removed using cdhit-454 software [71].  The most complex "generalized isolation-with-migration" model with nine parameters (GIM9) assumes that an ancestral population of size Ɵ a split at time t 0 into two descendant populations of sizes Ɵ 1 and Ɵ 2 interconnected by gene flow at rate M until t 1 , from which time onwards both species are at their current sizes Ɵ c1 and Ɵ c2 and gene flow occurs at its current rate M c . Note that all models described above are nested within the GIM9-model.
All models were ranked according to their AIC score and we also calculated the evidence ratio for each model, providing a relative measure of how much less likely a given model is compared to the best-fitting model, given the set of candidate models considered and the data [69]; Table 3. In addition, Likelihood Ratio tests were performed to compare pairs of nested models, where we have assumed that the use of the χ2 distribution with the appropriate number of degrees of freedom is conservative ( [38,39,74]).
Because the currently available implementation of the above models allows the analysis of only two species at once, we prepared separate datasets, each including two of the analyzed species with three rows -two rows corresponded to sites with successfully resolved allelic states in both compared species and the third row contained SNPs of the outgroup species C. strumicae. The coalescent models assume free recombination among loci but no recombination within loci and require three sets of data for each pairwise species analysis. These consist of the number of nucleotide differences between pairs of sequences sampled a) both from species 1, b) both from species 2, and c) one from each species. To estimate all model parameters simultaneously, data from all three types of sequence pairs must be included but each pair of sequences must come from a different, independent locus [37]. Therefore, we randomly divided the analyzed loci into three nonoverlapping subsets. The intraspecific data were simply calculated as the number of heterozygous positions at each locus, which in fact represents the number of nucleotide differences between a pair of alleles brought together by segregation into a sampled individual. The preparation of the third dataset was more complex since it requires the comparison of two haploid sequences from different species, while our sequences originated from diploid individuals often possessing multiple heterozygotic positions with unknown phase. Therefore, we extracted the longest possible alignment of SNPs where both compared individuals (species 1 and species 2) had at most one heterozygous position by trimming the per-locus alignments of SNPs (similar to the trimming procedure done in [75]). Subsequently we randomly phased the remaining heterozygous SNP markers into either one or the other base and compared the number of differences between both species (this procedure in fact represents the comparison of two alleles drawn from species 1 and 2 respectively while keeping the allelic states resolved in the outgroup on a shortened alignment). In order to incorporate intrapopulation variability into our estimates, we combined the data from two individuals of each species. In doing so, we randomly sampled each locus from either one or the other individual (assuming free recombination among loci) but we always kept the sequence of SNPs within each locus from a single randomly selected individual to avoid the introduction of false recombination.
The relative mutation rates at all loci were estimated by comparison with an outgroup species [76]. C. strumicae was used as outgroup, whose divergence time from the ingroup was set to 17.7 Mya according to the penalized likelihood time-calibrated tree of [35]. Computation of 95% confidence intervals was based on the profile likelihood as described in [39]. Because of the computational time required, such confidence intervals were provided only for the best-fitting model for each dataset. In some cases, the likelihood profiling broke down prematurely, since forcing one parameter too far from its maximum-likelihood estimate rendered the whole model impossible to fit.
In such cases we indicate the nearest value that could be obtained, keeping in mind that the true 95% CI should be wider (Table 3). 5. The general applicability of patterns revealed in Cobitis was tested by a comparative study of published sexual and asexual fish hybrids. Russell [5] investigated the correlation between the genetic divergence between hybridizing fish species (distances in the cytochrome b gene corrected with the K2P model) and the level of postzygotic isolation, which was categorized by an index ranging from 0 (both hybrid sexes fertile) to 4 (both sexes inviable). The index value 2 assigned the stage when both hybrid sexes are viable but infertile, therefore preventing any effective interspecific gene flow. We amended Russel's study by introducing additional value (5) of postzygotic isolation index to those fish hybrids that have been documented to transmit their genomes clonally (gynogenesis, androgenesis) or hemiclonally (hybridogenesis). Altogether, performed literature search lead us to 13 cases of fish asexual hybrids, which were added to the database of Russell [5]. In single case we modified Russell's [5] data since he assigned R. rutilus x A. bramma hybrids with postzygotic isolation index 0.5 but Слынько [77] showed that such hybrids produce clonal gametes and can reproduce via androgenesis. Therefore, R. rutilus x A. bramma hybrids were assigned with index value 5. In accordance with Russell's data, the cytochrome b gene divergence was calculated from available sequences of parental taxa using the K2P correction using the Mega 5.0 software [78] (Appendix S5).
The genetic distances of the group 5 was compared with the other types of hybrids using the t-test after the normality of data was evaluated with the Shapiro-Wilk test.
To avoid phylogenetic dependence, Russell [5] considered each species in one cross only.
Published cases of asexual hybrids concern non-overlapping species pairs with three exceptions. In Cobitis, C. elongatoides has been involved in at least two crosses leading to naturally occurring asexuals (C. elongatoides-taenia and C. elongatoides-tanaitica) but since the original C. tanaitica-like mitochondrion has been lost (see above), we considered C. elongatoides -C. taenia cross only. On the other hand, Hexagrammos octogramus and Poeciliopsis monacha produce asexual hybrids by mating with two (H. otaki and H. agrammus) and three (P. lucida, P. occidentalis, P. latipina) congenerics, respectively. Therefore, we performed the t-tests several times with only one such cross per species.
Supporting Information Table S1. Allelic profiles of parents and their progeny within crossing experiments on the basis of seven microsatellite markers.  Appendix S4. SNP dataset. Every contig is represented by ten lines, starting with its name (ref_acc) on the first line and followed by nine lines corresponding to nine individuals. Individuals' genotypes are indicated by letters as follows: s, C. strumicae; p, C. pontica; e, C. elongatoides; t, C. taenia; n, C.
tanaitica. Variant positions are expressed using IUPAC notation and invariant positions are excluded.
Appendix S5. Summary of genetic divergences of fish species producing asexual hybrids.