Probiotic consumption influences universal adaptive mutations in indigenous human and mouse gut microbiota

Probiotic intake commonly altered the genetic composition of gut microbial residents

To comprehensively understand the adaptive mutations in the resident gut microbiota due to probiotic intake, we first collected and curated publicly available metagenomic studies related to probiotics with the following criteria. (1) The study has a longitudinal design, which at least has a baseline and end time point for the probiotic consumption for a human or animal host subject. (2) The study does not use probiotics in combination with any other substance, such as medications, prebiotics, minerals, vitamins. (3) The study’s raw data were published and had detailed metadata. (4) The study’s sequencing data quality allows us to analyze at least species-level composition in the gut microbiome. (5) The study provided clear probiotics species/strains/product information. (6) The study has a clear statement on the dose and duration for probiotic intake. Finally, 11 high-quality metagenomic studies were included, among which seven were human cohorts (U.S., N = 1; Israel, N = 1; New Zealand, N = 2; China, N = 3), and five were animal cohorts (dog, N = 1; rat, N = 1; mice, N = 3). In total, 224 probiotic-treated individuals and 197 placebo controls (Tables 1 and 2) were included. The probiotic-administration duration for hosts in studies ranged from 1 week to 2 years and its median was 4 weeks. The median dose of probiotic administration was 9¹⁰ CFU/day, ranging from 10⁸ to 10¹⁰. Next, MetaPhlan2 was employed to identify the microbial compositions that have a relative abundance >0.5% for SNV profiling (Supplementary Data 1). The metagenomic reads were then mapped to the reference genomes of these selected species for SNV identification. We compare the SNVs against each reference genome for each host before and after probiotic treatment. A total of 16,901 SNVs were associated with probiotic administration (Supplementary Data 2). We first wondered how diverse were resident gut microbes that spontaneously mutated after probiotic consumption and if such diversity can be different from usual (the control group). Interestingly, the number of gut resident species occurring SNVs significantly decreased with hosts after the dietary intervention with Probio-Fit, Lactobacillus rhamnosus GG (L. rhamnosus GG) and Bifidobacterium lactis HN019, besides in the mice with Lactobacillus plantarum HNU082 (L. plantarum HNU082) (Wilcoxon rank-sum test, Fig. 1a). Next, raw SNV frequency might be not comparable across studies due to the inevitable sample/study-level disparity in the metagenome sequencing depth. Specifically, raw SNV frequency positively correlated with sequencing depth in both human and mice populations (Fig. 1b and Supplementary Fig. 1). To reduce this technical bias across studies, a sequencing-depth-normalized number of SNVs (nSNVs) was used for the following cross-study comparisons.

a The barplot indicates the number of gut resident species that had SNVs after probiotic intake. The bars represent the number of species that had SNVs. b The scatter plot shows the correlation between the raw number of SNVs and the sequencing depth of the target genomes. Green represents the mice (R = 0.66, P = 0.003), while orange represents the human cohort (R = 0.69, P < 0.001). c The barplot indicates the normalized number of SNVs identified in two time points (baseline and end point of probiotic intervention in each study). All error bars represent the SEM. The source data for graphs are available as Supplementary Data 7 or FigShare (figshare.com/projects/Probiotic_consumption_influences_universal_adaptive_mutations_in_indigenous_human_and_mouse_gut_microbiota/122447).

The consumption of probiotic L. plantarum HNU082, L. rhamnosus GG, and Bifidobacterium lactis HN019 significantly reduced the total frequency of SNVs (nSNVs) in the gut residents (Wilcoxon rank-sum test, Fig. 1c). Overall, these suggested that probiotic intake can significantly change the genetic composition of a wide range of indigenous gut microbiota that was often not assumed.

The nSNVs introduced by probiotics consumption were strain-specific

We next compared the nSNVs before and after probiotic intake in each of the studies. Overall, probiotic intake caused more SNVs in gut microbiota than the control group without any probiotic consumption (Fig. 2a). Furthermore, the SNVs in mice native gut microbiome outnumbered that in humans. Next, alpha diversities (Shannon and Simpson index) and beta diversity were calculated for each sample based on the profile of species-level SNVs. PERMANOVA was used to measure the effect size of probiotic intake on the SNV profiles at the species level (p < 0.001) (Fig. 2b−d and Supplementary Data 3). Our results suggested that the overall pattern of SNVs induced by probiotics was highly specific to probiotic strains. Furthermore, there is no significant correlation between nSNVs and experimental factors such as probiotics dose and duration of probiotics, observed from our investigation (Fig. 2e).

a A scatter plot shows the normalized number of SNVs (i.e., SNV number/sequencing depth of the target genome) induced by probiotic intervention for all included studies. Dashed lines connect the same probiotic strains in the animal experiment and human cohorts. The blue shadow represents the comparison of nSNVs produced by probiotics with that of the control group (NT-HNU082 and NT-Supherb Bio-25). b, c In human cohorts, alpha diversity indexes (Shannon and Simpson index) were calculated based on the species-level profile of SNVs induced by probiotic intake (P < 0.001). The shaded area indicates the full probability distribution of the variable. d In human cohorts, the genetic beta-diversity difference within and between studies was estimated based on the species-level SNV profile of native gut microbiomes. e No significant correlations between duration (R = 0.134, P = 0.205) and dose of probiotics (R = −0.112, P = 0.289) and the normalized number of SNVs induced by probiotics. f The study-dependent correlation between the Bray−Curtis distance of microbial species abundance profiles and the normalized number of SNVs induced by probiotics, including positive (AH1206 and HNU082), negative (LGG and Zhang), and uncorrelated (mixed probiotics). The source data for graphs are available as Supplementary Data 7 or FigShare (figshare.com/projects/Probiotic_consumption_influences_universal_adaptive_mutations_in_indigenous_human_and_mouse_gut_microbiota/122447).

To mitigate the potential effect of confounding factors, such as individuality in the gut microbiome, in our analyses, six probiotic studies were focused, including Bifidobacterium longum AH1206 (B. longum AH1206)⁷, Supherb Bio-25 ¹⁹, L. rhamnosus GG²⁰, Probio-Fit²¹, L. plantarum HNU082 ²² and Lactobacillus casei Zhang (L. casei Zhang) where host participants had paired/repeated microbiome measurements before and after the probiotic intervention. The correlation pattern between nSNVs and the beta diversity of gut microbiota was highly specific to what probiotic strains had been consumed (Fig. 2f). We found that nSNVs caused by B. longum AH1206 and L. plantarum HNU082 consumption had a positive correlation with Bray−Curtis distance of gut microbiota between baseline and post-probiotic intervention (B. longum AH1206, R = 0.43; L. plantarum HNU082, R = 0.607), while L. rhamnosus GG and L. casei Zhang had a negative correlation (L. rhamnosus GG, R = −0.346; L. casei Zhang, R = −0.367). No correlation between the nSNVs caused by mixed probiotics and Bray−Curtis distance of gut microbiota was found. These suggested that the overall pattern of nSNVs induced by probiotics was highly probiotic-strain-specific.

Universal adaptive mutations in indigenous gut microbes in response to probiotic intervention

We identified three bacterial gut residents that accumulated the convergent genetic changes in response to probiotic consumption in six human metagenomic studies, including Faecalibacterium prausnitzii (F. prausnitzii), Eubacterium rectale, and Roseburia intestinalis (Supplementary Data 2 and Fig. 3a−c). Interestingly, the probiotic interventions did not significantly alter the relative abundance of all these three species (Wilcoxon rank-sum test, p > 0.05), except for the LGG cohort (Wilcoxon rank-sum test, p < 0.05, F. prausnitzii and Eubacterium rectale). While ecological alterations in the gut microbiome were limited, the probiotic intervention led to widespread shifts in the genetic composition (detectable SNVs) of these individual gut residents (Fig. 4a, F. prausnitzii, p = 0.026). These suggested that evolutionary response might precede the ecological changes in the microbial communities under selection pressure.

The boxplots indicate the relative abundance and the normalized number of SNVs of a F. prausnitzii, b Eubacterium rectale, c Roseburia intestinalis in between time points in each study. Wilcoxon rank-sum test was used to compare the abundance or nSNVs between time points and considered the significance at 0.05 levels. For each comparison, T0 represents the baseline phase and T1 represents the end point of the intervention phase (time variable). The box represents the 25–75th percentile, whiskers represent the full range, and the line represents the median value. The source data for graphs are available as Supplementary Data 7 or FigShare (figshare.com/projects/Probiotic_consumption_influences_universal_adaptive_mutations_in_indigenous_human_and_mouse_gut_microbiota/122447).

a The details of 17 adaptive SNVs of F. prausnitzii. b The location of 17 unique SNVs in the genome of F. prausnitzii. All SNVs belonged to four types of mutations (A > G, A > C, G > A and G > T). c The crystal structure of Nitroreductase family protein and FprA family A-type flavoprotein. The red was the amino acid perssad. d The dN/dS ratios of genes with SNVs due to probiotics consumption (The dN/dS ratios can be performed when SNV has both non-synonymous and synonymous on the nine functional proteins.) The source data for graphs are available as Supplementary Data 7 or FigShare (figshare.com/projects/Probiotic_consumption_influences_universal_adaptive_mutations_in_indigenous_human_and_mouse_gut_microbiota/122447).

To investigate if different probiotic interventions can lead to similar genomic variations, candidate adaptive SNVs were explored, which can be commonly found in at least three out of six probiotics-intervention studies. Remarkably, F. prausnitzii ATCC 27768 had the most shared SNVs (N = 19) across independent studies (Supplementary Data 4), while Eubacterium rectale and Roseburia intestinalis also had two shared SNVs respectively. We next validated whether these candidate adaptive SNVs produced by probiotic intervention can also occur in the control group (null model, Supplementary Data 5). The four SNVs from Eubacterium rectale and Roseburia intestinalis can be also identified in Israel control cohorts (null model). Two SNVs from F. prausnitzii in the probiotics group were detected in the control group as well. Therefore, we pinpointed a total of adaptive 17 SNVs occurred in F. prausnitzii specifically adapted to probiotic intake and can be validated across distinct host cohorts (Fig. 4a).

Functional annotation of SNV-related genes of F. prausnitzii induced by probiotic intervention

Among those 17 adaptive SNVs due to probiotics consumption, 13 (76.5%) occurred in the gene coding regions of functional genes. Seven were non-synonymous mutations, while six were synonymous mutations. These mutations involved in nine functional proteins, including 30S ribosomal protein S5, phosphohydrolase, sensor histidine kinase KdpD, ferritin, fprA family A-type flavoprotein, nitroreductase family protein, ribonucleotide-diphosphate reductase subunitbeta, peptidase S24 and Type II toxin-antitoxin system PemK/MazF family toxin (Fig. 4a and Table 3), including four types of mutations A > G (n = 6), A > C (n = 6), G > A (n = 3) and G > T (n = 2) (Fig. 4b and Table 3). Given six protein-expressing genes contained non-synonymous mutations. Next, Phyre2 was employed to predict the protein structure before and after probiotic intake and further visualized how these non-synonymous genetic mutations significantly changed the protein structure via EZMOL. The predicted structure of nitroreductase family protein and fprA family A-type flavoprotein has been substantially modified (Fig. 4c), suggesting significant changes in the functional potential of the gut microbiome after probiotic exposure. The structures and amino acid sequences of other proteins have been provided in Supplementary Fig. 2.

To investigate how differentially functional genes responded to the gut selective pressure due to probiotic intake, the ratio of non-synonymous and synonymous (dN/dS) was calculated. The dN/dS ratio < 0.25 indicated the purifying selection acting on the genes, while the ratio >1 suggests that a gene was under positive selection for adapting to a new and or changing habitat^23,24. In our study, the dN/dS ratios in different probiotic interventions ranged from 0.15 to 2.0 or from 0.25 to 1 (Fig. 4d). This suggested that different functional genes of a gut microbial strain can have diverse evolutionary trends. Moreover, the same gene may present parallel evolutionary trends under the different interventions of probiotics. Specifically, the dN/dS ratio of nitroreductase family protein was >1 in probiotics B. longum AH1206, L. plantarum HNU082, and L. casei Zhang group. Phosphohydrolase was positively selected during the probiotic treatment with both L. plantarum HNU082 and mixed probiotics (Probio-Fit). Also, the same dN/dS ratios pattern for mixed probiotics (Probio-Fit) and a single-strain probiotic (L. plantarum HNU082) was exhibited in peptidase S24 and type II toxin-antitoxin system PemK/MazF family toxin. Nonetheless, different probiotic products may still have distinct patterns of evolutionary effect on a microbial functional gene of gut residents. Under the intervention of probiotic strain L. plantarum HNU082, the dN/dS ratios of ferritin and fprA family A-type flavoprotein were >1, while the mixed strains intervention was the opposite. Notably, only one gene, sensor histidine kinase KdpD, was under purifying selection (dN/dS < 0.25). It suggests that most genes in F. prausnitzii tend to be neutral by the new gut environment shaped by the probiotic ingestion. The above results illustrated the distinct evolutionary changes in the intestinal microbiota under the environmental pressure of different probiotic interventions.

The heritability of adaptive SNVs induced by probiotic intervention

To investigate whether or how long such adaptive mutations accumulated in the key gut residents, such as F. prausnitzii, can be inherited, an independent longitudinal microbiome study of probiotic intervention was conducted using L. plantarum HNU082 as a model strain (Fig. 5a). All six human participants in this validation study successfully completed two experimental phases: (I) continuous probiotic intervention for 7 days; (II) a long-term follow-up microbiome study (6 months after phase I). They volunteered to provide stool samples throughout all experimental phases as requested. Firstly, we identified 610 SNVs of F. prausnitzii at the end point of phase I, while a total of 1828 SNVs genomes were identified at phase II. Among those 610 SNVs identified from phase I, 317 (51.96%) were transient mutations that were not detectable at phase II, while 293 (48.04%) were retained on the F. prausnitzii genome at phase II (Fig. 5b). These suggested that probiotic intervention led to long-lasting yet often overlooked genetic changes in the gut residents. In the 293 heritable SNVs and 317 transient SNVs we observed, 129 functional genes were identified. Within the 129 functional genes, 39 were uniquely from heritable SNVs, 43 were uniquely from transient SNVs, and 47 overlap (Fig. 5c and Supplementary Data 6).

a Our experiment design of independent longitudinal microbiome study of probiotic intervention using Lp082 (n = 6). b Venn diagram indicates the unique and shared SNVs of F. prausnitzii induced by probiotic HNU082 in phase I (right after probiotic intervention) and II (6 months after probiotic intervention). c Venn diagram indicates the unique and shared proteins which inherited SNVs and transient SNVs involved. d The barplot indicates the SNV number of entirely SNV-inherited proteins. e The barplot indicates the SNV number of 20 protein products that had SNVs transiently occurred during the probiotic intervention. The source data for graphs are available as Supplementary Data 7 or FigShare (figshare.com/projects/Probiotic_consumption_influences_universal_adaptive_mutations_in_indigenous_human_and_mouse_gut_microbiota/122447).

We next characterized the functional genes with entirely inherited or transient SNVs induced by probiotic intervention from phase I to II. Sixteen entirely SNV-inherited proteins were identified firstly (Fig. 5d), which contained at least two consistent SNVs at both phases I and II. We next functionally annotated 20 protein products that have at least two transient SNVs at phase I whereas these two SNVs were not detectable at phase II (Fig. 5e and Supplementary Data 6). For example, one of those entirely SNV-inherited proteins, FprA family A-type flavoprotein, possesses ten SNVs induced by probiotic HNU082 that can inherit in an extraordinarily long period. Intriguingly, most transient-SNVs-related proteins are involved in carbohydrate transport and metabolism, such as carbohydrate ABC transporter permease, carbohydrate ABC transporter substrate-binding protein and carbohydrate-binding protein. These suggested that residents in the gut microbial communities tended to adaptively evolve carbohydrate-related proteins for the short-term probiotic invasion.