Capture and return of sexual genomes by hybridogenetic frogs provides clonal genome enrichment in a sexual species

Hybridogenesis is a reproductive tool for sexual parasitism. Hybridogenetic hybrids use gametes from their sexual host for their own reproduction, but sexual species gain no benefit from such matings as their genome is later eliminated. Here, we examine the presence of sexual parasitism in water frogs through crossing experiments and genome-wide data. We specifically focus on the famous Central-European populations where Pelophylax esculentus males (hybrids of P. ridibundus and P. lessonae) live with P. ridibundus. We identified a system where the hybrids commonly produce two types of clonal gametes (hybrid amphispermy). The haploid lessonae genome is clonally inherited from generation to generation and assures the maintenance of hybrids through a process, in which lessonae sperm fertilize P. ridibundus eggs. The haploid ridibundus genome in hybrids received from P. ridibundus a generation ago, is perpetuated as clonal ridibundus sperm and used to fertilize P. ridibundus eggs, yielding female P. ridibundus progeny. These results imply animal reproduction in which hybridogenetic taxa are not only sexual parasites, but also participate in the formation of a sexual taxon in a remarkable way. This occurs through a process by which sexual gametes are being captured, converted to clones, and returned to sexual populations in one generation.


Results
Microsatellite taxon assignments. Experimental crosses between P. esculentus males and P. ridibundus females (Fig. 1A) resulted in 34 families and 523 juveniles that reached metamorphosis (Supplementary  Table S1). Microsatellite analyses performed on 18 parental individuals and 220 progeny across 17 families with metamorphosed juveniles detected a total of 71 alleles, from which eight were known as lessonae-specific and 63 alleles as ridibundus-specific (Table 1). Microsatellites supported the taxon assignments of adults based on morphology and further assigned RR/RL genotypes to 60/160 juveniles, respectively (Table 1; Supplementary  Tables S2 and S3). The presence of one or two alleles per locus together with allele peak sizes (dosage effect) suggested that all individuals were diploid. Allele-dosage effects as indications of polyploidy were not observed. Nine out of 12 hybrid males had only sons (always of RL genotypes). Two males (M4, M11) had mixed progeny of daughters (always of RR genotypes) and sons (always of RL genotypes) with sex ratios of approximately 1:1. Hybrid male M12 had only three RR daughters (Table 2). Four juveniles had no fully developed gonads making it impossible to determine their sex.
Inheritance and relatedness of parental genomes. Here, we examined the heredity modes of RL hybrid males concerning the fate of their L and R genomes. We assumed that RL males are fixed F1´s and, from a definition of hybridogenesis, eliminate one parental genome from their germline and transfer only the second parental genome to gametes. To do so, we designed methodological approaches allowing us to trace DNA variants of progeny backward to their parentals and separated a genome of the father from that of the mother. Phased haploid genomes were then analyzed for their relationships within families, taxon assignments, and levels of recombination using Principal component analysis (PCA), Principal coordinate analysis (PCoA), fineRADstucture heatmap, and STRU CTU RE ancestry bar plots. We denoted sexual genomes of P. ridibundus as R and of P. lessonae as L, the hybridogenetic hybrid P. esculentus as RL‚ and the clonally inherited genomes as [R] and [L].
GenAlEx estimated a total of 204 multilocus genotypes (MLGs) based on microsatellite profiles of 220 phased offspring. A total of 11 MLGs were found within lessonae genomes of 160 RL sons. A total of 193 MLGs were found within ridibundus genomes of 160 RL sons and 60 RR daughters. Inspecting the allelic variability by filtering out variation caused by missing data showed RL sons shared the same MLG. This indicates the clonal inheritance of the L genome (further [L]) in all but one RL son (60-2013 JUV19), which inherited the R-specific allele at locus Re1Caga10 from the father (Supplementary Figures S1 and S2). In contrast to hybrid sons, the paternal ridibundus genomes in RR daughters exhibited three MLGs, corresponding to the allelic profiles of three P. esculentus fathers. Daughters and fathers shared the same MLG within each family, indicating that R genomes (further [R]) were inherited clonally as well (Supplementary Figure S1). PCoA run in GenAlEx and PCA performed with R package radiator supported two groups of microsatellite MLGs; the two principle coordinates accounted for 65.16%/14.60% (GenAlEx/radiator) for axis 1 (Fig. 1D). Of the 15,676 SNPs called from the RADseq data, 2219 could be phased and were present in at least both parents and 75% of offspring in one or more families. The dendrogram of phased haplotypes showed more variability in maternal recombined genomes than in [R] and [L] genomes inherited from P. esculentus fathers (Supplementary Figure S4). PCA was used again to summarise this phased dataset and discriminated between the L and R genomes. PC 1 (78,60% in radiator), PC 2 (3,10%) and   Figure S2), even though algorithms of both Lep-Map3 and MSTmap are not designed for data with non-Mendelian inheritance patterns. In general, we revealed a complete linkage of loci within the germline of RL males. The linkage maps of RR females consisted of several linkage groups of different lengths, indicating functional crossing-over and recombination (Supplementary Figure S2).

Discussion
Our study shows that asexual organisms possess much higher reproductive variation than previously thought. The results of genetic analyses carried out on the basis of SNP markers and microsatellites indicate the absence of recombination between the L and R genomes in P. esculentus males. We could not exclude, however, an unrecognized leakage of subgenomic DNA due to general marker limitations. Among other asexual taxa, hybridogenetic P. esculentus males are exceptional in their ability to generate individuals of their parental species (Fig. 1A)   www.nature.com/scientificreports/ in relatively substantial numbers. More importantly, this reproductive strategy indicates mechanisms of non-Mendelian inheritance that may change our view of the role of asexuality in the dynamics of mixed population systems in which most asexuals exist.
Clonal inheritance of both paternal genomes. Both the microsatellite and fineRADstructure analyses, and PCA of the phased genomes, revealed the variation of inheritance patterns in P. esculentus‚ allowing fathers to produce not only clonal [L] but also [R] gametes, or both ( Fig. 1; Supplementary Figure S3). Allelic profiles compared between mother and offspring verified the Mendelian inheritance (i.e., with recombination) of R genomes coming from P. ridibundus females (Supplementary Figure S1). The comparison of male and offspring haplotypes, however, showed the complete physical linkage between markers in haploid [L] and [R] genomes from P. esculentus males, indicating their non-Mendelian inheritance from the father. The level of relatedness between the clonal [L] genomes was similar among individuals and might represent a single hemiclonal lineage. Indeed, allelic variation compared with previously analyzed hybrid males 35 supported the sharing of lessonae-specific alleles within the [L] genomes, containing a Y chromosome or at least a male sex-determining region. Clonal [R] gametes from fathers and sexual R gametes from the mothers were characterized by low divergence. Clustering of the phased genomes showed a pattern in which a terminal branch length varied between the sexual R genomes while clonal ones were characterized by a much shorter branch length (Supplementary Figure S4).
Our data supported the absence of recombination between R and L genomes except one P. esculentus offspring, which inherited a single R microsatellite allele instead of an expected [L] allele from its father. Traces of occasional introgression, called genome leakage, have already been identified in water frogs 22,[35][36][37] . The recent example of introgression has been linked to unexpected paternal genomic elements in the celibate genome 38 . Whether the recombination of P. esculentus offspring was caused by disturbed segregation during meiosis, intrachromosomal recombination, or through paternal genomic elements remains unclear.
Hybrid amphispermy as a reproductive strategy. The phenomenon of hybrid amphispermy, when even a single P. esculentus male can form sperm with clonal genomes of two different species, seems to be widespread 22,32,[39][40][41][42][43] . A similar strategy is assumed in the Australian carp gudgeons, which also exhibit hybridogenetic heredity of both parental genomes, although the pattern has not been observed at the individual scale; a male-biased group clonally inherits the paternal genome, while a female-biased group, clonally inherits the maternal genome 11 .
In our study, at least one-quarter of water frog hybrids formed both types of progeny, P. esculentus sons, and P. ridibundus daughters. Besides a pattern confirming an XY sex determination system, the total number of P. ridibundus daughters produced by hybrid fathers was quite impressive: 67 out of 274 offspring. Previous experimental data support a broader scale of the phenomenon in which even one-third of progeny were P. ridibundus females 22 . In the case of vertebrate systems that mirror hybrid female asexuality, one would assume a selective evolutionary preference for males to produce 100% hybrid progeny, providing them a two-fold reproductive advantage compared with sexual propagation 44 .
Rather than a failure of a proper genome elimination of one parental genome, we propose that P. esculentus´ male amphispermy is a result of an evolutionary strategy providing benefits to asexual males in coping with mating challenges. Sperm-dependent asexuals may face some degree of sperm-limitation 2 . In spite of this, asexual females sharing sexual males do not need to produce males on their own as sexual males are usually willing to spawn repeatedly during the whole reproductive season. The opposite situation can be expected in the R-E system in which P. esculentus males compete with P. ridibundus males for P. ridibundus females. The ovarian cycle of water frog females in the temperate zone of Europe provides a finite number of oocytes that will be ovulated in the next spawning season 45 . Sexual females are thus a much more limited source of mates than sexual males. Under such circumstances, the production of P. ridibundus daughters by P. esculentus males may increase the number of females in mixed populations and, in turn, the chance of hybrid males to successfully find a mate.
Genomes jumping from sex to asex and back to sex. Our crossing experiments allowed us to investigate a remarkable mechanism of inheritance operating through the interaction between both non-sexual hybrids and sexual species. Pelophylax esculentus males themselves were originally formed through the fusion of ancestral R gametes provided by P. ridibundus and ancestral L gametes obtained from P. lessonae, both representing species with Mendelian inheritance. However, once trapped within a hybrid, the [R] and [L] genomes usually do not recombine during spermatogenesis. The recombinant R gamete from a sexual species is assumed to be removed during the DNA elimination process in the hybrid germline, and in this way to be permanently lost from the sexual population. According to current data, the recombinant R genome is only temporarily misplaced from the sexual population rather than eliminated. In fact, the R genome is essentially ´frozen´ in some P. esculentus males for a single generation before it is returned, through non-Mendelian inheritance, back to the sexual population in the form of P. ridibundus females. Presumably, as the non-Mendelian inheritance of these R genomes has only persisted for one generation, deleterious effects are negligible. Therefore, we might expect the fitness of these females to be comparable to those of homotypic crosses, although this fate remains unexplored.
Concluding remarks and outlook. Many asexual hybrids are reproductively dependent on sexual partners and can, therefore, be viewed as sexual parasites 3 . Despite the hundreds of crossing experiments investigating the inheritance of hybrids conducted in the past, the reproductive strategies of hemiclonal hybrids are not fully understood. A characteristic feature of hybridogenetic taxa is that one parental genome is perpetuated by hybrids in a clonal form, while the other is recombined in each generation by the Mendelian parental species. www.nature.com/scientificreports/ Whether the ability of P. esculentus males to form new sexual females is beneficial or rather disadvantageous for the hybrids‚ as the females are produced at the expense of hybrid genotypes, remains unknown. We hypothesize, however, that this may represent a way for hybrid males to increase the number of potential female partners in the population, and thus, a so-far unexplored strategy of asexual vertebrates to maintain themselves. Our study adds another aspect to the phenomenon related to sexual populations. If standard, a mechanism affecting the R genome that is captured, made clonal, then returned to the sexual population, may theoretically slow down a recombination-driven evolution of parasitized P. ridibundus populations in comparison to pure sexual populations with fully Mendelian inheritance. This situation suggests sexual parasites have outstanding potential to change the evolutionary trajectory of their host species.

Material and methods
Sampling and crossing experiment. A total of 26 adult water frogs that originated from three mixed R-E and two pure P. ridibundus populations were included in this study (Table 3; Supplementary Table S2). Each of 16 P. esculentus males was crossed with two or three of nine P. ridibundus females to explore their gamete production (Supplementary Table S1). One P. ridibundus male, as a control, was crossed with all females. Crossing experiments followed Berger et al. 46 with modifications following Pruvost et al. 47 . Fully metamorphosed froglets were anesthetized in a 2 mg/l MS-222 solution and dissected for a morphological sex determination.  (25-, 27-, 30-, 39-, and 41-2013), which differed in the genome clonally passed down by the parents. We scanned genomes for SNPs data of both parental genomes. ddRADseq libraries were prepared following the protocol of Peterson et al. 49 with modifications by Brelsford et al. 50 . Briefly, template DNA for each sample at equimolar concentrations (approx. 15 ng/μl) was digested using a combination of MspI and SbfI restriction enzymes. P1 and P2 adapters containing barcodes were then ligated‚ and these ligated fragments were amplified in four PCR replicates of 20 cycles, lowering the risk of polymerase-induced errors. The products were then pooled, and fragments of 350-500 bp were selected from 2% agarose gels, of which 135-139 bp belonged to the pair of Illumina primers, adaptors, and barcodes. The pool of 96 individual libraries was sequenced with a single-end run of 100 cycles (for 150 bp) on a HiSeq 2500 Illumina instrument at the Genomic Technologies Facility at the University of Lausanne, Switzerland.

Analyses of genome inheritance with phasing variants in hybrids.
The microsatellite alleles were scored with GeneMapper v. 3. 7 (Applied Biosystems, Zug, Switzerland). Raw genotypes of sexual species were checked for genotyping errors and null alleles with Micro-Checker v. 2.2.3 51 . Hybrid genotypes were converted to separate lessonae-and ridibundus-specific alleles, following Doležálková-Kaštánková et al. 35 . To distinguish clonally inherited and recombined genomes in backcrossed B1 individuals, we ran the analysis to separate multilocus genotypes (MLGs) using GenAlEx v. 6. 41 52 . To reveal potential gene flow between L and R genomes, we analyzed allele frequencies, heterozygosity, and polymorphism using GenAlEx v. 6. 41 52 . We performed PCoA in GenAlEx v. 6. 41 52 and, in parallel, PCA in radiator v 0.0.18 R package 53 with data standardization, which was based on the total variation of microsatellite allele frequencies. MLGs with more than 65% of missing data were removed from the analysis. ddRADseq reads were demultiplexed with the STACKS v 1.42 module 54 process_radtags (-renz_1 sbfI -renz_2 mspI -E phred33 -r -c -q -i fastq -adapter_1 ACA CTC TTT CCC TAC ACG ACG CTC TTC CGA TCT -adapter_2 GTG ACT GGA GTT CAG ACG TGT GCT CTT CCG ATCT -adapter_mm 2). To clip the adaptor sequence fragments, if present, we applied a custom script (Supplementary File S1). We generated 105,780,761 single-end reads of average read quality > 35, and after applying processing with the Stacks pipeline, 15,676 SNPs remained. Of these, 2265 were informative and were thus retained for the subsequent analyses. The Pelophylax lessonae draft genome assembly (3,453,250 scaffolds) was prepared for mapping purposes with a custom awk script, connecting series of 450 sequential scaffolds (> 190 bp) into an artificial sequence spaced with the fragments of 300 N's. Such an artificially stitched reference genome contained 4,605 scaffolds that were used for further mapping. The reads were mapped on a draft genome assembly of P. lessonae with bowtie2 (-local -very-sensitive), the mapped data were filtered with samtools v 0.1.19 for the mapping quality > 20. Further, we constructed a common catalog in Stacks v 1.42 using the ref_map.pl pipeline. Genotypes were outputted with the Stacks module populations treating each family separately by specifying a family-specific "population map " file and keeping parameters -p = 2 and -r = 0.95, having aimed to retain only the loci that were present in both parents and > 95% of the offspring in a given family. Each population catalog was whitelisted for the biallelic loci only. Ultimately, the five whitelists were merged into one (comprising 4,755 loci), and the population's module was rerun on the initial catalog for all five families together.
So far, known scripts phase parent-progeny duos and trios on data coming from model Mendelian species like humans and imply imputation, i.e., the probability of the mistake is calculated for a missing allele in a genotype, and, if it does not follow the Mendelian pattern, it is considered as an artifact of the data. We designed a script www.nature.com/scientificreports/ that separates non-Mendelian hybrid profiles to haploid ones according to species-specificity of alleles or SNP sites in parental genomes. Using the 15,676 SNP data set for five families, the R genomes were separated from L genomes using custom python scripts. Loci were phased based on their segregation patterns among mother, father, and offspring, separately for each family. For example, an R specific allele will be heterozygous in the father, homozygous in the mother, and heterozygous in all offspring of an R[L] x RR cross. The [L] allele, in this case, is the one which is present only in the father and the offspring, but not in the mother. Similar tests were performed and tailored specifically for all cross types. Such tests can be sensitive to data artifacts. For example, null alleles in one parent, so that the progeny would be heterozygous (null-allele/allele), while appearing as a homozygous. These artifacts were checked for in our scripts, and where appropriate, these loci were removed. For a detailed explanation of this phasing approach and all code used, see https:// github. com/ DanJe ffries/ Hybri dogen_ paper.

Analyses of the structure and relatedness between clonal and sexual haploid genomes. STRU
CTU RE analysis (v 2.3.4) 55 with K = 2 was made by running a chain of 30,000 with burn-in 10,000 under the admixture model with correlated allele frequencies of microsatellites. The alpha parameter was set to 0.5 following recommendations of Wang 56 , who showed the importance of setting the initial value of this prior being equal to 1/K for improved accuracy while analyzing unbalanced populations (as in our case, we had RR, RL, but not LL individuals). STRU CTU RE analysis of the phased dataset was run with the same settings. PCA on diploid and phased SNP datasets was performed with the Pegas (Paradis 2010), radiator 53 , and adegenet 57 R packages. To explore relationships between phased haplotypes, we constructed SNP dendrogram (additive branch length) with the R package SNPRelate v1. 16.0 58 using a set of 2219 SNPs, i.e., without filtering on presence/absence of loci to account for null allele variation between the samples (Supplementary Figure S4).
To investigate the traces of possible recombination within the germline of the RL males, we reconstructed the sex-wise linkage maps using diploid genotypes of the parents and progeny. In sexual females, which recombine, loci were expected to distribute along with the linkage groups; however, in hemiclonal males, the distance between markers was expected to be 0 cM. We reconstructed linkage maps using programs Lep-Map3 59 and MSTmap 60,61 . In order to assess coancestry and genomic variation among phased R-sexual and [R]-/[L]clonal haplotypes, we used the fineRADstructure pipeline 62 of the fineSTRU CTU RE Markov chain Monte Carlo (MCMC) clustering algorithm 63 . This method takes into account linkage information between loci and generates a coancestry matrix accounting for the closest neighboring haplotype, the data undergoes MCMC sampling, and ultimately closely related genetic clusters are inferred. Phased data in the form of a VCF file were converted into fineRADstructure input with a genomic converter module of radiator v 0.0.18 R package 53 . Since our dataset was mapped on a draft genome assembly that does not reflect physical linkage to a large extent, we corrected for linkage disequilibrium with the sampleLD R-script following the recommendations of fineRADstructure developers.
Ethical approval. All experimental protocols were approved by the Federal Commission for Animal Experiments (FCAE) and by the Standing committee on animal health under the Federal Food Safety and Veterinary Office FSVO, Switzerland, under permit nos 119/2013: TV 5113 and TH 103. We declare that all other manipulations with animals were performed in accordance with relevant guidelines and regulations. The study was carried out in compliance with the ARRIVE guidelines.

Data availability
All the microsatellite data supporting our results and conclusions, along with sufficient details, are included in the manuscript and Supplementary files. The initial P. lessonae genome assembly is available at https:// blast. natur kunde museum. berlin/ frog/. ddRADseq libraries are accessible under the BioProject ID PRJNA662869.
Received: 28 April 2020; Accepted: 4 January 2021 www.nature.com/scientificreports/ Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.