The phylogeny, phylogeography, and diversification history of the westernmost Asian cobra (Serpentes: Elapidae: Naja oxiana) in the Trans‐Caspian region

Abstract We conducted a comprehensive analysis of the phylogenetic, phylogeographic, and demographic relationships of Caspian cobra (Naja oxiana; Eichwald, 1831) populations based on a concatenated dataset of two mtDNA genes (cyt b and ND4) across the species' range in Iran, Afghanistan, and Turkmenistan, along with other members of Asian cobras (i.e., subgenus Naja Laurenti, 1768). Our results robustly supported that the Asiatic Naja are monophyletic, as previously suggested by other studies. Furthermore, N. kaouthia and N. sagittifera were recovered as sister taxa to each other, and in turn sister clades to N. oxiana. Our results also highlighted the existence of a single major evolutionary lineage for populations of N. oxiana in the Trans‐Caspian region, suggesting a rapid expansion of this cobra from eastern to western Asia, coupled with a rapid range expansion from east of Iran toward the northeast. However, across the Iranian range of N. oxiana, subdivision of populations was not supported, and thus, a single evolutionary significant unit is proposed for inclusion in future conservation plans in this region.

N. oxiana is subdivided into eastern and western populations, split by the Hindu Kush Mountains and desert zones of southern Afghanistan, southeastern Iran, and southwestern Pakistan (Wüster, 1990). However, a new population in its eastern distribution has been recently recorded in Himachal Pradesh, India, at an altitude of approximately 2,100 m (Santra et al., 2019). Based on morphological characters, Wüster (1990) suggested that the western population is relatively homogeneous but morphologically different from the eastern population. Additionally, western N. oxiana populations appear to lack cuneate scales, whereas eastern populations, similar to most cobras of Asiatic Naja, have a cuneate scale on each side.
The distribution, demography, biology, ecology, and conservation status of the Caspian cobra remain largely undescribed.
Even though some molecular studies have investigated Asian and African cobras (Kazandjian et al., 2020;Lin et al., 2008Lin et al., , 2012Lin et al., , 2014Ratnarathorn et al., 2019;Santra et al., 2019;Wallach et al., 2009;Wüster et al., 2007), the phylogenetic and phylogeographic status, as well as the evolutionary history and population structure of the Caspian cobra in the Trans-Caspian region, are still poorly known.
The only study conducted on the genetic structure and phylogeny of the Caspian cobra in Iran, using 589 base pairs (bp) of the mitochondrial D-loop region, revealed low genetic diversity and unstructured populations of the Iranian Caspian cobra (Shoorabi et al., 2017). Yet, it is unclear when or how the species dispersed into northeastern Iran, expanding its range from the dry and semidry mountainous habitats bordering Iran, Afghanistan, and Turkmenistan to the hot F I G U R E 1 The Caspian cobra (Naja oxiana), aka the Central Asian cobra. This member of the family Elapidae is a nonspitter. Credit: Ali Khani F I G U R E 2 Top: map showing the distribution of the 11 cobra species of Asiatic Naja, ranging from Trans-Caspia to southern and southeastern Asia. Geographic ranges were obtained from the IUCN Red List (IUCN, 2020) and Wüster (1996) and updated for some species based on new occurrence records. Blue circles represent N. oxiana occurrence points that lie outside the species' known range. White dotted line represents N. sumatrana's range boundary. Bottom: map of sampling localities of the 38 specimens from Iran, Turkmenistan, and Afghanistan. Colors correspond to sampling localities of N. oxiana, namely Golestan (purple), Semnan (orange), Northern Khorasan (light blue), Central Khorasan (dark blue), and Southern Khorasan (yellow) provinces of Iran, as well as Turkmenistan (red) and Afghanistan (green). Black triangles once again represent N. oxiana's occurrence points that lie outside the species' known range and humid plains of Golestan Province, southeast of the Caspian Sea. This species is known as rare and vulnerable in Iran (Darvish & Rastegar-Pouyani, 2012) and is listed in Appendix II of CITES; however, it has been designated as "data deficient" (DD) under the criteria of the International Union for Conservation of Nature (IUCN).
For nearly a century, this cobra has been intensively harvested for venom extraction by the largest Iranian venom manufacturing center, the Razi Institute for Serum and Vaccine Research (Darvish & Rastegar-Pouyani, 2012), which is perceived to be the main cause for the dramatic decline in its wild populations.
Here, we aimed to ascertain the first insight into the phylogeny and phylogeography of the Caspian cobra across its distribution range in the Trans-Caspian region, using the Cytochrome b (cyt b) and NADH dehydrogenase subunit 4 (ND4) markers. We sought to (a) demonstrate the phylogenetic relationships of the Caspian cobra and other members of Asian cobras (subgenus Naja), (b) uncover the phylogenetic processes, including the timing and mechanisms of the species' colonization from southwest Asia into Iran and subsequently to the plains of the Caspian Sea, and (c) delineate the evolutionary lineages of the Caspian cobra across its geographical range in the Trans-Caspian region and define the species' ESUs in Iran.   Taylor, 1922;three samples of N. samarensis Peters, 1861;two samples of N. sagittifera WALL, 1913;three samples of N. sputatrix F. Boie, 1827;nine samples of N. kaouthia LESSON, 1831; and three samples of N. mandalayensis Slowinski & Wüster, 2000) from GenBank (Appendix A).

| DNA isolation, amplification, and fragment analysis
For genomic DNA extraction, we used the standard phenol/chloroform protocol. Two fragments of the mtDNA containing 1,060 bp of the cyt b and 699 bp of the ND4 loci were amplified using four pairs of modified primers (ND4/Leu (Arevalo et al., 1994) and L14910/ H16064 (Burbrink et al., 2000; Table 1). For PCR mix preparation, a final volume of 25 µl consisting of 1 µl of primer, 12.5 µl of Taq Mix (2×), 5 ng of template DNA, and 5.5 µl of double-distilled H 2 O was used. Typical amplification conditions involved a 5-min initial denaturation step at 95°C and 35 cycles at 95°C for 30 s. Next, primer annealing was done at 58°C (cyt b) and 52°C (ND4) for 30 s, and 72°C for 1 min, followed by a 5-min primer extension step at 72°C. PCR products were then analyzed and verified by agarose electrophoresis. Sequencing of the purified products was performed in an ABI PRISM 3730xl automatic sequencer (Korea Genomics, Bioneer).

| Sequence analysis
The SeqScape version 2.6 software (Applied Biosystems) was used for editing the sequences, and ClustalW via MEGA version 5 (Tamura et al., 2011) was applied for generating a multiple sequence alignment. Sequences of protein coding genes were translated to search for possible stop codons due to pseudogene production. The substitution saturation test (Xia et al., 2003) was conducted using DAMBE version 6.0.4 (Xia, 2013). Genetic diversity (nucleotide and haplotype diversities), polymorphic sites, and haplotype numbers were measured using DnaSP version 5.0 (Librado & Rozas, 2009).
Genetic distances and composition of nucleotides were examined based on the uncorrected pairwise genetic distances with 1,000 bootstraps using MEGA version 5 (Tamura et al., 2011). Furthermore, we computed all pairwise distances, between and among groups, TA B L E 1 Modified sequences of the primers used for PCR and/or sequencing

Primers
Sequences Sources

| Phylogenetic analyses
We constructed a concatenated dataset of two fragments of mtDNA We used IQ-Tree version 1.6.8 (Nguyen et al., 2014) to build an ML tree using 1,000 nonparametric bootstrap replicates and estimate support of the tree topology (Hoang et al., 2018). Furthermore, a Bayesian phylogenetic reconstruction was performed using the selected scheme in MrBayes version 3.2.4 (Ronquist & Huelsenbeck, 2003). For Bayesian inference, four Markov chains Monte Carlo runs (one cold and three heated chains (MC 3 )) were simultaneously used for 40 million generations with two replicates, sampling the trees every 1000th generation and removing the first 25% burn-in trees. By combining the remaining trees, we obtained a 50% majority rule consensus tree.
Tracer version 1.7 was implemented to visualize the results and to measure stationarity and convergence of the chains (Rambaut et al., 2018). Bayesian posterior probabilities were computed to check for support of the Bayesian tree branches.

| Analyses of population genetic structure
To detect genetic clusters in the concatenated mtDNA dataset, we applied the Bayesian inference of genetic structures of populations in BAPS version 6.0 (Corander et al., 2013). A range of 1-20 was considered for values of the number of clusters (K). Moreover, two haplotype networks were constructed to reveal haplotype relationships among the 11 Asian cobras. To reconstruct the phylogenetic relationships among Asian cobras, SplitsTree version 4.6 (Huson & Bryant, 2006) was used and a neighbor-net network was generated using the uncorrected patristic distances with 1,000 bootstrap replicates, and to plot haplotype relationships within N. oxiana populations, TCS algorithm was used in PopArt (Leigh & Bryant, 2015).

| Molecular dating estimates
To determine divergence times, BEAST version 1.8.0 (Drummond & Rambaut, 2007) was implemented using three fossil calibration points recommended by Head et al. (2016) including (a) Viperinae: the divergence between Crotalinae and Viperinae dating back to at least 20 Mya (Szyndlar & Rage, 1990;Wüster et al., 2008). For this node, we used a lognormal prior with 20 Mya as zero offset, a standard de- fossil (Hoffstetter, 1939) which marks the divergence between Naja and Haemachatus dating back at 17 Mya modeled using lognormal prior with a 95% confidence interval of 17.52-31.08.
In search of the fittest partitioning schemes for our data and models of evolution for molecular dating analyses, PartitionFinder version 1.1.1 (Lanfear et al., 2012) was utilized. Also, a birth-death process was applied as it is more fitting for a multispecies sequence dataset (Drummond & Rambaut, 2007). Fitness of the molecular clocks (strict, exponential relaxed, and lognormal relaxed) was assessed based on the Bayes factor support value (Brandley et al., 2005)

| Neutrality tests and demographic analyses
The Extended Bayesian Skyline Plot (EBSP; Ho & Shapiro, 2011) was created in Beast version 2.4.7, allowing us to estimate the demographic history of N. oxaina using the concatenated mtDNA dataset.
Strict molecular clocks were used with a mutation rate of 0.0065/ site/myr (Lin et al., 2014). The MCMC chain was set to a total of 500 million steps, and the Markov chain was sampled every 10,000 steps. To check for convergence of MCMC runs and examine effective sample sizes (>200), Tracer version 1.5 was utilized. Bayesian skyline plots were created using the EBSP R function (Heled & Drummond, 2008 (Fu, 1997) and Tajima's D (Tajima, 1989) were estimated to test equilibrium of the population in Arlequin version 3.1 (Excoffier & Lischer, 2010).
Negative statistics signify an excess of low frequency alleles, which suggests size expansion and/or purifying selection of a population (Tajima, 1989).

| Phylogenetic reconstructions
Based on our phylogenetic reconstruction analyses, our concat- Ml) (Figure 3) is new.

| Haplotype networks
We detected 23 haplotypes for N. oxiana using 38 novel sequences.

| Divergence times
Analyses of divergence placed the basal cladogenesis among the

| Phylogenetic and phylogeographic patterns of the Asiatic Naja
The evolutionary relationships among lineages of the Asiatic Naja have not yet been clearly refined and addressing the issue of the evolution of spitting behavior in African and Asiatic cobras seems to be the core subject for generating a comprehensive and robust phylogenetic resolution for cobras (Wüster et al., 2007 (Ineich, 1995;Minton, 1986). However, Wüster

et al. (2007) argued against this hypothesis and proposed that Asiatic
Naja originated from a single invasion of Asia from Africa. Barbour (1922) Figure 2 (bottom) F I G U R E 6 Bayesian spatial clustering for groups of individuals of N. oxiana performed in BAPS. The mixture analysis was set for K = 1 (one cluster) Our reconstruction of the phylogenetic relationships among Asiatic cobras at interspecific levels recovered a generally well-supported phylogenetic structure for the Asiatic Naja. Our results continue to strongly confirm the monophyly of the Asiatic Naja group, in congruence with previous studies based on molecular and morphological evidence (Slowinski & Keogh, 2000;Szyndlar & Rage, 1990;Wüster, 1990;Wüster et al., 2007). Moreover, our phylogenetic reconstructions confirmed that N. naja formed a basal lineage relative to other Asiatic Naja species, a finding similar to previous studies (Ashraf et al., 2019;Wallach et al., 2009;Wüster et al., 2007). Excluding N. naja, which forms the basal clade for our phylogenetic tree, we could distinguish two broad groups of taxa: (a) southeastern Asian cobras including N. siamensis, N. sumatrana, N. philippinensis, N. samarensis, N. sputatrix, and N. mandalayensis, and (b) western, central, southern, and eastern Asian cobras including N. oxiana, N. kaouthia, N. sagittifera, and N. atra

F I G U R E 7
Chronogram resulting from dating analyses using the concatenated mtDNA dataset (cyt b + ND4) with 116 sequences (40 sequences from Naja oxiana, 47 sequences from the other 10 Asiatic cobras, and 29 sequences from 9 genera as out-groups), generated by BEAST version 1.8.0 (Drummond & Rambaut, 2007). Branch numbers display times of divergence. Colors correspond to lineage colors in Figure

| Phylogenetic and phylogeographic patterns of N. oxiana
In this study, we attempted to present a first view of the phylog- According to the EBSP results, our demographic history reconstruction exhibited that N. oxiana's western ranges expanded between 1,000 and 2,000 years ago during the last age (Meghalayan) of Holocene series (Walker et al., 2019). Furthermore, mismatch distribution provided strong evidence for N. oxiana's unimodal distribution, which is associated with a panmictic population undergone sudden demographic expansion in its evolutionary history (Rogers & Harpending, 1992;Slatkin & Hudson, 1991). In addition, the starshaped structure of the haplotype network evidently indicates a sudden expansion (Kerdelhué et al., 2009;Slatkin & Hudson, 1991). Palearctic species (Avise, 2000;Hewitt, 2001). Studies on snakes living at high altitudes highlight that mountain vipers have been modified during the glacial climates of Pleistocene and have diversified as a result of the progression and recession of glaciers (Behrooz et al., 2018;Ding et al., 2011). Our findings corroborated that the western population of N. oxiana has a genetically uniform structure belonging to a single genetically homogeneous population, while also acknowledging the small number of samples from Afghanistan and Turkmenistan in our study.
We posit that N. oxiana was not influenced by cold Pleistocene periods and therefore not restricted to glacial refugia, in that it typically inhabits arid and semiarid regions, or rocky, shrub or scrub-covered foothills and tends to avoid high-altitude mountain habitats which were subject to past glacial episodes in Iran (Moghimi, 2008(Moghimi, , 2010

| Conservation units and management propositions
For long, conservation managers have debated how to delineate the minimal units appropriate for conservation management purposes (Amato, 1991;Fraser & Bernatchez, 2001;O'Brien & Mayr, 1991). To fulfill this aim, the term evolutionarily significant unit (ESU) emerged to define special groups of taxa below the species level that warrant specific conservation due to their evolutionary originality (Moritz, 1994;Ryder, 1986). Yet, there exists no consensus on the definition of ESU, and here, we favored the one proposed by Fraser and Bernatchez (2001) as it is focused on isolated populations. They advocated that, for adaptive evolutionary conservation, ESUs are best represented by phylogenetic lineages that exchange gene flow within a species. Therefore, each isolated population of N. oxiana qualifies as a fitting ESU for conservation.
In addition to being classified as Data Deficient (DD) by the IUCN, the conservation status of N. oxiana within the complex group of Naja has not yet been assessed and its national conservation has been largely neglected.
Here, we suggest that all N. oxiana populations in northeastern Iran, Afghanistan, and Turkmenistan could be treated as a single ESU.

ACK N OWLED G M ENTS
We would like to thank Atefeh Asadi and Ali Norouzi for their contri-

CO N FLI C T O F I NTE R E S T
The authors declare to have no competing interests.

DATA AVA I L A B I L I T Y S TAT E M E N T
Dataset has been deposited in DRYAD (https://doi.org/10.5061/ dryad.r7sqv 9s97) and submitted to NCBI (MW172773-MW172810 for cyt b and MW145451-MW145488 for ND4).