Phylogeny and biogeography of Arabian populations of the Persian Horned Viper Pseudocerastes persicus (Duméril, Bibron & Duméril, 1854)

The Persian Horned Viper (Pseudocerastes persicus) is distributed from northeast Iraq through the Iranian Plateau to western Pakistan with isolated populations in the Hajar Mountains of south-eastern Arabia. Like the other members of the genus Pseudocerastes, P. persicus is a sit-and-wait ambush feeder with low vagility, a characteristic that often results in high levels of population differentiation. In order to clarify the level of genetic variability, phylogenetic relationships, and biogeography of the Arabian populations of P. persicus we sequenced 597 base pairs of the mitochondrial cytochrome b of four individuals from the Hajar Mountains in south-eastern Arabia and inferred their phylogenetic relationships including 10 samples of P. persicus from Iran and Pakistan, four P. urarachnoides and one P. fieldi downloaded from GenBank. The four Arabian samples are genetically very similar in the gene fragment analysed and are phylogenetically very closely related to populations of P. persicus from coastal south Iran. Biogeographically, it appears that colonisation of the Hajar Mountains by P. persicus took place from Iran very recently, most probably during the last glaciation, when most of the Persian Gulf was above sea level and did not represent a barrier for dispersal.

The distribution area of the Persian Horned Viper (P. persicus) extends from northeast Iraq through the Iranian Plateau to western Pakistan with isolated populations in southeastern Arabia (Oman and UAE; Fathinia & Rastegar-Puyiani, 2010;Sindaco et al., 2013). Arnold and Gallagher (1977) first reported the presence of the species in Oman and provided morphological data on one male and one female from the Jebel Akhdar (Hajar Mountains). Subsequently, Gasperetti (1988), Feulner (1999, 2014, van der Kooij (2001), Cunningham (2002) and Grossmann, Kowalski, Zwanzig, and Zilger (2012) provided additional records and substantially extended the known distribution of P. persicus in Oman and the UAE, which seems to favour rugged and rocky areas above 800 metres that are seldom visited by humans (see Fig. 1 and Supplementary Fig. S1). Despite these new records, P. persicus still ranks as one of the rarest reptiles of the Hajar Mountains (Gardner, 2013) and very little is known about its systematics, distribution, ecology and conservation in this region. It is categorized in the Arabian Peninsula as regionally vulnerable by the IUCN (Cox, Mallon, Bowles, Els, & Tognelli, 2012). Habitat destruction due to development and quarrying has been identified as a threat to the species within the region at present and in the foreseeing future.
Like the other members of the genus, Pseudocerastes persicus is a sit-and-wait ambush predator with low dispersal capacity, a characteristic that often results in high levels of population differentiation, something that has recently been corroborated by the presence of substantial genetic variability at the mitochondrial DNA level in the Iranian populations (Fathinia et al., 2014). The Arabian populations of P. persicus are distributed across the Hajar Mountains in the extreme southeastern Arabia and are therefore isolated from the Iranian populations by a sea barrier. Despite the obvious interest from a systematic and biogeographic point of view, the lack of samples of Arabian P. persicus for genetic studies have prevented any phylogenetic studies that could clarify the origin and taxonomy of these isolated populations. Recent molecular studies indicate that the diversity of reptiles in Arabia and especially in the Hajar Mountains of Oman and UAE is much higher than it was previously thought, with some species showing very high levels of genetic variability even across very short distances (Papenfuss et al., 2009;Carranza & Arnold, 2012;Šmíd et al., 2013;Badiane et al., 2014;Metallinou et al., 2015;de Pous et al., 2016). Besides its relevance from a biodiversity point of view, to know the level of genetic differentiation and to have a good systematic knowledge that correctly reflects the phylogenetic relationships have important implications in venomous snakes. Venom composition can vary even between closely related species and therefore a robust taxonomic framework is highly relevant for toxinological research and the production and administration of the correct antivenom (see e.g. Ali et al., 2015).
In this work we use DNA sequences of the cytochrome b mitochondrial gene from four individuals of P. persicus from the Hajar Mountains to infer the phylogenetic relationships and biogeography of the enigmatic Arabian population of P. persicus. ). The specimen was basking in the afternoon sun while partly hidden in the crevice of a large rock (Appendix 2C); (4) On 17 October 2015 a live specimen was collected by J.E. and S.J. on the Hajar mountain range along the east coast of Sharjah, UAE (25.03°N, 56.21°E, 815 m a.s.l.), very close to a locality previously published by Feulner (2014). The habitat was nearly absent of vegetation and rocky and the animal was in an ambush position slightly cratered in the soft substrate between smaller rocks. Sampling, DNA extraction and sequencing. A total of 21 individuals were included in the phylogenetic analyses: 4 Pseudocerastes urarachnoides, 14 P. persicus and one specimen of P. fieldi. Two individuals of Eristichophis macmahoni were used as outgroups based on published evidence (Pyron et al., 2013). A list of all specimens included in the molecular study with their taxonomic identification, sample codes, voucher references, geographical distribution data and Gen-Bank accession numbers is presented in Supplementary Table 1. For the four new samples of P. persicus from Arabia, genomic DNA was extracted from ethanol-preserved tissue samples using the SpeedTools Tissue DNA Extraction kit (Biotools, Madrid, Spain). Partial sequences of the mitochondrial cytochrome b (cytb) were PCR-amplified and sequenced in both directions using primers L14910 5´-GAC-CTG-TGA-TMT-GAA-AAC-CAY-CGT-TGT-3´ and H16064 5´-CTT-TGG-TTT-ACA-AGA-ACA-ATG-CTT-TA-3' from Burbrink, Lawson, and Slowinski (2000) using 40 PCR cycles with the following conditions: 94° 40'', 46° 30'', 72° 1'. Geneious v. R6 (Biomatters Ltd.) was used for assembling and editing manually the chromatographs. The cytb coding fragments were translated into amino acids and no stop codons were observed. DNA sequences were aligned using MAFFT v.7 (Katoh & Standley, 2013) applying parameters by default (Auto strategy, Gap opening penalty: 1.53, Offset value: 0.0). Uncorrected p-distances with complete deletion were estimated for the cytb mitochondrial fragments using MEGA v.5 (Tamura et al., 2011). The level of substitution saturation was assessed by using the substitution saturation test of the program package DAMBE v. 5.3.46 (Xia, 2013). Phylogenetic analyses. Phylogenetic analyses were performed using maximum likelihood (ML) and Bayesian inference (BI) methods. The best-fitting model was inferred using jModeltest v.2.1.3 (Guindon & Gascuel, 2003;Darriba, Taboada, Doallo, & Posada, 2012) under the Akaike Information Criterion (AIC) (Akaike, 1973). ML analyses were performed in RAxML v.7.0.3 (Stamatakis, 2006). A GTR+G model of sequence evolution was used with parameters estimated independently for each partition and the reliability of the ML tree was assessed by bootstrap analysis (Felsenstein, 1985) including 1000 replications. Bayesian analyses were performed with BEAST v.1.7.5 (Drummond & Rambaut, 2007) using the same dataset used in the ML analysis but without outgroups. Analyses were run twice for 1x10 7 generations with a sampling frequency of 5000 generations. Models and prior specifications applied were as follows (otherwise by default): model of sequence evolution TrN+I; Yule process of speciation; random starting tree; fix mean rate of molecular clock model to 1. Posterior trace plots and effective sample sizes (ESS) of the runs were monitored in Tracer v1.5 (Rambaut & Drummond, 2009) to ensure convergence. The results of the individual runs were combined in LogCombiner discarding 10% of the samples and the ultrametric tree was produced with TreeAnnotator (both provided with the BEAST package). Nodes in the phylogenies were considered strongly supported if they received ML bootstrap values ≥70% and posterior probability (pp) support values ≥0.95 (Huelsenbeck & Rannala, 2004;Wilcox, Zwickl, Heath, & Hillis, 2002).

Results
The dataset used for the phylogenetic analyses consisted of an alignment of 597 base pairs (bp) of cytb, of which 117 were variable and 60 were parsimony-informative. Ingroup sequences do not present substitution saturation (data not shown). The results of the BI and ML analyses of the genus Pseudocerastes are presented in Figure 2. The two trees were very similar, differing only in the position of P. persicus KF314711 (sister taxon to the clade formed by coastal south Iran and Omani P. persicus populations in the ML tree) and in the branching order of the two coastal south Iran populations (KF314712 sister taxon to an unsupported clade formed by KF314713 and the Arabian populations in the ML tree). Pseudocerastes fieldi is sister to a clade formed by P. urarachnoides and P. persicus, although this latter clade is only supported by the BI analysis. The three Omani specimens sequenced for this study are genetically identical in the mitochondrial cytb gene and they differ by only 1 mutation from the UAE specimen (uncorrected genetic distance of 0.2%; see Supplementary Table S2). In the BI phylogenetic tree from Figure 2 the Arabian specimens branch as a sister group to P. persicus from Saraj and Bandar-e-Abbas in coastal south Iran (Figure 1). The uncorrected genetic distances for the cytb gene between the coastal south Iran and the Arabian populations range between 0.4%-1.2% (see Supplementary Table S2).
Apart from the four new records of P. persicus that correspond to the sequenced specimens included in the present work, other recent unpublished sightings of P. persicus personally communicated to J.E. include: Jebel Hafeet (one specimen found by a naturalist in late 2011, 24.06°N, 55.76°E, 683 m a.s.l., another one on 19 March 2016 found by Jacobo Reyes-Velasco and Joseph Manthey, 24.05°N, 55.79°E, 487 m a.s.l.); east coast of Hajar Mountains, Sharjah, UAE (one specimen found by a naturalist in early October 2015, 25.02°N, 56.20°E, 460 m a.s.l.).

Discussion
The three records of Pseudocerastes persicus from Oman represent an interesting range extension to the eastern Jebel Akhdar and the isolated Jebel Kawr (see Supplementary Figure S1). The species likely has a continuous range in the higher elevations from the Ru'us al-Jibal in the Musandam Peninsula to Jebel Qahwan in the extreme eastern Hajar Mountains. As mentioned by Cunningham (2002), this species could be limited to the higher elevations due to competition with Echis carinatus and E. omanensis at lower elevations, the latter being much larger and aggressive than P. persicus (Gardner, 2013;pers. obs.). As an Asian member of the genus Pseudocerastes, P. persicus probably tolerates slightly lower temperatures than the mainly Arabian and African carpet vipers of the genus Echis (Arnold, Robinson, & Carranza, 2009) and can therefore adapt better to the higher elevations of the Hajar Mountains. Nevertheless, some recent records situate P. persicus at altitudes of 460-487 m a.s.l., contradicting Egan (2007), who suggested that the species is only found above 600 m a.s.l.. Field surveys carried out by J.E. and S.J. around the Hajar mountain range along the east coast of Sharjah, UAE has confirmed that P. persicus lives at close proximity to other snake species including Echis omanensis, Telescopus dhara and Psammophis shokari, although E. omanensis was present in lower densities compared to lower altitude localities with permanent water where the species is abundant predating on Sclerophrys arabica (J.E. and S.J., pers. observ.). In view of the new data presented here, further research is needed on the possible competition between P. persicus and E. omanensis across their distribution ranges in the Hajar Mountains.
The phylogenetic relationships agree with the current taxonomy of Pseudocerastes and confirm that the Arabian populations belong to P. persicus (Uetz, 2015). Judged by the low level of genetic divergence between the coastal south Iran and Arabian populations and the extremely low genetic variability between the four Arabian specimens (one single mutation in 597 bp) we suggest that Arabia was colonised from Iran very recently and that the expansion across the Hajar Mountains occurred relatively quickly. Colonisation could have happened by transmarine dispersal across the Persian or Oman Gulfs or directly by land dispersal across the Persian Gulf during the last glaciations. Water depths in the Persian Gulf do not generally exceed 100 m and the average depth is only about 35 m, thus it has been widely recognised that most of the Gulf has been above sea level during glacial time (Lambeck, 1996). Evidence suggests that until the Strait of Hormuz opened up as a narrow waterway approximately 14,000 years ago, the Persian Gulf was free of marine influence and continued like that until the flooding phase that reached the present shorelines 6,000 years ago (Lambeck, 1996). The deep basin of up to 3,400 m existing in the Gulf of Oman prevented its drainage even during the peak of the glaciations and therefore colonization across this area could have only been by transmarine colonisation.