Phylogenetic relationships of the populations of Iranolacerta brandtii (de Filippi, 1863) (Squamata: Lacertidae) recently found in Eastern Anatolia, Turkey

The Persian Lizard, Iranolacerta brandtii, was until recently considered to be restricted to north-western Iran (Azerbaijan and Esfahan provinces). However, two recent studies have revealed the existence of populations in Eastern Anatolia, extending the range of this species for about 230 km westwards. The fragmented distribution of this species has been considered to be a consequence of the climatic oscillations during the Pleistocene and Holocene, which created events of alternating contact and isolation of populations in distinct glacial refugia. According to our obtained genealogy derived from three mitochondrial fragments (12S rRNA, 16S rRNA and cytb), the Turkish specimens cluster together but form an independent clade, sister to the individuals from Tabriz in Iran. The separation of these two clades is concurrent with the cladogenesis between the Esfahan and Ardabil clades, estimated to have taken place during the late Holocene.


Introduction
The Persian Lizard, Iranolacerta brandtii, was described by de Filippi (1863) from Basmenj (near Tabriz) in Iran as Lacerta brandtii. Much later, Arnold, Arribas, and Carranza (2007) erected the genus Iranolacerta to include two species: the ground dwelling I. brandtii (de Filippi, 1863) and the saxicolous I. zagrosica (Rastegar-Pouyani & Nilson, 1998). Iranolacerta brandtii is distributed in north-western Iran, namely in Azerbaijan province and discontinuously in Esfahan province (Ahmadzadeh et al., 2013;Ahmadzadeh, Kiabi, Kami, & Hojjati, 2008;Nilson, Rastegar-Pouyani, Rastegar-Pouyani, & Andrén, 2003), and it has been recently recorded in eastern Anatolia (Van province) ( Fig.1) (Avcı, Ilgaz, Bozkurt, Üzüm, & Olgun, 2015;Yıldız & İğci, 2015), extending the distribution of the species for approximately 230 km westwards. Anatolia had already been predicted in Ahmadzadeh et al. (2013) as potentially suitable for I. brandtii, using species distribution modeling, due to its environmental features, namely, higher seasonal precipitation and lower temperatures than the surrounding areas. The three main mountain ranges in Iran, specifically, Azerbaijan in the northwest, Alborz in the north and Zagros in the west, constitute well-known geographic barriers for the dispersion of lacertid lizards (Ahmadzadeh et al., 2008;Arnold et al., 2007). Several studies on lacertids have thus reported the existence of fragmented ranges in this area, as well as in the neighbouring regions (Kapli et al., 2013;Nilson et al., 2003;Pavlicev & Mayer, 2009). In I. brandtii this fragmentation seems to be recent and a result of progressive Holocene aridification, which led to the differentiation of the species within several isolated pockets (Ahmadzadeh et al., 2013). Two subspecies have been described for I. brandtii, the nominal subspecies, comprising the north-western part of the species range, and I. b. esfahanica (Nilson et al., 2003), restricted to the Zagros Mountains (Esfahan province) in the isolated extreme south of the range. A phylogeographic study on Iranolacerta (Ahmadzadeh et al., 2013) showed no phylogenetic support for I. b. esfahanica, with the Esfahan populations clustering with those from Ardabil (Azerbaijan province) assigned to the nominal subspecies. However, no formal synonymisation of the two subspecies was attempted. In addition, the morphological characters analysed by both Avcı et al. (2015) and Yıldız & İğci (2015) did not provide sufficient evidence for subspecific ascription of the Turkish specimens. The new locality records together with the ecological models (Ahmadzadeh et al., 2013) suggest a fragmented distribution and point to the existence of a separate refugium in Turkey.
In the present study, we use mitochondrial molecular markers to determine the phylogenetic relationships between the specimens found in Turkey and the remaining individuals of I. brandtii from Iran, in order to draw conclusions about the historical biogeography of the species.

Material and methods
Sampling, DNA extraction and sequencing. Genomic DNA was extracted from ethanol-preserved tails of 10 individuals of I. brandtii from several localities in Turkey, using the DNeasy Extraction Kit from Qiagen, following the manufacturer's protocol. A total of 26 individuals was used in this study, belonging to all known geographic populations of I. brandtii. These comprise the 16 specimens used in previous studies (Ahmadzadeh et al., 2013;Harris, Arnold, & Thomas, 1998;Pavlicev & Mayer, 2009). The geographic location of each individual is shown in Figure 1 and detailed information in Annex 1.
For all individuals, the Polymerase Chain Reaction (PCR) amplification and sequencing of three mtDNA gene fragments, the 12S rRNA, cytb and the 16S rRNA were performed using the primers 12Sa/12Sb, cytb1/cytb2 and 16Sar/16Sbr from Kocher et al. (1989) and Palumbi (1996), respectively. We used the same general PCR protocol for all mitochondrial markers. Each PCR had a total volume of 25μl, using the MyTaqTM Red DNA Polymerase (BIOLINE), with 5μl of Buffer, 5mM of each dNTP, 15mM MgCl2, 1U of Taq DNA Polymerase, 0.4mM of each primer, and 10-25ng of DNA. The amplification protocol for all PCR reactions of the three mtDNA fragments were as follows: 95°C for 1 min of initial denaturation followed by 35 cycles of denaturation at 95°C for 15 s, annealing at 48°C for 15 s and extension at 72°C for 10 s. Amplification ended with a 10 min final extension step at 72°C. All PCR products were visualized in 2% agarose gels stained with GelRed nucleic acid stain (BIOTIUM) and successful amplifications were sent to a private company for sequencing (Beckman Coulter Genomics, UK). We used Geneious Pro v.5.5.9 (created by Biomatters, available from www.geneious.com) to assemble the contigs and MAFFT v7.0.17 (Katoh & Standley, 2013) to align each set of sequences, using the default parameters. Generated nucleotide sequences have been submitted to GenBank with the accession numbers from X to Y (to be added after acceptance of the manuscript). Phylogenetic analyses. The best-fit partitioning scheme and models of molecular evolution for the dataset were selected with PartitionFinder v.1.1.1 (Lanfear, Calcott, Ho, & Guindon, 2012) under Figure 1. Geographic location of all specimens of Iranolacerta brandtii used in this study and representation (in green) of the known distribution range of the species, according to IUCN (Tuniyev et al., 2009).
the Bayesian Information Criterion (Luo et al., 2010) and with the following parameters set: branch lengths: linked; models: mrbayes or raxml; search: greedy. For both phylogenetic models (mrbayes or raxml) the best partitioning scheme was to join all markers and apply the GTR+G nucleotide substitution model of evolution.
Maximum Likelihood (ML) analyses were done with RAxML 7.4.2 (Stamatakis, 2006) as implemented in raxmlGUI (Silvestro & Michalak, 2012), performing 100 ML inferences. Node support was assessed with bootstrap analysis (Felsenstein, 1985) of 1000 replicates. Bayesian Inference (BI) was implemented with the program MrBayes v3.2.4 (Ronquist & Huelsenbeck, 2003). The Bayesian posterior probabilities were estimated using a Metropolis-coupled Markov chain Monte Carlo sampling approach, with both runs starting with random trees running for 10 7 generations, saving one tree every 100 generations to produce 100,000 trees. Both convergence and appropriate sampling were confirmed by examining the standard deviation of the split frequencies between the two simultaneous runs and the Potential Scale Reduction Factor (PSRF) diagnostic. The first 25,000 trees of each run were included in the burn-in period and discarded, and a majority-rule consensus tree was generated from the remaining trees. Following Pavlicev & Mayer (2009) and Pyron, Burdrink, & Wiens (2013), both Iranolacerta zagrosica and Anatololacerta danfordi were used as outgroups in all phylogenetic analyses. Calculation of genetic distances was performed using Mega v.5 (Tamura et al., 2011). Haplotype Network. In order to investigate the haplotype diversity and structure, we built a haplotype network for the concatenated mtDNA dataset. A median-joining haplotype network was constructed using the software PopArt (Bandelt, Forster, & Röhl, 1999; http://popart.otago.ac.nz) with the parameter epsilon set to 0. given above the nodes and the star represents nodes with bootstrap values ≥ 95%. Iranolacerta zagrosica and Anatololacerta danfordi were used as outgroups. Colour code is the same as in Figure 1.

Results
A total of 1,003 bp of mtDNA, corresponding to 333 bp for 12S, 390 bp for 16S and 280 bp for cytb, was analysed. With minor discrepancies, all phylogenetic analyses generated similar topologies and phylogenetic relationships. The specimens from Anatolia formed an independent clade (Figure 2), sister to the individuals from Tabriz, although this relation is not highly supported (67% bootstrap). They are both separated by 2.9% uncorrected genetic p-distance for cytb. Regarding the remaining phylogenetic relationships obtained, they are the same as in Ahmadzadeh et al. (2013), with the individuals from Esfahan (Iran) clustering with the ones from Ardabil (Iranian Azerbaijan) and forming an independent clade.
According to the obtained haplotype network, two distinct mtDNA haplotypes within the Turkish clade exist, separated by a single mutation step, one exclusive to DB24965 and the other one shared by the remaining individuals. All other clades present the same number of mitochondrial haplotypes.

Discussion
Morphological studies conducted on the population of Iranolacerta brandtii in eastern Anatolia were inconclusive regarding their subspecific assignment. Although the subspecies arrangement of northern and southern populations of I. brandtii (Nilson et al., 2003) is not supported by the mitochondrial markers used here (Ahmadzadeh et al., 2013), the individuals found in Turkey are more closely related to the populations from Tabriz ascribed in the past to the subspecies I. b. brandtii, as indicated by Yıldız & Iğcı (2015) and based on pholidosis. According to the tree topologies, the separation of the Anatolian and Tabriz clades would be contemporary with the recent cladogenesis between the populations from Ardabil and the Zagros Mountains during the late Holocene (Ahmadzadeh et al., 2013).
During the Pleistocene climatic oscillations, North Africa and the Middle East experienced a shift from dry to humid periods (Le Houérou, 1991), which caused a series of expansions and contractions of climatic zones. These events led to the establishment of patterns of genetic subdivision currently observed in several animal groups (e.g. Fonseca, Brito, Paulo, Carretero, & Harris, 2009;Gonçalves et al., 2012). The last glacial period (Würm) reached its maximum around 18,000-21,000 years ago and in the Middle East corresponded to a period of higher aridity than at present (Ehlers & Gibbard, 2004). By contrast, once the glaciations finished the climate in this region during the Holocene Climatic Optimum (HCO; 9,000-5,000 years ago) began to become more humid than is currently the case (Kaufman et al., 2004). Therefore, one would expect that at first, during the glaciations, the populations of I. brandtii would have shifted their range to lower altitudes, but due to aridification and the establishment of unsuitable habitats between them they would have remained isolated. During the HCO, the increase in humidity would have favoured expansion and contact between populations, while the increase in temperature and aridification during the late Holocene would have caused them to retract to higher altitudes. As a result of these climatic oscillations, I. brandtii presents a fragmented geographic distribution with four mtDNA genetic lineages which have diverged over several ice ages within different glacial refugia (Ahmadzadeh et al., 2013). In fact, most of the Lacertini taxa diversified very rapidly and often have fragmented distributions, suggesting a reduction of the total range and relictual distribution for some of them (e.g. Ferchaud et al., 2015).
Such climatic shifts during the Pleistocene and the Holocene have also affected eastern Anatolia in the Van Lake area (Wick, Lemcke, & Sturm, 2003). The finding of a distinct clade of I. brandtii in this region is therefore not completely surprising and suggests that this area might have served as a glacial refugium for the Persian Lizard. In fact, recent studies on the genetic variation of several reptiles and amphibians acknowledge the potential role of the Anatolian region as a major refugium and a source of re-expansion for these taxonomic groups during the Pliocene and Pleistocene (e.g. Bellati, Carranza, Garcia-Porta, Fasola, & Sindaco, 2015;Kapli et al., 2013;Kornilios et al., 2011;Kornilios et al., 2012). Because of its geographic position and geological history through the Tertiary and Quaternary, Anatolia played an important role in the past, acting as both bridge and barrier for species dispersal between Asia, Europe and the Afrotropical region via the Middle East, by providing a natural pathway or by acting as a vicariant agent (Tchernov, 1992). through the European Regional Development Fund (ERDF). Yamanyurt specimens were collected thanks to financial assistance by TÜBİTAK (Project No. 113Z752). Gürpınar, Çaldıran and Özalp specimens were collected by MZY and NI during the National Biodiversity Inventory Project financially supported by The Republic of Turkey Ministry of Forestry and Water Affairs, General Directorate of Nature Conservation and National Parks.

Disclosure Statement
No potential conflict of interest was reported by the authors.

Supplementary Material
The Supplementary Annex 1 is available via the "Supplementary" tab on the articles online page (http://dx.doi.org/10.1080/09397140.2015.1101925).