Introduction

The indigenous peoples of Southeast Asia (SEA) are often referred to as “negritos” because of their differences in phenotypic appearances, i.e. on average shorter stature, frizzier hair, and darker skin color1,2. The phenotypic resemblances of SEA negritos and African pygmies initially suggested a separate origin in Africa, although convergent evolution now seems more likely (for reviews see Détroit et al.3 and Stock4). Hunting-gathering was the mode of subsistence traditionally practiced by most SEA negritos, whose communities are distributed in three regions: Andaman Islands, Peninsular Thailand and Malaysia, and the Philippines5. Several genetic studies have been conducted in some of these relic populations and indicated deep ancestry6,7,8,9,10,11,12,13 which could link these people with the ancestral groups who migrated from Africa into SEA during the late Pleistocene14. However, the hunter-gatherer groups of Thailand are relatively under-studied compared to those from other SEA regions.

There are two hunter-gatherer groups in Thailand, the Maniq and the Mlabri. The Maniq, also known as Sakai or Kensui, have a census size of ~250 individuals15,16,17 and resemble nearby hunter-gatherer populations in Malaysia, e.g. Jehai and Batek18. For example, the Maniq language belongs to the northern Aslian group of Austroasiatic languages, as do the languages of the Jehai and Batek19, although this most likely reflects a later introduction by agriculturalists20,21. A previous study of sequences of the first hypervariable region (HVR) of the mtDNA control region reported a relatively large divergence of the Sakai (Maniq) from southern Thailand from other Thai groups22. Analyses of mtDNA from prehistoric human remains in southern Thailand suggested that contemporary SEA negritos are descended from pre-Neolithic populations23.

The other hunter-gatherer group, the Mlabri, are characterized by a nomadic traditional lifestyle, in which they construct shelters from the leaves of banana trees that they move out when the banana leaves turn yellow. Their neighbors call them “Phi Tong Luang” which means the “Spirits of the Yellow Leaves”. Although they do not exhibit the typical “negrito” phenotypes, their mode of subsistence is hunting and gathering24,25. This, however, is disputed by Waters (2005)26 who believed that the Mlabri have had contact and trade with neighboring ethnic groups, and thus they are not hunter-gatherer in the strict sense. The conflicting theories as to how they sustain their livelihood make it interesting to study their origin and lineages. With a census size of only ~130 individuals16, the Mlabri live in a mountainous area of northern Thailand and speak a language belonging to the Khmuic branch of Austroasiatic languages. Khmuic languages are also spoken by their neighbors, Htin (of Mal and Pray subgroups) and Khmu, who are also the other hilltribes in northern Thailand. Preceding genetic studies reported extremely reduced genetic diversity in the Mlabri and suggested a cultural reversion from agriculture to hunting and gathering27 and genetic affinities between the Mlabri and the Mal, a Htin subgroup28.

Comparisons of variation in maternally-inherited mtDNA and paternally-inherited, non-recombining portions of the Y chromosome (NRY) have been useful in illuminating differences in the female vs. male demographic history of human populations29,30. However, such comparisons are complicated by the different molecular methods that have been typically employed to assay variation (e.g., HVR sequencing for mtDNA vs. genotyping of SNPs and/or microsatellites for the NRY). Advances in next-generation sequencing make it feasible to obtain mtDNA and NRY sequences that allow directly comparable, unbiased inferences about female vs. male demographic history31,32. Here, we present and analyze complete mtDNA genome sequences and ~2.364 Mbp of NRY sequence from the Maniq and the Mlabri (Fig. 1). These two groups have very distinct mtDNA and NRY lineages and contrasting patterns of variation. Overall, our results support an affiliation of the Maniq with other SEA negrito groups, and a recent origin of the Mlabri lineages.

Figure 1
figure 1

Map of the Mlabri and Maniq sampling locations (stars) and the Thai populations selected for NRY comparisons (comparative mtDNA data are from previous studies33,34). Further details concerning the populations are provided in Tables S1 and S2. There are four Maniq collection sites; Pabon District (site 1) and Kong-Ra District (site 2), Phatthalung Province; Palien District, Trang Province (site 3); and La-Ngu District, Satul Province (site 4). (Adobe Illustrator CS4 14.0.0. http://www.adobe.com/sea/).

Full size image

Results

Mitochondrial DNA

We generated 29 complete mtDNA sequences (18 from Mlabri and 11 from Maniq) with mean coverages ranging from 277× to 1342×. Extraordinarily low mtDNA diversity was observed in the Mlabri (haplotype diversity (h) = 0.29 ± 0.12; mean number of pairwise differences (MPD) = 0.29 ± 0.32; nucleotide diversity (π) = 0.00002 ± 0.00002; segregating site (S) = 1) with just one polymorphic site (nt 43) defining two haplotypes. With 5 haplotypes defined, the Maniq showed higher diversity values (h = 0.82 ± 0.08; MPD = 29.85 ± 14.16; π = 0.0018 ± 0.00097; S = 64) than the Mlabri. However, both of them have lower genetic diversity values than other populations from Thailand (Figure S1)33,34.

The two different mtDNA sequences found in the Mlabri belong to haplogroup B5a1b1 and have a coalescent age of ~735.48 ya (95% highest posterior density (HPD): 1.71–2321.88; Figure S2). The overall coalescent age of B5a1b1 is ~13.99 kya, which was similar to our previous study34. The network analysis (Fig. 2) also supports a single maternal ancestor for the Mlabri. The average number of mutations among the 18 Mlabri sequences is 0.29; with a mutation rate of one mutation in every 3,624 years35, this suggests an age for the Mlabri mtDNA variation of ~1,000 years, close to the coalescence age estimate. Overall, the mtDNA results are in keeping with previous analyses suggesting a single maternal ancestor for the Mlabri and a recent founding age within the past 1000 years27.

Figure 2
figure 2

Networks of mtDNA haplogroups: B5a1b1 is found in Mlabri and R21, M21a and M17a are found in Maniq.

Full size image

By contrast, there are three basal haplogroups (M17a, M21a and R21) observed in the Maniq. When combined with sequences from other populations, the respective estimated ages of M17a, M21a and R21 are ~23.52, ~18.44 and ~8.83 kya (Figure S3). The ages of M17a and M21a are much older than our previous estimations in a Thai dataset34 while R21 is newly-reported in populations from Thailand. The networks and maximum clade credibility (MCC) trees of M21a and R21 indicate a close relatedness of the Maniq sequences with those from the Jehai (an aboriginal Malaysian negrito group), while the Maniq M17a sequences cluster with sequences from the Laotian and Austroasiatic-speaking Blang (Fig. 2; Figure S3).

The multi-dimensional scaling (MDS) plot based on Φ st distances indicates that the studied groups have large genetic distances from each other and from most other populations (Fig. 3a). In the Neighbor Joining (NJ) tree (Fig. 3b), the Mlabri do not show any close relationship with other SEA hunter-gatherer or negrito groups, but instead are closest to other northeastern Thai populations (Seak, Soa and Bru). In contrast, the Maniq cluster with populations from West Malaysia, and are closest to the Jehai of Malaysia (Fig. 3b).

Figure 3
figure 3

The MDS plot (a) and NJ tree (b) based on the mtDNA Φ st distance matrix for 142 populations from Asia, New Zealand and Australia. Population details are provided in Table S1. TN1 is Htin subgroup Mal and TN2 and TN3 are Htin subgroup Pray. KA, SO, SK and BU are Khmu, So, Seak and Bru, respectively. Black indicates Thai/Lao populations; red indicates populations from Indonesia and Malaysia; green indicates populations from China.

Full size image

Y chromosome

We generated 14 sequences (10 from the Mlabri and 4 from the Maniq) for ~2.364 Mb of the NRY, with mean coverages ranging from 10× to 25×. In contrast to the mtDNA results, the Mlabri have higher NRY diversity (six different sequences; h = 0.7778 ± 0.1374; MPD = 17.46 ± 8.94; π = 7.39 × 10−6; S = 37) than the Maniq, which show just one haplotype (i.e., all four sequences were identical). It should be noted that these four unrelated sequences were collected from three different locations (sites 1, 3 and 4) (Fig. 1).

Among six Y chromosomal haplotypes in ten Mlabri males, an MCC tree shows two clades: one comprises five haplotypes of haplogroup O1b1a1a1b which diverged ~1.43 kya and the other is a single haplotype, found in five individuals, classified to haplogroup O1b1a1a1b1a1 (Figure S4), suggesting two paternal lineages for the Mlabri. When we combine the Y sequences of haplogroup O1b1a1a1b and its sublineages retrieved from previous studies (Table S2), the newly reported age of haplogroup O1b1a1a1b is ~6.89 kya with a signal of population expansion. The MCC tree shows that the O1b1a1a1b sequences of the Mlabri are clustered with sequences from the Htin subgroup Mal, while the O1b1a1a1b1a1 sequences are clustered with sequences from the Khmu and the Htin subgroup Pray (Figure S4). The network of haplogroup O1b1a1a1b shows an overall star-like structure, suggesting lineage expansion, and further supports two paternal lineages for the Mlabri (Fig. 4). The NRY sequence of the four unrelated Maniq males belongs to haplogroup K, and when compared to the two haplogroup K2b1 sequences32, the estimated coalescent age is ~41.10 kya (Figure S5). The Maniq sequence thus adds to a small set of basal Y chromosomes sequences.

Figure 4
figure 4

Network of NRY haplogroup O1a1a1b sequences.

Full size image

The Mlabri NRY sequences cluster with other northern Thai populations in the MDS plot and NJ tree (Fig. 5a,b), particularly with their linguistic relatives, Khmu and the Htin subgroup Pray. Although the Maniq are highly diverged from other populations due to their very low diversity, they show closest affinity with Papuan, island SEA (ISEA), and South Asian populations (Fig. 5b).

Figure 5
figure 5

The MDS plot (a) and NJ tree (b) based on the Y chromosomal Φ st distance matrix for 23 populations from Asia and Oceania. Black circles indicates Thai populations. Population details are provided in Table S2.

Full size image

Discussion

We here report and analyze sequences of mtDNA genomes and ~2.364 Mbp of the NRY from the Mlabri and Maniq, two enigmatic hunter-gatherer groups from Thailand. The large genetic divergence of the Mlabri and the Maniq from each other and from the other populations probably reflects genetic drift associated with isolation and very small population sizes. Notably, we observed contrasting patterns of variation in mtDNA vs. the NRY: mtDNA diversity is higher in the Maniq and lower in the Mlabri, while the reverse is true for the NRY.

The SEA negritos are of interest because they are characterized by basal genetic lineages that likely reflect the earliest dispersal of modern humans from Africa6,36. However, there is also evidence that at least some negrito groups have interacted and mixed with non-negrito groups10,11,12,37. Among the three basal mtDNA haplogroups (M21a, M17a and R21) observed in the Maniq (Fig. 2), two (M21a and R21) have been previously reported in negrito groups from Malaysia9,13. MtDNA lineages of the Maniq (also known as the Kensui) were mentioned as unpublished data in a previous study13, which also found M21a and R21 to be common lineages, along with several minor SEA-specific haplogroups (B5a, F1a1a, M7c3c and N9a6a), indicating some SEA influences. In contrast, we did not detect any SEA-specific mtDNA haplogroups (B, F and M7) in the Maniq, but instead we found M17a, which has been reported at low frequency in MSEA34,38 and ISEA39. Moreover, we detected Y chromosomal haplogroup K in the Maniq, which occurs elsewhere in populations of Papua New Guinea40,41, the Philippines30 and Australia42. Haplogroup K is an ancient lineages that is dated here to ~41.10 kya, which is younger than previous estimates in Aboriginal Australians (~48.40 kya)42. Overall, the Maniq display basal genetic lineages and closest genetic affinities with indigenous groups in Malaysia or ISEA, supporting a common ancestry of the Maniq with other SEA negrito communities. The reduced genetic diversity (and population size) of the Maniq probably occurred as a consequence of agricultural expansions to this region. The archaeological evidence in southern Thailand suggests that the Neolithic period of the Malay Peninsula started around 4 kya20 and after the expansion of agriculturalists who also brought Austroasiatic languages to this region, the Proto-Aslian groups were fragmented21, probably leading to gradual isolation.

The Mlabri present a different picture. Their mtDNA and NRY haplotypes are closely related to those of neighboring agricultural groups, rather than belonging to basal haplogroups that are characteristic of typical hunter-gatherer groups. Overall, the mtDNA and NRY results suggest a severe founder event in the Mlabri ancestors, involving just one mtDNA lineage and two NRY lineages that occurred ~1000 ya. These results are in good agreement with previous studies that suggested that the Mlabri originated from an agricultural group and reverted to hunting-gathering27. We further attempted to identify the closest living relatives of the Mlabri, as these might be their source population. We found close genetic relationships in the NRY lineages between the Mlabri and various Htin subgroups as well as the Khmu, in agreement with an oral tradition, which indicates that the Htin are the ancestors of the Mlabri27. A shared common ancestor of these Khmuic-speaking groups might be another explanation. However, the Mlabri mtDNA lineage is most closely related to the Seak, Soa and Bru groups from northeastern Thailand (Fig. 3). Historical evidence indicates that the Mlabri, Soa, Seak and Bru in Thailand migrated from the present-day area of Laos about 100–200 ya25. Although each group migrated independently, their close relatedness might reflect maternal gene flow among various groups in Laos. However, the Mlabri stand out among these groups in apparently having just one maternal ancestor (Fig. 2).

A limitation of our study is that the sample sizes are quite small, which reflects the difficulty in contacting these elusive and enigmatic groups. Nonetheless, a rigorous participant recruitment procedure, and overall agreement with previous studies, suggests that our results are indeed representative of these small populations. Notably, the Mlabri are probably not indigenous SEA hunter-gatherers, but rather originated recently from an agricultural group via an extreme founder event, reflected in their one maternal and two paternal ancestral lineages. The Htin and Khmu are their potential ancestors. By contrast, the Maniq are likely to be indigenous SEA hunter-gatherers whose ancient ancestry is not mixed with non-negrito populations; further study of additional Maniq samples from the other bands as well as autosomal data would provide more information on the history and relationships of this group.

Methods

Samples

Genomic DNA samples from 18 Mlabri who live in a village located in Wieng Sa District, Nan Province, Northern Thailand were obtained from a previous study43. Because the census size of the Mlabri is about 130 individuals, we collected all possible unrelated samples which could be assumed to be representative of this group. The Maniq still lead a highly mobile way of life in the forested Banthad Mountain of three southern Thai Provinces, i.e. Patthalung, Trang and Satun. They live in about 9 small bands or groups, with 8–50 individuals per group, and a total census size of about 250 individuals. All groups are quite closely related, and there is constant movement of family members between groups15,17,18. Written informed consent was obtained prior to the collection of Maniq saliva samples, and all volunteers were interviewed by interpreters and coordinators to screen for subjects unrelated for at least two generations. During the research, we protect the rights of participants and their identity and we confirmed that all experiments were performed in accordance with relevant guidelines and regulations based on the experimental protocol on human subjects which was approved by the Khon Kaen University Ethic Committee (Reference number: HE582391) and the Ethics Commission of the University of Leipzig Medical Faculty. We obtained saliva samples from 11 Maniq from four different bands in three provinces, i.e. Pabon District (site 1: 3 samples) and Kong-Ra District (site 2: 1 sample), Phatthalung Province; Palien District, Trang Province (site 3: 2 samples); and La-Ngu District, Satun Province (site 4: 5 samples) (Fig. 1). We extracted DNA using the QIAamp DNA Midi Kit (Qiagen, Germany) according to manufacturer’s directions.

Sequencing

mtDNA

Twenty-nine complete mtDNA sequences (18 Mlabri and 11 Maniq) were generated from genomic libraries prepared with double indices and enriched for mtDNA as described previously44,45; the libraries were sequenced on the Illumina Hiseq 2500. MtDNA consensus sequences were obtained as described by Arias-Alvis et al.46 with minor modifications. We used Bustard for Illumina standard base calling with a read length of 76 bp. We manually checked sequences with Bioedit (www.mbio.ncsu.edu/BioEdit/bioedit.html) and aligned the sequences to the Reconstructed Sapiens Reference Sequence (RSRS)47 using MAFFT 7.27148.

Y chromosome

We designed a capture probe set for the NRY following a previous study49 and using genome sequences incluldedin SGDP Panel A50 and Panel B51. Sites were kept that were included in the UCSC table browser CRG Align 100 mappability track of the human reference genome hg19 (9,655,499 bp). In order to obtain confident Y haplogroup calls later, we selected regions with at least one ISOGG SNP (version 8.78) and then applied the following filter criterion: a region was removed when its length was lower than 500 bp and its repeat count of 24-mers was higher than 10 in an average window of 120 bp using the CRG Align 24 mappability track from the UCSC table browser. A total of 2,264 regions with a combined length of 2,364,049 bp were included into the capture probe set, which included 2,071 ISOGG SNPs (version 11.144).

Genomic libraries were prepared for each sample using a double indices scheme50. We enriched the libraries for the NRY regions via in-solution hybridization-capture51 using our probe set and the Agilent Sure Select system (Agilent, CA, USA) and sequenced on the Illumina HiSeq. 2500 platform, generating paired-end reads of 125 bp length. Standard Illumina base-calling was performed using Bustard. We trimmed Illumina adapters and merged completely overlapping paired sequences using leeHOM52 and de-multiplexed the pooled sequencing data using deML53.

The sequencing data were aligned against the human reference genome hg19 using BWA’s aln algorithm54 and all sequence pairs that aligned to the NRY region which are previously defined were retained49. Duplicate reads were removed using PicardTools MarkDuplicates, indel realignment was performed using GATK IndelRealigner55 and base quality were re-calibrated using GATK BaseRecalibrator. We identified single nucleotide variants (SNVs) using GATK UnifiedGenotyper v3.3-0 across all samples simultaneously setting the parameter ploidy to 1 and using DBSNP build 138 as a prior position list. The identified SNVs were further filtered as previously described56. We imputed all samples with missing genotype information at any of these variant sites using BEAGLE57.

Statistical analyses

Haplogroup assignment was performed by Haplogrep58 for mtDNA and yHaplo59 for the NRY. Summary statistics for the studied populations and genetic distances (Φ st ) among the studied and compared populations were computed by Arlequin 3.5.1.360. Comparative data for the mtDNA genome sequences were selected from the literature encompassing East Asia33,34,38,61,62,63,64,65,66, Island Southeast Asia12,13,36,67,68,69, Taiwan70, India71, Andaman Islands6,72, Australia73,74,75 and New Zealand76. Further information on the populations used for mtDNA comparison are listed in Table S1. Overlapping NRY sequences of relevant populations from Thailand (Htin, Khmu, Soa, Bru and Seak) (unpublished data), East Asia, Oceania, South Asia and also Indian groups from the United Kingdom and the United States of America32,77,78 were included for comparison; additional details are provided in Table S2.

The Φ st distance matrices of both mtDNA and the NRY were used to construct MDS plots and NJ trees79 by STATISTICA 13.0 (StatSoft, Inc., USA) and MEGA 7.080, respectively. Median-joining networks81 of haplogroups without pre- and post-processing steps were performed by Network (www.fluxus-engineering.com) and visualized in Network publisher 1.3.0.0.

BEAST 1.8.0 was used to construct MCC trees by haplogroup, based on Bayesian Markov Chain Monte Carlo (MCMC) analyses82. We ran jModel test 2.1.783 for selecting the suitable models in each run during the creation of the input file of BEAST via BEAUTi v1.8.082. For mtDNA, the BEAST calculation was done with the data partitioned between coding and noncoding regions with mutation rates of 1.708 × 10−8 and 9.883 × 10−8, respectively35. For the NRY, we used a mutation rate of 8.71 × 10−10 84, and the BEAST input files were modified by an in-house script to add in the invariant sites found in our dataset. After BEAST runs, we used Tracer 1.6 to check the results and TreeAnnotator v1.6.2 and FigTree v 1.4.0 to assemble and draw the Bayesian MCC trees, respectively. The Bayesian MCMC estimates (BE) and 95% HPD intervals of haplogroup coalescent times were calculated using the RSRS for rooting the mtDNA tree and African sequences belonging to haplogroup B2 (GS000035246-ASM) for rooting the NRY tree32.

Data availability

The mtDNA dataset generated during the current study are available in GenBank (https://www.ncbi.nlm.nih.gov/genbank/) with accession numbers MG672483-MG672511 and the Y chromosomal dataset are available from the corresponding author on request.