Phylogenetic relationships amongst the snake-eyed lizards of the genus Ablepharus Fitzinger, 1823 (Sauria, Scincidae) in the Iranian Plateau based on mtDNA sequences

We recovered molecular phylogenetic relationships amongst species of the genus Ablepharus in Iran and Iraq. Partial sequences of three mitochondrial genes (cytochrome C oxidase subunit I – COI, 12S rRNA and 16S rRNA) were analysed. In addition, phylogenetic relationships and taxonomic evaluation of Ablepharus species in Cyprus, India, Greece, Turkey and Syria were performed using partial sequences of the 16S rRNA gene. Phylogenetic trees and estimated genetic distances showed that the Ablepharus populations of Iran and Iraq clustered into three distinct clades. One is found in northwest Iran (A. bivittatus in Ardabil, East and West Azerbaijan and Hamedan Provinces). The second clade, formed by A. chernovi, is found only in Uromia. The third and most heterogeneous clade is divided into two subclades, the first includes two lineages of Ablepharus in Khorasan Razavi and Semnan Provinces (A. pannonicus) and in eastern and south-eastern Iran (A. grayanus); the second subclade is distributed in the eastern part of Iraq and west and south-western Iran (Ablepharus sp.). Our analyses indicated that splitting of A. chernovi within the genus occurred in the early Miocene [about 22.5 million years ago (Mya)]. Ablepharus bivittatus diverged 15.2 Mya, in the middle Miocene. Ablepharus pannonicus diverged in the late Miocene (8.4 Mya) and A. grayanus separated in the late Miocene (6.7 Mya). The lineages of eastern Iraq and south-western Iran (Ablepharus sp.) diverged also in the late Miocene (7.0 Mya).


Introduction
The family Scincidae encompasses more than 25% of all living genera and species of lizards (Uetz and Hošek 2020) with a nearly worldwide distribution (Vitt and Caldwell 2013). Based on morphological characters, they consist of four subfamilies, Acontinae, Feylininae, Lygosominae and Scincinae (Greer 1970), but based on molecular analysis, they include three subfamilies Acontiinae, Scincinae and Lygosominae (Pyron et al. 2013). According to the molecular phylogenic studies, the genus Ablepharus is nested within Lygosominae and is sister (Gray, 1839) which are distributed in southern Europe (the whole Balkan Peninsula and islands of the Aegean Sea), Turkey, Syria to Egypt, Azerbaijan, Armenia, Caucasus, Tajikistan, Kazakhstan, Kyrgyzstan, Uzbekistan, Turkmenistan, Afghanistan, Iran, Iraq, United Arab Emirates, Pakistan and NW India (Fühn 1969a, b;Anderson 1999;Khan 2002;Vyas 2011). Skourtanioti et al. (2016), based on two mitochondrial genes (cytochrome b and 16S rRNA) and two nuclear markers (the melanocortin 1 receptor -MC1R and the natural killer-tumour recognition -NKTR), found that A. kitaibelii represented a species complex that includes species A. kitaibelii, A. budaki, Α. chernovi and A. rueppellii. Ablepharus bivittatus, A. chernovi, A. grayanus and A. pannonicus occur in Iran -our study area (Anderson 1999;Karamiani et al. 2015).
The two-streaked snake-eyed skink, A. bivittatus was first described as Scincus bivittatus from Perimbal, Talysch Mountains and Azerbaijan. Ablepharus bivittatus is distributed in eastern Turkey, northern Iran, Armenia and Azerbaijan (Baran and Atatür 1998;Anderson 1999;Ilgaz et al. 2007). It is expected that A. bivittatus occurs in higher elevations of the Zagros Mountains and Blanford (1876) mentioned an isolated record from the Lar area in northern Fars Province. A study of phenotypical variation of the populations of Hamedan and East Azerbaijan Provinces from Iran found there was not much difference between them (Karamiani et al. 2017).
Chernov's snake-eyed skink, A. chernovi, was regarded as a subspecies of Ablepharus kitaibelii from vicinity of the settlement Tkhit, Ashtarak Region, in the middle current of the River Razdan, Armenia. Eremchenko and Szczerbak (1986) wrote a comprehensive monograph about the genus Ablepharus. In their study, A. chernovi, which had previously been accepted as a subspecies (Fühn 1969a;Baran 1977), is considered a separate species (Schmidtler 1997). A molecular phylogenetic study confirmed A. chernovi to represent a genetically-distinct species (Poulakakis et al. 2005). Ablepharus chernovi is distributed in Razdan River Valley, Armenia, northern Syria and southern and central parts of Anatolian Turkey and Uromia in West Azerbaijan Province, Iran (Baran and Atatür 1998;Ananjeva et al. 2006;Sindaco and Jeremčenko 2008;Arakelyan et al. 2011;Karamiani et al. 2015).
The minor snake-eyed skink, A. grayanus was first described as Blepharosteres grayanus from Waggur District, northeast Kutch, India. Later, it was treated as a subspecies of A. pannonicus (Fühn 1969a). Ablepharus grayanus is distributed in northern and western India, through Pakistan and Afghanistan (Khan 2002) and was recorded from the eastern and south-eastern margins of the central Iranian Plateau (Fühn 1969a;Leviton and Anderson 1970), but authors, such as Anderson (1999), Rastegar-Pouyani et al. (2008) and Šmíd et al. (2014), have not listed A. grayanus for the Iranian herpetofauna. Karamiani et al. (2015) reported a range expansion in the Minor Snake-eyed Skink, A. grayanus, from eastern and south-eastern Iran, that, based on examination of pholido-sis and morphometric measurements, were different from A. pannonicus.
The Asian Snake-eyed Skink, A. pannonicus, is distributed in a vast area from west to south in the Zagros, central Iranian Plateau and Kopet-Dagh ranges (Anderson 1999;Šmíd et al. 2014;Karamiani et al. 2018a). This species, in its extensive distribution range, based on the morphological characters, has previously been divided into several species, but later, all of them have been treated as synonyms of A. pannonicus (Anderson 1999;Bauer et al. 2003). So far, molecular studies for phylogenetic relationships have not been performed in species of the genus Ablepharus occurring in the Iranian Plateau.
In this study, we evaluate the phylogenetic relationships amongst species, based on sequences of three mitochondrial genes: COI, 12S rRNA and 16S rRNA. In addition, divergence times within the genus Ablepharus are estimated and a hypothesis of the historical biogeography on the Iranian Plateau is presented.

Sampling and DNA extraction
The tissue samples of Ablepharus specimens were collected from various populations throughout the Iranian Plateau from April to October 2010 and 2019 (about 150 days in total). Moreover, we used voucher specimens in the Zoological Museum of Tehran University (ZUTC) and Zoological collection of Hakim Sabzevari University (ZCHSU). Collected specimens were preserved in 80% ethanol and were deposited in the Razi University Zoological Museum (RUZM, Kermanshah, Iran). The outgroup taxa, Plestiodon elegans (Boulenger, 1887) and Ophiomorus persicus (Steindachner, 1867) were chosen according to Pyron et al. (2013). For DNA extraction, tail tips were used, following Kapli et al. (2013). The quality of extracted DNA was measured using 1% agarose gels stained by 0.5 µl GreenViewer 6X and visualised under ultraviolet light.

PCR amplification and sequencing
Three mitochondrial genes were selected for molecular phylogenetic analysis: 1) a partial sequence (631 bp) of the protein-encoding COI gene, 2) a partial sequence  Table S1.
(360 bp) of the non-coding 12S rRNA and 3) a partial sequence (532 bp) of the non-coding 16S rRNA. The genes were amplified by the polymerase chain reaction (PCR) procedure using primers presented in Table 1. PCR amplification for each of the three primers was performed in the condition of 1.5 mM MgCl 2 and 36 cycles with different temperatures and times, also for COI gene using semi-nested runs (for the semi-nested PCR run, one of the primers used in the first run was used again and the other primer was within the target sequence) ( Table 2). The amplified products of three mtDNA genes were sequenced by an automated sequencer ABI 3730XL (Macrogen, Seoul, South Korea).

Alignment and genetic divergence
Alignment of the individual and concatenated COI, 12S rRNA and 16S rRNA sequences was performed with Clustal W (Thompson et al. 1994) as implemented in Bioedit version 7.0.0 (Hall 1999), then checked and corrected by eye. Sequences of COI gene were translated by the programme MEGA 7 (Tamura et al. 2013) into amino acids for checking unexpected stop codons. Uncorrected genetic distance (p-distance) and corrected genetic distance (K2p; Kimura 1980) were calculated using MEGA 7 (Tamura et al. 2013) for each gene in specimens from the Iranian Plateau (Tables 3-5); the concatenated COI,  (Posada 2008) for extracting the best-fitting evolutionary model with Akaike Information Criterion (AIC). The best fit model for the combined dataset (COI, 12S rRNA and 16S rRNA) was TVM+I+G and that for 16S rRNA gene was TIM+I+G.
Maximum Likelihood (ML) phylogenetic tree was recovered by raxmlGUI 2.0 (Edler et al. 2020) and significance was estimated by 1000 bootstrap repartitions. Bayesian Inference (BI) was performed with the programme MrBayes 3.2.5 (Huelsenbeck and Ronquist 2001) under default for the number of runs and the num-ber of chains for 10 7 generations and sampled every 1000 generations with 25% of the initial samples discarded as burn-in (Condamine et al. 2015). We used all snake-eyed lizard sequences for Reconstruct Ancestral State in Phylogenies (RASP) (Yu et al. 2012). The final trees were combined as a 50% majority-rule consensus tree to calculate posterior probabilities (PPs) for each node.

Estimation of divergence time
BEAST 1.8 (Heled and Drummond 2010) was used to estimate divergence times amongst the Ablepharus populations of this study using sequences of 16S rRNA only, because the rate for the other two genes was not available. The models and priors were applied as follows: random starting tree; clock models were set as log-normal relaxed clock with unlinked status; tree priors were set as coalescent and constant size. Finally, divergence times were assessed by the mean rate of molecular evolution for the mean rate under the uncorrelated log-normal relaxed molecular clock (ucld) priors for 16S rRNA (mean: 0.00457, standard deviation: 0.0025) (Poulakakis et al. 2005;Barley et al. 2015).

Phylogenetic relationships amongst populations of Ablepharus in Iran and Iraq
Of the 360 bp examined for 12S rRNA gene, 67 (18.6%) were variable and 25 (6.9%) were parsimony informative. From the 532 bp of 16S rRNA, 67 (12.5%) were variable and 43 (8.1%) were parsimony informative, of 631 bp of COI 217 (34.4%) were variable and 59 (9.4%) were parsimony informative. In total 1,523 nucleotides were unambiguously aligned in a combined dataset, including 416 (28.3%) variable and 147 (9.6%) parsimony informative sites. The nucleotide frequencies in the combined dataset are as follows: A = 30%, C = 28%, G = 17% and T = 25%. P-distances and corrected genetic distances (K2P; Kimura 1980) for the individual clades are listed in Table 6. Reconstructed phylogenetic trees using COI, 12S rRNA and 16S rRNA, as well as the     (Fig. 2). The largest clade in the tree, clade C, is divided into subclades (C1 and C2), with subclade C1 being subdivided into lineages (C1a, C1b; Fig. 2). Clade C is very heterogeneous, contains many populations and is distributed in a vast area, from south-eastern Iraq to southwest and west of Iran, deep along the Zagros and Alborz mountains to northeast and south of the country. Lineage C1a, attributed to A. pannonicus, includes populations from the foothill regions of Khorasan Razavi and Semnan Provinces, while the lineage C1b, attributed to A. grayanus, includes populations from eastern and south-eastern Iran. Lineage C1b was also divided into two well-supported sub-lineages (b1 and b2) of haplotypes that correspond to two geographical areas; the first (b1) consisted of the protected Parvand National Park (Sabzevar) and Kerman and the second (b2) of Sistan and Baluchistan Province. Finally, subclade C2 is distributed in the eastern part of Iraq and the west and southwest of Iran (Fig. 2). The phylogenetic relationship, produced from the molecular dataset, supports the putative new species predicted from the morphological classification of our populations, with the exception of A. pannonicus in Iraq and the west and southwest of Iran.
Of 536 bp examined for 16S rRNA gene, 143 (26.7%) were variable and 63 (11.8%) parsimony informative. The estimated nucleotide frequencies were: A = 33%, C = 23%, G = 20% and T = 24%. P-distances and K2P distances for the comparison amongst clades are listed in Table 7. Comparison of Ablepharus from Cyprus, India, Greece, Turkey and Syria with those from Iran and Iraq, based on BI in the 16S rRNA gene clearly showed that A. chernovi from Iran (clade B) is close to A. chernovi from Turkey and Syria with very low genetic distance (p-distance = 0.018). The resulting tree also indicated that populations of Khorasan Razavi and Semnan Provinces (C1a) are close to the grayanus complex. The grayanus complex forms two subclades with relatively low genetic distance (p-distance = 0.035): populations from Dostmohammad (border of Iran with Afghanistan), Kerman and west Khorasan which are closely related to the Indian single sample (b1) and populations of central areas of Sistan and Baluchistan (b2). The results showed that populations from south-eastern Iraq, western and south-western Iran are sister taxa to the remaining populations (Fig. 3). Genetic divergence between C2 and each of the other populations is more than 6% for 16S rRNA; therefore, we consider it as an undescribed species (Ablepharus sp.).

Biogeographical analysis
Estimated divergence times (Fig. 4) indicated that the ancestors of Ablepharus species diverged around 22.5 Mya and subsequently resulted in the evolution of two groups. The first group includes A. budaki from Cyprus, Syria (Alawit Mountain, Allepo, Lattakia); A. kitaibelii from Greece, Turkey (Izmir) and A. chernovi which is now distributed in a vast area in Syria (Homs), Turkey to north-western Iran (Uromia). The ancestor of A. chernovi diversified from A. kitaibelii about 13 Mya and the second group subdivided into four lineages. Divergence time of the lineages A (A. bivittatus), B and C is about 15.2 Mya (Fig. 4). With regard to genetic distance analysis, the most important divergence times are within the C populations. The split of the C2 (Ablepharus sp.) populations from the C1 populations occurred at about 10.59 Mya; about 8.41 Mya C1a (A. pannonicus) from C1b (A. grayanus complex); about 6.79 Mya C1b1 (Dostmohamad and India) from C1b2 (Kerman, Parvand Park and Sistan and Baluchistan) (Fig. 4).

Molecular phylogeny of the genus Ablepharus
The results of this study show a well-resolved molecular phylogeny and identify species, based on sequence divergence of the genus Ablepharus in Iran and Iraq and their relationship with species existing in Cyprus, India, Greece, Turkey and Syria. Different lineages within the genus Ablepharus have diverged in the early Miocene and have undergone different evolutionary histories. The phylogenetic analyses (ML and BI) reveal at least five lineages of Ablepharus within Iran and Iraq with high bootstrap values and posterior probabilities.
A molecular phylogenetic study confirmed A. chernovi to represent a genetically-distinct species (Poulakakis et al. 2005). It is distributed in Armenia, Syria and Turkey (Baran and Atatür 1998;Ananjeva et al. 2006;Arakelyan et al. 2011;Sindaco et al. 2013). This species was recorded for the first time in Iran from 40 km north of Uromia, Province of West Azerbaijan, Iran by Karamiani et al. (2018b). Our analysis of three mitochondrial genes confirmed the presence of A. chernovi (clade B) in Iran as a sister taxon to all samples in clade C. This study confirmed previous records of A. grayanus in Iran and indicated a high level of variation in populations of eastern and south-eastern Iran. In addition, analysis of the 16S rRNA gene revealed a close   relationship between populations of east Iran (5537 Dostmohamad) and India (CES09-858). Our molecular results showed that A. pannonicus is a paraphyletic taxon regarding the position of A. grayanus in the phylogenetic tree (Fig. 2). Therefore, to solve the taxonomic problems of A. pannonicus in its distribution area, the population of subclade C1a from northeast should be considered as A. pannonicus, because of its relatively short geographic distance to the type population. The populations of subclade C2 should be considered as a separate evolutionary species. Ablepharus bivittatus (clade A) and A. chernovi (clade B) are monophyletic units, supporting the taxonomic status of these species.
Although there was only one sample available for A. chernovi, it is quite divergent from the other clades in all analyses (Table 6).

Biogeography of the genus Ablepharus
Fühn (1969a) believed that the representative distribution is a good reason to accept the centre of origin of Ablepharus in SE Asia, but the fossil fragment of a jaw found in lower Pleistocene layer, extends the ancient range of Ablepharus to 800 km northwards from the Altai in Kazakhstan (Fühn 1969b;Anderson 1999). The past and current model for distribution and habitat suitability for Ablepharus confirmed that they prefer different climatic conditions across the Middle East, Central Asia (Tajikistan, Turkmenistan and Uzbekistan) to India (Karamiani et al. 2018a). Based on the fossil evidence from Altai (south of Russia, close to the border with Kazakhstan; Estes 1983), we suggest that Central Asia is the origin of the genus Ablepharus. After the rise of Kopet Dagh, Caucasus and Zagros Mountains, they have subsequently been dispersed into the Iranian Plateau. The route of its dispersal is suggested to be similar with Paralaudakia caucasia (Macey et al. 1998) and in the opposite direction as Eremias velox, which dispersed from the Iranian Plateau to Central Asia before the rise of the Caucasus Mountains (Rastegar-Pouyani 2009;Rastegar-Pouyani et al. 2012). According to Brandley et al. (2011Brandley et al. ( , 2012, Plestiodon, the outgroup taxon of this study, originated from Eastern Eurasia. Plestiodon was separated from the ancestor of Ablepharus soon after collision of the Indian Plate to Eurasia and the rise of the Himalaya Mountains and the Tibetan Plateau around Oligocene (32 Mya) (Tripathy-Lang et al. 2013). In addition, the collision of the Indian Plate with Eurasia resulted in changing the geological structure of Iran, so that the ocean between the Lut Block and the Afghan Block was closed (Aghanabati 2004). During the Eocene to the early Miocene (55 to 20 Mya), at the same time as the formation of the Red Sea, with the collision of Africa with Europe, a new Alpine orogeny was introduced into Europe and the new Tethys Sea disappeared (Darvishzadeh 2003;Berra and Angiolini 2014;Karamiani 2016). Based on our molecular clock, two groups of the snake-eye lizard had diverged about 22.5 Mya and subsequently dispersed into different areas. The first group occurs in south-eastern Europe, Greece, Turkey, north-western Iran, Cyprus and Syria and the second group on the Iranian Plateau and its adjoining eastern regions. Recently, fossil remains of skinks have been described (e.g. Ablepharus) from Greece (late Miocene and early Pliocene) and from the European part of Turkey (Georgalis et al. 2019(Georgalis et al. , 2021Loréal et al. 2020; (Anderson 1999;Karamiani et al. 2017). In the late Miocene, about 10 Mya, the rise of the Zagros Mountains (Sborshchikov et al. 1981) caused the splitting of clade C into two branches, subclade 1 and subclade 2. About 8 Mya, Helmand Zone and the eastern and northeast Iranian Plateau were the habitat of ancestors of subclade 1 and about 7 Mya, progenitors of subclade 2 occurred in the west and southwest Zagros Mountains. The subclade 1 has been dispersed in east and southeast Iran, Afghanistan, Pakistan, north and west India and Central Asia (the Pamir-Altai Mountains). Subclade 2 dispersed into the west, south and southwest Zagros Mountains (Anderson 1999;Khan 2002;Karamiani 2016). The Arabian Plate separated from Africa at the beginning of Pliocene (about 5 Mya; Girdler 1984) and this phenomenon accelerated uplifting of Caucasus and Kopet Dagh mountains in the north and the Zagros Mountains in the south of the Iranian Plateau. Almost coincident with the rise of the Zagros, the populations of Khuzestan, Iran and Basra, Iraq dispersed into other regions of west and southwest Iran. The distribution area of the lineage C1a (population of north-eastern Iran) reaches the Shahrood Fault, which is the northern boundary of the eastern and south-eastern part of the Iranian Plateau (Darvishzadeh 2003). The ancestor of sub-lineage C1b separated into two groups around 5.4 Mya in Kerman-Parvand and Saravan-Taftan (Sistan and Baluchistan) areas. The drying Neotethys Sea in the Pliocene (Aghanabati 2004) introduced divergence in populations of south-eastern Iran. The global climate in the late Miocene has changed from moist to dry (Janis 1993), causing the drying of Asia (Guo et al. 2004). Climate change and changing habitats may have a fundamental role in speciation and widespread extinction (Hewitt 1996;He et al. 2010). Considering the divergence time of clade A (A. bivittatus) from its central Asian ancestor (about 5 Mya), it may have occurred after orogeny of the Caucasus, Kopet-Dagh and Alborz mountains. It has subsequently dispersed to the northwest of Iran and then dispersed to the Zagros Mountains. The warm and humid climate (3.6-2.5 Mya; Webb and Bartlein 1992) led to the fragmentation of plant and animal populations (Bennett 1990). It can be assumed that, at this time, most populations of Ablepharus had invaded or sheltered in cool habitats (e.g. migration of A. bivittatus from northwest to lower latitudes, such as Hamadan; 3.25 Mya). The climate changed to cool and dry (Webb and Bartlein 1992) in the early Pleistocene (2.5-1.8 Mya). The snow line in the cold weather periods in Hamadan has been up to 1800 m (Wright 1962), so it is likely that the Hamedan population was sheltered at low altitudes. The effect of hot and humid climate (6.5-2.5 Mya) on subclade 2 caused divergence of Zagros (Baghmalek, north Khuzestan) from Abadan (south Khuzestan), Iran