Phylogeographical analysis shows the need to protect the wild yaks' last refuge in Nepal

Abstract The wild yak Bos mutus was believed to be regionally extinct in Nepal for decades until our team documented two individuals from Upper Humla, north‐western Nepal, in 2014. The International Union for Conservation of Nature (IUCN) seeks further evidence for the conclusive confirmation of that sighting. We conducted line transects and opportunistic sign surveys in the potential wild yak habitats of Humla, Dolpa, and Mustang districts between 2015 and 2017 and collected genetic samples (present and historic) of wild and domestic yaks Bos grunniens. We also sighted another wild yak in Upper Humla in 2015. Phylogenetic and haplotype network analyses based on mitochondrial D‐loop sequences (~450 bp) revealed that wild yaks in Humla share the haplotype with wild yaks from the north‐western region of the Qinghai‐Tibetan Plateau in China. While hybridization with domestic yaks is a major long‐term threat, illegal hunting for meat and trophy put the very small populations of wild yaks in Nepal at risk. Our study indicates that the unprotected habitat of Upper Humla is the last refuge for wild yaks in Nepal. We recommend wild yak conservation efforts in the country to focus on Upper Humla by (i) assigning a formal status of protected area to the region, (ii) raising awareness in the local communities for wild yak conservation, and (iii) providing support for adaptation of herding practice and pastureland use to ensure the viability of the population.

Kush Himalaya and across the QTP in China (Zhao & Gao, 1991 cited in Han, 2014). Schaller & Liu (1996) updated the distribution to include the Tibetan Plateau habitats in China, India, and Nepal, whereas (Harris & Leslie, 2008) considered the species to be present only in China and India, while regionally extinct in Bhutan and Nepal. At present, around 10,000-15,000 wild yaks roam QTP in an area of approximately 400,000 km 2 (Buzzard & Berger, 2016;Schaller, 1998;Zhang et al., 2020).
For many decades, only indirect evidence like horns, skulls, and pelts of presumed wild yaks indicated their past and/or current presence in Nepal. The Transhimalayan wild yak habitats in Nepal are contiguous with the Tibetan Plateau of the Tibetan Autonomous Region (TAR) in China as they represent parts of the western end of the Plateau. This allows seasonal movement of wild yaks from TAR into northern Nepal (Miller et al., 1994). However, due to the lack of evidence of live animals, Jnawali et al. (2011) assessed the species as "data deficient" and "possibly regionally extinct" in Nepal. Wild yaks are a protected priority species in the country (GoN, 1973) and The retaliatory killing of wild yaks (usually bulls) by livestock herders to prevent hybrid offspring and abductions of female domestic yaks Bos grunniens seems the most evident threat to wild yaks in Nepal (F. Tamang and G. Lama pers comm 2015). The herders and their domestic livestock threaten the wild yaks further through habitat encroachment and displacement (Harris, 2007), while the possibility of disease transmission between wild-domestic yak interfaces remains as a substantial threat (Buzzard & Berger, 2016;Schaller & Liu, 1996). Illegal hunting for meat and trophy adds to the threats (F.

Tamang and G. Lama pers comm 2015).
In 2014, we sighted two wild yaks in the remote Transhimalayan valleys of Upper Humla in north-western Nepal, leading to the rediscovery of the species in the country (Acharya et al., 2015). National and international experts identified the sighted wild yaks by referring to their morphology and behavior from our photographs and video footage. Some experts also suggested to perform genetic analyses for conclusive identification. As such, IUCN considers our sighting record as uncertain and seeks for an additional evidence for conclusive confirmation (Buzzard & Berger, 2016).
The mitochondrial D-loop or control region has been widely used for investigating intraspecific genetic variation, population structure, and demographic histories of animal domestication (Beja-pereira et al., 2004;Guo et al., 2006;Jansen et al., 2002;Lai et al., 2007;Larson et al., 2005;Luikart et al., 2001;Troy et al., 2001;Wolf, 1999). This is because D-loop is the most variable region of mtDNA having high nucleotide polymorphism compared to the rest of the regions in mtDNA. For example, the variable sites in the D-loop of yak mtDNA are seven times more than that observed in the entire yak mtDNA (Wang et al., 2010). The genetic diversity of domestic and wild yaks in the QTP has been inferred through phylogeographical studies based on the D-loop region (Guo et al., 2006;Wang et al., 2010;Ma et al., 2010). Wang et al. (2010) derived distinct phylogeographical differences between wild and domestic yaks based on the D-loop sequence data. We chose the mitochondrial D-loop data to discern the genetic identity of our wild yak sightings by comparing our samples with the reference data from the previous studies. However, mtDNA being maternally inherited is a haploid and nonrecombining locus and thus cannot present evidence of admixture or hybridization. But the presence of reference dataset (about wild versus domestic haplotypes) allows us to utilize the mitochondrial D-loop data to study phylogeographic patterns and domestication histories. Conservation Area (ACA). All the three study areas share an international border with the TAR of China. Elevation in the study areas ranged between 3,600 and 5,600 masl, and the vegetation is characterized by alpine grasslands and steppes interspersed with patches of shrubland (Miehe et al., 2016). Local communities, belonging to the Tibetan ethnic group, are mostly agro-pastoralists who graze domestic yaks, cattle Bos taurus, dzos/jhoppas (yak-cattle hybrids) Bos spp., horses Equus ferus coballus, goats Capra aegagrus hircus, and sheep Ovies aries in the alpine pastures during spring and summer seasons. They shift the livestock herds among pasturelands and bring the animals down to the villages during winter (Bauer, 2004).

| Noninvasive genetic sampling
We conducted line transects (Nichols & Karanth, 2005) and opportunistic sign surveys in the potential areas to look for any historical or present evidence of wild yaks. We noted geographical locations of all places (e.g., monasteries, local houses, mountain passes, and other natural features), with a GPS device (GPSMAP 64s), where we retrieved evidence (e.g., body parts such as horn, skull, pelt) of presumed wild yaks ( Figure 1). We collected genetic samples from these wild yak specimens (bone from horn and skull and hair with hair bulb from pelt) and noninvasive dung samples from live wild yaks. We also collected dung, hair, and bone samples from domestic yaks. We stored the dung samples in 2-ml cryo-vials containing DET buffer and the bone and hair samples in 15-ml sampling vials containing silica gels with a small layer of cotton atop the silica, before transferring them to the Intrepid Nepal Lab in Kathmandu, Nepal, where we performed the genetic analyses.
We collected 28 samples: six of presumed wild yak (three bone samples, two dung samples, and one hair sample), 21 of presumed domestic yak: (five dung samples, 14 hair samples, and two bone samples), and one unknown dung sample. We obtained genetic data from 25 samples, while two bone samples of suspected wild yaks and one hair sample of domestic yak failed possibly due to poor sample quality (see Table S1 for details). We treated the samples to be from wild yaks based on their morphology, that is, long hair and large skeleton (Wang et al., 2010), and the local history, that is, reports stating their origin from wild yaks.

| Genetic analysis
We extracted DNA off the genetic samples and generated par-  Table S2). This dataset represented all D-loop haplogroups of wild and domestic yaks generated by Wang et al. (2010) who inferred phylogeographical structure of the yaks across their range in QTP. We then compared the 25 sequences generated in this study (Table S1) with the reference sequence dataset by aligning them using MUSCLE (Edgar, 2004)  We also collected informal information on historical and current presence of wild yaks and traditional use of their body parts (see Acharya et al. (2015) for details).  We identified a total of 55 haplotypes from these 155 sequences, consisting of 48 haplotypes from published sequences of China and seven haplotypes newly identified from sequences of Nepal (Tables S1 and   S2). Of the 25 yak sequences from Nepal, 16 sequences belonged to eight previously known haplotypes and nine sequences belonged to seven new haplotypes (Figures 3 and 4, and Table S2). The sequences generated in this study have been deposited in NCBI GenBank under accessions MW048416 to MW048440. F I G U R E 1 Locations of yak samples collected in this study. Mustang lies within the Annapurna Conservation Area (ACA) and Dolpa within the Shey-Phoksundo National Park (SPNP); Humla lies outside the protected area system of Nepal. Numbers 1-25 are the serial numbers (SN) of the samples as they appear in Table S1. A zoomed-in view of the samples collected in Upper Humla is provided in the inset

| RE SULTS
The haplotype network showed that the domestic yak haplotypes of Nepal are mostly distributed among three major hap-  Figure 1) belonged to the same haplotype that was identical with a haplotype of wild yaks in north-western QTP. This haplotype associated with an endemic haplogroup of wild yak within lineage I (circled in Figure 3; haplotype "china21" in Figure 4). The haplotype F I G U R E 2 Photographs of wild yaks from Upper Humla. (a) The single wild yak seen in Upper Humla in July 2015. Its greyish-white muzzle, handle-bar horns, hump raised above the shoulder, long shaggy black-dark brown coat and thick tail are characteristic of wild yaks. Note the dung pile between its hind legs. Photo: Naresh Kusi. (b) A wild yak sighted with a domestic yak in Upper Humla in September 2020. The domestic yak lacks a greyish-white muzzle, has a nonuniform coat (with white patches on its forehead and shoulder), thinner and smaller horns, is smaller in size, and carries an inconspicuous hump. Photo: Bishnu Bahadur Lama  were found to be identical to a haplotype belonging to an endemic haplogroup of wild yaks sampled in north-western QTP by Wang et al. (2010). This particular haplotype is only five nucleotide mutations away from the next domestic animal, which is less than the variation observed among domestic haplotypes in haplogroups A and B (Figure 3). However, there are also many haplotypes that are unique to the wild yak gene pool which are at a relatively small genetic distance from the closest domestic haplotype. A broader and more comprehensive sampling and analysis involving a genomewide SNP-based assessment will probably provide clearer insights into the distinction between wild and domestic yaks. Interestingly the location where the haplotype "WYNepal_Humla" is found lies approximately 1,000 km south-west of north-western QTP where the "China21" wild haplotype was identified. Most of this region (between Upper Humla (Nepal) and north-western QTP) represent wilderness areas with no human settlements. So it is very unlikely that domestic yak herds can get transported between these points.
This further suggests that "WYNepal_Humla" is a wild haplotype.
Our data strongly suggests that the wild yak haplotype identified from Humla in this study is unique and rare among wild yaks and has a potential to get incorporated into domestic gene pool from hybridization. National and international experts agreed that the morphological features (greyish-white muzzle, long, and shaggy black-dark brown coat, thick tail, hump raised above the shoulder and handlebar horns) and behavior (they were very shy and ran away as soon as they saw humans) of the animals we sighted in 2014 are typical of wild yaks (see Acharya et al., 2015). A proper management intervention is urgent to formally protect the wild yak and its habitat in north-western Nepal.
The other two haplotypes of presumed (historical) wild yaks collected from Humla (hair from pelt, sample age: approx. 5 years) and Dolpa (bone from skull, sample age: approx. 10 years) are identical to contemporary haplotypes belonging to domestic yaks in hap- Phylogenetic tree of all domestic and wild yak D-loop haplotypes/haplogroups in western China and north-western Nepal (this study) constructed by Bayesian inference, rooted with Bison bison. The length of the alignment is 637 bp. The haplotypes from reference sequences are referred with prefix "china" while new haplotypes identified in this study are referred with prefix "nep." Branches with filled circles indicate haplotypes found only in wild yaks; filled diamond indicates haplotype found in unknown yak; open circles are haplotypes shared by domestic and wild yaks; branches without circles are haplotypes identified only in domestic yaks. The haplotypes representing 25 yak samples from this study are marked with an asterisk. The wild yak sighted in Humla belonged to haplotype "china21" which is a wild yak haplotype previously identified in north-western QTP in China. Support values at the nodes represent Bayesian posterior probabilities. Accession numbers are listed in Tables S1 and S2 their geographic locations. These results corroborate the previous findings of Guo et al. (2006) and Wang et al. (2010) (Figure 2b). This lack of reproductive barrier between wild and domestic yaks allows them to hybridize. Similar trends have been observed in Asian wild water buffaloes Bubalus arnee whose wild population face serious threats due to hybridization with domestic water buffaloes Bubalus bubalis (Flamand et al., 2003;Scherf, 2000).
Domestic yak bulls that mate with wild females can mediate genetic introgression into wild yak populations. From the perspective of genetic diversity, hybridization is potentially detrimental to wild yaks because frequent introgression will gradually change their genetic integrity and cause the wild gene pool to relinquish among the admixed ancestry in the longer term. However, the wild yak gene pool in Upper Humla, which consists of only a few individuals, is unlikely to get exposed to a severe introgression because the hybrid progeny (produced through mating of a wild male and domestic female) usually remains within the domestic herd that is a part of a bigger domestic population comprising thousands of animals. Importantly, herders do not prefer hybrid offspring as they are generally shyer and more difficult to handle (Buzzard & Berger, 2016). Extent of such hybridization footprints are best studied by using markers from mitochondrial as well as nuclear genomes.
However, we have utilized mitochondrial marker (that are maternally inherited) only in this study by referring to Guo et al. (2006)

| Importance of habitat protection for wild yak conservation
Almost 10  poaching and trade will support wild yak conservation in Nepal.

ACK N OWLED G M ENT
We thank the Department of National Parks and Wildlife Conservation