More evidence of cryptic diversity in Anatololacerta species complex Arnold, Arribas and Carranza, 2007 (Squamata: Lacertidae) and re-evaluation of its current taxonomy

. Genetic diversity is not always congruent with phenotypic heterogeneity, resulting in cryptic species complexes which cause a great struggle for scientists trying to deﬁne ‘species’ and describe relationships among taxa. Anatololacerta is a lizard genus distributed in southern and western Anatolia and some neighboring Aegean islands. Three morphospecies were recognized in Anatololacerta but a recent molecular study revealed the presence of cryptic diversity within the genus which led to the raise of a subspecies to species level. Currently the genus includes the species A. anatolica , A. danfordi , A. budaki and A. pelasgiana . Using a comprehensive sampling concerning both the number of specimens (218 specimens) and the genetic markers (3 nuclear and 3 mitochondrial), we performed phylogenetic analyses including tree reconstruction, species delimitation and divergence times estimation. The results revealed the occurrence of one more cryptic lineage which should be regarded as a separate species for which the name A. ibrahimi stat. nov . has priority. The existence of ﬁve well differentiated species with parapatric distributions in Anatololacerta is strongly supported. There is also evidence of recent and rapid radiation of the genus which probably causes phylogenetic relationships between these species to remain largely unresolved. At last, we proceeded to some nomenclatorial changes: The current name A. budaki was synonymized with A. pelasgiana because specimens of the type-locality of A. budaki are assigned genetically to A. pelasgiana . The genetic lineage including specimens currently assigned to A. budaki was named A. ﬁnikensis stat. nov . , raising the subspecies A. b. ﬁnikensis to species level.

This lack of phylogenetic resolution is also true for lower level taxa (Tamar et al., 2016;Psonis et al., 2017;Šmíd et al., 2017). Failure to reconstruct a phylogeny and recover the true historical relationships between any group of taxa could be the result of several methodological or evolutionary factors [e.g., sufficiency of reconstruction methods (Moran, Morgan and O'Connell, 2015), assumptions in the modeling of sequence evolution (Sullivan and Joyce, 2005), taxon or genome sampling (Rokas and Carroll, 2005), incomplete lineage sorting (Galtier and Daubin, 2008), hybridization (Degnan and Rosenberg, 2009), recent and/or rapid diversification (Philippe et al., 2011;Giarla and Esselstyn, 2015)].
The above study revealed that Anatololacerta is a very diversified taxon including cryptic variation. Yet the phylogenetic relationships of the four species were not fully resolved and taking into account the presence of one ambiguous clade, it is essential to further elucidate the phylogenetic history within the genus. Using a more comprehensive sampling both in terms of genomic sampling and geographic coverage ( fig. 1), we performed phylogenetic reconstructions using nuclear DNA (nuDNA) and mitochondrial DNA (mtDNA) markers. Gene tree approaches were combined with coalescentbased species tree methods, haplotype networks, species delimitation as well as divergence time estimation to clarify the speciation processes and the delimitation of operational taxonomical units within Anatololacerta, revising its current taxonomy.

Sample collection
A total of 218 Anatololacerta specimens were used in this study covering the distributional range of the genus. Sequences from several lacertid genera (Hellenolacerta, Iberolacerta, Lacerta, Parvilacerta and Phoenicolacerta) were used as outgroup and for age constraint comparisons. For the outgroup taxa, we used chimeric sequences (combination of several genetic data from multiple sources) in order to reach the maximum amount of genetic information for each species. All newly produced sequences have been deposited in GenBank. Information about all samples used in this study is presented in supplementary tables S1, S2 Figure 1. a) Map containing sampling localities of Anatololacerta specimens used in this study. Different colors represent the five main phylogenetic lineages (species) indicated in our results. A case of sympatry between clades B and C which was discovered in a sampling locality is shown in the map with a small black dot. b) Type localities of subspecies budaki and finikensis (sensu Eiselt and Schmidtler, 1986). According to current nomenclature clade C is named A. budaki. But, specimens of the type-locality of A. budaki (2.5 km N Çobanisa) are assigned genetically to A. pelasgiana. Therefore, the name A. budaki should be synonymized with A. pelasgiana. Since the name A. budaki is not available for clade C, the clade should be named Anatololacerta finikensis, which is the next available name. and supplementary text S1, while a map containing codes of each specimen and type localities of all subspecies is presented in supplementary fig. S1.

DNA extraction, amplification and sequencing
Total genomic DNA was isolated from ethanol-preserved museum specimens using a standard ammonium acetate extraction protocol from small tissue pieces (tail or leg muscle). Partial segments of three mtDNA genes [12S ribosomal RNA (12S rRNA), 16S ribosomal RNA (16S rRNA), cytochrome b (cyt b)] and three nuDNA genes [melanocortin 1 receptor (mc1r), oocyte maturation factor (c-mos), recombination activating gene 1 (rag1)] were selected for the phylogenetic analyses. Primers and conditions used in PCR amplifications and in cycle sequencing reactions are shown in supplementary table S3. Double stranded sequencing of the PCR products was performed using a Big-Dye Terminator Cycle sequencing Kit (v.3.1) on ABI 3730XL automated sequencer.

Sequence alignment and data analysis
Chromatographs were checked and sequences were corrected using the CodonCode Aligner v.8.0.2 (CodonCode Corporation). DNA sequences were aligned separately for each gene with MAFFT v.7 (Katoh and Standley, 2013) with default parameters. Ends of some sequences were cut to keep the same length in all. C-mos, cyt b, mc1r and rag1 gene fragments were translated into amino acids in order to ensure the absence of stop codons.
Sequence divergences were estimated in MEGA v.X (Kumar et al., 2018), for the main lineages of Anatololacerta revealed by the present phylogenetic analysis.
The sequence data were concatenated in three different matrices including 226 sequences. The first dataset included the three mtDNA fragments (mtDNA dataset), the second included the three nuDNA fragments (nuDNA dataset) and the third included both the mtDNA and the nuDNA fragments (combined dataset). All phylogenetic analysis were separately performed for each dataset as described below.

Phylogenetic analysis
Phylogenetic estimation was performed on the three different datasets using two methods; Maximum Likelihood (ML) and Bayesian Inference (BI).
Maximum Likelihood analysis was performed with RAxML v.8 (Stamatakis, 2014) under the Generalized timereversible + gamma (GTRGAMMA) model with 20 ML searches. The confidence of the branches of the best ML tree was further assessed based on 1000 rapid bootstrap replicates (using the GTRCAT model). Higher bootstrap values indicate stronger support for the particular branch/clade being analyzed.
Bayesian Inference analysis was performed in MrBayes v.3.2.1 (Ronquist et al., 2012), with six runs and five chains for each run for 2 × 10 7 generations sampling every 1000 generations (see table S4 for details). The evolutionary models used for the BI analysis were selected for each gene separately using jModelTest v.2.1.7 (Darriba et al., 2012) based on the Bayesian information criterion (BIC) (Schwarz, 1978), ignoring the models that included both gamma distribution and invariable sites (Yang, 2006). Several MCMC diagnostics were used to check for convergence and stationarity [the plot of the generation versus the log probability of the data (the log likelihood values), the average standard deviation of split frequencies, the average Potential Scale Reduction Factor (PSRF) and the minimum value of minimum Estimated Sample Sizes (ESS)]. The first 25% of the trees were discarded by default, as a conservative measure to avoid the possibility of including random sub-optimal 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, where probabilities 95% indicate significant support.
To infer and visualize phylogenetic relationships within Anatololacerta complex, we performed haplotype network reconstructions in the six nuclear and mitochondrial DNA loci studied as well as in the concatenated mitochondrial dataset. Networks were constructed in PopART (http:// popart.otago.ac.nz) (Leigh and Bryant, 2015) using the statistical parsimony algorithm TCS (Templeton, Crandall and Sing, 1992).

Species tree and species delimitation analysis
The multi-species coalescent-based method SVDQuartets (Chifman and Kubatko, 2014) implemented in PAUP v.4.0a165 (Swofford, 2003) for the combined dataset was used to infer the species tree. Using a coalescent model, this method infers the topology among randomly sampled quartets of predefined species (we assigned individuals to species according to the phylogenetic analysis results), and then a quartet method is used to assemble the sampled quartets into a species tree. The option of exhaustive search of quartet sampling was selected and the uncertainty in relationships was measured using non-parametric bootstrapping with 100 pseudo-replicates.
Species delimitation was performed using a single-locus coalescent-based method implemented in mPTP . The mPTP analysis method uses single-locus data, in this case the mtDNA dataset, to fit the branching events of each delimited species to a particular exponential distribution, using ML optimization. mPTP assumes that branching events within species will be more frequent than between species, with each substitution having a small probability of generating branching events. The input file for mPTP is a binary rooted phylogenetic tree and we used the one produced by RAxML and the mtDNA dataset. The presence of very similar sequences may cause the method to falsely over split, so we calculated the minimum branch length to correct for this potential type of error. Two MCMC chains of 5 × 10 7 steps each were used to assess the confidence of the ML species delimitation. Convergence of each chain was assessed by calculating the average standard deviation of delimitation support values (ASDDSV). The ASDDSV is computed by averaging the standard deviation of per-node delimitation support values across all MCMC chains (Ronquist et al., 2012;Kapli et al., 2017), with values closer to zero indicating that the independent chains are converging on the same distribution. Finally, support for the species delimitation was estimated by calculating average support values (ASV) across the two runs where values close to one indicate high support for the ML species delimitation.
For comparison reasons joint Bayesian species delimitation and species tree estimation was conducted using the program BPP v.4.1.3 (Yang and Rannala, 2010;Rannala and Yang, 2013;Flouri et al., 2018) for the nuDNA dataset only since nuclear and mitochondrial loci have very different characteristics, including different mutation rates and effective population sizes, while the latter by being inherited as a single unit, thus have a single evolutionary history. The method uses the multi-species coalescent model to compare different models of species phylogeny in a Bayesian framework, accounting for incomplete lineage sorting (ILS) due to ancestral polymorphism and gene tree-species tree conflicts. The guide-trees used were the ones produced by the gene-tree analysis but different guide-tree topologies were used as well to examine their impact in the result. The population size parameters (θ's) were assigned the inverse gamma prior IG(3, 0.004) with mean 0.004/3-1 = 0.002. The divergence time at the root of the species tree (τ 0 ) was assigned the inverse gamma prior IG(3, 0.004), with a mean of 0.002. To explore the impact of prior specifications to the posterior model probabilities we also used different values and combinations of θ's and τ 's [θ = IG(10, 0.018); τ = IG(10, 0.018)]. We ran the analysis for 25 × 10 4 generations with a burn-in of 8000. Each analysis was performed at least twice to confirm consistency between runs.

Chronophylogenetic analysis
We estimated relative divergence times using the *BEAST method (Heled and Drummond, 2010) implemented in BEAST v.2.4.7. For this analysis we used a subset of the ingroup combined dataset including five individuals from each of the main clades indicated by the phylogenetic analysis (individuals used are marked with an asterisk in supplementary table S2). As in BI analysis, nucleotide substitution models were estimated according to the BIC criterion in jModelTest (supplementary table S4). Other prior specifications applied were Birth Death model for speciation, the Constant Population model for population model and the Uncorrelated Lognormal for describing the molecular clock. Because of the lack of internal calibration points in Anatololacerta we included the outgroup individual of Parvilacerta parva (Boulenger, 1887) and used as an external calibration point the split between Anatololacerta and Parvilacerta which according to Mendes et al. (2016) occurred 10.91 Mya. The analysis ran for 1 × 10 8 generations with a sampling frequency of 1 per 5000 trees. Results were analyzed in Tracer v.1.7 (Rambaut et al., 2018) to assess effective sample sizes (ESSs) for all parameters. The -lnL was stabilized prior to 1 × 10 8 and the first 10% of the sampled generations were discarded. The final tree with divergence estimates and their 95% highest posterior densities (HPD) was computed in TreeAnnotator (included in the BEAST package).
The average of the pairwise genetic distances (p-distance) between and within the five major lineages of Anatololacerta revealed by the phylogenetic analysis and the outgroup specimen Parvilacerta are shown in table 2. As can be ascertained from table 2, the lowest divergence for all mitochondrial gene fragments is located between clades C and D which appear as sister clades in the phylogenetic analysis. The highest distances observed are those between Anatololacerta lineages and Parvilacerta, with the highest being the one between Parvilacerta and A. anatolica (clade A) reaching 20% in cyt b. Regarding the divergence of nuDNA, clades A and E seem to be the most distant.
All MCMC diagnostic metrics indicated that the iterations of BI analysis reached convergence and stationarity for all datasets. The average standard deviation of split frequencies was less than 0.01 in nuDNA/mtDNA datasets and 0.016 for the combined dataset (according to MrBayes 3.2 manual, a rough guide is that an average SD < 0.01 is very good indication for convergence, while values between 0.01 and 0.05 may be adequate depending on the purpose and therefore the runs were considered as converged), the plot of generation versus log-likelihood of the data had characteristic "white-noise" morphology after burnin. In addition, for all parameters, the Potential Scale Reduction Factor (PSRF) values were near 1.00 (range 0.999-1.000) and Estimated Sample Sizes (ESS) values were well over 100 for all parameters.
Maximum Likelihood (−lnL = 10 029.92) and BI (−lnL = 10 453.27) analysis of the concatenated mtDNA data produced identical topologies (full resolution tree in supplementary fig. S2). Both analyses identified the presence of five well supported lineages. The first [clade A: posterior probability (pp) = 1.00, and bootstrap support (bs) = 89], corresponds to the species A. anatolica. It includes all populations north of Büyük Menderes River in Turkey and the insular populations of the  Greek islands Samos and Ikaria. The second clade (clade B: pp = 1, bs = 100) corresponding to A. pelasgiana, includes the populations from Mugla Province, the mainland populations spread across the provinces of central Denizli and northwestern Antalya in Turkey and the insular populations inhabiting the Greek islands Symi and Rhodes. Clades C and D, inhabiting southcentral Turkey, appear as sister clades in the phylogenetic trees. The first one (clade C: pp = 1, bs = 99) corresponding to the taxon currently named A. budaki [or A. finikensis (Eiselt and Schmidtler, 1986) stat.nov., see Taxonomic rearrangements], includes populations distributed in the west part of Antalya Province in Turkey and the Greek islet Psomi. The second one is a new Anatololacerta genetic lineage (clade D: pp = 1, bs = 77) and includes populations from Burdur city to Göksu River in Turkey. This clade was also observed in the previous study (clade V of Bellati et al., 2015) but was not described as a separate species as its relationship with the rest of the clades was ambiguous probably due to poor sampling. Finally, the fifth clade (clade E: pp = 1, bs = 94) which corresponds to A. danfordi, includes all the populations east of Göksu River. Relationships between clades were poorly resolved except from the close relation of clade C and D (pp = 1, bs = 99) and the differentiation of clade E (pp = 1, bs = 100).
Phylogenetic analysis of the concatenated nuDNA data with ML and BI (trees not shown), produced unresolved trees with almost no support in any branch concerning Anatololacerta specimens. The only information obtained is the differentiation of A. anatolica (clade A: pp = 1, bs = 68) from the other branches.
Maximum Likelihood (−lnL = 16 815.47) and BI (−lnL = 18 005.33) analysis of the combined dataset produced identical topologies ( fig. 2). The analysis recovered the same five major monophyletic clades (A-E) as the analysis of the mitochondrial data but with a slightly different topology although with higher support values for the phylogenetic relationships among lineages. Clades C and D still appears to be the closest lineages. This sister-group relationship with clades E and B is not well resolved (pp = 0.77, bs = 32 and pp = 0.91, bs = 55 respectively). Clade A seems to be the first branch-off.
Haplotype network for Anatololacerta inferred from concatenated mtDNA is shown in fig. 3. Networks of the six individual loci are given in supplementary figs S3 and S4. Patterns of mtDNA mostly agree with the output of phylogenetic analysis. There is a case of allele sharing in 16S between clades A and D and in 12S between clades C and D. In nuDNA, very high levels of allele sharing were observed suggesting possible incomplete lineage sorting. Alleles of c-mos are shared across all species while in Regarding the outgroup lineages, taking into account the phylogenetic trees and p-distances, our study confirms the phylogenetic affinity of Anatololacerta and Parvilacerta indicated by other authors (Arnold, Arribas and Carranza, 2007;Hipsley et al., 2009;Mendes et al., 2016).

Species tree estimation and species delimitation analysis
Species tree estimation was not able to resolve relationships among the five main clades except from the sister relationship between clades C and D (bs = 85). Species delimitation analysis with BPP strongly supported the presence of five well differentiated species with posterior probability (pp) for presence of nodes equal to 1.0 in all nodes, when a 5-species assignment was used as input. Different prior combinations and guide-tree topologies had no impact on this result. mPTP analysis delimited 7 species as clade B was split into three species.

Divergence times
According to the results, the first split in Anatololacerta occurred in Early Pleistocene approximately 1.62 Mya with the separation of A. anatolica (95% HDP 1.06-2.25 Mya). This result is close to the date Bellati et al. (2015) proposed (2.29 Mya). Anatololacerta pelasgiana and A. danfordi diverged 1.52 and 1.15 Mya, respectively whereas the split between A. budaki and clade D seems to be very recent (0.56 Mya).

Discussion
In this study we used 173 newly generated and 45 previously published DNA sequences from 6 loci to validate the taxa included in the Anatololacerta species-complex. Our results reveal the presence of five well differentiated species, hence one more than the current taxonomy of the genus. This new species is distributed between the city of Burdur and the Göksu River (yellow color in fig. 1). According to Eiselt and Schmidtler (1986) this is an area where ranges of two taxa, ssp. bileki and ssp. ibrahimi, partially overlap. It is also an area where the sampling of the previous phylogenetic study of Anatololacerta (Bellati et al., 2015) was limited. The current study indicates the presence of only one distinct species in this area.
Our field sampling shows that ranges of the five species are non-overlapping. Only one case of coexistence of two species was observed between clades B and C indicating the presence of a contact zone between the two species ( fig. 1). Divergence within the species complex of Anatololacerta has been suggested to be allopatric (Bellati et al., 2015). The five lineages-species have probably been genetically isolated for a long time in multiple distinct refugia (Bellati et al., 2015). Bellati et al. (2015) suggested that geological factors such as wide sea-level changes and extensive tectonic uplifts of landmasses at the Plio-Pleistocene boundary (Kosswig, 1955;Davis, Harper and Hedge, 1971;Schmidtler, 1998) and Quaternary climatic oscillations resulting in southwards progression of ice sheets and intense aridification during phases of cooling (Avise, 2000;Hewitt, 2001Hewitt, , 2004 have affected the distribution of Anatololacerta ancestor, triggering the evolution and allopatric divergences of various lineages within Anatolia. The speciation of A. anatolica could be further explained by the presence of a physical barrier resulted from tectonic events in the area of Büyük Menderes River (Bellati et al., 2015). Populations north of the river have probably been isolated since then. A similar case could be true for the populations comprising A. danfordi the eastern-most distribution of the genus which is isolated east of the Göksu River. The excavation of the Göksu gorge, which is a deeply incised valley along the lower reaches of the Göksu River, started in the late Pliocene, some 2.5 Mya (Walsh-Kennedy et al., 2014). More research in the zones where two species meet regarding hybridization, the existence of contact zones or isolation is required in order to assess the evolutionary history and biogeography of the genus.
Intraspecific genetic structure in the sense of high support in the internal nodes within species, has been observed in all species. In A. pelasgiana for example, mPTP analysis delimited three species within this species. This does not mean that A. pelasgiana should be split into three species. Species should not be delimited using only genetic data because methods designed for this purpose are able to diagnose genetic structure, not species, and do not statistically distinguish structure associated with population isolation vs. species boundaries (Sukumaran and Knowles, 2017). In some cases (clades A, C and E) the genetic structure is strongly associated with spatial distribution (supplementary figs S1, S2). Further analysis of intraspecific diversity including population genetics, ecology and morphology would shed light into the relationships of taxa within species.
According to divergence times analysis speciation in Anatololacerta initiated in Early Pleistocene and could therefore be characterized recent. It is also rapid if we consider that the first three speciation events occurred in less than 0.5 million years. Rapid speciation has already been proposed for the entire family of Lacertidae (Fu, 2000;Pavlicev and Mayer, 2009;Mendes et al., 2016). A similar case was recently suggested for Lyciasalamandra that has a similar range of distribution and species/genetic diversity (Veith et al., 2016(Veith et al., , 2020. Recent and rapid radiation phenomena are often the reason of failure to recover a wellsupported phylogeny since little time exists for the evolved characters to sort to species. The amount of genetic data used in the present study did not provide enough information about phylogenetic relationships in species level. Obtaining data from as much of the genome as possible will increase the odds of including the few characters that provide phylogenetic resolution. Genome-wide SNP data are broadly used to address similar evolutionary problems (Giarla and Esselstyn, 2015;Leaché et al., 2015Leaché et al., , 2016Meiklejohn et al., 2016;Grummer et al., 2018;Blair et al., 2019;Dufresnes et al., 2019a,b) and can easily be recovered and analyzed using modern sequence technology (Elshire et al., 2011;Griffin, Robin and Hoffmann, 2011;Peterson et al., 2012) and sophisticated new analytical models (McCormack et al., 2013;Aberer, Kobert and Stamatakis, 2014;Kozlov, Aberer and Stamatakis, 2015). Genome-wide SNP data might be the only way, so far, to improve the resolution of the phylogeny.
Nevertheless, the existence of five welldifferentiated species is strongly supported by our results. Ideally, species should be identified using an integrative taxonomy approach involving data from multiple sources such as genetics, morphology, ecology, behavior and geography (Padial et al., 2010;Guillot et al., 2012). In order for the new lineage that appeared in our analysis not to remain unnamed and in line with the independently evolving species concept (De Queiroz, 2007) and the phylogenetic species concept, we recommend the following taxonomic rearrangements.
The easternmost populations of the genus [population group VII in figure S1 and population group "k" in Eiselt and Schmidtler (1986)], were assigned morphologically to Lacerta oertzeni ibrahimi by Eiselt and Schmidtler (1986). Phylogenetic analysis of the present study assigns these populations in Anatololacerta danfordi. More research regarding morphology of this population group is required in order to confirm our results.  Eiselt and Schmidtler, 1986. NMW 18355. Type locality: Approximately 3 km southwest of Finike, by the cave above the coastal road, Antalya Province, Turkey. Taxa included: Anatololacerta budaki finikensis (Eiselt and Schmidtler, 1986), Lacerta oertzeni ibrahimi partim Eiselt and Schmidtler, 1986. Distribution: Southwestern Anatolia in the west part of Antalya Province (Turkey) and Psomi islet (Greece). Remarks: According to current nomenclature, clade C is named A. budaki. But, specimens of the type-locality of A. budaki (2.5 km N Çobanisa) are assigned genetically to A. pelasgiana. Therefore, the name A. budaki should be synonymized with A. pelasgiana.
Since the name A. budaki is not available for clade C, the name Anatololacerta finikensis is the next available name and has priority. The name was first used in 1986 by Eiselt and Schmidtler and refers to its type locality which is near Finike.
Populations close to Antalya [population group V in supplementary fig. S1 and population group "i" in Eiselt and Schmidtler (1986)], were assigned morphologically to Lacerta oertzeni ibrahimi by Eiselt and Schmidtler (1986). Phylogenetic analysis of the present study assigns these populations in clade C = Anatololacerta finikensis. More research regarding morphology of this population group is required in order to confirm our results.

Conclusions
Five well differentiated species with distinct spatial distributions have been identified one of which is presented here for the first time. It is the second cryptic lineage revealed in Anatololacerta (the first was indicated by Bellati et al., 2015). Other cases of cryptic lineages in amphibians and reptiles of Anatolia have been reported in the past (Plötner et al., 2001;Fritz et al., 2007;Kyriazi et al., 2008;Akın et al., 2010), highlighting the importance of the mountains of the area as "hotspots" of biodiversity (Çiplak, 2003;Çiplak, 2004;Kornilios et al., 2011;Sindaco et al., 2014, Gül, Kumlutaş andIlgaz, 2015). Our findings support a Pleistocenic rapid radiation of the Anatololacerta species group. This might be the reason that phylogenetic relationships between species to remain largely unresolved. More research is needed regarding more genetic loci, morphology, behavior and ecology to confirm our results and achieve a robust delimitation and taxonomy of the species included as well as the reconstruction of their evolutionary history.