https://orcid.org/0000-0002-5394-0114
Abstract
We still have little knowledge concerning the phylogeography of amphibians and reptiles from the Balkan Peninsula compared with the other two Mediterranean peninsulas. This raises concerns for endemic taxa from these peninsulas, because it might interfere with further conservation efforts. Here we focus on the endemic Greek stream frog (Rana graeca) and reconstruct its biogeography and evolutionary history. Using four genetic markers (Cytb, 16S, COI and BDNF) in > 350 sequences covering the whole distribution range, we conducted phylogenetic, demographic and ecological niche analyses, which revealed the phylogeography of this species. Surprisingly, this examination of R. graeca reveals a very shallow level of intraspecific genetic variability through the Balkans, with two main, statistically supported lineages having a partly sympatric distribution. The most variable marker was Cytb, which showed 19 haplotypes in 123 analysed sequences in the whole species distribution area. Here presented genetic data, together with the environmental niche projection and demographic analyses suggest that R. graeca was probably affected only marginally by climatic oscillations, with the Hellenides as the most suitable area for the occurrence of the species in different geological periods. This is consistent with the observed genetic diversity, which is mostly related to these mountains. Although the species shows a certain level of phenotypic variability and ecological preferences, this might be related to species plasticity affected by the micro-climatic conditions in small areas, which merits further research. Comparing phylogeography of other amphibian and reptile species in the Balkans, we showed that the observed pattern represents a new view on the phylogeography of the Balkan herpetofauna.
INTRODUCTION
Biodiversity hotspots are areas with an exceptionally high ecosystem, species and genetic diversity (Hewitt, 2011). Examples of such regions in the Western Palaearctic are the Balkan Peninsula and/or overall Eastern Mediterranean. Their fauna has experienced dramatic shifts during past climatic and topographic changes, which has left a mark on their current species and genetic diversity (Schmitt & Varga, 2012; Salvi et al., 2013; Poulakakis et al., 2015; Jablonski et al., 2016, 2019; Kindler et al., 2018). The Balkans has a complex topography, with the presence of different mountain reliefs and numerous environmental niches that have preserved biodiversity during the Pliocene/Pleistocene climatic changes (Griffiths et al., 2004). Compared with the other two southern European Peninsulas (Iberian and Italian), the Balkan Peninsula is phylogeographically less studied, although it is an important source of species diversity and endemism as a consequence of its role as the late Miocene radiation centre (Kryštufek et al., 2007; Dufresnes et al., 2013, 2019a, b; Poulakakis et al., 2015). These endemics stayed located mainly in or around the Dinarides and Hellenides, which is why they mostly have small distribution ranges (Wielstra et al., 2014, 2018; Jablonski et al., 2016; Kotsakiozi et al., 2018). A good example is the Peloponnese (south continental Greece), which hosts several endemic taxa either at the species level (e.g. Anguis cephallonica, Algyroides moreoticus and Podarcis peloponnesiacus) or at the genus level (Hellenolacerta graeca; Valakos et al., 2008). Despite their small geographical ranges, some of these taxa are characterized by a higher level of intraspecific genetic diversity (Dufresnes et al., 2013; Pabijan et al., 2015; Psonis et al., 2017, 2018; Spilani et al., 2019).
Another endemic species in the Balkans is Rana graeca Boulenger, 1891 (Figs 1, 2), formerly considered also as a part of the Apennine frog fauna (Rana graeca italica Dubois, 1987). Picariello et al. (2002), however, showed that R. italica is not a subspecies of R. graeca but a different species. This was the opposite of the findings of Asimakopoulos (1994), who stated morphological uniformity between R. graeca and R. italica. Later, Veith et al. (2003) clearly confirmed the findings of Picariello et al. (2002) and showed that R. graeca forms a monophyletic clade with Rana dalmatina and Rana latastei and that this clade is sister to the clade of Rana species from the Anatolian and Caucasian region (Anatolian clade; Veith et al., 2003). Both R. graeca and R. italica are currently recognized as endemic taxa. Rana graeca occurs in almost the whole of continental Greece, south-western Bulgaria, North Macedonia, Albania, southern and central Serbia and Montenegro and extends up to north-western Bosnia and Herzegovina, close to the Croatian border (Asimakopoulos & Grossenbacher, 2014; Sillero et al., 2014; Šukalo et al., 2015; Šunje et al., 2017, 2018; Fig. 2). Interestingly, the species is missing in most of the large islands that are close to its mainland range (i.e. Corfu, Kefalonia, Lefkada and Zakynthos) except Thasos (North Aegean). This species often lives in valleys that have streams, springs or small rivers with cold, clear water. This species has a high vertical distribution, starting from sea level and reaching high mountain areas (≤ 2100 m a.s.l.; Stojanov et al., 2011; Asimakopoulos & Grossenbacher, 2014; Szabolcs et al., 2017). This might be an indication of the ecological plasticity of the species. This frog also shows a certain level of phenotypic variability (Bringsøe, 2011; Asimakopoulos & Grossenbacher, 2014; Fig. 1), which suggests that the external morphology could be interpreted as a genetic signal. We decided, therefore, to delve deeper into the evolutionary history of R. graeca based on molecular data, which are less known compared with other taxa of the genus (Veith et al., 2003; Carranza & Arribas, 2008; Vences et al., 2013; Canestrelli et al., 2014; Yuan et al., 2016; Teixeira et al., 2018). The phylogeography of R. graeca was studied recently by Šunje et al. (2018), based on one mitochondrial DNA (mtDNA) marker and including only 22 samples along the distribution range of the species. The authors stated that a wider sampling and more genetic markers were required to resolve the evolutionary history of the species (Šunje et al., 2018).
Individuals of Rana graeca from different geographical origins express different colour patterns. A, individual from the Theth valley, Albania. B, individual from Banjë, Albania. C, individual from Samarina, Greece. D, individual from Chora Getson, Peloponnese, Greece. Photographs by Daniel Jablonski.
Geographical distribution of samples of the Greek stream frog (Rana graeca) used exclusively in our study, with colours corresponding to the main recovered phylogenetic lineages (for locality details of particular samples, see Supporting Information, Table S1). For geographical coverage including published sequences, see Figure 4. The distribution range of the species is indicated in dark brown (Asimakopoulos, 1994; Asimakopoulos & Grossenbacher, 2014; Sillero et al., 2014; Šukalo et al., 2015; Šunje et al., 2017, 2018). The pictured specimen originates from the locality Lojmë, Albania.
Herein, we investigate the multilocus intraspecific phylogeny and phylogeography of R. graeca, combining phylogenetic trees (gene trees and species trees), network and demographic analyses with past and current species distribution modelling. Our objectives were as follows: (1) to evaluate the evolutionary history, genetic diversity and phylogeography of the species; (2) to estimate demography and species range fluctuations; and (3) to carry out a biogeographical comparison of the studied species with amphibians and reptiles that are endemic in the continental Balkans.
MATERIAL AND METHODS
SAMPLE COLLECTION
In total, 101 R. graeca tissue samples were obtained from 55 localities of four countries and used in the phylogenetic, phylogeographical and demographic analyses, covering the whole species range (Fig. 2; Supporting Information, Table S1). The Natural History Museum of Crete (NHMC) specimens (N = 54) have been collected during its research activities, for the most part the recording of the Greek Fauna through collections, from 1995 to the present, following the provisions of the Presidential Decree 67/81. Collected specimens, after a small abdominal dissection to allow ethanol to enter the body cavities, were stored initially in 96% ethanol, which was changed to fresh 96% ethanol upon return of the expedition members to the NHMC premises. For the present study, a small piece of muscle tissue from the thigh of the adult or, alternatively, from the tail of the tadpole was used. Other material was obtained from road-killed individuals (liver and muscle) or, alternatively, from living animals as blood droplets or a ting piece of a digit on the right hindleg.
To place our data into a phylogeographical context, we compiled an additional dataset supplemented by all known published sequences of R. graeca from the GenBank database (Veith et al., 2003; Ilić et al., 2016; Yuan et al., 2016; Šunje et al., 2018). These published sequences represent mostly the cytochrome b marker and come from Albania, Bosnia and Herzegovina, Greece, Montenegro and Serbia (Figs 2, 3; Supporting Information, Table S1).
A, the dated phylogeny of the genus Rana, showing the studied species. Numbers between clades of species show the estimated time of divergence (in millions of years ago). B, the Bayesian tree reconstructed from the concatenated dataset (2707 bp; Supporting Information, Table S1), rooted with the sequence of congeneric sister species (Rana dalmatina and Rana temporaria) and Pelophylax spp. (not shown). Numbers above the branches show posterior probabilities/maximum likelihood bootstrap support values. The scale bar corresponds to one substitution per 100 nucleotide positions. Each terminal branch represents a sample with its locality number (see also Fig. 2; Supporting Information, Table S1). Phylogenetic lineage colours correspond to those used in Figures 2 and 4.
Moreover, 13 R. dalmatina and two Rana temporaria from the southern Balkans (the NHMC collection) were also obtained and used in the phylogenetic and phylogeographical analyses. Additionally, six R. dalmatina, 15 R. temporaria, one Rana pyrenaica, one Rana arvalis, one Rana asiatica and one Rana macrocnemis were were retrieved from GenBank based on the availability of the genetic markers used in our study and used in phylogenetic analyses. Finally, 16 Pelophylax specimens (two Pelophylax cretensis, two Pelophylax kurtmuelleri, three Pelophylax cerigensis and nine Pelophylax cf. bedriagae) were also collected by the NHMC and used in the chronophylogenetic analyses (Supporting Information, Table S1), in which these specimens were used as external calibration points based on the results of previously published studies (Lymberakis et al., 2007; Poulakakis et al., 2013).
WET LAB AND SEQUENCING
Total genomic DNA from all sampled individuals was extracted using a standard ammonium acetate protocol (Bruford et al., 1998) from muscle, liver or blood. Double-stranded polymerase chain reaction (PCR) was used to amplify partial sequences of three mitochondrial genes (mtDNA) encoding cytochrome oxidase subunit I (COI) and cytochrome b (Cytb) and the large ribosomal subunit (16S rRNA), in addition to one nuclear gene (nuDNA) encoding brain-derived neurotrophic factor (BDNF). Primers and conditions used in PCR amplification and in cycle sequencing reactions are given in the Supporting Information (Table S2). Single-stranded sequencing of the PCR product was performed using the Big-Dye Terminator (v.3.1) Cycle Sequencing kit on an ABI3730 automated sequencer, following the manufacturer’s protocol and using the same primers as in the PCR. Sequences were viewed and edited using CodonCode Aligner v.3.7.1 (CodonCode Corporation). The authenticity of the sequences and the homology to the targeted genes were evaluated with a BLAST search in the NCBI genetic database (http://blast.ncbi.nlm.nih.gov/Blast.cgi). All coding gene sequences (Cytb, COI and BDNF) were translated (DnaSP v.6.0; Rozas et al. 2017) before further analysis, in order to check for the existence of stop codons.
All newly determined sequences have been deposited in GenBank (Supporting Information, Table S1). The alignment of the sequences was performed separately for each locus using the algorithm ClustalW as is implemented in MEGA X (Kumar et al., 2018).
Phylogenetic and phylogeographical analyses
In order to investigate the phylogenetic relationships, genetic diversity and demography of R. graeca, four datasets were constructed: (1) a concatenated DNA dataset with all four loci (16S, Cytb, COI and BDNF); (2) a concatenated DNA dataset with mitochondrial loci (16S, Cytb and COI); (3) a concatenated DNA dataset with 16S and Cytb, excluding COI, in order to access the effect of missing data (see Results section); and (4) a dataset with the nuclear marker (BDNF), all including R. graeca, R. dalmatina, R. temporaria, R. pyrenaica, R. arvalis, Rana macrocnemis, R. asiatica and Pelophylax spp. that were used as an outgroup.
The mtDNA alignment was partitioned into seven blocks, including three blocks for the first, second and third codon positions for the protein-coding locus Cytb, three blocks for the first, second and third codon positions for the protein-coding locus COI, and one block for 16S rRNA. Likewise, the nuDNA alignment (BDNF) was also partitioned into three blocks for its first, second and third codon positions (Supporting Information, Table S3). These initial partition schemes were loaded into PartitionFinder v.2.1 (Guindon et al., 2010; Lanfear et al., 2017) in order to calculate and select the best-fitting partitioning scheme and evolutionary models. The results were implemented in each program (MrBayes and RAxML), based on the Bayesian information criterion (BIC), the greedy algorithm and considering that the blocks of each alignment had linked branch lengths. The models that included both gamma distribution and invariable sites were ignored (Yang, 2006).
Phylogenetic reconstruction was conducted using two different approaches, maximum likelihood (ML) and Bayesian inference (BI). The same parameters were used for all datasets. Maximum likelihood analyses were performed using RAxML v.8.1.21 as implemented through raxmlGUI v.1.5 (Silvestro & Michalak, 2011). The best ML tree for each dataset was selected from 50 iterations, and the confidence of the branches was assessed further based on 1000 thorough bootstrap replicates. Bayesian inference analyses were conducted in MrBayes v.3.2.7 (Ronquist et al., 2012), with four runs and eight Markov chain Monte Carlo (MCMC) chains per run for 2 × 107 generations, sampling every 1000th generation. Several MCMC convergence diagnostics were used to check for convergence and stationarity [the plot of the generation vs. the log probability of the data (the log-likelihood values), the average standard deviation of split frequencies, the average potential scale reduction factor, and the minimum value of minimum estimated sample sizes (ESS)]. The lnL was stabilized after ~6 × 106 generations for the mitochondrial and the combined analysis and after 2 × 106 generations for the analysis of the nuclear dataset. The first 35% (7 × 106) trees were discarded as burn-in, as a measure to sample from the stationary distribution and avoid the possibility of including random, suboptimal trees. A majority-rule consensus tree was then produced from the posterior distribution of trees, and the posterior probabilities were calculated as the percentage of samples recovering any particular clade. Posterior probabilities ≥ 0.95 indicate statistically significant support (Huelsenbeck & Ronquist, 2001).
Given that a network approach can sometimes generate a more effective presentation of the intraspecific evolution than the tree-based phylogenetic approaches (Posada & Crandall, 2001), we also constructed haplotype networks for three mtDNA and one nuclear gene using the 95% limit of parsimony (TCS algorithm) as implemented in the software PopArt (http://popart.otago.ac.nz). The resultant phylogenetic pattern was then presented according to the frequency of a particular lineage in the haplotype network. For network analyses, we compared overall 382 sequences (107 sequences from 16S of 387 bp, 123 from Cytb of 322 bp, 53 from COI of 631 bp, and 99 from BDNF of 453 bp; Supporting Information, Table S1) with the same sequence length for the particular marker.
DIVERGENCE TIME ESTIMATIONS
The species tree and divergence times of Rana species in the southern Balkans were estimated using the dataset containing both mtDNA and nuDNA sequences using the Standard Template BEAST 2 package as implemented in BEAST 2 v.2.6.2 (Bouckaert et al., 2019). The input files (xml format) were created using BEAUti v.2.6.2. We used two unlink partitions, one for mtDNA (unlinked site and clock models and linked trees) and one for nuDNA. We used the model averaging tool bModelTest (Bouckaert & Drummond, 2017) to select the most appropriate substitution model. The Yule model for speciation and the uncorrelated lognormal model for describing the molecular clock were selected. Regarding the divergence time estimation, a calibration point retrieved from the study on Pelophylax phylogeny was used (Lymberakis et al., 2007): the split of P. cretensis from all other Pelophylax species at ~5–5.5 Mya, using normal distribution. This calibration point is based on a very well-known paleogeological event, which is the permanent isolation of the island of Crete from the Peloponnese (Meulenkamp, 1985; Dermitzakis, 1990). For this reason, we included in the analysis the corresponding sequences of P. cretensis, P. kurtmuelleri, P. cf. bedriagae and P. cerigensis. The MCMC analysis was run for 6 × 108 generations, saving the result every 5 × 103 generations. The first 25% of the saved trees were discarded after inspection of the log files with TRACER v.1.6 (Rambaut et al., 2014). The maximum clade credibility (MCC) tree that best represented the posterior distribution was identified using TreeAnnotator v.2.6.2 (also included in BEAST 2).
DNA polymorphism and demography
The sequence divergence (p-distances) was estimated using MEGA X. The program DnaSP v.6.0 (Rozas et al., 2017) was used to estimate the number of haplotypes (h), haplotype diversity (hd), number of segregating sites (S), nucleotide diversity (π) and Watterson’s theta per site (θW).
The past population dynamics were inferred using the Bayesian coalescent-based approach of the Bayesian skyline plot (BSP; Drummond et al., 2005), as implemented in BEAST 2. This method computes the effective population size through time directly from sampled sequences and does not require a specific a priori assumed demographic model. This method was applied for all R. graeca sequences of the Cytb dataset owing to the highest number of available sequences of this variable marker. We used the methodology of Vences et al. (2013), a uniform prior for the mean substitution rate with the initial value of 0.005 mutations per site/Myr. Preliminary analyses were run using both a strict molecular clock and the uncorrelated lognormal relaxed molecular clock. Given that the parameter of the standard deviation of the uncorrelated lognormal relaxed clock was close to zero, the final analyses were run enforcing the strict molecular clock model. Using PartitionFinder, all codon positions being treated together as one partition was selected as the best-fitting partitioning scheme, and the HKY substitution model was selected as the best-fitting model. The final BSP analysis was run in duplicates to check for consistency between runs, each run for 10 million generations and sampled every 1000 generations. Convergence, ESS > 200, stationarity and the appropriate number of generations to be discarded as burn-in (10%) were assessed using TRACER. The resulting BSP was also summarized in TRACER, with the maximum time as the median of the root height parameter. In addition, the mismatch distributions (MD) were calculated as the distributions of the observed pairwise nucleotide differences and the expected values under a growing- or declining-population model using DnaSP. The occurrence of historical demographic changes was assessed with neutrality-test statistics (Fu’s FS and Tajima’s D) implemented and calculated in DnaSP, with the estimation of the statistical significance using 10 000 coalescent simulations. Negative significant values of FS and D are expected when population expansion occurs (Tajima, 1989; Fu, 1997).
SPECIES DISTRIBUTION MODELLING
The occurrence dataset of the Greek frog (R. graeca) was constituted by 1067 points from field campaigns, citizen science databases (GBIF, 2020; Ueda, 2020), the literature or our own databases (Supporting Information, Table S4). The area of the modelling was based on the current species distribution (Šunje et al., 2018) along with implications for possible past range in the context of the Balkan Peninsula. The environmental data are based on bioclimatic variables from the CHELSA database (Karger et al., 2017) for the current climate and from the PaleoCLIM database for the past conditions, the Last Glacial Maximum (LGM) and the Mid-Holocene (Fordham et al., 2017; Brown et al., 2018), at 5 km pixel size. The spThin R package (Aiello-Lammens et al., 2015) and an occurrence thinner radius of 5 km were used for the minimization of any effect of sampling bias (Boria et al., 2014). The use of the USDM R package (Naimi et al., 2014) has been selected for calculation of the variance inflation factor (VIF) for the set of selected predictors and excludes the highly correlated variables from the set through a stepwise procedure (VIF values < 10; Dormann et al., 2013). The Wallace R package (Kass et al., 2018) was used for the modelling, allowing the fine-tuning (Hao et al., 2020) of the MaxEnt algorithm using the ENMeval R package (Muscarella et al., 2014). The ‘block’ geographical partition scheme of ENMeval was selected for all the analyses (Muscarella et al., 2014) because it split the dataset into four different independent datasets across longitudinal and latitudinal directions (Supporting Information, Table S5; Fig. S1) for training/validation. ENMeval allowed us to evaluate models using a geographical partitioning scheme and to ‘fine-tune’ two parameters of MaxEnt that affect model complexity and predictive power. These parameters are the regularization multiplier (RM) or beta values and the feature classes (FCs). The RM penalizes overly complex models, whereas the FCs are functions of the raw environmental data (Morales et al., 2017). All FCs (L = linear, Q = quadratic, H = hinge and P = product) were selected, and the RM was set between one and five with steps of 0.5, allowing for model complexity and model tuning. Basically, all predictor variable coefficients were shrunk progressively until some reached zero, when they dropped out of the model. Only those variables with the greatest predictive contribution remained in the mode (Supporting Information, Table S6). The model selection was based on the average test Area Under Curve value (avg.test.AUC) along with the lowest delta corrected Akaike information criterion (delta.AICc) as calculated for each model following the method of Warren & Seifert (2011). In total, 45 different models were built, run and tested. A threshold rule distinguishing area with ‘suitable’ vs. ‘unsuitable’ climatic conditions was applied; we adopted the commonly used ‘minimum training presence threshold’ (Radosavljević & Anderson, 2014), and a ‘stability’ map was also created by assembling all of the model projections, that is, depicting as ‘suitable’ only the pixels that were predicted as ‘suitable’ by all of the model projections. The final maps were designed in ArcGIS v.10.5 (Environmental Systems Research Institute, 2019).
COMPARISON WITH ENDEMIC SPECIES OF AMPHIBIANS AND REPTILES IN THE BALKANS
To compare the genetic diversity of R. graeca with other endemic species of amphibians and reptiles distributed in the continental Balkans (sensuSpeybroeck et al., 2016), we compared available mitochondrial protein-coding gene sequences that have a similar substitution rate (Johnson & Sorenson, 1998): Cytb and NADH dehydrogenase subunits II–IV (ND2, ND3 and ND4; overall, 631 sequences) for seven amphibians and 17 reptiles (however, data for two Algyroides species and Podarcis erhardii were unavailable; for details see Supporting Information, Tables S7 and S8). For 21 species, we calculated the nucleotide diversity (π), haplotype diversity (hd) and number of haplotypes using DnaSP. The threshold for describing genetic diversity was considered as follows (π): 0–1%, low level; 1–2.5%, medium; and > 2.5%, high. This diversity was also put into the context of approximate species ranges sensuSillero et al. (2014), Speybroeck et al. (2016). The area of the distribution range was estimated using QGIS (QGIS Development Team, 2020). Moreover, we gathered the expected time of divergence for these endemic species based on available published data (Supporting information, Table S8). We performed general linear model (GLM) analysis to investigate possible relationships between the size of the species distribution range and the genetic diversity and between the age of divergence and the genetic diversity. We tested both the haplotype (hd) and the nucleotide diversity π (expressed as a percentage). We used the GeoDa python application (Anselin et al., 2006) for normal ordinary least squares (OLS) and R (R Core Team, 2020) for GLM with Poisson family using the built-in GLM of core R, and we used the rsq R package (Zhang, 2020) for the calculation of R2 values for each combination.
RESULTS
PHYLOGENY AND PHYLOGEOGRAPHY
We obtained a total of 351 sequences (99 sequences from 16S, 101 from Cytb, 53 from COI and 98 from BDNF) of R. graeca and 121 outgroup sequences (Rana spp. and Pelophylax spp.), with no signal of contamination or sequences of nuclear genomic origin (Supporting Information, Table S1). Analyses of saturation tests for each gene separately did not show any evidence of saturation (data not shown). Together with outgroup sequences, a total alignment length of 2707 bp was analysed. The best-fitting partitioning schemes for each downstream analysis and the selected nucleotide substitution models are presented in the Supporting Information (Table S3). Both datasets (mtDNA and mtDNA + nuDNA) identified only two lineages with similar tree topologies after ML (lnL = −7985.72) and BI analysis (lnL = −8073.12; Figs 2, 3B). For the mtDNA datasets (with and without COI), the ML and BI showed similar topologies (see Supporting Information, Figs S2–S7), which were in accordance with the topologies of the complete dataset with lnL = −6727.10 and −6851.08 (mtDNA dataset) with all mitochondrial gene fragments and lnL = −4214.19 and −4345.71 (mtDNA dataset) with all mitochondrial gene fragments except COI, respectively. Unfortunately, for the BDNF marker, sequences were not available in GenBank for the species R. pyrenaica, R. arvalis, R. macrocnemis and R. asiatica, and the analyses were restricted to R. graeca, R. dalmatina and R. temporaria. For BDNF, ML and BI produced similar topologies, in which R. graeca appeared monophyletic (lnL = −1062.54 and −1111.07 for ML and BI, respectively), as for the mtDNA trees, but without any intraspecific differentiation.
All MCMC diagnostic metrics indicated that the iterations of BI analysis reached convergence and stationarity. The average standard deviation of split frequencies was 0.009 for the complete dataset and 0.007 for the mtDNA dataset (when this value approaches zero, the tree samples are more similar), and the plot of generation vs. log-likelihood of the data had a characteristic ‘white-noise’ morphology after burn-in. In addition, for all parameters, the potential scale reduction factor values were 1.000 in both datasets, and the minimum ESS values were > 200 for all parameters in both datasets. According to the chronophylogenetic analyses of the concatenated dataset (lnL = −6917.64, ESS > 1070), the temporal time of divergence of the European Rana was estimated as 7.9 Mya (Tortonion age in the Late Miocene), with the divergence of the R. temporaria group from its sister clade that includes R. graeca, R. macrocnemis and R. dalmatina. The divergence of R. graeca and R. dalmatina was estimated to have occurred during the Messinian Salinity Crisis (MSC) in the Messinian age (Late Miocene; 5.7 Mya), whereas the intraspecific differentiation of R. graeca was dated to the Late Pleistocene (0.27 Mya) (Fig. 3A).
The resulting topology of R. graeca represented an unexpected and very shallow divergence, with two main statistically supported lineages (I and II). The highest recorded p-distance between lineages of the concatenated dataset was 0.67%, and overall nucleotide diversity was only 0.46%. Lineage I had a nucleotide diversity of 0.410% (N = 82) and lineage II 0.077% (N = 17). Lineage I included specimens collected in the Dinarides (northern Albania, Bosnia and Herzegovina), Hellenides (southern Albania and northern Greece), the Macedonian–Thracian Massif (Rhodopes; eastern Bulgaria and north-eastern Greece) and Peloponnese. Lineage II was distributed in the central and southern Hellenides and northern Peloponnese. Lineages I and II were recorded in sympatry only in north Peloponnese (gr46–gr49; Fig. 2).
Haplotype networks showed a different level of differentiation depending on the DNA marker used (the highest level was in mitochondrial Cytb and the lowest in nuclear BDNF; Fig. 4). Both lineages were well recognized in Cytb, containing unique haplotypes that were separated from the central and most common haplotype by two mutation steps (Fig. 4). The geographically marginal populations from Bosnia and Herzegovina or Bulgaria were detected in the central haplotype of lineage I.
Seventeen haplotypes were recognized in 107 sequences of 16S, with hd = 0.53 and π = 0.19% (Table 1). Both lineages were revealed, with only two haplotypes from southern mainland Greece and the northern Peloponnese belonging completely to lineage II (Fig. 4).
Summary of DNA polymorphism for concatenated dataset of Rana graeca and particular genes
| Length (bp) | N | h | S | π ± SD (%) | hd ± SD | θW ± SD (%) | |
|---|---|---|---|---|---|---|---|
| Rana graeca | 2707 | 99 | 33 | 27 | 0.46 ± 0.027 | 0.93 ± .0.011 | 0.74 ± 0.23 |
| 16S | 387 | 107 | 17 | 15 | 0.19 ± 0.025 | 0.53 ± 0.055 | 0.75 ± 0.26 |
| Cytb | 322 | 123 | 19 | 15 | 0.78 ± 0.054 | 0.87 ± 0.017 | 0.87 ± 0.30 |
| COI | 631 | 53 | 9 | 18 | 0.23 ± 0.059 | 0.70 ± 0.045 | 0.64 ± 0.23 |
| BDNF | 453 | 99 | 2 | 1 | 0.009 ± 0.006 | 0.04 ± 0.027 | 0.04 ± 0.04 |
| Length (bp) | N | h | S | π ± SD (%) | hd ± SD | θW ± SD (%) | |
|---|---|---|---|---|---|---|---|
| Rana graeca | 2707 | 99 | 33 | 27 | 0.46 ± 0.027 | 0.93 ± .0.011 | 0.74 ± 0.23 |
| 16S | 387 | 107 | 17 | 15 | 0.19 ± 0.025 | 0.53 ± 0.055 | 0.75 ± 0.26 |
| Cytb | 322 | 123 | 19 | 15 | 0.78 ± 0.054 | 0.87 ± 0.017 | 0.87 ± 0.30 |
| COI | 631 | 53 | 9 | 18 | 0.23 ± 0.059 | 0.70 ± 0.045 | 0.64 ± 0.23 |
| BDNF | 453 | 99 | 2 | 1 | 0.009 ± 0.006 | 0.04 ± 0.027 | 0.04 ± 0.04 |
Abbreviations: h, number of haplotypes; hd, haplotype diversity; N, number of species; π, nucleotide diversity; S, number of segregating sites; SD, standard deviation; θW, Watterson’s theta per site.
Summary of DNA polymorphism for concatenated dataset of Rana graeca and particular genes
| Length (bp) | N | h | S | π ± SD (%) | hd ± SD | θW ± SD (%) | |
|---|---|---|---|---|---|---|---|
| Rana graeca | 2707 | 99 | 33 | 27 | 0.46 ± 0.027 | 0.93 ± .0.011 | 0.74 ± 0.23 |
| 16S | 387 | 107 | 17 | 15 | 0.19 ± 0.025 | 0.53 ± 0.055 | 0.75 ± 0.26 |
| Cytb | 322 | 123 | 19 | 15 | 0.78 ± 0.054 | 0.87 ± 0.017 | 0.87 ± 0.30 |
| COI | 631 | 53 | 9 | 18 | 0.23 ± 0.059 | 0.70 ± 0.045 | 0.64 ± 0.23 |
| BDNF | 453 | 99 | 2 | 1 | 0.009 ± 0.006 | 0.04 ± 0.027 | 0.04 ± 0.04 |
| Length (bp) | N | h | S | π ± SD (%) | hd ± SD | θW ± SD (%) | |
|---|---|---|---|---|---|---|---|
| Rana graeca | 2707 | 99 | 33 | 27 | 0.46 ± 0.027 | 0.93 ± .0.011 | 0.74 ± 0.23 |
| 16S | 387 | 107 | 17 | 15 | 0.19 ± 0.025 | 0.53 ± 0.055 | 0.75 ± 0.26 |
| Cytb | 322 | 123 | 19 | 15 | 0.78 ± 0.054 | 0.87 ± 0.017 | 0.87 ± 0.30 |
| COI | 631 | 53 | 9 | 18 | 0.23 ± 0.059 | 0.70 ± 0.045 | 0.64 ± 0.23 |
| BDNF | 453 | 99 | 2 | 1 | 0.009 ± 0.006 | 0.04 ± 0.027 | 0.04 ± 0.04 |
Abbreviations: h, number of haplotypes; hd, haplotype diversity; N, number of species; π, nucleotide diversity; S, number of segregating sites; SD, standard deviation; θW, Watterson’s theta per site.
In 123 sequences of Cytb spaced across the Balkans, we found 19 haplotypes, with hd = 0.87% and π = 0.78% (Table 1). Fifteen out of 19 haplotypes corresponded to lineage I, and lineage II contained four haplotypes from the central and southern Greece and the northern Peloponnese (Fig. 4).
Only nine haplotypes were revealed among the 53 sequences of COI, with hd = 0.70 and π = 0.23% (Table 1). Only lineage I was distinct in the network (Fig. 4), and lineage II was observed sharing haplotypes with the lineage across the Greek part of the Hellenides (Fig. 4).
In the highly conservative nuclear marker (BDNF), only two haplotypes were revealed among the 99 sequences analysed, with very low hd (0.04) and π (0.009%) (Table 1). One was present in both lineages, distributed from south-eastern Bosnia and Herzegovina to the southern Peloponnese. The only haplotype that corresponded completely with lineage I was found in one locality in north-eastern Greece (Fig. 4).
DEMOGRAPHY
The Bayesian skyline plot analysis on Cytb (Fig. 5A) showed a slight population growth (Ne) of R. graeca, which started at ~15 kya, close to the end of the last glacial period. The observed value of complementary mismatch distributions (MDs; Fig. 5B) mirrored the expected value for a growing- or declining-population model. The neutrality tests showed no significant evidence of an expansion (FS = −0.18520, P = 0.549; Tajima’s D = −0.067, P = 0.445).
![A, historical demography of Rana graeca estimated with Bayesian skyline plots. The central line shows the mean value of the population size [Ne × τ × μ; where Ne is the effective population size, τ is the generation length in units of time (substitutions/site), and μ is the mutation rate] on the logarithmic scale, and the shaded area represents the 95% highest posterior density. LGM indicates the time of the Last Glacial Maximum. B, mismatch distributions of observed frequencies (dashed line) compared with the expected frequencies (continuous line) under the expansion model of population size.](https://oup.silverchair-cdn.com/oup/backfile/Content_public/Journal/biolinnean/132/4/10.1093_biolinnean_blaa224/2/m_blaa224_fig5.jpeg?Expires=1748021184&Signature=T6Y7Yg1jh5E2x2Vh8jz0h-k2hCYUPSaqKUJEopcjXzlR-mn3pW2osCwylnKiE8eEpo0Lg8JLVfLPAMDqgbu0Gai6KMcev32pviai3j4fTxyDoj22FqMBfhonF8jfVZOr5zPTiicIcj1ei2lIQe6TGgKIcPwzSsDfmp6IzrrecPBtiWPopf3CECcMPevYrudFGqMKOE0Wjlt4xCYrvbyT2npoi8al2M1JcqiBIPgyh76~7dIwU4g-Qe6iBnQORnxeaEd1R32pzHnrEjfiUBSbh2j~o3c5Ne4NGDqLkaihrgdf3GQgVmhHgmp6M5psbRk-f~p-IPTpAQUsY428AEA2Sw__&Key-Pair-Id=APKAIE5G5CRDK6RD3PGA)
A, historical demography of Rana graeca estimated with Bayesian skyline plots. The central line shows the mean value of the population size [Ne × τ × μ; where Ne is the effective population size, τ is the generation length in units of time (substitutions/site), and μ is the mutation rate] on the logarithmic scale, and the shaded area represents the 95% highest posterior density. LGM indicates the time of the Last Glacial Maximum. B, mismatch distributions of observed frequencies (dashed line) compared with the expected frequencies (continuous line) under the expansion model of population size.
SPECIES DISTRIBUTION MODELLING
From the initial 1067 occurrence data points (656 unique locations), 337 occurrence data points were included in the final model (Supporting Information, Tables S4 and S5). From the predictors after the VIF analysis, we ended up with eight out of 19 initial predictors; these were the isothermality (bio_3), temperature seasonality (bio_4), maximum temperature of the warmest month (bio_5), mean temperature of the wettest quarter (bio_8), mean temperature of the driest quarter (bio_9), precipitation seasonality (bio_15), precipitation of the warmest quarter (bio_18) and precipitation of the coldest quarter (bio_19) (Supporting Information, Table S6). Forty-five models were run consecutively, and the model with the lowest delta.AICc was the one that used hinge (H) features along with a regularization multiplier of two (H_2). The AUC was high for both train. AUC (0.9) and avg.test.AUC (0.84); it is worth mentioning that almost all models had high AUC values (> 0.85). The projections in the present, LGM and Mid-Holocene climate periods showed a stabile overall spatial coverage (where the species might occur), but with a vast reduction of its distribution in the Mid-Holocene (with appearance in some new areas, i.e. east Greece and European Turkey) and a slow spread in the current period (Fig. 6). The congruence of models showed that the species distribution was confined to the territory of the Balkan Peninsula, with some new potential areas of occurrence, i.e. where the presence of the species has not yet been confirmed (the Dinarides in Croatia, the Balkanides in Bulgaria or European Turkey). The most stable conditions across all models were detected in the south-western Balkans and the Peloponnese, corresponding to the Hellenides. The most important variables were bio_3, bio_44, bio_18 and bio_19, related to temperature and to rainfall (Fig. 6; Supporting Information, Table S6).
The environmental niche model projection on climatic conditions at the Last Glacial Maximum (LGM), Mid-Holocene, the present and a congruence of all models (green colour) for Rana graeca among geological times. White circles represent georeferenced occurrence records obtained from our own field research, publications and citizen science databases. The maps were designed in ArcGIS v.10.5 using Min-Max as the stretching histogram and bilinear interpolation as smoothing for display needs only.
Rana graeca and other endemics
Our comparative dataset of 24 endemic species of amphibians and reptiles from the continental Balkan Peninsula showed that R. graeca had the second-lowest value of genetic/nucleotide diversity (19 haplotypes with π = 0.78%) among the endemic amphibians and sixth lowest of all species (Supporting Information, Table S8), with the expected time of divergence dated at 5.7 Mya. In contrast, R. graeca was the species with the largest distribution range among the examined endemic taxa (~220 000 km2). The available published data of expected times of divergence showed that speciation of inspected endemic taxa of the Balkans occurred between ~16 and 3 Mya (the Miocene to Pliocene). For three species, we are missing comparative molecular data. Finally, 15 of those endemic species are biogeographically associated mostly with the Hellenides, but only six with the Dinarides and three (including the focus species) with both of them (Supporting Information, Table S8). No statistical evidence has emerged using OLS or GLM with Poisson family (R2 < 0.005), and it cannot support any relationships between the examined variables.
DISCUSSION
EVOLUTIONARY HISTORY
Five species of the genus Rana (R. arvalis, R. dalmatina, R. graeca, R. latastei and R. temporaria), from two different evolutionary groups (Veith et al., 2003; Yuan et al., 2016), are found in the Balkan Peninsula (Ilić et al., 2016; Speybroeck et al., 2016). Our knowledge of their evolutionary history is, however, limited in comparison to that of other amphibians in the region (e.g. Ichthyosaura, Triturus, Bombina, Bufotes, Pelobates and Pelophylax; Lymberakis et al., 2007; Fijarczyk et al., 2011; Wielstra et al., 2012, 2013; Recuero et al., 2014; Vucić et al., 2018; Dufresnes et al., 2013, 2019a, b). The reason is that there was less sampling effort in previous studies, which did not allow the resolution of relationships among the Balkans and other populations of Rana species from Europe (Babik et al., 2004; Vences et al., 2013; Canestrelli et al., 2014; Dufresnes et al., 2020). Therefore, this needs further attention. As our data suggest, genetic variability among members of the genus in the Balkans is expected to be higher than is currently known (see R. dalmatina and R. temporaria; Fig. 3; Supporting Information, Figs S2–S7).
In our phylogenetic analyses, R. graeca is a monophyletic species showing close relationships to R. dalmatina and R. macrocnemis. This is consistent with previously published studies, in which R. graeca, R. dalmatina, R. macrocnemis and R. latastei form a monophyletic clade (Picariello et al., 2002; Veith et al., 2003; Yuan et al., 2016; Dufresnes et al., 2020). This clade has a sister relationship with the R. temporaria species group that includes R. arvalis, R. pyrenaica, R. parvipalmata, R. iberica and R. italica. From the chronophylogenetic point of view, our estimated time of divergence supported the split of the European Rana in the Late Miocene (7.9 Mya; Fig. 3A). This is relatively close to the estimated radiation presented by Veith et al. (2003), but much younger than estimations by Yuan et al. (2016). According to Yuan et al. (2016), the European Rana diverged in the early Miocene (19 Mya), whereas the time differentiation of R. graeca from its sister lineages (R. dalmatina) was in the middle Miocene (12 Mya). This major discrepancy is probably attributable to the different calibration strategy (species tree approach vs. gene tree approach) used in our study compared with the study of Yuan et al. (2016). Moreover, our calibration was based on previous age estimations vs. the combination of fossils and age estimations used by Yuan et al. (2016). However, our estimations seem more reliable because they are closer to the MSC, confirming the current hypothesis about numerous speciation events during the late Miocene period (e.g. Dufresnes et al., 2018, 2019a, b). Thus, R. graeca represents an old lineage, probably with a Miocene common ancestor, that became the endemic brown frog of the Balkan Peninsula, similar to the way in which R. italica, R. iberica, R. parvipalmata, R. pyrenaica or R. tavasensis became endemics for Apennine, Iberian or Anatolian Peninsulas (Veith et al., 2003; Yuan et al., 2016; Dufresnes et al., 2020).
PHYLOGEOGRAPHY
Rana graeca showed a low level of intraspecific genetic diversity in view of its large range size (the largest among other endemic amphibians and reptiles of the Balkans), as shown by the concatenated markers and by the single markers (Figs 3, 4). In general, the two genetically shallow lineages revealed in R. graeca showed a certain phylogeographical pattern, whereby lineage I was distributed in the major part of the species range, whereas lineage II was limited to central and southern parts (Greece). Concurrently, these lineages formed a partly sympatric distribution pattern (Figs 2, 4).
In the richest dataset of 123 sequences of Cytb throughout the Balkans, we found 19 haplotypes distributed from the southern Peloponnese to north-western Bosnia and Herzegovina and western Bulgaria, also with very low intraspecific genetic variation. This variation was even lower in the case of the remaining analysed datasets (16S, COI and BDNF; Fig. 4). In the Cytb network analysis, two lineages of the concatenated dataset were recognized, separated by two mutation steps and partly suggesting a star-like pattern. Despite extensive sampling, our data did not show any significant hotspot of deeper genetic diversity suggestive of geographically limited microrefugia, as was presented for other species of the genus in Europe (Canestrelli et al., 2008, 2014; Vences et al., 2013; Teixeira et al., 2018; Dufresnes et al., 2020). The observed genetic pattern was, in contrast, very similar to that of R. pyrenaica (representative of the mountain Pyrenean endemism), suggesting very rapid colonization of its present small range (Carranza & Arribas, 2008). Likewise, the lineages of R. graeca probably diverged in a very narrow spatial and temporal window with subsequent rapid colonization, and possible extinction of peripheral populations, resulting in the observed low divergence between detected lineages and haplotypes (p-distances ≤ 0.62%). This corresponds to the so-called ‘R’ species model (Recuero & García-París, 2011) characterized by reduced genetic diversity, which was also observed in some other European brown frogs (Canestrelli et al., 2008; Vences et al., 2013; Teixeira et al., 2018). Although low genetic diversity accompanied by a star-like phylogeographical structure often predict past population growth (Slatkin & Excoffier, 2012), this was not significantly supported here by neutrality tests (Fu’s FS and Tajima’s D). Therefore, we can hypothesize that the Hellenides, as such, were probably the major centre for the species divergence and further evolution and colonization (Figs 2–4, 6), supporting the biogeographical importance of these south European mountains (Poulakakis et al., 2015; Jablonski et al., 2016; Marzahn et al., 2016; Psonis et al., 2017, 2018).
Šunje et al. (2018) suggested that dispersal of R. graeca in the Pleistocene was probably ongoing during warm phases and was restricted by climatic oscillations and geographical barriers, e.g. mountain glaciers. Surprisingly, however, our data suggested the opposite results and are in concordance with the presumption of Teixeira et al. (2018) that Mediterranean populations of Rana are currently surviving in interglacial refugia. This is unexpected, because the Plio-Pleistocene climatic changes are considered to be factors behind the diversified phylogeographical patterns observed today (Fijarczyk et al., 2011; Pabijan et al., 2015; Jablonski et al., 2016; Dufresnes et al., 2019b). First, based on genetic data (Cytb), the colonization of northern, eastern and southern areas, outside of the Hellenides, was probably conducted by several haplotypes very close to those from the central range (Fig. 5). More precisely, northern populations (Bosnia and Herzegovina, Montenegro and Serbia) were formed by three haplotypes that are close to those occurring in Albania or Greece (Fig. 4). Second, the SDM projection showed a possible range expansion of the species during the LGM and contraction during the warmest middle Holocene, probably followed by another expansion that was subsequently observed in the present model and corresponds to the known species range (Figs 2, 6).
This is also supported by demographic analysis, showing slight growths in population from around the LGM until today, but no significant population expansion. Similar results have been observed in several taxa (e.g. Stewart et al., 2010), for which climatic conditions during the LGM were probably more suitable for population growth and/or range expansion compared with the warmer periods. Such a pattern was observed in very few cases of western Palaearctic amphibians (mainly cold-adapted species; Teixeira et al., 2018; Afroosheh et al., 2019). These species occur in lower-elevation areas during the LGM and then shift to higher altitudes. However, this cannot be fully applicable for R. graeca, because this frog is currently known from sea level or even thermal springs (Bringsøe, 2011; Šunje et al., 2017; D. Jablonski, personal observations) up to high mountains (e.g. Szabolcs et al., 2017). The observed long-term stability in the evolutionary history of R. graeca located in the south-western part of the peninsula therefore represents a new view of the phylogeography of the Balkan fauna with the following conclusions: (1) R. graeca is currently associated with regions of main mountain ranges of the Balkans, i.e. lineage I mostly in the Dinarides, Macedonian–Thracian Massif, northern and central Hellenides and lineage II in the southern Hellenides, including Peloponnese; (2) the species shows relative genetic homogeneity, with one geographically widespread lineage and a second, probably sympatric lineage inside the species range; and (3) the morphological and ecological differentiation observed in the species (Bringsøe, 2011; Asimakopoulos & Grossenbacher, 2014) is not mirrored in a phylogeographical signal of the markers studied. This means that R. graeca did not diversify in deeper phylogenetic lineages that could have a connection to certain historical microrefugia and/or be mirrored in the observed phenotypic variability. This suggests that the main influences were probably caused by microclimatic and local environmental conditions rather than past, long-term climatic cycles. Concurrently, the morphological variability (colour and pattern variation) and environmental adaptability of the species (type of the habitat, wide elevational range) are probably attributable to the phenotypic and ecological plasticity of the species, a fact which calls for further research (see also Šunje et al., 2017).
Endemism in the Hellenides
Rana graeca is one of 24 endemic species of amphibians and reptiles in the continental Balkans (Supporting Information, Table S8). Although it has the largest distribution range among all endemic herpetofauna living there, its genetic diversity is very low. This could probably be attributed to the evolutionary history of the species (rapid colonization with a limited level of genetic variability), ecology or its overall natural history (the species mainly follows river valleys). Two phylogenetically supported lineages found in R. graeca have different sizes of the detected range, whereby the smaller (lineage II) is restricted to the southern Balkans and should be considered as a separate unit from the conservation point of view. As shown by our comparison (Supporting Information, Table S8), different levels of genetic diversity, range size or estimated time of divergence could be found among amphibians and reptiles endemic to the Balkans, especially in the Hellenides, e.g. cold-tolerant species showing a low level of genetic diversity in small areas of distribution (Vipera graeca; Mizsei et al., 2016, 2017) or species with larger distribution ranges influenced by long-term climatic oscillations that resulted in a high level of genetic diversity (e.g. Anguis graeca, Mediodactylus kotschyi and Podarcis ionicus; Jablonski et al., 2016; Psonis et al., 2017, 2018; Kotsakiozi et al., 2018). We did not find relationship between the age of divergence and the level of currently observed genetic diversity or the size of the species range. However, this comparison is only general, and it needs further comprehensive research. Despite a wide current species distribution, the evolution of the specie′s is most probably associated mainly with the Hellenides. Despite the wide current distribution of R. graeca, the evolution of this frog has been most probably connected with the Hellenides. The Hellenides, partly shared with the Dinarides, represent 63% (15 species) of the overall continental Balkan herpeto-endemism (Supporting Information, Table S8) and are one of the most important biodiversity hotspots in the Eastern Mediterranean. This is definitely a topic for further evolutionary research and conservation priorities, in addition to the presented view of these southern European mountains as sources of species, genetic, or ecological diversity.
SUPPORTING INFORMATION
Additional Supporting Information may be found in the online version of this article at the publisher’s web-site:
Table S1. List of samples and sequences used for analyses.
Table S2. Primers and conditions used in polymerase chain reaction amplifications and in cycle sequencing reactions.
Table S3. Partitioning schemes and best-fitting models of sequence evolution selected in PartitionFinder (PF) for downstream analyses.
Table S4. Complete occurrence data for species distribution modelling.
Table S5. Selected and partitioned data for species distribution modelling.
Table S6. Eight main variables forming the species distribution modelling of the species.
Table S7. GenBank numbers of sequences used for DNA polymorphism of endemic species.
Table S8. Overview of DNA polymorphism and genetic diversity of all species of amphibians and reptiles endemic in the continental Balkans, based on mitochondrial DNA, and their species ranges (in square kilometres). Abbreviations: h, number of haplotypes; hd, haplotype diversity; N, number of species; π, nucleotide diversity. *The source with an estimated time of divergence.
Figure S1. The area of work, covering almost the whole of the Balkan Peninsula. Using coloured dots, the occurrence data of the Greek frog (Rana graeca) are split using the block method.
Figure S2. Phylogenetic relationships of Rana reconstructed using Bayesian inference of 16S, Cytb and COI sequences. The numbers above the branches represent Bayesian posterior probabilities showing the branch support. For details, see the Supporting Information (Table S1).
Figure S3. Phylogenetic relationships of Rana reconstructed using maximum likelihood of 16S, Cytb and COI sequences. The numbers above the branches represent bootstraps showing the branch support. For details, see the Supporting Information (Table S1).
Figure S4. Phylogenetic relationships of Rana reconstructed using Bayesian inference of 16S and Cytb sequences. The numbers above the branches represent Bayesian posterior probabilities showing the branch support. For details, see the Supporting Information (Table S1).
Figure S5. Phylogenetic relationships of Rana reconstructed using maximum likelihood of 16S and Cytb sequences. The numbers above the branches represent bootstraps showing the branch support. For details, see the Supporting Information (Table S1).
Figure S6. Phylogenetic relationships of Rana reconstructed using Bayesian inference of BDNF sequences. The numbers above the branches represent Bayesian posterior probabilities showing the branch support. For details, the Supporting Information (Table S1).
Figure S7. Phylogenetic relationships of Rana reconstructed using maximum likelihood of BDNF sequences. The numbers above the branches represent bootstraps showing the branch support. For details, see the Supporting Information (Table S1).
ACKNOWLEDGEMENTS
We would like to thank P. Balej, V. Gvoždík, D. Jandzik, M. Meszáros, P. Mikulíček, E. Mizsei, S. Papežíková, J. Poláková, M. Raffaj, E. Themeli and B. Vági for donation of samples, technical support or help in the field. We are also grateful to M. K. Lawson, S. R. Goldberg (who helped us with English) and Emina Šunje, in addition to three anonymous reviewers, for their comments, which significantly improved the first submitted version of the manuscript. This work was supported by the Slovak Research and Development Agency under the contract numbers APVV-15-0147 and APVV-19-0076.
REFERENCES
