New insights from Thailand into the maternal genetic history of Mainland Southeast Asia

Tai-Kadai (TK) is one of the major language families in Mainland Southeast Asia (MSEA), with a concentration in the area of Thailand and Laos. Our previous study of 1234 mtDNA genome sequences supported a demic diffusion scenario in the spread of TK languages from southern China to Laos as well as northern and northeastern Thailand. Here we add an additional 560 mtDNA genomes from 22 groups, with a focus on the TK-speaking central Thai people and the Sino-Tibetan speaking Karen. We find extensive diversity, including 62 haplogroups not reported previously from this region. Demic diffusion is still a preferable scenario for central Thais, emphasizing the expansion of TK people through MSEA, although there is also some support for gene flow between central Thai and native Austroasiatic speaking Mon and Khmer. We also tested competing models concerning the genetic relationships of groups from the major MSEA languages, and found support for an ancestral relationship of TK and Austronesian-speaking groups.


Introduction
The geography of Thailand encompasses both upland and lowland areas, and Thailand is one of the most ethnolinguistically-diverse countries in Mainland Southeast Asia (MSEA). With a census size of~68 million in 2015, there are 70 different recognized languages belonging to five different major language families: Tai-Kadai (TK) (90.5%), Austroasiatic (AA) (4.0%), Sino-Tibetan (ST) (3.2%), Austronesian (AN) (2.0%), and Hmong-Mien (HM) (0.3%) [1]. The majority of the people (29.72%) are called Thai or Siamese and speak a central Thai (CT) language that belongs to the TK family. Since it is the country's official language, the number of people speaking the CT language as their primary or secondary language is~40 million [1], or~68% of the population.
The recorded history of the CT people or Siamese started with the Sukhothai Kingdom, around the 13th century A.D [2]. However, before the rise of the TK civilization, Thailand was under the control of Mon and Khmer people [3,4]. Linguistic and archaeological evidence suggests that the prehistorical TK homeland was situated in the area of southeastern or southern China, and that they then spread southward to MSEA around 1-2 kya [5,6]. This process could have occurred via demic diffusion (i.e., a migration of

Sequencing
We generated complete mtDNA sequences from genomic libraries with double indices and mtDNA enrichment based on protocols described previously [14,15]. The libraries were sequenced on the Illumina Hiseq 2500. MtDNA consensus sequences were obtained as described by Arias-Alvis et al. [16]. except that Illumina standard base calling was performed using Bustard and the read length was 76 bp. Sequences were manually checked with Bioedit (www. mbio.ncsu.edu/BioEdit/bioedit.html). A multiple sequence alignment of the sequences and the Reconstructed Sapiens Reference Sequence (RSRS) [17] was obtained by MAFFT 7.271 [18]. The 560 mtDNA genomic sequences reported in this study have been deposited in NCBI GenBank under accession numbers MG272576-MG273135.

Statistical analyses
Haplogroup assignment was performed with the online tools Haplogrep [19] and MitoTool [20]. Arlequin 3.5.1.3 was used to obtain summary statistics [21]. For the population comparisons, we included an additional 1234 mtDNA genomes from 51 Thai/Lao populations from our previous study [7] (Supplementary Table S1), for a total of 1794 sequences from 73 populations (Fig. 1). The matrix of genetic distances (Φ st , pairwise difference), permutation test, analyses of molecular variance (AMOVA), and a Mantel test of the correlation between genetic and geographic distances were also carried out with Arlequin [21]. Three types of geographic distances were computed, as previously described [7]. To get a broad picture of population relationships in Asia, we included 1936 published mtDNA genomes from 61 Asian populations (Supplementary Table S1) and calculated the Φ st matrix by Arlequin [21].
The Φ st distance matrix was visualized by a multidimensional scaling plot (MDS) using STATISTICA 10.0 (StatSoft, Inc., USA). A Discriminant Analysis of Principal Components (DAPC) was employed using the dapc function within the adegenet R package [22]. Median-joining networks [23] of haplogroups without pre-processing and post-processing steps were constructed with Network (www.fluxus-engineering.com) and visualized in Network publisher 1.3.0.0.
Bayesian skyline plots (BSP) per population and maximum clade credibility (MCC) trees per haplogroup, based on Bayesian Markov Chain Monte Carlo (MCMC) analyses, were constructed using BEAST 1.8.0 [24]. BEAST input files were created with BEAUTi v1.8.0 after first running jModel test 2.1.7 in order to choose the most suitable model of sequence evolution [25]. BSP calculations per population and the BEAST runs by haplogroup were executed with respective mutation rates of 1.708 × 10 −8 and 9.883 × 10 −8 for data partitioned between coding and noncoding regions [26] and Tracer 1.6 was used to generate the BSP plot from BEAST results. The Bayesian MCMC An approximate Bayesian computation (ABC) approach was utilized to test different demographic scenarios concerning the origin of CT populations and the relationships between SEA language families. For the maternal origin of CT populations, we considered the same three demographic scenarios tested in our previous study for the origins of North/Northeastern Thai and Laos populations: [7] demic diffusion (Fig. 2a); an endogenous origin (with cultural diffusion of the TK language) (Fig. 2b); and continuous migration (Fig. 2c). For testing the genetic relationships of populations from the different SEA language families, we included populations speaking AA, AN, ST, and TK languages but excluded HM because of its low population size in SEA and limited mtDNA genome data. We analyzed five tree-like demographic histories based on linguistic data [8][9][10][11] for Model 1-Model 3 ( Fig. 3a-c) and based on the geographic distribution of these languages for Model 4 and Model 5 (Fig. 3d,e). Since the AA, TK, and ST are the languages spoken in MSEA while AN is the major language in ISEA, Model 4 and Model 5 propose a closer affinity of AA, TK, and ST and set AN as an outgroup. Model 4 postulates an AA-TK affinity while Model 5 is a trifurcation of AA, TK and ST. Because of the computational cost of simulating a large number of complete mitochondrial sequences, we utilized a novel approach [27] based on a machine learning tool called "random forests" (RF) [28]. Additional details concerning the ABC-RF analyses are described in Supplementary Text.
We observed 174 haplogroups among the 560 sequences (Supplementary Table S2); when combined with our previous study [7] of Thai/Lao populations, there are a total of 1794 sequences from 73 populations (Fig. 1). In total there are 1103 haplotypes and 274 haplogroups, of which 62 haplogroups were not observed in the previous study [7] (Supplementary Table S3). An analysis of haplotype sharing (Supplementary Figure S1) shows that all four Karen groups (KSK1, KSK2, KPW, and KPA) share haplotypes, indicating high gene flow among them. The Mon (MO6-MO7) shared haplotypes with several other ethnic groups, e.g., Yuan (YU) and Central Thai (CT), whereas most of the CT populations shared haplotypes more often with northeastern Thai than northern Thai groups (Supplementary Figure S1). The AMOVA revealed that overall, 7.10% of the genetic variation is among populations ( Table 2). Classifying populations by language family resulted in a slightly higher proportion of variation among groups (0.91%, P < 0.01) than a geographic classification (0.17%, P > 0.01), indicating that language family seems to be a better indicator than geography of the genetic structure of Thai/Lao populations, however there is much more variation among  [7]. Population abbreviations are in Supplementary Table S1  Table 2). Within each language family, the variation among AA groups (11.14%) was greater than that of ST (6.51%) or TK (4.59%) groups, indicating greater genetic heterogeneity of AA groups. Interestingly, we observed that the CT groups are the most homogenous of the TK groups,  with only 1.64% of the variation among groups. However, Lue groups had higher heterogeneity (7.26%) than the average for TK groups (4.59%). A Mantel test for correlations between genetic and geographic distances indicates no correlation for all three types of geographic distances, i.e. great circle distance (r = 0.0216, P > 0.01), resistance distance (r = −0.0996, P > 0.01) and least-cost path distance (r = 0.0459, P > 0.01), further supporting the limited impact of geography on the genetic structure of Thai/Lao populations. Furthermore, a DAPC analysis showed that clustering groups by language family resulted in more discrimination among groups than clustering by geographic criteria (Supplementary Figure S2). The MDS showed that the most differentiated groups were two H'tin groups (TN2 and TN1) and Seak (SK), as found previously [7] and in the central cloud of the plot it is difficult to see population clustering trends (Supplementary Figure S3). Genetic distance values exhibited larger significantly genetic difference of the AA populations than the TK populations, supporting greater genetic divergence of the AA groups (Supplementary Figure S4). After omitting these outlier groups (TN1, TN2, and SK), a 3-dimensional MDS provides an acceptable fit (Fig. 4a-c)

MtDNA haplogroups
Fourteen of the 174 haplogroups occur in at least ten individuals and together account for 33 Table S3).
In the combined Thai/Lao dataset, SEA specific haplogroups (B, F, and M7) are prevalent in almost all groups (overall frequency 55.18%), with the exception of some AA groups (i.e., Mon, Suay, Nyahkur, Khmer, and Lawa), Karen, and CT groups; these groups have other widespread haplogroups, e.g., D, M12-G, M (xM12-G, M7), A, C, and N (xN9a) (Fig. 1). Networks of common SEA specific haplogroups, e.g., B5a, F1a, F1f, and M7b, tend to exhibit star-like structures, indicative of population expansions (Supplementary Figure S6). Apart from F1a1a (xF1a1a1), other more-prevalent haplogroups of Karen (B6a1a and C7a1) do not show indications of population expansion, but rather sharing of sequences, suggesting population contraction (Supplementary Figure S6). Apart from B and F1, other lineages, that is, C7a1 and A17 and N8 which are sublineages of C, A, and N (xN9a), respectively are observed in the Karen (Fig. 1). Haplogroup C7 has a very high frequency in northeast Asia and eastern India [29] while haplogroup A was previously reported to be specific to North and Central Asia [30]. A high proportion of C and A lineages were previously observed in ST-speaking Barmar and Karen from Myanmar [31]. For the TK-specific haplogroups, i.e. B4 and M7c, there was no obvious signal of population expansion in the networks (Supplementary Figure S6).
For the combined dataset, we estimated coalescence ages of SEA haplogroups and their sublineages. We analyzed haplogroups that have additional sequences from the present study and have more than five sequences in total ( Table 3). The ages of major haplogroups are generally consistent with previous studies [7]. However, we obtained more data from several sublineages which were not dated previously, e.g., B4c1b, B6a1, C4, C7a, D4a, F1c, F1e, F1g, F2, F3, F4a2, and G2a (Table 3).
There are many lineages with ages older than 30 kya found in our Thai/Lao samples, e.g.  Table 3). Many of them are major lineages and distributed in our Thai/Lao samples as well as in other SEA populations, and have been previously discussed [7]. Here, we focused on some uncommon ancient lineages, i.e., M*, M17, M21, M71, M73, M91 and U; these  Supplementary  Table S1 New insights from Thailand into the maternal genetic...   [32]. Finally, several haplogroups associated with the AS expansion from Taiwan [33][34][35][36], namely B4a1a1a, M7b3, M7c3c, E1a1a, and Y2 were not observed, suggesting that this expansion had at most a limited impact on mtDNA lineages in MSEA.

Bayesian skyline plots
BSP of population size change over time were constructed for each group, and five typical patterns were observed (Fig. 5). The four Karen populations all showed different patterns: KSK2 (and also MO6 and LU4) displayed unchanged population size until~1-2 kya followed by sharp reductions (Fig. 5, pattern a); KSK1 was also constant in size, with a sudden increase in the last 1-2 kya (Fig. 5,  pattern b); KPA was basically constant in size over time (Fig. 5, pattern c); and KPW exhibited the most common pattern (also observed in MO7, KPW, TKH, LU1-LU2, YU3-YU6, CT6-CT7), consisting of population expansion between 50-60 kya, followed by a decrease in the last 5 kya (Fig. 5, pattern d). While recent reductions in population size could reflect recent bottlenecks, such changes during recent times should be interpreted cautiously as they may reflect a bias in sampling [37]. Finally, population growth without further change was found for LU3 and CT1-CT5 (Fig. 5, pattern e). The BSP plots for each individual population are depicted in Supplementary Figure S7.

Demographic models for the origin of central Thai people
In our previous study we used demographic modeling to show that northern and northeastern Thai groups most likely originated via demic diffusion from southern China [7]. Here we use the same approach to test three demographic scenarios concerning the origins of central Thai groups: (1) descent from the prehistorical Tai stock of southern China via demic diffusion, like their neighbors in the North and Northeast of Thailand (Fig. 2a); (2) local AA groups (Mon and Khmer) who changed their identity and language via cultural diffusion to become TK groups (Fig. 2b); or (3) descent from a migration from southern China that received gene-flow from the local Mon and Khmer people (Fig. 2c).
The estimated values of ABC-RF prior error rate with respect to the number of trees in the forest indicated that 500 trees were sufficient in our analyses (Supplementary   Table S4). This result is understandable, as the two scenarios only differ for the events of migration hypothesized between the AA and CT groups. However, the demic diffusion model had the highest posterior probability at 0.604 and was selected slightly more often among the classification trees (0.515) than the continuous migration model (0.404); both of them were selected much more often than the cultural diffusion model (0.081) (Supplementary  Table S5). Moreover, the estimated parameters for the continuous migration model indicate very low level of migration between CT and AA groups (Supplementary  Table S6, Supplementary Figure S10). We conclude that demic diffusion, possibly with a very low level of gene flow between CT and AA groups, is the most likely scenario for the origins of central Thai populations.

Genetic relationships of populations from different language families
We also used the demographic modeling approach to test different models for the genetic relationships of populations belonging to the four main SEA language families (TK, AA, AN, and ST). In doing so, it is important to keep in mind that we are not testing the relationships of these language families, as that would require linguistic data. However, determining the best-fitting model based on genetic relationships may help discriminate among hypotheses concerning the language family relationships that make predictions about the genetic relationships of populations speaking those languages. We tested five models of the language family relationships (Fig. 3). The observed data fall within the range of the simulated data in the LDA plot (Supplementary Figure S9) while the confusion matrix confirms our ability to distinguish amongst the different scenarios (Supplementary Table S7). The model that best fit the mtDNA genome data was Model 1, according to Starosta [11] (Fig. 3a). The posterior probability of this model is 0.657, and it was selected slightly more often among the classification trees (0.509) than Model 2 (0.311); the other models were much less often selected among the classification trees (0.037 for Model 3; 0.112 for Model 4; 0.031 for Model 5) (Supplementary Table S8). Because of the high selection frequency of Model 1 and Model 2, which have in common an ancestral relationship of TK and AN groups (Fig. 3a,b), we conclude that the TK and AN groups are descended from a common ancestral population.

Discussion
The present study adds to our previous study of Thai/Lao mtDNA genome sequences by including 22 additional groups from Thailand, including the AA-speaking Mon (MO), ST-speaking Karen, and several TK speaking groups, especially the CT. Similar to our previous mtDNA study [7], genetic heterogeneity among populations belonging to the same ethnic groups was still observed, especially for the Karen (the hill tribes) (Figs. 1, 4 and 5). Geographic isolation and matrilocal residence appear to be important factors influencing the genetic landscape of the highlanders.  Table 1. Each line is the median estimated maternal effective population size (y-axis) through time from the present in years (x-axis) For the remaining lowland groups, i.e., Mon and other TKspeaking groups, gene flow with other groups is a rather more important factor. The Mon, who were a previously dominant group in MSEA located in present-day southern Myanmar and central Thailand since the 6 to 7th century A.D. [38], have been reported to link with Indian populations with some haplogroups, i.e., W3a1b [7]. With data from two additional Mon groups, there is still support for a connection between India and the Mon in the distribution of M subhaplogroups characteristic of South Asia or the Near East [39][40][41], e.g. M6a1a, M30, M40a1, M45a, and I1b (Supplementary  Table S3). Genetic relationship analysis also reveals some Mon populations (MO1, MO5, MO6) clustering with Indian groups, although the other Mon groups were closely related to Thai/Lao populations (Supplementary Figure S5), possibly reflecting gene flow. Thus, based on the many older mtDNA lineages observed, the modern Mon from both Thailand and Myanmar could be an important group for further studies to reconstruct early SEA genetic history.
The Karen in Thailand are refugees who migrated from Myanmar starting from the 18th century A.D. due to the influence of Burmese [42]. However, the ancestors of the Karen probably migrated from some unknown location to Myanmar, as the Karen languages are thought to have originated somewhere in north Asia or in the Yellow River valley in China, i.e., the homeland of ST languages [43]. In agreement with previous studies of different Karen subgroups and/or different genetic markers [31,44,45], we find both northeast and southeast Asian components in the maternal ancestry of the Karen.
The present results emphasize the common maternal ancestry of CT and other TK speaking groups in MSEA, e.g., Laos and Southern China. Demic diffusion is still the most probable scenario for TK-speaking populations (Fig. 2a, Supplementary Table S5), possibly accompanied by some low level of gene flow with autochthonous Mon and Khmer groups. It seems that the prehistoric TK groups migrated from a homeland in south/southeast China to the area of present-day Thailand and Laos, and then split to occupy different regions of Thailand, expanding and developing their own history. During the migration and settlement period, genetic contact with the local AA people was certainly limited, but nonetheless the modeling results, haplogroup profiles and genetic diversity values all suggest some degree of admixture in the CT groups (Supplementary  Table S3, Supplementary Table S6, Table 1). However, historical records indicate that large proportions of the CT groups were taken to neighboring kingdoms as war captives in multiple episodes from 500 to 300 ya [4]. The presentday CT people are probably not solely descended from the prehistorical TK groups, and admixture with local groups might have occurred starting from this time. However, in sum, cultural diffusion did not play a major role in the spread of TK languages in SEA.
Finally, we used simulations to test hypotheses concerning the genetic relationships of groups belonging to different language families. We found that Starosta's model [11] provided the best fit to the mtDNA data; however, Sagart's model [9,10] was also highly supported. These two models both postulate a close linguistic affinity between TK and AN. Although genetic relatedness between TK and AN groups has been previously studied [7,46,47], to our knowledge this is the first study to use demographic simulations to select the best-fitting model. Our results support the genetic relatedness of TK and AN groups, which might reflect a postulated shared ancestry among the proto-Austronesian populations of coastal East Asia [48].
Specifically, the best-fitting model suggests that after separation of the prehistoric TK from AN stocks around 5-6 kya in Southeast China, the TK spread southward throughout MSEA around 1-2 kya by a demic diffusion process, accompanied by population growth but with at most minor admixture with the autochthonous AA groups. Meanwhile, the prehistorical AN ancestors entered Taiwan and dispersed southward throughout ISEA, with these two expansions later meeting in western ISEA. The lack of mtDNA haplogroups associated with the expansion out of Taiwan in our Thai/Lao samples has two possible explanations: either the Out of Taiwan expansion did not reach MSEA (at least, in the area of present-day Thailand and Laos); or, if the prehistoric AN migrated through this area, their mtDNA lineages do not survive in modern Thai/Lao populations. Ancient DNA studies in MSEA would further clarify this issue. Moreover, although mtDNA analyses are informative in elucidating genetic perspectives in geographically and linguistically related populations, they have an obvious limitation in that they only provide insights into the maternal history of populations. Future studies of Y chromosomal and genome-wide data will provide further insights into the genetic history of Thai/Lao populations and the role of factors such as post-marital residence patterns and migration in shaping the genetic structure of the region.

Compliance with ethical standards
Conflict of interest The authors declare that they have no conflict of interest.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this license, visit http://creativecommons. org/licenses/by/4.0/.