Introduction

Speciation is a continuous process, often progressing with decreasing gene flow between diverging populations. According to the concept of the speciation continuum ^3,4, hybridization between relatively closely related lineages typically results in fertile and recombining progeny and may even lead to an evolutionary meltdown or homoploid hybrid speciation ⁵. Increasing genetic divergence between hybridizing taxa may further promote the emergence of transgressive phenotypes ⁶ or hybrid vigour but typically lowers the mean fitness of hybrids, which consequently impacts the establishment of postzygotic reproductive isolating barriers ⁵. In species with well-differentiated sex chromosomes, incompatibilities particularly affect the heterogametic sex (Haldane’s rule) ^7–9. Ultimately, speciation is completed when the accumulation of incompatible alleles in differentiating populations induces hybrid sterility or even inviability.

However, hybridization is also known to give rise to fertile, yet asexually reproducing hybrids in various animal and plant taxa. Asexual hybrids exhibit a broad spectrum of cytological mechanisms for the production of unreduced, often clonal, gametes, which range from entirely ameiotic processes (apomixis) to those involving more or less distorted meiosis (automixis) ^10,11. A relatively common gametogenic alteration leading to asexuality is premeiotic genome endoreplication (Fig. 1), found in a variety of hybrid fishes, amphibians, and reptiles ^12–16. By auto-duplication, from each chromosome usually an identical copy is produced to pair with during the meiotic prophase. Premeiotic genome endoreplication thus not only ensures clonal reproduction, but also allows hybrids to overcome problems in chromosome pairing that would otherwise lead to their sterility ^15,17. Hybridization and asexuality thus represent important evolutionary phenomena ¹⁰.

Fig. 1.

To explain the apparent link between asexuality and hybridization, the “balance hypothesis” has been proposed ¹⁸, which assumes that incompatibilities, accumulated among parental genomes, may disturb gametogenesis in hybrids leading to the formation of unreduced gametes. Hybrid asexuality may thus emerge in a “favourable evolutionary window” of genetic divergence between hybridizing species, which should be large enough to trigger particular gametogenic aberrations, but no too large to compromise the hybrid’s fertility or even viability. The recently proposed concept of the “extended speciation continuum” ⁹ explicitly links the emergence of asexual reproduction to the speciation continuum. While forming species ultimately reach complete intrinsic reproductive isolation under increased divergence, fertile hybrids may be produced at earlier stages that exhibit aberrant gametogenesis, leading towards clonal reproduction. This not only facilitates the formation of allodiploid and allopolyploid asexual hybrid lineages, but asexuality per se also effectively restricts introgression between parental species due to non-recombinant production of gametes. Hybrid asexuality may thus be considered a particular type of Bateson-Dobzhansky-Muller (BDM)-incompatibilities, representing a form of intrinsic reproductive isolation during the speciation process ¹⁷.

While the evolutionary link between hybridization and asexuality has been empirically documented, it remains unclear how widespread such a link is and how commonly might asexual gametogenesis occur in hybrids. On the one hand, naturally occurring asexual hybrids exist in many major animal and plant lineages, but appear rare, suggesting that the genomic pre-conditions required for asexual reproduction may be rarely met during hybridization ¹⁹. On the other hand, most of the knowledge about the evolution and biology of asexual hybrids comes from studies of naturally occurring lineages, while laboratory crossing experiments have been rarely performed to directly address the rate of clonal gametogenesis in hybrids. The questions how easily asexual reproduction is induced and why it is so often linked with a hybrid constitution thus remains unclear.

To estimate the incidence of clonal gametogenesis in a radiation of freshwater fish, we performed extensive crossing experiments of sexual species and investigated the gametogenic pathways in their hybrids, using European loaches of the family Cobitidae as a model. Cobitidae presents a speciose group of freshwater fish of South Asian origin that spread over most of the Palearctic region since late Eocene ²⁰. They colonized Europe in three distinct lineages ²¹, which diverged from each other around 17-20 Mya ^21,22. The so-called Adriatic and Bicanestrinia lineages colonized the Balkan and Anatolian regions, while Central and Eastern Europe have been colonized by the Cobitis sensu stricto lineage (Fig. 1a-b). This lineage is composed of two sub-groups, the so-called C. taenia clade, involving several closely related species diverged during the last 1 Mya, and C. elongatoides, which diverged ∼9 Mya from the C. taenia clade ¹⁷. Such independent colonization events and subsequent range shifts gave rise to natural hybrid zones across the Palearctic and Europe in particular ²³. Diploid and polyploid hybrids have been reported within the family, including gynogenetic ‘asexual’ lineages. These lineages reproduce via gynogenesis, i.e. clonal eggs require sperm from sexual males to activate further development into genetically identical progeny. After fertilization, male pronucleus is usually eliminated from the eggs, without genetic contribution to the progeny ^24,25 (Fig. 1a-c). Such forms usually discard the sperms ‘genome after fertilization and are currently found in hybrids of the Asian Misgurnus anguillicaudatus complex ²⁶, the Cobitis hankugensis-longicorpa hybrids in Korea ²⁷, and in hybrids of the European Cobitis sensu stricto lineage (Fig. 1b), ²⁸. In these natural systems, only hybrid females are known to produce unreduced gametes, employing premeiotic genome endoreplication, while hybrid males are usually sterile due to aberrant pairing of chromosomes in meiosis ^15,29.

In the present paper, we report meiotic and premeiotic gametogenic stages of experimental hybrids of eight loach species, representing the three main European lineages. These species have been crossed in various combinations from closely related taxa to phylogenetically distant ones, which allowed us to test, how often unreduced hybrid gametes arise, whether their emergence is linked to genetic and karyotypic divergence between parental species, and how asexuality is related to hybrid sterility.

Results

2.1 Parental species and their karyotypes

The selected parental species represent three distinct phylogenetic lineages from Europe (Fig. 1a, Supplementary Fig. S1) ²¹, possessing diverse karyotypes with diploid chromosome sets (2n) between 48 – 50 chromosomes and variable numbers of meta-/submeta- and subtelo-/acrocentrics (Supplementary Fig. S2, Supplementary Table S1). Karyotype dissimilarity among our focal species was measured as the autosomal karyotype index (AKD) which evaluates the differences in chromosome numbers and fundamental numbers of their arms ³⁰. We found that AKD significantly increases with the genetic divergence among analyzed parental species, evaluated as the p-distances in SNPs in coding sequences (Mantel test, r = 0.338, p = 0.0261). However, the correlation was not always linear (Supplementary Fig. S3). This is because some phylogenetically distant species belonging to the Bicanestrinia and Adriatic phylogroups possess morphologically almost identical karyotypes (Supplementary Fig. S2, Supplementary Table S1).

2.2 Artificial crosses and the presence of germ cells in all F1 hybrids

We crossed the eight focal species in 12 combinations, including some reciprocal crosses (Fig. 2). To investigate the effect of increasing phylogenetic distance, we predominantly used EE males and females for crosses with most species belonging to the three sampled phylogroups, but species combinations without C. elongatoides were performed too. We obtained viable progeny of both sexes from all attempted species combinations, altogether resulting in 19 F1 hybrid families (Fig. 2, Supplementary Table S2 Supplementary Fig. S1). Here we present the abbreviation for every hybrid combination: C. elongatoides × C. taenia (ET), C. elongatoides × C. tanaitica (EN), C. elongatoides × C. taurica (EC), C. elongatoides × C. pontica (EP), C. elongatoides × C. strumicae (ES), C. elongatoides × C. bilineata (EB), C. elongatoides × C. ohridana (EO), C. taenia × C. bilineata (TB), C. ohridana × C. bilineata (OB), C. tanaitica × C. pontica (NP), C. pontica x C. strumicae (PS). We observed no significant deviations from equal sex ratios among the F1 (Binomial test, all p values > 0.25 after corrections for multiple testing; Supplementary Table S2) and -about 6 months after hatching -randomly selected 3-5 juveniles per family to examine the gametogenesis in F1 hybrids.

Fig. 2.

Whole-mount immunofluorescence staining of gonads with antibodies against Vasa protein showed that the analysed F1 hybrids contained germ cells organized in clusters, which were similar to germ cell clusters of sexual fish, used as a control group (Supplementary Fig. S4a,c,e). The immunofluorescence staining of SYCP3 protein was performed either on whole-mount staining of gonadal tissue or on pachytenic spreads, and revealed the formation of lateral components of synaptonemal complexes (SCs) in most controls and F1 hybrids, suggesting the presence of meiocytes (Fig. 3,4 Supplementary Fig. S4b,d,f, S7, S9, S10). The only exceptions were PS and TB hybrid juvenile males, whose gonads contained clusters of germ cells without meiocytes. The absence of meiocytes was later verified also in both, adult PS and TB males (Supplementary Fig. S4a,b).

Fig. 3.

Fig. 4.

2.3. Analysis of meiocytes

We then investigated chromosome pairing during the pachytene stage of the meiotic prophase by immunostaining of the SYCP3 and SYCP1 components of the synaptonemal complex, and we analysed lampbrush chromosomes of the diplotene meiotic stage in females. We analysed 24 parental individuals, including males and females (EE, TT, NN, PP, OO) and 79 individuals among all F1 families (Supplementary Tables S2, S3).

2.3.1. Sexual species produce reduced gametes with properly formed bivalents

Depending on their chromosome numbers, parental species always contained 24-25 synaptonemal complexes (Fig. 3a; Fig. 4a) in 852 analysed cells (Supplementary Table S3). Flow cytometry of testicular cell suspensions from parental species revealed the presence of haploid and diploid cells with a major peak corresponding to haploid cells (Supplementary Fig. S5a). Histological analysis showed many sperm cells organized in cysts in all sexual males investigated, containing properly developed Sertoli and Leydig cells (Supplementary Fig. S6a), suggesting normal development of testes.

2.3.2. Sterility of hybrid males is caused by both, the aberrant pairing in the meiotic prophase and the inability to proceed to meiosis

Flow cytometry revealed that the analysed hybrid males lacked the haploid population of cells, characteristic for spermatids and spermatozoa of sexual males. Hybrid males had only diploid cell populations (Supplementary Fig. S5b-h). Interestingly, TB males had a small population of diploid cells corresponding to pachytene stage. These males had morphologically underdeveloped gonads, indicating that here the sterility is probably caused by the inability of gonocytes to proceed to meiosis (Supplementary Fig. S4a,b; Supplementary Fig. S5h). Histological examinations demonstrated that analysed hybrid males, had defective testes with an asynchronous development of germinal cells in cysts (Supplementary Fig. S6b-g). Similar results were reported in naturally occurring ET hybrid males ¹⁵.

Furthermore, hybrid males from most experimental crosses showed aberrant chromosomal pairing during the pachytene in 244 analysed cells (Supplementary Table S3). Abnormal cells had variable numbers of univalents, bivalents and multivalents (Fig. 3c; Supplementary Fig. S7a-g). The only exception to aforementioned observations were the NP hybrids involving closely related parental species (C. tanaitica and C. pontica) where we detected fully-paired chromosomes during the pachytene, in 36 cells (Fig. 3b). Chromosomal spreads of NP males also clearly showed the presence of formed spermatid nuclei, similar to sexual males (Fig. 3b). Additionally, we counted SCs, based on the presence of crossing-over foci (Supplementary Fig. S10a-f) ^15,31.

2.3.4. Production of unreduced gametes by premeiotic genome endoreplication is the prevalent gametogenic pathway among hybrid females

Analysis of diplotenic oocytes of most hybrid females revealed no univalents or mispaired chromosomes, but instead showed the formation of 49 or 50 properly formed bivalents (Supplementary Fig. S8a-e). Similar to previous data on natural clones ^13,29, we therefore conclude that diplotene oocytes in the majority of hybrid females contain only cells with duplicated genomes, suggesting their ability to complete clonal gametogenesis (Fig. 4; Supplementary Fig. S8a-e). The only exceptions were observed in OB and NP hybrids, where diplotenic oocytes contained only 25 fully formed bivalents, suggesting normal meiosis with the pairing likely between orthologous chromosomes (Supplementary Fig. S8f).

However, the analysis of 665 pachytenic oocytes revealed different patterns (Supplementary Table S3). OB and NP females possessed 25 bivalents with delayed pairing of two chromosome pairs (Supplementary Fig. S9d), generally matching their diplotene configurations. However, most F1-females contained pachytenic cells (n = 535) with abnormally-paired chromosomes and various numbers of univalents, bivalents, and multivalents (Fig. 4c; Supplementary Fig. S9a-c, e-h). This suggests that most pachytenic oocytes in hybrids did not undergo genome endoreplication. Pachytenic oocytes with duplicated genomes and thus normal chromosomal pairing were only observed in females from three F1-families (EC, EO and ET hybrids). These showed sets of 50 or 49 bivalents with properly loaded lateral and central synaptonemal components and MLH1 loci on each bivalent (Fig. 4d).

Overall, the incidence of genome endoreplication was very low among pachytenic cells of any F1 hybrid female, occurring on average in 0.7% of cells only (Supplementary Table S3) and there were no significant differences between any combinations of parental species (ANOVA Chi square test of generalized linear mode (GLM); DF = 9, dev. = 7.5, p-value = 0.59). We also did not observe any significant effects of maternal species in those crosses, where both directions were performed, i.e. there were no obvious differences between ET and TE as well as between OE and EO cross types.

We therefore treated all experimental progenies as a single group and tested whether they differed from previously-analysed natural diploid ET-hybrids, whose incidence of duplicated gonocytes was higher, on average ∼6%; ²⁹. Here, the difference was significant (ANOVA Chi square test of GLM; DF = 1; dev = 22.07; p-value < 10^-5), suggesting that natural clones produce a higher proportion of oocytes with endoreplicated genomes than experimental F1 hybrids.

The contrasting incidence of oocytes with endoreplicated genomes between pachytene (approx. 6 %) and diplotene stages (100%) is consistent with previous results ²⁹ suggesting that premeiotic genome endoreplication is a rare event in all types of crosses. This allows cells to form normal bivalents and to enter the diplotene, while the majority of cells remain non-duplicated and cannot proceed beyond the pachytene checkpoint due to mispairing among orthologous chromosomes.

Importantly, we noted that incidence of duplicated gonocytes was significantly higher among F1 hybrids (treated as a single group), than in females of any tested parental species, where in fact no duplications have been observed among 792 analyzed oocytes (ANOVA Chi square test of GLM; DF = 1; dev = 7.506; p-value = 0.00615).

2.3.5. Divergence of parental karyotypes affects the formation of bivalents between orthologous chromosomes in hybrids but has different effects in males and females

We than focused exclusively on those pachytenic cells with unduplicated genomes and found that F1-hybrids largely vary in numbers of paired orthologous chromosomes, ranging from 3 detectable synaptonemal complexes in EO males to 25 fully-formed bivalents in OB females. To explain such a variability, we fitted the distribution of numbers of synaptonemal complexes by variables assuming a) sex of the hybrid individual, b) chromosomal divergence between its parental species measured through the AKD index and c) their genetic divergence measured in exome-wide SNP p-distance. All three variables significantly contributed to the data distribution with numbers of bivalents in pachytenic cells being negatively correlated with morphological dissimilarity of parental karyotypes (GLMM z value = -5.86; p < 10^-8) and with their genetic distance (GLMM z value = -3.73; p < 10^-3). Males also tended to have fewer bivalents than females (z = -4.54; p < 10^-5). The model also suggested that both slopes are steeper in males than in females, however, the interaction of AKD or SNP divergence with hybrids’ sex was not significant (Fig. 5). This suggests that the formation of bivalents in hybrids is negatively affected by chromosomal dissimilarity as well as overall genetic divergence between parental species, and that male offspring generally form less bivalents than female progeny.

Fig. 5.

2.4. Comparison with other crossing experiments within Cobitidae

Finally, to put our results in context with previously published data, we plotted the reproductive outputs of all known loach hybrids in relation to genetic divergence and chromosomal differences between their parental species. To express the amount of chromosomal dissimilarity among parental species, we used the AKD index calculated from published and here newly prepared karyotypes of all spined loaches known to hybridize. To express the amount genetic divergence among these species, we used their sequence distances in the nuclear RAG1 gene, as a proxy, because exome-based sequence data are not available for some of the previously examined crosses. Data on RAG1- and karyotype divergences were then correlated with reproductive outputs in their hybrids and presented on Fig. 6a.

Fig. 6.

Inviable or fully sterile hybrids were generally found among hybrids between the genetically most distant species, while sexual and/or clonal reproduction were observed in hybrids between more closely related species. However, the window of genetic divergences allowing for sexual as well as asexual reproduction was relatively wide (Fig. 6a). Interestingly, we noticed that asexuality tends to be observed among female hybrids between species with substantially diverged karyotypes, while sexual reproduction generally required more similar karyotypes.

Discussion

Using freshwater fishes, loaches (Cobitidae), as a model, the present study reveals a wide range of reproductive outcomes in F1 hybrids in relation to genetic and karyotype divergences between the parental species (Fig. 6). On one extreme, hybrid inviability has been reported in crosses between distant species ³². On the other extreme, the production of haploid gametes was found in crosses of closely related species, such as between C. lutheri and C. tetralineata ³³, and between C. pontica and C. taenia ¹⁷ or C. tanaitica species, whose divergence has been dated to ∼1 Mya ¹⁷, but also in OB female hybrids in a relatively distant species pair (Fig. 2, Fig. 6, Supplementary Fig. S9d, S10d,). However, the most common types of reproductive outputs were hybrid sterility and production of unreduced gametes.

In case of hybrids between relatively far-distant species pairs, like the intergeneric Cobitis-Misgurnus cross, or between Adriatic and Cobitis sensu stricto lineages (i.e. EB and TB hybrids), sterility occurred in hybrids of both sexes (Fig. 6). However, the most common outcome was asymmetric, with sterility affecting hybrid males, while F1 females produced unreduced oocytes. This pattern occurred in the majority of crosses, including species pairs, like SE, SP or EO, with divergences reaching ∼17 or even >20 Mya back ²² (Fig. 2, Fig. 6) and it has also been reported from several natural hybrid complexes ^13,15,29. An intermediate situation was observed in TP-hybrids between a closely related species that diverged ∼1 Mya ¹⁷, where males produced reduced gametes and females produced both, reduced and unreduced ones.

It is generally assumed that the emergence of asexuality is a rare phenomenon that may require complex changes in gametogenic pathways and reproductive modes ¹⁹. This opinion has been corroborated by the “twiggy” distribution of ‘asexual’ organisms on the tree of life ³⁴ and by many failures to re-generate naturally occurring asexual hybrids by experimental crosses of their parental species (e.g. in apomictic Poecilia fish and parthenogenetic Darevskia lizards ^19,35.

However, we observed the production of unreduced gametes in 7 out of 11 crossed species pairs of loaches and in 10 out of 18 species pairs when including literature reports and natural hybrids (Fig. 2, Fig. 6). Such a widespread emergence of asexual gametogenesis therefore appears not to depend on the presence of a genome from any particular parental species, but to arise relatively commonly among independent species pairs across the entire subfamily Cobitinae. Asexual gametogenesis has also been observed in experimental F1 progeny of several other vertebrate species ³⁶. This suggests that the formation of unreduced gametes, i.e. the necessary prerequisite for asexual reproduction, may be a much more common outcome of interspecific hybridization than it would appear from the rarity of asexual hybrid taxa that are successfully established in nature.

Some studies also suggested that the emergence of asexuality requires the evolution of specific genes within hybrid lineages after their origin ³⁷. However, our finding that asexual gametogenesis originates relatively commonly already in the F1 generation somewhat refines this hypothesis and suggests that it is a direct result of genome merging by interspecific hybridization. Importantly, all F1-females capable of clonal gametogenesis employed an identical cytogenetic mechanism, i.e. premeiotic genome endoreplication, as also observed in many natural asexual vertebrates ^{13,15,16,29,36,38,39} suggesting that gametogenic perturbations induced by hybridization are often canalised into a similar developmental pathway. The reasons for such common patterns are not clear. The regulation of cell cycle is particularly sensitive to stress and perturbations ⁴⁰ and interspecific hybridization may act as such stimulus ¹⁸. While asexual reproduction may employ various pathways ¹⁰, it is possible that aberrations of cell cycle leading to endoreplication are more common since they require less accumulated incompatibilities than other routes to the production of unreduced gametes. Alternatively, the occurrence of endoreplication may reflect some predispositions, like the spontaneous production of unreduced gametes in some sexual parental lineages, which would then ensure clonal reproduction in their hybrids. However, our data argue against this alternative, since we analysed hundreds of gonocytes of sexual females without a single incidence of genome endoreplication ^15,29. Moreover, during years of crossing experiments and population sampling, we genotyped hundreds of C. elongatoides and C. taenia and their progenies and never found an autotriploid specimen ^23,41, which would be expected if fertilisation of unreduced gametes existed in these species ⁴².

Our in-depth analysis of hybrid gametogenesis in Cobitis loaches further provided compelling evidence for a tight link between hybrid incompatibility, asexuality and sterility. Namely, the majority of meiocytes produced by both the male and female hybrids failed to form a full set of bivalents between orthologous chromosomes and were arrested at meiotic checkpoints (Fig. 3). However, unlike their sterile brothers, hybrid F1 females also usually contained a minor proportion of oogonia that underwent premeiotic genome endoreplication. In most our F1 hybrid females, only these rare oogonia were able to complete meiosis and form diploid gametes. Interestingly, similar observations were reported from several natural and experimental hybrid vertebrates, like ET and EN hybrid loaches, geckos and whiptail lizards (Aspidoscelis) Medaka ricefish (Oryzias) ^16,29,36,43. This suggests that even successfully established hybrid asexuals usually contain only a minor proportion of oogonia, which actually undergo clonal gametogenetic pathway.

Such congruent patterns across phylogenetically distant taxa thus conform to the concept of an “extended speciation continuum”. Similar to the classical concept of the “speciation continuum”, it assumes that speciation generally proceeds from early stages characterized by the production of sexual hybrid progeny to irreversibly diverged species, whose hybrids are either sterile or even inviable. At intermediate stages, however, hybrids between diverging evolutionary lineages may not only have a certain incidence of sterile gonocytes, but also tend to produce unreduced gametes (Fig. 6). Genome endoreplication thus appears as an effective way to overcome sterility in hybrid females, even at high levels of divergence. Asexuality may thus often represent the only reproductive pathway of hybrids between a given pair of species and, together with hybrid sterility, it may simultaneously contribute to reproductive isolation of parental species due to clonal genome propagation.

We observed that the fertility of female hybrids remains conserved due to sex-specific emergence of clonality, while males of the same cross are usually sterile. This may potentially indicate some analogy to other classical sex-related asymmetries in hybrids, like incompatibilities between sex chromosomes (Haldane’s rule) or other uniparentally inherited factors (Darwin’s corollary ⁴⁵). However, the sex determination systems are not yet known in most of the investigated species ⁹ and in any case, the reasons for such an asymmetry are likely more complex. For instance, meiocytes in hybrid males contained significantly lower numbers of properly-formed synaptonemal complexes than non-duplicated meiocytes in females (Fig. 5) and were more severely affected by decreasing similarity of parental karyotypes (see Results, 2.3.5.). This indicates some more fundamental differences in pairing mechanisms and affinity among orthologous chromosomes fundamentally between spermatogenesis and oogenesis, e.g. ⁴⁶. In fact, we recently found that spermatogonial cells taken from a sterile ET hybrid male may resume the ability of endoreplication when transplanted into female gonads ⁴¹. This suggests that initiation of premeiotic endoreplication may rely on female-specific gonadal environment rather than on genetic sex determination.

The genetic divergence among parents may affect early gametogenic stages in hybrids through various mechanisms causing e.g. depletion of primordial germ cells ⁴⁷, improper pairing of diverged chromosomes ⁴⁸, asynapsis between sex chromosomes ⁴⁹, or deleterious epistatic interactions between orthologous genes and their binding motifs on chromosomes ⁵⁰. In loaches, we documented at least two such mechanisms underlying hybrid sterility. TB and PS hybrid males had almost no spermatocytes, suggesting their gonocytes probably failed the transition to meiosis. The other type of sterility conforms to a chromosomal speciation model ^51,52, prevailing in most F1 hybrid males, in which spermatocytes developed but possessed aberrantly-paired chromosomes, preventing further progression to spermatids. Analogous patterns were observed in non-duplicated oocytes of most hybrid females, suggesting that problems in bivalent-formation have a detrimental impact on the meiosis in both male and female loach hybrids.

Karyotype evolution among parental species indeed plays a crucial role in hybrids’ reproductive capabilities ⁵². Our data suggest that it has a crucial, yet indirect, impact on asexual reproduction. Specifically, hybrids between species with diverged karyotypes had less bivalents in non-duplicated meiocytes than those combining similar karyotypes. Sometimes, even hybrids between phylogenetically distant, yet karyotypically similar species could produce full sets of bivalents (e.g. OB hybrids) (Fig. 5). On the other hand, the numbers of endoreplicated oocytes were low in any hybrid. Their proportion was relatively higher in natural clones than in experimental F1 Cobitis, suggesting that such a trait may recover during the evolutionary establishment of a natural asexual lineage ²⁹ and/or that selection among “freshly” emerged clones would favour those with the highest proportion of endoreplicated oocytes. Nevertheless, even the most widespread clones possessed only ∼6% of such oogonia, suggesting that karyotype similarity allows the completion of meiosis and hybrids would mostly reproduce sexually by reduced gametes, even if they had the capacity of producing unreduced gametes, too.

While hybridization is assumed to present an important trigger of asexuality in vertebrates and other organisms ¹⁰, the origins of clonal lineages are thought to be rare, and often can be traced back to a single or few hybridization events in a given species pair ^19,28. Multiple independent origins of clonal gametogenesis in spined loaches ^23,28, this study) in part contradicts this view since interspecific hybridization has been found to result in clonal gametogenesis relatively commonly up to a substantial divergence, where fertile or viable hybrids are no longer appearing (Fig. 6a). However, we emphasize that successful establishment of hybrid clones under natural conditions is likely restricted by several additional levels of selective pressure and thus filtering (Fig. 6b). On the one hand, the “balance hypothesis” ¹⁸ and its expansions and fine-tuned versions ^17,19, suggest that the production of unreduced gametes probably requires some level of genetic divergence among the parental species to alter the hybrid cell cycle towards genome endoreplication. On the other hand, it appears that most meiocytes of asexual hybrids anyway do not pass through endoreplication ^16,29,36,43. Consequently, if the parental karyotypes remain largely similar, then orthologous chromosomes may successfully pair even in such nonduplicated meiocytes, potentially, allowing them to accomplish normal meiosis and form reduced gametes. This would decrease the proportion of unreduced gametes produced by such hybrids and negatively impact on their likelihood to establish a stable asexual lineage.

We therefore propose, that for a cross to be successful in generating an asexual hybrid, it not only has to produce a sufficient number of unreduced gametes, but it also has to prevent the formation of reduced gametes from non-duplicated meiocytes. One way to achieve this may be if hybridization involves species with sufficiently diverged karyotypes, ensuring improper pairing of orthologous chromosomes in non-duplicated gametes. Increased divergence of parental karyotypes may therefore have a positive impact on the production of stable clonal hybrids and, for instance, the fixation of certain chromosomal rearrangements among diverging populations may be accelerated by certain geographical aspects of speciation, such as parapatry ^51,52. It would thus be attractive to test if the incidence of asexual hybrids correlates with increased rates of karyotype evolution among parental species.

The rarity of naturally occurring clonal organisms probably results from the simultaneous impact of the abovementioned selective processes. However, our data indicate that the production of unreduced gametes, which is the very basic conditio sine qua non for the evolution of asexuality, may be surprisingly “easy” to trigger, especially in hybrids of distantly related lineages. Of course, we may not rule out that the frequent emergence of premeiotic endoreplication in Cobitidae reflects some specific predispositions in their gametogenic machineries in comparison to other taxa, where asexuality appears rarer. Thus, our study demonstrates that conclusions about the evolution of asexuality should not only be drawn from natural asexuals, but should be based also on experimental crosses.

Acknowledgements

The authors would like to thank Sharka Pelikanova (Institute of Animal Physiology and Genetics), Daniel Kulik and Grzegorz Skórzewski (Museum of Natural History, University of Wroclaw) for helping with taking care of the fish; Ruxanda Cucerenco (Nicolae Testemițanu State University of Medicine and Pharmacy, Chisinau) and Anton Svinin (University of Tyumen) for helping with preparation of samples with immunofluorescent staining, Sean Burgess (College of Biological Sciences, University of California) for providing primary antibodies against SYCP1; Vladislav Vasiullin for the help with the preparation of illustrations; Stephen A. Schlebusch (Charles University, Prague) and Roger Butlin (University of Sheffield, Sheffield) for reading the first draft of the manuscript.

Funding

This study was supported by following grants:

Czech Science Foundation grant Nos. 19-21552S, 21-25185S, RVO 67985904, and Ministry of Education, Youth and Sports of the Czech Republic (grant No. 539 EXCELLENCE CZ.02.1.01/0.0/0.0/15_003/0000460 OP RDE. (AM, DD, KJ)

Ministry of Education and Research in Moldova grant DIVZOO 20.80009.7007.12 (AM)

Author contributions

Conceptualization: AM, DD, JB, KJ

Methodology: AM, DD, VSB, KH, LC, TT, AS, KJ

Investigation: AM, DD

Visualization: AM, DD

Funding acquisition: KJ, AM

Project administration: KJ

Supervision: KJ, DD

Writing – original draft: AM, DD

Writing – review & editing: KJ, MS

Competing interests

Authors declare that they have no competing interests.

Data and materials availability

All data are available in the main text and the supplementary materials and raw reads from sequence capture will be made available on the SRA (Sequence Read Archive) NCBI upon publication.

Materials and Methods

Selection of parental specimens and genotyping

All Cobitis specimens were obtained from the rearing facility stocks of the Laboratory of Fish Genetics, IAPG CAS CZ that had been collected during recent projects across Central Europe in accordance with environmental protection legislation. The Valid Animal Use Protocol was in force during the study at the IAPG (No. CZ 02386). All institutional and national guidelines were covered by the “Valid Animal Use Protocol” No. CZ 02386 of the Laboratory of Fish genetics, IAPG CAS. For the breeding experiments, we selected 38 individuals of 8 species. Taxonomic identification and genotyping of examined individuals were based on previously determined and routinely applied molecular markers involving two Sanger-sequenced nuclear markers (the intron in S7 gene and RAG1 gene), and one mitochondrial gene, cytochrome b (CytB), which were compared to previously published data to confirm their taxonomical identification using routine protocols described in ²³. Experimental design of crosses depended on available fish in our aquariums. So far, the most abundant, in males and females, were EE, which was bred with all other species. In the case of breeding of TT with PP, BB, crosses were unique or had only two replicates because of the limited number of available fish.

To obtain phylogenomic data and evaluate SNP, we applied the exome-capture approach to all experimentally crossed lineages ¹⁷. Briefly, isolated genomic DNA was sheared with Bioruptor™, tagged by indices, hybridized to custom-designed exome-capture probes ¹⁷ sequenced with 2*75 mode on Illumina NextSeq device. Raw data were processed (adapter and quality trimming) using bbduk.sh script from BBmap package ⁵³; ktrim=r k=23 mink=11 hdist=1 tpe tbo qtrim=rl trimq=10 maq=10 minlen=36). Processed reads were aligned to C. taenia reference transcriptome that was published and cleaned from potentially paralogous contigs by (Janko et al. 2018). Mapping was performed with bwa mem algorithm ⁵⁴ in default settings and SNP calling with GATK4 pipeline ⁵⁵ (indels were excluded).

Phylogenetic inference and estimation of pairwise genomic distances

Mitochondrial loci were eliminated from all downstream analyses of SNP, which focussed only on nuclear exon sequences. The VCF dataset was used to calculate pairwise p-distances in SNPs between all individuals using VCF2Dis v1.47 software ⁵⁶.

In order to reconstruct phylogenetic relationships among crossed species, we further used the GATK4 FastaAlternateReferenceMaker option to create locus-specific consensuses from each sample. Regions with low read depth in each new consensus were identified using the samtools software ⁵⁷ with depth option and the output file was rewritten into bed format where sites with depth < 10 were masked by ‘N’ using bedtools ⁵⁸ and maskfasta option. To mitigate the locus dropout in distant species, we than used custom R scripts (R Core Team 2020) with functionalities of the seqinr package ⁵⁹ to select the alignments where reads from all investigated species are present. Final phylogenetic analysis was thus based on 1960 loci with length > 750 bp & > 30 parsimony informative sites, where all species had correctly read sequence variants on more than 70% of sites. Individual Maximum Likelihood gene trees were reconstructed by IQ-TREE v. 2.0.3 ⁶⁰ using the extended model selection with free rate of heterogeneity in combination with 1000 ultrafast bootstrap replicates ^61,62. Consensus species tree was estimated by ASTRAL v. 5.5.6 ⁶³. Gene Concordance Factors (gCF) and site Concordance Factors (sCF) were estimated ⁶⁴ and resulting trees were processed using the ape package ⁶⁵ and plotted by DensiTree v2.0 ⁶⁶.

Crossing experiments

Artificial spawning was performed by hormonal stimulation of fish. Mature fish (males and females) were injected twice in the peritoneal cavity (24 and 12 hours before fertilization) with hormone Ovopel (Interfish Kf) ⁴¹. First injection was applied with a solution of one Ovopel pill per 20 mL of 0.9% NaCl. Second injection was done with a solution of 1 pill per 5 mL of 0.9% NaCl. The ratio of injected solution in both cases was 0.05 mL per 10g of fish weight. Fish eggs were gently squeezed 24 hours after first injection and transferred to a Petri dish (ø = 8.6 cm, 1.5 cm high) with a density not higher than 100 eggs per dish. Male sperm was also obtained by squeezing or by homogenization of male testes from euthanized fish and transferred in Hanks’ balanced salt solution (Sigma-Aldrich) or directly applied on loach eggs together with fresh water to activate spermatozoa. Incubation of eggs was performed at room temperature (23 – 25C). After hatching, free larvae were transferred into plastic pots with the following size: 25×25×15cm. Larvae were fed twice daily with nauplii of the brine shrimp Artemia. After 45 days from hatching, juveniles were fed with Tubifex worms once per day. Adult fish were fed with tablets of dry food when live food was not available. Starting from two weeks we randomly selected fish larvae for the whole mount investigation of fish gonads. For immunofluorescent staining we randomly selected juveniles with an age of two months after hatching. Analysis of diplotene chromosomes was performed on adult and subadult females older than 6 months. Adult and subadult males older than half a year were used for flow cytometry measurements of the cell suspension from their gonads, histology and pachytene chromosome analysis.

DNA flow cytometry

DNA content of the cell suspension from testes was measured by BD FACSAria™ flow cytometer. Dissected testes were collected into Versene solution (Thermo Fisher Scientific). In total, we investigated 22 hybrid males from 19 families and two sexual species (EE). To release and stain the cell nuclei, testes were homogenized and incubated with 0.1% Triton X100, 10 µ/ml DAPI and 15 mM MgCl2 at +4°C overnight. At least 10,000 events were analysed. The samples from testes of sexual diploid species EE, were used as an internal control for measurements. Data was further analysed by BD FACSDiva software (version 6.1.3). We assessed haploid cells to 1n1C peak, and diploid cells to 2n2C and 2n4C peak (2n2C cell population represent cells before DNA synthesis, while the 2n4C cell population represents cells after the S phase, possibly indicating meiocytes arrested during pachytene due to chromosomal mispairing in meiosis), where n represents a set of chromosomes, and C represents the number of chromatids.

Preparation of mitotic chromosomes

Cell suspensions with mitotic chromosomes from all specimens were obtained either from kidneys and/or regenerated caudal fins according to ²². Metaphase chromosomes were stained with Giemsa to check the morphology of chromosomes. Chromosome formula = 2n (48-50) metacentric, submetacentric/acro- and subtelocentric.

Pachytene chromosomes preparation and immunofluorescent staining

To obtain pachytene chromosomes we used protocols described in ^15,29 for males and females, respectively. Gonads were manually homogenized in 1× PBS. In the case of males, we incubated 1 μl of suspension in 30 μl of hypotonic solution (1/3 of 1× PBS), preliminarily dropped on SuperFrost® slides (Menzel Gläser) for 20 minutes. Subsequently slides were fixed in 2% paraformaldehyde (PFA) for 4 minutes. In the case of females, 20 μl of cell suspension was put on SuperFrost® slides (Menzel Gläser) followed by addition of 40 μl of 0.2 M Sucrose and 40 μl of 0.2% Triron X100 for 7 min and subsequently fixed in 2% PFA for 16 minutes. After fixation, slides with male and female pachytene spreads were air dried, washed in 1× PBS and used for immunofluorescent staining.

During immunofluorescent staining we visualized synaptonemal complexes (SC) of pachytene chromosomes using rabbit polyclonal antibodies (ab14206, Abcam) against SYCP3 protein (the lateral component of SC) and chicken polyclonal antibodies (a gift from Sean M. Burgess) against SYCP1 (the central component of SC). Crossing-over foci were identified using mouse monoclonal antibodies (ab15093, Abcam) against MLH1 protein, a mismatch repair protein. Slides were incubating with 1% blocking reagent (Roche) in 1× PBS and 0.01% Tween-20 for 20 min. Primary antibodies (dilutions as recommended by manufacturers) were added for 1h at RT followed by 3 times washing in 1× PBS. Afterwards, corresponding secondary antibodies were added for 1h at RT: (Alexa-594-conjugated goat anti-rabbit IgG (H+L) (Invitrogen), Alexa-594-conjugated goat anti-chicken IgG (H+L) (Invitrogen) and Alexa-488-conjugated goat anti mouse IgG (H+L) (Invitrogen). Slides were washed in 1× PBS with 0.05% Tween-20 and mounted in Vectashield/DAPI (1.5 mg/ml) (Vector, Burlingame, Calif., USA).

Diplotene chromosome preparation

Diplotene chromosomes (so called lampbrush chromosomes) were prepared from parental and hybrid females accordingly ^15,29. Ovaries with vitellogenic oocytes of size 0.5–1.5 mm in diameter were taken from females and placed in OR2 saline (82.5 mM NaCl, 2.5 mM KCl, 1 mM MgCl₂, 1 mM CaCl₂,1mM Na₂HPO₄, 5 mM HEPES (4-(2-hydroxyethyl)-1-piperazineethanesulfonic acid); pH 7.4). Individual oocytes were transferred to the isolation medium “5:1” (83 mM KCl, 17 mM NaCl, 6.5 mM Na₂HPO₄, 3.5 mM KH₂PO₄, 1mM MgCl₂, 1 mM DTT (dithiothreitol); pH 7.0–7.2) for microsurgical isolation of nuclei. Nuclear envelopes were manually removed in one-fourth strength “5:1” medium with the addition of 0.1% PFA and 0.01% 1M MgCl₂ in glass chambers attached to a slide. This procedure ensures the presence of a chromosomal set from an individual oocyte in each separate chamber. Slides with oocyte nuclei contents were centrifuged for 20 min at +4°C, 4000 rpm, fixed for 30 min in 2% PFA in 1× PBS, and post-fixed in 50% and 70% ethanol. After drying chromosomal spreads were mounted in Vectashield/DAPI (1.5 mg/ml) (Vector, Burlingame, Calif., USA).

Whole-mount immunofluorescence staining

Gonads of fish larvae were isolated and fixed in 2% PFA in 1× PBS followed by washing in 1× PBS. Tissues were stored until usage in 1× PBS with the addition of 0.02% sodium azide. Prior to immunofluorescent staining, gonadal fragments were permeabilized in a 0.5% solution of Triton X100 in 1× PBS for 4-5 hours at RT and washed in 1× PBS at RT. The following primary antibodies were used: rabbit polyclonal antibodies DDX4 antibody (C1C3, GeneTex) against vasa protein; rabbit polyclonal (ab14206, Abcam) against SYCP3 protein and chicken polyclonal against SYCP1 protein (a gift from Sean M. Burgess). After incubation for 1-2 hours in a 1% blocking solution (Roche) dissolved in 1× PBS, primary antibodies were added for 12 hours at RT. After three times washing in n 1× PBS with 0.01% Tween (ICN Biomedical Inc), secondary antibodies Alexa-488-onjugated goat anti-rabbit IgG (H+L) (Invitrogen) and Alexa-594-conjugated goat anti-chicken IgG (H+L) (Invitrogen) were added for 12 hours at RT. After three times washing in 1× PBS with 0.01% Tween (ICN Biomedical Inc) tissues were stained with DAPI (1 mg/ml) (Sigma) in 1× PBS at RT overnight.

Confocal laser scanning microscopy

Tissues fragments were placed in a drop of Vectashield antifade solution containing 1 mg/ml DAPI. Confocal laser scanning microscopy was carried out using a Leica TCS SP5 microscope based on the inverted microscope Leica DMI 6000 CS (Leica Microsystems, Germany). Specimens were analyzed using HC PL APO 40х objective. Diode and helium-neon lasers were used to excite the fluorescent dyes DAPI and Cy3, respectively. The images were captured and processed using LAS AF software (Leica Microsystems, Germany).

Wide-field and fluorescence microscopy and image processing

Mitotic chromosomes were examined by Olympus BX53 epifluorescence microscope and Axio Imager Z2 microscope equipped with CCD camera (DP30W Olympus) and CoolCube 1 b/w digital camera (MetaSystems, Altlussheim, Germany), respectively. Axio Imager Z2 epifluorescence microscope is equipped by Metasystems platform for automatic search, capture and image processing. Meiotic chromosomes were analysed using Olympus BX63 microscopes equipped with standard fluorescence filter sets, and sCMOS camera (Prime95B Teledyne Photometrics) using CellSense Dimension software (Olympus). The IKAROS imaging program (Metasystems, Altlussheim, Germany) were used to analyse grey-scale images. The captured digital images from FISH experiments were pseudocoloured (blue for DAPI, red for Alexa-594, green for Alexa-488) and superimposed using Microimage and Adobe Photoshop software, version CS5, respectively. Microphotographs were finally adjusted and arranged in Adobe Photoshop CS6 software.

Karyotype differences between species

Chromosomal data for each species was obtained from the literature and present work (Supplementary Fig. S2 and Supplementary Table S1). The inspection of karyotypes was done visually and we noticed the numbers of meta-/submeta- and subtelo-/acrocentrics for each species. The magnitude of the autosomal karyotypic differences between focal species was calculated using the autosomal karyotype index (AKD) ³⁰, which is calculated as the sum of absolute differences in diploid numbers of chromosomes (2n) divided by two and the absolute differences in the autosomal fundamental numbers of arms (NF) also divided by two.The data are presented in Supplementary Table S2.

Statistical analyses

To investigate whether karyotype divergence among species correlate with their genetic divergence, we used a Mantel test with 10,000 replicates to compare a matrix of pairwise chromosomal divergences expressed in pairwise AKD index and a matrix of genetic divergences expressed in pairwise exome-wide SNP distances.

We than investigated whether proportions of endoreplicated oocytes at the pachytenic stage differ among various types of F1 crosses. To do so, we compared the counts of duplicated and non-duplicated oocytes in every female using the generalized linear mixed effect model (GLMM) with individuals taken as a random factor and binomial error distribution. Model significance was evaluated with the Wald Z-test implemented in the R library CAR ⁶⁷.

Finally, we tested whether numbers of bivalents formed by F1 hybrids are explicable by the amount of karyotype differences measured by the AKD index of morphological differentiation among their parental species. To do so, we fitted the numbers of bivalents to hybrid’s sex and AKD as well as SNP distances between parental species by a GLMM model with individuals taken as a random factor and Poisson distribution of error in R library LME4 ⁶⁸. We performed the forward selection with each pair of nested model tested by analysis of deviance in GLMs, including their interaction with hybrids’ sex. Z values and Pr(>|z|) are reported as z and p, respectively.