Distribution of Enterobacteriaceae worldwide

To perform a comprehensive global characterization of human gut microbiome signatures linked to Enterobacteriaceae colonization, we retrieved 12,238 public human gut metagenomic samples from 65 studies across 45 countries (Fig. 1 and Supplementary Table 1). Samples were collected primarily from Europe (n = 4,284, 35%) and North America (n = 3,367, 27.5%), followed by Asia (n = 2,844, 23.2%) and Africa (n = 1,024, 8.4%). The majority of metagenomic datasets were from adults (n = 8,275, 67.6%) and healthy individuals (n = 7,606, 62.2%) (Extended Data Fig. 1a).

Fig. 1: Exploring the global ecological landscape of Enterobacteriaceae.

a, Geographic distribution of the 12,238 human gut metagenomic samples used in this study. b, Workflow developed to identify and functionally characterize Enterobacteriaceae co-excluders and co-colonizers. ML, machine learning; GSMM, genome-scale metabolic modelling; BGCs, biosynthetic gene clusters. c, Metadata distribution of the number of samples where no Enterobacteriaceae species was detected (Absent) or at least one species was detected (Presence).

Using the Unified Human Gastrointestinal Genome (UHGG) catalogue12, we employed a mapping-based approach to accurately detect the presence and abundance of 4,612 gut microbial species (113 from the Enterobacteriaceae family) in the 12,238 global metagenomes on the basis of their level of breadth, depth and expected coverage (see Methods); the chosen thresholds were further validated with an experimental mock community (Extended Data Fig. 1b). Applying these parameters to synthetic metagenomes showed the detection limit of our metagenomic approach to be within a relative abundance of 0.003–0.01% (Extended Data Fig. 1c). The overall prevalence observed for Enterobacteriaceae was 66%, which is in line with a previous culturing-based surveillance study of Escherichia coli found in stool13. The distribution of Enterobacteriaceae was generally well balanced across metadata categories (Fig. 1c) and studies (Extended Data Fig. 1d), with a median prevalence of 57% across continents, 67% across age groups and 71% across different health states. The genera Escherichia, Klebsiella and Enterobacter were the most prevalent across various age groups, health states and continents, with E. coli, Klebsiella pneumoniae and Enterobacter hormaechei representing the most frequent species in their respective genera. The prevalence of E. coli was highest among African samples (88%), infant metagenomes (74%) and rheumatoid arthritis patients (96%) (Fig. 2a).

Fig. 2: Distribution and diversity of the most prevalent Enterobacteriaceae species.

a, Prevalence and median abundance of representative species from the five most prevalent Enterobacteriaceae genera across different age groups, continents and health states. b, Upset plot showing Enterobacteriaceae co-colonization patterns. Vertical bars represent the proportion of samples by continent harbouring the species highlighted in the lower panel. Numbers below the bars indicate sample size. Horizontal bars in the lower left panel show the total number of samples in which each species was detected.

We further investigated the co-distribution of Enterobacteriaceae species across this gut metagenomic collection to characterize patterns of polymicrobial colonization (Fig. 2b). Although most samples were found to be uniquely colonized by E. coli, we identified a statistically significant co-colonization of E. coli with K. pneumoniae predominantly among samples from Asia (observed vs expected proportion of 16.2% vs 10.9%, binomial exact test, P < 0.001). Moreover, we detected a significant co-occurrence of E. coli, K. pneumoniae and Enterobacter hormaechei primarily in samples from Africa and Oceania (observed vs expected proportion of 5.3% vs 0.7%, binomial exact test, P < 0.001). Overall, these distinct geographic patterns of Enterobacteriaceae co-colonization might be reflective of variations in environmental conditions, dietary habits, lifestyle and/or healthcare practices.

As E. coli is the most prevalent Enterobacteriaceae in the gut, we performed a dedicated analysis of the strain diversity of this species. To reduce the effect of host state and better understand the subspecies diversity circulating asymptomatically in the human population, we focused on 5,128 metagenomic samples collected from healthy adults. Through a metagenomic multilocus sequence typing (MLST)14 we identified 585 distinct E. coli sequence types (STs) (Extended Data Fig. 2a). The most prevalent known STs were classified as ST10, ST95, ST131 and ST73, which represent dominant E. coli lineages worldwide13,15,16. However, 76.5% of detected strains belonged to unknown STs, including two STs (here labelled as 100024 and 100083) found to be among the 10 most frequent STs (Extended Data Fig. 2a). These unknown lineages were found to be overrepresented particularly among samples from Africa (Extended Data Fig. 2b). Given current reference biases towards E. coli clinical isolates and the extent of unknown subspecies diversity uncovered here, these results suggest that a substantial global diversity of E. coli remains uncharacterized.

Microbiome structure linked to colonization dynamics

Having access to a global collection of >12,000 human gut metagenomes enabled us to explore the relationship between the gut microbiome composition and Enterobacteriaceae colonization status. First, we built machine learning classifiers to distinguish samples with or without Enterobacteriaceae on the basis of the abundance and prevalence of the remaining non-Enterobacteriaceae microbiome species (Fig. 3a and Extended Data Fig. 3a). We tested three supervised learning methods (ridge regression, random forest and gradient boosting) using as outcome variables the colonization status of Enterobacteriaceae as a whole, or that of E. coli and K. pneumoniae in particular. Results across all methods and variables showed a consistently good performance (median area under the receiver operating curve, AUROC = 0.788), with gradient boosting outperforming ridge regression and random forest for the classification of Enterobacteriaceae (median AUROC = 0.812), E. coli (median AUROC = 0.797) and K. pneumoniae status (median AUROC = 0.773). Model performance was consistent between the entire metagenomic dataset and when only considering samples from healthy adults (Extended Data Fig. 3a). Given the variation in gut microbiome composition between geographic regions, we also investigated whether models were generalizable within and across different continents. Focusing on samples from healthy adults, models tested on a per continent basis had overall good performance (AUROC > 0.7 for all continents tested, Extended Data Fig. 3b), with the highest performance observed for samples from Africa. Cross-validation between continents also revealed that models specifically trained on metagenomes from Asia, Europe or North America performed well across other regions (AUROC > 0.7 in at least 2 other continents, Extended Data Fig. 3c), showing that models trained on larger sample sizes are more generalizable. Overall, these results indicate that the human gut microbiome harbours compositional differences linked to Enterobacteriaceae colonization, even across different health states and geographic locations.

Fig. 3: Gut microbiome composition is associated with Enterobacteriaceae colonization and abundance.

a, ROC curve of the machine learning results linking the gut microbiome composition with Enterobacteriaceae, E. coli or K. pneumoniae colonization status. AUROC values were obtained with gradient boosting applied to the 12,238 human gut metagenomes. b, Phylogenetic tree of the 306 bacterial species associated with Enterobacteriaceae colonization and abundance. Clades are coloured according to their affiliated order. Red and blue colours in the outer layer indicate the number of analyses (out of 12) in which each species was classified as a co-excluder (negative) or co-colonizer (positive). A maximum score of 12 denotes that the species was found to be consistently associated with Enterobacteriaecae, E. coli and K. pneumoniae co-colonization across all 12,238 datasets, as well as the subset of healthy adults. c, Top 10 gut microbiome species classified as co-excluders (negative effect size) or co-colonizers (positive effect size) coloured by their family affiliation.

We next performed diversity analyses to investigate overall community differences relating to Enterobacteriaceae colonization and abundance (Extended Data Fig. 4). Beta diversity estimates revealed higher pairwise distance among samples with Enterobacteriaceae compared with those without (Wilcoxon rank-sum test, P < 0.0001, Extended Data Fig. 4a). These differences were independent of alpha diversity and sample read depth, as we observed a low correlation between Enterobacteriaceae abundance and Shannon diversity estimates (Pearson’s coefficient of determination, R2 = 0.03, Extended Data Fig. 4b), even after subsampling to 500,000 mapped reads—the minimum depth needed for accurate beta and alpha diversity estimates17.

To identify the microbiome species associated with Enterobacteriaceae presence (co-colonizers) or absence (co-excluders), we performed a differential abundance analysis on the basis of the intersection of a generalized (ALDEx2 (ref. 18)) and mixed-effects model (MaAsLin2 (ref. 19)), while accounting for study, age group, health state, continent and read depth (see Methods for further details). In addition, beyond investigating the differential abundance of species according to Enterobacteriaceae colonization status at a binary level (presence/absence), we used a network-based approach20,21 to model Enterobacteriaceae–microbiome co-abundance patterns. We analysed all 12,238 human gut metagenomes, as well as the subset of 5,128 samples from healthy adults to further control for microbiome differences related to age and health state. This revealed 307 prokaryotic species (12% of prevalence-filtered species) significantly associated with Enterobacteriaceae, E. coli and/or K. pneumoniae colonization and abundance (Fig. 3b, Extended Data Fig. 5a and Supplementary Table 2): 172 were identified as Enterobacteriaceae co-colonizers and 135 as Enterobacteriaceae co-excluders.

At a taxonomic level, species from the orders Lachnospirales, Oscillospirales and Bacteroidales were overrepresented among the co-excluders (Fisher’s exact test, adjusted P < 0.05). In contrast, co-colonizers were significantly associated with the orders Lactobacillales, Veillonellales and Actinomycetales. Analysis of the 1,000 most prevalent species revealed that 17 bacterial orders contained neither co-excluders nor co-colonizers (Extended Data Fig. 5b), even though two taxa in particular (RF39 and Burkholderiales) were represented by >10 species each. At a species level, 89% of the 307 candidate species showed a consistent signal across datasets (all or healthy adults only) and taxa (Enterobacteriaceae, or E. coli and K. pneumoniae individually) (Extended Data Fig. 5c), showing that the identified microbiome signatures are robust to differences in host state and analysis resolution. Species from the Faecalibacterium genus were among the strongest co-excluders (Fig. 3c and Extended Data Fig. 5d), with one uncharacterized Faecalibacterium species (Faecalibacterium sp.900539885) identified as the top antagonistic candidate. Previous studies have shown that species from the Faecalibacterium genus (for example, F. prausnitzii) carry important beneficial functions in the intestinal tract such as the production of short-chain fatty acids (SCFAs)22. This in turn has been shown to directly inhibit the growth of Enterobacteriaceae species, including E. coli and K. pneumoniae23. Thus, the overrepresentation of Faecalibacterium species in gut microbiomes with low levels of Enterobacteriaceae might be reflective of a genus-wide mechanism of colonization resistance. Focusing on the co-colonization patterns, we found that members of the Intestinibacter, Veillonella and Enterococcus genera were the strongest candidates. E. faecalis has been previously shown to promote the growth and survival of E. coli in vitro and in vivo through the production of l-ornithine24. However, the relationship between species of the Intestinibacter and Veillonella genera with Enterobacteriaceae has not been previously explored and may potentially underlie uncharacterized mechanisms associated with Enterobacteriaceae colonization success and outgrowth.

To further evaluate the clinical relevance of our findings, we compared our results with an independent study that investigated the longitudinal dynamics of carbapenemase-producing Enterobacteriaceae (CPE)10. This study tracked for up to 12 timepoints a cohort of CPE-positive individuals that were later decolonized, as well as CPE-negative household controls (n = 46 participants; 361 samples). Comparison of differentially abundant species between CPE-positive individuals and household controls showed a statistically significant overlap (χ2 test, P = 0.0071) in relation to those we detected to be associated with Enterobacteriaceae in general (Extended Data Fig. 6a). However, the overlap was not significant when comparing CPE-positive to CPE-negative individuals that were previously colonized. This aligns with findings from the original CPE study10, which noted that individuals recently colonized by CPE are still undergoing microbiome recovery. Nevertheless, we suggest that the overlapping co-excluder and co-colonizer species here discovered (Extended Data Fig. 6b) are not only involved in Enterobacteriaceae colonization as a whole but may also be related to colonization of CPE lineages in particular.

Strain-specific patterns of Enterobacteriaceae colonization

Our analyses revealed a strong association between the microbiome species composition and Enterobacteriaceae colonization patterns. However, subspecies (that is, strain)-level differences could also be linked to Enterobacteriaceae–microbiome interactions. We therefore investigated strain-specific signatures of 39 gut microbiome species that were identified as either co-colonizers or co-excluders of Enterobacteriaceae among healthy adults, and that were represented by at least 10 genomes with CheckM25 statistics of >90% completeness and <5% contamination within the UHGG. Using the Unified Human Gastrointestinal Protein (UHGP)12 catalogue, we characterized the accessory genome (that is, genes detected in <90% of conspecific genomes) of these 39 selected species to identify subspecies populations associated with Enterobacteriaceae colonization. A total of 213 accessory genes were significantly associated with Enterobacteriaceae colonization status (207 positively- and 6 negatively-associated genes; adjusted P < 0.05). These were distributed across 15 of the 39 species and overrepresented in functions involved in nucleotide transport and metabolism (Fisher’s exact test, adjusted P = 7.93 × 10−5), with the majority belonging to the species Ruminococcus B gnavus and Faecalimonas phoceensis (Extended Data Fig. 7). We further investigated the phylogenetic similarity of strains with the highest number of accessory genes associated with Enterobacteriaceae (top 10% strains; Fig. 4). This revealed a much stronger population structure associated with Enterobacteriaceae co-colonization among F. phoceensis strains (permutational multivariate analysis of variance (PERMANOVA), R2 = 0.74, P < 0.001) compared with R. gnavus genomes (PERMANOVA, R2 = 0.02, P < 0.001; Fig. 4). Interestingly, only 2 Faecalimonas species among the 307 candidate species were identified as Enterobacteriaceae co-colonizers. In contrast to other genera such as Veillonella and Streptococcus, which harboured >10 candidate species each, these results suggest that there is a more specific, strain-level association between the diversity of Faecalimonas and Enterobacteriaceae co-colonization.

Fig. 4: Faecalimonas phoceensis exhibits strain-specific co-colonization patterns.

Core-genome phylogenetic tree of 200 Faecalimonas phoceensis (left) and 665 Ruminococcus B gnavus genomes (right). Clades are coloured on the basis of whether the genome is a metagenome-assembled genome (MAG) or an isolate. First outer layer denotes genome geographic origin and the second layer highlights in red those genomes with the highest number of significant accessory genes (top 10%). PERMANOVA was used to relate each species phylogenetic structure (pairwise cophenetic distances) with the number of significant accessory genes. Only genomes with >90% completeness were included in the analysis.

Co-colonizers are more functionally diverse

Complementing our taxonomic results, we analysed the functional capacity and diversity of the 245 (out of 307) non-Enterobacteriaceae species that consistently showed co-colonizing or co-excluding patterns across datasets and taxa. On the basis of the annotation of protein-coding sequences with KEGG26, we found that co-colonizers exhibited a greater functional diversity compared with co-excluders (Wilcoxon rank-sum test, P = 0.0049; Fig. 5a), which was found to be independent of annotation coverage and genome quality differences (Extended Data Fig. 8a). A previous study showed that metabolic independence drives gut microbial colonization and resilience of disease-associated species, particularly under inflammatory conditions27. Here we show that co-colonization of gut bacteria with Enterobacteriaceae is also associated with higher metabolic independence, even among healthy individuals. Statistical analysis of the KEGG Orthologs (KOs) associated with co-colonization or co-exclusion showed that functions involved in drug resistance and DNA regulation (for example, major facilitator superfamily transporter, 16S rRNA methyltransferase and HTH-type transcriptional regulators) are among the functions most strongly linked to co-colonization (Extended Data Fig. 8b and Supplementary Table 3). In contrast, genes related to iron metabolism and transport (for example, rubrerythrin, ferredoxin and ferrous iron transport protein B) show stronger association with colonization resistance. As iron is an essential nutrient for many human pathogens28, an enrichment of iron-utilizing genes among co-excluders suggests that there is competition for iron availability between Enterobacteriaceae and co-excluder species.

Fig. 5: Functional differences between co-excluders and co-colonizers.

a, Distribution of Shannon diversity values obtained among co-excluders (n = 129) and co-colonizers (n = 116), which did not belong to the Enterobacteriaceae family, based on the pattern of KOs detected per genome. Exact P values were calculated with a two-sided Wilcoxon rank-sum test. b, COG functional categories significantly associated with co-colonizers (positive effect size) or co-excluders (negative effect size). c, Primary metabolic pathways detected with gutSMASH differentially abundant between co-excluders and co-colonizers. d, Pairwise metabolic distances between co-excluders or co-colonizers compared to all Enterobacteriaceae species detected at >1% prevalence (co-excluders: n = 4,773 comparisons; co-colonizers: n = 4,292 comparisons). P values were calculated with a two-sided Wilcoxon rank-sum test. In a and d, box lengths represent the IQR of the data, the central line represents the median, and the whiskers depict the lowest and highest values within 1.5× the IQR of the first and third quartiles, respectively.

On a broader level, we characterized the main functional categories (Clustered Orthologs Groups, COGs) associated with Enterobacteriaceae co-colonization and co-exclusion (Fig. 5b). In general, functions involved in metabolism (for example, amino acids, nucleotides and inorganic ions) were enriched among co-colonizers (Fisher’s exact test, adjusted P < 0.05), which further supports the hypothesis that co-colonizers exhibit greater metabolic independence. On the other hand, co-excluders encoded a higher number of genes involved in signal transduction mechanisms, which includes functions such as sporulation, motility and quorum sensing. We also identified more genes with unknown functions among co-excluders, suggesting that they may carry uncharacterized mechanisms with potential roles in Enterobacteriaceae colonization resistance. These results were consistent even when only considering species from the Bacillota phylum, which contains the highest number of both co-excluders and co-colonizers (Extended Data Fig. 8c).

We further investigated the distribution of specialized primary metabolic pathways associated with Enterobacteriaceae colonization patterns using the gutSMASH29 algorithm. We found that co-excluders were primarily enriched in metabolic gene clusters involved in the production of the three major SCFAs (acetate, propionate and butyrate), as well as in other pathways involving the Rnf complex and 2-hydroxyglutaryl-CoA dehydratases (Fig. 5c). These results reinforce the genus-wide signal we observed for Faecalibacterium, one of the notable short-chain fatty acid producers in the gut. In addition, the Rnf complex signature further supports the importance of iron among co-excluders, as the Rnf complex consists of a ferredoxin:NAD+ oxidoreductase involved in energy production in anaerobic bacteria30. These analyses reveal that species co-occurring with Enterobacteriaceae are more functionally diverse, with iron metabolism and SCFAs potentially playing important roles in regulating gut environmental conditions and modulating Enterobacteriaceae colonization and abundance.

Gut colonization is largely driven by habitat filtering

To investigate the relationship between co-colonization patterns and interspecies metabolic interactions between Enterobacteriaceae species and the rest of the gut microbiome, we generated genome-scale metabolic models31 of all candidate co-colonizers and co-excluders (only considering those with a consistent signal across datasets and taxa), together with all Enterobacteriaceae species from the UHGG detected at >1% prevalence (n = 282). Using a phylogenetically adjusted quantification method32, we calculated metabolic competition and complementarity indices between all pairwise Enterobacteriaceae–microbiome combinations. Metabolic competition was calculated on the basis of the overlap between two given metabolic networks, while metabolic complementarity measured the potential of one species to utilize the metabolic output of another. Values were distributed across two main clusters based on high or low competition/complementarity indices (Extended Data Fig. 9a). However, within each cluster, there was a significant negative correlation between metabolic competition and complementarity (Cluster1: Pearson’s r = −0.85, P < 0.0001; Cluster2: Pearson’s r = −0.20, P < 0.0001). We combined both parameters to estimate a metabolic distance score and observed that co-colonizers showed a lower metabolic distance to Enterobacteriaceae species compared with co-excluders (Wilcoxon rank-sum test, P < 0.0001; Fig. 5d and Extended Data Fig. 9b). These results indicate that habitat filtering, a process that favours the coexistence of functionally similar species, is probably the main driver of colonization success and microbiome assembly among Enterobacteriaceae-associated species. In addition, metabolic comparisons were made both within and between the groups of co-excluders and co-colonizers. These analyses revealed that within-group metabolic distances were smaller than between-group differences (Wilcoxon rank-sum test, P < 0.0001; Extended Data Fig. 9c), indicating shared niche preference. Importantly, results were consistent even when simulating metabolic models under various media compositions reflecting differences in diets (Extended Data Fig. 9d). This supports a previous study showing that gut microbiome species with similar nutritional requirements tend to co-occur across individuals33.

To further explore functional differences between co-excluders and co-colonizers inferred from metabolic models, we compared the number and types of metabolite predicted from uptake and secretion fluxes within each species population. By simulating the models under a rich gut medium (M3)34, we observed that the number of metabolites predicted to be secreted was significantly higher among co-colonizers (Wilcoxon rank-sum test, P < 0.0001, Extended Data Fig. 10a), further supporting that higher functional diversity and versatility is associated with Enterobacteriaceae co-colonization. Statistical analyses identified 7 and 77 candidate metabolites from uptake or secretion fluxes, respectively, as differentially abundant between co-colonizers and co-excluders (Extended Data Fig. 10b). With regards to estimated uptake fluxes, l-serine and indole were the most significant metabolites enriched in co-colonizers and co-excluders, respectively. Interestingly, dietary l-serine has been previously described to provide a competitive fitness advantage to Enterobacteriaceae under inflammatory conditions35, while indole has been shown to alleviate intestinal inflammation through modulation of the gut microbiome composition36. Indole is also recognized as a signalling molecule among indole-producing bacteria, such as E. coli37. Therefore, a higher uptake of indole among co-excluders could impair intercellular signal communication of Enterobacteriaceae. In terms of secretion, we observed most notably an overrepresentation of the metabolites undecaprenyl phosphate/undecaprenyl-diphosphatase and thymine among co-excluders, concomitant with a higher secretion of oxidized glutathione and β-alanine among co-colonizers. Undecaprenyl phosphate is involved in the biogenesis of the bacterial cell wall38, while thymine is essential for DNA synthesis, repair and bacterial growth. In contrast, detection of oxidized glutathione among co-colonizers may be indicative of adaptation to an environment with higher oxygen tension and oxidative stress. Lastly, β-alanine, which we found to be overrepresented among co-colonizers, was previously identified as significantly increased in Crohn’s disease patients39. Overall, these results highlight metabolic differences between co-excluders and co-colonizing species that may reflect differences in colonization and adaptation to distinct gut niches. This further supports that modulation of the gut environment, for instance through diet, may affect susceptibility to Enterobacteriaceae colonization.

Co-excluders encode gene families involved in quorum sensing

As the production of secondary metabolites can influence bacterial fitness and interspecies ecological interactions, we performed a dedicated analysis of biosynthetic gene clusters (BGCs) among Enterobacteriaceae co-colonizers and co-excluders using the antiSMASH40 prediction tool. The BGCs identified as most significantly overrepresented among co-excluders belonged to a class of cyclic lactone autoinducers (Fisher’s exact test, adjusted P < 0.05; Fig. 6a). Autoinducers represent signalling molecules that play a role in bacterial communication through quorum sensing41. This is in line with the COG results which showed an enrichment of signal transduction mechanisms among co-excluders. By investigating the taxonomic distribution of all 147 autoinducer BGCs detected, we find that the majority were harboured by species from the Lachnospiraceae family (97/147, 66%). Further grouping of the BGCs by their genetic similarity (>50% nucleotide identity over >50% coverage; Supplementary Table 4) revealed 110 unique BGC families, segregated between co-excluders and co-colonizers (Fig. 6b). We compared all genes from each family to experimentally characterized BGCs in the Minimum Information about a Biosynthetic Gene Cluster (MIBiG) database42 and found an overall low amino acid identity to known clusters (interquartile range, IQR = 25–30%; Fig. 6c), suggesting that these autoinducer BGCs represent uncharacterized sequences. Quorum sensing molecules have been previously implicated as a means for members of the gut microbiome to provide colonization resistance against external pathogens43. Therefore, an enrichment of these autoinducer BGCs among co-excluders indicates that they could play a role in controlling the colonization and abundance of Enterobacteriaceae in the human gut.

Fig. 6: Co-excluders harbour biosynthetic gene clusters involved in quorum sensing.

a, BGCs detected with antiSMASH that were found to be overrepresented among co-excluders or co-colonizers (two-sided Fisher’s exact test, adjusted P < 0.05). b, Network of all cyclic lactone autoinducer BGCs detected among co-colonizers and co-excluders. BGC nodes are linked if they share >50% nucleotide identity over >50% alignment coverage. c, Distribution of amino acid identity values obtained by comparing each autoinducer BGC family against the MIBiG database (n values indicated in parenthesis next to each BGC family represent the total number of alignments against the MIBiG database). Box lengths represent the IQR of the data, the central line represents the median, and the whiskers depict the lowest and highest values within 1.5× the IQR of the first and third quartiles, respectively.