Repeated hybridization of two closely related gazelle species (Gazella bennettii and Gazella subgutturosa) in central Iran

Abstract Interspecific hybridization increasingly occurs in the course of anthropogenic actions, such as species translocations and introductions, and habitat modifications or occurs in sympatric species due to the shortage of conspecific mates. Compared with anthropogenically caused hybridization, natural hybridization is more difficult to prove, but both play an important role in conservation. In this study, we detected hybridization of two gazelle sister species, Gazella bennettii (adapted to dry areas) and Gazella subgutturosa (adapted to open plains), in five habitat areas, where G. bennettii naturally occur in central Iran. The hybrids have a nuclear genomic identity (based on two introns), habitat preference, and phenotype of G. bennettii, but the mitochondrial identity (based on cyt b) of G. subgutturosa. We suggest that natural hybridization of female G. subgutturosa and male G. bennettii happened twice in central Iran in prehistoric times, based on the haplotype pattern that we found. However, we found indications of recent hybridization between both species under special circumstances, for example, in breeding centers, due to translocations, or in areas of sympatry due to the shortage of conspecific mates. Therefore, these two species must be kept separately in the breeding centers, and introduction of one of them into the habitat of the other must be strictly avoided.

In the genus Gazella, the only record of hybridization until now is the introgression of mitochondrial DNA between G. marica and G. subgutturosa in eastern Turkey (Murtskhvaladze, Gurielidze, Kopaliani, & Tarkhnishvili, 2012). Morphologically, intermediate individuals, probably natural hybrids between the two nominal forms (Groves, 1997;Groves & Harrison, 1967;Mallon & Kingswood, 2001), occur in a vast area between the Tigris/Euphrates valley and the Zagros Mountains, but the only molecular study that included samples from this region (five captive individuals from the Rutba region in Iraq -Wacher et al., 2010) only detected mitochondrial sequences of G. subgutturosa. Detection of natural or anthropogenic hybridization is often difficult, but important for conservation actions because of the complex situation of policies and management decision for hybrids (Allendorf et al., 2001;Ellstrand et al., 2010;Stronen & Paquet, 2013;Trouwborst, 2014).
Iran is home to three, possibly four gazelle species and can be regarded as a hotspot for gazelle diversity. Arabian mountain gazelles (Gazella arabica dareshurii) are restricted to Farur Island in the Persian Gulf. In southwestern Iran, populations of sand gazelles (G. marica) exist based on mitochondrial cytochrome b (cyt b) sequences (Fadakar et al., 2019), but morphologically these individuals are intermediate in size and morphology between G. marica and G. subgutturosa and might be natural hybrids (Groves & Harrison, 1967), which remains to be confirmed using nuclear DNA. Goitered gazelles (G. subgutturosa) and chinkara (G. bennettii) occur in large parts of the country. Their ranges meet in central Iran, where both species occur in neighboring areas (e.g., G. subgutturosa in Kalmand-Bahadoran Protected Area and G. bennettii in Darre Anjir Wildlife Refuge).
Gazella subgutturosa is associated with open plains and is widely distributed in all steppes or semi-deserts, occurring in a relatively high number of individuals (Firouz, 2005;Karami, Hemami, & Groves, 2002). Its sister species, G. bennettii, adapted to dry deserts with the ability to survive with very little brackish water in hot weather in the deserts of central Iran, such as Dasht-e Lut and Dasht-e Kavir deserts, with a relatively low number of individuals (often less than 50) in each area (Akbari, Moradi, Sarhangzadeh, & Esfandabad, 2014).
The two species can easily be differentiated, as females of G. subgutturosa are usually hornless (some females may have short and deformed horns), while G. bennettii females have long, slender horns. In males, the horns of G. subgutturosa are much more out bowed, while males of G. bennettii have relatively straight horns. In both sexes, the body size of G. bennettii is much smaller than G. subgutturosa.
Due to the highly sculptured landscape, intraspecific patterning was detected in G. subgutturosa (Khosravi et al., 2018), where two subspecies were found to occur in Iran .
G. bennettii was also supposed to split into at least two subspecies or even species in Iran (Groves & Grubb, 2011), but this has not been confirmed using molecular data yet. For a study on the genetic diversity of G. bennettii, we therefore collected fecal and tissue samples of this species throughout the country (Fadakar, unpublished data). In a first step, these samples were sequenced for mitochondrial DNA (cyt b) to confirm species identification.
Surprisingly, however, some of them were identified as belonging to G. subgutturosa, despite being collected in habitats of G. bennettii and/or from individuals that were morphologically identified as belonging to G. bennettii (see Methods). Therefore, we hypothesize that these samples with G. subgutturosa mitochondrial identity belong to hybrids of G. subgutturosa and G. bennettii. Using additional sequences for two nuclear introns, we aimed at investigating this hypothesis. Furthermore, we asked whether hybridization occurred repeatedly or only once, and whether this is an ongoing process or a prehistoric event.

| Sampling
In total, 32 samples, collected in G. bennettii habitats and stored in 96% ethanol, are included in this study (Table 1 and  in Bahabad desert. The desert is located between Bahabad city in the south, Naybandan WR in the east, Robat-e Posht-e Badam in the north, and Ardakan city in the west. One tissue sample was collected from a dead animal in Darre Anjri WR (Figure 2g), and one from a carcass found in Bazman HPA (southeastern Iran). All fresh feces were TA B L E 1 List of gazelle samples used in this study, with cyt b haplotype (Hap -see Table S1), GenBank accession numbers for cyt b and nuclear intron (ZNF618 and CHD2) sequences, sample location, and type of sample   The thermocycling 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.
Sequences were edited for correction with SeqScape v.2.6 software (Applied Biosystems). All new sequences were preliminarily identified using a BLAST search against known gazelle sequences on GenBank and were submitted to GenBank (MT811607-MT811638, Table 1).

| Nuclear markers
For a phylogenetic analysis of the genus Gazella, Lerp et al. (2016) published a new set of nuclear intron markers. These had only limited variation, but showed very high consistency in the delineation of species. Only two of the six markers (chromodomain-helicase-DNA-binding protein 2 (CHD2) and zinc finger protein 618 (ZNF618)) were able to differentiate between G. subgutturosa and G.  (Thompson, Higgins, & Gibson, 1994) implemented in Mega v.5 (Tamura et al., 2011). New sequences were submitted to GenBank (MT822208-MT822269, Table 1).

| Haplotype network for cyt b
A median-joining (MJ) network was constructed for 295 sequences,

| Phylogenetic analysis
For the phylogenetic analysis of mitochondrial cyt b sequences, we used sequences of all haplotypes collected from the habitats of G. bennettii, as well as one previously published sequence of G. bennettii (NC020703) and all haplotype sequences of G. subgutturosa from GenBank. Four sequences from closely related sand gazelle (G. marica) and slender-horned gazelle (Gazella leptoceros) were used as outgroup representatives (Table S1). Sequences were aligned using the Clustal W algorithm (Thompson et al., 1994) implemented F I G U R E 2 Hybrid samples in the breeding center of Naybandan WR (a-d), typical Gazella subgutturosa of the breeding center of Naybandan WR (e), typical pure Gazella bennettii in Bahou Kalat WR at southeast Iran (f), and hybrid sample from Darre Anjir WR which died due to snakebite (g). (e and f) are presented here to show phenotypical trait differentiation between typical pure G. bennettii and G. subgutturosa in Mega v.5 (Tamura et al., 2011), and final adjustments were made by eye.
The first 25% of the sampled trees and estimated parameters were discarded as burn-in. Convergence of the model parameters was monitored using the program Tracer v.

| Sequencing and mitochondrial cyt b
Extracting and amplifying mitochondrial and even nuclear DNA sequences from fecal samples was unproblematic, but from some tissue samples from Naybandan WR we obtained no results, presumably because of improper storing conditions (in a shed) for a long time. Using cyt b, 16 of 32 samples of putative pure G. bennettii were identified as G. subgutturosa in a BLAST search against known gazelle sequences on GenBank. These are therefore hypothesized to belong to hybrids of G. bennettii and G. subgutturosa in the subsequent analyses.

| Nuclear intron markers
Sequencing of CHD2 (652 bp) and ZNF618 (676 bp For CHD2 the only difference between G. bennettii and G. subgutturosa is a mutation at position 221 ("C" in G. subgutturosa, G. marica, G. leptoceros, and G. cuvieri; "T" in G. bennettii and all other species). This is also the case in our pure G. subgutturosa samples, and in all pure G. bennettii samples (see Table 1). All sequences from putative hybrids were identified as G. bennettii for CHD2. five samples are homozygotic for the indel, two are heterozygotic (nine unknown -see Table 2). Heterozygotic indels are inferred from the raw sequencing reads: If the forward and reverse sequences can clearly be read up to the indel sequence but not across, the indel is heterozygotic (the two alleles produced PCR products with different lengths, and therefore, the sequencing results are unreadable from the point onwards where they start to differ).

| Haplotype network
The reconstructed MJ network based on complete cyt b provides an overview of the haplotype distribution and relationships within the two G. subgutturosa subspecies and the hybrid haplotypes ( Figure 3).
Samples from Bahabad HPA, Darre Anjir WR, Ariz HPA (all three are located in Yazd Province), and Khabr NP (Kerman Province) are connected with (or even identical to) H54, the supposed ancestral haplotype of G. subgutturosa that belongs to the nominate subspecies G. subgutturosa subgutturosa . Samples from Naybandan WR (South Khorasan Province) derived from H1, the most frequent haplotype of the G. s. yarkandensis subspecies ( Figure 3). Of the eight haplotypes found in hybrids, only one has been detected in pure G. subgutturosa specimens: the ancestral haplotype H54 is present in hybrids in Khabr NP. All other hybrid haplotypes are new (H68-H74).

| Phylogeny
The

| Hybridization and chromosome numbers
Hybridization in gazelles is probably restricted to very closely related species, as Robertsonian translocations are very frequent in the genome of gazelles (Vassart, Granjon, & Greth, 1995). These are major chromosomal changes where two acrocentric chromosomes fuse to form one biarmed chromosome. If individuals with two different chromosome configurations interbreed, the offspring may have reduced fertility as the chromosomes cannot properly be segregated during meiosis, leading to unbalanced gametes (Baker & Bickham, 1986;Benirschke & Kumamoto, 1991). The ancestral condition in bovids probably is an autosome number of 2n = 58 (not including sex chromosomes - Wurster & Benirschke, 1968). While most of the ancestral autosomes are still identifiable in gazelles using banding techniques, many of them have fused to form biarmed chromosomes (including a common X-to-autosome-5 fusion), so the chromosome numbers in the genus Gazella are as low as 2n = 30-35 (Groves & Grubb, 2011;Vassart et al., 1995 and references therein).

| Hybridization as a prehistoric event
Hybridization can be a prehistoric phenomenon that might lead to mitochondrial capture of a whole species, for example, as demonstrated for polar bears (Ursus maritimus -Hailer et al., 2012) or European bison (Bison bonasus) which is very closely related to American bison (Bison bison) based on nuclear DNA, but seems to have adapted the mitochondrial identity of some extinct bovine species from Europe (Hassanin, An, Ropiquet, Nguyen, & Couloux, 2013).
Also, natural hybridization can occur repeatedly, for example, among Gazella subgutturosa experienced a strong bottleneck in Iran during the last decades (Hemami & Groves, 2001;Khosravi et al., 2018), so it is possible that the haplotypes involved in hybridization have been extirpated from the pure G. subgutturosa populations. Furthermore, small populations of G. subgutturosa still exist in unprotected and protected areas not sampled so far. Some of these are relatively close to the areas where hybrids were confirmed, for example, Biduiyeh PA (Kerman Province, approximate population size = 150) from where individuals might be able to migrate to Khabr NP (Khosravi et al., 2018). It could therefore be possible that hybrid haplotypes are present in small G. subgutturosa populations not sampled so far. This should be tested by sampling more individuals from hitherto unsampled small populations close to the hybrid populations.
With the current data, however, we conclude that the populations evolved independently after the initial hybridization events and that there was no back-crossing with G. subgutturosa (at least not in the maternal line). Especially in northeastern Iran, it is interesting to see that the hybrid specimens have haplotypes (H68 and H69) that are only two to three mutational steps away from the proposed ancestral haplotype of the G. s. yarkandensis subspecies (H1, not present in Iran but reported from China ), but four to five steps away from the pure G. s. yarkandensis individuals (H49-H53) also living in Khorasan Province (Figure 3). So the natural hybridization occurred relatively early, probably when the subspecies were not even properly distinct.
With regard to G. bennettii, it might even be concluded that the hybridization took place before the chromosome numbers started to change drastically in this species. G. subgutturosa and G. bennettii are sister species, so it can be assumed that their chromosome numbers were initially identical, which would facilitate hybridization. If

| Recent natural hybridization
Hybridization of G. bennettii and G. subgutturosa might still be possible in the wild, in areas where both species are sympatric with low individual numbers. In Bahabad HPA, for example, the total number of gazelle individuals (G. bennettii and G. subgutturosa) probably is below 10. Under these conditions, mating between the two species could take place as conspecific mates are not present or hard to find, based on the "desperation hypothesis" (Hubbs, 1955). This could potentially lead to strong outbreeding depression (Randi, 2008), as hybrids might have reduced fertility. The hybrids in Khabr NP could potentially be explained by the "desperation hypothesis": Khabr NP is located far away from the other G. subgutturosa populations sampled by Fadakar et al. (2020), but previously there has been a small population of G. subgutturosa in the northern part of Khabr NP. In order to increase the population of G. bennettii, the core area of the National Park was fenced about three decades ago. After fencing, the population of G. bennettii increased sharply (from about 70 to over 1,000), while the G. subgutturosa population decreased dramatically because their habitat was outside the fenced area of the National Park and severely under the pressure of illegal hunting, competition with livestock and habitat destruction which eventually led to their extinction. However, within the fenced area, remnant individuals of G. subgutturosa might have mixed with the abundant population of G. bennettii. This would have been a rather recent hybridization, which could explain why there is still a recent haplotype of G. subgutturosa (H54) present in the hybrids. However, as the G. bennettii population has been strongly increasing, there seems to be no negative effect of this hybridization on the reproduction and viability of the mixed population, which would be expected according to the very different haplotype numbers of the two parent species. Other reported cases of interspecific hybridization in bovids have shown that even if the parent species have different chromosome numbers, female offspring can be fertile, for example, in an Alcelaphus buselaphus × Damaliscus lunatus cross reported by Robinson et al. (2015).
So maybe one fertile female hybrid was enough to keep the subgutturosa haplotype in the G. bennettii population. In any case, an investigation in chromosome numbers would be desirable especially for individuals in this mixed population.
If, however, recent hybridization occurs repeatedly, we expect to see a negative effect on the population viability of the (hybrid-)G. bennettii population. One potential case of this is Darre Anjir WR, and other protected areas in central Iran. Although the level of protection has increased, there, the population numbers of G. bennettii (and "old" hybrids), have not increased (Akbari et al., 2014).

| Recent hybridization due to the anthropogenic actions
There has been at least one instance of recent hybridization of the two species due to anthropogenic action. In Shir Ahmad PA (Razavi Khorasan Province, Sabzewar), three females and one male G. bennettii, which were previously kept in a 25 ha fenced area in Shir Ahmad breeding center, were released to the protected area in 2006. This area is a prime habitat for G. subgutturosa in northeastern Iran with around 650 individuals of that species living there. Three of the four G. bennettii individuals died in the same year and only one female survived that was frequently observed by the game wardens (A. Khani, personal communication). In the following year that female was seen with a young. As no natural populations of G. bennettii live anywhere close to Shir Ahmad PA, the juvenile very likely is the offspring of the G. bennettii female and a G. subgutturosa male. Other G. bennettii individuals from the Shir Ahmad breeding center and their offspring were translocated to Salami Breeding Center (South Khorasan Province, Khaf), from where they were introduced to South Khorasan (Ferdous) in 2016. It is possible, that these gazelles were also in contact with G. subgutturosa, as some individuals of the latter could have been placed within the fenced area in Shir Ahmad (for short periods of time) for recovery from illness or after having been confiscated. Therefore, introductions and translocations, which can substantially increase the rate of hybridization in mammals in general (Allendorf et al., 2001;Rhymer & Simberloff, 1996;Vonlanthen et al., 2012) might also have led to hybridization between G. bennettii and G. subgutturosa.

| CON CLUS ION
Finding gazelles with mixed ancestry of G. subgutturosa (mitochondrial DNA) and G. bennettii (nuclear DNA) was very unexpected, as hybridization between these two sister species was never hypothesized before. The hybrid populations are located in central Iran, in a contact zone between the two species. G. subgutturosa is predominantly living in larger herds and migrating between habitats, while G. bennettii is more sedentary with animals living in small groups (Lerp, Wronski, Butynski, & Plath, 2013). The hybrids are phenotypically and ecologically identical with G. bennettii and most likely belong to the same population as pure G. bennettii in these areas.
Especially, Khabr NP with more than 1,000 individuals is a prime habitat of G. bennettii in Iran and acts as a source for smaller populations of G. bennettii in the surrounding habitats. Although we propose that the hybridization goes back to two separate prehistoric events, as no shared haplotypes exist between hybrid and neighboring pure populations of G. subgutturosa, it is not entirely impossible that hybridization between G. bennettii and G. subgutturosa still occurs in areas with very low population numbers. The two species now have very different chromosome numbers, so it is possible that hybrids from recent hybridization events have reduced fertility. Therefore, keeping G. subgutturosa and G. bennettii individuals in the same fenced area in a breeding center and also introducing G. bennettii to the habitat of G. subgutturosa such as Shir Ahmad PA should be avoided under all circumstances.
For future captive breeding programs, the knowledge of these newly identified wild hybrid populations is very important. As long as the genetic makeup of these animals, especially the chromosome numbers, is not known, we strongly advise against using them in G. bennettii conservation programs.

ACK N OWLED G M ENTS
We thank the Iranian Department of Environment for sampling authorizations. We also thank Ali Khani (former manager of Shir

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