Thirteen P. lessonae, 69 P. ridibundus and 83 P. esculentus males and females were collected at 23 sample sites (Supplementary Table 1). Frogs were taxon-determined by phenotypic characters [18, 19] (Fig. 1A–C). DNA was extracted from an interdigital forelegs membrane using a commercial Tissue DNA Isolation Kit (Geneaid Biotech, Taipei, Taiwan) following a manufacturer protocol.

Fig. 1.

Morphological variation (shape of the metatarsal tubercle, dorsal/lateral/ventral colouration) in three Pelophylax taxa and examples of early development in progeny. (A) P. ridibundus (RR). (B) P. lessonae (LL). (C) P. esculentus (RL). (D) Stages of regular early development in progeny from RR $\times$ RL crosses: (a) morula. (b) 2cell morula. (c) 4cell morula. (d) 8cell morula. (e–g) morula. (h) gastrula. (i) neurulation. (j,k) neurula. (l–n) larvae. (E) Variants of embryonic/larval malformations detected in RL $\times$ RR crosses: (a–d) unequal divisions during embryogenesis. (e–h) dying apoptotic embryos in later phases. (i–l) malformed neurulas/larvae.

The experimental procedures followed directives of the State Veterinary Administration of the Czech Republic under the Ethical Committee of the Faculty of Science, Charles University, Prague, permit number 34711/2010-30 issued by the Ministry of Agriculture of the Czech Republic. Frogs were collected under permit no. 358/2011, 278/2011 provided by the Agency for Nature Conservation and Landscape Protection of the Czech Republic. The permits for the crossing experiments were obtained from the Swiss authorities (experimental permit 119/2013: TV 5113 and TH 103).

To explore gamete production and early development in water frogs, we made 211 crosses, including 63 P. esculentus males from the R-E populations, six hybrid males and 14 females from the L-E populations, 22 P. ridibundus males and 47 females, eight P. lessonae males and five females. The schema of crossing design is given in Supplementary Fig. 1 and Supplementary Table 2.

The artificial fertilization procedure followed the methodology from Berger et al. [20] with some modifications. Males were euthanized in a buffered (pH = 7.0) 2 mg/L MS-222 solution (Sigma A-5040, St. Gallen, Switzerland); their testes were removed and stored in a Petri dish with Holtfreter’s solution (pH = 7.4) before use. Females were triggered to ovulate by an injection of salmon luteinizing-releasing hormone (LHRH, Sigma L4897, Prague, Czech Republic). For this purpose, 2 mg of the hormone were diluted in 100 mL of Holtfreter’s solution; per 10 g of body mass, 0.1 mL of this solution was injected into the abdominal cavity. After 16–18 hours, ovulation was checked by pressing the female belly carefully between the thumb and the index finger of the left hand and opening the cloaca with curved forceps. Ovulation was indicated when some eggs were released from the cloaca. Females that did not ovulate received a second hormone injection.

Testes were sliced and crushed to release sperm in a new Petri dish containing aged tap water. Eggs from one female were gently stripped into the sperm solutions and covered with aged tap water. The fertilization success was indicated by egg rotation that turned the black animal pole to the top within 10–40 min after fertilization. On the second day, eggs were checked under a microscope for the presence of the fertilization membrane, indicating successful penetration of a sperm cell into the egg. All eggs were photographed and transferred to 1.5 L plastic boxes (20 $\times$ 12 $\times$ 7 cm), which contained 0.8 L of aged tap water. Unfertilized eggs, egg yolk, and aborted embryos were carefully removed on the third day of the experiment and repeated every second day to avoid bacterial and fungal infections causing embryonic mortality (EM). About 12 days later, when the embryos started to reach the free-swimming stage 25 [21], the boxes were photographed, and tadpoles were transferred to 500 L outdoor containers filled with pond water and steamed straw. Crosses with more than 70 tadpoles were evenly distributed into two containers to reduce inter-individual competition during development. Tadpole mortality (TM) was estimated from the number of tadpoles that died after stage 25. About 50 days later, tadpoles started to metamorphose. Fully metamorphosed froglets with completely reabsorbed tails were anaesthetized in 2 mg/L MS-222 solution and dissected or reared till the maturity (2–4 years) to be crossed. Their sex was determined through the morphology of gonads. The juveniles with malformed or insufficiently developed gonads were marked as “0” sex. All gonads were photographed.

For the following analysis, we divided crosses into five groups that differed in mother/father genotype combinations according to mating pairs known from natural water frog populations: R-E (female $\times$ male; RR $\times$ RL), L-E (LL $\times$ RL, RL $\times$ RL), primary hybridization (LL $\times$ RR, LL $\times$ RR), homospecific parental crosses (RR $\times$ RR, LL $\times$ LL) and other crosses (RR $\times$ RL from L-E, RL $\times$ RR, RL from L-E $\times$ RL from R-E).

We analysed data in R 4.1.2 [22] using primarily Generalized Estimating Equations Models (GEEGLM) with a binomial distribution (logistic link) from the library “geepack” [23]. As a dependent variable, we used (i) presence/absence of egg fertilization, (ii) presence/absence of embryos survival till stage 25, (iii) presence/absence of tadpole survival till metamorphosis. As one cluster in the model, we defined one clutch of eggs from one crossing. For the analysis, we selected only combinations of males and females with more than six cross IDs. In order to test, which factor best explains the successful development rate of embryos, we used as possible explanatory variables: female genotype, male genotype, mitochondrial information of mother, compatibility of mitochondrial and nuclear information, population type of female, population type of male, to create set of models. We compared these models by calculation of quasi-likelihood under the independence model criterion (QIC) for GEEGLM from the library “MuMIn” [24]. We performed the stepwise selection up to the third step, including also the interaction between tested variables. Models were compared based on QIC also to a null model.

We then used GEEGLMs and a set of explanatory variables useful to test particular hypotheses: (i) Are there differences in the proportion of fertilized eggs/embryos survival/tadpole survival among various types of crosses (L-E system, R-E system, primary hybridization, homospecific parental crosses)? (ii) Which gametotype (father’s/mother’s) influence the proportion of embryos survival among various type of crosses? (iii) Are there differences in embryonic survival among particular progeny in interaction with mitochondrial information? (iv) Are there differences in embryonic survival among individual mother $\times$ father combinations? (v) Are there differences in the proportion of embryo/tadpole survival between families fathered by RL males from L-E and R-E systems? In addition, we compared the fertility of females by the generalized linear mixed model (GLMM) with a negative binomial distribution (logarithmic link) from library “lme4” [25]. We used the number of eggs of each female as a dependent variable, the female's genotype as an explanatory variable, and her population type and mitochondrial information as random terms. The significance of the model was tested by comparison to a null model. The plots were created using library “sciplot” [26].

The type of mitochondrion may play a role in a frog sensitivity to oxygen deficiency and its survival [27, 28]. We, therefore, amplified and sequenced the mitochondrial ND2 gene following Plötner et al. [29]. PCR products were commercially Sanger-sequenced by SeqMe s.r.o. (Prague, Czech Republic). Sequence editing was performed and aligned in BioEdit v.7.0.9.0 (Bioedit Company, Raleigh, USA) [30], and variable sites in sequences were evaluated using Mega v 5.1. (Tempe, Arizona) [31].

Parentals were genotyped at ten microsatellite loci using two multiplex PCR sets. Multiplex 1: RlCA1b5 [32], Ga1a19 [33], Rrid013A [32, 34], Res14 [35]. Multiplex 2: Res22 [35], Rrid169A [36], Re1Caga10, Re2Caga3 and RlCA1b6 [33], Rrid082A [36]. PCR protocol was based on a study by Christiansen and Reyer [36]. Individual genotypes were based on species-specifity of amplified alleles described in Doležálková-Kaštánková et al. [16]. Fragment-length analyses were performed on the ABI 3730 Avant capillary sequencer (Applied Biosystems, Foster City, California, USA) with an internal size standard (GeneScan-500 LIZ, Thermo Fisher Scientific, Waltham, MA, USA); the alleles were scored with GeneMapper v. 3. 7 (Applied Biosystems, Zug, Switzerland).

To verify the preliminary phenotypic determination of the adult frogs, we run Principal Coordinate Analysis (PCoA) via covariance matrix with data standardization in GenAlEx v. 6. 41 (Canberra, Australia) [37]. PCoA was performed on multi-allelic microsatellite profiles of adult frogs used as parents in crosses.

Results of Principal Coordinate Analysis (PCoA) on microsatellite data from 143 adult individuals supported morphological taxon identification (Supplementary Fig. 2) of 81 P. esculentus (12 females, 69 males), 52 P. ridibundus (42 females, ten males) and 10 P. lessonae (four females, six males).

The cytoplasmatic background of the females was tested for the possible influence of the mito-nuclear compatibility on the type of produced gametes. Analysis of mitochondrial variation based on 55 ND2 sequences (partial length 669 bp) detected 86 variable sites (Supplementary Table 4) corresponding to three species-specific mtDNA profiles when compared with the NCBI BLAST results (GenBank numbers: MN864876, MN808439, AM749716); P. ridibundus mtDNA was found in 15 RR females; P. lessonae mtDNA contained 24 RR females, 10 RL females and four LL females; P. kurtmuelleri (KK) mtDNA was found in two RR females (Supplementary Table 3).

To begin with, we studied differences among hybrid male gonads in terms of functional development and the ability to reproduce. Pelophylax esculentus males differed in testis size and shape (Supplementary Fig. 3), with the right testis usually being smaller. Moreover, we observed atypical morphological structures with evidence of segmentation (Supplementary Fig. 3).

Experimental crosses between RR, RL and LL males and females resulted in 211 families (Supplementary Fig. 1), with 2939 juveniles undergoing metamorphosis (Supplementary Tables 5,6). The fertilization success (FS) in RR $\times$ RL (R-E) backcrosses ranged from 39.6% to 100%. The only case of 100% FS was observed in a cross 20-GH fathered by an R-E hybrid male 12-2016 AD5. Generally, higher values of FS were recorded in homospecific parental crosses RR $\times$ RR (65.1–99.9%) and RR $\times$ LL crosses (87.4–100%), indicating a better fertilization ability in recombined sexual sperms. On the other hand, one of the lowest values of FS (2.9 and 4.7%) was observed in the eggs of two P. lessonae females and hybrid males from the L-E system (crosses 14-BC and 23-BC). The same hybrid male (15SK2WF3M) had low FS (4.1%) in combination with another P. ridibundus female (cross 9-BC), suggesting a low quality of sperm. Detailed information about FS is given in Supplementary Table 5, Fig. 2a.

Fig. 2.

Proportion of (a) unfertilized eggs, (b) undeveloped embryos, (c) undeveloped tadpoles, (d) overall survival in relation to progeny genotype at different population systems (mean $\mathrm{\pm }$ SE). Detailed informations (n, clusters, df, $\chi$${}^2$, p-value, QIC${}_{\mathrm{\Delta }\text{null.m}}$, W) are described in the main text, section 3.5.

Female fecundity based on the quality and quantity of ovulating eggs was higher in RR females in comparison with RL and LL females (df = 57, $\chi$${}^2$ = 12.62, p = 0.002, QIC${}_{\mathrm{\Delta }\text{null.m}}$ = –8.60). The biggest egg clutch (4328 eggs) was found in P. ridibundus female (16CZ2WF2F), finally crossed with seven hybrid males and FS ranging from 50.2 to 86.1% (Supplementary Table 5).

For further analysis, we used families with fertile females (FS was higher than 5% and EM was lower than 90%). Within 192 families (Supplementary Table 5), there were 43028 fertilized eggs and 15417 embryos reaching stage 25, and 2939 tadpoles completing metamorphosis (Supplementary Table 6). Although most eggs cleaved, many zygotes did not get beyond the blastula stage.

R-E crosses: The number of RR $\times$ RL embryos decreased from 24708 to 5068 by stage25 during their development. The EM ranged between 4.5–100%, with an average of 79.5% (see Fig. 2, Supplementary Tables 5,6). The progeny fathered by 19 R-E males (Albrechtičky, Košatka) had an EM of 100% (24 families). In four R-E males (13cz3wf35m, 16cz4wf1m, 12-2016 AD1, 19cz1wf6m), we did not observe any developing fertilized eggs even in combinations with fertile females, implying that they were sterile. The average TM was 66.7% in total, with 1690 metamorphosed froglets (Fig. 2, Supplementary Table 6). The maximum TM of 100 % was observed in 18 families fathered by 14 R-E males (Košatka, Lébus, Karsibor). A comparison of EM and TM in these families detected a higher average EM (79.5%) than average TM (66.7%). Detailed information is given in Supplementary Table 5, Supplementary Fig. 4 and Fig. 2.

L-E crosses: Embryonic survival rates differed significantly in families with P. esculentus backcrossed to P. lessonae. Lower EM (35.4%) was observed in progeny mothered by hybrid females than in progeny fathered by hybrid males (EM 79.2%), see Supplementary Table 6, while tadpole mortality was similar (75.5/79.7%). Hybrid x hybrid crosses (resulting in RR progeny) showed extreme mortality rates (EM 99.1%, TM 100%), see Fig. 2, Supplementary Fig. 4.

Homospecific parental crosses and primary hybridizations: In parental combinations (RR $\times$ RR, LL $\times$ LL), the average EM was much smaller (36.9%) than in hybrid backcrosses (RR $\times$ RL, RL $\times$ RR) (58.0%), resulting in RR progeny. The survival rate in F1 hybrid embryos (RR $\times$ LL) was higher in families mothered by P. ridibundus (low EM 15.6%) than in reciprocal families mothered by P. lessonae (high EM 63.8%), see Supplementary Table 5. A high EM (78.4%) was detected in hybrid $\times$ hybrid crosses that resulted in RL progeny (up to 100% EM in the family 11-HI). Eight families (two RR $\times$ RR, three RR $\times$ RL and three RL $\times$ RL) exhibited 100% TM. The EM and TM in 11 homospecific RR $\times$ RR crosses showed a higher average TM (93.6%) than average EM (35.2%), Fig. 2.

There were significant differences among types of crosses (Fig. 2) in proportion of fertilized eggs (n = 48402, clusters = 166, df = 3, $\chi$${}^2$ = 590, p 0.001, QIC${}_{\mathrm{\Delta }\text{null.m}}$ = –1120), developed embryos (n = 37416, clusters = 166, df = 3, $\chi$${}^2$ = 117, p 0.001, QIC${}_{\mathrm{\Delta }\text{null.m}}$ = –4009), developed tadpoles (n = 7111, clusters = 112, df = 3, $\chi$${}^2$ = 124, p 0.001, QIC${}_{\mathrm{\Delta }\text{null.m}}$ = –267) and in consequence in overall survival (n = 41736, clusters = 144, df = 3, $\chi$${}^2$ = 215, p 0.001, QIC${}_{\mathrm{\Delta }\text{null.m}}$ = –356).

Proportion of fertilized eggs differed among types of crosses significantly (n = 48402, clusters = 166, df = 3, $\chi$${}^2$ = 590, p 0.001, QIC${}_{\mathrm{\Delta }\text{null.m}}$ = –1120), with RR (W = 1.27, p 0.001), LE system (W = 21.70, p 0.001), and RE system (W = 36.00, p 0.001) having higher mortality than LL. Proportion of developed embryos differed among types of crosses significantly (n = 37416, clusters = 166, df = 3, $\chi$${}^2$ = 117, p 0.001, QIC${}_{\mathrm{\Delta }\text{null.m}}$ = –4009), with RR not differing (W = 0.32, p = 0.570), but with LE system (W = 1.96, p 0.001) and RE system (W = 9.74, p 0.001) having higher mortality than LL. Proportion of developed tadpoles differed among types of crosses significantly (n = 7111, clusters = 112, df = 3, $\chi$${}^2$ = 124, p 0.001, QIC${}_{\mathrm{\Delta }\text{null.m}}$ = –267), with RR (W = 60.10, p 0.001), LE system (W = 29.60, p 0.001), and RE system (W = 34.50, p 0.001) having higher mortality than LL. Overall survival differed among types of crosses significantly (n = 41736, clusters = 144, df = 3, $\chi$${}^2$ = 215, p 0.001, QIC${}_{\mathrm{\Delta }\text{null.m}}$ = –356), with RR (W = 4.92, p 0.001), LE system (W = 7.35, p 0.001), and RE system (W = 9.24, p 0.001) having higher mortality/lower survival than LL.

Proportion of developed embryos (Fig. 3) depends primarily on male gametotype (n = 41763, clusters = 188, df = 2, $\chi$${}^2$ = 77.10, p 0.001, QIC${}_{\mathrm{\Delta }\text{null.m}}$ = –8447), and female genotype (n = 41763, clusters = 188, df = 2, $\chi$${}^2$ = 8.90, p = 0.012, QIC${}_{\mathrm{\Delta }\text{null.m}}$ = –258), and on their interaction (n = 41763, clusters = 188, df = 6, $\chi$${}^2$ = 94.70, p 0.001, QIC${}_{\mathrm{\Delta }\text{null.m}}$ = –8963).

Fig. 3.

Proportion of undeveloped embryos (EM) in relation to father genotypes (mean $\mathrm{\pm }$ SE). Detailed information (n, clusters, df, $\chi$${}^2$, p-value, QIC${}_{\mathrm{\Delta }\text{null.m}}$) are described in the main text, section 3.5.

Proportion of developed embryos (Fig. 4) differed among expected progeny types (n = 37644, clusters = 167, df = 3, ${\chi }^2$ = 102.20, p 0.001, QIC${}_{\mathrm{\Delta }\text{null.m}}$ = –6787), and among mitochondrial types (n = 37644, clusters = 167, df = 2, ${\chi }^2$ = 33.70, p 0.001, QIC${}_{\mathrm{\Delta }\text{null.m}}$ = –2950). RL progeny (W = 0.62, p = 0.430) and RR progeny (W = 0.11, p = 0.740) had similar mortality of embryos as LL progeny, while RR/RL progeny had much higher mortality (W = 65.26, p 0.001). Progeny with mitochondria from Central Europe, LL-like mtDNA (W = 31.57, p 0.001) and RR-like mtDNA (W = 27.12, p 0.001) had lower mortality than progeny bearing the South-European KK-like mtDNA.

Fig. 4.

Proportion of undeveloped embryos concerning their genotype and type of inherited mitochondrion (mean $\mathrm{\pm }$ SE). Detailed information (n, clusters, df, $\chi$${}^2$, p-value, QIC${}_{\mathrm{\Delta }\text{null.m}}$, W) are described in the main text, section 3.5.

Proportion of developed embryos (Fig. 5) differed among mother $\times$ father combinations (n = 37976, clusters = 161, df = 7, $\chi$${}^2$ = 408.00, p 0.001, QIC${}_{\mathrm{\Delta }\text{null.m}}$ = –9029). Homospecific crosses RR $\times$ RR had similar mortality as LL $\times$ LL crosses (W = 0.32, p = 0.571), from primary hybridizations, LL $\times$ RR had similar mortality (W = 2.39, p = 0.122) and RR $\times$ LL had lower mortality (W = 4.21, p = 0.040), from LE system RL $\times$ LL had similar mortality (W 0.01, p = 0.981) but RL $\times$ RL had higher mortality (W = 177.20, p 0.001) as well as LL $\times$ RL combination (W = 4.68, p = 0.030), and RE system represented by RR $\times$ RL combination had the higher mortality (W = 39.47, p 0.001). Thus, the model can be simplified on four groups: homospecific parental crosses, primary hybridization, LE system and RE system. There were significant differences among these groups (n = 37976, clusters = 161, df = 3, $\chi$${}^2$ = 62.10, p 0.001, QIC${}_{\mathrm{\Delta }\text{null.m}}$ = –7721), however this simplification leaded to worse QIC (QIC$\mathrm{\Delta }$ among models = +1308). In comparison to homospecific parental crosses, primary hybridizations had similar mortality (W = 1.70, p = 0.193), but LE system tendet to have higher mortality (W = 3.10, p = 0.078) and RE system had much higher mortality (W = 41.89, p 0.001).

Fig. 5.

Proportion of undeveloped embryos in relation to mother x father combinations (mean $\mathrm{\pm }$ SE). Detailed information (n, clusters, df, $\chi$${}^2$, p-value, QIC${}_{\mathrm{\Delta }\text{null.m}}$, W) are described in the main text, section 3.5.

Proportion of developed embryos (Supplementary Fig. 4) did not differ between RL males from LE and RE systems significantly (n = 28836, clusters = 135, df = 1, $\chi$${}^2$ = 1.84, p = 0.175, QIC${}_{\mathrm{\Delta }\text{null.m}}$ = +42.95), and did not differ among females (n = 28836, clusters = 135, df = 2, $\chi$${}^2$ = 1.52, p = 0.468, QIC${}_{\mathrm{\Delta }\text{null.m}}$ = +162.60), but there was strong interaction (n = 28836, clusters = 135, df = 4, $\chi$${}^2$ = 107.00, p 0.001, QIC${}_{\mathrm{\Delta }\text{null.m}}$ = –385) with LE males crossing with LL females significantly being the group with highest mortality (W = 13.27, p 0.001). Proportion of developed tadpoles differed between RL males from LE and RE systems significantly (n = 6248, clusters = 92, df = 1, $\chi$${}^2$ = 4.17, p = 0.041, QIC${}_{\mathrm{\Delta }\text{null.m}}$ = –33.70), with males from LE system with higher mortality, did not differ among females (n = 6248, clusters = 92, df = 2, $\chi$${}^2$ = 1.83, p = 0.400, QIC${}_{\mathrm{\Delta }\text{null.m}}$ = +4.62), and there was strong interaction (n = 6248, clusters = 92, df = 4, $\chi$${}^2$ = 2046.00, p 0.001, QIC${}_{\mathrm{\Delta }\text{null.m}}$ = –268) with RE males crossing with RL females significantly being the group with lowest mortality (W = 1.74, p = 0.187).

During the initial stages of embryogenesis, we captured phases of blastula cell divisions, gastrulation and neurulation. We sorted them sequentially (Fig. 1D). In some cases (especially RR $\times$ RL crosses), we registered atypical morphology as malformations or heterogeneous cell divisions (Fig. 1E) that led to embryo death (apoptosis). In later phases (after neurulation), larval malformations occurred not only in RR $\times$ RL crosses but also in other RR $\times$ RR or RL $\times$ RL crosses.

Survival rates in pure parental species pairs. Overall, the best values in the proportion of developing eggs, embryos, tadpoles and overall survival showed the progeny of homospecific LL $\times$ LL crosses (Fig. 2). This species, P. lessonae, typically occurs in low numbers within commonplace mixed populations with P. esculentus [15]. Better performance in early developmental stages may give P. lessonae a competing advantage in survival rates against the reproductive dominance of parasitizing P. esculentus that is not only higher than P. lessonae in numbers but also in clutch sizes (Supplementary Table 5). Still, a relatively high proportion of fertilized eggs and developed embryos were observed in RRxRR crosses. However, the proportion of undeveloped tadpoles was even higher than the values counted for L-E and R-E systems with included hybrids (Fig. 2).

We assume that later when the embryo finishes neurulation, and the arising tadpole reaches stage 25, ecological factors may impact TM more than genetic ones. Omitting the biotic factors like infections observed in some semi-natural experiments [38, 39, 40], factors such as competition or density of tadpoles may influence the amount of successfully metamorphosed juveniles of amphibian taxa. Reading and Clarke [41] found positively correlated evidence of proportionally higher tadpole mortality in higher tadpole densities. The negative effect of density in other amphibians delayed tadpole growth due to the excretion of specific substances by tadpoles [42] and resulted in smaller sizes at the metamorphosis stage [43]. Toxic defences were observed in Bufo toads [44], Litoria [45] and other Rana tadpoles [46]. Laboratory-raised tadpoles of P. shqipericus showed high sensitivity to a temperature that strongly affected developmental and survival rates [47]. We do not know the exact causation for the low P. ridibundus tadpole performance constantly observed in our study. If environmental factors are functional players in tadpole performance, the ecological disadvantage of this species will provide positive prospects from a conservation perspective. Pelophylax ridibundus has recently been expanding and invading many areas throughout Europe [48, 49, 50, 51]. Any ecological handicap in relation to P. esculentus or P. lessonae may prevent unwanted population replacements of L-E systems by pure P. ridibundus.

Total proportions of fertilised eggs and developed embryos were relatively similar between L-E and R-E systems but much lower than homospecific crosses in parental taxa. However, once the early stages passed this period, tadpoles mostly of hybrid hemiclonal origin performed even better than pure P. ridibundus progeny. Survival rates were not affected by the hybrid genotype as such but were defined by the parent’s gender. Greater fertility and a relatively high proportion of developed embryos, in values comparative to LL and RR females from homospecific crosses, were shown by hybrid females backcrossed to P. lessonae males. Their average EM of 35% significantly differed from reciprocal crosses between hybrid males and P. lessonae females, with an EM of almost 80% (Fig. 5). Even higher EM (99.5%) was recorded in progeny fathered by a hybrid male from Poland [52]. Combining two parental hybrids provided the highest EM over mother x father pairs (Fig. 5). The low performance of embryos in hybrid x hybrid crosses have also been previously observed [53, 54, 55]. However, comparing mother $\times$ father and father $\times$ mother mate pairs (Figs. 3,5), we found that a low proportion of developing embryos results from the hybrid male sex.

Low EM was not only a feature of hybrid males from L-E populations but also of conspecific males from R-E populations. The mortality of progeny fathered by such P. esculentus males was significantly higher in the initial phases of embryogenesis (EM) than during the later phases of tadpole development (TM, Supplementary Fig. 4). Our results correspond with the research of Kawamura and Nishioka [56], who presented mortality rates close to our results in two B1 families of amphispermic hybrid males from Western Germany. Some authors observed the opposite pattern i.e., a lower EM but higher TM [15, 52]. Variation in survival rates in progeny fathered by males from the R-E system may have been caused by variations in sperm production. The occurrence of different germ cells in the testes may lead to the formation of only R sperm, only L sperm, or both sperm types at once, resulting in variable progeny [16, 17, 52, 57].

One theory predicts that high embryonic mortality might be caused by the accumulation of deleterious mutations in clonally transmitted genomes [55]. A study of all-hybrid P. esculentus populations revealed post-zygotic selection against parental genotypes, P. lessonae and P. esculentus, explained by a role of low genetic diversity in clonal gametes and fixation of deleterious mutations [10]. However, this phenomenon alone cannot explain contrasting values of survival rates in early progeny between P. esculentus fathers (low values) and P. esculentus mothers (high values). We assume two alternative phenomena might have played a key role in the early development of hybridogenetic water frogs. Only maternal genes (mRNAs and ribosomes) are expressed during the initial stages of early embryonic development, whereas the expression of paternal genes starts in the late blastula [58]. The subsequent incompatibility between L and R genomes in the blastula stage may result in the interrupted or malformed development of cleaved embryos and high EM. Indeed, we detected various morphological abnormalities in hybrid offsprings such as deformed eggs, unhatched embryos, and tadpoles with irregular body shapes (Fig. 1E), correlating with earlier research [55, 59]. We, therefore, conclude that the high embryonic mortality in the hybrid offspring might have been influenced by the genomic incompatibility between the clonal genome from a father and the sexual genome from a mother. Additionally, previous authors have noted that disturbances and irregularities during hybrid spermatogenesis have led to lower fertility [15, 60]. We have no direct data on male spermatogenesis, but our observations of the resulting reproduction patterns suggest either a quality difference among the sperm produced or a generally low-fertile sperm produced by hybrid P. esculentus.

As well as laboratory-produced mate pairs typical for natural populations within L-E and R-E systems, we crossed P. ridibundus and P. lessonae to form new P. esculentus F1s. The successful development of newly born F1 was also reported by Berger and Günther [52] and Hotz et al. [61]. Here we analyzed reciprocal crosses and found significantly unequal survival rates in progeny in which crosses of RR mothers $\times$ LL fathers had EM of only 15.6%. In contrast, crosses of LL mothers $\times$ RR fathers had EM reaching 63.8% (Fig. 5). The natural origin of F1 hybrid progenies has been assigned to primary hybridizations between bigger-sized P. ridibundus females and smaller-sized P. lessonae males [18, 62, 63] due to assumed behavioural constraints in reciprocal amplexus between the large P. ridibundus male and smaller P. lessonae female [61]. Our study shows that in addition to preferential matings due to behaviour constraints, it is also important which parental species is genetically the father and the mother. Behavioural and genome compatibilities between parental taxa favour primary unidirectional hybridizations between P. ridibundus females and P. lessonae males. Both phenomena might have affected the historical origin of P. esculentus. Still, Pelophylax water frogs are one of few known cases because the successful primary crosses between sexual taxa in other animal systems were in most cases apparently unidirectional with respect to sex rather than reciprocal [64]. Possible mechanisms of significant differences in offspring survival in reciprocal crosses remain unclear to us. We may only hypothesize that P. esculentus F1s with P. lessonae mtDNA may have lower mito-nuclear compatibility and performance in survival rates than F1s with P. ridibundus mtDNA. Wild-caught P. esculentus and even P. ridibundus bearing P. lessonae mtDNA are, however, common in nature [29, 65].

Finally, we observed a significant difference in embryonic survival rates depending on the type of inherited mtDNA in progeny. Due to the matrilinear heredity, assumed preferential primary hybridizations between P. ridibundus females and P. lessonae males in nature formed P. esculentus F1s with a P. ridibundus mtDNA type. However, hybrid males producing R gametes and mating with P. lessonae females may have also resulted in the transfer of P. lessonae cytoplasm to hybrid cells [66]. Indeed, not only P. esculentus but also many P. ridibundus bear P. lessonae type mtDNA. The greater effectiveness of the enzymes encoded by P. lessonae mtDNA may give some advantage to P. ridibundus and probably P. esculentus in the northern parts of their ranges [29]. Pelophylax lessonae is known for better cold tolerance reaching higher latitudes and longitudes, while P. ridibundus is a cold-sensitive species from lowlands that recolonized the deglaciated northern parts of Europe later than P. lessonae [29]. The laboratory progeny of P. ridibundus and P. esculentus showed very similar embryonic survival rates in having P. ridibundus or P. lessonae mtDNA (Fig. 4), supporting the compatibility of P. lessonae cytoplasm with the nuclear DNA background of P. ridibundus genomes collected in Central Europe. A striking contrast can be seen however, in the significantly higher numbers of undeveloped embryos in P. ridibundus and P. esculentus progeny when mtDNA were of the P. kurtmuelleri type (Fig. 4). Although Pelophylax kurtmuelleri mtDNA occurs at low levels in Central European water frogs, it may represent an opposite case to P. lessonae mtDNA, as P. kurtmuelleri is endemic to the warmer Balkans in Southern Europe. Its mtDNA may not be fully compatible with Central European, more cold-adapted, Pelophylax genomes.