All sampling procedures in this study were carried out in accordance with the recommendations of “Guide for Care and Use of Agricultural Animals in Research and Teaching, 3rd edition” published by the Federation of Animal Science Societies (ISBN: 978-1-884706-11-0) and also “Guidelines for sampling and handling animals for purposes of Veterinary and Other researches” of Livestock sector of Mongolia, with full approval from Committee Ministry of Agriculture and Food of Mongolia (ref. A67/2010, in Mongolian). The protocol approved by Committee of Ministry of Agriculture and Food of Mongolia was used in this study. Animal owners had informed consent for inclusion of their animals in the study, and they involved sampling procedures.

Sampling was conducted across Mongolian territory in 2017 and 2018. In this study, a total of 158 unrelated genetic samples that consisted of 97 blood (2017), 10 ear skin, and 51 tail hair (2018) were gathered from 11 different populations (Figure 1A). We collected different samples during sampling years due to availability of field equipment. Of 11 populations, two populations were recognized as separate native strains (Zalaa Jinstiin White; MZG and Gobi-Gurwan Saikhan Black; MGG), while the other nine populations were believed to be Mongolian indigenous native goats (MNG) (Table 1, Figure 1A). Geographical distances between sampling sites ranged from ~65 km to ~1,450 km. When selecting animals from each investigated population, we involved owners to avoid sampling from related animals. Genomic DNA was extracted from punches of filter paper (FTA classic cards Whatman: Lot No.9767339) using a MagMax DNA Multi-Sample Ultra Kit, while genomic DNA was extracted from ear tissue and hair samples using a DNeasy Blood and Tissue Kit (QIAGEN, Hidden, Germany) following the recommended manufacturer’s instructions.

To study genetic diversity and phylogenetic analyses of MG, we amplified a 452 bp mtDNA fragment of the caprine partial HV1 region in D-loop (positions 15,678 to 16,130 on the complete goat mitochondrial sequence of reference NC_ 005044) from158 animals using a set of primers: gd-loop F (5′-CAC AAA CTT CCC ACT CCA CA -3′) and gd-loop R (5′-AGC GTG TTT AAA ACG GTG GT -3′). We designed these two internal primers using the online tool Primer3Plus (online available at: http://www.bioinformatics.nl/cgi-bin/primer3plus/primer3plus.cgi) based on a reference sequence (accession number: NC_005044). Polymerase chain reaction (PCR) amplifications were performed on a total of 20 μL final volume that included 2.0 μL genomic DNA (20 ng/μL), 0.8 μL of each primer (10 pmoles/μL), 2.0 μL 10× reaction buffer, 1.6 μL dNTP, 0.2 μL prime Taq DNA polymerase, and 12.6 μL distilled water to adjust the final PCR product volume. The reaction conditions for amplification consisted of initial denaturation at 94°C for 3 min, followed by 35 cycles with denaturation for 30 s at 94°C, annealing at 55°C for 30 s, extension for 1 min at 72°C, and final extension at 72°C for 5 min. All PCRs ran with negative controls and were performed by a thermal cycler C1000 system (Bio-Rad laboratories Inc., Hercules, CA, USA). After amplification, approximately 5 μL of PCR product was run on 1.5% agarose gel at a constant voltage of 150 V for 30 min in 1×TAE buffer prior bands and size checked under UV light. A service provider company (GenoTech Corp. Deajeon, Korea) performed direct sequencing of purified PCR products using a Sanger DNA sequencing approach in an ABI Prime 3730 DNA sequence machine. All 158 novel sequences of Mongolian goats in this study were deposited in the Genbank public database under accession numbers MK165491–MK165648.

A total of 158 sequences representing 11 MG populations from this study were merged with the previously described 16 sequences (considered as separated populations) for total of 12 MG populations. To examine the phylogenetic status and the origin of MG, 174 sequences of Mongolian goats were initially edited with BioEdit v.7.2.5 [23] and partial sequences aligned positions using the ClustalW tool [24] together with an additional 329 reference sequences of 24 foreign breeds from different geographic regions in the Old World and three wild caprine species, as well as representative sequences of the previously defined six caprine haplogroups (A, B, C, D, F, and G) (accession numbers of all reference sequences in Supplementary Table S1).

At the individual population level, as well as all MG together, DNA Sequence Polymorphism v5.10.01 (DnaSP [25]) was used to calculate the following diversity parameters of mtDNA: number of haplotypes (h), haplotype diversity (Hd), number of segregating sites (S), nucleotide diversity (Pi), and average number of nucleotide differences (K).

We constructed a neighbor-joining (NJ) tree (using 500 bootstraps to determine evolutionary distances computed based on the Tamura-Nei distance model; [26]) to examine possible haplogroups in MG populations and phylogenetic relationships among MG populations, as well as other foreign breeds using Molecular Evolutionary Genetic Analysis version.7.0.21 (MEGA [27]). The genetic distance (Tamura-Nei) between 12 MG populations were also assessed using MEGA. The Roehl network data (*.rdf) files for network analysis were generated by DnaSP. Then, we drew the Median-Joining (MJ) network using NETWORK v.5.0.0.1 [28] to trace genetic relationships among identified haplotypes within MGs and between other breeds from different regions.

Arlequin v3.5 [29] was used to estimate pairwise genetic distance (significant values were accepted when p<0.05), average number of pairwise differences between and within populations, and analysis of molecular variance (AMOVA) to assess the population genetic structure (with 1,000 permutations). We used, and explored in Arlequin v3.5, two main approaches: i) Fu’s (1997) neutrality tests of Fs and ii) mismatch distribution analysis [30], to assess signs of historical population expansions. Two tests of goodness of fit, sum of squared deviation (SSD) and Harpending’s raggedness index (r) [31], were used to evaluate the significance of deviations within observed patterns of mismatch from the simulated model of demographic expansion (p value accepted at <0.05).

Comparison of 174 sequences of the HVI region of goats from 12 different Mongolian goat populations showed high polymorphism. There were 109 polymorphic sites over 452 bp, while only two insertion/deletion (InDel) mutations were detected. Estimated diversity parameters for MG populations are given in Table 1. Among these polymorphic sites, we identified 20 singletons and 89 parsimony informative sites that defined 137 haplotypes (three of these haplotypes defined by two indels) with a relatively high Hd (overall; 0.997±0.001). Estimated values of Hd ranged from 0.933 (Erdenet; ERD) to the maximum value of 1.000 (MGG, MZG, and Ref), while nucleotide diversity (Pi) varied between 0.0199 (Lun Tsag; LTS) and 0.0384 (Altai; ALT) (overall, Pi = 0.0283) (Table 1). The mean K was 12.75, and the highest K-value was observed for the ALT population (17.34). In addition, an average base composition was A = 31.0%, T = 30.3%, G = 16.5%, and C = 22.2% in a 452 bp HVI region. Thus, the percentage of A+T (61.3%) nucleotide pairs was higher than the nucleotide pairs of C+G (38.7%).

Phylogenetic relationships of goat populations within Mongolian and between 36 breeds (MG plus foreign breeds) were initially investigated in this study based on the high polymorphic HVI region in mtDNA. Because previous studies [6,11] have described two origins for haplogroups of domestic goat (one for B, remaining one same for other five haplogroups), haplogroups in MG were identified through an un-rooted NJ tree using representative sequences of six haplogroups to discern possible origin for MG in this study. More in detail, the NJ tree clearly divided all MG sequences into four different haplogroups: A, B, C, and D (Figure 1B). The findings suggested two origins for the MG population. Among these, haplogroup A was highly predominant and found in all populations of MG (100% for ERD, LTS, Baganuur; BAG, and MZG). In contrast, relatively low numbers of animals belonged to haplogroups C (6/12 populations) and D (5/12), and only one for haplogroup B (1/12) (Table 1, Supplementary Table S2). Thus, the frequencies of haplogroups A, B, C, and D in MG populations were 90.8%, 0.57%, 3.4%, and 5.2%, respectively. Furthermore, ALT, Ikh Nart (INT), and Mankhan (MAN) populations were found to contain three haplogroups (A, C, and D), while there was no population with all four haplogroups (Table 1). To obtain further insight into the phylogenetic status of MG, a NJ tree was constructed once again for 12 MG populations and a total of 503 sequences that consisted of 36 populations (24 foreign and 12 MG), as well as three species of wild caprine using Tamura-Nei distance model (Figure 2, Supplementary Figure S1). Supplementary Table S3 presents Tamura-Nei genetic distance matrix (within/between) of MG. The results of NJ tree and genetic distance matrix showed relatively close relationships of MG populations, and indicated an absence of population-specific clusters in the NJ tree (Supplementary Figure S1). In addition, goats from China, Korea, and Central Asian (Kazakh and Kyrgyz), and some from Arabia had relatively close relationships compared to MG, while goats from Europe, Africa, and South Asia (except some Indian goats) clustered separately (Figure 2). Among the MG populations, Tamura-Nei distance within population reached a maximum estimated value in ALT (0.0402), while LTS accounted for the lowest value (0.0204). The genetic distance values between populations varied between 0.0389 (LOG versus Ref.) and 0.0207 (LTS vs BGN) (mean value was 0.0302) (further information see Supplementary Table S3).

To further discern relationships among identified haplotypes in MG populations and in the total dataset (503 sequences), a MJ network analysis was constructed separately. Results obtained from this analysis also supported four well-resolved haplogroups of MG (Figure 3). The topology of the MJ network illustrated relationships of 355 haplotypes of 503 sequences. As mentioned above, 137 of 355 haplotypes belonged to MG, and the majority of these 137 haplotypes clustered into haplogroup A (160 total haplotypes), with only one haplotype for haplogroup B, six haplotypes in haplogroup C, and seven haplotypes in haplogroup D (Figures 1B, 3). Among MG, 121 of 137 haplotypes were unique for a single population, whereas only 16 (11.6%) haplotypes were shared between two or three MG populations, and H19 and H60 were observed as the most frequent haplotypes with four animals (Figure 3, Supplementary Tables S2, S4). In addition, among 355 haplotypes, the three haplotypes that clustered into haplogroup A (H60, H88, and H188) were shared by Mongolian goats and goats from China (seven goats) and Korea (one goat) (Figure 3). Interestingly, Figure 3 also revealed lack of haplotypes that shared between multiple populations or breeds from different regions.

Additionally, we calculated Wright’s F-statistic (pairwise Fst) and average number of pairwise differences between (PiXY) and within (PiX) population values among 12 MG populations, while we also calculated Fst values for comparison among all 36 populations. Estimates of values of genetic differences are indicated in Tables 2 and Supplementary Table S5, respectively. The pairwise Fst values between 12 MG populations revealed relatively low genetic differentiation and ranged between −0.065 and 0.150 (mean Fst = 0.0168). The MGG population was closely related genetically to Nomgon (NOM; −0.065), followed by the distance between NOM and Ref (−0.0543), while comparison of LOG and BGN populations had a greater distance value (0.150) (Table 2). There were no significant relationships between genetic (Fst) and geographical (km) distances among MG populations (R² = 0.04, d.f = 53, p>0.05) (Figure 1C). Table 2 also presents the average number of pairwise differences between and within MG populations, and mean values were relatively similar, 13.8 and 13.6, respectively. The maximum PiXY value was obtained from a comparison between LOG versus Ref (18.08), and the minimum value was observed between LTS and BGN (9.326). Meanwhile, PiX values ranged from 18.44 (ALT) to 9.21 (LTS). The low genetic differentiation of MG populations was also evidenced by AMOVA analysis in this study. AMOVA revealed poor geographical structure among populations, 98% of total variation was within populations and only <2% was due to differences among populations (Table 3).

In addition, pairwise Fst values from comparisons of 36 breeds from MG and foreign countries and the majority of all possible pairwise Fst comparisons (234 of 288) indicated that MG populations were significantly (p<0.05) different from foreign populations, in particularly those from Africa, Europe, and some of southern Asia (Iran and Loas) (Supplementary Table S5). However, some populations were closely related to MG populations, in particularly a negative Fst value was observed between Kazakh and Mongolian goats (−0.190), followed by low positive values 0.026 for Indian and MGs, and goats from Arabian countries (including the domestication center).

We initially studied past demographic dynamics of MG based on two standard approaches of mismatch distribution (number of pairwise differences) and the neutrality test of Fu’s [32] Fs statistic. Because of the small sample size at the individual breed level (<25), these statistics were performed for all datasets of MG and at the haplogroup level (but, haplogroups B, C, and D were precluded due to small sample sizes). All MG and haplogroup A displayed multimodal and unimodal curves that mismatched distributions, respectively (Figure 4). In particularly, all MG showed two major curves at 11 and 23 pairwise differences with one minor curve at 39 pairwise differences. These results indicate there were at least two population expansion events that occurred in the past. Fu’s [32] Fs statistic provided more evidence for past population expansion of MG; this value was significantly negative for all MG and haplogroup A (p<0.001) (Figure 4). In addition, SSD and r values were not significant for all MG goats (p>0.05; 0.004 and 0.003, respectively), indicating that our total dataset (174 MG) showed relatively good fit to the population expansion model [31].

To date, no study has explored genetic diversity and the origin of MGs in depth. This study conducted the first assessments of genetic population structure, genetic diversity, phylogenetic status, and the origin of Mongolian goats using the most polymorphic HVI region in mtDNA. This highly polymorphic region (indeed, the complete D-loop) has been routinely used in previous studies that examined genetic diversity, phylogenetic relationships, and especially the maternal origin of most domesticated animals, including goats [14,33–37].

Similar to previous results [17,38,39], goats from Mongolia investigated in this study showed high mtDNA diversity. The assessed overall Hd and Pi of MG were 0.997 and 0.0283, respectively (Table 1). Our results also indicated that Hd in Mongolian goats were slightly higher (not significant) than that found in goat breeds from Africa, Egypt, and south Asia, but significantly higher than in Chinese and European breeds (Supplementary Figure S2A) [17,38,39]. Whereas, Pi values obtained from this study were significantly higher than that found in other regions including domestication center (except Chinese and southern Asian breeds) (Supplementary Figure S2B) [17,38,39]. The high level of mtDNA diversity in the world’s goat populations was previously explained by multiple maternal effects of wild ancestors [5,6]. Similarly, the high genetic diversity in MG likely resulted from the presence of an abundant haplogroup and crossbred/mixed with foreign populations from the Union of Soviet Socialist Republic (USSR) (Mountain Altai, Don and Unjuul) to improve the quality of their fibers [20].

As the origin of MG was only initially investigated in this study, these results are far from clear. Mongolians traditionally believed that their goats descended directly from wild Siberian ibex, a species widely distributed in Mongolia [19]. Yet, genetic and archeological studies widely agree that domestic goats in Mongolia originated from C. aegagrus in Fertile Crescent region of the NE [1,2,6,7]. Meanwhile, Lin et al [22] emphasized that haplogroup B or at least sub-group B1 originated in some region of Southeast Asia. Therefore, all animals that clustered into haplogroup B belong to east or south Asian countries, in particularly sub-group B1 accounted for only Chinese and Mongolian goats in this study (Figures 1B, 3).

The results of our phylogenetic analysis revealed four distinct maternal haplogroups (A, B, C, and D) in MG populations (Figure 1B). Consistent with previous studies [6,14,15,22,40], haplogroup A was predominant (90.8%) and found in all MG populations. Haplogroup B was found mostly in Asia, but only one individual from the Ref population was found in this haplogroup in this study (Figure 3). Alternatively, haplogroups C and D were found from both Asia and Europe, and goats from G and F had very limited distributions in Sicily and the NE region, respectively [6]. Some previous studies [6,14] reported that three haplogroups (A, B, and C) occurred in Mongolia based on a very limited sample size. More recently, Lin et al [22] included some Mongolian goats (population’s locality unknown) in their diversity analysis. They reported the four haplogroups that we found in this study. To obtain further insight into the origin of MG, NJ tree, MJ network, and estimation of pairwise Fst values were performed based on a combination of goats from Mongolia and those from different regions, including the domestication center (Figures 2, 3). The close relationships of Mongolian and Chinese (from Huanghuai and Leizhou) goats were obtained from MJ network analysis; they shared three haplotypes (H60, H88, H118, Figure 3, Supplementary Table S4) and occurred in NJ trees’ clusters. Interestingly, goat populations (n = 1) from Kazakhstan were very closely related genetically to MG (mean, Fst = −0.016), and followed by Fst values (mean = 0.103) obtained from comparisons of goat populations (n = 5) from some Arabian countries (including the domestication center; i.e. Turkish goat and Baladi breed from Jordan) and MG populations (n = 12), this value was lower than that found in comparisons of Chinese (n = 5) versus MG populations (0.128), as well as southern Asia (n = 4, mean Fst = 0.275) (Supplementary Table S5). These close relationships among goats from different regions likely resulted from a weak phylogeographic structure of domestic goats, which is likely explained by a high rate of gene flow among populations, even between continents [38].

According to findings from this and previous studies [6,14, 22,39], goats from the four haplogroups in Mongolia arrived from two possible domestication centers through ancient human migration or commercial trade routes. We hypothesize that the first domestication center was the NE region for MG from haplogroups A, C, and D whose ancient migration routes could be explained by a migration that started in the NE region and passed through central Asia (i.e. Kyrgyzstan and Kazakhstan, Fst values versus MG was 0.106 and −0.016, respectively), prior to arriving in Mongolia. This migration may have involved ancient sheep population, because since ancient times Mongolian nomads reared sheep and goat together.

Second, a domestication center was located close to Mongolia somewhere in Southeast Asia [22]. We believe that a smaller group of ancient goats from haplogroup B (B1 subgroup) arrived in Mongolia from this center. It is difficult to determine the date of the first arrival of domestic goats into Mongolia. Nevertheless, archeological findings based on petroglyphs in Mongolia demonstrate that the first goat populations arrived at least 3,000 to 3,500 years ago [6,41]. If the first goats arrived with sheep in Mongolia, this dated could be pushed back to 5,000 to 7,000 years ago [10].

AMOVA analysis of all MG dataset showed very weak phylogeographic structure among 12 populations, and more than 98% of the total variation occurred within populations. This weak structure was also evidenced by limited genetic differentiations (mean pairwise Fst = 0.0168). Our results were consistent with previous studies [6,38,42]. The traditional nomadic culture of Mongolia may have contributed to this weak geographical structure and limited genetic differentiations among MG populations. Under nomadic pasturing, there is a high rate of animal exchanges of the five main types of livestock animals in Mongolia (sheep, horse, cattle, camel, and goat) and often long-distance movements even cross the country.

Mismatch distribution and the neutrality test of Fu’s [33] Fs statistics are recognized as suitable approaches to examine past population expansion and were routinely used in previous studies [6,36,38,39]. Mismatch distribution is a frequency graph of pairwise differences between alleles, and it has two common patterns, including multimodal (populations at demographic equilibrium) and unimodal (populations passed through a recent demographic expansion) [29,43]. Fu’s Fs statistic is based on the probability of having a number of alleles equal to or greater than the observed number in a sample drawn from a stationary population [32]. Both approaches revealed at least two population expansions of Mongolian goats in the past: i.e. mismatch distribution pattern was multimodal for all MG while unimodal for haplogroup A, and highly significant and negative Fu’s Fs values were observed for all MG (−24.06) and haplogroup A (−24.39). Relatively similar patterns of mismatch distribution and significant and negative Fu’s Fs values were obtained in previous studies [37, 39,42]. More recently, the number of goats in Mongolia has increased fivefold within the last 30 years, more specifically the goat population grew from 5.1 million in 1990 to 27.3 million by 2017 [18].