Origin of the Bunun Indigenous People of Taiwan, a Review of Published Material Using Y-Chromosome and Mitochondrial DNA Gene Systems

: Anthropological and linguistic studies place the ﬁrst settlement of Austronesian speaking Taiwanese (AN_Tw) in the mid-Holocene era. However, geneticists have revealed exclusive diversity among the Bunun indigenous people, implying that their ancestral origin needs further study. The mitochondrial DNA (mtDNA) polymorphism of the Bunun shows a homogeneous relationship with other AN_Tw. However, the Y-chromosome polymorphism shows two major haplogroups, O1a2-M50 (60.7%), also seen to a lesser extent among the Northern AN_Tw, and O1b1a1a1a1a1-M88 (37.5%), scarce among all other AN and non-AN groups in Taiwan, but prevailing in Southeast Asia (SEA) and Mainland Southeast Asia (MSEA). While the present-day mtDNA proﬁle of the Bunun typiﬁes the long-term demographic standard for all AN_Tw since the Neolithic era, their Y-chromosome proﬁle suggests an arrival of male settlers in the last two to three millennia from SEA or MSEA, who mixed exclusively with the Bunun indigenous people.

To this day, the Babuza, Basay, Hoanya, Ketagalan, Luilang, Makatao, Pazeh/Kaxabu, Papora, Qauqaut, Siraya, Taokas, and Trobiawan groups are not recognized as indigenous people by the central government. They are often referred to as the Taiwan Western Plain tribes or more formally as the Pingpu. Together, the Pingpu represent 0.5% of the Taiwan population and speak Mandarin, since their original languages are extinct.
Ethnologists and archaeologists have distinguished several ancient indigenous cultures in Taiwan. They generally propose that these cultures trace their origin with the Ta-Pen-Kang culture (TPK) in the mid-Holocene period (7000 to 5000 kya) [2][3][4]. However, the discovery of Paleolithic artifacts of a tool industry at the "Changbin culture" site (28 kya) near Taitung in Southeast Taiwan, and ancient human remains (bones and teeth of the Tso-Chen Man) lead archaeologists to propose that Taiwan entered human history in the Paleolithic era and that these people were later replaced by settlers of the Ta-Pen-Kang culture [5,6]. The currently dominant Bunun language, Isbukun, is spoken in the south. Takbunuaz and Takivatan dialects are spoken by the Central Bunun, and Takituduh and Takibaka are spoken in the north. In Bold: Location of Taiwan, officially recognized indigenous people (AN_Tw). In Grey: Not officially recognized indigenous groups (also referred to as Pingpu). Reference data of other groups in East and Southeast Asia were obtained from the literature. The currently dominant Bunun language, Isbukun, is spoken in the south. Takbunuaz and Takivatan dialects are spoken by the Central Bunun, and Takituduh and Takibaka are spoken in the north. In Bold: Location of Taiwan, officially recognized indigenous people (AN_Tw). In Grey: Not officially recognized indigenous groups (also referred to as Pingpu). Reference data of other groups in East and Southeast Asia were obtained from the literature. Furthermore, linguistic studies trace all AN_Tw languages back to a common ancestral language, the proto-Austronesian language (PAN), whose filiations are shared by extant groups of Mainland Southeast Asia (MSEA), Peninsular Asia (Malaysia), Polynesia, and Madagascar [7].
The Bunun, with approximately 51,000 members represent approximately 12% of the indigenous population (or 0.23% of Taiwan population), and are the fourth-largest indigenous group in Taiwan [1,8]. In terms of family, kinship and, descent principles, the Bunun are classified, along the Tsou and Saisiyat, as a patrilineal clannish society [9]. Similar to most AN_Tw, they live in family groups, often in mono-ethnic villages that are widely dispersed in the central mountain range. In the past, their staple food was foxtail millet, and females were responsible for gathering and preparing food. On the other hand, males were very individualist, responsible for a slash and burn agriculture, hunting, and are known for having been head-hunter and fierce warriors [10].
The beginning of Taiwan's recorded history started with the arrival of the Dutch, in 1624 who established a colonial regime to expand their deerskin trade, and the Spanish later who established a base in Northern Taiwan [11,12]. Along with these new settlers was the entrance of missionaries, who translated the Bible into Bunun language and indigenized the Christian conceptual universe. By 1940, many AN_Tw were converted to Christianity. The Bunun lost many traits of their cultural heritage, and as a result, revived present-day  cultural features have become poorly supported and their genuineness debatable. Under Japanese rules, from 1895 to 1945, for security reasons against the slash and burn leading to deforestation, 62% of Bunun people were forcibly resettled in villages of the Western lowlands. Traditional hunting was also prohibited, and despite Bunun Chief Raho Ari who embarked on guerrilla warfare against the Japanese authorities, rice cultivation was imposed as the preferred alternative [10]. In October 1945, the Chinese Nationalist Kuomintang entered Taiwan and encouraged the use of standard Mandarin. This generated drastic changes in their culture. The family set up, that the Bunun were accustomed to, ended and led to a mode of living in individual units. However, the Bunun people endlessly pursue the revival of their cultural heritage. They are best known for the ear-shooting ritual that recreates their ancient practice for young man passage to adulthood, and for the "Pasibutbut", a singing prayer using an uninterrupted polyphonic song performed by eight singers still performed during the annual Millet Harvest Festival [13]. Linguistic factors are well characterized and well identifiable from other groups. Presently, the Bunun have five distinct communities whose extant dialects have been divided into three main branches [14]: The Northern Bunun (Takituduh and Takibaka), the Central Bunun (Takbanuaz and Takivatan also spoken in Nantou and Hualien counties), and the dominant Southern Bunun dialect (Isbukun). Although many Bunun under the age of 30 today are monolingual Mandarin speakers [15], the revival of their strong ethnic heritage by a large number of Bunun individuals is using Isbukun as language re-establishment [9].
Genetic studies aiming at characterizing the genetic heritage of the Taiwanese started in 1994 with Umetsu et al. [16]. The authors analyzed the polymorphisms of complement components I and C1R in nine Aboriginal Taiwanese populations, and for the first time, suggested long-term genetic and geographic isolation. Within a short time, other genetic studies within the Taiwanese Han and indigenous populations followed. The ABO and Rhesus blood groups [17], the genetic polymorphisms of the promoter region of dopamine D2 receptor and dopamine transporter [18], and a comprehensive study of six Alu insertions [19], contributed to the genetic profile of the Han and Indigenous Taiwanese.
Furthermore, a study using the Histocompatibility Leucocyte Antigens (HLA-A, -B and -DRB1), a highly polymorphic diploid gene system on chromosome 6, showed that Bunun was a member of the genetic cluster composed uniquely of AN_Tw, which separated clearly from NAN_Tw and other east Asian groups [20][21][22].
Other studies, using the uni-parental gene system [22][23][24][25][26] showed that, except for the moderate prevalence among the Bunun of mtDNA haplogroups B4b1a2, B5a2a2b, F4b1, and E1a1a, the Bunun appeared as a group maternally closer to the Northern AN_Tw. Conversely, the low Y-chromosome diversity and the high prevalence of Y-chromosome haplogroup O1a2-M50 and O1b1a1a1a1a1-M8 make the Bunun a paternally distinct outlier to all other Taiwanese groups. However, these characteristics have been attributed to founder effect, bottleneck, and long isolation. Most studies divided Taiwan into two groups, the Han Taiwanese and the Taiwan Indigenous people, and rarely [22,27,28] concentrated on tracing the genetic heritage of a single group. As a result, the genetic heritage of Bunun is still believed to be identical to the other indigenous groups in Taiwan and remains poorly characterized.
The analysis proposed here will combine evidences from two uni-parental gene systems, the non-recombining Y-chromosome (NRY) and the mitochondrial DNA (mtDNA) to compare the genetic structure of the Bunun with other Formosan speaking people, and other groups of East Asia and Island Southeast Asia (ISEA). While the mtDNA findings support previous HLA findings, the extant paternal profile suggests that the Bunun males likely reached Taiwan later in the mid-Neolithic separately from other Austronesian speaking groups. For this analysis, 235 Northern Bunun partial mtDNA sequences comprising segments of the HVS-I and coding regions were compared to 4423 mtDNA sequences from population groups of SEA, MSEA, ISEA, and Taiwan (Table 1, Supplementary Text S1). Using haplogroup assignments according to mtDNA Haplotree-Family Tree 17 [29], 699 different haplogroups were found. Compared to other ethnic groups of East Asia, SEA, and MSEA, the mtDNA polymorphism (the number of haplogroups) seen in each of the three Bunun groups (Northern, Central, and Miscellaneous) appeared to be reduced, but was in range with other indigenous groups of Taiwan. The Miscellaneous Bunun, a self-declared group collected in Taitung city [30], showed 20 haplogroups (Supplementary Text S1 and Table S1) from which only six were exclusive to AN_Tw, and fourteen were of Han origin. However, the table shows that 60.36% are shared with AN_Tw and 10.1% are shared with Han. This suggests that the remaining 20.54% is exclusive to the miscellaneous group or represents a mixture with a group other than the putative parents. This pattern is similar for the Northern Bunun, whereas the Central Bunun appears to be the closest group to the AN_Tw putative parental population. Pairwise comparisons ( Figure 2) between the Northern Bunun and the Miscellaneous Bunun group did not show significant differences (p > 0.3). However, these two groups showed striking differences from the Central Bunun group (p < 0.03) lacking mtDNA subtypes of haplogroup F4b1 and showed a lower prevalence of haplogroup E. On the contrary, haplogroups B5a1c, D5b3a, F3b1, and M7c1 were not seen in the Northern and Miscellaneous Bunun groups, indicating that the genetic history of the Northern and Miscellaneous groups may be different from the Central Bunun.
Pairwise comparisons (Figure 2) between the Northern Bunun and the Miscellaneous Bunun group did not show significant differences (p > 0.3). However, these two groups showed striking differences from the Central Bunun group (p < 0.03) lacking mtDNA subtypes of haplogroup F4b1 and showed a lower prevalence of haplogroup E. On the contrary, haplogroups B5a1c, D5b3a, F3b1, and M7c1 were not seen in the Northern and Miscellaneous Bunun groups, indicating that the genetic history of the Northern and Miscellaneous groups may be different from the Central Bunun.
The genetic diversity of mtDNA haplogroups in the Northern Bunun (h = 0.99) was high, but in range with other indigenous groups in Taiwan (0.989 to 1.00) ( Supplementary  Table S1). Similarly, the nucleotide diversities between groups were not significantly different (π ~ 0.044) (data not shown). The three Bunun groups exhibited multimodal mismatch distribution patterns (Supplementary Figure S1) with the mean number of nucleotide-pairwise-differences (MNND) between 12.5 and 14 as well as significant "sum-of-square differences" (SSD, p ≤ 0.05) and raggedness values (p < 0.01). These results indicate a long stationary maternal population history [51]. Furthermore, the Tajima D values of the three Bunun groups were not significant (Supplementary Figure S1d), suggesting the groups were close to or under the neutrality expectation. However, a small Tajima D difference was observed between the Northern Bunun (0.3; p = 0.7) and the Central Bunun (−1.28; p = 0.086). This suggests a minor expansion event of the Central Bunun or possibly gene flow from neighboring groups. This observation is also supported in the Mismatch Distribution Analysis by the presence of two short peaks (left in Supplementary Figure S1C), possibly due to the same sampling size or the result of the admixture. This is not surprising since the Central Bunun group [33], similar to the miscellaneous group [38], were self-declared Bunun indigenous individuals (collected in the south of Taiwan in the Kaohsiung and Taitung counties, respectively).

MtDNA Polymorphism
When using only major mtDNA haplogroups, the AN_Tw can be divided into three homogeneous groups, Northern, Central, and Southern group; however, advances in the sequencing techniques have allowed for the determination of discrete differences within the groups [33,49,52]. Some of these differences can be seen throughout ISEA, sometimes at higher frequency, suggesting drift, bottleneck, and/or long isolation, and reflecting the complexity encountered when linking genetics, languages, and demography throughout ISEA.
In this analysis, the Northern Bunun profile is limited to sub-lineages of four mtDNA para-haplogroups, namely B4b1, B5a2, F4b1, and E1a1a ( Figure 2, Supplementary Table S1), that represent 80% of the Northern Bunun genetic structure, a common characteristic of the Northern and Central AN_Tw. Interestingly, these haplogroups do not represent clades that participated in the "out-of-Taiwan" [53].
The most common subtype of B4b1 in Taiwan, B4b1a2g, is seen in the Northern Bunun and Thao groups (24% and 27%, respectively) and has an estimated expansion age The genetic diversity of mtDNA haplogroups in the Northern Bunun (h = 0.99) was high, but in range with other indigenous groups in Taiwan (0.989 to 1.00) ( Supplementary  Table S1). Similarly, the nucleotide diversities between groups were not significantly different (π~0.044) (data not shown).
The three Bunun groups exhibited multimodal mismatch distribution patterns (Supplementary Figure S1) with the mean number of nucleotide-pairwise-differences (MNND) between 12.5 and 14 as well as significant "sum-of-square differences" (SSD, p ≤ 0.05) and raggedness values (p < 0.01). These results indicate a long stationary maternal population history [51]. Furthermore, the Tajima D values of the three Bunun groups were not significant (Supplementary Figure S1d), suggesting the groups were close to or under the neutrality expectation. However, a small Tajima D difference was observed between the Northern Bunun (0.3; p = 0.7) and the Central Bunun (−1.28; p = 0.086). This suggests a minor expansion event of the Central Bunun or possibly gene flow from neighboring groups. This observation is also supported in the Mismatch Distribution Analysis by the presence of two short peaks (left in Supplementary Figure S1C), possibly due to the same sampling size or the result of the admixture. This is not surprising since the Central Bunun group [33], similar to the miscellaneous group [38], were self-declared Bunun indigenous individuals (collected in the south of Taiwan in the Kaohsiung and Taitung counties, respectively).

MtDNA Polymorphism
When using only major mtDNA haplogroups, the AN_Tw can be divided into three homogeneous groups, Northern, Central, and Southern group; however, advances in the sequencing techniques have allowed for the determination of discrete differences within the groups [33,49,52]. Some of these differences can be seen throughout ISEA, sometimes at higher frequency, suggesting drift, bottleneck, and/or long isolation, and reflecting the complexity encountered when linking genetics, languages, and demography throughout ISEA.
In this analysis, the Northern Bunun profile is limited to sub-lineages of four mtDNA para-haplogroups, namely B4b1, B5a2, F4b1, and E1a1a (Figure 2, Supplementary Table S1), that represent 80% of the Northern Bunun genetic structure, a common characteristic of the Northern and Central AN_Tw. Interestingly, these haplogroups do not represent clades that participated in the "out-of-Taiwan" [53].
The most common subtype of B4b1 in Taiwan, B4b1a2g, is seen in the Northern Bunun and Thao groups (24% and 27%, respectively) and has an estimated expansion age of 4.2 ± 3.5 kya [20]. Excluding B4b1a2a, which is principally seen in SEA and MSEA, all other B4b1 subtypes in Taiwan are subclades of B4b1a2. They prevail in the Northern tribes (Saisiyat, Atayal, and Taroko) and groups of the central mountain range, Bunun, Tsou, and Thao [22,33]. They find their origin in SEA and north MSEA, likely along the Mekong river since they are absent from Vietnam and the Cham Austronesian speaking groups of Vietnam [22,33,45,46,54]. Haplogroup B4b1a2g is the most common mtDNA B subtype in Taiwan. It prevails in the Northern Bunun and Thao groups (25% and 27%, respectively). It has an estimated expansion age of 4.2 ± 3.5 kya and most likely expanded exclusively in Taiwan [22].
Haplogroup B5a is common in South and Southeast China. Subtypes of B5a are seen in the three Bunun groups and the Tsou (Supplementary Table S1). They are primarily represented by B5a2a2b, at a lower frequency by B5a2a2a, and are scarce out of the central mountain range. The type seen among Bunun is B5a2a2b1 and likely evolved in Taiwan among the Bunun (the Tsou) 4.3 ± 2.2 kya.
Haplogroup F4b1 is determined by the presence of nucleotide position np 10,097 and 513. It is a common characteristic of the Northern AN_Tw (Supplementary Table S1). With an expansion estimate of 4.3 kya (CI 1.3-7.3 kya), it is found to be a daughter haplogroup of F4b, which is common in Southeast Asia, and represents a mid-Holocene dispersal from this region in Taiwan [22,55].
Haplogroup E1a1a1 is a descendent of haplogroup E1. Commonly seen in the Saisiyat, the Tsou, and the Bunun tribes, its ancestral lineages trace back to haplogroup pre-E1, first described in the Liangdao man [33,47]. All E-lineages are branches of haplogroup M9 seen along the east coast of China and Southeast Asia [33,47]. Since no E-lineages have been reported in Mainland Asia, haplogroup E likely evolved in Island Southeast Asia. With an average founder age of 6.25 ± 1.6 kya, haplogroup E1a1a1 relates to the early dispersal of the Austronesian people to Taiwan.
Finally, lineages B4b1a2, B5a2a2b, and F4b1 prevail in Northern Bunun. Their molecular variation was used to estimate the time of gene flow events in that group. Convergence estimates showed 4.2 kya (CI 0.6-7.8 kya), 4.3 kya (CI 1.3-7.3 kya), and 4.3 kya (CI 1.2-7.5 kya), respectively. While their confidence interval encompasses the time of the Austronesian expansion, it is also likely that these lineages represent later gene flow events of the past four millennia.
Further mtDNA haplogroups with considerably lower frequencies were identified. These were R9c, M7b1a2a, and F4b1c in Northern Bunun, A5a, B4a2a3, B5a1c, D5b3a, F3b1, and M7c1c in Central Bunun, and D5b1c and F4b1a in Miscellaneous Bunun. Except for A5a, a fingerprint of Han admixture, these haplogroups were all associated with the profile of Taiwan indigenous people. Their presence elsewhere was limited to the Philippines and Western Indonesia (Supplementary Table S1). Their low frequencies were likely the result of drift.

Principal Component Analysis
As expected from the pairwise comparison (Figure 2), the principal component ( Figure 3) characterizes the Northern Bunun, the Thao, and the Tsou as outgroups. Moreover, the Principal Component Analysis (PCA) plot separates the Northern indigenous groups of Taiwan (Atayal, Taroko, and Saisiyat) seen in the lower right from the Bunun, the Thao, and the Tsou. Furthermore, the plot separates the Southern indigenous groups in the left lower pane from MSEA, the Philippines, Indonesia, and the Taiwan Han top seen on the top of the lower left quadrant.

Y-SNP
Y-SNP typing revealed only two Non-Recombining Y-Chromosome (NRY) haplogroups among 56 Northern Bunun males ( Table 2). Almost 61% belonged to haplogroup O1a2-M50 and 37.5% to haplogroup O1b1a1a1a1a1-M88 (Supplementary Table S1). The frequency of the last haplogroup is the highest seen in Taiwan and Western ISEA. One single individual was bearing haplogroup p-M45. Another individual bearing the same haplogroup was seen among the Amis. This presence of P-M45 in Taiwan is interesting since it was associated with the Negrito groups of south Luzon in the Philippines (Aeta from Santa Juliana and Agta from Umiray) [56], and its ancestry has been traced back to the upper Paleolithic [57].   (Table 2). Almost 61% belonged to haplogroup O1a2-M50 and 37.5% to haplogroup O1b1a1a1a1a1-M88 (Supplementary Table S1). The frequency of the last haplogroup is the highest seen in Taiwan and Western ISEA. One single individual was bearing haplogroup p-M45. Another individual bearing the same haplogroup was seen among the Amis. This presence of P-M45 in Taiwan is interesting since it was associated with the Negrito groups of south Luzon in the Philippines (Aeta from Santa Juliana and Agta from Umiray) [56], and its ancestry has been traced back to the upper Paleolithic [57]. Haplogroup O1a2-M50 seen in the Northern Bunun is a subclade of O1a-M119 [58]. It was reported in Cambodia [59], and a few Daic-speaking groups of Southeast Asia [32,60]. Its coalescence is estimated to be 15.8 kya (CI 17.5-14.1 kya) [61,62], but its origin in SEA is not well defined. The expansion of O1a2-M50 likely took place in the South of Taiwan, where it prevails at approximately 16.7%. Its presence in any insular groups out of Taiwan represents an out-of-Taiwan signature. From Taiwan, O1a2-M50 reached the admiralties islands (17%). In the south, it expanded to the Philippines, Borneo, Java, the South Nias Island off the coast of Sumatra (60% to 80%), Malaysia, and finally Madagascar [60,[63][64][65][66].

Y-STR Haplotype Networks
The Y-STR haplotype networks of O1a2-M50 and O1b1a1a1a1a1-M88 were constructed to evaluate their early dispersal ( Figure 4). However, only the Northern Bunun NRY data were available. All Bunun O1a2-M50 haplotypes clumped within shallow clusters on the outskirt of the network away from the ancestral core. Genetic studies have shown [60] that Taiwan is the most likely region of expansion of O1a2-M50 and it is now generally accepted by linguists, archaeologists, and geneticists [2,60,73], as well as O1a2-M50 is a male signature of the Austronesian expansion. Here, it appears that the absence of NRY haplogroup O1a1-P203 and the high prevalence of O1a2-M50 in Bunun are likely the result of drift. Similarly, the exclusive occurrence of haplogroup O1b1a1a1a1a1-M88 in the Northern Bunun (in purple on the right network of Figure 4) showed a shallow structure, supporting a paternal founding or drift event and isolation. Moreover, we notice that the haplogroup pair (O1a2-M50 and O1b1a1a1a1a1-M88) is seen along the Taiwan-Philippines-Borneo-Sumatra pathway, and can be traced further toward Madagascar and the Comoros islands [63,64,74,75]. While this distribution may suggest that O1a2-M50 and O1b1a1a1a1a1-M88 dispersed as a pair from Taiwan, the low Bunun maternal contribution to ISEA and the exclusive presence of O1b1a1a1a1a1-M88 in the Northern Bunun suggest that O1b1a1a1a1a1-M88 reached Taiwan directly from SEA and that the dispersal of O1a2-M50 from Taiwan and O1b1a1a1a1a1-M88 coming from MSEA merged in Indonesia.
The Y-STR variation of haplogroups O1a2-M50 and O1b1a1a1a1a1-M88 was then calculated based on the seven most used Y-STR loci in the NRY scientific literature. Corroborating the ISSOG dating [61,62] O1a2-M50 and O1b1a1a1a1a1-M88 coalesce at approximatelỹ 11.8 + 4.1 kya and 16.7 + 2.7 kya, respectively, suggested a late Paleolithic Southeast Asian at the beginning of their expansions [59,76]. Furthermore, estimation of the Y-STR variation exclusive to the Northern Bunun provided a TMRCA of 2.4 ± 1.6 kya for O1a2-M50 and 4.1 ± 3.1 kya for O1b1a1a1a1a1-M88. These results contrast with the SSOG TMRCAs obtained for the whole of Asia Oceania using the software Yfull, version 8.04.00 (i.e., 5.1 and 5.6 kya, respectively) [61,62,77], and suggest a separate migration from MSEA and later expansion in Taiwan among the Bunun Indigenous people.  [42] (e.g., Tonga, Samoa, Fiji).

Principal Component Analysis from NRY Haplogroups
The two dimensions of the principal component analysis (PCA) using NRY haplogroup frequencies represent 55.1% of the polymorphism variation ( Figure 5). It separated the Northern and Southern AN_Tw groups (red dots) from each other. The southern AN_Tw group shows more affinity with Mainland Southeast Asia (yellow squares) and suggests admixture with Taiwan plain dwellers Pazeh and Siraya groups (light blue dots). The Northern Bunun groups are equidistant from MSEA and Western Indonesia [25,32].  [42] (e.g., Tonga, Samoa, Fiji).

Principal Component Analysis from NRY Haplogroups
The two dimensions of the principal component analysis (PCA) using NRY haplogroup frequencies represent 55.1% of the polymorphism variation ( Figure 5). It separated the Northern and Southern AN_Tw groups (red dots) from each other. The southern AN_Tw group shows more affinity with Mainland Southeast Asia (yellow squares) and suggests admixture with Taiwan plain dwellers Pazeh and Siraya groups (light blue dots). The Northern Bunun groups are equidistant from MSEA and Western Indonesia [25,32].

TreeMix
Finally, a Maximum-Likelihood phylogram was constructed using the TreeMix package [78] to infer mtDNA and NRY gene-flow events (Figure 6), which contributed to the genetic history of Bunun.
In concordance with the previous study [79], the mtDNA phylogram (left on Figure  6) inferred gene flow to the Northern and Central tribes (Saisiyat, Taroko, and Bunun), most likely associated with the E1a clade from the Philippines (Supplementary Table S1). Similarly, gene flow associated with F4b1a was seen between Saisiyat, Atayal, Taroko, Thao, and Bunun. Finally, the mtDNA phylogram showed a series of relationships restricted to the Central AN_Tw, Thao, Tsou, and Bunun, likely associated with the high prevalence of subtypes of B4b1a2 in these groups. These results indicate gene flow between Central and Northern AN_Tw and support the findings of the PCA (Figure 3).

TreeMix
Finally, a Maximum-Likelihood phylogram was constructed using the TreeMix package [78] to infer mtDNA and NRY gene-flow events (Figure 6), which contributed to the genetic history of Bunun.
In concordance with the previous study [79], the mtDNA phylogram (left on Figure 6) inferred gene flow to the Northern and Central tribes (Saisiyat, Taroko, and Bunun), most likely associated with the E1a clade from the Philippines (Supplementary Table S1). Similarly, gene flow associated with F4b1a was seen between Saisiyat, Atayal, Taroko, Thao, and Bunun. Finally, the mtDNA phylogram showed a series of relationships restricted to the Central AN_Tw, Thao, Tsou, and Bunun, likely associated with the high prevalence of subtypes of B4b1a2 in these groups. These results indicate gene flow between Central and Northern AN_Tw and support the findings of the PCA (Figure 3). Figure 6 also shows results for the NRY-STR dataset (right). Taiwan samples were analyzed along with ISEA, MSEA, Oceania, and Madagascar populations [63,64]. For clarity, the arrows in Figure 6 only show the gene flow corresponding to migration weight, which is greater than 0.5. Most of the gene flow associated with Bunun finds its origin in MSEA/SEA, here represented by the heterogeneous ethnic groups of Lao [81] and primarily indicating the origin of haplogroup O1b1a1a1a1a1-M88. The coalescence age estimate of O1b1a1a1a1a1-M88 in Bunun was~2.5 kya (95% CI 0-5 kya) ( Figure 4). Furthermore, the highest occurrence of haplogroup O1a2-M50 in Taiwan is found in Bunun (60%). One must note that O1a2-M50 is also seen among the non-Austronesian speakers of Taiwan (TwH) and shows signs of local expansion.
Comparisons of the Bunun mtDNA and Y-chromosome diversity to other ethnic groups in Taiwan show that women moved between the groups. Conversely, the almost exclusive NRY profile of the Bunun suggests that 37.5% of the men had a phylogeographic history distinct from other AN_Tw. Most Taiwan indigenous groups are patrilineal. Patrilocality might have been the principal social element that maintained the bearers of O1b1a1a1a1a1-M88 within the Bunun community until the present day. Furthermore, genetic variation in the Bunun shows high within-group diversity for the mtDNA (Supplementary Table S1) and a low Y-chromosome diversity ( Table 2) with large between-group distances [25], which is in range with a patrilocal residence expected pattern. Finally, the mtDNA gene flow detected by TreeMix ( Figure 6) was principally restricted to the gene flow between Taiwan tribes, whereas the Y-chromosome gene flow covered remote geographic regions. This suggests that the male dispersal pattern in the past covered longer geographical distances compared to the female dispersal pattern.  Table S1, and mitochondrial studies [27,80]. Righ Gene flow of NRY-STR for haplogroups O1b1a1a1a1a1-M88. Note: Haplogroups shown on a yellow background are inferences from Table 2 and Y-chromosome studies [25,68]. TreeMix detected gen flow for mtDNA from closer geographic regions compared to the Y-chromosome. This suggests tha male dispersal patterns cover longer geographical distances compared to female patterns. Figure 6 also shows results for the NRY-STR dataset (right). Taiwan samples wer analyzed along with ISEA, MSEA, Oceania, and Madagascar populations [63,64]. For clar ity, the arrows in Figure 6 only show the gene flow corresponding to migration weigh which is greater than 0.5. Most of the gene flow associated with Bunun finds its origin i MSEA/SEA, here represented by the heterogeneous ethnic groups of Lao [81] and primar ily indicating the origin of haplogroup O1b1a1a1a1a1-M88. The coalescence age estimat of O1b1a1a1a1a1-M88 in Bunun was ~2.5 kya (95% CI 0-5 kya) (Figure 4). Furthermore the highest occurrence of haplogroup O1a2-M50 in Taiwan is found in Bunun (60%). On must note that O1a2-M50 is also seen among the non-Austronesian speakers of Taiwa (TwH) and shows signs of local expansion.
Comparisons of the Bunun mtDNA and Y-chromosome diversity to other ethni groups in Taiwan show that women moved between the groups. Conversely, the almos  Table S1, and mitochondrial studies [27,80]. Right: Gene flow of NRY-STR for haplogroups O1b1a1a1a1a1-M88. Note: Haplogroups shown on a yellow background are inferences from Table 2 and Y-chromosome studies [25,68]. TreeMix detected gene flow for mtDNA from closer geographic regions compared to the Y-chromosome. This suggests that male dispersal patterns cover longer geographical distances compared to female patterns.

Conclusions
The Bunun indigenous people hold significant cultural, anthropological, historical, and linguistic diversity, which differentiates them from other groups in Taiwan. While previous HLA and mitochondrial DNA studies posited the Bunun as a group distinct to the Southern Taiwan indigenous groups, their Y-Chromosome profile showed a mixture of ancestral origin. First, haplogroup O1a2-M50, commonly seen among other AN_Tw and Southeast Asia, showed evidence of expansion exclusive to the Bunun (at 1.6 to 2.4 kya). Second, haplogroup O1b1a1a1a1a1-M88 (O2a1a-M88), frequently seen among Daic and Hmong-Mien populations [32,82,83] of Southeast Asia and Mainland Southeast Asia, was primarily seen in the Bunun group and scarcely elsewhere in Taiwan. This suggests a mid-Neolithic (1.6 to 4.1 kya) intrusive migration event of men from Mainland Southeast Asia (Figure 4) into the Bunun groups. This migration event contributed to the modern genetic pool of the Bunun. In this case, it is likely that the strong cultural practices of the Bunun, such as animism, their strict sexual division of labor, hunting, cultivation, language, and social prestige were placed aside, even replaced by those of the newcomers. However, it had a significant motivating impact on the Bunun of today to revive their initial Austronesian heritage. Moreover, this uni-parental genetic review allowed for the determination of the dual origin of Bunun people.  Table S1: Polymorphisms of mtDNA exclusive to the Northern Bunun and comparison with other populations, Supplementary Text S1: Materials and Methods. References [22,24,25,27,29,31,33,34,37,38,45,48,50,51,56,63,64,66,68,75,78,80,[84][85][86][87][88][89][90][91][92][93] are cited in the supplementary materials.
Funding: This research received no funding.

Institutional Review Board Statement: Not applicable.
Informed Consent Statement: Not applicable.

Data Availability Statement:
This study is a Review of published material; the data presented are available from articles described in Table 1.