Revisiting genetic structure of Wild Buffaloes Bubalus arnee Kerr, 1792 (Mammalia: Artiodactyla: Bovidae) in Koshi Tappu Wildlife Reserve, Nepal: an assessment for translocation programs

Koshi Tappu Wildlife Reserve (KTWR) has the last remaining Nepalese population of the Endangered Asiatic Wild Buffalo (Bubalus arnee Kerr, 1792). Individual animals protected inside KTWR may be of purely wild, domestic or hybrid origin, and the wild population is under potential threat due to habitat loss and genetic introgression from feral backcrosses. Identification of genetically pure wild individuals is important for identifying animals for translocation to other areas within their former range. In this study we have sequenced a highly variable 422bp region of the Cytochrome b gene of 36 animals, and added 61 published sequences of both River and Swamp Buffalo from Italy and some southern Asian countries including India. The haplotype diversities ranged from 0.286-0.589 with slightly higher diversities in domesticated individuals. The AMOVA analysis revealed that 97.217% of the genetic variation was contained within groups and 2.782% occurred among groups. An overall fixation index (FST) was found to be 0.02782 (p>0.05). Phylogenetic relationships derived through a reduced median network and maximum parsimony analyses reconfirmed the ancestral nature of the Wild Water Buffalo. Moreover, this study has reviewed recent achievements of molecular research in wild buffalo, assessed the technical capacities of research institutes in Nepal to conduct molecular research required for identifying pure wild individual in KTWR and more importantly initiated DNA bank and DNA sequence library of buffalos, which will enable an international collaboration for advanced molecular research in the future.


INTRODUCTION
The domesticated Water Buffalo Bubalus bubalis Kerr, 1792 is one of the most important dairy and draft animals in southern Asia. Buffaloes are broadly categorized into two general 'breeds' sometimes characterized as subspecies: river (Bubalus bubalis bubalis) and swamp buffalo (Bubalis bubalis carabanesis). Despite having distinct morphological and behavioural traits and different karyotypes between these two categories of buffalo, they interbreed easily and produce progeny with intermediate chromosomes (Mishra et al. 2015). Swamp and river buffalo have different purposes and are found in different geographical areas. The swamp type, with wide-spreading horns and some white markings, is more similar to Wild Buffalo Bubalus arnee (Kumar et al. 2007;Mishra et al. 2015).
The Government of Nepal (GoN) established Koshi Tappu Wildlife Reserve (KTWR), an IUCN Category IV protected area (Heinen 1995), in 1976 primarily to conserve the last Nepalese population of Wild Buffalo (Heinen 1993). Wild Buffalo cohabit the reserve with highly backcrossed feral buffalo thought to have been released in the area in the 1950s (Dahmer 1978). Wild Water Buffalo Bubalus arnee (Kerr, 1792) is considered Endangered globally (Kaul et al. 2019), with isolated populations in KTWR and selected areas of Bhutan, India, Thailand and possibly Myanmar and Cambodia (Groves 1996;Heinen & Srikosamatara 2003;Choudhury & Barker 2014). In 2016, 433 Wild Buffaloes were counted in KTWR (2016).
Despite the possibility of interbreeding between wild and domesticated buffalo, it is important to assign extant individuals to wild and other types where possible to broaden our understanding of the genetic structure of different types and to maintain genetic fingerprinting of wild breeds for their conservation. The introgression from domestic to wild population is female-mediated, therefore mitochondrial DNA (mtDNA) sequencing is likely to be helpful in the identification of group-specific mitotypes. The presence of wild-specific or domesticspecific haplotypes in either group would allow us to identify hybrids (Lau et al. 1998;Flamand et al. 2003). Furthermore, mtDNA sequence variations have been widely applied in mammals to study inter-and intraspecies phylogenetic relationships (Kikkawa et al. 1997;Lau et al. 1998;Conroy & Cook 1999;Kuwayama & Ozawa 2000;Kumar et al. 2007).
Conservation decisions on translocation should be based on putative wild, feral and domestic genetic assignments reliably performed through standard and widely accepted techniques. Therefore, selecting individuals for translocation programs, identification of wild individuals through detailed molecular study of the buffalo population protected in the reserve is a high priority for the Nepal government. In addition, understanding the genetic makeup of Wild Buffalo could be used as the basis for genetic improvement of domestic stock. National capacity building to conduct advanced molecular studies should be initiated from collecting blood and faecal samples, creating a DNA reference library and carrying out genetic research on various aspects such as population genetics, breeding behaviours among different buffalo types, disease dynamics, and food habits of buffalo population in the reserve. We present results of DNA sequence variation in the partial but variable cytochrome b gene among purely wild, feral and domesticated individuals and future prospects for advancing genetic research on Wild Buffaloes inhabiting KTWR in eastern Nepal.

Identification of organizations
Nepal Government, Ministry of Forest and Environment is working together with local communities and national and international conservation partners to protect wild animals and restore their natural habitats. With the aim of establishing viable populations of endangered species in different areas and safeguarding them from poaching and natural calamities such as flood, fire and epidemics, the Ministry of Forests and Environment, Department of National Parks and Wildlife Reserve (DNPWC) has been regularly translocating species to their original natural habitats, viz.: Onehorned Rhinoceros and Wild Water Buffalo. In 2016, 18 individuals of Wild Water Buffaloes were translocated from Koshi Tappu Wildlife Reserve (KTWR) to Chitwan National Park. Translocation was carried out by a team of 60 people including three veterinarians and 12 wildlife technicians led by DNPWC with support from the World Wildlife Fund Nepal, the USAID-supported Hariyo Ban Program, National Trust for Nature Conservation, Biodiversity Conservation Centre (NTNC-BCC), and the Zoological Society of London (Nepal). A series of consistent identification criteria that are based on phenotypic and behavioural characteristics were used to choose Wild Buffaloes from different herds for translocation. Given the availability of highly polymorphic genetic markers, the expert team strongly recommended to adopt genetic translocation as an effective and reliable management strategy. For this management strategy to be implemented effectively the team further emphasized the need for institutional strengthening and capacity building of national laboratories to conduct genetic research on buffalo before translocating them from the reserve in future. Wildlife veterinarians, biologists and research scholars from DNPWC and NTNC-BCC collected blood and faecal samples of buffaloes and initiated a genetic study in the wildlife laboratory of NTNC-BCC and molecular biotechnology laboratory of the Nepal Academy of Science and Technology.

Sampling and blood collection
A total of 42 blood and faecal samples (Table 1) were collected mainly from individuals residing in and around KTWR located in the Terai of southeastern Nepal (Figure 1). Animals were provisionally divided into three classes: domestic (D, n=11), hybrid (H, n=11), and wild (W, n=20) based on the consistent phenotypic and behavioural criteria (Dahmer 1978;Heinen & Paudel 2015) and location of herds sampling (domestic buffalo were sampled from the villages nearby KTWR while wild and feral were sampled using a location map of natal herds prepared by the reserve) and behavioural and anatomical phenotypic traits. All animals classified as domestic were river type buffalo with black bodies and curled horns (as in the Murrah breed of river buffalo), while those classified as wild had white chevrons, socks and tail tips, and larger, relatively straight, pale-coloured horns (Image 1) similar to swamp buffalos (Heinen 2002). Hybrid animals had intermediate phenotypes and may be first generation crosses or the result of various levels of backcrossing to either wild or domestic.

DNA Extraction, Polymerase Chain Reaction and Sequencing
Genomic DNA was extracted from blood using a Qiagen DNEasy Blood kit, according to the manufacturer's protocol. For the faecal samples, a Qiagen QIAMP DNA Stool Mini Kit was used following the manufacturer's instructions. Extracted DNA samples were stored at 4°C until they were used for molecular analyses. Aliquots of extracted DNA were used for PCR and sequencing. The mitochondrial partial Cytochrome b gene of 422bp was amplified using primer pairs (L14724: 5'-CGAAGCTTGATATGAAAAACCATCGTTG-3' and H15149: 5'-AAACTGCAGCCCCTCAGAATGATATTTGTCCTCA -3') (Kocher et al. 1989). PCR was carried out with 3µl template DNA, 15µl of Hot Start Taq 2X Mastermix (New England Biolab, UK), 1µl of each primer and 7µl of nuclease free water in a total reaction volume of 30µl using an ABI Veriti TM Thermal Cycler (Model no. 9902). Of

14946
the 42 samples, only 36 were amplified successfully and rest of the six didn't amplify even in multiple attempts due to poor DNA quality. The PCR conditions were an initial denaturation at 94°C for 10 minutes, followed by 35 cycles of denaturation at 94°C for 30s, annealing at 55°C for 30s, elongation at 72°C for 45s and a final extension at 72°C for 10 minutes. The PCR products were electrophoresed at 100 volts for 30 minutes in 1.5% agarose gels, viewed in Gel Doc (Syngene InGenius) after staining with Sybr safe and photographed. Amplified DNA fragments were purified using ExoSap-IT Express PCR Product Clean up (Affymetrix Inc., Santa Clara, CA, USA) following a cycle of 37°C for 15 minutes and 80°C for 1 minute in a thermo-cycler. High quality purified PCR amplicons were subjected to Dideoxy sequencing in a total volume of 10µl containing 1µl purified PCR product, 1µl Big Dye terminator sequencing mixture (V3.1) (BIGDYE Terminator Cycle Sequencing Kit, Applied Biosystems, Foster, CA, USA), 1.5µl sequencing buffer and 1.5µl primer (10µm). Sequencing was done at Nepal Academy of Science and Technology, Molecular Biotechnology Laboratory, in an ABI 3500XL automated DNA sequencer (Applied Bio-systems, Forster City, CA, USA). Sequencing for the majority samples was performed for one time but in the case sequence quality was low and/or polymorphism was observed then resequencing was done for confirmation. Successful sequences of 36 samples were deposited in GenBank (Accession no. MH718851-85)
All sequences were aligned and mitochondrial haplotypes were defined in DnaSP v5, and haplotype (h) and nucleotide diversity (π) were estimated using DnaSP v5 (Librado & Rozas 2009). Genetic differentiation within and between buffalo groups (wild, hybrid and domestic) was estimated by an analysis of molecular variance (AMOVA) using 10,000 permutations in Arlequin v3.5 (Excoffier & Lischer 2010). Similarly, pairwise genetic divergences between groups (F ST ) was also calculated and significances tested using 10,000 permutations in Arlequin v3.5.
Phylogenetic relationships among haplotypes were

Image 1. Translocation of Wild Buffalo and collection of blood samples from both domestic and Wild Buffalo: a-b-blood sample collection from domestic buffalo | c-Wild Buffalo Bubalus arnee | d-e-blood collection and translocation of Wild Buffalo from Koshi Tappu Wildlife
Reserve to Chitwan National Park, Nepal. Taxa  derived through a reduced median network v4.6 (Bandelt et al. 1999). To identify phylogenetic lineages, maximum parsimony (MP) analysis was performed in PAUP 4.0b10 (Swofford 2002) using the heuristic search option with 1,000 random additions and the tree bisectionreconnection (TBR) swapping and the MULTrees option on. Branch support was provided by a bootstrap analysis of 10,000 replicates of heuristic searches, with the MULTrees option on and TBR swapping off. Consistency indices (CI) and retention indices (RI) were obtained in PAUP. In addition, a neighbor-joining (NJ) tree (Saitou & Nei 1987) was produced using MEGA7 (Kumar et al. 2016).

Sequence Variation and Divergence
The lengths of the 36 partial cytochrome b sequences including primers at both ends for all the three types of buffalo were 486bp. When primers on either end were removed the sequence lengths were 422bp. The aligned matrix without primers contained two variable sites (Table 2), however, when our data set was compared with the accession data of Zhang et al. (2016), Kumar et al. (2007), and Kikkawa et al. (1997), base substitutions at 13 nucleotide positions (variable sites) were obtained and, among them, three nucleotide positions were specific for river buffalo and 10 positions were specific for swamp buffalo. The sequence divergence within buffalo of Nepal and India was 0.24 to 0.49 %; however, when compared with river buffalo sequence of Kikkawa et al. (1997), the divergence was slightly higher: 0.24 to 0.74 %. Sequence divergence within swamp buffalo was 0.24 to 0.98%. The sequence divergence between swamp and river buffalo was calculated to be 1.49 to 2.49%.

Haplotype Identification, Differentiations, and Phylogenetic relations
The 97 partial cytochrome b sequences, including 36 from this study, showed a total of 13 variable sites (Table 2, see Figure 2), which defined nine haplotypes. Nepalese buffalo had either haplotype H1, H2 or H3, but the three classes domestic (D), hybrid (H) and wild (W) were represented among both H1 and H2 (H1: D=23,W=18, H=17; H2: H=8, D=7, W=2 and H3: W=2). The most common haplotype (H2) was widely distributed among groups and represented by 66% of sequences (Nepal: D=23, W=18 and H=17; one of Kumar et al. (2007); and five of Kikkawa et al. (1997). The second most common haplotype (H1) was represented by 20% of the samples (H=8, D=7, and W=2) and one each from Kumar et al. (2007) and Kikkawa et al. (1997). Three sequences of Nepalese samples, two from this study and one each from Zhang et al. (2016) and Kumar et al. (2007) were restricted to the third haplotype (H3). Of the remaining six haplotypes, one was reported by Kumar et al. (2007;H4) while the remaining five (H5 to H9) were defined by Kikkawa et al. (1997). Haplotypes 4, 6 and 7 were found in river animals from other countries (i.e., not Nepal) and haplotypes 5, 8 and 9 were specific As above As above As above

14948
to swamp buffalo. Within groups, haplotype diversity was highest in hybrids followed by domesticated buffalo, and slightly lower in Wild Buffalo (Table 3). Overall haplotype and nucleotide diversities were 0.403 and 0.00105, respectively. Haplotypes were divided into two branches corresponding to river and swamp buffalo by six nucleotide mutations. River buffalo were found to be less diverse genetically than swamp buffalo. The AMOVA conducted for 78 Nepalese buffaloes under three groups showed highest variation (97.21%) partitioned within groups and very little variation (2.782%) partitioned among them with an overall fixation index F ST of 0.0278 (p > 0.05) ( Table 4). The pairwise genetic variations between groups also revealed non-significant low F ST scores (wild vs hybrid, F ST -0.055, p = 0.126; wild vs domestic, F ST -0.002, p = 0.416 and hybrid vs domestic, F ST -0.020, p = 0.546). The maximum parsimony analysis of the 97 sequences revealed four distinct clades (clades A through D) with moderate to high bootstrap values ( Figure 3). The three most parsimonious trees (CI = 1.00, RI = 1.00, length = 15 steps) were recovered. Clade A contained all H1 sequences plus one H4 and one H6, Clade B contained all 19 H2 sequences plus one H7, Clade C contained all four sequences of Haplotype H3 and Clade D, with a high bootstrap value (99%), included all swamp animals (Haplotypes H5, H8, H9). The strong separation of swamp and river buffalo is also shown by the median joining network, with six mutations separating H8 and H2. Neighbour-joining (NJ) phylogenetic analysis showed the same topology (Supporting information, Appendix 1).

Conservation status of Wild Buffalo population in Nepal
Nepal is home to a population of the Endangered Asiatic Wild Buffalo Bubalus arnee, the progenitor of domesticated water buffalo (Lei et al. 2007). Recent   (Khatri et al. 2013 ). There were only 63 individuals at the time of establishment of KTWR in 1976 (Dahmer 1978). The historic range of this species extended further west within Nepal, and at least as far as Chitwan National Park, however, Wild Buffalo have been restricted to KTWR for an estimated 60 years and they are under constant threat of extirpation from floods, habitat deterioration, hybridization, and the potential for diseases and parasites transmitted by domestic livestock (Heinen & Kandel 2006). Since Wild Buffalo have been eliminated from the greater part of their former  range, the Nepalese population is very important for the survival of the species globally (Beyers et al. 1995;Hedges 1995Hedges , 2001Choudhury & Barker 2014). Given the precarious existence of Wild Buffalo within KTWR, several wildlife conservationists have emphasized the need to translocate a sufficient number of individuals to sites within their indigenous range. Chitwan National Park had this species at least until the 1950s (Spillet & Tamang 1966;Aryal et al. 2011) and has extensive grassland areas, and much larger riverine habitats with sufficient upland areas that are not prone to flooding, compared to KTWR (Heinen & Paudel 2015). For these reasons 18 Wild Buffaloes were translocated to Chitwan National Park from KTWR recently, and more need to be moved in the near future (Shah et al. 2017;Kandel et al. 2018).

Identification of wild individuals in KTWR
Translocation of endangered species can restore species, protect populations from threats and reinstate the local ecosystem functions (Tarszisz et al. 2014).
Adequate morphological and genetic studies should be carried out to distinguish putative wild from feral backcrossed animals for translocation programs. During translocation enough purely wild individuals must be moved to assure genetic variation in the founding population. Although the recent selection of individuals for translocation to Chitwan was based on phenotypic and behavioural characteristics widely recommended by Dahmer (1978), many individuals of mixed wilddomestic ancestry may not be correctly distinguished from wild animals (Flamand et al. 2003).
On the basis of the partial cytochrome b sequences, we were able to define only three haplotypes in Nepalese buffalo. None of the haplotypes were specific to wild, domestic and hybrid types identified here. These haplotypes had been identified by Zhang et al. (2007Zhang et al. ( , 2016, Kumar et al. (2007), and Kikkawa et al. (1997) in their Nepalese, Indian, and wider samples (Pakistan, Bangladesh, Thailand and Sri Lanka). Although we have used the partial sequence of cytochrome b gene, the complete length (1,120bp) of this gene sequence reported from other studies (Kikkawa et al. 1997;Kumar et al. 2007) including 42 Nepalese samples from Zhang et al. (2016) did not show consistent distinctions between wild, hybrid, and domestic group of buffalos. Non-significant genetic variability observed between three groups of buffaloes and the highest, i.e., 97.21% total variations found within group in AMOVA analysis, clearly shows an evidence of gene flow between groups. In KTWR, most of the time wild, feral and domestic herds share grazing areas, where crossbreeding between groups occurs frequently. Moreover, low genetic variation between groups is also attributed to local the farmers' practice of crossbreeding domestic females with wild males (Heinen 2001).
NJ and MP analysis performed (results not shown) with the complete length (1,120bp) cytochrome b sequences of Zhang et al. (2016, 42 Nepalese sequences), Kumar et al. (2007, four river buffalo haplotypes), and Kikkawa et al. (1997, seven river and eight swamp haplotypes) provided essentially the same topography of the tree but an addition of one more haplotype represented by Genbank accession KR009945NP_D07 alone. Lau et al. (1998) using partial cytochrome b sequence (303bp) and D-loop (158bp) sequence suggested that Wild Asian Water Buffalo (Bubulus arnee) in Assam, Nepal and Indo-China is the possible ancestor of both river and swamp buffalo. The study by Tanaka et al. (1996) also supports this hypothesis. Nepal's Wild Buffalo show swamp type phenotypic characteristics but Zhang et al. (2011) found them genetically closer to river buffalo. Our study is consistent with Zhang et al. (2011) in that our MP, NJ and Network analyses of swamp buffalo showed distinct variation with a bootstrap value of almost 100% and six nucleotide differences between these groups. Similar results were obtained by Mishra et al. (2015) in upper Assamese and Chilika populations of domestic buffalo, which show phenotypic similarities to swamp buffalo but are also genetically closer to wild-type buffalo in the region. Our multilineage MP and NJ trees and several previous studies (Tanaka et al. 1996;Lau et al. 1998;Zhang et al. 2011Zhang et al. , 2016Mishra et al. 2015) inferred the ancestral nature of Wild Water Buffalo including the remnant population in Nepal. In this context, to determine genetically pure wild individuals in KTWR detailed and advanced genetic research is necessary.

National capacity building and close collaborations for genetic translocation
Phylogenetic analysis of partial cytochrome b region revealed overlapping clusters of wild, feral and domestic buffalo types residing in KTWR. The findings of this study along with previous genetic studies on water buffalo populations distributed globally including Wild Buffaloes of KTWR suggest that, single marker or partially overlapping markers are not sufficient to show the level of admixture and introgression of domestic in wild stock at the local level. Genome level assessments of source population offfers specific criteria and objective means for translocation of an appropriate  Heinen & Paudel (2015) should be represented from a source population of KTWR to the translocated population in the native area such as Chitwan, Bardiya or other appropriate sites (Heinen 2002) in Nepal. Founding genetic diversity of translocated population will be determined by the number of genetically variable wild individuals, the proportion of diverse pure stocks, those that contribute genetically to the next generation and number and frequency of polymorphic alleles that represent whole genomes of the source population. Translocated populations are mostly small in size therefore they are prone to loss of genetic diversity very rapidly through genetic drift (Frankham et al. 2012). Genetic assessment of source populations in advance of translocation (pretranslocation) helps to guide translocation plans and inform post translocation assessment or monitoring of genetic diversity in the founders (Groombridge et al. 2012). In addition to geneticists, active involvement of conservation biologists, wildlife experts, wildlife veterinarians, ecologists, physiologists and local people in the translocation program can ensure longer-term welfare and wellbeing of the re-introduced population.
Before embarking on a genetic translocation program for the buffalo of KTWR, Nepal should upscale its laboratory facilities, design population-based advanced genetic research and take the initiative to build a DNA bank of all possible individuals counted in 2016. The DNA bank, reference DNA sequences and genotype database are crucial for research and conservation efforts, to enhance our understanding of genetic effects of introgression, the study of dietary patterns on different buffalo types and assess the status of pathogens affecting the buffaloes with different genetic backgrounds. Using the same blood samples collected during this study we have reported the prevalence of malaria parasites for the first time in buffaloes of KTWR (Kandel et al. 2019). Given the lack of highly technical laboratories and trained manpower in Nepal and the urgency to identify wild individuals reliably, collaborative research with international universities, research institutes and conservation partners on advanced molecular studies are to be jointly conducted. In conclusion this research sets a baseline to develop well defined action-oriented strategies that guide pretranslocation genetic study of wild buffaloes in KTWR and their monitoring through post-translocation genetic studies. Key actions highlighted in this paper such as collaboration between partners, establishment of DNA bank of all extant individuals in KTWR, involve experts from different disciplines, upscale and strengthen present laboratories and build capacity of available human resources for genomic level data management are important steps to be taken by the Ministry of Forest and Environment, Department of National Parks and Wildlife Conservation and its national and international conservation partners for genetic translocation of Wild Buffalo including other threatened species of Nepal.