Genetic diversity of the Caucasian Parsley Frog, Pelodytes caucasicus (Anura: Pelodytidae)

Caucasia is a global biodiversity hotspot, rich in amphibians, including several endemic species. We sequenced samples from Parsley frogs (genus Pelodytes) across their Anatolian range to generate a barcode reference database and to assess patterns of genetic diversity in the species. Different species delimitation methods (ABGD, ASAP, GMYC and PTP) were applied to assess species diversity in the genus Pelodytes based on published and newly obtained mtDNA sequences. A majority of the species delimitation tests (ABGD, GMYC and ASAP) recovered four taxonomic units corresponding to currently accepted taxonomy (P. atlanticus, P. caucasicus, P. ibericus and P. punctatus). PTP, on the other hand, recovered only two taxonomic units, one combining the three Iberian taxa (P. atlanticus, P. ibericus, and P. punctatus), and the other, P. caucasicus. In Anatolia, individuals from Giresun and Trabzon were found to be genetically closer to each other compared to those from Rize and Artvin, based on genetic distances and phylogenetic and haplotype network analyses.

. Locations from where Pelodytes material was studied. Sample numbers correspond to those given Table S1 (red diamond: own data; red dots: data derived from BOLD). across all species of Pelodytes, including the Iberian taxa. This study additionally assesses species diversity in the genus using various species delimitation methods applied to COI sequences.

Material and Methods
Material. We analysed 41 COI sequences belonging to all species in the Pelodytidae family. Of these, 15 sequences were newly generated for this study from specimens collected in the provinces of Giresun, Trabzon, Rize, and Artvin (Turkey) during the reproductive period (1-30 July 2020 and 1-30 July 2021). The remaining 26 sequences were retrieved from the Barcode of Life Data Systems (BOLD) database (Table S1, Figure 1 DNA extraction and sequencing. The total DNA was extracted from tissue samples using a Nucleospin DNA isolation kit (Macherey Nagel, Germany). A universal primer mix (comprising of equal proportions of forward primers: COI-C01 and COI-C02, and reverse primers COI-C03 and COI-C04), designed by Che et al. (2012), was used for polymerase chain reactions (PCRs). Only COI-C01 and COI-C03 primers were used for sequencing. The following PCR protocol was used for amplification: initial denaturation for 5 min at 95°C; 36 cycles of denaturation for 1 min at 95°C, annealing for 1 min at 46°C, an extension for 1.15 min at 72°C; and final extension for 10 min at 72°C. Data analyses. The sequences were checked and edited using the software BioEdit (Hall, 1999) and then trimmed to 642 bp. The new sequences were submitted to both GenBank and BOLD systems (Table S1). The intra-and inter-specific genetic distances were calculated with Mega X (Kumar et al., 2018), based on p-distance parameters (Table 1). DNASP (Librado & Rozas, 2009) was then used to calculate the haplotype diversity (Hd), and nucleotide diversity (pi). The operational taxonomic units (OTUs) were additionally inferred through Refined Single Linkage (RESL) analyses on BOLD (Ratnasingham & Hebert, 2013).
Four different species delimitation methods were used to assess the species diversity in Pelodytidae: Poisson-Tree-Processes (PTP), Generalized Mixed Yule Coalescence (GMYC), Automatic Barcoding Gap Discovery (ABGD) and Assemble Species by Automatic Partitioning (ASAP). Although a multiple-threshold version of GMYC has been developed (Monaghan et al., Table 2. Inter-specific p-distances between individuals of different species of the genus Pelodytes based on COI sequences. P-cau = P. caucasicus; P-atlan = P. atlanticus; P-iber = P. ibericus; Ppunc = P. punctatus.

P-cau (Artvin)
- 2009), only the single-threshold version of GMYC was used in this study, which has been shown to outperform the multiple-threshold version (Fujisawa & Barraclough, 2013;Talavera et al., 2013). From each of the different species delimitation methods, the number of Operational Taxonomic Units (OTUs) was derived. For PTP species delimitation, a maximum likelihood phylogenetic tree was used, which was reconstructed with PhyML (Guindon et al., 2010) with 1000 bootstrap replications (Zhang et al., 2013). The analysis was performed on the PTP web server (http://species.h-its.org and http://mptp.h-its.org/#/tree). An ultrametric tree was generated using BEAST 1.8.2 (Drummond et al., 2012) before conducting the GMYC analysis. The analysis was carried out for 50 million generations with a sampling frequency of 10,000. All runs were checked with Tracer 1.6 to assess convergence. The maximum clade credibility tree was reconstructed using TreeAnnotator 1.8.2 (Drummond et al., 2012). The resulting ultrametric tree was imported into R 3.1.3 (R Core Team, 2013), and the single threshold ST-GMYC analysis was carried out using the R packages 'splits' (Ezard et al., 2009) and 'ape' (Paradis, 2004). ABGD (Puillandre et al., 2012) was used to test for a barcode gap in the COI sequences of the Pelodytidae (Relative gap = 1, Pmin=0.001, Pmax=0.1, Steps=10, NBins=10, K2P). ASAP is a recently developed method that infers the number of OTUs by implementing a hierarchal clustering algorithm based on pairwise genetic distances from single-locus sequence alignments (Puillandre et al., 2021).

Results
Among the 642 bp sequences which were analysed, 490 sites were conserved and 152 were variable. The average nucleotide frequencies were: G (17.0%), C (23.5%), T (33.7%), and A (25.7%). The overall transition/transversion ratio (R) was 2.1. The sequences of P. caucasicus were recovered as a single OTU according to cluster analysis (RESL), with Iberian species comprising three additional OTUs. The haplotype diversity (Hd) was found to be 0.9 and nucleotide diversity (Pi) 0.008. Pairwise genetic distances (p-distances) ranged from 0.002 to 0.015 in P. caucasicus. The highest p-distances were found between sequence pairs from Russia-Artvin, Russia-Rize, Trabzon-Artvin and Trabzon-Rize. Minimum p-distances were found between sequence pairs from Giresun-Trabzon, Giresun-Russia and Artvin-Rize (Table 1). Accordingly, the specimens from Giresun, Trabzon and Russia were found to be more closely situated than specimens from Rize and Artvin in phylogenetic trees and haplotype networks (Figures 2-3). In Iberian Pelodytes, the maximum p-distances were 0.04 for P. punctatus-P. ibericus and P. ibericus-P. atlanticus, and 0.03 for P. punctatus-P. atlanticus.
The PTP analysis recovered two OTUs in Pelodytes: one comprising the three Iberian species combined and the second, all sequences of P. caucasicus. In contrast, ABGD, GMYC and ASAP analyses recovered four OTUs corresponding to the currently recognised species: P. caucasicus and P. ibericus, P. atlanticus and P. punctatus. The GMYC clusters were delimited using the single threshold model, and the likelihoods of the null and GMYC models were 260.6498 and 262.4722, respectively. The total GMYC analysis was represented by four ML entities (CI=1-12). The ABGD gap analysis recovered four groups of sequences corresponding to the currently accepted Pelodytes species and revealed a clear gap between intra-and inter-specific distances, between 0.05 and 0.22, respectively. According to the detected prior intraspecific divergence (P) between 0.001 and 0.01, 2, 4 and 9 OTUs were recognised in initial partitions and recursive partitions, respectively. This is because the initial partitions were stable on a wider range of prior values and usually close to the hypothetical species number described by taxonomists (Puillandre et al., 2012). ASAP also identified four OTUs (ASAP-Score = 1, threshold distance = 0.02) ( Figure S1).

Discussion
Moderate levels of genetic diversity were observed in P. caucasicus, with nine different haplotypes in the 19 individuals analysed, in addition to evidence for regional population differentiation.
The abundance of P. caucasicus across its range is limited by the lack of appropriate breeding sites. One of the key threats for the species is habitat loss and fragmentation associated with forest logging (Litvinchuk & Kidov, 2018). Further studies about the genetic structure of P. caucasicus should be undertaken in order to assess general patterns of population connectivity at finer spatial scales. These studies should incorporate additional molecular markers from the nuclear genome to obtain more further robust assessments of current patterns of genetic diversity and structure in the species.
A majority of the results of the species delimitation analyses carried out in this study are consistent with the results of Dufresnes et al. (2020), which, based on hybrid zone analyses, support the existence of three Iberian taxa: P. atlanticus, P. punctatus and P.
ibericus. An exception was PTP, which tends to provide more conservative results than other delimitation tests (Zhang et al., 2013). In contrast, the results of this study support a single taxonomic unit in P. caucasicus. The genetic distances are consistent with this taxonomic arrangement, with evidence for a "barcoding gap" separating intra-and interspecific genetic distances.

Supplementary Material
Supplementary Material is given as a Supplementary Annex, which is available via the "Supplementary" tab on the article's online page.