Genetic evidence supports linguistic affinity of Mlabri - a hunter-gatherer group in Thailand

Background The Mlabri are a group of nomadic hunter-gatherers inhabiting the rural highlands of Thailand. Little is known about the origins of the Mlabri and linguistic evidence suggests that the present-day Mlabri language most likely arose from Tin, a Khmuic language in the Austro-Asiatic language family. This study aims to examine whether the genetic affinity of the Mlabri is consistent with this linguistic relationship, and to further explore the origins of this enigmatic population. Results We conducted a genome-wide analysis of genetic variation using more than fifty thousand single nucleotide polymorphisms (SNPs) typed in thirteen population samples from Thailand, including the Mlabri, Htin and neighboring populations of the Northern Highlands, speaking Austro-Asiatic, Tai-Kadai and Hmong-Mien languages. The Mlabri population showed higher LD and lower haplotype diversity when compared with its neighboring populations. Both model-free and Bayesian model-based clustering analyses indicated a close genetic relationship between the Mlabri and the Htin, a group speaking a Tin language. Conclusion Our results strongly suggested that the Mlabri share more recent common ancestry with the Htin. We thus provided, to our knowledge, the first genetic evidence that supports the linguistic affinity of Mlabri, and this association between linguistic and genetic classifications could reflect the same past population processes.


Background
The Mlabri are a hill tribe in northern Thailand, inhabiting a dispersed area along the border with Laos [1,2]. Today, they are a small population of nomadic huntergatherers, unusual in a region of almost entirely agricultural economies [3]. The modern population size is estimated at around 300 individuals, with some estimates being as low as 100 [4]. The name Mlabri is a Thai/Lao alteration of the word Mrabri, which appears to derive from a Khmuic term for "people of the forest" -in Khmu, mra means "person" and bri "forest". They are also known locally as Phi Tong Luang or "spirits of the yellow leaves", apparently because they abandon their shelters when the leaves begin to turn yellow with the onset of the dry season.
Little is known about the origins of the Mlabri and most evidence comes from linguistic studies. The Mlabri language is classified as a Khmuic language, a subgroup of the Mon-Khmer language in the Austro-Asiatic language family [5]. The available linguistic evidence suggests that the present-day Mlabri language most likely arose from Tin, a Khmuic language [2,6]. However, so far there is no genetic evidence supporting this idea. A recent study suggested Mlabri was founded recently from an agricultural group, thus representing a typical example of cultural reversion [7]. This work, although very interesting, was criticized for not including any of populations neighboring the Mlabri, such as the Htin, Hmong, and northern Thai. As a result, these authors were unable to demonstrate any similarities in the genetic and linguistic affinity of the Mlabri, and so made little comment on the possible source population (s) from which the Mlabri originated [8].
In this study, we analyzed populations samples from throughout northern Thailand, including the Mlabri as well as several neighboring groups, including the Htin, 1 Hmong, Yao, and other populations speaking Austro-Asiatic and Tai-Kadai languages. Four HapMap population samples, representing Altaic, Sino-Tibetan, Indo-European and Niger-Congo language speakers, were also included in this study. We conducted a genome-wide analysis on these samples using 50K SNPs, to investigate the genetic affinity of the Mlabri, examine the concordance of genetic and linguistic affinities, and further explore probable origin(s) of this enigmatic hunter-gatherer group.

Genetic Characteristics of Mlabri
Since this is the first genome-wide genetic study of this enigmatic population, we calculated several population genetic parameters, including SNP diversity, haplotype diversity and linkage disequilibrium (LD).

Reduced genetic diversity in the Mlabri
Expected heterozygosity for SNPs (HS e ) were calculated based on allele frequencies of 55,561 SNPs and the results were shown in Figure 1A. The HS e in Mlabri (0.197) is lower than that of any of other populations in which HS e is at least 0.250 (HM). The expected heterozygosity for haplotypes (HH e ) were calculated based on haplotypes in 500-kb genomic regions (Methods) and the results are shown in Figure 1B. The HH e in Mlabri (0.666) is also much lower than that of any of other populations in which HH e is at least 0.820 (TN). The HH e comparison obtained from larger size of genomic  regions (1 Mb) show the similar results (see Additional file 1, Figure S1). All the above comparisons are statistically significant (t-test, p < 10 -5 ). We also compared genetic diversity among populations using the cumulative proportion of the genome given the number of haplotypes (see Methods). The number of haplotypes was estimated for two different window sizes (500-kb or 1-Mb) respectively, with adjustment for sample size difference among populations (see Methods). Again, we found that the genetic diversity was significantly lower in Mlabri than in other populations for both 500-kb segments ( Figure 2A) and 1-Mb segments ( Figure 2B), respectively. For example, in Mlabri, 99% of the 500 kb segments across the genome carry 17 or less haplotypes in Mlabri, and it is much larger than those in other East Asian populations (52%~68%), CEU (48%), and YRI (20%).

Increased linkage disequilibrium in the Mlabri
The significantly reduced genetic diversity in Mlabri was also reflected by its extent of linkage disequilibrium (LD). We assessed the extent of LD among markers with minor allele frequency (MAF) ≥ 0.05 ( Figure 3A, B) and ≥ 0.1 ( Figure 3C, D). The LD extended substantially longer in Mlabri than all the other populations, measured as the fraction of SNP pairs with r 2 ≥ 0.5 ( Figure  3A, C) or r 2 ≥ 0.8 ( Figure 3B, D). For marker pairs with moderate LD (r 2 ≥ 0.5), we observed this fraction to be 1.6-to 12.3-fold higher in Mlabri than in all the other Asian populations for the distance range above 10-kb to 200-kb. For those marker pairs with strong LD (r 2 ≥ 0.8), the fraction in Mlabri is 2.2-to 31.3-fold higher in Mlabri than in all the other Asian populations, and 6to 259-fold higher than in YRI. Furthermore, LD of r 2 ≥ 0.8 extended more than 1 Mb in Mlabri, whereas in all the other populations, such strong LD extended only up to 200 kb.

Genetic Affinity of Mlabri
The genetic characteristics obtained from above analysis, such as significantly increased LD and extremely reduced haplotype diversity are both consistent with the view from a previous study [7] that the Mlabri were recently founded from a very small number of individuals. The available linguistic evidence suggests that the present-day Mlabri language arose from a Khmuic language, most likely Tin [2,6,7]. To search for the group that gave rise to the founders of Mlabri and to examine if the genetic affinity is consistent with linguistic affinity, we further investigated the genetic relationship of Mlabri and other populations. The rational is that the group with closest genetic relationship with Mlabri, if also consists with linguistic relationship, is most likely the genetic and linguistic founder source.

Individual-based clustering analysis
We first studied the clustering relationships among 446 individuals representing 13 populations in Thailand and the CHB and JPT from the HapMap project (YRI and CEU samples were not included in this analysis). We used an allele sharing distance (ASD) [9] as the genetic distance between individuals and reconstructed an individual tree ( Figure 4) using the Neighbor-Joining algorithm [10]. There are several clear clusters on the tree which coincide with individual linguistic or ethnic affiliations, for example, as denoted in Figure 4, JPT, CHB, Hmong-Mien, Tai-Kadai, Austro-Asiatic, Htin and Mlabri. Notably, all the Mlabri and Htin individuals cluster together tightly (100 per cent bootstrap) although there is a bifurcation between clusters of Mlabri and Htin, indicating that the Mlabri have a closer relationship with the Htin than any of other populations studied.
The above clustering relationships among individuals were also confirmed by principal components analysis (PCA) at the individual level [11]. As shown in a 2dimentional plot of first two PCs ( Figure 5A), individuals tend to cluster with other members of their linguistic or ethnic affiliations. Again, Mlabri showed a closer relationship with the Htin for PC1, which explains 21.8% of variation represented by the first ten PCs. The closer relationship between Mlabri and Htin is even more pronounced in the 2-dimentional plot of PC1 and PC3 ( Figure 5B).
Since the closer relationship between Mlabri and Htin could be due to recent gene flow from Htin to Mlabri or vice versa, we further performed Bayesian cluster analysis as implemented in the STRUCTURE algorithm [12] to examine the ancestry of each person. This analysis considers each person's genome as having originated from K ancestral, but unobserved, populations whose contributions are described by K coefficients that sum to 1 for each individual [13]. Individuals are posited to derive from an arbitrary number of ancestral populations, denoted by K. We ran STRUCTURE from K = 2 to K = 18, with results at K = 8 showing the greatest posterior probability (see Additional file 2, Figure S2). Estimated individual membership fractions in K genetic clusters are shown in Figure 6A. At K = 3, the three clusters correspond with Asian, European and African ancestry, respectively. At K = 4, the new cluster corresponds to a Mlabri specific component, which is exclusively shared by all Mlabri individuals with 100 percent membership fractions and this pattern persisted for all choices of K>3. Similar analyses were also performed using the program frappe [14] which implements a maximum likelihood method. The results obtained from frappe ( Figure 6B) showed a general concordance with that of STRUCTURE; but slight differences were also observed, such as the order with which new clusters emerge at K = 5 and K = 6, and the estimated individual membership fractions for all K>3. Notably, both analyses showed that all Asian populations shared some proportion of the major Mlabri component at K = 4 and K = 5. However, this sharing pattern, unless it is an artifact, is more likely to be explained by shared common ancestry rather than recent gene flow, because it appears highly unlikely that the Mlabri received (or contributed) nearly identical amounts of gene flow from (or to) all Asian populations, and with similar proportions, in every instance. Therefore, the close relationship between Mlabri and Htin is most likely the result of a considerable degree of common ancestry.

Population-and component-based clustering analyses
Because the analyses discussed above were all consistent in showing that individuals from the same population cluster together, it is meaningful to evaluate the genetic relationships among populations. A maximum likelihood tree of populations [15], based on 55,561 SNPs showed that Mlabri (MA) and Htin (TN) have the closest relationship, and this topology was supported by 100% of bootstrap replicates ( Figure 7A). However, the Htin showed signs of admixture in both STRUCTURE and frappe analyses ( Figure 6A, B). This raised the concern that whether the close relationship between Mlabri and Htin was confounded by external immigrants from other populations, given that about half of components of Htin are also found in both Austro-Asiatic and Tai-Kadai populations at Ks>4 in STRUCTURE results ( Figure 6A). We thus further investigated this potential confounding effect by reconstructing the phylogenetic relationships of those clusters inferred from STRUCTURE and frappe (referred to as the "component tree"). The rationale is that the component tree, given the statistical independence of the components, should reveal an evolutionary history that is  Table 1 and Figure 1. Left plot: population structure inferred by STRUCTURE; right plot: population structure inferred by frappe. For both STRUCTURE and frappe results, the figure shown for a given K is based on the highest probability run of ten runs at that K.
less perturbed by recent gene flow and admixture than is a population phylogeny. At K = 8, both STRUCTURE and frappe identified a cluster predominant in the Htin, and with each of the other seven clusters easily associated with a predominant linguistic or ethnic group. We therefore refer to the eight clusters (or components) by their representative linguistic or ethnic group as follows: Altaic/Sino-Tibetan, Hmong-Mien, Tai-Kadai, Austro-Asiatic, Mlabri, Htin, European and African. The component tree was reconstructed based on allele frequencies in each cluster inferred from the STRUCTURE analysis ( Figure 7B). We found that the Mlabri specific and Htin specific component clustered tightly on the tree (supported by 100% of bootstrap replicates), strongly indicating once again that the Mlabri share a more recent ancestry with the Htin than with any other group in our sample.

Discussion
In this study, we analyzed genome-wide SNP data on the Mlabri, as well as several neighboring populations and HapMap population samples. The Mlabri population shows several substantial differences from the other populations: significantly increased LD, extremely reduced haplotype diversity and small effective population size (29), all of which are consistent with the view that the Mlabri were recently founded from a very small number of individuals of an agricultural group but subsequently adopted their current hunting-gathering lifestyle, as proposed by a recent study based primarily on mtDNA and Y chromosome data [7]. Although an alternative scenario could also explain the above genetic characteristics of Mlabri, i.e. the Mlabri are an ancient hunter-gatherer group and maintain their hunting-gathering lifestyle from the very beginning but experienced a severe bottleneck event in the history, the results from the clustering analyses do not favor this scenario. If the Mlabri are an ancient hunter-gatherer group, we expect Mlabri is outside of the clade of all Asian populations and close to the root of Asian clade, but Mlabri is actually inside of Asian clade with Austro-Asiatic group outside on both population tree ( Figure 7A) and component trees ( Figure 7B, C) where no signal of admixture was found have disturbed tree topology.
Both model-free and model-based clustering analyses strongly suggest that the Mlabri share a degree of common ancestry with the Htin, a group speaking Tin language. In this caseas is the general rule in many human populationsthe genetic affinity of these populations is consistent with its linguistic affinity. This result, to our knowledge, is the first genetic evidence supporting the linguistic affinity of the Mlabri and Tin languages. Cavalli-Sforza and colleagues showed an apparent congruence between linguistic phyla and genetic clusters, and they proposed that this congruence indicates "considerable parallelism between genetic and linguistic evolution" [16]. Subsequent studies using diverse scales and methodologies have found variable degrees of association between linguistic and genetic classifications [17][18][19][20][21][22].
Some typical examples of exceptions are populations with language replacement [23][24][25][26] or recent admixture between divergent populations [27,28]. However, human genetic and linguistic diversity have been proposed to be generally correlated, either through a direct link, whereby linguistic and genetic affiliations reflect the same past population processes, or an indirect one, where the evolution of the two types of diversity is independent but conditioned by the same geographic factors [29].
Hunting and gathering was presumably the subsistence strategy employed by human societies for more than two million years, until the end of the Mesolithic period. Contemporary hunter-gatherer groups are often thought to serve as models of an ancient lifestyle that was typical of human populations prior to the development of agriculture. However, there has been complex interaction between hunter-gatherers and non-huntergatherers for millennia. There are contemporary huntergatherer peoples who, after contact with other societies, continue their ways of life with very little external influence. There are also contemporary groups usually identified as hunter-gatherers do not have a continuous history of hunting and gathering, and in many cases their ancestors were agriculturalists and/or pastoralists who were pushed into marginal areas as a result of migrations, economic exploitation, and/or violent conflict [30]. Our current data are not sufficient to distinguish the two scenarios, but in case cultural reversion occurred in the history of Mlabri, the Htin is most likely the source population from which the Mlabri genetically originated. The Htin samples in this study speak Mal language, represent only one of the two varieties (Mal and Prai) of Tin language [31,32], it is possible to further determine which variety the Mlabri language originated from by comparing the genetic relationships between the Mlabri and populations speak the two Tin varieties, although such evidence is indirect and would only make sense when the assumption hold that the genetic origin of the Mlabri was not earlier than the divergence of the two language varieties and there was no language replacement.

Conclusions
In summary, our results strongly suggested that the Mlabri share more recent common ancestry with the Htin, a group speaking a Tin language. This result, to our knowledge, is the first genetic evidence supporting the linguistic affinity of the Mlabri and Tin languages. We proposed that Htin is most likely the source population from which the Mlabri genetically originated in case cultural reversion occurred in the history of Mlabri.

Populations and Samples
Samples from the Mlabri as well as other 12 populations were collected in Thailand. The sample information, including sample size, ethnic and linguistic information is shown in Table 1, and the sampling locations are shown in Additional file 3, Figure S3. These samples were also described previously [33]. In this study, eight Mlabri samples were not included because they were identified as close relatives (IBD > 0.2) of one of the rest samples. Four population samples (60 YRI, Yoruba from Ibadan, Nigeria; 60 CEU, Utah residents with ancestry from northern and western Europe; 45 CHB, Han Chinese in Beijing; and 44 JPT, Japanese in Tokyo) obtained from the database of the International HapMap Project [34] were also included in this study.

Data Sets
Genotype data of 13 Thailand population samples generated using Affymetrix Genechip Human Mapping 50K Xba array were obtained from the Pan-Asian SNP Initiative [33]. Detailed information about data filtration and data quality control was described elsewhere [33].

Statistical analysis Haplotype inference
Haplotypes of 22 autosomes were inferred for each individual from its genotypes with fastPHASE [37] version 1.2. "Population labels" were applied during the model fitting procedure to enhance accuracy. The number of haplotype clusters was set to 20, the number of random starts of the EM algorithm (-T) was set to 20, and the number of iterations of EM algorithm (-C) was set to 50. This analysis was used to generate a "best guess" estimate of the true underlying patterns of haplotype structure [37]. We run fastPHASE for 55,561 SNPs shared by 17 populations, and only unrelated individuals were included.

SNP heterozygosity
Heterozygosity for each SNP (HS e ) was calculated based on allele frequencies.

Haplotype heterozygosity
To calculate heterozygosity for haplotypes (HH e ), the genome was divided into 500-kb regions, with each region having roughly 14 SNPs. HH e were calculated for each region using haplotype frequencies [38]. Considering the substantial variation of recombination across human genome [39,40], we adopted a slide window strategy and let the sliding window move 100 kb each time. For each population, HH e were averaged over all windows.

Number of haplotypes and its cumulative proportion of the genome
The number of haplotypes was obtained by counting the number of haplotypes for a given window size, i.e. 500kb or 1-Mb, respectively, for each population. The same sliding-window scheme as mentioned before was employed. Since this measurement could be affected by sample size, we sampled 36 chromosomes (equal to the sample size of Mlabri) without replacement in each population. Note that Mlabri has the smallest sample size in all the populations studied. For a population with sample size larger than 36 chromosomes, the sampling was repeated 100 times for each segment and the average of the number of haplotypes of all replications was taken as the number of haplotypes.
The cumulative proportion given a number of haplotypes was obtained by estimating the proportion of the sliding-windows across the genome carrying equal or less haplotypes.

LD calculation
Linkage disequilibrium (LD) between SNPs were measured using r 2 following Hill and Weir [41] and calculated from haplotype data.

Principal component analysis for individuals
Principal component analysis (PCA) was performed at individual level using EIGENSOFT version 2.0 [42].

Genetic distance for individuals
We used an allele sharing distance (ASD) [9,43] as a measure of genetic distance between individuals and a 454 × 454 inter-individual genetic distance matrix was generated according to genotypes of 55,561 autosomal SNPs.

Tree reconstruction
The tree of individuals was reconstructed based on ASD distance and using Neighbor-Joining algorithm [10] with the Molecular Evolutionary Genetics Analysis software package (MEGA version 4.0) [44]. Trees of populations as well as components were reconstructed using maximum likelihood method [15] with CONTML program in PHYLIP package [45].

STRUCTURE analysis
Ancestry of each person was inferred using a Bayesian cluster analysis as implemented in the STRUCTURE program [12,46]. We ran STRUCTURE from K = 2 to K = 18 and repeated 10 times for each single K. All STRUCTURE runs used 20,000 iterations after a burnin of length 30,000, with the admixture model and assuming that allele frequencies were correlated [46].

frappe analysis
The program frappe [14] implements a maximum likelihood method to infer genetic ancestry of each individual. As in STRUCTURE analysis, this analysis considers each person's genome as having originated from K ancestral, but unobserved, populations whose contributions are described by K coefficients that sum to 1 for each individual [13]. The program was run for 10,000 iterations from K = 2 to 18 and repeated 10 times for each single K.
Additional file 1: Contains Figure S1 -Haplotype heterozygosity (HH e ) in 17 populations. In table at the bottom of each plot displayed the average and the standard deviation of HH e in each population sample. The sample information of each population is shown in Table 1. HH e were calculated from haplotypes of 1-Mb windows, SD denotes standard deviation of the HH e values across windows. Figure S3 -Probability Estimations for the Number of Clusters, with Ten Repeats for Each K. The ordinate shows the Ln probability corresponding to the number of clusters (K) shown on the abscissa. A: showing maximal probability estimation of ten runs at each K (from K = 2 to K = 9); B: showing probability estimation at K = 2 to K = 18 in all ten runs.

Additional file 2: Contains
Additional file 3: Contains Figure S2 -Geographical distribution of Thailand population samples. Red dots on the map indicated sampling locations. Information of population IDs can be found in Table 1.