Diversification and subspecies patterning of the goitered gazelle (Gazella subgutturosa) in Iran

Abstract Goitered gazelles, Gazella subgutturosa, exist in arid and semiarid regions of Asia from the Middle to the Far East. Although large populations were present over a vast area until recently, a decline of the population as a result of hunting, poaching, and habitat loss led to the IUCN classification of G. subgutturosa as “vulnerable." We examined genetic diversity, structure, and phylogeny of G. subgutturosa using mitochondrial cytochrome b sequences from 18 geographically distant populations in Iran. The median‐joining network of cyt b haplotypes indicated that three clades of goitered gazelles can be distinguished: a Middle Eastern clade west of the Zagros Mountains (and connected to populations in Turkey and Iraq), a Central Iranian clade (with connection to Azerbaijan), and an Asiatic clade in northeastern Iran (with connection to Turkmenistan, Uzbekistan, and other Asian countries as far as northeastern China and Mongolia). Based on our results, we argue that Iran is the center of diversification of goitered gazelles, due to the presence of large mountain ranges and deserts that lead to the separation of populations. In accordance with previous morphological studies, we identified the Asiatic clade as the subspecies G. s. yarkandensis, and the other two clades as the nominate form G. s. subgutturosa. The new genetic information for goitered gazelles in Iran provides the basis for future national conservation programs of this species.


| Taxonomy of goitered gazelles
In recent taxonomic literature, following the revision by Groves and Grubb (2011), some of the previously recognized subspecies of goitered gazelle were elevated to full species status based on morphological data: the Yarkand gazelle (G. yarkandensis) in China and Mongolia (including the previously recognized subspecies G. s. hillieriana, G. s. sairensis, G. s. reginae, and G. s. mongolica), the Turkmen goitered gazelle (G. gracilicornis) in Turkmenistan, Kazakhstan, Tajikistan, and Uzbekistan, and the nominate subspecies, the Persian gazelle (G. subgutturosa) in Iran (including a previously recognized subspecies, the Seistan gazelle, G. s. seistanica, from eastern Iran), Iraq, Azerbaijan, and Turkey ( Figure 1). Recent molecular studies show that differences between these species are only detectable when using fast-evolving mitochondrial markers, preferably control region sequences (e.g., Abduriyim, Zibibulla, Eli, Ismayil, & Halik, 2018b), so a distinction on species level might not be appropriate, and the groups are referred to as subspecies throughout this manuscript. However, molecular data were able to support subspecific patterning within G. subgutturosa. Sorokin, Soldatova, Lukarevskiy, and Kholodova (2011) found a marked distinctness of G. s. subgutturosa from Azerbaijan compared to G. s. gracilicornis from Turkmenistan and Uzbekistan. Abduriyim, Zibibulla, et al. (2018) distinguished three subspecies in Asia: G. s. gracilicornis (referred to as G. s. subgutturosa in their study, but sequences from Iran and Azerbaijan were hardly included) stretching from Turkmenistan to eastern Kazakhstan, G. s. hillieriana (possibly including G. s. reginae and G. s. sairensis) in northern China and Mongolia, and G. s. yarkandensis in southern Xinjiang Uygur Autonomous Region (XUAR) in China. The majority of the recent studies, however, focused on population genetics of specific regions within the range of G. subgutturosa, for example XUAR in China (Abduriyim, Nabi, & Halik, 2018a;Abduriyim, Zibibulla, et al., 2018;Dong et al., 2016) or Iran (Fadakar et al., 2019;Khosravi et al., 2019;Mirzakhah et al., 2015;Zachos et al., 2010) without addressing taxonomic questions.

| Conservation and population status of goitered gazelles in Iran
Knowing more about the population genetics of G. subgutturosa in Iran is important to improve their conservation. The species is associated with open plains near hilly escape terrain that decreases their susceptibility to poachers (Farhadinia et al., 2009). It is widely distributed in all steppes or semideserts in Iran, except in the far northwest, along the Caspian Sea, and in the southeast (Firouz, 2005;Karami et al., 2002). Even though many areas exist with high numbers of resident G. subgutturosa, for example, Mooteh Wildlife Refuge (>7,000) and Qamishlou National Park (>3,000) in F I G U R E 1 Distribution of subspecies of G. subgutturosa in Asia based on the morphological data, cyt b sequences, and D-loop sequences. Blue (G. s. subgutturosa), yellow (G. s. gracilicornis), and red (G. s. yarkandensis) polygons correspond to the distribution based on the morphological data (Groves & Grubb, 2011) which are modified from distribution map of G. subgutturosa (IUCN SSC Antelope Specialist Group, 2017), large circles represent samples sequenced for cyt b, small circles represent samples sequenced for D-loop. Genetically identified subspecies (based on this study) are represented with blue (G. s. subgutturosa) and red (G. s. yarkandensis) circles. The locations of sequences of China and Mongolia are approximate. Different shades of gray represent different altitudes Isfahan Province, populations have declined considerably during the last four decades (Hemami & Groves, 2001;Karami et al., 2002).
The species has even disappeared from some protected areas, for example, Borouieh Wildlife Refuge (WR) in Central Iran (Fadakar et al., 2013), as well as many nonprotected areas, for example, the Chenaran plain of Razavi Khorasan Province and Turkmen Sahra of Golestan Province in northeastern Iran.
In the early 1990s, about 100,000 individuals of G. subgutturosa occurred in Asia, but now the species is threatened in many parts of its natural range (Mallon & Kingswood, 2001), which led to the IUCN classification "vulnerable" in 2008 (IUCN SSC Antelope Specialist Group, 2017). Hunting and habitat loss have caused a recent decline of more than 30% in many populations (IUCN SSC Antelope Specialist Group, 2017). In Iran, intensified hunting and habitat destruction due to overgrazing, conversion of natural gazelle habitat to agricultural land, urbanization, road development, and mining account for the dramatic decline of G. subgutturosa in deserts and plains (Fadakar et al., 2013(Fadakar et al., , 2019Karami et al., 2002;Mallon & Kingswood, 2001;Mirzakhah et al., 2015).
One recent study showed that a reintroduced population of G. subgutturosa in Dimeh Protected Area (PA) is actually a mixed population of G. subgutturosa and G. marica (Fadakar et al., 2019). These were considered to be conspecific based on external appearance and translocated from two source populations without prior knowledge on their genetic identity. Also, reintroduced goitered gazelle individuals from eastern Turkey to eastern Georgia showed haplotype identity of G. marica (Murtskhvaladze, Gurielidze, Kopaliani, & Tarkhnishvili, 2012). These illustrate how important it is to analyze genetic markers before carrying out further translocations of gazelles, especially in Iran where three species, G. subgutturosa, G.
bennettii, and G. marica, occur partially in sympatry (Groves, 1997;Groves & Harrison, 1967), that is, the geographic ranges of G. marica and G. subgutturosa overlap in southwestern Iran, and the geographic ranges of G. subgutturosa and G. bennettii meet in Yazd Province, Central Iran, where both species occur in neighboring areas (G. subgutturosa in Kalmand-Bahadoran PA, G. bennettii in Darre Anjir WR). Also, a resident population of G. marica in Mond PA is in close proximity to the habitat of G. bennettii in Nayband National Park (NP) northwest of the Persian Gulf.
Another important aspect for the in situ conservation of G. subgutturosa in Iran is to understand the connectivity of (sub-) populations for assessing the impact of inbreeding, or for evaluating which areas are most important for maintaining a healthy population. The Zagros mountain range that stretches through the country from the northwest to the southeast separates suitable gazelle habitats, and significant morphological differences between gazelle populations east and west of the mountain range are detectable (Hayatgheib et al., 2011). Contrastingly, in the central Iranian Plateau large groups of G. subgutturosa exist, that migrate between protected areas (Khosravi et al., 2018), for example, from Mooteh WR to Qamishlou NP in winter, and back to the Mooteh WR in summer (Fakheran Esfahani & Karami, 2005) and from Mooteh WR to northern areas such as Haftad Gholle PA. Therefore, genetic connectivity can be assumed between those populations.
With our study, we want to broaden the taxonomic picture for G. subgutturosa by sequencing the complete cytochrome b (cyt b) gene for this species in Iran. The new sequences are combined with published cyt b sequences of G. subgutturosa from Asia (Abduriyim, Nabi, et al., 2018;Dong et al., 2016) and other published mitochondrial sequences available on GenBank to study the phylogeny and subspecific patterning of G. subgutturosa. We aim to infer for the first time a genetic framework for Iranian G. subgutturosa populations living in almost all parts of the country, in order to understand the phylogenetic relationships among them, and with G. subgutturosa in other parts of Asia. Therefore, we analyzed the mitochondrial sequence variation of cyt b from 18 sampling sites. We hypothesized that (1) all supposed Iranian G. subgutturosa populations belong to this species and no further maternal introgression from G. marica occurred, and (2) the Zagros mountain range acts as a geographical barrier between the gazelle populations that occur east and west of the mountains, as was proposed based on morphological studies (Geptner, Nasimovich, & Bannikov, 1961;Hayatgheib et al., 2011).

| Sampling
This study was conducted with permission by the Iranian Department of Environment that authorized access to all sampling locations. In total, 60 fecal samples from G. subgutturosa were collected from 18 different localities in Iran ( Figure 2). Information regarding the origin, collector, and kind of material are summarized in Table S1 for each sample (Supporting Information). Fresh feces were collected in the field, after observing the animals from a distance to allow for species identification. Samples were stored in 96% ethanol.

| DNA extraction, amplification, and sequencing
Whole genomic DNA was extracted from fecal samples using The thermocycling conditions for L14724 and H15915 primers was performed as follows: initial denaturation (3 min at 95°C), followed by five cycle steps of 60 s at 94°C (denaturation), 90 s at 45°C (primer annealing) and 90 s at 72°C (elongation), then 40 cycle steps of 60 s at 94°C, 60 s at 50°C and 90 s at 72°C, and lastly, a final extension step (10 min at 72°C) (Lerp et al., 2011).
For CYTB_F and CYTB_R primers, we used the following protocol (Rezaei et al., 2010): 10 min at 95°C followed by 35 cycles of 30 s at 95°C, 30 s at 55°C, and 60 s at 72°C, and finally followed by 7 min at 72°C.
Sanger sequencing was performed using the BigDye Terminator Cycle Sequencing kit v.3.1 (Applied BioSystems), and electrophoresis of the purified sequencing product was carried out on an ABI PRISM 3730xl automatic sequencer.

| Alignments
For subsequent analyses, we worked with different alignments.
Alignment 1 is restricted to the 60 new Iranian sequences generated in this study (i.e., "Western group" in Tables 1 and 2). Alignment 2 (i.e., "Total group" in Tables 1 and 2) additionally includes 219 cyt b sequences (48 haplotypes) referred to as G. subgutturosa from GenBank (see Table S1 for accession numbers). The majority  Table S1]. The final dataset of complete cyt b (Alignment 2) comprises 279 sequences of G. subgutturosa and 8 outgroup sequences. Alignment 3 further includes 32 partial cyt b sequences of G. subgutturosa and is therefore restricted to only 390 bp.
The additional sequences come from our study (in cases where no complete cyt b fragment could be obtained) and from GenBank, restricted to samples with known geographic provenance (Fadakar et al., 2013;Khosravi et al., 2019;Mirzakhah et al., 2015). Sequences were aligned using the Clustal W algorithm (Thompson, Higgins, & Gibson, 1994) implemented in Mega v.5 (Tamura et al., 2011), and obvious misalignments were corrected by eye.

| Haplotype network
Based on Alignment 1, a median-joining (MJ) network was constructed for all 60 sequences from Iran using PopART v.1.7 (Leigh & Bryant, 2015) with the default settings. A second MJ network was generated, using 279 sequences (Alignment 2), including 60 sample sequences from Iran (i.e., Alignment 1 with 19 Haplotypes) plus 219 published sequences (48 haplotypes) for whole Asia. For this haplotype network, we had to work around the problem that 13 sequences (from Abduriyim, Nabi, et al., 2018) were only 995 bp in length. As we did not want to lose the information from the remaining 145 bp that show some mutations mostly in Iranian samples, we filled in the missing information for the shorter sequences using the consensus sequence of all full-length Asian samples (from Dong et al., 2016). Only three of these 32 sequences showed one unique single point mutation each in these parts of the sequence, so the consensus most likely reflects the true sequence information of the shorter sequences. The same was done for the sequence from Mongolia (Lerp et al., 2016) that has a length of 1,083 bp.
Models were selected by the Bayesian information criterion (BIC). We found the optimal partitioning scheme includes three partitions (optimal models are shown in brackets) 1st codon (K80), 2nd codon (HKY + I), and 3rd codon (GTR + G).
Phylogenetic analysis using Bayesian inference was carried out in MrBayes v.3.2.7a (Ronquist et al., 2012) with two independent runs of four Markov chains (one cold and three heated) over 10,000,000 generations and sampling every 1,000 generations. The first 25% of the sampled trees and estimated parameters were discarded as burn-in. Convergence onto the stationary distribution was monitored with average standard deviation of split frequencies (below 0.01), the potential scale reduction factor (close to 1 for all parameters) in MrBayes v.3.2.7a (Ronquist et al., 2012), and the effective sample size (ESS) value (above 200) in Tracer v.1.7.1 (Rambaut, Drummond, Xie, Baele, & Suchard, 2018 (Nei, 1987)), nucleotide diversity (π, the average number of nucleotide differences per site), and the mean number of pairwise differences within a group (k).
We calculated mismatch distributions separately for the phylogenetic subgroups and the total group to test if their frequency graph shows a chaotic/multimodal pattern characteristic for populations in demographic equilibrium, or a unimodal profile which is found in populations that experienced recent geographic expansion (Hey & Nielsen, 2004

| Additional sequencing of mitochondrial D-loop sequences
During the course of this study, we found it quite dissatisfying that for some already sampled populations of G. subgutturosa, especially in Azerbaijan, Turkmenistan, and Uzbekistan, cyt b sequences were not available. Restricting the study to this marker would limit the analysis of subspecies patterning of G. subgutturosa, as potentially an entire subspecies-G. s. gracilicornis-would be left out. We therefore decided to sequence ten of our samples for mitochondrial D-loop, as this marker has been used in other studies of G.
subgutturosa (Abduriyim, Zibibulla, et al., 2018;Sorokin et al., 2011;Zachos et al., 2010) (Tajima, 1989), and Fu's F s (Fu, 1997) F I G U R E 3 Median-joining network based on the complete cyt b gene depicting the relationships among the three clades described for G. subgutturosa in Iran (Asiatic, Central Iranian, and Middle Eastern), delimited by dashed lines. Circle sizes are proportional to haplotype frequencies; colors refer to sampling areas using numbers and abbreviations (compare Figure 2) Cen tr al Ir a n i a  Phylogenetic analysis was carried out similarly to the analysis of cyt b sequences (see above). A MJ network was also constructed using PopART v.1.7 (Leigh & Bryant, 2015) with the default settings.

| Haplotype network of cyt b
The were also detected in Central Iran, and H54 even in one location in northern Iran, east of the Caspian Sea (Golestan NP, Figure 2). In the Central Iranian clade (H62-H67), H62 was the most frequent haplo-  Figure 1).
When looking at the short segment (390 bp) of cyt b (Alignment 3), three additional haplotypes were found by Khosravi et al. (2019) in Central Iran, while one previously published short haplotype is identical with H54 of the complete cyt b, and another is identical with H55. In total, only 29 haplotypes can be distinguished based on the short fragment, so we will not discuss these results further.

| Phylogenetic analysis of cyt b
All G. subgutturosa sequences form a monophyletic clade (posterior probability (PP) = 1; Figure 5) which is placed as sister species

| Population differentiation and structure
In order to obtain an overview of the genetic diversity of G. subgutturosa, basic molecular diversity indices were determined for several subgroups of the species, based on geographic patterns (Western group and Eastern group) and phylogenetic clades according to the haplotype network and phylogeny ( Mismatch distributions were calculated separately for the Asiatic clade, for the combination of Central Iranian and Middle Eastern clades, and for the Total group that includes all sequences ( Figure 6).
The mismatch distributions for the Asiatic clade and the Total group were unimodal and fully consistent with a population expansion, but the combined Central Iranian and Middle Eastern clades showed a pattern that is interpreted to represent demographic equilibrium. A similar picture is found using R 2 statistic (Table 2), which was significant for the Eastern group (p < .01) and the Total group (p < .01), but not for the Western group. When looking at the clades, only the Asiatic clade shows values indicative of population expansion (p < .01). The same picture emerges when looking at Tajima's D (Table 2), but here a significant negative value indicative of population expansion was also found for the Central Iranian clade. Fu's F s values were highly significant also for the Eastern group, Total group, and Asiatic clade, and much smaller though still significant for the other groups.

| Phylogenetic analysis of D-loop sequences
Eight out of the ten samples were successfully sequenced for D-loop, including representatives for each northeastern Iranian haplotype and all six sequences from Turkmenistan form a monophyletic clade (PP = 1) that is separated by at least nine mutational steps.

| Geographic patterning and origin of G. subgutturosa based on cyt b
The haplotype networks (Figures 3 and 4)

G . s . s u b g u t t u r o s a G . s . y a r k a n d e n s i s
expansion is also reflected in the high number of rare alleles in the Asiatic clade that are arranged in a star-like pattern ( Figure 4) with only one or two mutational steps distance around a very frequent central haplotype (Slatkin & Hudson, 1991 Sea. This area was strongly affected by the Alpine orogeny (Tchernov, 1988) and especially by the glacial cycles during the Pleistocene when the first ancestors of G. subgutturosa might have emerged (Lerp et al., 2016). Rapid habitat changes might have promoted the subsequent range expansion in this group. The presence of a third subspecies, G. s. gracilicornis, was evaluated using D-loop sequences, as only this marker is currently available for samples from Turkmenistan and Uzbekistan where it is supposed to occur. The haplotype network ( Figure S1) and the phylogeny (Figure 7) show a clade of unique haplotypes from

| Population differentiation and taxonomic implications
Turkmenistan that might be interpreted as representing G. s. gracilicornis. However, these haplotypes were found in sympatry with other haplotypes that can clearly be assigned to G. s. yarkandensis as they are deeply nested within the Chinese haplotypes. So it is possible that Turkmenistan and Uzbekistan were colonized by G. subgutturosa several times. The first immigrants stem from the nominate form, G. s. subgutturosa, and seem to have been isolated from their conspecifics for some time, therefore evolving the unique haplotypes that were only found in Turkmenistan. However, at

| Population expansion
The mismatch distribution analysis clearly supports a demographic expansion for the whole G. subgutturosa species, and particularly for the G. s. yarkandensis subspecies (Figure 6). For G. s. subgutturosa the picture is not very clear, as the curve is relatively flat and does not show a marked peak.
R 2 statistics and Tajima's D support this interpretation ( has not significantly increased its range compared to the historic distribution and still occurs in the geographic area where it originated.
This picture might change when sequences from Azerbaijan are included, as it has been shown that the population in Shirvan Steppe Reserve inherits several unique haplotypes for mitochondrial control region sequences compared with sequences from Iran (Sorokin et al., 2011, Figure S1).

| Geographic barriers
Iran is a country with very diverse landscapes. The Zagros moun-

G . s . y a r k a n d e n s i s G. s. subgutturosa
spreading of G. subgutturosa from Iran further into Asia has happened a long time ago, when probably migration was possible due to more favorable climatic conditions at that time.
It has been proposed before that G. subgutturosa populations east and west of the Zagros mountain range are separate subspecies (Hayatgheib et al., 2011;Hemami, 1994;Karami et al., 2002). Hayatgheib et al. (2011) compared skulls of animals from both sides of the mountain range and found differences in at least 10 variables.
However, the three specimens (two males and one female) that they included from east of the Zagros mountain range were all collected in Khorasan Province in the northeast of Iran. Based on our genetic results (Figures 2 and 3   This is corroborated by morphological differences between the two groups (Geptner et al., 1961;Hayatgheib et al., 2011). We found no evidence to support the existence of a third subspecies, G. s. gracilicornis, in Turkmenistan and Uzbekistan based on D-loop sequences and conclude that gazelles in these areas belong to G. s. yarkandensis.

| Implications for conservation
For every future gazelle translocation action, for setting up captive breeding programs, or restocking wild populations, these clades need to be taken into account: Source populations and target populations should belong to the same clade. Most importantly, G.
subgutturosa from northeastern Iran must not be mixed with the remaining Iranian populations. As pure or hybrid populations of G.
marica could be present west of the Zagros mountain range, it is especially important to not relocate animals from west of Zagros to Central Iran. As long as the presence of G. marica in Iran and Iraq is not studied in detail, conservation efforts for G. subgutturosa should focus on Central Iranian populations. The viability of these populations will depend on increased persecution of illegal hunting, inside and outside of protected areas, so that stepping-stone habitats are maintained and the connectivity of all areas is improved to allow for gazelle migration in a protected area network. Gareth Jenkins, whose comments and suggestions helped to improve this manuscript, and especially encouraged us to include D-loop sequences as well. This study was supported by Gorgan University of Agricultural Science and Natural Resources, and the project funding 2018 of the German Society of Mammalian Biology (DGS).

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