Goats (Capra hircus), domesticated in the Fertile Crescent approximately 10,000 years ago, are bred in most parts of the world today1,2,3. They are essential livestock for human life, providing milk, meat, hair, and other products. The genetic information of goats is crucial for understanding their genetic relationships, the origins of domestication, and the history of their propagation. Studies of goat genetics have primarily focused on mtDNA and the SRY gene on the Y chromosome4,5,6,7,8.

Previous studies mainly based on sequences of hypervariable mtDNA control region (481 bp) have identified six maternal haplogroups (A, B, C, D, F, and G)1,2,9,10,11. Haplogroup A is the most divergent globally, while the minor haplogroups C, D, F, and G are found sporadically at low frequencies across Asia, Europe, the Near East, and North Africa, respectively. The mtDNA haplogroups have been reported to show a weak phylogenetic structure due to extremely high mtDNA variation across regions and continent1,2,12,13. In contrast, haplogroup B is found at a high frequency in Southeast Asia. Given the unique phylogeographical structure of haplogroup B, it is postulated that the Baluchistan region of Pakistan may be the origin of haplogroup B1,11. However, the details of the origin and propagation routes of haplogroup B remain unclear.

As informative paternal genetic variations, five haplotypes (Y1AA, Y1AB, Y1B, Y2A, and Y2B) have been reported in previous studies using the Sex-determining region Y chromosome (SRY) 3’UTR sequence8,14,15,16,17. Y1A is further classified into Y1AA and Y1AB; however, Y1A in certain parts of Europe and Africa has not been divided into these haplotypes. The available data suggest a geographic contrast, with Y1AB predominating in northern China and Y1AA along with Y2B in the south. Y1B, Y2A, and Y2B are primarily observed in Europe, Africa, and East Asia, respectively17. The distinct distribution of SRY haplotypes across regions indicates that SRY haplotypes exhibit a strong phylogeographic structure, in contrast to the weak phylogenetic structure seen in the major mtDNA haplogroups17.

The Philippines and Indonesia, located in Southeast Asia, are archipelagic nations consisting of numerous islands. The Philippines, situated at the crossroads of past human migrations in the Asia-Pacific region, is believed to have never been connected to the Asian continent, even during the major marine regressions of the Pleistocene. In contrast, Indonesia was part of Sundaland, which was connected to present-day Mainland Southeast Asia (MSEA) by land during the Pleistocene and is located in the southeasternmost part of MSEA. In Island Southeast Asia (ISEA), there are historical records of trade and Islamic missionary work by Muslim merchants beginning in the 9th century, as well as colonial rule by Spain in the Philippines and by the Netherlands in Indonesia from the 16th century onward18,19,20.

It is generally thought that ISEA goats were propagated from the domestication center via MSEA and Maritime Silk Road. These small goats, called “Kambing Katjang”, are well known to be raised over a broad area of Southeast Asia17,21. Given geographical and historical background of ISEA, Island Southeast Asian goats may have distinct genetic influences compared to those in Mainland Southeast Asia. Supporting this, several previous studies on Southeast Asian goats have reported that the SRY haplotype Y2A, primarily observed in Africa, is found in ISEA but not in MSEA17,22,23. However, the genetic structure and introgression patterns of goats in ISEA remain unclear due to limited geographic coverage and sample numbers in previous studies.

Therefore, this study aimed to (1) provide a comprehensive understanding of the genetic structure of goats by sampling from a wide range of ISEA locations, (2) estimate the genetic differences among maternal and paternal lineages, and (3) infer the processes of propagation and gene introgression in ISEA goats, using sequences of mtDNA HV1 (hypervariable segment 1) and the SRY gene on the Y chromosome.

Results

We sequenced the mtDNA HV1 region (481 bp) of 176 Philippine goats and 72 Indonesian goats and analyzed these sequences in conjunction with previously published mtDNA sequences. Supplementary Figures S1 and S2 show the alignment of the D-loop sequence with the reference goat sequence (Accession number AF533441.1). The sequence analysis identified 59 and 45 variant sites, as well as 32 and 19 mitochondrial haplotypes, in Philippine goats and Indonesian goats, respectively (Supplementary Figures S1 and S2). Subsequently, neighbor-joining trees were constructed using 30 previously published Philippine goat sequences and 22 reference goat sequences2,22. Consequently, 206 Philippine goats were classified into haplogroups B (n = 152/206) and A (n = 54/206), while all 72 Indonesian goats belonged to haplogroup B (Supplementary Figure S3).

Figure 1 shows the median-joining networks of mtDNA haplotypes for Philippine and Indonesian goats with previously published Philippine and Indonesian goat sequences22,23. In Philippine goats, the median-joining network showed that the haplotypes were classified into haplogroups A and B, consistent with the neighbor-joining tree. In haplogroup A, two moderately frequent haplotypes (ISEA1 and ISEA2) were observed, differing by 19 base pairs. In contrast, haplogroup B showed a star-like appearance with haplotype ISEA12 being the most frequent (n = 95/152) as a center (Fig. 1a).

Fig. 1
figure 1

Median joining network of 481 bp in the mtDNA HV1 sequence. Separator lines on the Network indicate mutation positions. Blue and yellow colors are shown haplogroup A and B, respectively. The size of nodes is proportional to the number of goats present in the node. Each short bar on the edges represents a nucleotide mutation. (a) Philippine goats; (b) Indonesian goats.

Full size image

In Indonesian goats, the median-joining network showed that haplotype ISEA14 (n = 31/72) was the most frequent, followed by ISEA34 (n = 9/72) and ISEA33 (n = 7/72) (Fig. 1b). Compared with Philippine goat haplotypes, ISEA14 and ISEA13 were observed also in Philippine goats. However, the predominant haplotype ISEA12 found in Philippine goats was absent in Indonesian goats.

Nucleotide diversity (π) and haplotype diversity (Hd) were calculated for Southeast Asian goat populations. Table 1 shows the nucleotide and haplotype diversities for haplogroups A and B in Southeast Asia. The Philippine goat population showed higher genetic diversity in haplogroup A (π = 0.01390, Hd = 0.813) compared to haplogroup B (π = 0.00224, Hd = 0.582). In the Indonesian goat population, the nucleotide and haplotype diversities for haplogroup B were 0.00286 (π) and 0.789 (Hd), respectively. In ISEA, Indonesian goats show higher genetic diversity in haplogroup B than Philippine goats. The nucleotide diversity for haplogroup A in Philippine goats was relatively low compared to other Southeast Asian goat populations, while haplogroup B in the Philippines showed the lowest nucleotide diversity across Southeast Asia. The nucleotide diversity in Indonesia was also relatively low compared to other Southeast Asian goat populations.

Table 1 Genetic diversity of haplogroup A and B in six Southeast Asian countries and eight Southeast Asian countries, respectively.
Full size table

Figure 2 shows the distribution of mtDNA haplogroups in the Old World and seven regions of ISEA. This map highlights the increase in the frequency of haplogroup B toward the southeast, as previously described23. In ISEA, haplogroup B was particularly prevalent in Mindanao Island (n = 29/30) and South Sulawesi (n = 110/112).

Fig. 2
figure 2

Geographic distribution of caprine mtDNA haplogroups in the Old World and Island Southeast Asia. The size of nodes is proportional to the number of goats present in the node. The used data were shown in Table S5. The map was generated using Microsoft PowerPoint 2019 based on a free map website by 3kaku-K (https://www.freemap.jp/item/world/world1.html).

Full size image

Table 2 presents the SRY haplotype frequencies for each region in ISEA, while Fig. 3 shows the distribution of SRY haplotypes in the Old World and seven regions of ISEA. This distribution reveals that goats in ISEA possess haplotypes Y1AA, Y1B, Y2A, and Y2B. Haplotype Y1B was observed exclusively on Mindoro Island, Philippines, likely due to the introduction of modern European breeds by the Philippine Department of Agriculture. Haplotype Y2A, which is absent in MSEA, was observed in all ISEA regions, with the highest frequency in Mindanao Island, Philippines (n = 6/16).

Table 2 SRY haplotype frequencies in Island Southeast Asia.
Full size table
Fig. 3
figure 3

Geographic distribution of caprine SRY haplotypes in the Old World and Island Southeast Asia. The size of nodes is proportional to the number of goats present in the node. The used data were shown in Table S6. The map was generated using Microsoft PowerPoint 2019 based on a free map website by 3kaku-K (https://www.freemap.jp/item/world/world1.html).

Full size image

Approximate Bayesian Computation (ABC) analysis results indicated posterior probability values for Scenario 1, Scenario 2, and Scenario 3 of 0.1762 (95% CI 0.1703–0.1820), 0.3586 (95% CI 0.3523–0.3650), and 0.4652 (95% CI 0.4586–0.4717), respectively. Scenario 3, with the highest value (0.4652), indicates two propagation routes: one from China/Taiwan and another from MSEA via Indonesia to the Philippines, followed by admixture. However, the value of Scenario 3 is close to that of Scenario 2 (0.3586), so that the possibility of non-admixture history cannot be denied. In both cases, DIYABC analysis supports two propagation routes: from China/Taiwan to the Philippines and from MSEA to Indonesia.

Discussion

The mtDNA haplogroup frequencies analyzed in this study revealed that haplogroup B in ISEA goats was more frequent than in MSEA (Fig. 2), although its genetic diversity was relatively lower compared to that of MSEA (Table 1). This lower genetic diversity may be attributed to a bottleneck effect in the propagation process due to the geographical location of ISEA. Haplogroup B was observed at higher frequencies particularly in the southeastern regions of ISEA, such as Mindanao and South Sulawesi. Additionally, the genetic diversity in the Philippines was lower than that in Indonesia, suggesting a migration event from South Sulawesi to Mindanao.

Haplogroup B, predominantly observed in Southeast Asia, has also been reported at low frequencies in South African region: South Africa (n = 3/15), Namibia (n = 2/4), and Tanzania (n = 22/627)2,24. The absence of haplogroup B in African regions north than Tanzania suggests a potential propagation route by maritime trade in the Indian Ocean. The Indian Ocean has a long history of human maritime activity, with Islamic merchants trading and spreading Islam throughout the region since the 9th century, and European countries using the African coast as a transit point for maritime trade to Southeast Asia from the 16th century. Moreover, the migration of Austronesians from Indonesia to East Africa and the Arabian Peninsula due to this maritime trade has been documented25,26. Therefore, the presence of haplogroup B in the South African region is likely associated with human migration and trade in the Indian Ocean since the Middle Ages.

A previous study analyzing SRY haplotypes worldwide found that Y1AA, Y1AB, and Y2B spread from the Fertile Crescent to East and Southeast Asia, while Y1B and Y2A emerged during the expansion process in Europe and Africa, respectively17. This study revealed that haplotype Y2A is prevalent across the entire ISEA (Fig. 3). Haplotype Y2A, a major haplotype in Africa and southern Europe, has not been observed in MSEA. Therefore, haplotype Y2A in ISEA has been genetically influenced by goats introduced from southern Europe and Africa, regions that served as base camps and relay ports during the Age of Discovery. Moreover, our results revealed regional differences in the frequency of haplotype Y2A in ISEA, which may reflect the gene flow process of Y2A. A worldwide SRY analysis also identified the East Asian-specific haplotype Y2B in northern Madagascar17.

Studies on the migration routes of canines and swine provide useful insights into the propagation of goats to Island Southeast Asia27,28. Zhang et al.27 suggested two possible routes for canine propagation: one through Taiwan to the Philippines and ultimately to Australia, and another through MSEA to Indonesia and then to Australia. These propagation are estimated to have occurred between 5000 and 2000 years ago. Similarly, Layos et al.28 proposed two possible routes for swine propagation: one through Taiwan to the Philippines (during the Neolithic) and another through MSEA to the Philippines (during the prehistoric times). The propagation of domestic animals is thought to have accompanied the migration of the Austronesian language family, estimated to have occurred after 5000 years ago in the Philippines and between 4000 and 3000 years ago in Indonesia29. The Philippines, located at the crossroads of past human migrations in the Asia-Pacific region, has shown genetic effects from several times of human migration30. Indonesians are suggested to have experienced parallel introductions of MSEA- and Austronesian-related ancestry sources between 1000 and 2500 years ago31,32.

Considering the migration and expansion of humans, dogs, and pigs, goats may have migrated through similar routes. The distribution of mtDNA haplogroups in ISEA showed that the frequency of haplogroup B tended to increase toward the southeast (Fig. 2). Within ISEA, haplogroup B is particularly prevalent in South Sulawesi and Mindanao, suggesting possible propagation between these regions due to their geographical proximity. Taking these results into consideration, a plausible hypothesis for goat propagation in ISEA is that haplogroup A arrived in the Philippines via China/Taiwan, while goats with haplogroup B flowed from the Malay Peninsula to the Philippines via Indonesian islands. Simulation results from DIYABC indicated a low probability for a single propagation route to ISEA via Indochina (Fig. 4). Although Scenario 3, with the highest posterior probability value, suggests genetic admixture in the Philippines, the actual occurrence of this admixture event remains unclear, given the close posterior probability values between Scenario 3 (0.4652) and Scenario 2 (0.3586).

Fig. 4
figure 4

Evolutionary scenarios for the propagation route of ISEA goats. Scenario 1: ISEA goats were introduced from MSEA. Scenario 2: Philippine and Indonesian goats were independently introduced from China/Taiwan and MSEA, respectively. Scenario 3: Philippine goats were introduced by two routes, from China/Taiwan and from MSEA via Indonesia, with admixture occurring in the Philippines. Upper figures show scenarios indicated by the trees; lower figures show geographic regions on maps. Colors of the tree correspond to the geographic regions on the map (Green: Philippines; Red: Indonesia; Light blue: China; Pink: MSEA; Yellow: Outgroup). The map was generated using Microsoft PowerPoint 2019 based on a free map website by 3kaku-K (https://www.freemap.jp/item/world/world1.html).

Full size image

Paternal analysis revealed that haplotype Y2A, which is not found in MSEA, was observed in ISEA. The ISEA has historically been an important trading hub, especially in the spice trade, linking Europe, China, India, and the Arab world. This was due to the activities of Muslim merchants from the nineth century, followed by the Spanish and Dutch colonization of the Philippines and Indonesia, respectively, from the 16th century19,33,34. Therefore, it is possible that some exotic goats were introduced to the ISEA via maritime trade from the Middle Ages onwards in the Indian Ocean, from Europe, Africa, and India. Furthermore, the Maluku Islands in eastern Indonesia, the original source of clove and nutmeg, and South Sulawesi and Mindanao, which are located nearby, were served as trading centers by European countries. This may explain the relatively high frequency of Y2A in these regions (Fig. 3). The presence of mtDNA haplogroup B and the SRY Y2B haplotype in the South African region and northern Madagascar further supports gene introgression between the South African region and ISEA through Indian Ocean maritime trade.

This hypothesis is supported by previous research on chickens using mtDNA, which demonstrated that human migrations across the Indian Ocean led to gene flow of chickens from Southeast Asia to East Africa35. There are also historical records of chickens being transported aboard ships as a source of protein during voyages in the Medieval Period36,37. Similarly, goats were used as a protein source on maritime expeditions during the Medieval Period36, and comparable gene flow may have occurred in goats. Therefore, the genetic structure of domestic goats in the ISEA may reflect the influence of human maritime trade during the Medieval Period, suggesting that their genetic structures reflect the extent of maritime trade activities.

In this study, we estimated the genetic structure and distribution patterns of goats in the Old World using maternal and paternal analyses. Our results suggested (1) two main migration routes for ISEA goats, and (2) reciprocal gene flow between the ISEA and the South African region through transoceanic trade across the Indian Ocean. The genetic information obtained from this study provides important insights into the conservation of genetic resources and diversity of goats in the Old World, particularly in the ISEA.

Materials and methods

Ethical statement

This study is reported in accordance with ARRIVE guidelines (https://arriveguid.elines.org). This study was approved by the University of the Philippines, Los Banos, under the attestations the H23Oct285/2–9 and the Hassanuddin University, Indonesia, under the attestations the 15/25Jun2018. All experiments were conducted in accordance with the Kobe University Animal Experimentation Regulations, and all protocols were approved by the Institutional Animal Care and Use Committee of Kobe University and by Association for the Promotion of Research Integrity (Approval Number: AP0000436777). All blood sample collections were approved by the animal owners with informed consent.

Animals

A total of 248 blood samples were collected from indigenous goats in the Philippines and Indonesia (Fig. 2 and Supplementary Fig. S4). In the Philippines, 176 blood samples were collected from seven islands and six regions: 34 from Luzon, 26 from Mindoro, 30 from Samar and Leyte, 26 from Panay, 30 from Palawan, and 30 from Mindanao (Table S1). In Indonesia, 72 blood samples were collected from South Sulawesi (Table S2). Efforts were made to avoid sampling from related individuals. DNA was extracted from the blood samples using the standard phenol/chloroform method38.

Sequencing

The amplification and sequencing of the mtDNA HV1 region (481 bp) and the SRY 3’UTR region (478 bp) were performed following the procedures described previously8,11,17,39,40 (Tables S3 and S4). All sequences have been deposited in the DDBJ database (accession numbers: LC849824–LC850069, LC856305, LC856306).

Statistical analysis

Sequence alignment of the HV1 region was performed according to a previous study39. The sequence of the mtDNA HV1 region was aligned using the ClustalW program41 in MEGA 7.042, based on the mtDNA reference sequence (DDBJ accession no. AF533441.1)43. An unrooted neighbor-joining phylogenetic tree was constructed using the Tamura-Nei distance method42. The distance computation and phylogenetic tree construction were performed in MEGA 7.0. The phylogenetic trees were assessed by bootstrap percentages computed with 1000 replications44. For phylogenetic tree construction, 22 reference sequences from six goat haplogroups were used as references2 (DDBJ accession nos. AY155721, EF618134, EF617779, EF618200, EF617945, EF617965, AB044303, EF617706, AJ317833, DQ121578, AY155708, AJ317838, EF618413, DQ188892, AY155952, EF617701, DQ188893, DQ241349, DQ241351, EF618084, EF618535, EF617727). Median-joining networks45 were constructed using Network 5.0.0.1 (https://www.fluxus-engineering.com/sharenet_rn.htm). Nucleotide diversity (π) and haplotype diversity (Hd) were calculated for HV1 sequences, in conjunction with sequences from Southeast Asian native goats deposited in public databases (Table S5) using DnaSP 6.12.0346.

To infer relationships between mtDNA haplogroups in the ISEA and other regions, we analyzed the distribution of maternal haplogroups in the Eurasian and African continents using previously published mtDNA data from 5939 goats. The used data were shown in Table S5.

The distribution of paternal haplotypes in Eurasia and Africa was analyzed using SRY data from 2001 goats deposited in public databases (Table S6). The SRY haplotype data for native goats from the Philippines and Indonesia were previously published in Vargoats Consortium et al.17. In this study, we calculated the frequency of paternal haplotypes in seven regions within the ISEA and inferred relationships between the ISEA and other regions.

Evolutionary scenarios on the propagation route of ISEA goats were evaluated using Approximate Bayesian computation (ABC) with the DIYABC program ver. 2.1.047, based on the mtDNA HV1 region. Three alternative scenarios were devised to assess the evolutionary relationships among four populations: the Philippines, Indonesia, China, and MSEA. Scenario 1: ISEA goats were propagated from MSEA. Scenario 2: Goats in the Philippines and Indonesia propagated from China/Taiwan and MSEA, respectively. Scenario 3: Philippine goats propagated via two routes—China/Taiwan and MSEA via Indonesia—and admixture occurred in the Philippines.

The following parameters were set: Ta, the time before goat domestication (10,000–470,000 years ago); Tb, the time when goat domestication occurred (10,000–15,000 years ago); Tc and Td, the times when goats propagated from MSEA to Indonesia and from China/Taiwan to the Philippines, respectively. We assumed Tc to be from 4000 years ago and Td to be between 5000 and 10,000 years ago. Other parameters were left as default. The mitochondrial genome sequence of haplogroup F of the Bezoar (Capra aegagrus) was used as the outgroup7. The mutation rate was estimated using the dataset from Colli et al.7 and was set between 3.75 × 10− 8 /site/generations and 3.75 × 10− 6 /site/generations. Simulations were conducted 3 million times with the HKY + I + Γ model under the three scenarios. To identify the most likely scenario, relative posterior probability with 95% confidence intervals (95% CI) was estimated for each scenario. Model validity was checked using the perform model checking function in DIYABC. Principal component analysis in DIYABC showed that the observed dataset belongs to the point cluster of posterior distribution (Supplementary Fig. S5), indicating that the model-posterior combination fits the observed data.