Research Article |
Corresponding author: Daniel Jablonski ( daniel.jablonski@uniba.sk ) Academic editor: Raffael Ernst
© 2023 Petr 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.
This is an open access article distributed under the terms of the Creative Commons Attribution License (CC BY 4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.
Citation:
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) Comparative mitochondrial phylogeography of water frogs (Ranidae: Pelophylax spp.) from the southwestern Balkans. Vertebrate Zoology 73: 525-544. https://doi.org/10.3897/vz.73.e95220
|
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 (
The Epirus water frog, Pelophylax epeiroticus (Schneider, Sofianidou & Kyriakopoulou-Sklavounou, 1984) is distributed from southern Albania to northwestern Peloponnese (west of the Pindos range;
One of the most endangered European amphibian species, the Albanian pool frog, P. shqipericus (
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 (
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.
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
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 (
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 relationships were inferred using maximum-likelihood (ML) and Bayesian interference (BI) by RAxML v. 8.1.12 (
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) (
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 (
For a better presentation of the intraspecific evolution, the parsimony network analysis with the algorithm of TCS (
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 (
Demographic history was analysed in BEAST 2, using Bayesian Skyline plots (BSP,
With the combination of literature and our own field data, we altogether collected 101 records for P. epeiroticus (
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
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 |
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.
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.
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.
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.
According to phylogenetic analyses and molecular dating (Fig.
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
Pelophylax shqipericus is represented by 69 sequences of ND2 with 25 unique haplotypes, hd = 0.91 and π = 0.92% (Table
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
Bayesian skyline plot (BSP) for P. epeiroticus (Fig.
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.
In P. shqipericus, population expansion based on BSP started approximately at ~25 Kya. MSD and neutrality tests, however, showed no significant population growth (Table
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.
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
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.,
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 (
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;
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
It seems that during the LGM, the southwestern Balkans hosted favourable conditions for water frogs and their distribution did not significantly change (Fig.
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.
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.
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
On the other hand, the kurtmuelleri lineage is known to be able to live for a certain time out of the water (
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.,
Since the description of P. kurtmuelleri, the taxon was considered the Balkan endemics based on the morphometry (
Apart from the natural range, haplotypes and alleles of kurtmuelleri lineage were recorded in several European countries in the last decades (e.g.,
Populations of amphibians are in decline on a global scale with various factors influencing their extinction (reviewed by
In the southwestern Balkans, P. epeiroticus and P. shqipericus represent endemic species with a limited range of several thousand square kilometers (
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).
Tables S1, S2 and Figures S1–S4
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á.