We present the first molecular analysis of the Bufo bufo (Linnaeus, 1758) complex population in the Levant. This complex contains four species distributed through the Western Palearctic. Based on the analysis of the widely used mitochondrial marker 16S rRNA, populations from the Levant belong to B. verrucosissimus (Pallas, 1814) and have a close affiliation to populations from the Transcaucasian region rather than southern Turkey. Such identified relationships may imply past rapid colonisation from north to south and support the need for further research.

common toad, Bufonidae, mtDNA, biogeography, taxonomy, Middle East

The Caucasian Toad, Bufo verrucosissimus (Pallas, 1814), is a member of the Western Palearctic B. bufo (Linnaeus, 1758) complex along with three other species: B. bufo, B. eichwaldi and B. spinosus (García-Porta et al. 2012; Recuero et al. 2012; Arntzen et al. 2013). Some authors consider B. verrucosissimus as a subspecies of B. bufo due to overlapping distribution, allele-sharing and possible hybridisation events (García-Porta et al. 2012). However, individual-based analyses of enzyme data, presented by Arntzen et al. (2013), support B. verrucosissimus as a full species until better evidence for its taxonomic position is presented. This species is larger in size than other members of the complex (70–190 mm snout-vent length), has a rugged skin texture, different colouration and marked sexual dimorphism in body size (Litvinchuk et al. 2008). The Caucasian Toad prefers wet and shaded sites in mountain coniferous, mixed and deciduous forests (Dufresnes 2019) in the Transcaucasian region (Armenia, Azerbaijan, Georgia, southern Russia and north-eastern Turkey). However, recent molecular data also confirmed this species in western and southern Turkey (García-Porta et al. 2012; Arntzen et al. 2013).

Although information about B. bufo complex in the Levant region (Syria, Lebanon) is absent in recent phylogeographic studies (García-Porta et al. 2012; Recuero et al. 2012), its occurrence had been documented there in the past (Disi and Böhme 1996; Hraoui-Bloquet et al. 2001). Potential isolated distribution of the complex within southern Turkey to northern Israel was also demonstrated through paleoclimatic reconstructions by García-Porta et al. (2012; fig. 7, 125 p.). However, recent information concerning their distribution is very incomplete due to ongoing political turmoil and the existence of fragmented and low density populations of these toads in the region.

We, therefore, document the species’ status in the Levant, specifically in Lebanon. The presence of B. bufo complex in this country was first recorded from the locality Jannat-Kartaba near Nahr Ibrahim (Abraham River) valley (Hraoui-Bloquet et al. 2001; Suppl. material 1: Fig. S1) and represents the southernmost known distribution limit for the complex in the eastern part of the range. However, for nearly two decades, no information about these toads was reported in literature as they are rarely encountered and are restricted to several river valleys. Although Kuzmin (2008) expected the presence of B. verrucosissimus in Lebanon due to its close proximity to known populations in the Anatolian part of Turkey, the taxonomic status (possible independent phylogenetic clade), genetic affiliation and origin (possible existence of the relict population) of this toad population in the Levant are unknown. Therefore, we conducted the first molecular-based study analysing the species/population affiliations of this toad.

Total genomic DNA was extracted from eight blood samples (see Forzán et al. 2012 for details on blood sampling) collected on 29 March 2018 in one population from south-central Lebanon (Moukhtara, around Awali river, 33.66N, 35.61E, 689 m a.s.l.; voucher specimens were not collected; Suppl. material 1: Fig. S2), using the Qiagen DNeasy Blood and Tissue Kit, following the manufacturer’s protocols. The majority of sequences for the B. bufo complex in GenBank (Suppl. material 1: Table S1) was available for the mitochondrial DNA (mtDNA) fragment of the commonly used 16S ribosomal RNA (16S rRNA). We used this standard phylogenetic marker in amphibian DNA barcoding and species classification (Vences et al. 2012) for DNA amplification of the Lebanese population (primers 16Sar-L and 16Sbr-H, using the laboratory protocol of Palumbi et al. 1991). PCR products were sent to Macrogene Inc. (Amsterdam, The Netherlands) for sequencing. For basic species affiliation, the new sequences were checked via BLAST (https://blast.ncbi.nlm.nih.gov/Blast.cgi) and combined with available 16S rRNA GenBank sequences of B. verrucosissimus published in previous studies by Kutrup et al. (2006), Van Bocxlaer et al. (2009), Recuero et al. (2012), García-Porta et al. (2012), Litvinchuk et al. (2012) and Liedtke et al. (2017; Suppl. material 1: Table S1 and see the comment and the map of fig. 3, 1206 p. in Arntzen et al. 2013). The 16S rRNA fragment (~550 bp) was aligned using Clustal W algorithm (Thompson et al. 1994) as implemented in BioEdit (Hall 1999). The alignment was checked by eye and low quality ends trimmed. Ambiguously aligned regions/gaps were ignored for the analysis. We used a network approach (Posada and Crandall 2001) to infer inter-population relationships. We constructed a haplotype network for B. verrucosissimus and allopatric representatives of its related species (JN647044 – B. bufo, JN647239 – B. eichwaldi, JN647207 – B. spinosus) as outgroups, using the 95% limit of parsimony (TCS algorithm; Clement et al. 2000). This was implemented and visualised in the software PopArt (http://popart.otago.ac.nz). The DnaSP 5.10 (Librado and Rozas 2009) was used to estimate the number of haplotypes (h), haplotype diversity (Hd), number of segregating sites (S), average number of nucleotide differences (k), nucleotide diversity (π) and Watterson’s theta per site (θW) for all available sequences of the species. The same programme was used to estimate uncorrected p-distances between haplotypes. The eight novel sequences were deposited in GenBank under accession numbers MN094870-MN094877.

Overall 32 sequences (16S rRNA) covering the species range in Turkey and Transcaucasia were combined with eight new ones from Lebanon (Fig. 1A, B; Suppl. material 1: Table S1). The network analysis supports the inclusion of the Lebanese population to B. verrucosissimus which is the first molecular-based confirmation of this species in Lebanon and in the Levant. As indicated by the network analysis (Fig. 1B), the sequences obtained from eight specimens form three different haplotypes in the studied Lebanese locality. The maximal distance amongst these haplotypes is two mutational steps (Fig. 1B). These haplotypes are only one to three mutation steps distant from the geographically most common haplotype (e.g. FJ882807) of the species occurring in Georgia and Russia (0.2 to 0.6% in p-distances) and surprisingly two to five steps (0.4 to 1.0%) from haplotypes detected in localities of southern Turkey (JQ348811 from Alanya and AY840240-1 from Mersin; Fig. 1A, B). Our analysis showed 12 detected haplotypes overall, forming a star-like pattern with a nucleotide diversity of 0.34% (Fig. 1, Suppl. material 1: Table S2).

Figure 1. 

(A) Map displaying the localities of genetically investigated (16S rRNA) populations of Bufo verrucosissimus. The expected range of the species, according to Kutrup et al. (2006), García-Porta et al. (2012), Litvinchuk et al. (2012), Recuero et al. (2012) and Arntzen et al. (2013), is highlighted in orange. The newly investigated population is represented by a red circle with a black dot in the middle. Geographically closest (Alanya and Mersin, Turkey) and furthest (Kyurdzhinovo, Russia) investigated localities, relative to Lebanon, are marked by GenBank accession numbers. The question mark represents parts where species distribution needs further investigation. (B) Haplotype network inferred from a 16S rRNA showing the intraspecific relations of B. verrucosissimus. Black dots inside the red circle represent population from Lebanon and arrows following geographically closest and furthest localities given to Lebanon as presented in Fig. 1A. Sequences of other species of the B. bufo complex were used as the outgroup. The pictured specimen originates from the locality Moukhtara, Lebanon.

The Levant region is known as a speciation centre with endemic phylogenetic lineages and cryptic taxa hidden within morphologically similar sister species that also have fully or partially isolated distributions (e.g. Van Riemsdijk et al. 2017, Dufresnes et al. 2019a, b). Our results are, therefore, surprising in that the Levant population (i) does not form any significantly divergent lineage as would be expected from their likely isolated distribution there, (ii) is genetically very close to populations from the Transcaucasian region rather than the population from southern Turkey (Fig. 1). The northernmost (Kyurdzhinovo, Russia; JQ348764 and JX218095) and southernmost genetically investigated localities (Moukhtara, Lebanon) are separated by approximately 1300 km (straight line distance). However, the genetic distance between these localities is only one or two mutations steps (average 0.4%), similar to some other haplotypes inside the Transcaucasian region (Fig. 1). In the context of geographical distance, one mutation step between the most common Transcaucasian haplotype (accession number, e.g. FJ882807) and the haplotype from Alanya, Turkey (JQ348811) is also noteworthy. The star-like pattern of haplotypes thus suggests a probable origin in Transcaucasia and species radiation from northern to western and southern regions of the current range through Anatolia and the Levant. Such results imply most likely past rapid but different colonisation scenarios. As we used here only mtDNA barcoding, that does not resolve strictly maternal inheritance, the lack of recombination or possible massive mitochondrial introgression (Ballard and Whitlock 2004), denser sampling in Anatolia and Transcaucasia with a combination of genome-wide multilocus analyses and compilation of phenotypic datasets is needed to accurately delineate the phylogeography of B. verrucosissimus and origin of populations in the Eastern Mediterranean.

Lebanon is a mountainous country with deep, wet valleys that could regionally provide suitable microrefugia for B. verrucosissimus. This species prefers wet and forested areas (Kuzmin 2008) and has a clear preference for river basins with particularly rocky riparian habitats. However, related Bufotes toads are also found here but in various different habitats. Bufo verrucosissimus is probably less adapted to the drier climate of the Levant than genus Bufotes. This could be a reason for its rarity in the region and its restriction to the western slopes of Mt. Lebanon which are wetter than other regions in Lebanon. Additional distribution research is needed to identify other localities for B. verrucosissimus in the Levant (see Tron and Rocha 2005) in order to establish local conservation priorities for these southernmost populations of the species.

We thank Daniel Gruľa for his help with the fieldwork in Lebanon, Jana Polakova for laboratory work, Chace Holzheuser for language improving and anonymous reviewers for their important suggestions. This study was supported by the Slovak Research and Development Agency under contract no. APVV-15-0147.