Origin and Genetic Diversity of Nepalese Indigenous Goats (Capra hircus)

Nepal has a sizeable indigenous goat population inhabiting from the high altitude of the Himalayas to lowland approaching the sea level. The Terai goat is well-adapted to the topography and climate of the Terai region (<610 masl). Khari goat, a prolific meat breed, is well-adapted to the mid hill region (610-4877 masl) throughout the country from East to West. Chyangra and Sinhal are adapted to the high hills (>4877 masl). None of the studies conducted till date has been thorough in estimating the origin of these goats.
       
Recent studies based on mitochondrial DNA (mtDNA) sequences (Wu et al., 2009) succeeded in the identification of several maternal haplogroups suggesting that the goat was domesticated from different populations of the wild bezoar goat (Capra aegagrus). A comprehensive analysis of domestic goat revealed a total of seven different monophyletic mtDNA haplogroups named as A, B, C, D, E, F and G (Naderi et al., 2007). Haplogroup A, which encompasses more than 90% of individuals, is the most predominant and widely distributed across all continents and corresponds to the initial domestication event. It is followed by haplogroup C and B (Luikart et al., 2001). Haplogroup D is considered to be as old as haplogroup A or even older (Sultana et al., 2003 Chyangra and Sinhal goats might have genetic introgression from Tibet, China as Nepalese goat and sheep flock of these areas share common pasture with the Tibetan flocks as they follow migratory system to escape from the extreme cold during winter (Gorkhali et al., 2006). The goat breeds in the mid hills (Khari) and lowlands (Terai) could have been influenced by the Indian goat breeds as these two countries share free border and easier accessibility. Apart from a single molecular level study on Khari goat which revealed three distinct strains within the population (Kunwar et al., 2000), very little attention has been paid to the genetic diversity of these indigenous goat resources. The main objective of this study was to examine genetic diversity and phylogeography of the six different Nepalese goat populations based on the analysis of mtDNA D-loop hypervariable (HVI) region. A comprehensive analysis of collected molecular and geographical data was conducted to address the issues regarding genetic diversity and distribution of mtDNA haplogroups followed by comparisons with domestic and wild goat populations in the region.

Population sampling
 
Blood samples were taken from 93 individuals belonging to (a) three populations of Khari breed based on Eastern, Western and Central hill regions, (b) three other distinct breeds- Chyangra, Sinhal and Terai, selected to have unrelated parents and grandparents on the basis of the information provided by the owners and cross-checked with the neighbours. The details of breeds, geographical regions and sample sizes are given in Table 1, Fig 1 and 2.
 

 

 

 
DNA amplification and sequencing
 
Total genomic DNA was extracted from the whole blood following the instructions of the manufacturer (TIANamp Blood DNA kit, Tiangen Biotech, Beijing, China).
       
625 bp belonging to the Mitochondrial DNA HVI region of goats were PCR amplified using the primers: forward 5'- CAT TAC ACC GCT CGC CTA C- 3' and reverse 5'- GGG CTG ATT AGT CAT TAG T- 3' following the procedure of Wu et al., (2009), using 50 µl reaction mixtures comprising of  primer (1µl of a 10 umol/L solution), dNTPs (4 µl of a 2.5 mmol/L solution), 10X buffer (5 µl) and of Taq DNA polymerase (5U/µl, 1.5 µl) (Tiangen Biotech, Beijing, China) was taken for each primer. The PCR process comprised of an initial denaturing step at 95oC for 5 minutes followed by 35 amplification cycles (94^oC for 30s, 56.2^oC for 30s and 72^oC for 30s) and a final extension at 72^oC for 10 min in a Programmable Thermal Controller. The mtDNA HVI region sequences generated in this study were deposited in the GeneBank under accession numbers KM198594- KM198686.
 
Analysis of sequence data
 
The mtDNA HVI region from the experimental samples were edited using Chromas version 2.23. The published mtDNA HVI region sequences from domestic goats from Asian region as well as wild goats were included in the analysis. Sequences were aligned by the Cluster W method and neighbor-joining (NJ) tree was constructed for these mtDNA HVI region sequences with all reference sequences based on the Kimura 2-parameter model of MEGA version 5.2.2 (Tamura et al., 2011). The robustness of internal branches was estimated based on 1000 bootstrap replications. Twenty two goat mtDNA control region recommended by Naderi et al., (2007) were also included in the analysis, to facilitate the recognition of haplogroup status of each individual. Haplotype diversity (h) and nucleotide diversity (p) (Naruya and Nei, 1987) for each goat breed/population were also estimated by using DnaSP 4.10 (Rozas et al., 2003). Median Joining (MJ) network (Bandelt et al., 1999) was constructed by using Network 4.2 (http://www.fluxusengineering.com/sharenet.htm). In order to examine whether there are genetic differences among different populations, analyses of molecular variance (AMOVA) (Excoffier et al., 2005) was performed using Arlequin version 3.11 (http://anthropologie.unige.ch/arlequin/).

Haplogroup classification and phylogenetic analysis
 
The mtDNA control region HVI sequences from all the 93 samples were found to be highly polymorphic with 102 variable sites across the 625 bp of the alignment. Out of the total variable sites, 100 variants were transitions and only 2 variants were transversions. There were no insertions/deletions observed. 
       
Sixty-six different haplotypes were identified out of the 93 mtDNA HVI region sequences revealing high haplotype diversity. Three haplotypes were found to be shared between breeds: one between Khari Bandipur and Khari Salyan and two between Khari Bandipur and Sinhal. Individuals of different breeds from different geographical regions such as Khari Bandipur and Sinhal also share haplotypes. This median joining network depicts clearly that breeds from different geographic regions did not cluster together and then separated from other regions. This result indicates that there is no correspondence between the geographic regions of origin and relationships among breeds.
       
Four haplogroups (A, B, C and D) were identified, when the NJ tree of 93 mtDNA HVI region sequences together with the 22 reference sequences belonging to the known six mitochondrial haplogroups (Naderi et al., 2007) (Fig 3) was constructed. Among them, haplogroups A was the predominant (74%) group followed by B (20%), both of which were present across all six goat populations (Table 2). There was no specific haplogroup distribution pattern in breeds/populations or among the different ecological regions (Table 1 and Fig 2). Only one breed (Chyangra goat) contained all four haplogroups (A-D). Khari Bandipur goats (NGKB) contained haplogroups A, B and D. Other breeds/populations contained only two haplogroups (A and B). When the ecological distribution was considered, high hill regions contained all the four haplogroups, mid hill regions contained three haplogroups (A, B and D) and lowlands contained only two haplogroups (A and B). Specifically for haplogroup B, all the breeds had sub-haplogroup B1except for Chyangra (one individual), which had sub-haplogroup B2 (Table 1). Haplogroups B2, D and C showed southward gene flow from Northwestern China down to alpine region and hills of Nepal and never reached to the lowlands, whereas B1 is found in all over the country except to the alpine region.
 

 

       
The NJ tree drawn against wild goat populations revealed that none of breeds were close to any of the existing wild goat populations.
 
Genetic diversity
 
Number of haplotypes found in each breed/ population ranged from 6 to 17 depending on the difference in number of samples and the diversity ranging from 0.86 to 0.99. Both indices revealed a relatively low mtDNA variation in Khari Ilam (NGKI) and Sinhal (NGS) while Chyangra (NGCh) goats were highly diversified. The genetic variations were found to be distributed at 92.43% and 7.54% within breeds and among breeds, respectively as per AMOVA. According to the ecological distribution, the samples from the high hill area had the highest haplotype diversity (0.98±0.02) and nucleotide diversity (0.04±0.01).
       
Genetic diversity for the four haplogroups were also estimated (Table 2). Among them, haplogroups A had the highest and haplogroup D had the lowest nucleotide diversity. As haplogroup C had only two individuals and belonged to the same haplotype, diversity could not be calculated for this group.
 
Phylogeographic analysis/ Origin of Nepalese goat
 
It was seen that both Khari Salyan and Terai goat sequences shared with Indian goats, while Khari Salyan shared with Pakistani and Chinese goats, indicating same origin or gene flow between goat populations (Table 3).
 

       
As we found the strong gene flow among different breeds of neighboring countries and considered China as the origin of Haplogroup B (Liu et al., 2006), we further presented MJ network to locate an interesting route of gene flow in the region (Fig 5). The network contained two sub-clades with star-like phylogeny, in which two high frequency haplotypes H_6 (78 individuals, 38.6% of the total B) and H_19 (28 individuals, 13.86% of the total B) were located in the center. In between H_6 and H_19, H_21 (11 individuals) with one mutation difference from H_19, was shared with one individual from Mongolia and with 10 individuals from China. In addition, one Nepalese goat was only one mutation away from H_21. The haplotype H_19 was shared with individuals only from China whereas haplotype H_6 was shared among individuals from China (65 individuals), Laos (2 individuals), Malaysia (6 individuals), Pakistan (1 individual), India (3 individuals) and Nepal (1 individual). Along with haplotype H_38, there could be four different domestication events in Asian population HapG B.
 

 

       
Furthermore, graphical distribution of haplogroups showed maximum gene flow under haplogroups level in goats among different Asian countries (Fig 6). The green arrow shows that sub-haplogroup B2 found in the Tibetan goats influenced Chyangra goat, found in the alpine region and did not cross the boundary of Himalayas; whereas in case of sub-haplogroup B1, gene flow (maroon arrow) revealed the southward and south-west direction from Southern China, a speculated domestication centre for HapG B, down under the trans-himalayan region to Nepal.
 

       
mtDNA sequences of the studied samples showed notable haplotypes diversity: sixty-six haplotypes were identified in the studied samples. The wide variability might have been due to the multiple maternal origin of the population, as reported by other researchers (Liu et al., 2006) studying on other goat populations. The result is corroborated by the NJ tree drawn against wild goat species which revealed that none of the Nepalese breeds are close to any of the expected wild ancestors. Multiple maternal origins in goats might be due to their unique adaptation characteristics in new environment which in turn cause gene flow easily among different populations of goats (When et al., 2004). Some of these haplotypes were shared by individuals of different breeds from completely different geographical regions suggesting that there is no association between the geographical regions of origin and the breeds. The low genetic variation appeared among the breeds/ population (7.54%) and the high rate of genetic variation within breeds/populations (92.43%) also reinforced the above statement. Furthermore, this view is supported by the network profile which showed that populations from different geographic regions in Nepal are intermixed.
       
In the studied samples, neither haplogroup specific distribution pattern in breeds/populations nor among the different ecological regions was observed. However, the overall frequency distribution of the four haplogroups in all the Nepalese goat breeds studied was consistent with the world scenario reported in previous studies (Naderi et al., 2007), in which haplogroups A (74.14%, 69/93) and B (20.41%, 19/93) were the main components of the Nepalese goats. Even though the haplogroup A is the largest group, it does not indicate any clear geographic pattern in the network profile which is complicated by the abundance of sequence homoplasmy and is consistent with the low bootstrap support in the phylogenetic tree.
       
While taking into account only haplogroup B, as the Eastern Asia is claimed to be centre of emergence of the Asian goat (Wu et al., 2009), a MJ network with the regional goat mtDNA sequences illustrates an interesting route of gene flow and shows at least two different origins for development of Nepalese breeds. On one hand, it shows that China from the northern border might be the origin at least for Chyangra goat in alpine region while on the other, the southern Nepalese goat (Khari from mid hills and Terai from lowland) shared haplotypes with Indian, Bhutan and Pakistan goat breeds. These results suggest a strong gene flow among goat populations promoted by the traditional seasonal pastoralism, annual long distance migrations and the ancient trade between Tibet and India via Nepal accounts for the pattern discerned in the regional goat gene pool. This prospective is again reinforced by the molecular evidence of some indigenous people of Nepal who show genetic relationship with Tibet and also have the common ancestry for the Tibeto-Burman language (Gayden et al., 2009).
       
Sinhal goat, which does not share with any of Nepalese and regional goats, might be different from other indigenous breeds and considered as native for Nepal. However, the finding of AMOVA reflects low differentiation in the maternal lineage among the different breeds. Interestingly, beside other haplogroups, haplogroup D is found to be trait-specific as this haplogroup is observed in all the hairy goats from Pakistan, Nepal, Tibet, Xinjiang, Liaoning and Inner Mongolia irrespective of geographically distant.
       
All the above evidences and the results of the current study support that the genetic diversity and structure in mtDNA genome among indigenous Nepalese goats were shaped not only by the intensive and continuous gene flow among the goats distributed in middle and lowland in Nepal and geographical vicinity in India but also by the exchanges between goats found in high hill of Nepal (e.g. the B2 haplotype present in Chyangra goats) and Tibetan goats in China.

The results indicated extensive and high mtDNA haplotype diversity among Nepalese goat breeds/populations. The complex mtDNA diversity and structure identified among indigenous Nepalese goats can be explained by the gene flow through ancient trading and current ‘free’ movement of goats from/to the geographic vicinities in India and China.

The authors declared no conflicts of interest.

This work was funded by “IFAR fellowship” of World Bank through International Livestock Research Institute (ILRI) and logistic facilities were provided by CAAS-ILRI Joint Laboratory on Livestock and Forage Genetic Resources, Institute of Animal Science, Chinese Academy of Agricultural Sciences, Beijing, China.