Industrialization drives convergent microbial and physiological shifts in the human metaorganism

Abstract

Understanding how host lifestyle and industrialization shape the human gut microbiome and intestinal physiology requires multimodal analyses across diverse global host contexts. Here, we generate multivariate data from the Global Microbiome Conservancy cohort, including gut microbiome, IgA-sequencing, host genotyping, diet, lifestyle and fecal biomarker profiles, to investigate host-microbiome interactions across gradients of industrialization and geography. We show that industrialization is associated with homogenized microbial compositions, reduced microbial diversity, and lower community stability, independent of host confounders. We further show that industrialization is linked to elevated markers of gut stress, increased IgA secretion, and altered patterns of IgA-bacteria interactions. Finally, we show that microbiome-based disease predictors trained on industrialized populations lose accuracy in less industrialized cohorts, highlighting limited cross-population transferability. Together, our results suggest profound restructuring of host-microbiome interactions due to industrialized lifestyles, and emphasize the need for inclusive, globally representative data to improve translational microbiome applications across diverse human populations.

Figure 1 –. Multimodal sampling in the GMbC cohort to investigate factors shaping human metaorganisms worldwide

I. Collected samples and metadata from GMbC participants (n = 1,015). Gut shotgun metagenome, human genotype, and cultured isolate genomic short read data were generated. Levels of fecal biomarkers (IgA, IgM, calprotectin and chromogranin) were measured. Lifestyle and diet metadata were collected. J. Countries (n = 12) and sampling locations (n = 35). Locations are colored based on industrialism and urbanism status of recruited populations. The number of participants per location is shown. K. Distribution of age, BMI and sex parameters in the GMbC cohort. L. Distribution of pairwise geographic distances across sampled participants (Haversine distance, assuming a spherical Earth). M. PCA of human genotypes of the GMbC cohort (green circles) and the 1000 Human genomes project (plain points, colors represent superpopulations). N. Left plot: Factor analysis of mixed data (FAMD) applied to lifestyle metadata factors (see Methods). Participants are colored by country. Top right: same plot, participants colored by industrialism status. Bottom right: same plot, participants colored by urbanism status. O. Correlation between individual lifestyle factors and the first ten principal components (PCs) of the Lifestyle PCA. Statistically significant associations (adjusted p-val < 0.05) are shown in plain circles. P. Principal component analysis (PCoA) of diet metadata (UniFrac-based distances of FFQ data, see Methods). Participants are colored by country as in F, left plot.

Figure 2 –. Effect of industrialized and urban lifestyles on community-level variations in the global human gut microbiome

G. Phylogenetic diversity (Faith’s PD index) of metagenomes (species level) as a function of the first PC of the Lifestyle PCA. PC1 Lifestyle is used as a continuous proxy for industrialization (Fig. 1F&G). Participants are colored by country (left plot), industrialization status (top right) or urbanism status (bottom right). The blue line shows the linear regression. H. PCoA of unweighted UniFrac compositional dissimilarities (species level). Participants are colored by country, industrialism status and urbanism status as in A. I. Inter-individual compositional dissimilarities across quantiles of PC1 Lifestyle, used as a continuous proxy for industrialization (Fig. 1F&G), increasing from first to fifth quantile. Industrialized and urban populations exhibit increased homogeneity of bacterial compositions (quantiles #4 & #5). J. Cumulative explained variance (as measured by stepwise dbRDA, forward multivariate model selection) of species-level microbiome compositions by individual host factors. The baseline covariates (latitude, longitude, age, sex and BMI) and the first 10 PCs for Lifestyle, Diet and Genetics were included in the model. K. Variance partitioning of alpha diversity of taxa (species level) (Faith’s PD and Shannon diversity) and KEGG KO functions (Richness). Marginal and maximum variances (R²) of alpha diversity were calculated for the following factors: Lifestyle (10 first PCs), Diet (10 first PCs), Genetics (10 first PCs), Non-Additive Effects of Lifestyle and Diet, and Non-Additive Effects of Lifestyle, Diet and Genetics. Statistical significance is represented as plain square and circle symbols. Marginal and maximum effects are calculated using a nested model approach where the null model contains baseline host covariates (age, sex, BMI and geography), and the full model contains all tested host variables (see Methods). L. Variance partitioning of PCs of beta diversity for taxa (species level) (unweighted UniFrac and Aitchison distance) and KEGG KO functions (Canberra distance). The same approach as in panel E was used. Colors and symbols for statistical significance are as in panel E.

Figure 3 –. Impact of industrialized and urban lifestyles on individual taxonomic and functional features of the gut microbiome

G. Meta-analysis of country-level associations between gut bacterial species abundance and urbanism status (urban vs. rural). The “META-ANALYSIS” column shows statistical results of the cross-country meta-analysis performed using an inverse-variance weighted fixed-effects approach. Because Canada included only rural participants in our cohort, it was grouped with the USA under the label ‘North America’. Species are shown in rows. The top 25 species most enriched in each category are shown. H. Effect of host covariates on the association between bacterial species abundance and industrialization (PC1 Lifestyle) (x axis: species; y axis: effect size). Two models were compared – a first model with Lifestyle PCs alone as predictors, and a second model with host covariates (baseline covariates and Diet and Genetic PCs). Significant associations are shown as plain symbols. I. Variance partitioning of species abundance to calculate marginal and maximum variance of baseline, Lifestyle, Diet and Genetic factors, along with Non-additive Effects. The same statistical framework and figure design as in Fig. 2E&F were used. “Environment” is defined as the overall contribution of Lifestyle and Diet (see Methods). Features are grouped by majority factors that explain >50% of variance. Data for KEGG KOs is shown in Supp. Fig. 3C and Supp. Table 3. J. Ternary plots of relative marginal effects on species (left) and KEGG KOs (right) profiles. Marginal effects for KEGG KOs were calculated as for species shown in panel C (see Supp. Fig. 3). Relative contributions of Genetics vs. Environment vs. Non-additive effects of Diet, Lifestyle and Genetics are shown in left plots, while contributions of Genetics vs. Diet vs. Lifestyle effects are shown in right plots. K. Identification of lifestyle-associated KEGG KOs across countries. Lifestyle-associated KEGG KOs were identified using the following criteria: (i) significant association with either Lifestyle alone or Environment, based on marginal variances calculated as in panel C, (ii) effect size (marginal R2) for Lifestyle higher than for Genetics and Diet, (iii) significant association with industrialization status (indus vs. non-indus, q-value < 0.05), (iv) significant association (q-value < 0.05) with PC1 Lifestyle and (v) significant within-country meta-association with PC1 Lifestyle. L. KEGG module enrichment analysis reveals Cobalamin biosynthesis as enriched in the industrialized group (left panel). The analysis is based on the average country correlation of KEGG KOs with PC1 Lifestyle. KEGG modules enriched in the industrialized group and highlighted in purple are described in further panels. Panels for modules M00924, M00925, M00122 and M00741 show the distribution of within-country correlation coefficients between all individual KEGGs of these modules and PC1 Lifestyle. KEGGs with plain symbols are part of the set of Industrialization-associated KEGGs identified in panel E.

Figure 4 –. Decreased gut microbiome stability among populations living in industrialized societies

G. Co-abundance networks reconstructed from metagenomic profiles of non-industrialized (left) and industrialized (right) samples. Only significant compositionally-corrected abundance correlations with |coefficient| >= 0.3 are retained (edges). Each node represents a bacterial species. Co-abundance networks were reconstructed by matching the two groups for alpha diversity (n=460 individuals in each group, see Methods and Supp. Fig. 4). Networks controlling for age, sex, BMI and geographic distribution are shown in Supp. Fig. 4. H. Increased node degree among co-abundance networks of more industrialized populations. Normalized node degree was calculated at the aggregate level for both non-industrialized and industrialized groups (left panel, Wilcoxon test), as well as for each individual locality (n participants >= 30) (right panel, linear regression with PC1 Lifestyle, coefficients shown at the bottom right). Average position along PC1 Lifestyle was calculated for each locality. I. Natural connectivity measurement of co-abundance networks based on industrialization status. Nodes were progressively and randomly removed (x axis) to simulate perturbations and measure the stability of the networks. J. Toy schematic illustrating the concept of a lineage as a series of matched co-abundance modules tracked across ordered localities. Eight synthetic locality-specific co-abundance networks are shown for eight localities ordered by PC1 Lifestyle. The highlighted module (colored edges and filled nodes) depicts a lineage negatively associated with PC1 Lifestyle, with module size diminishing along the PC. Curved grey connectors link the centroids of the highlighted module. Edge color indicates correlation sign. K. Top-10 module lineages with positive/negative trends along PC1 Lifestyle. Localities (x-axis) are ordered by mean PC1 Lifestyle; lineages (y-axis) are labeled (see Supp. Table 4). Points show the presence and size of the module in a given locality. Modules were matched across adjacent localities using taxonomic compositional overlap (Szymkiewicz–Simpson), with one-to-one assignment (Hungarian) and a single-gap lookahead to accommodate for module absence. Colors indicate trend category: green = negatively associated with PC1 Lifestyle (markers of less industrialized contexts), purple = positively associated (markers of more industrialized contexts). Top 10 lineages per trend were selected by module size and persistence over localities. L. Backbone co-abundance networks for two representative lineages. For each lineage, we identify the largest module occurrence (peak size) and build a concise correlation backbone over this set of taxa. Node size reflects prevalence of the taxon across localities where the lineage is present; edges represent the mean pairwise correlation across those localities. To improve readability, loosely connected nodes were removed from the co-abundance backbone. Top: lineage negatively correlated with PC1 Lifestyle; bottom: lineage positively correlated with PC1 Lifestyle.

Figure 5 –. Industrialized and urban lifestyles are associated with elevated levels of fecal markers of inflammation and immune response

A. Fecal IgA and chromogranin levels based on industrialization status (Wilcoxon tests, p-val = 3.6e-18 and p-val = 4.2e-30, respectively). B-E. Distribution of fecal IgA, chromogranin, calprotectin and IgM levels across countries, sorted by average position along PC1 Lifestyle, in ascending order. F. Percentage of participants with intermediary (50 to 200 ug/g) or high (>200 ug/g) levels of fecal calprotectin, across localities. The red dash line indicates the mean percentage of participants with calprotectin > 50 ug/g across all localities. Countries with a significant enrichment of participants with moderate to high levels of calprotectin (> 50 ug/g) are shown with asterisks (logistic regression with a GLM and a binomial distribution; Cameroon: p-val = 0.04; Ghana: p-val = 0.02; Nigeria: p-val = 0.03). Countries are colored as in panels B-E. G-H. PCoA of species-level microbiome compositions (as in Fig. 2B), with participants colored by level of IgA and chromogranin. I. Variance partitioning of IgA, Chromogranin, Calprotectin and IgM levels. The same statistical framework as in Fig. 2 and 3 was used, adding PCs of microbiome compositions to the models (dark green) to account for and measure the contribution of the microbiome to fecal biomarkers. J. Identification of bacterial taxa associated with fecal IgA levels (left panel) and fecal calprotectin levels (right panel). Both GLM and Bayesian-based (BIRDMAn) approaches were used to detect taxa being consistently differentially abundant based on IgA or calprotectin. Two regression models were tested with each approach. In the first one, the model controls for baseline host covariates (age, sex, BMI and geography) (triangle symbols). In the second model, host Lifestyle, Diet and Genetic PCs are added to further control for the effect of host covariates (circle symbols). Top 10 bacteria with the strongest positive and negative effect sizes based on GLM runs are represented (See Supp. Table 5 for the full data).

Figure 6 –. Shifts in IgA-gut bacteria interactions in the gut microbiome of industrialized populations

F. Correlation of OTU-specific average IgA Coating Index (ICI) between individuals from industrialized and non-industrialized groups (Pearson correlation of log ICI values). Comparisons were made on OTUs with sufficient representation and balanced sample size between lifestyle categories (n = 166, see Methods). G. Phylogenetic distribution of average ICI values in industrialized and non-industrialized groups (blue and red color gradients, respectively). OTUs are colored by phylum. H. Comparison of IgA-coating between industrialized (n = 51 samples) and non-industrialized (n = 232 samples) groups, based on: the frequency of IgA-coated OTUs (ICI > 10) per individuals (left); the comparison of observed vs. expected IgA-coating events in the industrialized group, based on the frequency of IgA-coating events in the non-industrialized group (p-value calculated with the Poisson distribution) (middle); and the Kullback–Leibler (KL) divergence between the IgA+ and the total 16S relative abundance profiles (right). I. Correlation between IgA coating index or relative abundance and industrialization status of the host. Compositionally-corrected (CLR-transformed) abundance values are used. OTUs with sufficient representation across lifestyles were included in the analysis (n = 197). Correlations between ICI and industrialization were calculated with both PC1_Lifestyle (ICI ~ PC1_Lifestyle) and industrialization as a binary variable (ICI ~ industrialized/non-industrialized). The correlation between the CLR-corrected abundance of OTUs and industrialization were also measured (Abundance ~ PC1_Lifestyle & Abundance ~ industrialized /non-industrialized). The correlation between ICI and relative abundance of the OTUs is also shown (ICI ~ Abundance). *: p-val < 0.05; **: p-val < 0.01; ***: p-val < 0.001. J. Significant associations between ICI and industrialization shown for six different OTUs.

Figure 7 –. Low portability of microbiome-based biomarkers of health and disease

D. Distribution of Gut Microbiome Health Index (GMHI) across GMbC countries. GMbC participants do not have known chronic diseases. The GMHI was also calculated on an external validation cohort of inflammatory bowel disease (IBD) patients from the USA (in grey; see Methods). Median GMHI values are shown as black lines. Positive values indicate healthy microbiomes, while negative values indicate dysbiotic microbiomes. Median values for Tanzania, CAR, Cameroon, Rwanda, Ghana and Thailand is 0, suggesting an absence of information to call for healthy or dysbiotic microbiomes. E. GMHI values across admixture clusters. Admixture clusters with a minimum of 50 individuals were included in the analysis. F. Percentage of samples that do not contain any of the health marker taxa that are informative for the GMHI, across admixture clusters. The lack of GMHI accuracy observed in panels A and B is due to the absence of key maker taxa of health defined from cohorts of majority populations.