Research Article
Research Article
Comparative mitochondrial phylogeography of water frogs (Ranidae: Pelophylax spp.) from the southwestern Balkans
expand article infoPetr Papežík, Peter Mikulíček, Michal Benovics§, Monika Balogová|, Lukáš Choleva#, Marie Doležálková-Kaštánková¤, Petros Lymberakis«, Edvárd Mizsei», Simona Papežíková, Nikos Poulakakis«˄, Enerit Saçdanaku˅, Márton Szabolcs», Radek Šanda¦, Marcel Uhrin|, Jasna Vukićˀ, Daniel Jablonski
‡ Comenius University in Bratislava, Bratislava, Slovakia
§ Masaryk University, Brno, Czech Republic
| P. J. Šafárik University in Košice, Košice, Slovakia
¶ Laboratory of Fish Genetics, Institute of Animal Physiology and Genetics CAS, Liběchov, Czech Republic
# University of Ostrava, Ostrava, Czech Republic
¤ Laboratory of Non-Mendelian Evolution, Institute of Animal Physiology and Genetics CAS, Liběchov, Czech Republic
« University of Crete, Heraklion, Greece
» Institute of Aquatic Ecology, Centre for Ecological Research, Debrecen, Hungary
˄ Institute of Molecular Biology and Biotechnology (IMBB), Foundation for Research and Technology - Hellas (FORTH), Heraklion, Greece
˅ Univeristy of Tirana, Tirana, Albania
¦ Department of Zoology, National Museum, Prague, Czech Republic
ˀ Charles University, Prague, Czech Republic
Open Access


The genus Pelophylax (water frogs) includes relatively common, widely distributed, and even invasive species, but also endemic taxa with small ranges and limited knowledge concerning their ecology and evolution. Among poorly studied species belong endemics of the southwestern Balkans, namely Pelophylax shqipericus, P. epeiroticus and P. kurtmuelleri. In this study, we focused on the genetic variability of these species aiming to reveal their phylogeographic patterns and Quaternary history. We used 1,088 published and newly obtained sequences of the mitochondrial ND2 gene and a variety of analyses, including molecular phylogenetics and dating, historical demography, and species distribution modeling (SDM). We revelated the existence of two mitochondrial lineages within P. epeiroticus and P. shqipericus that diverged at ~ 0.9 Mya and ~ 0.8 Mya, respectively. Contrarily, no deeply diverged lineages were found in P. kurtmuelleri. Pelophylax kurtmuelleri also shows a close phylogenetic relationship with widely distributed P. ridibundus, suggesting that both represent one evolutionary clade called here P. ridibundus/kurtmuelleri. The estimated split between both lineages in the clade P. ridibundus/kurtmuelleri date back to ~ 0.6 Mya. The divergence between the ridibundus and kurtmuelleri lineages on the ND2 gene is thus lower than the divergence between the two lineages found in P. epeiroticus and P. shqipericus. According to haplotype networks, demographic analyses, and SDM, endemic water frogs survived the last glacial maximum (LGM) in Balkan microrefugia, and their distribution has not changed significantly or even retracted since the LGM. Haplotypes of the kurtmuelleri lineage were also found in northern parts of Europe, where haplotype diversity is however much lower than in the Balkans, suggesting the possible hypothesis of their postglacial expansion to the north.


Biogeography, conservation, Eastern Mediterranean, endemism, evolution, Rana, taxonomy


The southwestern part of the Balkan Peninsula is characterized by heterogeneity in topography, habitats, and microclimate (Griffiths and Frogley 2004; Reed et al. 2004; Savić 2008). This region is demarcated by two mountain massifs, the Dinarides and Hellenides, and extends from the Neretva River to Peloponnese (Savić 2008). The region is also characterized by steep mountains and deep valleys that form natural barriers for latitudinal dispersal, and thus, may affect species distribution and diversity (Podnar et al. 2014; Jablonski et al. 2016; Psonis et al. 2018). The current genetic diversity and species richness in the southwestern Balkans are results of the complex paleogeographical history of the region that was accompanied by radiation and speciation of animal taxa during the Miocene and the Pliocene. A series of paleoclimatic fluctuations during the Pleistocene have led to subsequent allopatric differentiation of populations that resulted in high intraspecific diversity (Podnar et al. 2014). The southwestern Balkans is well known as a centre of endemism (e.g., Radea et al. 2017; Sfenthourakis and Hornung 2018; Allegrucci et al. 2021; Mermygkas et al. 2021) including amphibians and reptiles. Endemic species or unique evolutionary lineages are reported in salamanders (Recuero et al. 2014; Pabijan et al. 2017), toads (Fijarczyk et al. 2011; Dufresnes, Mazepa et al. 2019a), frogs (Dufresnes et al. 2013; Jablonski et al. 2021), lizards (Jablonski et al. 2016; Psonis et al. 2017; Kiourtsoglou et al. 2021; Strachinis et al. 2021), or snakes (Guicking et al. 2009; Musilová et al. 2010; Mizsei et al. 2017; Jablonski et al. 2019). The southwestern Balkans is inhabited by three endemic water frog taxa of the genus Pelophylax. Although we have preliminary views regarding their evolutionary relationships (e.g., Lymberakis et al. 2007; Plötner et al. 2012; Hofman et al. 2015), the range-wide intraspecific genetic variability and phylogeography have never been studied.

The Epirus water frog, Pelophylax epeiroticus (Schneider, Sofianidou & Kyriakopoulou-Sklavounou, 1984) is distributed from southern Albania to northwestern Peloponnese (west of the Pindos range; Sofianidou et al. 1987; Sofianidou and Schneider 1989; Schneider and Haxhiu 1992; Oruçi 2008) with the type locality “the eastern shore of Lake Ioannina ca. 5 km south of Amphithea”, Greece (Schneider et al. 1984). Recently, populations of P. epeiroticus were also discovered east of the Pindos Mountains (Strachinis 2021). The species belongs to the ridibundus/bedriagae group that comprises phylogenetically closely related Pelophylax species from the eastern Mediterranean (see Lymberakis et al. 2007 for details). Although its phylogenetic position is still questionable (see Lymberakis et al. 2007 vs. Akın et al. 2010; Hofman et al. 2015), it seems that P. epeiroticus diverged earlier than other taxa within the ridibundus/bedriagae group (Plötner et al. 2012; Hofman et al. 2015).

One of the most endangered European amphibian species, the Albanian pool frog, P. shqipericus (Hotz et al., 1987), is native to the coastal lowlands of Albania and Skadar Lake basin in southern Montenegro, from where it was originally described (Günther 2004; Plötner 2005). Its type locality lies in “Virpazar, Skadarsko Jezero, Crna Gora, Yugoslavia”. This species forms a phylogenetically distinct clade together with P. lessonae and P. bergeri (the lessonae group; Lymberakis et al. 2007; Plötner et al. 2012; Hofman et al. 2015; Vucić et al. 2018).

The last and taxonomically disputed species is the Balkan water frog P. kurtmuelleri (Gayda, 1940). Populations belonging to this species were described as Rana balcanica Schneider & Sinsch, 1992 based on bioacoustics and morphometry (Schneider and Sinsch 1992; Schneider et al. 1993). Rana balcanica is now a junior synonym of P. kurtmuelleri with the type locality “Arsen presso Petrelli” = Erzen river near Petrelë, Albania (Dubois and Ohler 1994). Although this species is considered to be endemic to the Balkans (Valakos et al. 2008), mitochondrial haplotypes and alleles of some genes affiliated with this taxon were recently discovered also in central and eastern Europe (Kolenda et al. 2017; Litvinchuk et al. 2020, Svinin et al. 2021). The species status of P. kurtmuelleri is controversial for a long time (Speybroeck et al. 2010, 2020), mostly due to its low mtDNA divergence and shared interspecific polymorphism with P. ridibundus (e.g., Lymberakis et al. 2007; Akın et al. 2010; Plötner et al. 2012; Hofman et al. 2015; Litvinchuk et al. 2020).

Considering the lack of comprehensive information on the genetic diversity and phylogeography of these three water frog taxa inhabiting the Balkan Peninsula, we decided to analyse a dense mitochondrial dataset to (1) describe their genetic diversity based on the mitochondrial ND2 gene, (2) reconstruct their Quaternary history, (3) compare these data with phylogeographical patterns of other amphibians and reptiles to infer general biogeographical scenarios in the southwestern Balkans, and (4) to deduce implications for conservation genetics and species protection.

Material and methods

Study species and sampling

We obtained a total of 996 samples of Pelophylax spp. from 230 localities across the Balkan Peninsula, with a focus on Albania and Greece (see Table S1 and Fig. S1). We sourced DNA from drops of blood or toe clips of varying ages and sources, which had been stored in 96% ethanol. To identify the species in the field, we used the morphological traits described in Speybroeck et al. (2016) and Papežík et al. (2021). Identification of individuals from areas sympatric for two species was subsequently confirmed through microsatellite loci (for details see Supplementary Material 1).

Lab work and sequence alignments

Genomic DNA was extracted from blood or toe clips using NucleoSpin® Tissue kit (Macherey-Nagel, Düren, Germany) following the manufacturer’s protocol. Complete NADH dehydrogenase 2 (ND2) gene 1,038 bp, was amplified using combinations of primer pairs L1 and H2 (Plötner et al. 2008). The following PCR program was used: 2 minutes denaturing step at 94°C, followed by 35 cycles consisting of denaturing for 30 seconds at 94°C, primer annealing for 20 seconds at 63°C and elongation for 1 minute at 72°C, with a final elongation step for 10 minutes at 72°C (modified from Plötner et al. 2008). Purification of PCR products was performed using the ExoSAP-IT enzymatic clean-up (USB Europe GmbH, Staufen, Germany) following the manufacturer protocol. PCR products were one-directionally sequenced with L1 primer using Macrogen Europe commercial services (Amsterdam, The Netherlands). The newly generated sequences were deposited in GenBank under accession numbers OP794049OP794075 for P. epeiroticus, OP794076OP794097 for P. shqipericus, OP794098OP794108 for P. ridibundus and OP820591OP820697 for P. kurtmuelleri (Table S1).

For genetic analyses, newly obtained sequences were combined with those already deposited in GenBank (Table S1). New sequences were manually aligned, checked, and translated to amino acids to detect possible stop codons using Geneious Prime v. 2020.2.4 (Biomatters, Auckland, New Zealand). Five different datasets were prepared. Three ND2 datasets of P. epeiroticus, P. kurtmuelleri/P. ridibundus, and P. shqipericus were prepared for the evaluation of their intraspecific genetic diversity. The fourth ND2 dataset, which includes the unique haplotypes of studied species, with sequences of other western Palaearctic water frog species (P. perezi, P. saharicus, P. lessonae, P. bergeri, P. cretensis, P. cerigensis and P. cf. bedriagae), and P. plancyi, P. nigromaculatus and Rana graeca as outgroups, was used for reconstruction of their phylogenetic relationships. Finally, the fifth ND2 dataset involves two sequences per phylogenetic lineage of each one of the studied taxa and sequences of the other western Palearctic water frogs was used for the estimation of the divergence times.

Phylogenetic analyses and DNA polymorphism

Phylogenetic relationships were inferred using maximum-likelihood (ML) and Bayesian interference (BI) by RAxML v. 8.1.12 (Stamatakis 2014) and MrBayes v. 3.2 (Ronquist et al. 2012), respectively. PartitionFinder 2 (Lanfear et al. 2017) was used to select the best-fit codon-partitioning schemes and substitution models using a greedy algorithm and the Akaike information criterion (AIC). The best-fit substitution model with each codon position treated separately for the BI and ML analysis was as follows: HKY+G (first position), TrN+I+G (second position), GTR+G (third position) for BI and GTR+G (first position), GTR+I+G (second position), GTR+G (third position) for ML. The best ML tree with the highest likelihood was selected from 20 generated trees. The clade support in the ML was assessed by 1,000 bootstrap pseudoreplicates. The BI analysis was set as follows: two separate runs with four chains in each run, 10 million generations sampled every 100 generations. The convergence of both runs was confirmed by the average standard deviation of split frequencies and potential scale reduction factor diagnostics (P < 0.01). The first 20% of trees were discarded as the burn-in after inspection for stationarity of log-likelihood scores of sampled trees in Tracer v. 1.7 (Rambaut et al. 2018). A consensus tree was constructed from the post-burn-in samples.

Molecular dating was performed using the calibration of the known split of P. cretensis that the most likely occurred at the end of the Messinian salinity crisis (~5.5–5.3 Mya) (Beerli et al. 1996; Lymberakis et al. 2007, Dubey and Dufresnes 2017). Divergence times were estimated using 28 sequences of ND2, representing two sequences per phylogenetic lineage that were divergent more than 1% on p distances of the studied species (P. epeiroticus, P. kurtmuelleri, P. ridibundus, P. shqipericus) and represent thus deeper intraspecific structure. To provide well-resolved phylogeny, selected sequences of the genera Pelophylax and Rana from GenBank were used (see Table S2 for GenBank accessions numbers). BEAST 2 package (Bouckaert et al. 2014) was used to show time-measured phylogeny and its credibility intervals. The above-mentioned calibration point (5.4 Mya; 95% HPD 5.9–4.8) was tested together with the model of nucleotide substitution in PartitionFinder 2 under the AIC. Two distinct runs were carried out, each one with four independent chains of 20 million iterations. Initial two million iterations were excluded as burn-in. Results were inspected using Tracer to check the convergence and effective sample size. For the visualization of the final time-calibrated phylogeny we used TreeAnnotator included in BEAST 2 package.

Inter- and intraspecific genetic structure were also analyzed using Principal Component Analysis (PCA) implemented in ‘adegenet’ package (Jombart et al. 2008) as a part of R statistical environment, v. 4.1.3 (R Core Team 2022).

For a better presentation of the intraspecific evolution, the parsimony network analysis with the algorithm of TCS (Clement et al. 2000) implemented in PopArt v. 1.7 ( was used.

The nucleotide diversity index (π), the number of haplotypes (h), the haplotype diversity (hd), the number of segregating sites (S), and Watterson’s theta per site (θW) were used to estimate the genetic diversity of each species and intraspecific lineages using DnaSP v. 6.00 (Rozas et al. 2017). The same software was used for average values of uncorrected p distances among species and intraspecific lineages.


Demographic history was analysed in BEAST 2, using Bayesian Skyline plots (BSP, Drummond et al. 2005). The strict molecular clock with the mean substitution rate of the initial value 0.016 substitution/site/Mya (based on the molecular dating calibration, see above) was set for each detected taxon or intraspecific lineage. PartitionFinder 2 (Lanfear et al. 2017) was used for the best-fit codon-partitioning schemes selection and substitution models under the (AIC). The best-fitting substitutional model was HKY for all tested alignments, i.e., P. epeiroticus and P. shqipericus, and P. kurtmuelleri and P. ridibundus, both tested separately and together as one evolutionary clade P. ridibundus/kurtmuelleri (see results). Each BSP analysis was run twice to check the consistency between the runs. The number of generations varied significantly among analyzed datasets, with values from 20 million to 300 million generations, and were sampled every 1,000 generations. The first 10% were discarded as a burn-in. The final BSP, effective sample sizes (ESSs) of parameters >200 and convergence of parameter estimates across runs were prepared and checked with Tracer v. 1.7. with the maximum times as the median of the root height parameter. Historical population expansions in the whole dataset and within detected lineages were tested using Tajima’s D (Tajima 1989), Fu’s Fs (Fu 1997) and Ramos-Onsins, Rozas’s R2 statistics (Ramos-Onsins and Rozas 2002), and pairwise mismatch distribution (MSD) analyses (Rogers and Harpending 1992) in DnaSP v. 6.00. The analyses were carried out with 10,000 coalescent simulations of Tajima’s D, Fu’s Fs and Ramos-Onsins and Roza’s R2.

Species distribution modelling

With the combination of literature and our own field data, we altogether collected 101 records for P. epeiroticus (Sofianidou and Schneider 1989; Haxhiu 1994; Beerli et al. 1996; Sofianidou 1996; Guerrini et al. 1997; Mellado et al. 1999; Casola et al. 2004; Oruçi 2008; Akın et al. 2010; Hofman et al. 2015; Sagonas et al. 2019), 53 for P. shqipericus (Hotz et al. 1985, 1987; Hotz and Uzzell 1982, 1983; Haxhiu 1994; Schneider and Haxhiu 1994; Sinsch and Schneider 1996; Guerrini et al. 1997; Spasić-Bosković et al. 1999; Bucci et al. 1999; Ragghianti et al. 1999, 2004; Casola et al. 2004; Akın et al. 2010; Jablonski 2011; Radojičić et al. 2015; Polović and Čađenović 2014; Hofman et al. 2015; Frank et al. 2018; Vucić et al. 2018) and 212 records for Balkan populations of P. kurtmuelleri (own field data). The presence coordinates for each species were filtered to remove potentially inaccurate records. Coordinates outside the Balkan Peninsula (sensu Džukić and Kalezić 2004) were discarded. Records outside of the 90% range of a kernel distribution fitted to the data were also discarded using the adehabitatHR package (Calenge 2006). Finally, a spatially balanced resampling was performed to account for potential spatial bias by iteratively (n = 1,000) selecting records according to the non-overlapping buffer of 0.05° retaining the maximum sample size. To build the SDMs we used the BIOCLIM variables as environmental response variables, obtained from the WorldClim v. 2.1 database for current conditions (derived from meteorological data recorded between 1970–2000) at a resolution of 2.5 arc minutes (~3.5×4.6 km raster cell resolution at the study area) (Fick and Hijmans 2017). For the LGM (the Last Glacial Maximum) past climate conditions, we used the downscaled climate data at the same resolution as for the current climate (2.5’) from simulations with three Global Climate Models (GCMs): (i) CCSM4, (ii) MIROC-ESM and (iii) MPI-ESM-P, which were produced using the bias-corrected WorldClim 1.4 as baseline current climate-based data recorded between ~1960–1990 (Hijmans et al. 2005). Preparation of environmental data and masking to the study area was done using the raster package (Hijmans 2022). Variable selection for each species was performed to reduce multicollinearity and variance inflation in the modelling by calculating the pairwise correlations among the variables, and among each pair of variables correlated above the threshold r>0.7. We excluded the variable with either the highest variance inflation factor (VIF) or the least significant or least informative bivariate relationship with the species presence/pseudo-absence data using the fuzzySim package (Barbosa 2015). We built SDMs for each studied Pelophylax species using the biomod2 package, which allowed us to build ensemble models from several different modelling algorithms, and to project that ensemble model in space and time (Thuiller et al. 2022). We created four presence/pseudo-absence datasets for each species by generating 5,000 random points as pseudo-absences. We used five linear and machine learning model algorithms to reduce uncertainties in different modelling approaches: Generalized Linear Models (GLM), Generalized Boosting Models/Boosted Regression Trees (GBM), Surface Range Envelop (SRE), Flexible Discriminant Analysis (FDA), Maximum Entropy (MaxEnt). We ran four modelling replicates for each presence/pseudo-absence dataset with a data split to 70% training and 30% testing data. To evaluate the individual model replicates we used the True Skills Statistic (TSS) metric, models were included in the ensemble model above TSS>0.7. Models building were done using current environmental conditions and ensemble models were projected for the LGM by projecting ensemble models for each GCM separately. Finally, the LGM projections for each species were combined by calculating the mean of the projections and were set to binary by the threshold of 0.8 suitability index. All steps of generating SDMs were done in R v. 4.1.3.


Phylogeny and phylogeography

Combining newly produced and available GenBank data, we obtained a dataset of 1,088 ND2 sequences (alignment length 1,038 bp) from four water frog taxa: P. epeiroticus (127 sequences), P. kurtmuelleri (780), P. ridibundus (112), and P. shqipericus (69). Both trees, ML and BI have identical topologies regarding the three major clades (P. epeiroticus, P. shqipericus and P. ridibundus/kurtmuelleri) with lnL = –7188.48 and –7744.64, respectively. Phylogenetic trees as well as haplotype networks and PCAs (Figs 1, S1) showed deep differentiation in both P. epeiroticus (genetic distance 2.1%) and P. shqipericus (genetic distance 1.6%) with two intraspecific lineages, named, due to their geographic positioning, Northern and Southern (Figs 2, 3). Moreover, possible third lineage of P. epeiroticus were detected in both phylogenetic network and PCA (Figs 1, 2). In contrast, the genetic differentiation within and between kurtmuelleri and ridibundus lineages (Fig. 4) was low and formed closely related clades in phylogenetic trees. Genetic differentiation of these clades was lower than between the Northern and Southern lineages in P. epeiroticus and P. shqipericus. (Table 1, Fig. 4). Thus, we will further call these three focal clades as P. epeiroticus, P. shqipericus and P. ridibundus/kurtmuelleri.

Table 1.

Summary of DNA polymorphism and demographic statistics for three clades of Pelophylax spp.

N h S π + SD (%) hd + SD θW + SD (%) D p Fs p R2 p
P. epeiroticus 127 31 50 1.143 ± 0.030 0.870 ± 0.019 0.919 ± 0.246 0.681 0.803 –0.574 0.502 0.113 0.818
Northern lineage 54 13 15 0.221 ± 0.019 0.790 ± 0.037 0.328 ± 0.120 –0.988 0.175 –3.853 0.052 0.068 0.142
Southern lineage 73 18 24 0.232 ± 0.036 0.718 ± 0.049 0.492 ± 0.158 –1.629 0.030 –7.692 0.003 0.046 0.036
P. shqipericus 69 25 36 0.912 ± 0.050 0.905 ± 0.026 0.769 ± 0.227 0.607 0.785 –2.231 0.275 0.124 0.814
Northern lineage 42 11 12 0.214 ± 0.022 0.767 ± 0.059 0.278 ± 0.117 –0.702 0.270 –2.870 0.087 0.087 0.262
Southern lineage 27 14 18 0.428 ± 0.031 0.934 ± 0.026 0.466 ± 0.179 –0.290 0.437 –3.818 0.050 0.109 0.368
P. ridibundus/kurtmuelleri 892 141 120 0.631 ± 0.016 0.923 ± 0.005 1.542 ± 0.003 –1.721 0.009 –135.069 0.000 0.027 0.024
ridibundus lineage 112 15 15 0.052 ± 0.009 0.416 ± 0.058 0.299 ± 0.103 –1.984 0.001 –32.132 0.000 0.021 0.006
kurtmuelleri lineage 780 126 102 0.430 ± 0.007 0.911 ± 0.006 1.393 ± 0.278 –2.274 0.000 –17.176 0.000 0.028 0.049
N the number of individuals, h the number of haplotypes, S the number of segregating sites, π the nucleotide diversity, hd the haplotype diversity, θW the Watterson theta, D Tajima’s D statistic, Fs Fu Fs statistic, R2 Ramos-Onsins, Rozas’s R2 statistic
Figure 1. 

Principal component analysis (PCA) of Pelophylax spp. ND2 dataset from southwestern Balkans with the inset PCAs for P. epeiroticus (orange), P. shqipericus (purple), and P. ridibundus/kurtmuelleri (green and blue) showing genetic diversity of selected clades. The oval outlines in PCAs represent 95% confidence intervals. The color scheme of the studied taxa corresponds with the one used in the mitochondrial phylogeny presented in Fig. S1. Insets photographs by Daniel Jablonski.

Figure 2. 

The phylogenetic position, haplotype network, and distribution of detected lineages and haplotypes of Pelophylax epeiroticus. Symbol sizes in the haplotype network reflect sequence frequencies. Small black circles are missing node haplotypes, each line connecting two haplotypes corresponds to one mutation step, if not otherwise indicated the number along the line. The shaded area represents the species range. The colour scheme of the haplotype networks and species ranges corresponds with the one used in PCAs (Fig. 1) and the molecular phylogeny (Fig. S1).

Figure 3. 

The phylogenetic position, haplotype network, and distribution of detected lineages and haplotypes of Pelophylax shqipericus. Symbol sizes in the haplotype network reflect sequence frequencies. Small black circles are missing node haplotypes, each line connecting two haplotypes corresponds to one mutation step, if not otherwise indicated the number along the line. The shaded area represents the species range. The colour scheme of the haplotype networks and species ranges corresponds with the one used in PCAs (Fig. 1) and the molecular phylogeny (Fig. S1).

Figure 4. 

Geographical range of the clade P. ridibundus/kurtmuelleri, haplotype network, distribution of detected lineages, and their haplotypes. For the detailed phylogeographic structure of kurtmuelleri haplotypes see Fig. S2. Symbol sizes in the haplotype network reflect sequence frequencies. Small black circles are missing node haplotypes, each line connecting two haplotypes corresponds to one mutation step. An uncertain contact zone between detected lineages/haplotypes is indicated by question marks. The black arrow represents the affiliation of the introduced population of the P. ridibundus/kurtmuelleri clade in southern Italy based on 517 bp only (MK116532). The colour scheme of haplotype networks and species ranges corresponds with the one used in PCAs (Fig. 1) and the molecular phylogeny (Fig. S1).

According to phylogenetic analyses and molecular dating (Fig. 5), studied taxa of the genus Pelophylax formed two main groups that diverged at 9.8 Mya (11.4–8.2 of 95% HPD). The first group was formed by species P. shqipericus, P. bergeri and P. lessonae, the second group by species P. epeiroticus, P. cretensis, P. bedriagae, P. cypriensis, P. ridibundus/kurtmuelleri. The Northern and Southern lineages in P. shqipericus diverged at 0.8 Mya (1.2–0.5 of 95% HPD), the Northern and Southern lineages in P. epeiroticus at 0.9 Mya (1.3–0.6 of 95% HPD), and the estimated split between kurtmuelleri and ridibundus lineages dated back to 0.6 Mya (0.9–0.4 of 95% HPD).

Figure 5. 

The dated phylogeny of the genus Pelophylax with the intraspecific structure of studied species. Numbers between clades show the estimated time of the divergence (Mya; blue bars indicates 95% highest posterior densities of the estimated node ages). The duration of the Messinian Salinity Crisis (MSC) is highlighted in yellow. The outgroup taxa (P. nigromaculatus and Rana graeca) are not shown. Photos of water frogs represent detected lineages (photographs by Daniel Jablonski).

In P. epeiroticus 31 unique haplotypes in 127 sequences with hd = 0.87 and π = 1.14% (Table 1) were recognized. The Northern lineage is distributed in southern Albania, north-western Greece, down to the Ambracian Gulf (Fig. 2) and it is formed by 13 haplotypes (hd = 0.79 and π = 0.22%; Table 1). The Southern lineage is distributed from the Ioannina Lake, Ambracian Gulf, southwestern Greece, and along the north-western coast of the Peloponnese to Kaiafa Lake (Fig. 2). This lineage is represented by 18 haplotypes (hd = 0.72 and π = 0.23%; Table 1). Uncorrected mean p distance between the lineages is 2.09%.

Pelophylax shqipericus is represented by 69 sequences of ND2 with 25 unique haplotypes, hd = 0.91 and π = 0.92% (Table 1). Its Northern lineage is distributed around Lake Skadar and along the northern Adriatic coast of Albania (Fig. 3) with the southernmost locality detected in Nishaj. This lineage contains 11 haplotypes (hd = 0.77, π = 0.21%; Table 1). The Southern lineage is distributed from Velipojë, south to Orikum where is the southern edge of the species range (Fig. 3), although one individual of this lineage was detected also in Virpazar, Montenegro. The Southern lineage is represented by 14 haplotypes (hd = 0.93 and π = 0.43%; Table 1). Both lineages were found in sympatry at the following localities: Nishaj, Velipojë, Virpazar, and Torovicë. The uncorrected mean p distance between the two lineages is 1.59%.

In the P. ridibundus/kurtmuelleri clade 141 unique haplotypes were identified among 892 sequences, with hd = 0.92 and π = 0.63%. In the kurtmuelleri lineage, we identified 126 haplotypes among 780 sequences (hd = 0.91, π = 0.43%). We detected haplotypes of the kurtmuelleri lineage in Albania, Greece, Northern Macedonia, Bulgaria, Kosovo, Montenegro, Serbia, Croatia, Bosna and Hercegovina, Romania, Ukraine (collected in Kyiv and confirmed by J. Plötner, pers. comm. 2022, GenBank accession number AM900651), Germany, Lithuania, Latvia, Russia, Italy, and France. Only a weak phylogeographic structure was observed in the kurtmuelleri lineage (Figs 4, S2). In contrast, ridibundus lineage showed much lower genetic variability with only 15 haplotypes identified in 112 sequences, and with hd = 0.42 and π = 0.05%. Haplotypes of ridibundus lineage were detected in Greece, Bulgaria, Serbia, Croatia, Romania, Slovakia, Poland, Russia, Italy, and France (Fig. 4). Uncorrected mean p distance between haplotypes representing kurtmuelleri lineage and ridibundus lineage was 1.37%.


Bayesian skyline plot (BSP) for P. epeiroticus (Fig. 6) showed population expansion started ~25 Kya. On the other hand, the mismatch distribution (MSD) and the neutrality tests did not show significant evidence for population growth (Table 1). When demographic analyses were carried out for the Northern and the Southern lineage separately, the Northern lineage of P. epeiroticus showed population growth starting at ~25 Kya according to BSP (Fig. 6) but population stability according to MSD and neutrality tests. In the Southern lineage, BSP (Fig. 6) showed moderate expansion in population size (starting at ~75 Kya) supported also by MSD neutrality tests.

Figure 6. 

Historical demography of Pelophylax spp. distributed in the southwestern Balkans as estimated with Bayesian skyline plots (BSP) and mismatch distributions (MSD). The color scheme of BSPs shading corresponds with the one used in PCAs (Fig. 1) and the molecular phylogeny (Fig. S1). The LGM line with BSP graphs indicates the period of the Last Glacial Maximum.

In P. shqipericus, population expansion based on BSP started approximately at ~25 Kya. MSD and neutrality tests, however, showed no significant population growth (Table 1). In the Northern lineage, BSP showed population expansion starting at ~35 Kya, however, MSD and neutrality test did not corroborate the expansion model. Similarly, BSP of the Southern lineage assumed recent expansion, starting at ~25 Kya, MSD and neutrality tests were non-significant.

Analysing both ridibundus/kurtmuelleri lineages together, the populations started the expansion ~30 Kya. The kurtmuelleri lineage alone started the expansion at ~30 Kya, followed by another steep expansion at ~12 Kya according to BSP (Fig. 6). This result is further supported by MSD and neutrality tests that both showed statistically significant evidence of population expansion (Table 1). The BSP of the ridibundus lineage showed population expansion starting ~17 Kya, which is consistent with MSD, where observed values mirrored expected values in the growing or declining population model (Fig. 6). In addition, also the neutrality test showed significant evidence of population expansion (Table 1).

Current and past distribution models

From the initial 366 occurrence data points (27 unique localities for P. epeiroticus, 20 for P. shqipericus and 198 unique localities for P. kurtmuelleri), 266 (58 for P. epeiroticus, 41 for P. shqipericus, and 167 for P. kurtmuelleri) occurrence data points were included in the final models. From the predictors after the variance inflation factor (VIF) analysis, we ended up with six, seven, and six, respectively, out of 19 initial predictors. For P. epeiroticus these were bio1 (annual mean temperature), bio2 (mean diurnal range – max/min temperature, monthly average), bio3 (isothermality), bio4 (temperature seasonality), bio8 (mean temperature of the wettest quarter), and bio13 (precipitation of wettest period) (Fig. S3). For P. shqipericus bio1, bio3, bio4, bio8, bio9 (mean temperature of the driest quarter), bio14 (precipitation of driest period), and bio15 (precipitation seasonality). For P. kurtmuelleri bio1, bio2, bio8, bio9, bio15, bio19 (precipitation of coldest quarter). The models under current climate conditions had robust evaluations metrics for both P. epeiroticus and P. shqipericus, but in the case of P. kurtmuelleri the algorithms not performed well (Fig. S3). The spatial projection of ensemble models corresponds well with the known range of the studied species (Figs 1, 2, 4). All three models of the LGM provided similar projections with the most suitable areas for P. epeiroticus and P. shqipericus in SW Balkans, and in the central and southern Balkans for P. kurtmuelleri. The LGM models for all taxa exhibit a potentially larger area of distribution than in the current model and observed occurrence (Fig. 7). The currently observed non-overlapping ranges were suitable for studying taxa also during the LGM according to the model predictions.

Figure 7. 

Species distribution modelling (SDM) based on four climatic models for P. epeiroticus, P. shqipericus and kurtmuelleri lineage projected only for the territory of the Balkans. The SDM corresponds to the current time (upper) and the last glacial maximum (LGM) with three Global Climate Models (lower).


Water frogs of the genus Pelophylax are among the most studied European amphibians (e.g., Plötner 2005). Although many studies focused on their unique reproduction mode associated with hybridogenesis, available data on their genetic diversity, or evolutionary history on a broad-geographic scale, are regionally limited, scattered in the literature, or missing. This is also a case of endemic Pelophylax taxa from the Balkan Peninsula, as they were studied only marginally from the range-wide phylogeographic point of view (Lymberakis et al 2007; Hofman et al. 2015; Vucić et al. 2018).

Historical biogeography

Our single gene molecular dating showed that the divergence of studied water frog taxa occurred between 9.8 and 2.5 Mya and was probably related to geological and climatic changes in the Balkans. During the Miocene, the opening of the mid-Aegean trench (~9–12 Mya) and the Messinian Salinity Crisis (MSC; ~5.3 Mya) were both well-discussed events as major factors of the evolution of local biota, including herpetofauna (Akın et al. 2010; Lymberakis and Poulakakis 2010). These events lead to the isolation of Crete from the mainland, i.e., Peloponnese (Paulakakis et al. 2015) and the local biota which occurred there. Although some authors hypothesized the isolation of Crete already during the upper Miocene (~9 Mya) and used it as a calibration point for molecular dating (Akın et al. 2010), the isolation of the island associated with the end of the MSC is well reviewed (Dermitzakis 1990; Poulakakis et al. 2015) and used in similar studies associated to the evolution of Mediterranean vertebrates (Dubey and Dufresnes 2017; Dufresnes et al. 2018a; Spilani et al. 2019; Kiourtsoglou et al. 2021). We thus decided to select the end of MSC as a calibration point isolating P. cretensis from other congeners. Applying this calibration point, we showed that intraspecific divergence in P. shqipericus, P. epeiroticus and P. ridibundus/kurtmuelleri occurred in Pleistocene (0.9 to 0.6 Mya) and might have been connected to the Early-Middle Pleistocene transition (EMPT; 1.4–0.4 Mya). The EMPT was a period of significant climatic turnover in the western Palearctic (Head and Gibbard 2005, 2015; Maslin and Brierley 2015) leading to the increased amplitude of climatic oscillations during the Pleistocene. These changes resulted in an overall widespread temperature decrease, ice volume increase, sea-level fluctuations, aridification, and habitat changes (Kahlke et al. 2011), especially during the so-called “0.9 Mya event” (Strani et al. 2019). The climatic shift reached its peak at approximately 0.9 Mya with an extremely dry climate during both glacial and interglacial periods. At the late stages of the EMPT (0.9–0.4 Mya) intensive cold periods with lowland glaciation occurred, which were followed by warmer periods, with increasing humidity leading to a variety of habitats and the development of forested landscapes (Suc et al. 1995; Bertini 2000; Kahlke et al. 2011; Head and Gibbard 2015). The impact of EMPT on the faunal composition is well-known for mammals (e.g., Kahlke et al. 2011; Strani et al. 2019), but poorly studied in other vertebrates, such as amphibians (Blain et al. 2018).

Another paleoclimatic event that has significantly affected the genetic diversity and range dynamics of the biota has been the Last Glacial Maximum (LGM; ~21 Kya; Hewitt 1999, 2011). The long-term in situ survival of populations in refugial areas has led to high genetic diversity and population differentiation. Subsequent rapid postglacial recolonization of the northern parts of Europe from refugia resulted in lower genetic diversity of the colonizing populations. This pattern of “the southern richness” and “the northern purity” (Hewitt 2004) is apparent also in studied water frogs. Based on the dispersal ability of a species and the opportunities for post-glacial expansion to the north, a lineage or species can be characterized either as a post-glacial re-colonizers (Hewitt 1999, 2000) or as a refugial endemics with limited expansion potential (Bilton et al. 1998). While the Balkan species/populations of Pelophylax show higher genetic diversity, ridibundus lineage, which colonized most of Europe, seems to be genetically uniform according to our data (Figs 14). Similarly, the haplotype diversity of kurtmuelleri lineage in the Balkans seems enormous compared to its haplotype diversity found in northerly located European populations (Figs 4, S2). Pelophylax epeiroticus and P. shqipericus thus represent refugial endemics, while P. ridibundus/kurtmuelleri might be characterized as post-glacial re-colonizers.

Based on the distribution of the haplotype diversity in the Balkans, it could be hypothesized that the microrefugial model might be applied to the kurtmuelleri lineage. Some haplogroups of recent origin have limited distribution and therefore we can assume they could have survived the LGM in Dalmatia (“orange” haplogroup), Peloponnese (“pink” haplogroup), coastal Albania (“green” haplogroup) or Albania and western Greece (“red” and “blue” haplogroups). Unfortunately, the geographic pattern of several haplogroups is mixed, a fact that complicates the detection of microrefugia. The Pleistocene history of the most widely distributed “yellow” haplogroup remains also unclear (Figs 4, S2). If the “yellow” haplogroup had limited distribution during the LGM, it had to expand rapidly (as suggested by demographic analyses) not only to other parts of the Balkans but also to vast areas of western, eastern, central, and northern Europe.

It seems that during the LGM, the southwestern Balkans hosted favourable conditions for water frogs and their distribution did not significantly change (Fig. 7). Considering details from haplotype networks of P. epeiroticus, the south-western continental Greece and northern Peloponnese were probably the suitable areas for the Southern lineage and north-western Greece for the Northern lineage (Fig. 2). In P. shqipericus, we can estimate the location of these areas for the Northern lineage based solely on the haplotype network, lying in the Skadar Lake lowland (Fig. 3). The Skadar Lake system represents the area of local endemism and harbours unique lineages in several animal groups (Grabowski et al. 2014; Pešić and Gloër 2018; Jabłońska et al. 2020) and thus, we might expect that this area had an impact also on the evolutionary history of the Northern lineage. However, the suitable area of the Southern lineage remains questionable. Due to the limited range of the Southern lineage, the whole Adriatic coast of Albania might represent refugium itself for this population. Although this uncertainty might be also a result of sampling bias or extinction of populations with peripheral haplotypes. On the other hand, populations of the kurtmuelleri lineage diverged probably in several microrefugia scattered across the peninsula, especially in the southwestern part (Fig. S2), as mentioned above. Observed resolution may suggest historical population admixture that is unprecedented in the observed phylogeographic patterns of studied Balkan amphibians (e.g., Sotiropoulos et al. 2007; Fijarczyk et al. 2011; Dufresnes et al. 2013; Jablonski et al. 2021). It calls, however, for further research where the influence of the rapid natural dispersion/invasion with a possible human-mediated introduction (see Dufresnes et al. 2018b) should be tested (see Conservation genetic implications). In the case of the ridibundus lineage, insufficient sampling covering only limited parts of the whole range and shallow mitochondrial structure (single haplogroup) through the species range suggest possible rapid colonization, but localization of glacial refugia makes rather impossible.

Except for the populations of the P. ridibundus/kurtmuelleri clade, the endemic species showed range stability or even retraction during the LGM (see also Fig. 7). In contrast, according to BSP of P. epeiroticus, both lineages showed population growth starting at approximately between 20–25 Kya. Based on neutrality tests, the Southern lineage assumes population expansion, whereas no sign of population expansion in the Northern lineage was found (see also Fig. 2). The Northern lineage is endemic to the area between southern Albania and the Ambracian Gulf, including Corfu Island which was several times connected and separated from the mainland during the last 1 Mya (Perissoratis and Conispoliatis 2003). Thus, we expect that this area could represent a long-term microrefugium of the Northern lineage. If we admit the possibility of expansion of the Southern lineage (according to BSP and neutrality tests), it would be probably along coastal areas of the Peloponnese and southern continental Greece up to Ioannina Lake. In the Ioannina Lake and Ambracian Gulf, both lineages were detected in sympatry. For P. shqipericus, population growth started approximately between 35 and 25 Kya in the Southern and Northern lineages, respectively. The statistical sign of the expansion is evident in the Southern lineage forming a contact zone with the Northern lineage in northern Albania and possibly Skadar Lake (see Fig. 3). The Northern lineage is distributed mostly in the vicinity of the Skadar Lake, including the type locality of the species in Virpazar. The Southern lineage covers most of the Adriatic coast of Albania, i.e., from Orikum to the border region between Albania and Montenegro. One individual of the Southern lineage was surprisingly recorded also further north, near Virpazar, Skadar Lake. If this record is a result of natural dispersion, denser sampling in the future should detect the mixed pattern of these two lineages around Skadar Lake. Our results suggest mostly pre-LGM population growth in both P. epeiroticus and P. shqipericus, with the evidence of further rapid post-glacial growth. As the expected times of population growth are close to the LGM period (21 Kya), results might reflect favourable conditions during the LGM for the frog population growth, although conditions differ significantly. On the contrary, southern parts of the Balkan were probably unaffected by glacials and might host favourable conditions for the population growth of water frogs even when environmental conditions in other parts of Europe were not suitable.

In the kurtmuelleri lineage, the first population growth started at approximately 30 Kya (possibly inside the Balkans only), following by the rapid expansion after the LGM (rest of the currently uncovered European range; Fig. 4 and S2). The ridibundus lineage also shows population growth and expansion corresponding approximately with the post-LGM period. The scenario of the expansion of the kurtmuelleri lineage is supported by the presence of kurtmuelleri haplotypes in the northern part of the peninsula and rarely also in central, northern, and eastern Europe. However, its frequency is decreasing northwards as is evident from own and published data (Herczeg et al. 2017; Kolenda et al. 2017; Litvinchuk et al. 2020). Kolenda et al. (2017) predict the occurrence of kurtmuelleri haplotypes and alleles in the geographic region between the northern Balkans and Poland. This is now confirmed by the presence of kurtmuelleri haplotypes in Ukraine and Germany (Figs 4 and S2). According to Litvinchuk et al. (2020), after the LGM, haplotypes associated with kurtmuelleri lineage possibly expanded north along lowlands in Romania, Serbia, or Hungary. During this expansion, haplotypes and alleles of the kurtmuelleri lineage dispersed to large distances across central, eastern, and northern Europe and met with ridibundus-specific haplotypes and individuals with different mtDNA and nDNA possibly hybridized. However, we cannot exclude this pattern could be also influenced by local human-mediated introductions, similar to western Europe or even northern Africa (e.g., Plötner 2005; Dubey et al. 2014; Dufresnes et al. 2017).

Contrasting genetic diversity

Comparisons among three water frog clades (P. epeiroticus, P. shqipericus and P. ridibundus/kurtmuelleri) of the southwestern Balkans revealed significant differences in their intraspecific genetic diversity. For example, both P. epeiroticus and P. shqipericus formed two, diverged and well-supported lineages (Figs 2, 3). Moreover, the presence of a possible third, separated lineage in P. epeiroticus (needing further attention) was detected particularly using PCAs in northern and central Greece (Fig. 1). However, due to the limited sample size of these individuals, we include them in the Southern lineage. The high nucleotide diversity was observed in both endemics, i.e., species with much smaller ranges compared to kurtmuelleri lineage (Figs 2, 3). However, it supports the hypothesis that the southwestern parts of the Balkans are sources of genetic diversity of endemic species/evolutionary lineages showing limited geographic dispersions (e.g., Marzahn et al. 2016; Psonis et al. 2018). For example, endemics to the southwestern Balkans, Triturus macedonicus or Anguis graeca, show a higher level of genetic diversity compared to other taxa of these genera from temperate Europe (Jablonski et al. 2016; Wielstra and Arntzen 2020). This contrast could be connected to the position of areas of long-time survival where lineages or haplotypes evolved, as well as to the natural history of taxa and their capacity for ecological and geographical dispersions or adaptability. A similar pattern could be found also in widespread species in the Balkans (see Sotiropoulos et al. 2007; Pabijan et al. 2014; Jablonski et al. 2016, Wielstra and Arntzen 2020). In P. kurtmuelleri, populations in northern and central Europe fall clearly into a single haplogroup (“yellow”), while in the Balkan Peninsula it forms a much more diverged pattern (Fig. S2). Most of the members of the genus Pelophylax are well known to be associated with water bodies which is probably a limiting factor of their dispersion to the landscape (Mikulíček and Pišút 2012). For example, it is very rare to observe P. epeiroticus out of the water (Sofianidou 1996; Sofianidou and Schneider 1989). Moreover, at least P. epeiroticus has a limited vertical distribution that is up to 500 m a. s. l. (Sofianidou and Schneider 1989). Similarly, the range of P. shqipericus is also limited to lowland areas of Montenegro and Albania. In this sense ca 70 km wide distribution gap between both studied endemic clades in southern Albania is notable. In that gap, the Ceraunian Mountains (up to 2000 m a. s. l.), rising directly from the seacoast, are situated and support the hypothesis of environmental limitations for the population expansion of P. epeiroticus and P. shqipericus due to natural barriers. Thus, we expect that these specific traits may be behind the distribution and the intraspecific genetic difference observed in our data, i.e., limited dispersion and long-term in situ persistence of populations.

On the other hand, the kurtmuelleri lineage is known to be able to live for a certain time out of the water (Sofianidou and Schneider 1989) and is also observed in high elevations (e.g., more than 2,000 m; Szabolcs et al. 2017). It suggests a significantly better dispersion capacity of this lineage (cf. SDM of the kurtmuelleri lineage vs. both endemic taxa, Fig. 7). According to our data, the kurtmuelleri lineage has rich haplotype diversity (126 haplotypes vs. 31 in P. epeiroticus and 25 in P. shqipericus) that geographically present a complex pattern (Figs 4, S2). Overall genetic (nucleotide) variability of the lineage is, however, shallow (see also Vucić et al. 2018) with low nucleotide diversity reflecting isolation in short-lasting time and space windows.

MtDNA, taxonomy and shared polymorphism in P. ridibundus/kurtmuelleri

Although evolutionary reconstructions based merely on mtDNA could face various disadvantages (e.g., hybridization and introgression of mtDNA), in the case of western Palearctic amphibians, mitochondrial data provide a comparable view to biparental markers in sense of the historical biogeography and even taxonomy (e.g., Dufresnes et al. 2019b). Moreover, the wide use of mtDNA in molecular studies on the genus Pelophylax (Cyt b: Dubey et al. 2014; Kolenda et al. 2017; Dufresnes et al. 2017, 2018b; Cyt b + 16S rRNA: Lymberakis et al. 2007, Radojičić et al. 2015; ND2: Litvinchuk et al. 2020; Svinin et al. 2020; ND3: Dubey et al. 2014; Vucić et al. 2018; ND2+ND3: Plötner et al. 2008, 2010, 2012; Akın et al. 2010) creates space for wider-data comparisons with dense sampling in these frogs from different aspects.

Since the description of P. kurtmuelleri, the taxon was considered the Balkan endemics based on the morphometry (Schneider et al. 1993), bioacoustics (Schneider and Sinsch 1992; Schneider et al. 1993; Lukanov et al. 2015) and supposed inability to induce hybridogenesis in interspecific hybrids (Hotz et al. 1985; Berger et al. 1994; Sinsch and Eblenkamp 1994). On the other hand, P. kurtmuelleri shows only weak morphometrical (Papežík et al. 2021), and genetic differentiation from P. ridibundus in mtDNA (Lymberakis et al. 2007; Hofman et al. 2015; Vucić et al. 2018) as well as nDNA (SAI-1 gene; Kolenda et al. 2017; Papežík et al. 2021). The comparison of full mitogenome sequences showed the same conclusions (about 1% of divergence) where the species level distance among other taxa is more than 5% (Hofman et al. 2015). Thus, it is not surprising that the taxon is considered conspecific for a long-time with P. ridibundus (Speybroeck et al. 2010) and Balkans populations currently recognized under the name P. kurtmuelleri presumably representing a closely related lineage. This is also supported by our data as the genetic distance between these two lineages is less (1.4%) than the same distance between lineages detected inside both studied endemic species, P. epeiroticus (2.1%) and P. shqipericus (1.6%). Low levels of genetic divergence are reflected also in PCAs and molecular dating (see above and Figs 1 and 5). Our improved view thus supports the hypothesis that the independent evolution of the ridibundus and kurtmuelleri mitochondrial lineages seems to be young, starting probably in the Pleistocene (Lymberakis et al. 2007). We could conclude that the mitochondrial variability (our data and Hofman et al. 2015) does not support the species status of P. kurtmuelleri and its taxonomic distinctiveness from P. ridibundus. Nevertheless, the description of hybrid zones and the rate of genetic admixture between these evolutionary lineages should be in the scope of further research (see Dufresnes and Mazepa 2020).

Apart from the natural range, haplotypes and alleles of kurtmuelleri lineage were recorded in several European countries in the last decades (e.g., Herczeg et al. 2017; Kolenda et al 2017; Bruni et al. 2020; Litvinchuk et al. 2020). Kolenda et al. (2017) proposed three hypotheses for the origin of P. kurtmuelleri haplotypes in central and northern Europe: (1) recent introduction followed by hybridization of autochthonous and allochthonous populations, (2) incomplete lineage sorting or, (3) hybridization before or during the post-glacial expansion. Some of the above-mentioned reports represent well-documented introductions of individuals from non-native kurtmuelleri lineage, especially to western Europe (Dubey et al. 2014; Dufresnes et al. 2017, 2018b; Bruni et al. 2020). In contrast, such introductions to multiple countries in central, northern, and eastern Europe are not reported. Also, ancestral polymorphism, characteristic by admixture of mtDNA between two taxa regardless of distance from secondary contact of both taxa should carry both types of mtDNA (Toews and Brelsford 2012). However, to the best of our knowledge, no individuals of kurtmuelleri lineage with ridibundus specific mtDNA are known. On the contrary, in the case of hybridization and introgression, one of the involved taxa should keep its mtDNA, while the rate of this mtDNA in another taxon decreases with distance from secondary contact (Toews and Brelsford 2012). This seems to be a case of kurtmuelleri lineage haplotypes and alleles as they seem to be quite rare in central, northern, and eastern Europe. Our data thus rather support the natural colonization of temperate Europe along lowlands in eastern Balkans and eastern Europe as proposed by Kolenda et al. (2017) and Litvinchuk et al. (2020).

Conservation genetic implications

Populations of amphibians are in decline on a global scale with various factors influencing their extinction (reviewed by Hayes et al. 2010). Even though genetic factors affect population extinction on a long-term scale (Frankham 2005), the loss of genetic diversity can increase extinction risk by insufficient adaptability to environmental changes and by fitness reduction (Spielman et al. 2004; Allentoft and O’Brien 2010). Here, we emphasize that the southwestern Balkan region is an important centre of genetic diversity of amphibians in the Palearctic with implications for conservation.

In the southwestern Balkans, P. epeiroticus and P. shqipericus represent endemic species with a limited range of several thousand square kilometers (IUCN SSC Amphibian Specialist Group 2020a,b; Jablonski et al. 2021). Both species are known to be threatened by habitat loss and degradation, and increasing levels of pollution of water bodies due to human overpopulation (IUCN SSC Amphibian Specialist Group 2020a,b; Pafilis and Marangou 2020; Saçdanaku et al. 2020; Crnobrnja-Isailović et al. 2022). This applies mostly to P. shqipericus which is present in the Adriatic coastland and lowlands of Albania that were transformed from swamps and marshes into the desiccated extensive agricultural landscape during the last decades (Haxhiu 1994; Frank et al. 2018; Saçdanaku et al. 2020). Moreover, both species are also collected for consumption by restaurants, the food industry, and by local people, even during the breeding season (Fig. S4.; Frank et al. 2018, Pafilis and Marangou 2020). Due to over-collecting for frog legs trade, P. shqipericus was even considered for listing in CITES (Gratwicke et al. 2010). Pelophylax epeiroticus is also still consumed in local restaurants. However, the assessment of both species in the IUCN Red List category was recently lowered to Near Threatened in P. epeiroticus and Vulnerable in P. shqipericus, respectively (Dufresnes 2019; IUCN SSC Amphibian Specialist Group 2020a,b). With respect to detected endemic genetic diversity discovered in both species and their possible increasing vulnerability, we appeal to a re-evaluation of their IUCN Red List assessment. Similarly, both species are listed in Appendix III of the Berne Convention, and they received only a lower level of protection but lack specific law protection in Albania and Greece. Thus, our study may also serve as a breakpoint for species conservation and management. This is especially true for Albania, where both endemic species occur. Limited ranges of both species together with the genetic intraspecific variability may lead to the increased threat of irreversible loss of their genetic diversity. Or such loss would reduce intraspecific diversity and could lead to an increased risk of possible local extinctions. Habitat alteration and anthropogenic pressure are strong factors in areas where these taxa occur. Thus, the formation of specific areas where these taxa and their genetic variability will be protected is necessary. Therefore, this new view should be considered by the conservation legislation, especially in Albania and Greece.


We would like to thank many of friends and colleagues that provided information, photographs, donated tissue samples, or assisted us during the fieldwork or in the laboratory, especially to Jana Poláková (Comenius University in Bratislava, Slovakia), Nuria Viñuela Rodríguez and Eva Ašenbrenerová (National Museum Prague, Czech Republic), Ivana Buj and Zoran Marčić (University of Zagreb, Croatia), Ivan Bogut (Josip Juraj Strossmayer University of Osijek, Croatia), Dario Marić (Dobrič b. b., Bosnia and Herzegovina), Spase Shumka (Agricultural University of Tirana, Albania), Denik Ulqini (University of Shkodra “Luigj Gurakuqi”, Albania), and Stamatis Zogaris (Hellenic Centre for Marine Research, Greece). We are also grateful to an anonymous reviewer for valuable comments that improved the submitted version of the manuscript. This work was supported by the Academy of Sciences of the Czech Republic (Grant Number RVO 67985904), the Czech Science Foundation (Grant Number GA19-24559S), the Slovak Research and Development Agency under the contract APVV-19-0076 and by the grant of the Scientific Grant Agency of the Slovak Republic VEGA 1/0286/19 and 1/0298/19. The research permits were provided by the Directorate of Forest Management, Ministry for Environment and Energy of the Hellenic Republic (154073/823/7-4-2017, 173857/1638/18-9-2018, 181012/807/28-3-2019, 183226/1246/11-6-2019, 123199/3356/22-12-2020, and 123199/3356/01-02-2021), and the National Agency of Protected Areas, Ministry of Tourism and Environment of the Albanian Republic (No. 480/2019).


  • Akın Ç, Bilgin CC, Beerli P, Westaway R, Ohst T, Litvinchuk SN, Uzzell T, Bilgin M, Hotz H, Guex GD, Plötner J (2010) Phylogeographic patterns of genetic diversity in eastern Mediterranean water frogs have been determined by geological processes and climate change in the Late Cenozoic. Journal of Biogeography 37: 2111–2124.
  • Allegrucci G, Rampini M, Chimenti C, Alexiou S, Di Russo C (2021) Dolichopoda cave crickets from Peloponnese (Orthoptera, Rhaphidophoridae): Molecular and morphological investigations reveal four new species for Greece. The European Zoological Journal 88: 505–524.
  • Allentoft ME, O’Brien J (2010) Global amphibian declines, loss of genetic diversity and fitness: A review. Diversity 2: 47–71.
  • Berger L, Uzzell T, Hotz H (1994) Postzygotic reproductive isolation between Mendelian species of European water frogs. Zoologica Poloniae 39: 209–242.
  • Bilton DT, Mirol PM, Mascheretti S, Fredga K, Zima J, Searle JB (1998) Mediterranean Europe as an area of endemism for small mammals rather than a source for northwards postglacial colonisation. Proceedings of the Royal Society B 265: 1219–1226.
  • Blain HA, Cruz Silva JA, Jiménez Arenas JM, Margari V, Roucoux K (2018) Towards a Middle Pleistocene terrestrial climate reconstruction based on herpetofaunal assemblage from the Iberian Peninsula: State of the art and perspectives. Quaternary Science Reviews 191: 167–188.
  • Bouckaert R, Heled J, Kühnert D, Vaughan T, Wu CH, Xie D, Suchard MA, Rambaut A, Drummond AJ (2014) BEAST 2: A software platform for Bayesian evolutionary analysis. PLoS Computational Biology 10: e1003537.
  • Bruni G, Mirabella I, Domeneghetti D, Fasola M, Bellati A (2020) Will there be a second extinction? Molecular identification of multiple alien water frogs (Pelophylax ridibundus sensu lato) in Tuscany, Central Italy, reveals genetic pollution within a unique hybridogenetic system. Herpetological Journal 30: 147–158.
  • Bucci S, Ragghianti M, Mancino G, Petroni G, Guerrini F, Giampaoli S (1999) Rana/Pol III: A family of SINE-like sequences in the genome of western Palearctic water frogs. Genome 42: 504–511.
  • Casola C, Marraci S, Bucci S, Ragghianti M, Mancino G, Hotz H, Uzzell T, Guex GD (2004) A hAT-related family of interspersed repetitive elements in genome of western Palaearctic water frogs. Journal of Zoological Systematics and Evolutionary Research 42: 234–244.
  • Crnobrnja-Isailović J, Adrović A, Bego F, Čađenović N, Hadžiahmetović Jurida E, Jablonski D, Sterijovski B, Jovanović Glavaš O (2022) The importance of small water bodies’ conservation for maintaining local amphibian diversity in the western Balkans. In: Pešić V, Milošević D, Miliša M (Eds) Small Water Bodies of the Western Balkans. Springer, Cham, 351–387.
  • Dermitzakis DM (1990) Paleogeography, geodynamic processes and event stratigraphy during the late Cenozoic of the Aegean area. Atti Convegni Lincei 85: 263–288
  • Drummond AJ, Rambaut A, Shapiro B, Pybus OG (2005) Bayesian coalescent inference of past population dynamics from molecular sequences. Molecular Biology and Evolution 22: 1185–1192.
  • Dubey S, Leuenberger J, Perrin N (2014) Multiple origins of invasive and ‘native’ water frogs (Pelophylax spp.) in Switzerland. Biological Journal of the Linnean Society 112: 442–449.
  • Dubois A, Ohler A (1994) Frogs of the subgenus Pelophylax (Amphibia, Anura, genus Rana): A catalogue of available and valid scientific names, with comments on name-bearing types, complete synonymies, proposed common names and maps showing all type localities. Zoologica Poloniae 39: 139–204.
  • Dufresnes C (2019) Amphibians of Europe, North Africa & the Middle East. A Photographic Guide. Bloomsbury Wild Life, London, 224 pp.
  • Dufresnes C, Denoël M, di Santo L, Dubey S (2017) Multiple uprising invasions of Pelophylax water frogs, potentially inducing a new hybridogenetic complex. Scientific Reports 7: 6506.
  • Dufresnes C, Leuenberger J, Amrhein V, Bühler C, Thiébaud J, Bohnenstengel T, Dubey S (2018b) Invasion genetics of marsh frogs (Pelophylax ridibundus sensu lato) in Switzerland. Biological Journal of the Linnean Society 123: 402–410.
  • Dufresnes C, Lymberakis P, Kornilios P, Savary R, Perrin N, Stöck M (2018a) Phylogeography of Aegean green toads (Bufo viridis subgroup): Continental hybrid swarm vs. insular diversification with discovery of a new island endemic. BMC Evolutionary Biology 18: 67.
  • Dufresnes C, Mazepa G, Jablonski D, Oliveira RC, Wenseleers T, Shabanov DA, Auer M, Ernst R, Koch C, Ramírez-Chaves HE, Mulder KP, Simonov E, Tiutenko A, Kryvokhyzha D, Wennekes PL, Zinenko OI, Korshunov OV, Al-Johany AM, Peregontsev EA, Masroor R, Betto-Colliard C, Denoël M, Borkin LJ, Skorinov DV, Pasynkova RA, Mazanaeva LF, Rosanov JM, Dubey S, Litvinchuk SN (2019a) Fifteen shades of green: The evolution of Bufotes toads revised. Molecular Phylogenetics and Evolution 141: 106615.
  • Dufresnes C, Strachinis I, Suriadna N, Mykytynets G, Cogălniceanu D, Székely P, Vukov T, Arntzen JW, Wielstra B, Lymberakis P, Geffen E, Gafny S, Kumlutaş Y, Ilgaz Ç, Candan K, Mizsei E, Szabolcs M, Kolenda K, Smirnov N, Géniez P, Lukanov S, Crochet PA, Dubey S, Perrin N, Litvinchuk SN, Denoël M (2019b) Phylogeography of a cryptic speciation continuum in Eurasian spadefoot toad (Pelobates). Molecular Ecology 28: 3257–3270.
  • Dufresnes C, Wassef J, Ghali K, Brelsford A, Stöck M, Lymberakis P, Crnobrnja-Isailović J, Perrin N (2013) Conservation phylogeography: Does historical diversity contribute to regional vulnerability in European tree frogs (Hyla arborea)? Molecular Ecology 22: 5669–5684.
  • Džukić G, Kalezić ML (2004) The biodiversity of amphibians and reptiles in the Balkan Peninsula. In: Griffiths HI, Kryštufek B, Reed JM (Eds) Balkan Biodiversity. Kluwer Academic Publishers, Dodrecht, 167–192.
  • Fick SE, Hijmans RJ (2017) WorldClim 2: New 1km spatial resolution climate surfaces for global land areas. International Journal of Climatology 37: 4302–4315.
  • Fijarczyk A, Nadachowska K, Hofman S, Litvinchuk SN, Babik W, Stuglik M, Gollmann G, Choleva L, Cogălniceanu D, Vukov T, Džukić G, Szymura JM (2011) Nuclear and mitochondrial phylogeography of the European fire-bellied toads Bombina bombina and Bombina variegata supports their independent histories. Molecular Ecology 20: 3381–3398.
  • Frank T, Saçdanaku E, Duda M, Bego F (2018) Amphibian and reptile fauna of the Vjosa river, Albania. Acta ZooBot Austria 155: 323–336.
  • Grabowski M, Bącela-Spychalska K, Pešić V (2014) Reproductive traits and conservation needs of the endemic gammarid Laurogammarus scutarensis (Schäferna, 1922) from the Skadar Lake system, Balkan Peninsula. Limnologica 47: 44–51.
  • Gratwicke B, Evans MJ, Jenkins PT, Kusrini MD, Moore RD, Sevin J, Widt DE (2010) Is the international frog legs trade a potential vector for deadly amphibian pathogens? Frontiers in Ecology and the Environment 8: 438–442.
  • Griffiths HI, Frogley MR (2004) Fossil ostracods, faunistics and the evolution of regional biodiversity. In: Griffiths HI, Kryštufek B, Reed JM (Eds) Balkan Biodiversity. Kluwer Academic Publishers, Dordrecht, 261–272.
  • Guicking D, Joger U, Wink M (2009) Cryptic diversity in a Eurasian water snake (Natrix tessellata, Serpentes: Colubridae): Evidence from mitochondrial sequence data and nuclear ISSR-PCR fingerprinting. Organisms Diversity & Evolution 9: 201–214.
  • Günther R (2004) Rana shqiperica. In: Gasc JP, Cabela A, Crnobrnja-Isailović J, Dolmen D, Grossenbacher K, Haffner P, Lescure J, Martens H, Martinez Rica JP, Maurin H, Oliveira ME, Sofianidou TS, Veith M, Zuiderwijk A (Eds) Atlas of Amphibians and Reptiles in Europe. Muséum national d’Histoire naturelle, Paris, 156–157.
  • Hayes TB, Falso P, Gallipeau S, Stice M (2010) The cause of of global amphibian declines: A developmental endocrinologist’s perspective. Journal of Experimental Biology 213: 921–933.
  • Haxhiu I (1994) The herpetofauna of Albania. Amphibia: Species Composition, Distribution, Habitats. Zoologische Jahrbücher, Abteilung für Systematik, Geographie und Biologie der Tiere 121: 321–334.
  • Head MJ, Gibbard PL (2005) Early-Middle Pleistocene transitions: An overview and recommendation for the defining boundary. In: Head MJ, Gibbard PL (Eds) Early-Middle Pleistocene Transitions: The Land-Ocean Evidence. Geological Society, London, 1–18.
  • Herczeg D, Vörös J, Christiansen DG, Benovics M, Mikulíček P (2017) Taxonomic composition and ploidy level among European water frogs (Anura: Ranidae: Pelophylax) in eastern Hungary. Journal of Zoological Systematics and Evolutionary Research 55: 129–137.
  • Hijmans RJ, Cameron SE, Parra JL, Jones PG, Jarvis A (2005) Very high resolution interpolated climate surfaces for global land areas. International Journal of Climatology 25: 1965–1978.
  • Hofman S, Pabijan M, Osikowski A, Litvinchuk SN, Szymura JM (2015) Phylogenetic relationships among four new complete migenome sequences of Pelophylax (Amphibia: Anura) from the Balkans and Cyprus. Mitochondrial DNA 27: 1–4.
  • Hotz H, Uzzell T (1982) Biochemically detected sympatry of two water frog species: Two different cases in the Adriatic Balkans (Amphibia, Ranidae). Proceedings of the Academy of Natural Sciences of Philadelphia 134: 50–79.
  • Hotz H, Uzzell T (1983) Interspecific hybrids of Rana ridibunda without germ line exclusion of the parental genome. Experientia 39: 538–540.
  • Hotz H, Mancino G, Bucci-Innocenti S, Ragghianti M, Berger L, Uzzell T (1985) Rana ridibunda varies geographically in inducing clonal gametogenesis in interspecific hybrids. Journal of Experimental Zoology 236: 199–210.
  • Hotz H, Uzzell T, Günther R, Tunner HG, Heppich S (1987) Rana shqiperica, a new European water frog species from the Adriatic Balkans (Amphibia, Salientia, Ranidae). Notulae Naturae the Academy of Natural Sciences of Philadelphia 468: 1–3.
  • Jabłońska A, Wrzesińska W, Zawal A, Pešić V, Grabowski M (2020) Long-term within-basin isolation patterns, different conservation units, and interspecific mitochondrial DNA introgression in an amphipod endemic to the ancient Lake Skadar system, Balkan Peninsula. Freshwater Biology 65: 209–225.
  • Jablonski D (2011) Reptiles and amphibians of Albania with new records and notes on occurrence and distribution. Acta Societatis Zoologicae Bohemicae 75: 223–238.
  • Jablonski D, Jandzik D, Mikulíček P, Džukić G, Ljubisavljević K, Tzankov N, Jelić D, Thanou E, Moravec J, Gvoždík V (2016) Contrasting evolutionary histories of the legless lizards slow worms (Anguis) shaped by the topography of the Balkan Peninsula. BMC Evolutionary Biology 16: 99.
  • Jablonski D, Nagy ZT, Avcı A, Olgun K, Kukushkin OV, Safaei-Mahroo B, Jandzik D (2019) Cryptic diversity in the smooth snake (Coronella austriaca). AmphibiaReptilia 40: 179–192.
  • Jablonski D, Gkontas I, Poursanidis D, Lymberakis P, Poulakakis N (2021) Stability in the Balkans: Phylogeography of the endemic Greek stream frog, Rana graeca. Biological Journal of the Linnean Society 132: 829–846.
  • Kahlke RD, García N, Kostopoulos DS, Lacombat F, Lister AM, Mazza PPA, Spassov N, Titov VV (2011) Western Palaearctic palaeoenviroment conditions during the Early and Middle Pleistocene inferred from large mammal communities, and implications for hominin dispersal in Europe. Quarternary Science Reviews 30: 1368–1395.
  • Kiourtsoglou A, Kaliontzopoulou A, Poursanidis D, Jablonski D, Lymberakis P, Poulakakis N (2021) Evidence of cryptic diversity in Podarcis peloponnesiacus and re-evaluation of its current taxonomy; insight from genetic, morphological, and ecological data. Journal of Zoological Systematics and Evolutionary Research 59: 2350–2370.
  • Kolenda K, Pietras-Lebioda A, Hofman S, Ogielska M, Pabijan M (2017) Preliminary genetic data suggest the occurrence of the Balkan water frog, Pelophylax kurtmuelleri, in southwestern Poland. AmphibiaReptilia 38: 187–196.
  • Lanfear R, Fradsen PB, Wright AM, Senfeld Calcott B (2017) partitionfinder 2: New methods for selecting partitioned models of evolution for molecular and morphological phylogenetic analyses. Molecular Biology and Evolution 34: 772–773.
  • Litvinchuk SN, Ivanov AY, Lukonina SA, Ermakov OA (2020) A record of alien Pelophylax species and widespread mitochondrial DNA transfer in Kaliningradskaya Oblasť (the Baltic coast, Russia). BioInvasions Records 9: 599–617.
  • Lukanov S, Tzankov N, Simeonovska-Nikolova D (2015) A comparative study of the matting call of Pelophylax ridibundus and Pelophylax kurtmuelleri (Anura: Ranidae) from syntopic and allotopic populations. Journal of Natural History 49: 257–272.
  • Lymberakis P, Poulakakis N, Manthalou G, Tsigenopoulos CS, Magoulas A, Mylonas M (2007) Mitochondrial phylogeography of Rana (Pelophylax) populations in the Eastern Mediterranean region. Molecular Phylogenetics and Evolution 44: 115–125.
  • Marzahn E, Mayer W, Joger U, Ilgaz Ç, Jablonski D, Kindler C, Kumlutaş Y, Nistri A, Schneeweiss N, Vamberger M, Žagar A, Fritz U (2016) Phylogeography of the Lacerta viridis complex: Mitochondrial and nuclear markers provide taxonomic insights. Journal of Zoological Systematics and Evolutionary Research 54: 85–105.
  • Mellado PV, Valakos ED, Gil MJ, Guerrero F, Lulch J, Navarro P, Maragou P (1999) Herpetological notes from mainland and insular Greece. British Herpetological Society Bulletin 67: 33–38.
  • Mermygkas D, Zikos A, Constantinidis T (2021) Biological traits, habitat preferences and endemism in the flora of Peloponnisos, Greece. Flora Mediterranea 31: 37–52.
  • Mikulíček P, Pišút P (2012) Genetic structure of the marsh frog (Pelophylax ridibundus) populations in urban landscape. European Journal of Wildlife Research 58: 833–845.
  • Mizsei E, Jablonski D, Roussos S, Dimaki M, Ioannidis Y, Nilson G, Nagy ZT (2017) Nuclear markers support the mitochondrial phylogeny of Vipera ursinii-renardi complex (Squamata: Viperidae) and species status for the Greek meadow viper. Zootaxa 4227: 75–88.
  • Musilová R, Zavadil V, Marková S, Kotlík P (2010) Relics of the Europe’s warm past: Phylogeography of the Aesculapian snake. Molecular Phylogenetics and Evolution 57: 1245–1252.
  • Oruçi S (2008) Data on geographical distribution and habitats of the Rana epeirotica in Albania. Natura Montenegrina 7: 419–423.
  • Pabijan M, Zieliński P, Dudek K, Chloupek M, Sotiropoulos K, Liana M, Babik W (2014) The dissection of a Pleistocene refugium: Phylogeography of the smooth newt, Lissotriton vulgaris, in the Balkans. Journal of Biogeography 42: 671–683.
  • Pabijan M, Zieliński P, Dudek K, Stuglik M, Babik W (2017) Isolation and gene flow in a speciation continuum in newts. Molecular Phylogenetics and Evolution 116: 1–12.
  • Pafilis P, Marangou P (2020) Atlas of Amphibians and Reptiles in Greece. Broken Hill Publisher Ltd. Nicosia, Cyprus (in Greek).
  • Papežík P, Kubala M, Jablonski D, Doležálková-Kaštánnková M, Choleva L, Benovics M, Mikulíček P (2021) Morphological differentiation of endemic water frogs (Ranidae: Pelophylax) from the southwestern Balkans. Salamandra 57: 105–123.
  • Perissoratis C, Conispoliatis N (2003) The impacts of sea level changes during latest Pleistocene and Holocene times on the morphology of the Ionian and Aegean seas (SE Alpine Europe). Marine Geology 196: 145–156.
  • Pešić V, Gloër P (2018) The diversity and conservation status of the Molluscs of the Lake Skadar/Shkodra. In: Pešić V, Kamaran G, Kostianoy AG (Eds) The Skadar/Shkodra Lake Environment. Springer, Cham, 295–310.
  • Plötner J (2005) Die westpaläarktischen Wasserfrösche: von Märtyrern der Wissenschaft zur biologischen Sensation. Laurenti Verlag, Bielefeld, 160 pp.
  • Plötner J, Baier F, Akın Ç, Mazepa G, Schreiber R, Beerli P, Litvinchuk SN, Bilgin CC, Borkin L, Uzzell T (2012) Genetic data reveal that water frogs of Cyprus (genus Pelophylax) are an endemic species of Messinian origin. Zoosystematics and Evolution 88: 261–283.
  • Plötner J, Uzzell T, Beerli P, Akın C, Bilgin CC, Haefeli C, Ohst T, Köhler F, Schreiber R, Guex GD, Litvinchuk SN, Westaway R, Reyer HU, Pruvost N, Hotz H (2010) Genetic divergence and evolution of reproductive isolation in eastern Mediterranean water frogs. In: Glaubrecht M (Ed.) Evolution in Action. Case Studies in Adaptive Radiation and the Origin of Biodiversity. Springer, Berlin, 372–403.
  • Plötner J, Uzzell T, Beerli P, Spolsky C, Ohst T, Litvinchuk SN, Guex GD, Reyer HU, Hotz H (2008) Widespread unidirectional transfer of mitochondrial DNA: A case in western Palaearctic water frogs. Journal of Evolutionary Biology 21: 668–681.
  • Podnar M, Mađarić BB, Mayer W (2014) Non-concordant phylogeographical patterns of three widely codistributed endemic endemic Western Balkans lacertid lizards (Reptilia, Lacertidae) shaped by specific habitat requirements and different responses to Pleistocene climatic oscillations. Journal of Zoological Systematics and Evolutionary Research 52: 119–129.
  • Polović L, Čađenović N (2014) The herpetofauna of the Great Ulcinj Beach area including Ada Island (Montenegro). Turkish Journal of Zoology 38: 104–107.
  • Poulakakis N, Kapli P, Lymberakis P, Trichas A, Vardinoyiannis K, Sfenthourakis S, Mylonas M (2015) A review of phylogeographic analyses of animal taxa from the Aegean and surrounding regions. Journal of Zoological Systematics and Evolutionary Research 53: 18–32.
  • Psonis N, Antoniou A, Karameta E, Leaché AD, Kotsakiozi P, Darriba D, Kozlov A, Stamatakis A, Poursanidis D, Kukushkin O, Jablonski D, Crnobrnja-Isailović J, Gherghel I, Lymberakis P, Poulakakis N (2018) Resolving complex phylogeographic patterns in the Balkan Peninsula using closely related wall-lizard species as a model system. Molecular Phylogenetics and Evolution 125: 100–115.
  • Psonis N, Antoniou A, Kukushkin O, Jablonski D, Petrov B, Crnobrnja-Isailović J, Soutiropoulos K, Gherghel I, Lymberakis P, Poulakakis N (2017) Hidden diversity in the Podarcis tauricus (Sauria, Lacertidae) species subgroup in the light of multilocus phylogeny and species delimitation. Molecular Phylogenetics and Evolution 106: 6–17.
  • Radea C, Louvrou I, Bakolitsas K, Economou-Amilli A (2017) Local endemic and threatened freshwater hydrobiids of Western Greece: Elucidation of anatomy and new records. Folia Malacologica 25: 1.
  • Radojičić JM, Krizmanić I, Kasapidis P, Zouros E (2015) Extensive heteroplasmy in hybrid water frog (Pelophlylax spp.) populations from Southeast Europe. Ecology and Evolution 5: 4529–4541.
  • Ragghianti M, Bucci S, Guerrini F, Mancino G (1999) Characterization of two repetitive DNA families (RrS1 and Rana Pol/III) in the genomes of the Palaearctic green frogs. Italian Journal of Zoology 66: 255–263.
  • Rambaut A, Drummond AJ, Xie D, Baele G, Suchard MA (2018) Posterior summarization in Bayesian phylogenetics using Tracer 1.7. Systematic Biology 67: 901–904.
  • Recuero E, Buckley D, García-París M, Arntzen JW, Cogălniceanu D, Martínez-Solano I (2014) Evolutionary history of Ichtyosaura alpestris (Caudata, Salamandridae) inferred from the combined analysis of nuclear and mitochondrial markers. Molecular Phylogenetics and Evolution 81: 207–220.
  • Reed JM, Kryštufek B, Eastwood WJ (2004) The physical geography of the Balkans and nomenclature of place names. In: Griffiths HI, Kryštufek B, Reed JM (Eds) Balkan Biodiversity. Kluwer Academic Publishers, Dordrecht, 9–22.
  • Ronquist F, Teslenko M, van der Mark P, Ayres DL, Darling A, Höhna S, Larget B, Liu L, Suchard MA, Huelsenback JP (2012) MrBayes 3.2: Efficient Bayesian phylogenetic inference and model choice across a large model space. Systematic Biology 61: 539–542.
  • Rozas J, Ferrer-Mata A, Sánchez-Del Barrio JC, Guirao-Rico S, Librado P, Ramos-Onsins SE, Sánchez-Garcia A (2017) DnaSP 6: DNA sequence polymorphism analysis of large data sets. Molecular Biology and Evolution 34: 3299–3302.
  • Saçdanaku E, Mesiti A, Bytyçi N, Halilaj M (2020) Species Action Plan: “Measures for Protection of Albanian frog (Pelophylax shqipericus) in the Bay of Vlora”. EcoAlbania, Tirana.
  • Sagonas K, Karameta E, Kotsakiozi P, Poulakakis N (2019) Cross-species testing of nuclear markers in Pelophylax water frogs in Greece and examination of their power to detect genetic admixture. AmphibiaReptilia 41: 1–7.
  • Savić IR (2008) Diversification of the Balkan fauna: Its origin, historical development and present status. In: Makarox SE, Dimitrijević RN (Eds) Advances in Arachnology and Developmental Biology. SASA Publishing, Belgrade, 57–79.
  • Schneider H, Haxhiu I (1994) Mating-call analysis and taxonomy of the water frogs in Albania (Anura: Ranidae). Zoologische Jahrbücher, Abteilung für Systematik, Geographie und Biologie der Tiere 121: 248–262.
  • Schneider H, Sinsch U (1992) Mating call variation in lake frogs referred to as Rana ridibunda Pallas, 1771. Taxonomic implications. Zeitschrift für Zoologische Systematik und Evolutionsforschung 30: 297–315.
  • Schneider H, Sofianidou TS, Kyriakopoulou-Sklavounou P (1984) Bioacoustic and morphometric studies in water frogs (genus Rana) of Lake Ioannina in Greece, and description of a new species (Anura, Amphibia). Zeitschrift für Zoologische Systematik und Evolutionsforschung 22: 349–366.
  • Sinsch U, Eblenkamp B (1994) Allozyme variation among Rana balcanica, R. levantina, and R. ridibunda (Amphibia: Anura): Genetic differentiation corroborates bioacoustically detected species status. Zeitschrift für Zoologische Systematik und Evolutionsforschung 32: 35–43.
  • Sinsch U, Schneider H (1996) Bioacoustic assessment of the taxonomic status of pool frog populations (Rana lessonae) with reference to a topotypical population. Journal of Zoological Systematics and Evolutionary Research 34: 63–73.
  • Sofianidou TS, Schneider H (1989) Distribution range of the Epirus frog Rana epeirotica (Amphibia: Anura) and the composition of the water frog populations in Western Greece. Zoologischer Anzeiger 223: 13–25.
  • Sofianidou TS, Schneider H, Asimakopoulos B (1987) Distribution of the water frog Rana epeirotica (Amphibia, Anura) in Greece. In: van Gelden JJ, Strijbosch H, Bergers PJM (Eds) Proceedings of the Fourth Ordinary General Meeting of the Societas Europaea Herpetologica. Faculty of Sciences Nijmegen, Nijmegen, 365–367.
  • Sotiropoulos K, Eleftherakos K, Džukić G, Kalezić ML, Legakis A, Polymeni RM (2007) Phylogeny and biogeography of the alpine newt Mesotriton alpestris (Salamandridae, Caudata), inferred from mtDNA sequences. Molecular Phylogenetics and Evolution 45: 211–226.
  • Speybroeck J, Beukema W, Bok B, Van Der Voort J, Velikov I (2016) Field Guide to the Amphibians and Reptiles of Britain and Europe. Bloomsbury, London, 432 pp.
  • Speybroeck J, Beukema W, Crochet PA (2010) A tentative species list of the European herpetofauna (Amphibia and Reptilia) – an update. Zootaxa 2492: 1–27.
  • Speybroeck J, Beukema W, Dufresnes C, Fritz U, Jablonski D, Lymberakis P, Martínez-Solano I, Razzeti E, Vamberger M, Vences M, Vörös J, Crochet PA (2020) Species list of the European herpetofauna – 2020 update by the Taxonomic Committee of the Societas Europaea Herpetologica. AmphibiaReptilia 41: 139–189.
  • Spielman D, Brook BW, Frankham R (2004) Most species are not driven to extinction before genetic factors impact them. Proceedings of the National Academy of Sciences 101: 15261–15264.
  • Spilani L, Bougiouri K, Antoniou A, Psonis N, Poursanidis D, Lymberakis P, Poulakakis N (2019) Multigene phylogeny, phylogeography and population structure of Podarcis cretensis species group in south Balkans. Molecular Phylogenetics and Evolution 138: 193–204.
  • Strachinis I, Poulakakis N, Karaiskou N, Patronidis P, Patramanis I, Poursanidis D, Jablonski D, Triantafyllidis A (2021) Phylogeography and systematics of Algyroides (Sauria: Lacertidae) of the Balkan Peninsula. Zoologica Scripta 50: 282–299.
  • Strani F, DeMiguel D, Alba DM, Moyà-Solà S, Belluci L, Sardella R, Madurell-Malapeira J (2019) The effects of the “0.9 Ma event” on the Mediterranean ecosystems during the Early-Middle Pleistocene transition as revealed by by dental wear patterns of fossil ungulates. Quarternary Science Reviews 210: 80–89.
  • Suc JP, Bertini A, Combourieu-Nebout N, Diniz F, Leroy S, Russo-Ermolli E, Zheng Z, Bessais E, Ferrier J (1995) Structure of West Mediterranean vegetation and climate since 5.3 Ma. Acta Zoologica Cracoviensis 38: 3–16.
  • Svinin AO, Dedukh DV, Borkin LJ, Ermanov OA, Ivanov AY, Litvinchuk JS, Zamaletnidov RI, Mikhaylova RI, Trubyanov AB, Skorinov DV, Rosanov YM, Litvinchuk SN (2021) Genetic structure, morphological variation, and gametogenic peculiarities on water frogs (Pelophylax) from northeastern European Russia. Journal of Zoological Systematics and Evolutionary Research 59: 1–17.
  • Szabolcs M, Mizsei E, Jablonski D, Vági B, Mester B, Végvári Z, Szabolcs L (2017) Distribution and diversity of amphibians in Albania: New data and foundations of a comprehensive database. AmphibiaReptilia 38: 435–448.
  • Thuiller W, Georges D, Gueguen M, Engler R, Breiner F, Lafourcade B (2022) biomod2: Ensemble Platform for Species Distribution Modeling. R package version 4.0.
  • Valakos ED, Pafilis P, Sotiropoulos K, Lymberakis P, Maragou P, Foufopoulos J (2008) The Amphibians and Reptiles of Greece. Chimaira, Frankfurt am Main, 463 pp.
  • Vucić M, Jelić D, Klobučar GIV, Prkljačić B, Jelić M (2018) Molecular identification of species and hybrids of water frogs (genus Pelophylax) from Lake Skadar, Southeast Adriatic drainages (Amphibia: Ranidae). Salamandra 54: 147–157.
  • Wielstra B, Arntzen JW (2020) Extensive cytonuclear discordance in a crested newt from the Balkan Peninsula glacial refugium. Biological Journal of the Linnean Society 130: 578–585.

Supplementary material

Supplementary material 1

Tables S1, S2 and Figures S1–S4

Papežík P, Mikulíček P, Benovics M, Balogová M, Choleva L, Doležálková-Kaštánková M, Lymberakis P, Mizsei E, Papežíková S, Poulakakis N, Saçdanaku E, Szabolcs M, Šanda R, Uhrin M, Vukić J, Jablonski D (2023)

Data type: .zip

Explanation note: Table S1. List of individuals included in this study with locality, haplotypes, species identification according to mitochondrial or microsatellite data, and GenBank accession numbers. The table is available as a separate Excel file. — Table S2. List of Pelophylax spp. GenBank ND2 sequences used in the time calibrated tree analysis. Sequences used as an outgroup are shaded in gray. — Figure S1. Maximum likelihood tree and phylogenetic relationships of the Pelophylax spp. haplotypes and map of sampled localities with indicated mtDNA of sampled individuals. — Figure S2. Phylogeographic structure and distribution of detected haplogroups in the kurtmuelleri lineage. — Figure S3. The plots under current climate conditions show evaluation metrics for P. epeiroticus, P. kurtmuelleri, and P. shqipericus. — Figure S4. Dozens of specimens of the kurtmuelleri lineage were killed for frog legs consumption along a drainage channel in Shënepremte, central Albania in 2018. Several other piles were found at the same site. Photo by Simona Papežíková.

This dataset is made available under the Open Database License ( The Open Database License (ODbL) is a license agreement intended to allow users to freely share, modify, and use this Dataset while maintaining this same freedom for others, provided that the original source and author(s) are credited.
Download file (1.55 MB)
login to comment