Understudied regions and messy taxonomy: Geography, not taxonomy, is the best predictor for genetic divergence of the Poecilimon bosphoricus species group

The complex and dynamic history of the Anatolian Peninsula during the Pleistocene set the stage for species diversification. However, the evolutionary history of biodiversity in the region is shrouded by the challenges of studying species divergence in the recent, dynamic past. Here, we study the Poecilimon bosphoricus (PB) species group to understand how the bush crickets' diversification and the regions' complex history are coupled. Specifically, using sequences of two mitochondrial and two nuclear gene segments from over 500 individuals for a comprehensive set of taxa with extensive geographic sampling, we infer the phylogenetic and geographic setting of species divergence. In addition, we use the molecular data to examine hypothesized species boundaries that were defined morphologically. Our analyses of the timing of divergence confirm the recent origin of the PB complex, indicating its diversification coincided with the dynamic geology and climate of the Pleistocene. Moreover, the geography of divergence suggests a history of fragmentation followed by admixture of populations, suggestive of a ring species. However, the evolutionary history based on genetic divergence conflicts with morphologically defined species boundaries raising the prospects that incipient species divergences may be relatively ephemeral. As such, the morphological differences observed in the PB complex may not to be sufficient to have prevented homogenizing gene flow in the past. Alternatively, with the recent origin of the complex, the lack of time for lineage sorting may underlie the discord between morphological species boundaries and genetic differentiation. Under either hypothesis, geography—not taxonomy—is the best predictor of genetic divergence.


INTRODUCTION
Classically, an association between species richness and geographic and/or climatic factors (e.g., topographic complexity, geographic barriers and climatic stability) has been used to provide insights into the drivers of diversification (Ahmadi et al., 2021;Rahbek et al., 2019).
However, with the ubiquity of DNA sequencing, the association between geographic and climatic factors with the geography of genetic structure and divergence among constituent taxa has become a predominate framework in studies of the factors contributing to diversification.For groups with historically unstable taxonomies for which estimates of species richness vary widely, genetic data offers an especially invaluable window into the drivers of diversification because the genetic data can also serve as a tool for corroborating hypotheses of species boundaries (Knowles & Carstens, 2007;Tonzo et al., 2019).For example, differences in species richness in Hercules beetles across the Isthmus of Panama suggested an uptick in diversification in the beetles; however, genetic analyses showed that regional differences in beetle species richness reflected differing taxonomic practices rather than process affecting diversification rates (Huang & Knowles, 2016).This dual functionality of application of genetic data is likewise especially valuable for regions with complex geologic and climatic histories and/or regions with relatively little biogeographic investigations (Bellati et al., 2015;Riemsdijk van et al., 2017).For such understudied regions, a lack of information may limit a priori predictions about what factors might affect species divergence.Likewise, for areas associated with climate-induced distributional shifts, it is unclear that the extent to the consequences of such shifts for patterns of genetic divergence can be generalized across taxa.For example, while there is evidence that shifting species distributions (e.g., Knowles & Massatti, 2017;Marske et al., 2020) and/or species persistence along ice sheet margins or within glacial refugia (e.g., Bemmels et al., 2019;Çıplak, 2008;Hampe & Petit, 2005;Hewitt, 1996;Ortego & Knowles, 2022) might have promoted genetic divergence, the effects of such factors can be species-specific, complicating generalizations about any particular region (see Çıplak et al., 2015;Chobanov et al., 2017;Kaya & Çıplak, 2016Kaya & Çıplak, , 2017;;Massatti & Knowles, 2016;Papadopoulou & Knowles, 2016;Uluar et al., 2023;reviewed in McGaughran et al., 2022).
Here, we use a genetic framework to investigate the geography of divergence in a species-rich group of bush crickets (Poecilimon bosphoricus [PB] species group) from Anatolia and adjoining Caucasian and Balkans regions, a relatively understudied region with a complex geologic history (Bozkurt, 2001;Popov et al., 2004).The PB species group is the most species-rich (with 27 recognized species; Borissov et al., 2023;Chobanov et al., 2015;Kaya et al., 2012;Ünal, 2012) of the diverse genus Poecilimon, which contains more than 170 species arranged into 24 species groups (up from estimates of 68 species in eight groups; Ramme, 1933), with two species not assigned to any species group (Borissov et al., 2023).The genus is primarily distributed throughout the Anatolio-Balkan region of the Western Palearctic, with a few species also occurring in Italy, the Levant, the Southern Caucasus and other parts of the Black Sea Basin (Çıplak, 2004a).The PB species group has a ring-shaped distribution around the Black Sea and the Marmara Sea, with a distributional gap along the dry Romanian and adjoining Bulgarian and Ukrainian Coasts.The topographically complex landscape of the PB range area contrasts to the few narrow lowland areas along the coasts of the geologically young Black Sea and the Marmara Sea, which connects to the Aegean Sea in the southwest (Elmas, 2003;Gökaşan et al., 2010;Görür et al., 2000;Krijgsman et al., 2019).However, these coast lines were not stable over time because of climate-induced cycles of expansion and regressions of sea levels, which were pronounced at times.Consequently, what constitutes aquatic barriers today may have become less of a barrier with drops in sea level.Moreover, climate-induced distribution shifts during the Pleistocene may have affected population connectivity as taxa occupying this temperate environment (or Mediterranean climate to the southwest) track shifts in habitat availability (e.g., Ahmadi et al., 2021;Chobanov et al., 2017;Kaya & Çıplak, 2016).Although such distributional shifts have driven diversification, including in other orthopteran groups (e.g., Çıplak et al., 2015;Kaya et al., 2015;Knowles & Massatti, 2017;Ortego & Knowles, 2022;Uluar et al., 2023), the high diversity of the PB species group may also reflect the relative stability of biodiversity during the dynamic Pleistocene.Specifically, the ranges of the PB species group extend to below 52 degree latitude, the margin of glacial/permafrost during glacial cycles of Pleistocene in the Western Palearctic (Hewitt, 2000) and a potential Balkan/Anatolian/Caucasian glacial refugia (Çıplak, 2008;Çıplak et al., 2015;Korkmaz et al., 2014).
In addition to the species richness of the PB group, with its biodiversity centre in the geologically complex and understudied region of Anatolia, there are several other factors that make the PB species group a particularly compelling group to study diversification.The species are defined by a genitalic structure-the species-specific shape of the male cercus.Species description based on these cercal forms are also frequently accompanied by bioacoustics signals of male calling song (e.g., Chobanov et al., 2015;Kaya et al., 2012).Both characters are postulated to play a significant role in mate recognition (see Eberhard, 1985;Heller et al., 2015;Ragge & Reynolds, 1998 and references there in).Such character-based species definitions imply a direct link to a mechanism of reproductive isolation, and such characters are, therefore, a prime candidate for understanding the bush cricket's diversification (Heller et al., 2015;Heller & Hemp, 2020).
However, sometimes there are individuals of intermediate morphologies between species, which has fuelled taxonomic discrepancies (Borissov et al., 2023;Chobanov et al., 2015;Kaya et al., 2012;Ramme, 1933;Tarbinsky, 1932;Ünal, 2012;Zhantiev & Korsunovskaya, 2015).Moreover, variation in the cerci (namely, in length, curvature, thickness, apical structure and the number/orientation of distal teeth) and variation in temporal parameters of the three types of male calling song (see Chobanov et al., 2015;Kaya et al., 2012) may represent a breakdown of reproductive barriers, or the ephemerality of incipient species divergences (Rosenblum et al., 2012).Yet, the parapatric and sympatric distributions of constituent taxa, especially in the Marmara Sea Basin, suggests taxa in the PB species group may have acquired sufficient reproductive isolation to co-exist, although syntopic occurrence is rare or absent (Chobanov et al., 2015;Kaya et al., 2012;Ünal, 2005, 2010;Zhantiev & Korsunovskaya, 2015; personal observation during collecting the material used in this study).Alternatively, potentially porous species boundaries, especially those in species rich areas such as the area surrounding the Marmara Sea, suggest that species diversity may be threatened by loss via hybridization (Quilodrán et al., 2020;Todesco et al., 2016), rather than representing the pinnacle of a biodiversity hotspot.Genetic studies offer a critical framework to test for a correspondence between morphologically delimited taxonomic units and genetic divergence, to corroborate the hypothesized role of male sexual characters in reproductive isolation and diversification.Likewise, with a genetic window into the history of diversification, tests of the hypothesis that species within the Marmara Sea region are each-others closest relatives would provide support for a model of in situ diversification, whereas a model of diversity accumulation via dispersal would be supported if geographically circumscribed species are related to more geographically distant taxa.However, currently genetic studies on the genus are either limited to small species groups (e.g., Borissov et al., 2020Borissov et al., , 2021;;Boztepe et al., 2013;Kaya et al., 2015;Kociñski, 2020) or are constrained by limited number of sequenced individuals per species (e.g., Borissov et al., 2023;Ullrich et al., 2010) that constrains tests about putative species boundaries and tests about the geography of divergence.
In our study, high geographic coverage and sampling per population and per putative taxon of the PB species group provide the requisite framework to test hypotheses about the diversification of a fascinating group of bush crickets, including insights into the processes related to diversification across the Black Sea Basin.Specifically, with a phylogenetic framework estimated from two mitochondrial genes and two nuclear gene fragments, we are positioned to test a variety of hypotheses about the timing and mechanisms of diversification, as well as the evolutionary ephemerality of the constituent taxa and their species boundaries, through tests (i) for a correspondence between morphologically defined species and genetic lineages, (ii) for a coincidence between distributions of constituent taxa and their evolutionary relationships and (iii) for the relative timing of species divergence and the approximate time of the groups diversification.As one of the first genetic studies with comprehensive geographic sampling of the PB species group, the study's findings also offer insights into the evolution of Euro-Siberian biota in the Black Sea Basin.

Taxon sampling and identification
Individuals from almost all of the currently recognized species of the PB species group (Orthoptera, Tettigonidae; Phaneropterinae) (Borissov et al., 2023;Chobanov et al., 2015;Kaya et al., 2012) were collected across the total geographic ranges of the species group (Figure 1) during 2015-2019 (see Table S1 for sample site details).The only exception is for P. athos that is isolated to the Athos Peninsula, Greece.
However, uncertainty about the identification of some collected individuals from portions of their purported range reduced the taxonomic representation in our study.Specifically, individuals of P. similis, P. geoktschajkus, P. djakonovi and P. oligacanthus (see Borissov et al., 2023;Kaya et al., 2012) from the northern and southern Caucasus collected from east of the Çoruh River were all identified as P. similis.In addition, P. similis-like individuals collected from the Black Sea Basin of Anatolia west of the Çoruh River show affinities with P. similis but also exhibit morphological differences.We, therefore, treat these two groups separately and designate them as P. similis-sp1 and P. similis-sp2, respectively, in this study.Because of this uncertainty, we were able to identify 23 morphospecies (see Table S1).In addition, several out group taxa were also included in our molecular analyses: specifically, P. lodosi, P. zonatus, P. inflatus, P. glandifer and P. hamatus, which are members of different Poecilimon species groups, as well as Isophya major.
We also note there is morphological variation across the ranges of some species in areas of parapatry.Specifically, P. bischoffi individuals from the same region as P. similis sp2 show a mix of morphological traits, but these individuals are morphologically distinguished by cerci typical of P. bischoffi in this region.Likewise, parapatric individuals of P. miramae, P. bidens, P. bosphoricus and P. istanbul from the Northern Marmara Basin (but not individuals from other parts of the species respective ranges) show a mix of morphological characteristics, as do individuals of P. sureyanus, P. diversus and P. similis-sp2 from the Western Black Sea and Southern Marmara Basins, and individuals of P. boldyrevi, P. scythicus, P. tauricus and P. pliginskii from Crimea.

DNA extraction and sequencing
DNA was extracted from the muscle of the hind femur with a proteinase K digestion and following the salt/isopropanol protocol (Aljanabi & Martinez, 1997).We amplified and sequenced the mitochondrial genes cytochrome C oxidase subunit I (COI) and nicotinamide adenine dinucleotide subunit 2 (ND2), and two nuclear ribosomal internal transcribed spacers of 5.8S rDNA, hereafter referred to as ITS; for details, see File S1.Both strands were sequenced on a 23 ABI 3730XL DNA analyser by Macrogen Europe (Macrogen Inc., Amsterdam, the Netherlands).GenBank accession numbers of sequences are given in Table S2.
Contigs from the forward and reverse sequences were visualized using SEQUENCHER v.4.1.4(Gene-Codes Corp.) and aligned using MAFFT v.7.245 (Katoh & Standley, 2013) online version (http://align.bmr.kyushuu.ac.jp/mafft/online/server/) with the default setups (FFT-NS-i strategy, scoring matrix 200PAM (k = 2), gap opening penalty = 1.53).The numt probability of the COI and ND2 sequences was assessed following the criteria of Kaya and Çıplak (2018).Because the mtDNA and nuclear loci were not sequenced in each individual, datasets for each of the mitochondrial (COI and ND2) and nuclear loci (ITS) were analysed separately; concatenating across loci is not possible without a significant reduction in individuals and geographic coverage (e.g., some localities were sequenced for COI, but not ND2, or some individuals had only ITS sequences).For example, 565, 849 and 678 individuals were sequenced for COI, ND2 and ITS, respectively, versus 381 individuals with sequences of COI and ND2 (for details about sequences see Tables S1 and S3).We recognize this is not ideal; however, given our study is focused on the geography of divergence our priority is on including all individuals so analyses were then run separately for each gene.Moreover, concatenating the data with the systemic pattern of missing data would introduce an unwanted bias that is undesirable.Unique haplotypes and their frequencies were identified by DNASP v.5 (Librado & Rozas, 2009) and the nucleotide composition, the number of variable sites and indels were calculated with MEGA v.X (Kumar et al., 2016) for each locus.Saturation probabilities of single gene datasets were assessed using DAMBE v.7 (Xia, 2018).

Phylogenetic inference and divergence time estimation
Phylogenetic relationships were estimated by maximum likelihood (ML) and Bayesian methods using RAxML (Stamatakis, 2006) via CIPRES with 1000 replicates and MRBAYES v.3.2.2 (Ronquist et al., 2012).Two independent runs and four Markov chains were used for Bayesian analysis of the COI, ND2 and nuclear ITS loci.Five million generations were used in analyses of the mitochondrial datasets and 15 million generations for the nuclear ITS dataset, sampling every 1000th generation; the first 25% of trees were discarded as burn-in and a majority-rule consensus tree was generated from the remaining trees.Substitution models for the genes were obtained using PARTITIONFINDER2 with the greedy algorithm (Guindon et al., 2010;Lanfear et al., 2012Lanfear et al., , 2017)).
To estimate divergence times, branch lengths were estimated using a molecular clock in BEAST v.2.6.1 (Bouckaert et al., 2019).Each BEAST analysis was run with a relaxed lognormal clock, linked site GTR and gamma model using four discrete gamma categories, linked Yule tree model for 100 million generations sampling every 10,000 generations.The maximum clade credibility trees were built using TREEANNOTATOR implemented in BEAST, discarding the initial 25% samples as burn-in samples.The effective sample size parameters for each BEAST and MRBAYES analysis were monitored by TRACER v.1.7 (Rambaut et al., 2018), and trees were visualized using FIGTREE v.1.4.2.The BEAST analyses were calibrated using substitution rates estimated for Pholidopterini from the Tettigoniidae family (Orthoptera): specifically, 0.0187 substitutions per site per million year for COI and 0.018 for ND2 (Çıplak et al., 2022).This contrasts with the substitution rates of ITS in two subfamilies of Tettigoniidae, namely Phaneroprerinae and Tettigoniinae (Chobanov et al., 2017;Çıplak et al., 2015;Uluar & Çıplak, 2020), which also produced extremely different nodal ages when calibrating the ITS dataset (Çıplak et al., 2022).Thus, BEAST analysis of ITS dataset was calibrated by fixing the ancestral node of PB species group with the date from the COI chronogram.

Evaluation of species status
We assessed the species status of the putative taxa within the PB species group using the automated barcode gap discovery (ABGD) using the web platform https://bioinfo.mnhn.fr/abi/public/abgd/abgdweb.html.The ABGD uses pairwise distances to group sequences into potential species based on detecting gaps in the variation between purported intra-and inter-species groups (barcode thresholds) (Puillandre et al., 2012).We focus on analyses of the COI and ND2 datasets given their breadth of geographic coverage and levels of genetic variation (see the Results section) with a Jukes-Cantor (JC69) model to measure genetic distance and a P min = 0.001, P max = 0.1 with a relative gap width X = 0.1-1.5.Note this approach was used rather than a coalescent-based model to avoid sampling gaps that exist if we limited our analyses to individuals with multiple sequenced loci (e.g., see Becheler & Knowles, 2022;Chambers & Hillis, 2020).

Characterization of geographic structuring of genetic variation
Genetic structure was assessed in each gene by a discriminant analysis of principal components (DAPCs) (Jombart et al., 2010) using the 'adegenet' (Jombart, 2008;Jombart & Ahmed, 2011) and 'seqinr' (Charif & Lobry, 2007) packages in the R environment (R Core team, 2020).A Bayesian Information Criterion (BIC) was used to determine the number of genetic clusters for each DAPC analysis with the number of the optimal principal components estimated using a-score optimisation-spline interpolation ('optim.a.score').Graphical representation of the clusters and membership probabilities were determined using the 'scatter' function.We adopted this approach rather than a method that characterizes genetic structure under an evolutionary model (e.g., analyses with structure; Pritchard et al., 2000) due to the hierarchical genetic structure evident in the phylogenetic trees (see the Results section).

RESULTS
For the 23 morphospecies from 103 localities (Figure 1 and Table S1), we obtained sequences for COI, ND2 and ITS (for details about sequences, see Tables S1 and S3).Of the 2701 bp sequenced across the three genes, 1588 sites are polymorphic; DAMBE suggested no saturation probability for any of the genes (not shown here).There is some haplotype sharing across species: specifically,10 COI, 12 ND2 and 24 ITS haplotypes are shared among morphospecies; among the 23 species, only P. kocaki, P. heinrichi and P. scythicus did not share haplotypes with any other species (note the numbers of samples were limited in the last two species).Shared haplotypes tended to occur among related species, with a few exceptions (see Figures 2   and S1).

Phylogenetic relationships and timing of divergence
The topologies of the ML and BI trees are largely congruent.Likewise, there is a general congruence of regionally sampled areas between the COI and ND2 trees (Figures 2a,b and S1A,B).However, we note that because of differences in the geographic representation of individuals with COI versus ND2 sequences from different regions and species (see Table S1 for details), it is difficult to evaluate congruence at fine geographic scales, and that the non-overlapping representation of individuals between the two mitochondrial genes may obscure potential concordance at regional scales (see the Discussion section).
As such, we focus here on specific regional and species relationships, noting difference among the mitochondrial datasets, as well as the ones estimated from nuclear ITS.Note that most of the obvious phylogenetic discord across loci is between the mitochondrial versus nuclear ITS datasets (Figures 2b,c and S1B,C).Below, we considered the genetic structure based on the mitochondrial trees, because of apparent lack of resolution of ITS tree due to lower mutation rate of the nuclear gene, and consequently relative short branch lengths that characterize the estimated tree (see Figure 1c).Note also that the ABGD genetic-based delimitation analysis also only delineated two taxa based on ITS.
The monophyly of the PB species group is fully supported across all analyses (Figures 2 and S1).We note that the genus Poecilimon does not appear to be monophyletic, but limited representation of taxa in our outgroup cautions against drawing conclusions about Poecilimon monophyly; tests of monophyly of the Poecilimon genus were not an objective of our study, hence the limited representation of taxa.In contrast, none of the taxa (at least those with appreciable sampling) are monophyletic (Figures 2 and S1) nor were they delineated individually by the species delimitation analyses, although groups of taxa form well-supported clades that are delimited by genetic species delimitation.That is, there is no strict correspondence between morphological or acoustic species boundaries and geneticbased delimited boundaries.
Throughout the mitochondrial trees, individuals from geographically contiguous ranges (colour-coded shaded areas in Figure 1) often tend to be closely related to each other, including those sampled from areas separated by narrow sea passages (e.g., the connection between the Black Sea and the Marmara Sea; Figure 1).In the ITS tree, such geographic structuring is less obvious with the F I G U R E 2 Chronograms for the diversification of the P. bosphoricus species group throughout the Pleistocene; the tree was calibrated using the substitution rate of 0.0187 substitutions per site per million years for (a) COI, 0.018 substitutions per site per million years for (b) ND2 and (c) ITS calibrated by fixing the ancestral node of PB species group with the date from the first chronogram.The shaded portion demarcates the mid-Pleistocene Transition from 1.1 to 0.8 myr ago when duration of glacial cycles shifted from 42 to 100 kyrs.Individuals are colour-coded and represented by species codes, as in Figure 1; the sampled site number are given at the tips (see Table S1 for details) followed by the haplotype number per species.Shared haplotypes across taxa are shown with multiple colour-coded circles.Of the two bars along the external margin of the tree, the outer one marks hypothesized species boundaries based on the ABGD analysis and the inner one to the numbers of clade and the genetic clusters obtained from DAPC analysis (see Figure 3).Note that cluster C5 of ND2 (B) consists of two separate clades (named as C5a and C5b).

Geography of divergence within the P. bosphoricus complex and genetically delimited taxa
The most pronounced split among the putative species geographically is circumscribed species is the split between the putative species in Crimea and those on the north shore of the Black Sea (containing three species: P. tauricus, P. pliginskii and P. scythicus) along with those from the Caucasus along the eastern shore of the Black Sea (namely, P. similis-sp1) from the rests of the species along the south shore of the Black Sea and surrounding area of the Marmara Sea (i.e., clade 1; Figure 2a,b).In the ITS tree, the former group includes some additional taxa (i.e., P. roseoviridis, P. proximus and P. similis-sp2) distributed along the south shore of the Black Sea but given the lower level of structure observed across the ITS tree (Figure 2c) compared with the mitochondria trees (Figure 2a,b), this discord may reflect low resolution of the nuclear genes as noted above.This well supported split represents the oldest divergence within the PB species group (around 2 million years ago; Figure 3).Within clade 1 (C1; Figures 2a,b and S1), two wellsupported distinct clades are also recognized by the genetic species delineation analyses that correspond to taxa from the north (i.e., P. tauricus, P. pliginskii and P. scythicus) versus eastern shore of the Black Sea (i.e., P. similis-sp1 [SM1] and P. similis-sp2 [SM2]).
Despite this clear phylogenetic structure, a few individuals of SM1 group phylogenetically with individuals from the more distantly related clade 5 from the south shore of the Black Sea (C5; The discrete genetic clusters identified by DAPC and determined by the highest BIC values for COI (upper) ND2 (middle), and ITS (lower) with the different clusters shown as different colours along with the membership of each genetic cluster based on the species codes (see Figure 1), followed by the number of individuals in parentheses.Species coded by two letter abbreviations: SR: P. sureyanus; TU: P. turciae; TR: P. turcicus; CN: P. canakkale; WR: P. warchalowskae; KC: P. kocaki; BD: P. bidens; BF: P. bischoffi; BS: P. bosphoricus; CR: P. cervus; DM: P. demirsoyi; NS: P. naskrecki; IS: P. istanbul; MR: P. miramae; PX: P. proximus; RS: P. roseoviridis; SM1: P. similis-sp1; SM2: P. similis-sp2; HR: P. heinrichi; SH: P. scythicus; PL: P. pliginskii and TA: P. tauricus.In the ITS tree, clade 4 is part of a larger clade that contains a wellsupported clade 5 from the mitochondrial trees.Clade 5 in both COI and ND2 is comprised of taxa from the south shore of the Black Sea (e.g., P. similis-sp2; see Figures 2 and S1) and taxa north and south of the Marmara Sea (e.g., SR; see Figure 2) as does clade 6 (Figure 4).The delimitation of clade 5 differs between COI and ND2: specifically, two taxa are delineated for clade 5 for ND2, but clade 5 is grouped with clade 6 in COI (Figure 2).The well-supported clade 6 in COI, ND2 and ITS is comprised of P. cervus, P. demirsoyi, P. similis-sp2, P. miramae, P. naskrecki and P. bosphoricus (Figure 2) from the south shore of the Black Sea and along both the north and south shores of the Marmara Sea (Figure 4).Finally, clade 7 contains members from the south-east coast of the Black Sea (e.g., P. bischoffi) and along the south shore of the Black Sea (e.g., P. similis-sp2); this clade is not delimited in ITS.
It is noteworthy that species divergences taking place during the Mid-Pleistocene (Figure 3) largely correspond to the clades recognized by the ABGD species delimitation based on a barcode gap, whereas the more recent divergences tend to go unrecognized by the delimitation analyses (see Figures 2a,b and S1).What also stands out for all the clades, including clade 1, is that they all contain a mixture of putative taxa.That is, none of the currently recognized taxa are monophyletic.

Inferred genetic-clusters from DAPC analyses
The number of clusters determined by the highest BIC-values in DAPC analyses identified K = 7 clusters for COI and K = 6 clusters for ND2 and ITS, with BIC-values of 1271.88, 1452.421 and 1728.445,respectively.The number of genetic clusters identified in the DAPC analyses differed from the number of taxa delineated by the genetic delimitation analyses.The relationships among these genetic clusters, or clades with reference to estimated phylogenetic relationships (Figures 2 and S1) show a strong correspondence.However, note that caution should be applied when comparing the relative position of genetic clusters across loci (e.g., the relative position of C1, C2 and C6 for COI compared to ND2), especially with the seventh genetic cluster of COI dominating the variation along the PC axes (Figure 3).

DISCUSSION
The recent origin of the PB group corresponds to the Pleistocene (Figure 2a,b), and this timing of divergence coincides with geological and/or climatic events (discussed below) suggestive of plausible drivers of this radiation.Specifically, because the topography of the area, except the Marmara and intermontane Caucasian Basin plus Crimea (see below), was relatively stable during Pleistocene (Görür et al., 2000;Krijgsman et al., 2019;Meijers et al., 2018), climate driven latitudinal and altitudinal range shifts may have been important events structuring diversification of the PB, especially along the southern Black Sea Basin.The entire range of the PB group in Anatolia, plus adjoining Balkans and Caucasus, constitutes a complex glacial refugial system (Çıplak, 2008;Korkmaz et al., 2014); the northern range margin of Anatolia is considered to be the border of glacier/permafrost (app.52 0 latitude; see Hewitt, 2000).However, as we discuss below, range shift dynamics of the PB group seem to be considerably different than proposed hypotheses of a stable rear/leading edge population (Çıplak, 2004b;Hampe & Petit, 2005;Hewitt, 1996;Uluar et al., 2023).Based on the geography and timing of divergence, we develop hypotheses about the joint role of climate-induced distributional shifts and ephemeral barriers in the structuring of the PB group diversification.

Historical bridges and barriers of Anatolia and adjacent regions
Changes in land/sea configuration, along with climate-induced range shifts, have likely impacted the PB group radiation.There are two water channels that could potentially act as regional barriers that could have structured the geography of species divergence and distributions: namely, one to the south in the Marmara Region and one to the north separating the Greater Caucasus plate and its connection/ disconnection with Anatolia (Figure 1).However, the two do not appear to have acted as a barrier in similar ways.For example, the geographic separation of the east and west coasts of the Black and Marmara Seas by narrow water channels (Dardanelles and/or Bosphorus) has not acted as a barrier to past movement.All the regional clades with individuals distributed in the area (i.e., C2, C3, C4, C5 and C6), except C7, contain putative species/populations that are found on both sides of the channel (Figure 4).Indeed, the geological evidence highlights the very dynamic land/sea configuration in the Marmara Region throughout the Pleistocene, with a water corridor established during the early Pleistocene 1.8 mya (Alpar & Yaltırak, 2002;Elmas, 2003;Görür et al., 2000;Yaltırak et al., 2000;Yılmaz et al., 2010).However, this channel was closed on both the Black Sea and Aegean Sea sides until the end of the Middle Pleistocene Transition around 0.8-0.7 myr ago, though based on geologic evidence it is unclear if ephemeral water channels may have formed archipelagolike landscape (Yılmaz et al., 2010).Within the PB species group, the development of a complete water channel 0.8-0.7 corresponds to the estimated divergence of P. canakkale and P. warchalowskae from each other, and divergences within C3 and C5 (Figure 2a,b) suggesting its role as a barrier.Likewise, some individual putative species are restricted to one or the other side of the channel (e.g., see P. turciae, P. proximus [or P. similis proximus according to Ünal, 2012], P. kocaki, P. heinrichi, P. canakkale, P. naskrecki, P. roseoviridis and P. diversus; Figure 1), suggesting there has been little movement in the recent past.However, these restricted geographic distributions may also reflect limited time for expansion such that it is unclear if the narrow channel is now acting as a barrier (but, see Korkmaz et al., 2014), especially given there are also several cases of individual putative taxa occurring on both sides of the channel (e.g., see P. similis, P. sureyanus, and P. bosphoricus; Figure 1).Moreover, the relationships of local endemics with other putative taxa suggest that movement across the barrier (rather than in situ diversification on one or the other side) took place, which is consistent with geologic evidence of frequent temporary connections (Badertscher et al., 2011).Episodic channels are also supported by other phylogeographic studies that suggest the region acted as both a barrier and a distribution corridor between Europe and Anatolia (termed the "Anatolio-Balkan phylogeographic fault"; Chobanov et al., 2017; see also Kaya & Çıplak, 2017;Korkmaz et al., 2014).
In contrast to most of the putative taxa that belong to clades with mixed taxonomic affinities that occur in the Marmara Region and the south east coast of the Black Sea, the oldest divergence in the history of the PB species group diversification (see clade 1; Figure 2a,b) occurs in the northern range of the PB species group-specifically, the regions of Crimea and the western Caucasus and involves the putative taxa P. tauricus, P. pliginskii, P. scythicus and P. similis-sp1 (Figure 1).Moreover, the early Pleistocene was also the time that the putative taxonomic pair P. tauricus and P. pliginskii from Crimea diverged from P. similis-sp1 distributed across the western Caucasus (e.g., Georgia and Armenia).The formation of these relatively distinct putative taxa is consistent with the relative geologic stability and periods of connection followed by isolation of Crimea (by establishment of a water channel between Black Sea and Caspian Sea), a terrestrial island plate that has existed since the Middle Miocene (Popov et al., 2004).Specifically, Crimea had intermittent connections with the Russian plate and a notable connection existed around 1.8 myr ago between Crimea and the Lesser and Greater Caucasus followed by a loss of connections at least 0.8 myr ago (Krijgsman et al., 2019;see Figures 3 and 7 therein).This lack of connection appears to have structured the F I G U R E 4 Geographic distribution of genetically defined clusters from the DAPC analyses of COI, ND2 and ITS in Figure 4.Note all the clades (except clade 1) overlap in distribution to varying degrees, especially along the southern coast of the Black Sea and the eastern margins of the Marmara Sea).Also shown are pictures of (or in a few cases sketches of) the species-specific male secondary-sexual character-specifically, the male cercus (most of the cercal illustrations produced during this study; some modified from Kaya et al., 2012 andChobanov et al., 2015).Species coded by two letter abbreviations: SR: P. sureyanus; TU: P. turciae; TR: P. turcicus; CN: P. canakkale; WR: P. warchalowskae; KC: P. kocaki; BD: P. bidens; BF: P. bischoffi; BS: P. bosphoricus; CR: P. cervus; DM: P. demirsoyi; NS: P. naskrecki; IS: P. istanbul; MR: P. miramae; PX: P. proximus; RS: P. roseoviridis; SM1: P. similis-sp1; SM2: P. similis-sp2; HR: P. heinrichi; SH: P. scythicus; PL: P. pliginskii and TA: P. tauricus.
geography of divergence of clade 1 and has most likely contributed to its relative isolation from the rest of the PB species group taxa, as observed in some other lineages with limited dispersal ability (Neiber & Hausdorf, 2015).The restricted distribution of putative taxa from clade 1 (see Figures 2a and 4) might also reflect potential river barriers, namely the Coruh River.The correspondence of genetic divergence with this river contrasts with other river systems in the range of the PB species group that tend to not be associated with genetic differentiation (e.g., the Kura, the Yeşilırmak and the Kızılırmak).The Sakarya River separates some of the species to the east and west (e.g., P. turcicus, P. diversus, P. turciae, P. bidens, P. miramae, P. sureyanus, P. kocaki and P. istanbul to the west of the river, and P. proximus, P. naskrecki and P. demirsoyi to the east. However, this geographic separation could also be explained by the water corridor associated with the narrow channel of the Marmara Sea (discussed above).

Porous species boundaries versus incomplete lineage sorting of recently diverged species
The PB species group has a history of taxonomic instability (see Borissov et al., 2023;Chobanov et al., 2015;Kaya et al., 2012;Ramme, 1933;Tarbinsky, 1932;Ünal, 2012;Zhantiev & Korsunovskaya, 2015).This instability raises the prospect of potential inflation of taxonomic diversity and obscured species boundaries.However, this instability has been viewed in large part as a problem of varying sets of characters being used to define the taxa in the genus Poecilimon (e.g., Boztepe et al., 2013;Kociñski, 2020), but our analyses also indicate that the species are not clearly defined genetically as well (Figure 2a,b).That is not to say the PB species group lacks genetic structure.To the contrary, there are genetic clades recognized by both the species delimitation tests and the DAPC analyses (see Figures 2a,b and 3).But none of the delimited species' correspond to those based on morphologically defined.The lack of corroborative genetic results raises the question that too many species have potentially been recognized morphologically.However, the timing of divergence itself suggest some alternative potential explanations for the lack of correspondence between morphologically defined and genetically delimited taxa in which species diversity may not be significantly overestimated.
One alternative explanation relates to the recency of divergence and the observation that genetically recognized units are relatively much older than the divergences associated with individuals from a given morphological taxon.These older divergences separating the clades recognized by the DAPC and species delimitation analyses also demarcate divergences occurring at a regional geographic scale, which are dramatically larger geographically than the range of most of the so-called endemics recognized morphologically.If the rate of speciation is very high in the PB species group, it could lead to many recent divergences that may not have had enough time to reach reciprocal monophyly (Knowles & Carstens, 2007).However, this explanation alone would not explain why the morphologically defined taxa correspond to much more recent events than the genetically recognized clades (Figures 2 and S1).This would require another assumption about constraints on species divergence during earlier periods in the diversification history of the PB species group.Such constraints could relate to a dynamic period in history of low population persistence (Dynesius & Jansson, 2013) that resulted in an insufficient time to reach reproductive isolation after initial species divergences (e.g., Ortego & Knowles, 2022).Such ephemeral or incomplete speciation (Rosenblum et al., 2012) could explain the lack of older speciation events.An alternative explanation for the lack of correspondence between morphological and genetically delimited taxa is homogenizing gene flow.With species distributions overlapping (Figure 1), there is the potential for interspecific gene flow.There is also at least one instance in which gene flow likely explains the clustering of individuals in two distantly related clades (i.e., a few individuals of SM1 from the taxonomically exclusive clade 1 group phylogenetically with individuals from the more distantly related clade 5; Figure 2a,b).However, the maintenance of regional genetic structure despite geographic overlap of the clades (Figure 4), as well as the co-distribution of taxa within a region that occur in genetically distinct groups (Figures 3 and   4) argues against widespread gene flow.Likewise, high levels of gene flow would be inconsistent with the observed species-specific phenotypic differences in the male cercus-a secondary sexual trait hypothesized to be under sexual selection and plays a role in male-female reproductive interactions (Arnqvist, 1998;Arnqvist & Rowe, 2002;Eberhard, 1985;Hosken & Stockley, 2004;Knowles, 2000;Knowles & Otte, 2000;Marquez & Knowles, 2007).
These alternative explanations have very different ramifications for interpreting diversity estimates of the PB species group across Anatolia and adjoining Caucasian and Balkan parts.For example, rampant incomplete lineage sorting, which is consistent with the recent divergences of the morphologically recognized taxa, would still support the hypothesis of Anatolia as a biodiversity hotspot (Şekercio glu et al., 2011).Likewise, rapid speciation with constraints on successful speciation in the past (thus explaining the lack of species originations during the early diversification history of the PB species group [Figure 2]) would also support the hypothesized high species diversity.
However, if species boundaries remain porous enough to allow the observed levels of genetic mixes within regional clades (notwithstanding the caveats that challenge this explanation), then species divergence would be severely overestimated.
In the future, the additional resolution and the analytical tools genomic data would avail (e.g., coalescent-based model testing; Kubatko & Knowles, 2023) will help distinguish among the hypotheses.Such information is beyond the scope of the current study, but we are hoping to be able to collect it in the not-too-distant future.

CONCLUSIONS
Our data indicate that climate-induced range shift dynamics are different in the northern and southern sampled areas of Anatolia.Specifically, the genetic data suggest how channels connecting the Marmara and Black Seas, and the Black Sea and the Caspian Seas were associated with land bridges (in the former) versus a barrier (in the latter).
Differences in the transgressions and relative isolation over geologic history between these two regions appear to have resulted in different distributions of biodiversity.The genetic data also shows there is little association between putative species boundaries based on morphology and the phylogenetic relationship of the putative taxa.
The geographic overlap of genetic clusters that have remained relatively distinct argues against widespread gene flow via porous species boundaries.Instead, the recent origin of most putative taxa suggests the lack of genetic monophyly of the taxa may reflect widespread incomplete lineage sorting coupled with geologic constraints that limited species origins in the more distant past of the PB species group.However, future genomic-level data are needed to statistically distinguish between alternative explanations observed in the gene trees.

Figure
Figure 2a,b) in both the COI and ND2 trees, which is suggestive of gene flow.The taxonomic membership of the remaining six clades in COI (C2, C3, C4, C5, C6 and C7) and five clades in ND2 (C2, C3, C4, C5 and C6) are largely congruent; these clades also correspond to the genetic clusters from the respective DAPC analyses (see Figure 3 for clade membership across the genetic clusters).However, the relationships among these remaining clades differ between the genes.Unlike clade 1, there is also some haplotype sharing across closely related clades, although the haplotype sharing varies among clades.For example, clade 2 consists almost exclusively of P. canakkale and P. warchalowskae haplotypes (across all genes), which are distributed on the east and west side of the Dardanelles, the narrow passage connecting the Marmara Sea and the Aegean Sea (Figure 4), with a few individuals of P. turcicus, P. demirsoyi, P. kocaki and P. sureyanus that are much more abundant in other clades (Figures 2a,b and S1).Likewise, clade 3 is strongly supported and contains almost exclusively P. sureyanus, P. turcicus, P. turciae and P. diversus distributed along the north and south shore of the Marmara Sea.Both clade 2 and clade 3 are delineated by the species delimitation analyses as well.Taxa from both the east and west sides of the narrow channel separating the Black Sea from the Marmara Sea comprise a well-supported clade 4 (namely, P. proximus, P. istanbul, P. roseoviridis and P. bidens) that is delineated by the genetic delimitation analyses in both COI and ND2.