Description of populations

Our participants included 54 individuals from four Himalayan groups, including Chepang (N = 14), Raji (N = 9), Raute (N = 11), and Tharu (N = 20), with median age of 40 years (SD ± 14 years) from rural villages in Nepal (Fig 1 and S1 Table). These four populations are long-term residents of the Himalayan foothills (altitude less than 1,000 m). Ethnographic, linguistic, and cultural data suggest that these populations are of East Asian ancestries, they speak closely related languages, and their cultural practices are similar to one another [32–34]. Although all four of the Himalayan populations in this study were foragers until recently [35–38], habitat loss due to rapid deforestation, population expansions of non-native groups, establishment of new settlements, and construction of modern highways led to settlement of these groups at various time points in the last 300 years. Historical records indicate that the Tharu gradually transitioned into agrarian lifestyles beginning in the late 18th century (250–300 years ago) [38]. They have fully transitioned into farming and are virtually completely disengaged from foraging practices. The population size of the Tharu is approximately 1.5 million, and they are distributed throughout the Terai plains in Nepal [39]. Historically, the Raute, Raji, and Chepang were seminomadic foragers, and their diets included native tubers, greens, and fruits from the jungle and wild honey, fish, and occasional game [36,37,40]. The Raute and Raji abandoned their foraging lifestyles in the 1980s [35,36]. While the Raute have settled in the remote hills in western Nepal, the Raji have settled in the Terai plains, which are relatively more urbanized. The current census size of Raute and Raji are approximately 650 and approximately 3,750, respectively [39]. The Chepang were fully nomadic at least until 1848 [41] and began supplementing their foraging practices with subsistence agriculture less than a century ago [37]. The Chepang population size is approximately 48,500 [39]; however, they exist as fragmented tribes in small, geographically isolated villages of a few hundred individuals deep within the hills of lower Himalaya. The Chepang in this study currently inhabit a remote village that is devoid of modernity, with no electricity, running water, irrigation, fertilizers, modern machines, or marketplaces. They still practice slash and burn agriculture and are completely dependent on rainwater for farming. Because yields from such traditional farming are low, their daily diet consists of wild plants such as sisnu (nettles) that are foraged from the forests.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 1. Sampling locations and habitats of the Himalayan populations in Nepal.

(A) Map displaying the geographical locations of sampled villages in southern Nepal (altitudes <1,000 m above sea level, latitude 26.97–29.15 °N). The Tharu are geographically most distant from the Raute and Raji and reside closer to the Chepang. (B) Habitats of each population. From top left in a clockwise direction, the remote Chepang village, the Raute village, the Tharu harvesting rice, and the Raji village. The census population sizes of the Raute, Raj, Chepang, and Tharu are 650, 3,758, 48,476, and 1.5 million, respectively.

https://doi.org/10.1371/journal.pbio.2005396.g001

Lifestyle gradients in the Himalayan populations

We conducted surveys to assess the extent of lifestyle change as these seminomadic populations transitioned to farming in the last few hundred years. The survey questionnaire included questions pertaining to current dietary practices, traditional and modern medicines, and several environmental factors, including sources of drinking water, types of cooking fuel, alcohol use, and tobacco consumption (N = 53, S2 Table). We also surveyed presence of parasites in our participants microscopically.

Supervised learning using a Random Forest classifier model on the survey data (including intestinal parasite load) assigned the individuals to their respective populations with high overall accuracy (94%, AUC = 0.997, S1 Fig). The Chepang, Raute, and Tharu were classified with 100% accuracies, indicating these populations have distinct lifestyles (S3 Table). 67% of the Raji individuals were classified accurately as Raji while the remainders were classified as Tharu. A correspondence analysis (CA) of the survey data (including intestinal parasite load) also revealed lifestyle differences between these populations (Fig 2A). The first CA dimension (CA1) explained 15.8% variation in the data and was strongly correlated with lifestyle gradients. Along CA1, samples progressed from the Chepang foragers at one extreme, to the Raute and Raji transitioning populations, and then to the Tharu farmers at the opposite extreme (Fig 2B). Despite the geographical distance between them, the Raji lifestyle appears to be more similar to that of the Tharu farmers, consistent with the Raji settlement occurring in a more urbanized setting compared to the Raute. Similarly, the Raute reside in geographical proximity to the Raji, although their lifestyle partitions between the Raji and the Chepang, indicating geographical proximity is not driving the lifestyle differences.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 2. CA based on survey questionnaires and parasite assessment in the Himalayan populations.

(A) First two dimensions of the CA and the amount of variation explained are shown. Each circle represents an individual, and colors represent the populations. (B) Distribution of populations along the primary CA1 axis shows patterns of separation by lifestyles. Chepang foragers (red) and Tharu farmers (blue) are on two extreme ends of CA1. In between the two are the Raute (yellow) and Raji (cyan), the two communities that are transitioning from foraging to farming. (C) Factors in gold are those that have more than expected eigenvalues and thus contribute most to the top two dimensions in the CA. The data underlying this figure can be found in S4 Table. CA, correspondence analysis.

https://doi.org/10.1371/journal.pbio.2005396.g002

A total of 10 variables contributed highly to the first two CA dimensions, and most of them are strongly associated with dietary differences and modernity (Fig 2C). These differences are described in detail in S2 Fig. Briefly, foraged plants such as sisnu (nettles) and jaand, a slushy alcoholic beverage made from fermenting millet or corn, are staples of the Chepang diet. In contrast, sisnu and jaand consumption was minimal among the Raute, Raji, and Tharu. Also, perceived food scarcity was higher in the Chepang and Raute relative to the Raji and Tharu. Although meat consumption was low across all four populations, the Tharu consumed animal products such as yogurt more frequently than the other three populations. Furthermore, the Tharu and Raji also showed increased signs of modernity. For example, they have installed tube wells at their homes, enabling access to underground water for drinking. In contrast, the Chepang and Raute still fetch drinking water from rivers and streams. Also, use of solid biomass fuel was lower in the Tharu and Raji, while the Chepang and Raute are still completely dependent on burning firewood for cooking. Although we detected low overall levels of intestinal parasites across the participants, Ascaris, Entamoeba, Trichuris, Hymenolepis, and Coccidia were detected in some, and most of the infected were the Chepang. Together, the diet and lifestyle assessments provide unbiased support that the four populations represent a gradient from traditional to increasingly agrarian and urban lifestyles.

Gut microbiome composition varies by lifestyles

In order to assess whether the gut microbiome varies across lifestyles, we characterized the gut bacterial composition of these populations using the Illumina MiSeq to sequence the V4 region of 16S ribosomal RNA (rRNA) gene obtained from a total of 79 stool samples (including technical replicates), with an average of 11,570 (±4,653) high-quality reads per sample (S3 Fig and S4 Table). Since flash freezing of the samples was not possible in the remote sampling areas in the Himalaya, we used commercially available DNAGenotek OMNIgene kits to collect stool samples from the four populations (N = 54). We also collected stool samples from 10 Americans of European descent using OMNIgene kits and compared them with freshly frozen samples to evaluate whether the preservation method affected the microbiome profile. The 16S rRNA profiles of the same samples stored by flash freezing or by OMNIgene were remarkably similar, with reproducible differences in minor taxa (Euryarcheota and Cyanobacteria), demonstrating the reliable preservation of microbiome composition with the OMNIgene kits (S4 Fig). Due to the reproducible, albeit minor, differences between the two collection methods, we used the OMNIgene data from the Americans for consistency in subsequent comparative analyses. The American samples provide a thoroughly investigated population as an industrialized reference for the Himalayan data.

Comparison of the community structure in the five study populations using unweighted UniFrac distances, a measure of compositional similarity that includes the phylogenetic relatedness between microbiomes, showed that the gut microbial composition varies across populations (P < 2.2 × 10⁻¹⁶, Kruskal–Wallis test, S5 Table). Within Himalaya, the Chepang foragers were closest to the Raute, and the distance between the two was significantly smaller than the Chepang–Tharu distance and marginally smaller than the Chepang–Raji distance (false discovery rate [FDR] adjusted P = 5.7 × 10⁻⁴ and 0.057, respectively; Dunn’s posthoc test). The Raute, Raji, and Tharu were equidistant from one another (FDR adjusted P = 0.99 for all pairwise comparisons, Dunn’s posthoc test). Similar results were also observed with weighted UniFrac and Bray–Curtis distances, both of which take the taxa abundance into account (S5 Table). These results suggest that differences in traditional lifestyles as these groups transition from foraging to farming influence their gut microbiomes. Comparison of these traditional Himalayan populations with industrialized Americans showed that all four Himalayan populations exhibited much larger distances from the Americans than when compared to one another (P < 1.3 × 10⁻⁵ for all pairwise comparisons, Dunn’s posthoc test, S5 Table). The Chepang were the most distant from the Americans, followed by the Raute, while the Raji and Tharu were equally close to the Americans.

Visualization of these distances using a Principal Coordinates Analysis (PCoA) revealed separation of populations along the top two dimensions (P = 1 × 10⁻⁵, permutational multivariate analysis of variance [PERMANOVA], Fig 3A). Furthermore, gradients in lifestyles were reflected by the distribution of populations along the primary axis (PCoA1, Fig 3B). These distributions remained consistent when using Bray–Curtis and weighted UniFrac distances as well (P = 1 × 10⁻⁵ for both, PERMANOVA, S5 and S6 Figs).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 3. Gut microbiome compositions show gradients corresponding to lifestyles.

(A) PCoA of the unweighted UniFrac distances of the 16S rRNA data colored by populations. Each dot represents an individual, and colors indicate the populations. Chepang foragers (red), Raute (yellow), and Raji (cyan) communities that are transitioning from foraging to farming; Tharu farmers (blue); and Americans (orange). (B) Distributions of populations along the PCoA1 axis show patterns of separation by lifestyles. (C) Gut microbial composition of the Himalayan populations represented by the primary dimension of the unweighted UniFrac distance (PCoA1) strongly correlates with lifestyle differences represented by the top dimension of the corresponding analysis performed on the survey data (CA1, Spearman’s Rho = 0.44 and P value = 0.001). Correlation between CA2 and PCoA1 was not statistically significant. The data underlying this figure can be found in S1 Data. CA, correspondence analysis; PCoA, Principal Coordinates Analysis; rRNA, ribosomal RNA.

https://doi.org/10.1371/journal.pbio.2005396.g003

When American microbiomes were eliminated from the PCoAs, the gradient between the Himalayan populations remained pronounced (P = 1 × 10⁻⁵, PERMANOVA, S7 Fig). Among the four Himalayan populations, the strongest separation was observed between the Chepang foragers and the Tharu farmers.

A random forest classifier based on the 16S rRNA-defined amplicon sequence variant (16S ASV) data assigned the Chepang, Tharu, and American individuals to their respective source populations with 79%, 100%, and 100% accuracies (overall accuracy = 66%, AUC = 0.9, S8 Fig and S3 Table). The classification accuracy for the Raute and Raji, the two populations that recently transitioned from foraging to farming, were relatively poor (<10%). While some of the individuals from these groups were classified as the Chepang, others were classified as the Tharu. However, none of the Himalayan individuals were classified as American. These results show that the gut microbiome of the Chepang foragers differs from that of the Tharu farmers, while that of the Raute and Raji reflect a transitional state in their lifestyles. They also indicate that the gut microbiome compositions of the Himalayan populations are distinct from those of the Americans. Therefore, these findings collectively indicate that the transition from foraging to farming is accompanied by noticeable shifts in gut microbiome, which may be further exacerbated in industrial populations.

To formally evaluate whether variation in gut microbiota reflects lifestyle differences within Himalaya, we assessed associations between the respective primary dimensions from the lifestyle questionnaire, parasite analysis (CA1), and gut microbial composition analysis (PCoA1 calculated using the four Himalayan populations) (Fig 3C and S5 Fig). We found that the CA1 was strongly correlated with the PCoA1 obtained from all of the three distance matrices (Spearman’s Rho = 0.47, 0.44, and 0.28 for Bray–Curtis; unweighted UniFrac; and weighted UniFrac distances, respectively; P values = 4.5 × 10⁻⁴, 1.1 × 10⁻³, and 0.05; correlation test). The CA1 was also correlated with PCoA2 of all three distance matrices (Spearman’s Rho = 0.26, 0.44, and 0.39; P values = 0.06, 0.001, and 0.004 for Bray–Curtis; unweighted UniFrac; and weighted UniFrac distances; correlation test). Conversely, no significant correlations were detected between CA2 and either of the PCoA axes from all three distances (P value > 0.05, correlation test). Notably, CA1 but not CA2 is associated with lifestyle gradient in Himalaya (Fig 2). Strong and consistent correlations between CA1 and PCoA axes indicate that gut microbiome compositions of the Himalayan populations mirror their lifestyles.

Gut bacterial diversity (alpha diversity) does not vary across lifestyles

Previous studies have suggested that elevated species diversity in the gut microbiome is a hallmark of traditional populations [19,24]. We assessed the alpha diversity in the five study populations using species richness and Shannon’s H at various rarefaction depths ranging from 10–6,500 reads (Fig 4). Species richness measures the presence and absence of taxa, whereas Shannon’s H additionally accounts for the relative abundances of each taxon within each population. We compared alpha diversity across the five populations at a rarefaction depth of 3,000 to include all 64 samples and at a higher rarefaction depth of 6,500, which included 61 samples. Regardless of the rarefaction depth, species richness was not significantly different between any of the five populations (P > 0.05, Kruskal–Wallis test). We did find marginally significant differences in Shannon’s H between these populations (P = 0.01 and 0.03 at 3,000 and 6,500 rarefaction depths, Kruskal–Wallis test). A posthoc pairwise comparison of all five populations showed that only the alpha diversity in the Tharu was slightly lower than that in the Americans (FDR adjusted P = 0.02 and 0.045 at the two rarefaction depths, respectively; Dunn’s posthoc test). Next, we evaluated association between the 10 factors that differentiate the lifestyle of the Himalayan populations (Fig 2) with the two alpha diversity measures. Neither species richness nor Shannon’s H were significantly associated with any of these factors at either rarefaction depth (P > 0.05, nonlinear mixed effects model). Finally, we assessed additional metrics of diversity (Fisher’s alpha, Simpson’s D), which similarly fail to differentiate populations (S9 Fig). These results indicate that lifestyle differences among the Himalayan populations or between these populations and Americans have little effect on the alpha diversity of their gut microbiome.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 4. Alpha diversity across lifestyles.

Rarefaction curves showing two commonly used measures of alpha diversity—species richness (top) and Shannon’s H (bottom) calculated by subsampling 10–6,500 reads per sample. No significant differences in species richness was detected between the five study populations at a lower depth of 3,000 reads per sample, which included all 64 samples, or at 6,500 reads per sample, which included 61 samples. Shannon’s H was significantly lower in the Tharu relative to the Americans at both rarefaction depths. No differences in any of these two alpha diversity metrics were observed between the Chepang, Raji, Raute, and the Americans. Population labels are colored to indicate the range of different lifestyles (red, foragers; yellow and cyan, former foragers; blue, farmers; orange, industrialists). The data underlying this figure can be found in S1 Data.

https://doi.org/10.1371/journal.pbio.2005396.g004

Bacterial taxa are associated with lifestyle transitions

Although lifestyle differences have little effect on the alpha diversity, gut microbiome compositions of the Himalayan populations reflected the gradient in their lifestyles. To identify taxa driving the differences in the gut microbiomes across lifestyles, we compared the abundance of individual phyla across the five populations using a negative binomial generalized linear model (GLM), as implemented in differential expression analysis for sequence count data version 2 (DESeq2) [42]. Differential abundances were detected for six out of 10 phyla (FDR adjusted P values are shown in S6 Table), and four of the six phyla reflect a traditional western lifestyle gradient. The Himalayan populations were characterized by higher abundance of Proteobacteria, while abundances of Actinobacteria, Firmicutes, and Verrucomicrobia were highest in the Americans, intermediate in the farmers (Tharu, Raji, and Raute), and lowest in the Chepang foragers (Fig 5A). Higher levels of Proteobacteria and lower levels of Actinobacteria and Verrucomicrobia are common features of many traditional human gut microbiomes around the world [19,21,24,29].

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 5. Distinctions in the gut microbiome across lifestyles.

(A) Phyla with most significant differences in abundances between the five populations. Abundances of Firmicutes, Verrucomicrobia, and Actinobacteria reflect gradients of traditional industrialized lifestyles. Proteobacteria distinguishes rural Himalayan populations from the Americans. (B) Heatmap displaying 52 genera with significantly different abundance across the five populations. Bars on the top represent the grouping of individuals in the heatmap columns by their populations or lifestyles. Genera labels in rows are colored by their phylum. Purple, Actinobacteria; dark blue, Bacteroidetes; light red, Elusimicrobia; orange, Firmicutes; light blue, Proteobacteria; magenta, Spirochaetes; light pink, Tenericutes; brown, Verrucomicrobia. Heatmap colors reflect relative abundances of each genus. Among the Himalayan populations, Ruminobacter, Campylobacter, unknown Veillonellaceae genus, Bulleidia, Weissella, Treponema, Barnesiella, Odoribacter, Alistipes, and Bifidobacterium differed significantly. The data underlying this figure can be found in S6 and S7 Tables.

https://doi.org/10.1371/journal.pbio.2005396.g005

To characterize the taxonomic differences between populations at a finer taxonomic level, we repeated the above analysis at the genus level and identified 52 out of 116 genera that showed significant differences in abundance across the five populations (Fig 5B, FDR adjusted P values are shown in S7 Table). The majority of these genera show consistent differences along the lifestyle gradient within the Himalayan samples (Fig 5B). For example, among the Himalayan populations, the Chepang foragers were enriched for Ruminobacter, Campylobacter, and Treponema relative to the Tharu farmers (S10 Fig). Although we did not detect significant differences in the abundance of Bacteroidetes phylum across these populations, several members of this phylum distinguished the Himalayan and American populations. The rural Himalayan communities were enriched for Prevotella, Alloprevotella, and Anaerophaga and significantly depleted in Bacteroides, Alistipes, Butyricimonas, Odoribacter, and Barnesiella. 29 genera belonging to Firmicutes differed significantly across the five populations, and their distribution was complex across these populations (S11 Fig). The Himalayan populations were enriched for Clostridium sensu stricto, Catenibacterium, Lactobacillus, Bulleidia, Sarcina, Enterococcus, Eubacterium, Oribacterium, Mogibacterium, Mitsuokella, Allisonella, Weissella, Papilbacter, and two unknown genera of Erysipelotrichaceae and Veillonellaceae families. Alternatively, abundances of several Clostridium genera, Oscillibacter, Blautia, Butyriciococcus, Anaerostipes, and Flavonifractor were elevated in the Americans. The Americans also showed highest abundances of Bifidobacterium (Actinobacteria) and Akkermansia (Verrucomicrobia), both of which were extremely low in the Chepang foragers and intermediate in the Tharu farmers. Elevated abundances of Treponema and Prevotella with reduction of Bacteroides and Bifidobacterium is a characteristic feature of gut microbiomes of foraging communities [19,21,24,29].

Microbiome structure across lifestyles

In addition to the individual taxa that differ across subsistence strategies, we wanted to determine whether microbial networks are also associated with lifestyle differences [24,43,44]. To understand how the gut microbiome network structure varies across these populations, we calculated the correlations between all pairs of bacterial genera in the gut using Sparse Correlations for Compositional data (SparCC) [45]. Clustering based on these correlations revealed seven bacterial coabundance groups (CAGs, S12 Fig). The dominant genera that defined these CAGs are Ruminococcus, Bacteroides, Roseburia, Escheria/Shigella, Suturella, Prevotella, and Dialister (Fig 6A and S13 Fig). These seven CAGs showed two antagonistic clusters: one cluster contains CAGs defined by Ruminococcus, Bacteroides, and Roseburia, and the second cluster contains CAGs defined by Prevotella, Escheria/Shigella, and Dialister (S13 Fig). Notably, the CAG dominated by Prevotella is most prominently represented in the Chepang and Raute, while members of the CAG dominated by Bacteroides are elevated in the Raji and Tharu (Fig 6B). Within the Prevotella CAG, Treponema and Ruminobacter are characteristic of the Chepang foragers. Conversely, the American gut is highly depleted of the Prevotella CAG and is dominated by the Bacteroides CAG. The results suggest how these changes in the microbiome that accompany lifestyle transitions may be viewed both at the level of individual taxa as well as higher-order community structure. Transitions from foraging to farming in Himalayan populations show changes in gut microbial networks, which appear to become more profound in industrialized societies.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 6. Microbial coabundance networks across lifestyles.

(A) Visualization of the occurrence patterns of bacterial genera using the Fruchterman–Reingold force–directed layout algorithm. Nodes (circles) represent bacterial genera, node colors represent the seven CAGs, and node sizes represent genus abundance. Only the most dominant genera in each CAG are labeled. Edges represent the significant and positive correlations between genera. Members of the red, blue, and yellow CAGs are tightly correlated to one another and mostly negatively correlated with the members of cyan, magenta, gold, and green CAGs. Labels with “x__unk” indicate taxa with unknown classification level. (B) The relative proportions of these CAGs vary across the lifestyle gradient. Chepang foragers show elevated proportion of the magenta CAG, which is dominated by Prevotella, Succinivibrio, Ruminobacter, and Treponema. This CAG decreases in the Raute, Raji, and Tharu farmers with concurrent increase in the blue CAG, which is dominated by Bacteroides, Faecalibacterium, and Bifidobacterium. The American gut is dominated by the blue CAG and highly depleted of the magenta CAG. The data underlying this figure can be found in S1 Data. c, class; CAG, coabundance group; f, family; g, genus; o, order.

https://doi.org/10.1371/journal.pbio.2005396.g006

Factors associated with gut microbiome composition in the Himalaya

We next assessed whether any of the 10 dietary and environmental factors that differentiate the Himalayan populations (from Fig 2) correspond to the variation in gut microbiome composition. A canonical correspondence analysis (CCA) revealed that the 10 factors collectively explain 28% of the gut microbiome variation within Himalaya, while 72% of the variation remained unexplained. Of the 10 variables, the source of drinking water and use of solid biomass fuel were significantly associated with the gut microbiome composition in the Himalayan populations (P value = 0.009 and 0.028, respectively; permutation test). Both of these factors contributed most to the first CCA axis (CCA1), which distinguished the Chepang and Raute individuals who drink river water and exclusively burn solid biomass fuel for cooking from the Raji and Tharu who drink underground water and use biogas for cooking (Fig 7). As an alternative approach, we assessed associations between the gut microbiome composition of the Himalayan individuals (using Bray–Curtis, unweighted UniFrac, and weighted UniFrac) and the 10 lifestyle-associated factors by performing a PERMANOVA (S14 Fig). These analyses also revealed that drinking water was strongly associated with the gut microbiome variation within Himalaya (P = 0.001 for all three distances; effect sizes = 0.096, 0.1095, and 0.088 for Bray–Curtis; unweighted UniFrac; and weighted UniFrac distances, respectively). Individuals who drank river water had higher abundances of Treponema, and those who drank underground water had elevated levels of Fusobacterium (FDR adjusted P value = 0.01 and 0.003, respectively; Mann–Whitney test). Although cooking fuel was significantly associated with overall composition, none of the individual genera reached statistical significance after correcting for multiple testing.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 7. Gut microbiome composition is associated with environmental factors in Himalaya.

(A) The two primary CCA axes and the proportion of constrained variance they explain are shown. Triangles represent individuals and circles represent genera. Individuals and genera are color coded by their respective drinking water sources and phyla. Drinking water and cooking fuel contributed most to CCA1, and sisnu (nettles) contributed most to CCA2. Genera labeled in grey contribute to the top two CCA axes. Among these, Fusobacterium and Treponema were significantly associated with drinking water. (B) PCoA of the unweighted UniFrac distances. Each dot represents an individual, colors indicate the two drinking water sources, and shapes represent different populations. Gut microbiomes of the Chepang (circles) and Raute (diamonds) who drink water from rivers and streams vary significantly from those of the Raji (squares) and Tharu (triangles). Statistical significance was assessed using PERMANOVA using the 10 variables that differentiate their lifestyles (P value = 0.0001). The data underlying this figure can be found in S1 Data and S4 Table. CCA, canonical correspondence analysis; PCoA, Principal Coordinates Analysis; PERMANOVA, permutational multivariate analysis of variance.

https://doi.org/10.1371/journal.pbio.2005396.g007

To assess whether the association between gut microbiome and drinking water extend beyond the Nepali populations, we reanalyzed an independent 16S rRNA amplicon data set from Hadza hunter–gatherers from the Hukamako camp (N = 60) [21]. In the late dry season, the Hadza use water from two distinct sources—springs (N = 22) and streams (N = 38). We used a CCA to assess the associations between the gut microbiome of the Hadza and several dietary and environmental factors, including 72-hour recall of baobab, berries, honey, maize, meat, and tuber consumption; alcohol and cigarette use; as well as differences in drinking water sources (S8 Table). These variables collectively explained 16.5% of the gut microbiome variation in the Hadza gut microbiome. Among the variables used in the CCA, difference in drinking water source was most strongly associated with the Hadza gut microbiome composition followed by honey consumption (P = 0.0001 and 0.03, respectively; permutation test; S15 Fig). We also performed a PERMANOVA to assess associations between the gut microbiome composition of the Hadza individuals and their dietary and environmental factors using Bray–Curtis, unweighted UniFrac, and weighted UniFrac distances. All three analyses revealed the association between drinking water source and the Hadza gut microbiome (P = 0.001, 0.002, and 0.003, PERMANOVA; effect sizes = 0.051, 0.051, and 0.059 for Bray–Curtis, unweighted, and weighted UniFrac distances, respectively; S15 Fig). Therefore, our results in the Hadza and the Himalayan populations suggest that drinking water is strongly associated with the human gut microbiome and emphasize the need for additional work to elucidate the mechanisms by which drinking water may influence the gut microbiome.

Ethics statement

This work was approved by the Ethical Review Board of the Nepal Health Research Council (NHRC) as well as by the Stanford University Institutional Review Board.

Study sites, participating individuals, and sample collection

Stool samples were collected with informed consent from 56 genetically unrelated adult participants (over 18 years old with different grandparents) from four indigenous Himalayan populations from Nepal and 10 adult Americans of European descent. Indigenous populations from Nepal included Chepang (N = 14), Raji (N = 10), Raute (N = 12), and Tharu (N = 20) inhabiting Chitwan, Bardia, Dadeldhura, and Sarlahi districts, respectively. The samples were collected in winter of 2016 (March and April) with consent from all participants.

In addition to collecting the fecal samples, we also obtained ethnolinguistic, demographic, environmental, and dietary data from the Himalayan participants using a survey questionnaire specifically designed for this study. The survey questionnaire assessed participants’ age, gender, diet, health status, use of medication, and behavioral practices such as tobacco and alcohol consumption, along with several environmental variables (S2 Table). In addition, we also visually inspected the stool samples of each individual under the microscope for the presence of intestinal parasites (triplicate slides per individual). Participants’ responses to survey data questionnaires are included in S4 Table.

DNA extractions

Freshly produced stool samples from the Himalayan participants were collected on a clean OMNIgene gut accessory collection paper (OM-AC1). About 500 mg of the stool samples was transferred to the OMNIgene gut kit collection tube containing the stabilizing buffer using the clean spatula provided with the kit. The tubes were shaken hard in a back and forth motion until the fecal samples were completely homogenized. Tubes were transported at room temperature within 48–72 hours of collection to the Tribhuvan University Institute of Medicine, Kathmandu, Nepal, where they were transferred to −80 °C until DNA extraction. DNA was extracted using a MolBio Power Soil Kit according to the manufacturer’s protocol. Extracted DNA was shipped to Stanford University on dry ice and stored at −20 °C until sequencing. Samples from Americans were collected from volunteers at Stanford University in 15-ml centrifuge tubes and transported to the laboratory on ice. Half of each sample was immediately frozen at −80 °C. From the other half, 500 mg stool was transferred to OMNIgene collection tubes and kept at room temperature for 48–72 hours after which they were stored at −80 °C. DNA was extracted from both sets of samples simultaneously using the MolBio Power Soil Kit according to the manufacturer’s protocol and stored at −20 °C until sequencing.

16S rRNA gene amplicon sequencing and analyses

The V4 region of the 16S rRNA gene was PCR amplified using the primers and protocols described previously [55]. The amplified DNA fragments were multiplexed and subjected to paired-end sequencing using Illumina MiSeq. Of the 66 samples, one yielded very low levels of DNA and another failed the paired-end sequencing. After discarding these two samples, the final data set included 64 individuals (14 Chepang, 9 Raji, 11 Raute, 20 Tharu, and 10 Americans). The amplification primers and barcodes used for multiplexing are described in S4 Table.

Paired-end reads were processed using DADA2 [56] and subsequently analyzed in R using phyloseq [57]. In order to identify high quality sequences, reads were trimmed to 150 bp. Sequences with N nucleotides and/or >2 expected errors were discarded (maxN = 0, maxEE = 2, truncQ = 2), and sequence variants were inferred by pooling reads from all samples (pool = TRUE). Sequence tables were then created by merging paired-end reads. A naïve Bayesian classifier method [58] implemented in DADA2 algorithm was used to assign taxonomy using the RDP v14 training set [59]. Multiple alignment was conducted using DECIPHER [60] package in R, and a maximum likelihood phylogenetic tree was constructed using phangorn [61], with a neighbor-joining tree as the starting point.

A total of 1,183,760 merged reads passed quality control, and 1,630 taxa were initially identified. After removing chimeric sequences, which constituted 22% of the reads, 921,345 merged reads remained. Further elimination of low-abundance phyla—Synergistetes and Deferribacteres—that were observed only once across all samples resulted in 883 taxa in the data set. After quality control, mean (±SD) sequencing depth per sample was 11,570 (±4,653). We performed three technical replicates of the frozen sample for one individual and a total of five replicates for two additional individuals for the OMNI samples. Since we did not observe marked differences in the technical replicates (S3 Fig), we retained the sample with highest coverage for these individuals. After removing the replicate samples, 64 individuals and 875 taxa remained in the final data set.

Hadza sample collection and 16S sequencing

Stool samples and dietary recall from 60 Hadza individuals were conducted in the field at the time of sampling with the aid of an interpreter. Following informed consent, each participant provided a list of the plants, animals, and animals products consumed over the previous 72 hours, including alcohol and cigarette use. Location and type of water source was also recorded. Although the Hadza consume water primarily from a single source near their camp, foraging activities often take subjects several kilometers away from camp where their water source may vary. Raw 16S reads from the Hadza were previously published [21] and were processed using DADA2, as described above. This data set included 1,038,333 nonchimeric reads from 60 individuals that were assigned to 1,511 taxa.

Random forest classifier model

A random forest classifier with 5,000 trees was constructed using all 35 variables (S3 Table) from the survey data. The R-package randomForest [62] was used to build the trees, and its “tuneRF” function was used to assess the optimal number of variables randomly sampled as candidates at each split. (“mtry” parameter, mtry = 6 for survey data). We also repeated these analyses on the 16S data using the RSVs as features and using mtry = 29 as determined by the tuneRF. The R-package “pROC” [63] was used to calculate and plot the area under the receiver operating characteristic curves (AUC) for each of the populations. In addition, brier skill scores (BSSs) in R-package “verification” [64] was used to assess the calibration of the RF model.

Statistical analyses

Intestinal pathogens and all 35 variables recorded from the Himalayan populations were used to perform CA of the survey data using “FactoMineR” package in R [65]. Associations between rows and columns in the correspondence analysis were evaluated by performing a chi-squared test (P = 3.9 × 10⁻⁶). The contributions of each factors to the top two dimensions of CA were visualized using the “fviz_contrib” function in R-package factoextra [66]. The expected contribution to the top two dimensions under a uniform model was determined, and factors that contributed more than the expected were considered important in differentiating lifestyles.

CCA was performed in the Himalayan populations using the 10 variables that differentiated lifestyles in the correspondence analysis by calling the “cca” function from vegan package [67] via phyloseq. For the Hadza, CCA was performed using nine variables including six dietary variables (baobab, berries, maize, tubers, honey, and meat consumption in the last 72 hours), alcohol and cigarette consumption, as well as source of drinking water. RSV counts were used as features of gut microbiomes in both the Himalayan and the Hadza populations. Permutation tests with 10,000 permutations were performed to evaluate the significance of each CCA model and terms using “anova.cca” function in “vegan.” For all CCA models, the P values from the permutation tests were less than 0.05, indicating that the CCA model explained more variance of the gut microbiome in the Himalayan and the Hadza populations than expected by chance. American samples were excluded for both CA and CCA analyses.

Phylogenetic diversity was computed by rarefying the samples to various depths starting from 10–6,500 sequences per sample. Alpha diversity was measured using species richness, Shannon’s H, Simpson’s D, and Fisher’s alpha calculated as the mean values from 100 iterations at each depth. Kruskal–Wallis tests were used to assess the significance of differences in each of the alpha diversity metrics between populations at each rarefaction depth. Dunn’s posthoc test was performed to assess pairwise differences between populations. Differences in rarefaction depth did not alter significance of the observed differences. A generalized linear mixed effect model was used to evaluate associations between the 10 dietary and environmental factors and the two metrics of alpha diversity. Four models were created, each with the four metrics of alpha diversity (observed species, Shannon’s H, Simpson’s D, and Fisher’s alpha) as the response variables; the 10 factors were treated as explanatory variables with fixed effects; and each individual had random effect. Beta diversity was assessed using Bray–Curtis as well as unweighted and weighted UniFrac distances calculated by log transformation of the nonrarefied 16S count data. PERMANOVA was performed using the vegan package in R [68]. For all PERMANOVA analyses, 10,000 randomizations were performed to assess the statistical significance. In order to identify differentially abundant taxa at the phylum and genus levels, we first agglomerated the taxa abundance (counts) at each taxonomic level, respectively. The differences in taxa abundance (counts) were then assessed using the DESeq2 package [42].

SparCC was used to assess correlations between bacterial genera, as described previously [45]. SparCC is specifically designed to measure correlations in microbiome data and computes compositionality robust correlations by averaging multiple iterations of data. The statistical significance of the inferred correlations is then assessed using a bootstrap procedure. First, a large number of simulated data sets, in which all components are uncorrelated, are generated. Then, correlations are inferred from each simulated data set with the same parameter setting as is used for the original data. Finally, for each component pair, pseudo P values are assigned to be a proportion of simulated data sets for which a correlation value is at least as extreme as the one computed for the original data. We computed bacterial correlations for all pairs of genera after removing genera with less than 2 reads in at least 5% of samples (3 individuals and 124 genera). Correlations were computed from 100 iterations of the data, and we repeated the iterative procedure 100 times to compute the P values. P values < 0.05 after multiple testing correction were considered significant. Bacterial networks were visualized using the Fruchterman–Reingold force–directed layout algorithm implemented in the igraph package in R [69].

All multiple testing corrections were performed by computing FDRs using the Benjamini–Hochenberg method, and adjusted P values < 0.05 were considered statistically significant. A phyloseq object containing the 16S data and metadata as well as the analyses protocols used in this work are included in the supplementary data.