Is there low maternal genetic variation in West Asian populations of leopard?

Persian leopards Panthera pardus saxicolor have been extirpated from over 84% of their historic range and are now limited to rugged landscapes of West Asia and the Caucasus. Understanding and maintaining genetic diversity and population connectivity is important for preventing inbreeding and genetic drift, both of which can threaten population viability. All previous analyses of intraspecific genetic variation of West Asian leopards based on the NADH dehydrogenase subunit 5 gene have reported low mitogenomic diversity. In the current study, we sequenced 959 bp of the mtDNA cytochrome b gene to describe the spatial genetic structure of 22 wild Persian leopards across Iran, which hosts most of the subspecies extant range. The findings based on phylogenetic trees and median-joining network indicated that leopards from Iran formed a distinct subclade, i.e., P. p. saxicolor. The AMOVA analysis showed significant differentiation (88.55%) between the subclades of Persian leopards and other Asian leopards. The lowest levels of haplotype (0.247) and nucleotide (0.00078) diversity were estimated in Persian leopards from Iran. Mitochondrial genome sequencing revealed only two closely related haplotypes. There was no evidence for recent sudden demographic expansion scenario in Persian leopards. The low diversity in cytochrome b gene could potentially be brought about by selective pressure on mitochondria to adapt to oxidative stress and higher metabolic rates in cold environments.


Introduction
As a large-bodied, wide-ranging felid, the leopard Panthera pardus has evolved to persist across a diverse range of habitat types (Jacobson et al. 2016). Adapting to extreme heterogeneity across its vast global range has resulted in remarkable morphological variation, encouraging taxonomists to describe up to 27 putative subspecies (Ellerman and Morrison-Scott 1966;Herrington 1986), recently revised to only 8 subspecies (Kitchener et al. 2017). Leopards also show higher genetic diversity than many other big cats such tigers P. tigris (Luo et al. 2004), jaguars P. onca (Wultsch et al. 2016) and snow leopards P. unica (Janecka et al. 2017). Nonetheless, their heterozygosity is not spatially homogenous across their global range (Uphyrkina et al. 2001).
West Asian leopards have been the subject of controversies surrounding their taxonomy, phylogeography, population genetic structure, and variability (Khorozyan et al. 2006;Kitchener et al. 2017). A total of seven putative subspecies has been reported for the region (Herrington 1986), although many of them were recently identified as synonyms based on molecular analyses (Uphyrkina et al. 2001;Rozhnov et al. 2011;Farhadinia et al. 2015). Currently, two subspecies P. p. saxicolor = P. p. dathei = P. p. sindica = P. p. ciscaucasica = P. p. tulliana and P. p. nimr = P. p. jarvisi are proposed for West Asia (Kitchener et al. 2017). Several recent studies have provided understanding of taxonomy, phylogeography, spatial population structuring, and genetic variability based on ancient samples (Paijmans et al. 2018) as well as recent mitogenome sequencing in West Asian leopards, such as Iran (Farhadinia et al. 2015), the Caucasus (Rozhnov et al. 2011), andPakistan (Asad et al. 2019).
Nonetheless, all previous analyses of intraspecific genetic variation of West Asian leopards have been limited to only one mitochondrial gene, the NADH dehydrogenase subunit 5 (NADH-5), which exhibits relatively high rates of mutation in leopards (Miththapala et al. 1996;Uphyrkina et al. 2001), accounting for 44 of 50 variable sites found in mtDNA regions of the leopards (Uphyrkina et al. 2001). Surprisingly, low mitogenomic diversity was reported in all studied populations of Persian leopards, from Pakistan to Iran and the Caucasus (Rozhnov et al. 2011;Farhadinia et al. 2015;Asad et al. 2019). Mitochondrial analyses of nearly 100 wild-born leopards revealed no more than 8 haplotypes across West Asia.
Mitochondrial markers have been extensively used to reveal the stochastic processes of coalescence and lineage sorting of leopards (Uphyrkina et al. 2001(Uphyrkina et al. , 2002Rozhnov et al. 2011;Farhadinia et al. 2015;Ropiquet et al. 2015;Paijmans et al. 2018). The objective of our study was to clarify if the low mitgenomic variation in Persian leopards is only seen on NADH-5 gene or other genes, with high variability in other leopard populations, will also show similar patterns in Persian leopards. Accordingly, the present study aims to evaluate the genetic variability of leopards using another mtDNA gene which has never been tested for West Asian populations. Within the leopard mitogenome, cytochrome b (cyt b) coding genes exhibit high levels of both inter-and intrapopulation variability (Uphyrkina et al. 2001;Ropiquet et al. 2015). We use a spatially representative sample of Persian leopards in Iran, which hosts 86% of the subspecies extant range (Jacobson et al. 2016). With the demand for leopard capturing and translocation in order to resolve the threat of severe livestock raiding (Weise et al. 2015;Farhadinia et al. 2018) or to restore the historic range of the species (Rozhnov et al. 2011;Breitenmoser et al. 2014), our findings are helpful for understanding spatial genetic structure and phylogeography of leopards in West Asia.

Sample collection and DNA extraction
We obtained 22 tissue samples, from both muscle and skin, from wild-born dead leopards confiscated from poachers or live individuals captured for a satellite telemetry program, all with known geographic origin in Iran ( Fig. 1  . Samples were collected from three major regions of leopard range associated with major mountain chains, namely, Zagros, Alborz, and Kopet Dag (Fig. 1). All samples were stored in 96% ethanol until further analysis. Total DNA was extracted from leopard tissue samples using a standard phenol/chloroform method described in Sambrook et al. (1989), employing negative controls for each batch of extractions. No detectable DNA (as visualized on agarose gel) or polymerase chain reaction (PCR) product was obtained from any of the negative controls.

DNA sequencing
We used leoF: GACYAATGATATGAAAAACCATCGTTG/ Leo-R: GTTCTCCTTTTTTGGTTTACAAGAC for a fragment of 959 bp cyt b amplification (Ropiquet et al. 2015). For each reaction, 35 cycles were performed with 5 min initial denaturation at 95°C followed by denaturation at 94°C for 0.5 min, 1.5 min annealing at 50°C, 1 min extension at 72°C, followed by a 5 min final extension at 72°C. The PCR products were checked in 1.5% agarose gels. Purification was carried out using column-based purification kits (Millipore) using a vacuum for filtering. Sequences were visualized using an ABI-3730xl genetic analyzer (Applied Biosystems; http://www.appliedbiosystems.com).

Data analysis
We aligned sequences using the progressive pairwise alignment implemented in Geneious® 11.0.3 to a 959 bp segment of mtDNA cyt b. New sequences were combined with 50 GenBank entries (Supplementary Table 2) to create the data set in order to clarify the phylogenetic relationships among leopard (P. pardus) populations. Lion (P. leo) (KC701376 and KM374704) and tiger (P. tigris) (AF053053, AF053054, and DQ151551) were used as outgroups, following Ropiquet et al. (2015).
Substitution saturation was analyzed using DAMBE6 (Xia 2017). The best-fit substitution models were evaluated using the corrected Akaike Information Criterion (AICc) and the Bayesian information criterion (BIC) calculated in PartitionFinder 2.1.1 (Lanfear et al. 2017). The Bayesian inference was reconstructed using MrBayes 3.2 (Ronquist et al. 2012) in four Monte Carlo Markov chains for 50,000,000 generations and tree sampling every 100 generations. IQ-TREE 1.6 (Nguyen et al. 2015) was used for inferring the maximum likelihood tree with 10,000 replicates. The median-joining network was inferred using Network 10.0 (available at http://www.fluxusengineering.com/sharenet.htm), in order to reconstruct the relationships among haplotypes.

Results
Totally, 22 samples of Persian leopard across Iran were sequenced for the 959 bp segment of the mtDNA cyt b. A total of two haplotypes were detected within the 22 individuals from Iran, namely, IRAN4 (N = 21 individuals) and IRAN5 (N = 3 individuals), which have been deposited in GenBank with accession numbers MK028952 and MK028953. The nucleotide differences between the two haplotypes were three transitions.
Haplotype IRAN4 occurred commonly across the study area, whereas the haplotype IRAN5 was detected only in the Alborz region (Fig. 1). Additional sequences from global geographic range of the species extracted from GenBank were used in the phylogenetic analyses, where the final dataset contained 72 sequences and 22 haplotypes ( Table 3). Result of substitution saturation indicated that base substitutions did not reach saturation, so the dataset can be used in phylogenetic analyses.
The phylogenetic analyses (Bayesian and maximum likelihood trees) resulted in identical topologies included two clades: Asian and African. Each of these clades was subdivided into two subclades (Fig. 2). These clades and subclades are very well supported by bootstrap values and posterior probabilities. The leopard living in Iran formed a distinct subclade that fell into the Asian clade. The medianjoining network (Fig. 3) also demonstrated the subclades and clades detected in the phylogenetic trees. The pairwise Φ ST values indicated significant differences among the identified subclades (ranging from 0.732 to 0.945). The lowest distance (0.732) occurred between two African subclades, with Nm value of 0.183. The highest distance (0.945) estimated between the Asian subclade A (Persian subclade) and African subclade A, with Nm value of 0.029 (Table 1). Persian leopards exhibited smaller distance to other Asian leopards (India and East Asia) than African leopards. The AMOVA results revealed significant differentiations between the two identified subclades within each of the Asian (genetic variations between subclades = 88.55%) and African (genetic variations between subclades = 73.17%) clades ( Table 2).
The lowest levels of haplotype (0.247) and nucleotide (0.00078) diversity were estimated in Persian leopards (Asian subclade A), whereas the highest values of haplotype (0.942) and nucleotide (0.00795) diversity were calculated in the African subclades (Table 3). The results of the neutrality tests (Fu's F S = 2.272, Tajima's D = − 0.260, R 2 = 0.1234) were not significant for Persian leopards, where negative value was only calculated for Tajima's D test. Thus, the sudden population expansion scenario can be rejected in Iran. However, the neutrality tests for the African subclade A were significantly negative (Fu's F S = − 5.873, p < 0.01; Tajima's D = − 1.21; R 2 = 0.835, p < 0.05). Fig. 2 Phylogenetic relationships of Persian leopards (Asian subclade A) with other identified subclades of leopard (P. pardus) derived from a sequence of 959 bp mtDNA cyt b and rooted with P. leo and P. tigris. The numbers on the branches are bootstrap values and posterior probabilities in the maximum likelihood and Bayesian analysis, respectively. IRAN4 and IRAN5 denote detected haplotypes in the present study. For detailed information on sequences, see Table S1 Discussion We found low mitogenome diversity for Persian leopards in Iran without substantial population subdivision, based on the mtDNA cyt b gene. Our data supported a common maternal genealogy for leopards in West Asia and the Caucasus (Rozhnov et al. 2011;Farhadinia et al. 2015) without evidence of demographic expansion. The unstructured matrilines suggest that Iran's leopards currently form a single gene pool, although nuclear markers are needed to confirm this hypothesis.
Previous studies found low haplotype diversity in leopards living in West Asian mountains, such as Pakistan, Iran, Turkmenistan, and the Caucasus based on NADH-5 gene (Rozhnov et al. 2011;Farhadinia et al. 2015;Asad et al. 2019). Our analyses using another mitochondrial marker, i.e., cyt b confirmed the low variability in mitogenomic sequences of leopards. Similarly, snow leopards, another Panthera species roaming Asia's rugged mountains (Janecka et al. 2017), show low haplotype diversity. In contrast, other wide-ranging large carnivores in Asian highlands such as brown bear Ursus arctos and Himalayan gray wolf Canis himalayensis show remarkable haplotype diversity (Murtskhvaladze et al. 2010;Ashrafzadeh et al. 2016;Werhahn et al. 2018). Felids, especially Panthera species, have the most recent divergence time and consequently a lower molecular evolution rate compared with Ursids and Canids (Wayne et al. 1991). This has resulted in little intraspecific genetic diversity and shorter inter-specific genetic distances (Slattery et al. 2000;Kim et al. 2016).
The mitochondrial genomes of leopards are highly variable across most of their global range (Uphyrkina et al. 2001;Ropiquet et al. 2015;Anco et al. 2017). There are two potential hypotheses for low haplotype diversity in Persian Haplotypes inside dark red circles denote to IRAN4 and IRAN5 obtained from Persian leopards in this study  (Ellegren and Galtier 2016), similar to Amur tiger P. t. altaica (Russello et al. 2004), Amur leopard P. p. orientalis (Uphyrkina et al. 2002), and Asiatic cheetah Acinonyx jubatus venaticus (Charruau et al. 2011). All these have experienced a recent or ongoing demographic bottleneck with a small founding population. Although Persian leopards survived widespread eradication attempts in the Caucasus during the 1900s, leaving a few dozen (Khorozyan and Abramov 2007), no contemporary demographic bottleneck is known from Iran, with an estimated population of between 550 and 850 individuals and wide distribution by the end of the 1900s (Kiabi et al. 2002) with some areas harboring high leopard densities (Farhadinia et al. 2019). It is therefore less likely that demography can comprehensively explain the current genetic pattern in Persian leopards. Alternatively, it might be associated with the selective pressure on endothermic vertebrates to cope with the combined effects of hypoxia and cold stresses at high altitudes (Cheviron and Brumfield 2012). There, they need adaptive evolution of the mitochondrial oxidative phosphorylation system to regulate oxygen usage and energy metabolism (da Fonseca et al. 2008;Hassanin et al. 2009;Werhahn et al. 2018). Similarly for another montane-living large felid, Janecka et al. (2017) suggested that the snow leopard mitogenome may have undergone selective pressure in hypoxic environments, resulting in low mitochondrial diversity. Future researches can test this hypothesis using leopards sampled at different elevations analyzed for candidate hypoxia genes (Zhang et al. 2014;Werhahn et al. 2018).
These findings support a distinction between Asian and African leopard populations. We found additional strong support for the divergence between Persian leopards and other Asian leopards. In parallel with previous works (Uphyrkina et al. 2001;Rozhnov et al. 2011;Farhadinia et al. 2015), our study confirms a monophyletic group for the Iranian female lineage. The pattern can be expanded to leopards from other West Asian countries within the range of Persian leopard to share similar monophyletic group with Iran (Rozhnov et al. 2011;Asad et al. 2019). Lack of substantial population subdivision suggests that Iran's leopards currently form a single gene pool.
We acknowledge that our results are limited by the use of mtDNA, and consequently single locus data. However, our spatially representative samples from the majority of Persian leopard range along with previous sequencing from other range countries in West Asia will provide the basic understanding in mitochondrial variability, partitioning, and phylogeographic patterns of leopards in West Asia. Nonetheless, further studies are encouraged to apply multi-locus sampling to reveal fine-scale population structure and variability. Also, leopard samples from central Asian countries such as Afghanistan, Turkmenistan, and Kazakhstan can further enhance our understanding about the phylogeorgraphy and taxonomy of Persian leopards. Importantly, in most countries in  West Asia and the Caucasus, the majority of the remaining leopard range is along international borders, and thus in most countries, conservation of these subspecies is dependent on transboundary collaboration (Farhadinia et al. 2020). Therefore, a trans-national genetic analysis with samples representing the range of the P. p. saxicolor to understand the gene flow and spatial differentiation is desirable.
Availability of data and material All sequences have been deposited in GenBank with accession numbers MK028952 and MK028953.
Code availability All software are cited in the manuscript.

Compliance with ethical standards
Conflict of interest The authors declare that they have no conflict of interest.
Ethics approval For dead samples; no ethical approval was needed. Sampling was conducted based on permissions number 31/12630 and 93/16258 issued by the Iranian Department of Environment. We also obtained samples from live individuals captured for a satellite telemetry study in Tandroueh National Park. The Iranian Department of Environment reviewed all sampling, trapping, and handling procedures and approved permits for the work conducted (93/16270). The trapping and handling protocol was also approved by the University of Oxford's Ethical Review Committee (BMS-ERC-160614).
Consent to participate All authors agreed to participate in this study and co-authorship.
Consent for publication All authors agreed with the content and that all gave explicit consent to submit.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.