Mitochondrial DNA analyses revealed low genetic diversity in the endangered Indian wild ass Equus hemionus khur

Abstract The Indian wild ass Equus hemionus khur, belonging to ass-like equid branch, inhabits the dry and arid desert of the Little Rann of Kutch, Gujarat. The E. h. khur is the sole survivor of Asiatic wild ass species/subspecies in South Asia. To provide first ever insights into the genetic diversity, phylogeny, and demography of the endangered Indian wild ass, we sampled 52 free-ranging individuals from the Little Rann of Kutch by using a non-invasive methodology. The sequencing of 230 bp in cytochrome b (Cyt b) and displacement loop (D-loop) region revealed that current ∼4000 extant population of Indian wild ass harbours low genetic diversity. Phylogenetic analyses confirmed that E. h. khur, E. h. onager, and E. h. kulan belong to a single strict monophyletic clade. Therefore, we suggest the delimitation of the five E. hemionus subspecies in vogue to a single species E. hemionus. The application of molecular clock confirmed that the Asiatic wild ass had undergone diversification 0.65 Million years ago. Demographic measurements assessed using a Bayesian skyline plot demonstrated decline in the maternal effective population size of the Indian wild ass during different periods; these periods coincided with the origin and rise of the Indus civilization in the northwest of the Indian subcontinent during the Neolithic. In conclusion, maintaining high genetic diversity in the existing isolated population of 4000 Indian wild asses inhabiting the wild ass sanctuary is important compared with subspecies preservation alone.


Introduction
Equids, belonging to the family Equidae, include asses, zebras, extant horses, and extinct horse-like species. Despite their vast fossil records and them being a classic example of macroevolution over many years, members of Equidae remain debatable amongst researchers (MacFadden 2005). Currently, IUCN has recognized eight extant species in the genus Equus: Horse (Equus caballus), Przewalski horse (E. przewalskii), African wild ass (E. africanus), kiang (E. kiang), Asiatic wild ass (E. hemionus), plains zebra (E. quagga), mountain zebra (E. zebra), and Grevy's zebra (E. grevyi) (Moehlman 2002). Because of its rapid radiation and recent divergence, the genus Equus has been extensively studied for addressing questions mainly related to its phylogeny. Most previous studies on equids were based on few samples obtained mainly from captive individuals from zoos and other sources having limited information regarding their exact geographical origin (Oakenfull et al. 2000;Steiner & Ryder 2011;Steiner et al. 2012). Among all equid branches, the ass-like branch of equids comprising of the Asiatic and African wild asses remains largely unresolved and requires attention. Rapid climate change and increased anthropogenic pressure on their habitats have fragmented the current distribution of all Asiatic wild asses. On the basis of their morphology and karyotype, two main species of the Asiatic wild ass, E. hemionus and E. kiang, have been described (Groves & Maz ak 1967;Ryder & Chemnick 1990); however, molecular data failed to distinguish these species (McCue et al. 2008;Vilstrup et al. 2013 (Groves & Maz ak 1967;Schaller 1998;Shah 2002). Similarly, three subspecies of E. kiang, E. k. kiang, E. k. holdereri, and E. k. polydon, have been recognized (Shah 2002;Grubb 2005). This study focuses on one such Asiatic wild ass population from western India called E. h.khur. The E.h.khur inhabits the highly fragmented, dry, and arid desert of the Little Rann of Kutch, Gujarat, India ( Figure 1). Moreover, it is the sole survivor of the once widespread Asiatic wild ass population, showing continuous distribution from the Arabian peninsula to Manchuria, and declared endangered by the IUCN (Grubb 2005). Equus hemionus khur is the flagship species of the costal desert in western India, a unique ecosystem spread into an area of 9000 km 2 , encompassing the Great and Little Rann of [68][69][70][71]. In the past, E. h. khur was distributed from southern Pakistan (Sindh and Baluchistan provinces) and Afghanistan to southeastern Iran. However, the Indian wild ass sanctuary (4900 km 2 ), located in the Little and Great Rann of Kutch Gujarat, is currently the last refugia of the Asiatic wild ass population in southern Asia. Although construction of the Indian wild ass sanctuary in this saltimpregnated desert has provided shelter to these endangered Asiatic wild asses, but the constant anthropogenic pressure such as salt production, grazing, habitat destruction, and dam construction in and around the sanctuary has threatened their existence (Goyal et al. 1999;Srivastav & Nigam 2010). Recent Gujarat forest reports have confirmed an increase in both population (over 4000) and distribution range of E. h. khur in western India; however, in 1967, a dreaded protozoic disease called Surra caused by Trypanosoma evansi reduced the population of the wild equid to mere 362 individuals (Gujarat Forest Department 1967;Singh 2000;Srivastav & Nigam 2010). Most previous studies on this free-ranging Indian wild ass were confined to demographic or ecological measurements and behaviours without assessing their genetic diversity (Gee 1963;Shah & Qureshi 2007). The only study that evaluated the genetic diversity of E. h. khur was conducted by Srivastav and Nigam (2010) with few captive individuals from a zoo. Recent improvements in the non-invasive conservation genetics methods have enabled researchers to study the diversity of the elusive and protected wild species worldwide (Beja-Pereira et al. 2009;Allendorf et al. 2010). We conducted the first pilot study for (1) measuring the level of genetic diversity underlying the natural Indian wild ass population in the Little Rann of Kutch, (2) confirming the position of E. h. khur in the ass-like branch of equid phylogeny, and (3) estimating the divergence and demography of these endangered wild asses by using a non-invasive sampling approach.

Samples
Over 12 months, 52 faecal samples of the free-ranging Indian wild ass E. h. khur were collected from three localities in the sanctuary ( Figure 1). The samples were placed in individual zippered bags, dried naturally, and stored at room temperature until further processing. Permission for collecting faecal samples in the wild ass sanctuary was obtained from the state of Gujarat through a letter (No. WLP/28/C/574-76/2013-14; Date: 18/12/2013) issued by the Principal Chief Conservator of Forest (Wildlife), Gujarat.

DNA extraction, PCR amplification, and sequencing
Genomic DNA was extracted using the modified protocol of a spin column tissue DNA extraction kit (GENETIX Miniprep Kit, GENETIX, Rome, GA). Modification to the standard protocol included long-term incubation of the fecal samples in suitable concentrations of lysis buffer, along with the addition of Proteinase K and use of the inhibitEX V R tablet (QIAGEN, Valencia, CA) for preventing PCR inhibition. Khurcytb FOR (5 0 -GAC ACA ACA ACC GCC TTC TC-3 0 ) and Khurcytb REV (5 0 -GAC TGT TGC TCC TCA GAA GGA T-3 0 ) were used to amplify approximately 230 bp of the cytochrome b (Cyt b) fragment in 33 samples, whereas KhurDL FOR (5 0 -TCC CCA TGT GCT ATG TCA GT-3 0 ) and KhurDL REV (5 0 -GAT ARG CGT GTT GAC TGG AAA-3 0 ) were used to amplify approximately 230 bp of displacement loop (D-loop) in 52 samples. PCR reactions were performed using 20 lL of total reaction mixture composed of 3 lL of deionized water, 8 lL of Dream Taq Mastermix (Thermo Scientific, Waltham, MA), 1.5 lL of 2.5 mM MgCl 2 , 1 lL of 0.3% BSA, 0.5 lL (5U/lL) Taq polymerase (Thermo Scientific, Waltham, MA), 0.8 lM of each primer, and variable amounts of genomic DNA. PCR cycling for both loci was performed in a PCR Verity 96-well thermal cycler (Applied Biosystems, Waltham, MA). For the Cyt b gene, the PCR mixture was subjected to initial denaturation at 94 C for 15 min, followed by 45 cycles at 94 C for 50 s, 60 C for 50 s (annealing), and 72 C for 45 s (extension), and final extension at 72 C for 20 min. For D-loop amplification, identical cycling conditions were used, except for the annealing temperature, which was 61.7 C. After purification with ExoSAP-IT V R (Affymetrix), amplicons were sequenced for both strands on an automated DNA sequencer ABI PRISM V R 377 (Applied Biosystems, Waltham, MA).

Molecular phylogeny, dating, and demography
Chromatogram files were edited and aligned using Geneious 8.0.5 (Applied Biosystems, Waltham, MA). All 33 Cyt b and 52 D-loop partial sequences generated were deposited in Genbank with accession nos. KJ490647-KJ490651 and KT221803-KT221829, KU342004-KU342014, respectively. To clarify the relative position of E. h. khur D-loop and Cyt b sequences generated using the available wild ass species/subspecies sequences in the database (Supplementary Data), we constructed phylogenetic trees. The best-fit model of evolution was selected using jModelTest (version 2.1.3) (Guindon & Gascuel 2003;Darriba et al. 2012) according to the Akaike information criterion for both D-loop and Cyt b datasets. The phylogenetic trees were inferred using the Bayesian approach as implemented in MrBAYES (version 3.2.1) (Ronquist et al. 2012). To obtain the Bayesian-inferred trees, two independent analyses were performed, and four MCMC chains were run for 20 million generations with sampling at every 1000 generations. After checking convergence, 25% trees were discarded as burn-in. The molecular divergence of E. h. khur was measured using BEAST (version 2.3.1) by estimating time to the most common ancestor among equid species for both datasets (Drummond et al. 2012). A relaxed uncorrelated lognormal clock model and a Yule process of speciation were used. The Asiatic wild ass branch was calibrated using several internal calibration points drawn from the study of Vilstrup et al. (2013) and Rosenbom et al. (2015) for extant equid species. Apart from E. hemionus and E. kiang, monophyly was constrained for species represented in the analyses by more than one individual. While constraining monophyly, normal prior with HKY þ G þ I and GTRþ G þI models were used for Cyt b and D-loop datasets, respectively. Parameters were sampled at every 1000 generations, convergence was viewed using TRACER (version 1.5) (Rambaut & Drummond 2007), and four MCMC chains over a total of 20 million generations were run with 25% generations discarded as burn-in. The program BEAST (version 2.3.1) was used to measure the present and past demography of E. h. khur by using the Bayesian skyline model that estimates the posterior distribution of population sizes (Drummond et al. 2005). A strict molecular clock with normal distribution prior with a clock rate of 3. 6 Â 10 À 8 per generation, based on the studies of the mitochondrial DNA control region of domestic donkeys and African wild asses (Kimura et al. 2011;Rosenbom et al. 2015) was applied.

Diversity and phylogeny
The alignment of Cyt b and D-loop regions revealed only one haplotype, indicating low genetic diversity in the extant Indian wild ass population. Homology searching using BLAST 2.2.32 (Zhang et al. 2000;Aleksandr et al. 2008) showed that the E. h. khur haplotype shares 98% identity with the two E. h. kulan sequences (JX312728.1 and JF718886.1). The Bayesian-inferred tree based on Cyt b and D-loop haplotypes demonstrated two well-supported clades within the Asiatic wild ass E. hemionus. In Cyt b haplotype-based BI tree, E h. khur cluster with E. h. kulan with high posterior support while in another clade which include, Equus kiang, E. h. hemionus and E. h. onager, we failed to recover the relatedness among these species/subspecies owing to resulted polytomy (Figure 2(A)). In addition, similar to the Cyt b-inferred tree, the tree constructed on the basis of the Dloop sequences also show two well-supported clades, where Indian wild ass, E. h. khur cluster with clade formed by E. h. kulan and E. h. onager while in another well-supported clade, the monophyly of Equus kiang with E. h. hemionus was illustrated (Figure 2(B)). The measurement of time to the most common ancestor in the presence of the Indian wild ass haplotype demonstrated that the Asiatic wild ass branch diverged approximately 0.65 million years ago (Figure 2(B)). The estimation of demographic dynamics by using the Bayesian skyline model revealed that the Indian wild ass maintained a steady maternal effective population size (N*g) until 10 kya, after which a population decline was detected (Figure 3).

Discussion
The occurrence of a single haplotype in the sampled population strengthens the 1967 Gujarat forest reports, revealing that severe decline in the E. h. khur population in 1962 caused by the protozoic disease Surra in the state. Phylogeny reconstruction of the Asiatic wild ass in the presence of the firstreported Indian wild haplotype is in agreement with the phylogeny inferred by Rosenbom et al. (2015), where two major clades of Asiatic wild asses were identified without any further distinction of subspecies within these clades. These findings were not consistent with the results of other previous studies, where a single monophyletic clade was reported for the entire E. hemionus species/subspecies in which E. kiang haplotypes formed the monophyletic clade within a larger variation of E. hemionus (Steiner et al. 2012;Vilstrup et al. 2013). Therefore, our phylogenetic analyses confirmed that the Indian wild ass, E. h. khur is genetically similar to the endangered E. h. kulan and E. h. onager, commonly called as Turkmenian kulan and Persian onager, respectively. However, the surviving population of approximately 4000 individuals of the Indian wild ass E. h. khur in the Indian wild ass sanctuary is the relic of once widespread but currently critically endangered population of E. hemionus or Persian onager in Asia (http://www.IUCN redlist.org). The successful reintroduction of E. h. kulan to Kazakhstan and Uzbekistan during the Soviet Union times, followed by its reintroduction in Israel, where it is forming hybrids with Persian onager in wild, strengthens this finding (Renan et al. 2015). The time depth of 0.65 million years estimated for the Asian wild ass branch is in agreement with the recent studies on wild asses conducted by Rosenbom et al. (2015) and Vilstrup et al. (2013). These results corroborate that climatic events during the Pleistocene were the major driving force in the differentiation processes in wild asses. A decrease in the maternal effective population size of the Indian wild ass approximately after the 10 kya coincides with the rise of 8000-year-old Neolithic human settlement in Mehrgarh in the northwest of the subcontinent eventually culminating into the 5000-year-old Indus civilization (Possehl Gregory 1996;Kenoyer 1998). However, the expansion of the human population along with hunting and habitat loss must have continued to impose a decreasing trend until recently. A similar decline in effective population sizes has been reported for E. hemionus and wild horses that showed a decline in their genetic diversity after the Last Glacial Maximum (LGM) (Der Sarkissian et al. 2015;Rosenbom et al. 2015).

Conclusion
By using a non-invasive methodology, we demonstrated that the Indian wild ass E. h. khur is genetically similar to E. h. kulan and E. h. onager, two subspecies of the Asiatic wild ass. Despite the current population of over 4000 individuals, the maternal effective population size of the Indian wild ass has declined, resulting in low genetic diversity. In addition to invoking the revision of systematic of Asiatic wild ass species/ subspecies, this pilot study can be a stepping stone for the genetic conservation of the Indian wild ass.