Prehistoric global migration of vanishing gut microbes with humans

Commensal gut bacteria play a complex part in human biology¹, but the evolutionary history of these microbial relationships remains poorly characterized⁹. Identifying long-term associations between humans and microorganisms is important for understanding their co-evolution^7,10,11, and this could reveal microbial species and functions integral to human health and physiology.

Direct observations of historical microbial relationships have been obtained by sequencing palaeofaeces¹². This study suggested that ancient gut microbiomes were more similar to those of contemporary non-industrialized populations than to those of industrialized ones (the term industrialized here refers to post-industrial populations as well as those in advanced and ongoing stages of industrialization). However, the comparatively recent ages of existing palaeofaeces (around 1,000–2,000 years old), combined with the challenges of analysing degraded ancient DNA, makes it difficult to determine whether these microbial lineages have been faithfully transmitted over long evolutionary timescales (of around 10,000–100,000 years) or were acquired or transmitted more recently.

Comparative genomics methods provide an alternative approach to inferring evolutionary history from the patterns of genetic diversity among contemporary microbial genomes^{8,13,14,15,16,17,18,19}. Building on previous work in non-human primates^2,8,19,20,21, advances in strain-resolved metagenomics have made it possible to use this approach on gut microbial species sampled from human populations worldwide. Previous work has identified phylogenetic correlations between gut microbial genomes and their human populations of origin⁷, potentially consistent with long-term associations^18,22,23. However, interpreting these correlations remains challenging because they require genome comparisons at the subspecies, or strain, level. Conventional phylogenetic methods²⁴ often fail in such settings because intra-species gene transfer and homologous recombination between strains produce mosaic patterns of ancestry that cannot be represented by a single phylogenetic tree^25,26,27,28. These high rates of recombination make it difficult to distinguish ancient co-migration from other sources of population structure, such as geography or host lifestyle variation^22,23. Moreover, existing strain-level analyses have been biased towards gut microbial species in industrialized populations²⁹. These sampling biases have hindered our ability to infer the evolutionary history of the microbial taxa that are the most likely candidates for ancestral co-migration.

To address these challenges, we conducted deep metagenomic sequencing on faecal samples from the Tsimane horticulturalists of Bolivia (previously characterized using 16S ribosomal RNA sequencing⁶) and compared their microbiomes with those of the Hadza hunter-gatherers of Tanzania³. These two cohorts represent unique and geographically distant contemporary populations whose ancestors separated thousands of years ago, providing a rare opportunity for examining human–microorganism co-migration patterns across a wide range of candidate ancestral microbial taxa. We identified more than 1,200 microbial species shared between the two populations, most of which (60%) were rare in or absent from industrialized populations. Using population genetic analyses that explicitly account for microbial recombination, we provide evidence of strain-level reproductive isolation that is stronger than that observed in microorganisms from comparable industrialized populations. These findings suggest the transgenerational maintenance of gut microbial communities throughout ancient human migrations, providing insights into the roles of horizontal gene transfer and host–microorganism co-evolution in the long-term structuring of the human microbiome.

The Tsimane people are Indigenous forager-horticulturalists living in the Bolivian Amazon. The individuals in this cohort live in six villages along the Maniqui River and three villages in a nearby forest region⁶. We performed deep metagenomic sequencing on 133 faecal samples collected from 85 Tsimane adults in 2009 and 2012–2013 (median depth = 31.9 gigabases; Extended Data Fig. 1a). We recovered 12,746 metagenome-assembled genomes (MAGs) representing 1,408 bacterial and archaeal species (Fig. 1a, Methods and Supplementary Table 1). We then compared these genomes with an analogous collection of 32,034 MAGs previously obtained from deeply sequenced microbiome samples of 137 Hadza hunter-gatherers in Tanzania³ (Extended Data Fig. 1b) to quantify how the gut microbial species of these populations have diverged over the past 100,000 years.

Despite their large geographical separation and lack of historical direct contact, 87.4% (1,231) of the species in the Tsimane cohort were also found in Hadza individuals (Fig. 1b). Some of these shared species (for example, Ruminococcus bromii) are also prevalent worldwide, but most (60.2%; 848 species) are rare in or absent from the faecal microbiomes of industrialized populations (Extended Data Figs. 2 and 3 and Supplementary Table 2). Similar results were obtained when controlling for sequencing depth outliers between populations (Extended Data Fig. 2d) or when using a mapping-based analysis to rule out common assembly biases (Extended Data Fig. 3). However, differences in strain-level heterogeneity and/or sequencing depth distributions between cohorts may have contributed to the differential recovery of some rare MAGs. Species sharing between distant populations could result from ancient co-migration or more-recent acquisition from the environment or human contact. We sought to distinguish between these scenarios by examining the genetic variation in each of these shared species and comparing these patterns with expectations from historical human migrations.

As a first step, we calculated the average nucleotide identity (ANI) between all pairs of MAGs in each of the 636 shared species with at least four MAGs each in both the Hadza and Tsimane cohorts (Extended Data Fig. 4). These ANI values provide an upper bound on the time to the most-recent common ancestor for each MAG pair. Estimates of the molecular clock in gut bacteria suggest that mutations accumulate at a rate of at least 10⁻⁷ per site per year (Extended Data Fig. 5 and Supplementary Methods 1). The typical ANI between Tsimane and Hadza MAGs (around 98%; Fig. 1c) therefore suggests that they shared a common genetic ancestor in the past 100,000–200,000 years (Extended Data Fig. 5). The fact that most species present in Tsimane samples share a common ancestor with those in the Hadza in this timeframe (Extended Data Fig. 6) places strong constraints on their previous evolutionary history and suggests that the ancestors of the Tsimane did not acquire a completely new suite of gut microorganisms during or after the population’s migration to the Americas.

Further insights can be obtained by comparing the ANI values between Tsimane–Hadza genome pairs and their corresponding within-population counterparts (Tsimane–Tsimane, Hadza–Hadza). As expected, these data revealed that pairs of conspecific genomes from the same human population were, on average, more similar to each other than to genomes in different populations (P < 0.001, analysis of variance, Fig. 1c). However, the effect sizes were smaller than the variation within each population (mean fixation index (F_ST) = 0.10; Extended Data Fig. 6). The broad overlap between these distributions indicates that most of the strain-level diversity in the Tsimane microbial populations pre-dates the onset of geographical isolation.

Although the average ANI values showed little geographical differentiation, we observed a much stronger trend in the high-ANI tail of the distribution (>99.75%), indicating a marked enrichment of pairs within populations compared with between populations (P < 2.2 × 10⁻¹⁶, Wilcoxon rank-sum test; Fig. 1c, inset). This within-population enrichment suggests that there was restricted strain transmission between these populations in recent evolutionary history, consistent with divergence driven by historical separation rather than ongoing exchange. However, estimating the timing of this separation can be challenging, because these high ANI values are often strongly influenced by recombination with other strains^25,28,30. Previous work has shown that these high-ANI pairs represent nearly clonal strains, which have genomes that have not yet been overwritten by recombination^25,28,31 (Extended Data Fig. 5 and Supplementary Methods 3). Their clonal ancestry can be inferred from the large stretches of nearly identical DNA sequences, which are interspersed with recombined segments from other strains. We estimated the total amount of clonal ancestry for each pair of MAGs using the fraction of core genes with identical DNA sequences (Fig. 2a). By calibrating these estimates with the rate of mutation accumulation in the vertically inherited genome regions, we determined that pairs of MAGs with more than 10% identical genes correspond to clonal strains that shared a recent common ancestor within the past approximately 5,000 years (Fig. 2a,b, Extended Data Fig. 7, Supplementary Table 4 and Supplementary Methods 3). Using this threshold, we identified all such pairs of clonal strains among the 636 shared species in Fig. 1c to investigate the landscape of recent strain transmission between the Tsimane and the Hadza.

Aggregated comparisons across all 636 species revealed that, in total, 78,746 MAG pairs (7.8%) were clonally related by this criterion, implying that they shared a most-recent common ancestor at most around 5,000 years ago (5 ka; Fig. 2b,c). Most of these clonally related pairs (93.8%) were found in the intra-population comparisons, indicating that most recent relationships occurred in the same human population (P < 10⁻⁴, permutation test). Moreover, the few clonal strain pairs shared between the Tsimane and Hadza were unevenly distributed across species, with 545 of the 636 species (85.6%) showing no recent strain sharing between the two populations (Fig. 2d; P < 10⁻⁴, permutation test). These results indicate that, in most shared species, very limited strain sharing has occurred between the Tsimane and Hadza over the past approximately 5,000 years.

However, there were some notable exceptions to this trend of limited inter-population strain sharing. Some of the species had levels of Tsimane–Hadza strain sharing comparable to their within-population baselines, consistent with a more recent spread. These species, which include Akkermansia muciniphila and Bacteroides ovatus (Fig. 2e, bottom), tend to be enriched in industrialized populations. To explore this pattern further, we expanded our analysis to include an additional 8,855 MAGs from Europe, Asia and North America^32,33,34. The results revealed that these recently shared Tsimane–Hadza species also had closely related MAGs distributed globally (Fig. 2c–e and Extended Data Fig. 8). This suggests that species with high Tsimane–Hadza strain-sharing rates have been introduced into the Tsimane and Hadza through more-recent interactions with other human populations.

Notably, the comparison with industrialized MAGs also revealed another class of species that is prevalent across lifestyles but showed little recent strain transmission between the Tsimane and Hadza (Fig. 2e and Extended Data Fig. 8). These taxa, which include common species such as Agathobacter rectalis and Ruminococcus bromii, indicate that Tsimane and Hadza strains can remain isolated despite the global prevalence and recent sharing of these species in industrialized populations.

One limitation of strain-sharing metrics is that the total number of clonal pairs is shaped by several factors, such as host contact networks^35,36 and genome-wide selective sweeps³⁷, that can be difficult to compare across populations. To further corroborate our observation of microbial isolation between the Tsimane and Hadza (Fig. 2), we used an orthogonal approach that included information from the remaining pairs of MAGs (that is, those with <10% identical genes). Although their original clonal ancestry has been overwritten by recombination with other strains, recent recombination events between the ancestors of these more-distantly related strains generate shorter segments of nearly identical DNA against a backdrop of lower genome-wide ANI (95–99%) (Fig. 3a). Subsequent mutation and recombination events gradually degrade these segments over time, creating a relationship between the length of an identical stretch of DNA and the age of the transfer event. The timing of these events can be used to measure the rates of historical microbial gene flow within versus between human populations (Fig. 3a)—the conventional genetic signature indicating that they were part of the same microbial population³⁸.

We quantified this signal by identifying the longest tract of identical synonymous sites between each pair of MAGs for each of the 636 shared species (representative species are shown in Fig. 3b). We excluded clonal pairs (>10% identical genes) to focus on gene flow rather than on DNA acquired through common ancestry (as analysed in Fig. 2). Differences in tract lengths for within-population versus between-population comparisons were most pronounced in the longer tracts, coinciding with recent recombination events (Fig. 3b and Supplementary Methods 4). Thus, to quantify this trend systematically, we calculated the 99th percentile of tract lengths (L99 length) for each of the 636 species.

Consistent with our other results, most species had shorter L99 lengths for comparisons between the Tsimane and Hadza than for comparisons within populations, suggesting genetic isolation in the recent past (Fig. 3c; P < 10−4, Wilcoxon rank-sum test). However, some exceptions—including Escherichia coli—had similar L99 lengths across populations, reflecting their high levels of more recent strain-sharing (Fig. 2). We repeated this analysis for MAGs from Asia, Europe and North America (Fig. 3c). The differences in L99 lengths for within-population versus between-population comparisons were smaller for these industrialized populations than for the Tsimane and Hadza populations (Wilcoxon rank-sum test; Asia–North America, P = 0.01; Europe–North America, P = 0.41). This suggests that there is a stronger barrier to gene flow between the Tsimane and Hadza populations than between the industrialized ones. The same trend was also observed in species that were prevalent across lifestyles (Extended Data Fig. 9, P = 0.018, Wilcoxon rank-sum test), suggesting that it is not driven by systematic differences in the underlying microbial taxa being compared.

Next we used the observed L99 lengths (Fig. 3c) to estimate the timing of genetic isolation between the Tsimane and Hadza microbiomes. For each species, we added simulated mutations to pairs of within-population MAGs until their L99 lengths matched their observed between-population levels. By calibrating these in silico population splits with the estimated mutation rates of gut commensals²⁸, we can estimate the number of years since the gut strains in a given pair of populations last exchanged DNA (Fig. 3d, Supplementary Table 5 and Supplementary Methods 4). These estimates revealed that genetic exchange was considerably more recent between industrialized populations (Europe–North America: median = 1,109 years, s.e.m. = 377 years; Asia–North America: median = 4,637 years, s.e.m. = 1,267 years) than between non-industrialized populations (Tsimane–Hadza: median = 17,090 years, s.e.m. = 863 years, P < 10⁻⁴, two-tailed Mann–Whitney U-test). Notably, the estimated onset of genetic isolation between the Tsimane and Hadza microbiomes roughly aligns with the timeframe between the initial migration of humans out of Africa (around 60–70 ka)^39,40 and the settlement of the Americas (around 16–23 ka)^41,42, considering the uncertainties in our molecular clock estimates (Supplementary Methods 1).

To obtain an additional independent estimate of microbial split times, we used an existing demographic inference method⁴³. Rather than using stretches of identical DNA, this method uses an entirely orthogonal approach, analysing the distribution of single-nucleotide variant frequencies across populations (Fig. 4a,b, Supplementary Methods 5). Although applying this model to microorganisms has challenges, such as small sample sizes and complex population structures (for example, the presence of clades or subspecies), we identified 15 species with sufficient data and an appropriate population structure to estimate isolation times (Supplementary Table 6; see Supplementary Methods 5 for detailed selection criteria). For each species, the inferred models broadly align with those based on L99 lengths (Fig. 3d and Extended Data Fig. 10) and clearly support an ancient population split in microbial demography rather than two populations with no isolation (Fig. 4b,c and Supplementary Methods 5). Together, these results indicate that many gut microbial species have cohabited with humans for at least as long as humans have been migrating worldwide.

Using several complementary analyses, we show that hundreds of commensal species inhabited the human gut for tens of thousands of years. These commensals have been transmitted across generations in Africa and during human migration worldwide. The findings shed light on some of the historical dynamics that have shaped the present-day geographical distribution of gut microbiome species. Many of these species were previously identified as VANISH taxa (volatile and/or associated negatively with industrialized societies of humans) based on their enrichment in Hadza hunter-gatherers^3,44,45. Our genetic analyses demonstrate that many VANISH species were also present in the common ancestors of the Tsimane and Hadza, suggesting that widespread losses might explain the absence of these species from industrialized populations. We also observed examples of BloSSUM species (bloom or selected in societies of urbanization or modernization)⁴ with extensive strain-sharing across industrialized and non-industrialized populations; this probably resulted from species transmission from industrialized populations to the Tsimane and Hadza. Finally, we identified a third class of species (stable independent of lifestyle (StIL)) that is prevalent across lifestyles but has signatures of Tsimane–Hadza isolation consistent with prehistoric co-migration. Notably, in this subset of StIL species, the degree of Tsimane–Hadza isolation was generally stronger than that observed between similarly separated industrialized populations (Fig. 2e and Extended Data Fig. 8). These observations suggest that ongoing exchange in industrialized populations has probably degraded signals of co-migration; this may partly explain why previous observations of these signals have been limited.

Our results highlight that, by focusing on the microbiomes of contemporary populations with limited exposure to industrialization, we can better understand long-term human–microbiome associations that extend to early human history. Although these are contemporary populations, the genetic relationships between their gut bacteria can shed light on their past evolutionary history, similar to the benefits of including diverse populations in human genome studies⁴⁶. As we expand our knowledge of microbiomes in rural and Indigenous populations, it is crucial to adhere to ethical research practices and foster collaboration with local communities, respecting their autonomy throughout the research process while ensuring mutual benefit, fostering open dialogue and building trust. Notably, the loss of VANISH taxa with industrialization, which is linked to an increased risk of several chronic metabolic and inflammatory disorders, raises the possibility that reintroducing these species may help to restore microbial functions that co-evolved with human biology and improve people’s health⁴⁷. This endeavour will require thoughtful and inclusive dialogue between researchers, ethicists and Indigenous communities.

Ethics and inclusion statement

Informed consent was obtained from the Tsimane at three levels: the Gran Consejo Tsimane (Tsimane governing body), community leaders and study participants. All study protocols were approved by the institutional review boards of the University of California Santa Barbara (3-21-0652) and the University of New Mexico (07-157, 15-133, 17-230). In Bolivia, study protocols were approved by the Universidad Mayor de San Simón in Cochabamba, Bolivia, and the Gran Consejo Tsimane. For specimen export, approval was provided by the Gran Consejo Tsimane, the Universidad Mayor de San Simón Medical School (Cochabamba, Bolivia) and the National Center for Tropical Diseases (CENETROP) in Santa Cruz, Bolivia. The binational (Bolivia–US) Tsimane Health and Life History Project (THLHP) team has worked in Bolivia for more than two decades and has trained dozens of Tsimane research assistants employed by the team. All of our research was conducted in collaboration with local Tsimane communities and the Gran Consejo Tsimane. The original impetus for collecting faecal samples was to analyse their microbial composition. However, as a benefit to participants, samples were also analysed by field microscopy for applicable parasitic species that could be treated with antihelminthic or anti-protozoan medications. The aim of the study was explained verbally to the participants. The THLHP has also engaged in capacity building through education (partnering with the One Pencil Project, a US non-governmental organization that provides educational supplies to primary schools and scholarships for Tsimane students pursuing university or medical school). The THLHP has invested in health infrastructure at the municipal hospital, has employed more than a dozen Bolivian physicians over the years, providing direct medical care during more than 65,000 patient visits, and has trained local health officials, built an official medical and dental clinic and provided expertise, personnel and equipment to two local hospitals.

We also analysed previously published metagenomic data³ from stool samples that were collected from Hadza hunter-gatherers in 2013–2014. The samples were legally and ethically collected in consultation with Tanzanian field guides after obtaining informed consent from Hadza participants who were given a translated verbal description of the research project. A material transfer agreement with the National Institute for Medical Research in Tanzania specifies that the collected stool samples will be used solely for academic purposes. Permission for the study was obtained from the National Institute of Medical Research and the Tanzania Commission for Science and Technology (MR/53i 100/83 and NIMR/HQ/R.8a/Vol.IX/1542). The results of a recent study³ were converted into an infographic using the local language and conveyed to the Hadza by an anthropologist and a Tanzanian field guide.

Sample collection

The Tsimane stool samples were part of a previously published study⁶ and, as part of the current work, were subjected to additional sequencing. In brief, the two sets of samples (2009 and 2012–2013) were collected from participants across nine villages under the auspices of the THLHP⁵¹, which collects health and demographic data from Tsimane participants and provides primary health care to their communities. The set of stool samples used in this study was from Tsimane individuals who were at least 3 years old (mean age = 28.4 years, s.d. = 13.2 years). Field specimen collection protocols were devised considering the logistical challenges and were consistent across the 2009 and 2012–2013 cohorts. During data collection, none of the families in the participating villages had plumbing, used pit toilets or nappies, or had access to refrigeration⁶. Faecal samples were collected in sterile urine specimen cups given to the participants the day before sample collection. Researchers visited the participants’ homes between 7 a.m. and 9 a.m. to collect the specimens. Samples were transported to a field laboratory in coolers with reusable ice packs within 1–2 h of collection. Samples were homogenized in the collection cup and then divided into 2-ml sterile cryotubes using non-sterile wooden tongue depressors. Cryotubes were immediately stored in liquid nitrogen and then transferred to −20 °C freezers before being transported to the USA on dry ice. Samples were stored in −80 °C laboratory freezers in the USA until analysis. Sample collection information for the Hadza cohort is provided in the previous study⁴⁵. In brief, at the time of sample collection, when available, we recorded the age (either provided by the individual or estimated based on their recall), gender, weight and related information. After retrieval from the participants, the samples were immediately placed into liquid-nitrogen containers in the field, transported frozen to Arusha and then mailed by air on dry ice to the USA, where they were maintained at −80 °C until further processing.

Library preparation and sequencing

DNA was extracted from stool samples using MoBio PowerSoil kits (Qiagen) following the manufacturer’s instructions. Libraries were prepared using Nextera Flex kits with a target of 10 ng of DNA with 12 bp unique dual-indexed barcodes for 12 cycles to minimize amplification bias. Paired-end sequencing (2 × 140 bp) was performed on a NovaSeq 6000 using S4 flow cells at the UCSD Microbiome & Metabolomics Core. Samples were randomized across runs and sequenced repeatedly until the target depth was reached. The minimum target depth for each sample was 100 million paired-end reads (around 28 gigabases). A total of 5.25 tera-base pairs of metagenomic data were generated.

Metagenome quality control and assembly

Raw sequencing reads were demultiplexed and processed using a custom processing pipeline (https://github.com/MrOlm/nf-genomeresolvedmetagenomics). In brief, adapters were trimmed using FastP⁵². Reads that aligned with the human and PhiX genomes were removed using Bowtie2 (ref. ⁵³). FastQC was used to ensure read quality⁵⁴. Metagenomes were assembled using SPAdes⁵⁵, and reads were mapped against the assemblies using Bowtie2. Contig coverage was determined using CoverM⁵⁶. Assembly size and contig metrics were evaluated using QUAST⁵⁷, and assemblies were filtered to contigs of ≥1,500 bp for all downstream analyses.

Metagenome-assembled genome recovery

Genome binning was performed using METABAT2⁵⁸, and genome bin quality was assessed using CheckM⁵⁹. Genomes with ≥50% completeness and <10% contamination according to CheckM were retained, in accordance with MIMAG standards⁶⁰. Genomes with ≥50% completeness and <10% contamination were classified as medium-quality genome bins, and those with ≥90% completeness and <5% contamination as high-quality genome bins. We next used GTDB-Tk⁶¹ (r220) to determine the taxonomic classification of each medium- and high-quality genome bin. Genome bins derived from publicly available studies^32,33,34 were also re-run using GTDB-Tk r220 to ensure up-to-date taxonomic classification.

Mapping-based prevalence calculations

To determine species presence across samples, quality-filtered metagenomic reads were mapped against the species-level representative genome database using Bowtie2. Alignment files were processed using inStrain quick_profile (v.1.2.14) and CoverM (v.0.4.0) to calculate per-genome coverage statistics. A species was considered present in a given sample if it achieved a breadth of coverage of ≥0.5 (at least 50% of the representative genome was covered by mapped reads). Prevalence for each species was then calculated as the percentage of metagenomes in which it was detected.

Sequence divergence and clonal strain-sharing analyses

To determine the all-versus-all ANI in each GTDB-assigned species, we used dRep^62,63 (v.3.5.1) with the following parameters: (dRep compare --S_algorithm goANI --SkipMash). This command ensures that all genomes in a GTDB-assigned species are compared using the goANI algorithm^64,65. In brief, goANI identifies open reading frames using Prodigal (v.2.6.3)⁶⁶, aligns their nucleotide sequences with NSimScan⁶⁷ and calculates ANI as the average sequence identity of all aligned genes. We used these ANI values to estimate an upper bound on the corresponding divergence times using molecular clock analysis described in Supplementary Methods 1.

To evaluate the statistical significance of the ANI differences within versus between populations, we performed a one-way analysis of variance for each species comparing the ANI distributions of intra- and inter-population strain pairs. We then applied Benjamini–Hochberg corrections to the P values for each species to obtain a false discovery rate-corrected q value for each species. We next used Fisher’s sumlog method to combine these q values across species into a single final P value. We also calculated the fixation index F_ST (a standard measure of population differentiation) using the procedure described in Supplementary Methods 2.

To perform the strain-sharing analysis (Fig. 2), we first calculated the number of identical genes shared by each genome pair. The set of genes aligned between each genome pair was filtered to include only those with at least 500 bp aligned. The percentage of identical genes was then calculated as the number of gene pairs with 100% nucleotide identity divided by the total number of aligned gene pairs with at least 500 bp aligned.

We defined strains as clonally related when they shared more than 10% identical genes. To estimate the divergence time between such clonal pairs, we used a hidden Markov model-based method²⁸ to identify recent recombination events and obtain recombination-corrected divergence time estimates. Further details on this analysis and the corresponding inference procedure are provided in Supplementary Methods 3.

We used a permutation test to determine whether within-population genome pairs are more clonally related than are between-population genome pairs. For each species, we computed the test statistic (p_HH+ p_TT)/2 − p_HT, where p_XX denotes the proportion of genome pairs with >10% identical genes in each population combination. We generated a null distribution by permuting population labels across genomes within each species (n = 10,000 permutations) and calculated empirical P values as (n* + 1)/(n + 1), where n* is the number of permutations with a test statistic of at least the observed value. For species-specific permutation tests, P values were corrected for multiple comparisons using the Benjamini–Hochberg procedure. For the global test, we permuted labels within species and pooled pairs across all species to calculate a single test statistic and P value.

SNV catalogues and gene flow analysis

For each species sufficiently represented in the Hadza and Tsimane, we generated a catalogue of SNVs present in the core genomes of that species (core genomes are defined as the core genes present in >90% of MAGs in a species). To create an SNV catalogue for a given species, we first chose one genome in the species as the representative genome based on the highest score determined using the following formula: Completeness − 5 × Contamination + 3 × log₁₀(N₅₀) + Length/10⁶ (N₅₀ is the contig length at which 50% of the total assembled sequence length is contained in contigs of that length or longer). We then used a previously published pipeline to annotate the sites in the representative genome that had SNVs compared with all of the other genomes in the species⁶⁸ (https://github.com/zjshi/snv_analysis_almeida2019). We next filtered this whole-genome SNV catalogue to include only sites located in core genes. To identify core genes in the pangenome of each species, we called open reading frames for all genomes in the species using Prodigal (v.2.6.3)⁶⁶, clustered all open reading frames in the pangenome at 90% identity using MMseqs2⁶⁹ and determined which gene clusters were present in at least 90% of the genomes. On the basis of these core genes, we created an additional mask that selected fourfold degenerate sites, thus focusing on sites where all mutations would be synonymous. The resulting SNVs were used to perform the gene flow analysis shown in Fig. 3. Further details of these analyses, including the underlying assumptions and motivations, are provided in Supplementary Methods 4.

Demographic inference using Moments

We used the SNV catalogues to perform the demographic inference shown in Fig. 4. For each species, we computed the frequency of each SNV in the Tsimane and Hadza populations to obtain a joint frequency spectrum across all fourfold degenerate sites. We then used the Moments software package⁴³ to fit these data to a demographic model of a historical population split. Further details of these analyses, including the quality-control criteria used to select the species in Fig. 4, are provided in Supplementary Methods 5.

Reporting summary

Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article.

1. Sommer, F. & Bäckhed, F. The gut microbiota — masters of host development and physiology. Nat. Rev. Microbiol. 11, 227–238 (2013).

   Article  CAS  PubMed  Google Scholar 

2. Sanders, J. G. et al. Widespread extinctions of co-diversified primate gut bacterial symbionts from humans. Nat. Microbiol. 8, 1039–1050 (2023).

   Article  CAS  PubMed  PubMed Central  Google Scholar 

3. Carter, M. M. et al. Ultra-deep sequencing of Hadza hunter-gatherers recovers vanishing gut microbes. Cell 186, 3111–3124 (2023).

   Article  CAS  PubMed  PubMed Central  Google Scholar 

4. Sonnenburg, J. L. & Sonnenburg, E. D. Vulnerability of the industrialized microbiota. Science 366, eaaw9255 (2019).

   Article  CAS  PubMed  Google Scholar 

5. Blaser, M. J. The theory of disappearing microbiota and the epidemics of chronic diseases. Nat. Rev. Immunol. 17, 461–463 (2017).

   Article  CAS  PubMed  Google Scholar 

6. Sprockett, D. D. et al. Microbiota assembly, structure, and dynamics among Tsimane horticulturalists of the Bolivian Amazon. Nat. Commun. 11, 3772 (2020).

   Article  ADS  CAS  PubMed  PubMed Central  Google Scholar 

7. Suzuki, T. A. et al. Codiversification of gut microbiota with humans. Science 377, 1328–1332 (2022).

   Article  ADS  CAS  PubMed  PubMed Central  Google Scholar 

8. Moeller, A. H. et al. Cospeciation of gut microbiota with hominids. Science 353, 380–382 (2016).

   Article  ADS  CAS  PubMed  PubMed Central  Google Scholar 

9. Amato, K. R. & Carmody, R. N. Gut microbial intersections with human ecology and evolution. Annu. Rev. Anthropol. 52, 295–311 (2023).

   Article  Google Scholar 

10. Nyholm, S. V. & McFall-Ngai, M. J. The winnowing: establishing the squid–vibrio symbiosis. Nat. Rev. Microbiol. 2, 632–642 (2004).

    Article  CAS  PubMed  Google Scholar 

11. Moran, N. A., McCutcheon, J. P. & Nakabachi, A. Genomics and evolution of heritable bacterial symbionts. Annu. Rev. Genet. 42, 165–190 (2008).

    Article  CAS  PubMed  Google Scholar 

12. Wibowo, M. C. et al. Reconstruction of ancient microbial genomes from the human gut. Nature 594, 234–239 (2021).

    Article  ADS  CAS  PubMed  PubMed Central  Google Scholar 

13. Moodley, Y. et al. The peopling of the Pacific from a bacterial perspective. Science 323, 527–530 (2009).

    Article  ADS  CAS  PubMed  PubMed Central  Google Scholar 

14. Linz, B. et al. An African origin for the intimate association between humans and Helicobacter pylori. Nature 445, 915–918 (2007).

    Article  ADS  PubMed  PubMed Central  Google Scholar 

15. Falush, D. et al. Traces of human migrations in Helicobacter pylori populations. Science 299, 1582–1585 (2003).

    Article  ADS  CAS  PubMed  Google Scholar 

16. Thorpe, H. A. et al. Repeated out-of-Africa expansions of Helicobacter pylori driven by replacement of deleterious mutations. Nat. Commun. 13, 6842 (2022).

    Article  ADS  CAS  PubMed  PubMed Central  Google Scholar 

17. Mah, J. C., Lohmueller, K. E. & Garud, N. R. Inference of the demographic histories and selective effects of human gut commensal microbiota over the course of human history. Mol. Biol. Evol. 42, msaf010 (2025).

    Article  CAS  PubMed  PubMed Central  Google Scholar 

18. Moeller, A. H., Sanders, J. G., Sprockett, D. D. & Landers, A. Assessing co-diversification in host-associated microbiomes. J. Evol. Biol. 36, 1659–1668 (2023).

    Article  PubMed  PubMed Central  Google Scholar 

19. Nishida, A. H. & Ochman, H. Captivity and the co-diversification of great ape microbiomes. Nat. Commun. 12, 5632 (2021).

    Article  ADS  CAS  PubMed  PubMed Central  Google Scholar 

20. Amato, K. R. et al. Evolutionary trends in host physiology outweigh dietary niche in structuring primate gut microbiomes. ISME J. 13, 576–587 (2019).

    Article  CAS  PubMed  Google Scholar 

21. Ochman, H. et al. Evolutionary relationships of wild hominids recapitulated by gut microbial communities. PLoS Biol. 8, e1000546 (2010).

    Article  PubMed  PubMed Central  Google Scholar 

22. Perez-Lamarque, B. & Morlon, H. Distinguishing cophylogenetic signal from phylogenetic congruence clarifies the interplay between evolutionary history and species interactions. Syst. Biol. 73, 613–622 (2024).

    Article  PubMed  Google Scholar 

23. Good, B. H. Limited codiversification of the gut microbiota with humans. mBio 17, e03727-25 (2026).

    Article  PubMed  PubMed Central  Google Scholar 

24. Washburne, A. D. et al. Methods for phylogenetic analysis of microbiome data. Nat. Microbiol. 3, 652–661 (2018).

    Article  CAS  PubMed  Google Scholar 

25. Sakoparnig, T., Field, C. & van Nimwegen, E. Whole genome phylogenies reflect the distributions of recombination rates for many bacterial species. eLife 10, e65366 (2021).

    Article  CAS  PubMed  PubMed Central  Google Scholar 

26. Lawson, D. J., Hellenthal, G., Myers, S. & Falush, D. Inference of population structure using dense haplotype data. PLoS Genet. 8, e1002453 (2012).

    Article  CAS  PubMed  PubMed Central  Google Scholar 

27. Yahara, K. et al. Chromosome painting in silico in a bacterial species reveals fine population structure. Mol. Biol. Evol. 30, 1454–1464 (2013).

    Article  CAS  PubMed  PubMed Central  Google Scholar 

28. Liu, Z. & Good, B. H. Dynamics of bacterial recombination in the human gut microbiome. PLoS Biol. 22, e3002472 (2024).

    Article  CAS  PubMed  PubMed Central  Google Scholar 

29. Abdill, R. J., Adamowicz, E. M. & Blekhman, R. Public human microbiome data are dominated by highly developed countries. PLoS Biol. 20, e3001536 (2022).

    Article  CAS  PubMed  PubMed Central  Google Scholar 

30. Dixit, P. D., Pang, T. Y., Studier, F. W. & Maslov, S. Recombinant transfer in the basic genome of Escherichia coli. Proc. Natl Acad. Sci. USA 112, 9070–9075 (2015).

    Article  ADS  CAS  PubMed  PubMed Central  Google Scholar 

31. Milkman, R. & Bridges, M. M. Molecular evolution of the Escherichia coli chromosome. III. Clonal frames. Genetics 126, 505–517 (1990).

    Article  CAS  PubMed  PubMed Central  Google Scholar 

32. Qin, J. et al. A metagenome-wide association study of gut microbiota in type 2 diabetes. Nature 490, 55–60 (2012).

    Article  ADS  CAS  PubMed  Google Scholar 

33. Li, J. et al. An integrated catalog of reference genes in the human gut microbiome. Nat. Biotechnol. 32, 834–841 (2014).

    Article  CAS  PubMed  Google Scholar 

34. The Human Microbiome Project Consortium. Structure, function and diversity of the healthy human microbiome. Nature 486, 207–214 (2012).

    Article  ADS  Google Scholar 

35. Beghini, F. et al. Gut microbiome strain-sharing within isolated village social networks. Nature 637, 167–175 (2025).

    Article  ADS  CAS  PubMed  Google Scholar 

36. Valles-Colomer, M. et al. The person-to-person transmission landscape of the gut and oral microbiomes. Nature 614, 125–135 (2023).

    Article  ADS  CAS  PubMed  PubMed Central  Google Scholar 

37. Yu, X. A. et al. Genome-wide sweeps create fundamental ecological units in the human gut microbiome. Nature 655, 202–209 (2026).

    Article  CAS  PubMed  PubMed Central  Google Scholar 

38. Arevalo, P., VanInsberghe, D., Elsherbini, J., Gore, J. & Polz, M. F. A reverse ecology approach based on a biological definition of microbial populations. Cell 178, 820–834 (2019).

    Article  CAS  PubMed  Google Scholar 

39. Haber, M. et al. A rare deep-rooting D0 African Y-chromosomal haplogroup and its implications for the expansion of modern humans out of Africa. Genetics 212, 1421–1428 (2019).

    Article  PubMed  PubMed Central  Google Scholar 

40. Soares, P. et al. The expansion of mtDNA haplogroup L3 within and out of Africa. Mol. Biol. Evol. 29, 915–927 (2012).

    Article  CAS  PubMed  Google Scholar 

41. Potter, B. A. et al. Current evidence allows multiple models for the peopling of the Americas. Sci. Adv. 4, eaat5473 (2018).

    Article  ADS  PubMed  PubMed Central  Google Scholar 

42. Bennett, M. R. et al. Evidence of humans in North America during the Last Glacial Maximum. Science 373, 1528–1531 (2021).

    Article  ADS  CAS  PubMed  Google Scholar 

43. Jouganous, J., Long, W., Ragsdale, A. P. & Gravel, S. Inferring the joint demographic history of multiple populations: beyond the diffusion approximation. Genetics 206, 1549–1567 (2017).

    Article  PubMed  PubMed Central  Google Scholar 

44. Fragiadakis, G. K. et al. Links between environment, diet, and the hunter-gatherer microbiome. Gut Microbes 10, 216–227 (2019).

    Article  CAS  PubMed  Google Scholar 

45. Smits, S. A. et al. Seasonal cycling in the gut microbiome of the Hadza hunter-gatherers of Tanzania. Science 357, 802–806 (2017).

    Article  ADS  CAS  PubMed  PubMed Central  Google Scholar 

46. Ragsdale, A. P. et al. A weakly structured stem for human origins in Africa. Nature 617, 755–763 (2023).

    Article  ADS  CAS  PubMed  PubMed Central  Google Scholar 

47. Li, F. et al. Cardiometabolic benefits of a non-industrialized-type diet are linked to gut microbiome modulation. Cell 188, 1226–1247 (2025).

    Article  CAS  PubMed  Google Scholar 

48. Medina-Muñoz, S. G. et al. Demographic modeling of admixed Latin American populations from whole genomes. Am. J. Hum. Genet. 110, 1804–1816 (2023).

    Article  PubMed  PubMed Central  Google Scholar 

49. Ragsdale, A. P. & Gravel, S. Models of archaic admixture and recent history from two-locus statistics. PLoS Genet. 15, e1008204 (2019).

    Article  CAS  PubMed  PubMed Central  Google Scholar 

50. Pico, T., Mitrovica, J. X. & Mix, A. C. Sea level fingerprinting of the Bering Strait flooding history detects the source of the Younger Dryas climate event. Sci. Adv. 6, eaay2935 (2020).

    Article  ADS  CAS  PubMed  PubMed Central  Google Scholar 

51. Gurven, M. et al. The Tsimane Health and Life History Project: integrating anthropology and biomedicine. Evol. Anthropol. 26, 54–73 (2017).

    Article  PubMed  PubMed Central  Google Scholar 

52. Chen, S. Ultrafast one-pass FASTQ data preprocessing, quality control, and deduplication using fastp. iMeta 2, e107 (2023).

    Article  CAS  PubMed  PubMed Central  Google Scholar 

53. Langmead, B. & Salzberg, S. L. Fast gapped-read alignment with Bowtie 2. Nat. Methods 9, 357–359 (2012).

    Article  CAS  PubMed  PubMed Central  Google Scholar 

54. Babraham Bioinformatics. FastQC: a quality control tool for high throughput sequence data. v.0.12.0 http://www.bioinformatics.babraham.ac.uk/projects/fastqc/(Babraham Institute, 2023).

55. Prjibelski, A., Antipov, D., Meleshko, D., Lapidus, A. & Korobeynikov, A. Using SPAdes DE Novo Assembler. Curr. Protoc. Bioinform. 70, e102 (2020).

    Article  CAS  Google Scholar 

56. Aroney, S. T. N. et al. CoverM: read alignment statistics for metagenomics. Bioinformatics 41, btaf147 https://doi.org/10.1093/bioinformatics/btaf147 (2025).

    Article  CAS  PubMed  PubMed Central  Google Scholar 

57. More, S. & More, A. Assessment the quality of genome assemblies by using QUAST tool for metagenomics. Int. J. Recent Technol. Eng. 8, 4253–4259 (2020).

    Google Scholar 

58. Kang, D. D. et al. MetaBAT 2: an adaptive binning algorithm for robust and efficient genome reconstruction from metagenome assemblies. PeerJ 7, e7359 (2019).

    Article  PubMed  PubMed Central  Google Scholar 

59. Parks, D. H., Imelfort, M., Skennerton, C. T., Hugenholtz, P. & Tyson, G. W. CheckM: assessing the quality of microbial genomes recovered from isolates, single cells, and metagenomes. Genome Res. 25, 1043–1055 (2015).

    Article  CAS  PubMed  PubMed Central  Google Scholar 

60. Bowers, R. M. et al. Minimum information about a single amplified genome (MISAG) and a metagenome-assembled genome (MIMAG) of bacteria and archaea. Nat. Biotechnol. 35, 725–731 (2017).

    Article  CAS  PubMed  PubMed Central  Google Scholar 

61. Chaumeil, P.-A., Mussig, A. J., Hugenholtz, P. & Parks, D. H. GTDB-Tk: a toolkit to classify genomes with the Genome Taxonomy Database. Bioinformatics 36, 1925–1927 (2019).

    Article  PubMed  PubMed Central  Google Scholar 

62. Olm, M. R., Brown, C. T., Brooks, B. & Banfield, J. F. dRep: a tool for fast and accurate genomic comparisons that enables improved genome recovery from metagenomes through de-replication. ISME J. 11, 2864–2868 (2017).

    Article  CAS  PubMed  PubMed Central  Google Scholar 

63. Jain, C., Rodriguez-R, L. M., Phillippy, A. M., Konstantinidis, K. T. & Aluru, S. High throughput ANI analysis of 90K prokaryotic genomes reveals clear species boundaries. Nat. Commun. 9, 5114 (2018).

    Article  ADS  PubMed  PubMed Central  Google Scholar 

64. Varghese, N. J. et al. Microbial species delineation using whole genome sequences. Nucleic Acids Res. 43, 6761–6771 (2015).

    Article  CAS  PubMed  PubMed Central  Google Scholar 

65. Olm, M. R. et al. Consistent metagenome-derived metrics verify and delineate bacterial species boundaries. mSystems 5, e00731-19 (2020).

    Article  PubMed  PubMed Central  Google Scholar 

66. Hyatt, D. et al. Prodigal: prokaryotic gene recognition and translation initiation site identification. BMC Bioinform. 11, 119 (2010).

    Article  Google Scholar 

67. Novichkov, V., Kaznadzey, A., Alexandrova, N. & Kaznadzey, D. NSimScan: DNA comparison tool with increased speed, sensitivity and accuracy. Bioinformatics 32, 2380–2381 (2016).

    Article  CAS  PubMed  Google Scholar 

68. Almeida, A. et al. A new genomic blueprint of the human gut microbiota. Nature 568, 499–504 (2019).

    Article  ADS  CAS  PubMed  PubMed Central  Google Scholar 

69. Steinegger, M. & Söding, J. MMseqs2 enables sensitive protein sequence searching for the analysis of massive data sets. Nat. Biotechnol. 35, 1026–1028 (2017).

    Article  CAS  PubMed  Google Scholar 

70. Carter, M. M. Prehistoric global migration of vanishing gut microbes with humans [Data set]. Zenodo https://doi.org/10.5281/zenodo.21765848 (2026).

71. Carter, M. M. Prehistoric global migration of vanishing gut microbes with humans [Computer software]. Zenodo https://doi.org/10.5281/zenodo.21206858 (2026).

Download references

Author notes

1. These authors contributed equally: Matthew M. Carter, Zhiru Liu, Matthew R. Olm

Authors and Affiliations

1. Department of Microbiology and Immunology, Stanford University School of Medicine, Stanford, CA, USA

   Matthew M. Carter, Matthew R. Olm, David A. Relman, Erica D. Sonnenburg & Justin L. Sonnenburg

2. Department of Applied Physics, Stanford University, Stanford, CA, USA

   Zhiru Liu & Benjamin H. Good

3. Department of Integrative Physiology, University of Colorado Boulder, Boulder, CO, USA

   Matthew R. Olm & Parsa Ghadermazi

4. Department of Anthropology, University of California Santa Barbara, Santa Barbara, CA, USA

   Melanie Martin & Michael Gurven

5. Department of Anthropology, University of Washington, Seattle, WA, USA

   Melanie Martin

6. Department of Microbiology and Immunology, Wake Forest University School of Medicine, Winston-Salem, NC, USA

   Daniel D. Sprockett

7. School of Human Evolution and Social Change, Arizona State University, Tempe, AZ, USA

   Benjamin C. Trumble

8. Center for Evolution and Medicine, Arizona State University, Tempe, AZ, USA

   Benjamin C. Trumble

9. Institute of Human Origins, Arizona State University, Tempe, AZ, USA

   Benjamin C. Trumble

10. Economic Science Institute, Chapman University, Orange, CA, USA

    Hillard Kaplan

11. Toulouse School of Economics, Toulouse, France

    Jonathan Stieglitz

12. Universidad Mayor de San Simón, Cochabamba, Bolivia

    Daniel Eid Rodriguez

13. Department of Medicine, Stanford University School of Medicine, Stanford, CA, USA

    David A. Relman

14. Infectious Diseases Section, Veterans Affairs Palo Alto Health Care System, Palo Alto, CA, USA

    David A. Relman

15. Department of Biology, Stanford University, Stanford, CA, USA

    Benjamin H. Good

16. Chan Zuckerberg Biohub, San Francisco, CA, USA

    Benjamin H. Good & Justin L. Sonnenburg

17. Center for Human Microbiome Studies, Stanford University School of Medicine, Stanford, CA, USA

    Justin L. Sonnenburg

Contributions

Conceptualization: M.M.C., Z.L., M.R.O., B.H.G. and J.L.S. Methodology: M.M.C., Z.L., M.R.O., P.G., B.H.G. and J.L.S. Software and analysis: M.M.C., Z.L., M.R.O., P.G. and B.H.G. Writing, review and editing: M.M.C., Z.L., M.R.O., M.M., D.D.S., B.C.T., H.K., J.S., D.E.R., D.A.R., E.D.S., M.G., B.H.G. and J.L.S. Visualization: M.M.C., Z.L., M.R.O., B.H.G. and J.L.S. Administration of the THLHP: B.C.T., H.K. and M.G. Project administration: B.H.G. and J.L.S. Supervision: B.H.G. and J.L.S. Funding acquisition: D.A.R., B.H.G. and J.L.S. All authors read and approved the manuscript.

Corresponding authors

Correspondence to Benjamin H. Good or Justin L. Sonnenburg.

Competing interests

The authors declare no competing interests.

Peer review information

Nature thanks the anonymous reviewers for their contribution to the peer review of this work. Peer reviewer reports are available.

Open Access This article is licensed under a Creative Commons Attribution-NonCommercial-NoDerivatives 4.0 International License, which permits any non-commercial use, sharing, 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 licence, and indicate if you modified the licensed material. You do not have permission under this licence to share adapted material derived from this article or parts of it. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence 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 licence, visit http://creativecommons.org/licenses/by-nc-nd/4.0/.

Reprints and permissions