Universal adaptive evolution strategy of probiotic Lactobacillus plantarum to gut selection pressure from humans, mice and zebrash

Improving the probiotic engraftment in humans requires a thorough understanding of in vivo gut-adaptive strategy of probiotics in diverse contexts. Here, we exposed the probiotic Lactobacillus plantarum HNU082 (Lp082) to gut microbiota of healthy humans, mice, and zebrash for four weeks. Independent of host choice, Lp082 established and adapted the gut by acquiring highly consistent single-nucleotide mutations, which modulated carbohydrate utilization, and acid tolerance, and signicantly promoted its in vivo competitive tness. In turn, resident gut microbial strains, especially competing strains with Lp082 (e.g. Bacteroides spp. and Bidobacterium spp.), actively respond to Lp082 engraftment by accumulating 10-70 folds more evolutionary changes than usual. Human gut microbiota exhibited a higher ecological and genetic stability than that of mice. Collectively, the highly convergent adaptation strategy of Lp082 across host environments lay the foundation for leveraging the animal model for ex vivo engineering of probiotics for better engraftment outcomes in humans.


Introduction
Probiotics have been long used as gastrointestinal therapeutics under many health conditions 1 . Unlike other therapeutics, probiotics are live micro-organisms and typically leverage different strategies to adapt within the human gut under high selective pressure 2,3 , e.g., spontaneous adaptive mutations 4,5 . The autonomous gut-adaptive behaviors of probiotics can confer them su cient tness advantages to engage in interactions with gut residents and host factors 6,7 . On the other side, it also provided novel opportunities for leveraging gut selective forces for genetic engineering of probiotics in vivo 3 . Unfortunately, for most probiotic strains, these in vivo genetic processes are still poorly characterized, while they can be subject to a few known factors, e.g., the individualized genomic characteristics of probiotics, the complexity or speci c members in indigenous gut microbiomes, and host-ltering mechanisms 8,9 . Therefore, more studies are urgently needed to comprehensively characterize the adaptive evolutionary behaviors of individual probiotic strains under the diverse schemes (e.g., microbiome complexity, diet and host factors 6,7 ). Among these, the impact of host factors on probiotic adaptive behaviors have been challenging to investigate 9,10 . Until now, the interactions between probiotics and indigenous gut microbiota/host phenotypes have been mainly established using animal models rather than humans. A wide range of animal models (including mice, ies, and zebra sh) have been employed to investigate how an engrafted probiotic (such as Lactobacillus plantarum) adapted within the gut 11 , and modulated host physiology 12 . However, the applicability of these ndings from animal models to humans has not yet been validated. Each host model has a speci c selection mechanism by which maintains the resident microbial compositions 13 and exogenous microbes in the gut. Given the characteristics of each probiotic strain is highly individualized, it is less clear that the adaptive behaviors can be always host-dependent. But, if a particular probiotic applies a universal adaptive strategy across hosts, more novel strategies on genetic engineering of probiotics can be proposed for human use. It is especially intriguing if the tness advantages (adaptive mutations) of a probiotic gained from the mouse gut can be stably transferred to human hosts resulting in a better probiotic engraftment outcome. Unfortunately, no study has bridged these knowledge gaps, which impede translational research and hamper the performance improvement of probiotics in human populations.
On the other hand, the probiotics engraftment can turn to be ecological and evolutionary forces reshaping the indigenous microbial communities. Many metagenomics studies explored the ecological impact of probiotic ingestion on gut microbiota, yet marginal changes in the composition of gut microbiota have been noted 10,[14][15][16] . But it is under-reported that the genetic compositions in the human gut microbiome are constantly evolving under intrinsic (such as aging) or external environmental disturbance (such as diet) 8,17,18 regardless of observable ecological changes. It is likely that evolutionary responses buffered a variety of environmental changes and further exerted a long-term effect on the population dynamics of evolving species 17,19 . While probiotics were proposed to modulate the gut ecosystem for digestive health 3,6,20 , they can also have a notable impact on the in vivo evolutionary trajectories or functions of those residents regardless of whether probiotics colonize or transiently pass through the gut 14 . However, few studies have systematically assessed risks in the evolution of native gut microbiota under the selection of ingested probiotics.
To address these questions, we employed Lactobacillus plantarum HNU082 (Lp082), which was isolated from traditional fermented food 21 and whole-genome sequenced (PRJCA000348, PRJNA637783), as a model probiotic strain 22 . It has gained increasing attention as its conventional characteristics of probiotic Lactobacillus plantarum 12,22 and speci c functions such as hyperlipidemia prevention 23 and regulation of neurotransmitter secretion disorder. Here, we explored the effects of host-derived selection pressures (humans, mice, and zebra sh) on the genetic stability of probiotics and, in turn, its ecological and evolutionary impact on the indigenous gut microbiota using shotgun metagenomic sequencing and isolate resequencing methods.

Results
The adaptive evolution of probiotics within the gut of distinct hosts We employed the standard reference-based approach to explore the genetic changes of the consumed probiotics Lp082 under the in vivo natural selection in the gut of multiple hosts (humans, mice and zebra sh) (Fig. 1A). We sequenced the complete genome of this reference probiotic strain, including one chromosome and four plasmids. Next, we isolated this strain from the feces of hosts at different time points for whole-genome sequencing. To identify and quantify the putative genetic mutations (such as single nucleotide variants) in the host-adapted strains, we compared the genome of all isolates with the reference genome. In total, 109 bacterial strains of Lp082 were isolated from feces or intestine content of three hosts, out of which 77 isolates were from humans, 25 from mice, and only 7 from zebra sh for the duration of the whole experiment (Fig. 1B, 1C, Table S1). A total of 71 putative single nucleotide variants (SNPs) and 2 mobile genetic elements (Fig. S1B) were identi ed and annotated from genome sequencing data of in vivo-adapted probiotic strains in human, mouse, and zebra sh models, out of which only 22 SNPs could be experimentally veri ed using PCR (Fig. 1D). By contrast, under the in vitro condition, no SNP was annotated in Lp082 incubated in the de Man, Rogosa and Sharpe (MRS) agar within 28 days (Fig. S1A), suggesting these adaptive mutations only occurred during the gut passage of this probiotic.
We next sought to characterize in vivo adaptive mutations of Lp082 across host models. First, despite the radically distinct host selection pressures, we did not identify any host-speci c SNPs in the genome of Lp082. Furthermore, no between-host difference in the accumulated mutations, that is, average number of SNPs at the last time point was observed (human=9.93; mouse=10; zebra sh=9; p=0.926, Kruskal-Wallis test) as well as mutation types (Fig. 1C).
We next pro led and compared the temporal evolutionary dynamics of Lp082 over a four-weeks sampling period between the humans and mice (Fig 1E, 1F, S1C-1G). As for all hosts, the mutational frequency of Lp082 steadily increased in the initial two weeks of gut colonization, suggesting the tness advantages of lineages were carrying them relative to their ancestors. Since relatively few new SNPs occurred, the mutation rate tends to stabilize after the third week of colonization (Fig. 1C). Based on these adaptive mutations, we employed the phylogeny to identify ve lineages of the probiotic Lp082, which commonly emerged in humans and mice (Fig. 1E). We further sought to investigate how similar the evolutionary history of these co-existing probiotic lineages can be between two hosts (Fig. 1E, F). Overall, the temporal dynamics pattern of the probiotic lineages throughout the sampling period was highly conservative between humans and mice (p=1e-4, R=0.86, Mantel test based on Jaccard distance). Despite its occupation in the ecologically and spatially distinct niches, Lp082 acquired highly similar mutations in almost the same order and timing for the gut adaptation. At the initial stage of colonization, the sublineage E1 rstly emerged by acquiring a mutation (SNP41) on the gene encoding IgG binding protein (Gene 2601). In the following sampling period, this lineage co-existed with other descendants with a consistently high frequency in the Lp082 population (>50%) (Fig. 1E, F, Table 1), suggesting that Lp082 should persistently resist the host immunity for establishing and maintaining its ecological niche in the gut. Next, the lineages E1-A and E1-B appeared and maintained their proportion in the Lp082 population in the gut throughout the sampling period. The most notable adaptive mutations (SNP10, 47, and 32) from these lineages were involved in the genes (Gene1717, 0658 and 2804) encoding inner membrane protein response to acid pH, transcriptional activator for 3-phenylpropionic acid catabolism, and transpose, indicating Lp082 was actively developing the acid-tolerance capability in vivo by adaptive mutations. Around 7 days after probiotic ingestion, E1-B-1 emerged and carried the representative mutations (e.g., SNP 23) that can enhance the capability of rhamnose utilization (Gene 1257).
Interestingly, this sub-lineage vanished at the end of our sampling period, and survived slightly longer in mice (Day 28) than in human gut (Day 21), which is the only observable difference in the probiotic evolutionary history in the gut between two hosts. We next wondered if or what genes in the probiotic genome had a similar pattern of adaptive mutations across independent lineages or host subjects. Firstly, we focused on those genes that have at least two SNPs (and their distance on the genome is no more than 2000bp) in the independent lineages from distinct host subjects or host species, indicating that they underwent parallel evolution. Such adaptive mutations identi ed from a gene are not likely a consequence of the genetic drifts, indicating these genes are under natural selection. Among all 12 mutated genes, we identi ed ve such genes from all 109 isolates (encoding trehalose 6-phosphate phosphorylase, transposase, and ABC transporter ATP-binding protein) ( Fig. 2A). Remarkably, these genes carried consistent multiple mutations across most host subjects of humans, mice and zebra sh ( Fig. 2A, Table S3). Next, we identi ed the other seven genes as single-mutation genes involved in bacterial alpha-L-rhamnosidase 6 hairpin glycosidase, immunoglobulin G-binding protein, glutathione reductase, inner membrane protein response to acidic pH, and transcriptional activator for 3-phenylpropionic acid catabolism (Fig. 2B). Collectively, the presence/absence patterns of adaptive mutations in each of 12 genes are highly consistent within and between host species (Fig. 2B), strongly demonstrating a universal adaptive strategy of this probiotic strain under varying selection pressures.
We next tested if those SNPs conferred the tness advantage of Lp082 in the host gut using the in vitro model. We selectively validated the phenotypic changes of the host-adapted isolates related to the rhamnose utilization and acid tolerance (Fig. 2C, S2B). Firstly, we focused on a non-synonymous SNP in Gene 1257 (annotated as bacterial alpha-L-rhamnosidase 6 hairpin glycosidase) which was identi ed in 21 host-adapted isolates. Intriguingly, all these mutants exhibited a higher capability in utilizing rhamnose in vitro than the original strain. While those mutants rapidly responded to the rhamnose substrate in the medium, the reference strain did not grow in the early 3 hours after its inoculation, suggesting this SNP had conferred a signi cant tness improvement of Lp082 after the gut passage. Furthermore, mutational isolates from humans and mice exhibited a highly consistent growth pattern.
The human-derived strains just have a slightly longer period of exponential phase than did mouse-derived strains. Next, we tested the acid tolerance of 78 mutational isolates with a synonymous SNP in Gene 1717 (encoding inner membrane protein response to acidic pH) independently. Among those 78 mutants, 16 (20.5%) exhibited a signi cant higher survival rate under each of low pH conditions (p<0.05, pH=3.5, 2.5, 2.0 respectively) than the reference strain. Notably, the in vitro experiments of these isolates were conducted 9 month after their passage through the gut of hosts, suggesting that these gut-adaptive mutations and related phenotypic changes can be stably inherited and preserved in a relatively long time.
Compositional and functional alterations of resident gut microbiota in response to probiotic engraftment We next examined how the resident gut microbiota ecologically responded to probiotic engraftment. Compared to the placebo group, the ingestion of Lp082 introduced a notable uctuation in the taxonomic structure of indigenous microbial communities in both humans and mice (Fig. 3A, S3A, S3B). As for humans, the probiotic intake did not signi cantly alter the overall taxonomic and functional composition of gut microbiome in any time points compared to baseline (R^2=0.088, p=0.979; PERMANOVA, Fig. 3A). Such a low probiotic-derived impact on the microbial compositions of the human gut may result from a higher level of individuality or resilience in the human gut microbiota prior to ingestion. By contrast, we observed a strong impact of probiotics on the structure and function of mouse gut microbiota (R^2=0.28, p=0.005; PERMANOVA, Fig. 3A, S3C, S3D). Compared to the baseline, the structure of the mouse gut microbiome signi cantly changed on Day 7, 14 and 28 after probiotic ingestion (R^2=0.28, p=0.005 on Day 14; R^2=0.19, p=0.039 at Day 21; R^2=0.59, p=0.004 at Day 28; PERMANOVA). All these observations were con rmed by both Aitchison (Fig. 3) and Bray-Curtis distance metrics (Fig. S3A, S3B).
Next, we sought to identify metagenomic markers associated with probiotic intake. Speci cally, we de ned a metagenomic marker according to the following conditions: 1) it signi cantly changed over time (e.g., from Day 0 to 14) in the probiotic group; 2) it maintains stable abundance during the experiment in the placebo group; 3) it was signi cantly different in abundance at a given time point (such Day 14) between the two groups. In the human population, no organismal or functional markers were identi ed based on these criteria which may be due to large inter-individual variations in the gut microbiome. In the mouse model, Faecalibaculum rodentium and 10 species from Bacteroides (including Bacteroides caccae, Bacteroides ovatus, Bacteroides vulgatus, Bacteroides xylanisolvens, etc.) increased signi cantly, while other 12 species including Prevotella dentalis, Prevotella oris and Pseudomonas uorescens decreased sharply after probiotic ingestion (Fig. S3C). Meanwhile, we also identi ed a total of 21 metabolic pathways (including D-galacturonate degradation, pyruvate fermentation to acetate and lactate, starch degradation, tetrapyrrole biosynthesis, L-citrulline biosynthesis, glycogen biosynthesis and L-isoleucine biosynthesis) enriched in the time points during probiotic colonization (Fig. S3D).
In order to compare microbiome responses to probiotic colonization across host species, we here focused on 19 microbial species commonly found in both mammalian hosts, which taxonomically belong to the genus of Bacteroides or Bi dobacterium (Fig. 3C). We found that these 19 microbial species residing the human and mouse gut had distinct temporal dynamics in response to probiotics ingestion. This suggests that the ecological response of resident microbes to probiotic engraftment highly depends on the host or gut microbiota context (Fig. 3B).
Co-occurrence relationship between Lp082 and resident microbes in the gut of humans and mice We next constructed a co-occurrence network in each host species to identify the correlation between resident microbes and the ingested probiotic. The co-occurrence networks were constructed using SpiecEasi 24 based on the species-level taxonomic pro les derived from shotgun metagenomic sequencing data (Fig. 3C). In the human model, we identi ed 10 species whose abundances strongly correlated with that of Lp082: Alistipes negoldii and Odoribacter splanchnicus positively correlated with Lp082, while several Bacteroides spp. and Bi dobacterium spp. (Bacteroides intestinalis, Bacteroides salanitronis, Bacteroides uniformis, Bacteroides thetaiotaomicron, Bi dobacterium longum, Bacteroides xylanisolvens, and Bi dobacterium adolescentis) negatively correlated with Lp082. In the mouse model, we identi ed 4 species (Bacteroides caccae, Bacteroides vulgatus, Bacteroides ovatus, and Bacteroides xylanisolvens) strongly correlated with Lp082, and all the 4 species negatively correlated with Lp082. We found that no common resident gut microbes across host species showed consistent correlation with Lp082. Except for Bacteroides xylanisolvens, there were no common negatively correlated bacterial species found between the human and mouse models (Fig. C). While most microbial competitors (with the exception of Bacteroides xylanisolvens) of Lp082 in the mouse model were also identi ed in the human gut, none of them exhibited the same co-occurrence relationship in the network analysis.
Profound evolutionary changes within the resident gut microbiota in response to probiotic engraftment Probiotic intake can give rise to dramatic changes in the genetic compositions in the resident species, yet it was surprisingly underreported. Here, we studied the within-genome evolution of approximately 37 prevalent microbial strains in the human or mouse gut. Firstly, we observed a striking difference in the overall number of adaptive mutations in the resident microbiota between placebo and probiotics groups (Fig. 4A). The average number of mutations that occurred in the placebo group is only 1.92 (or 1.69) for humans (or mice), respectively. By contrast, probiotic intake gave rise to a dramatic change in mutation frequency in the resident gut microbes, which reached up to 16.90 (or 78.02) per species in the human (or mouse) model, respectively (Fig. 4B, Fig. S4A). It indicated that a resident strain can accumulate in average 4.33 (or 7.67) SNPs per day in humans (or mice), respectively, which was remarkable given the high inter-individual variation in the gut microbiota of both humans and mice. Conversely, the genome of Lp082 was relatively stable: only 8 and 10 SNPs were detected on Day 14 and 28 respectively (Fig. 4B). We next asked if these adaptive behaviors of residents might give rise to compositional or functional changes in the communities. In the mouse gut, the more mutations occurred in resident gut microbiota, the more pronounced shifts can be observed in the taxonomic or functional pro les after probiotic ingestion (Mantel test; Fig. S4B). This suggested that adaptive mutations could be a strong driving force reshaping the resident gut microbiota.
We next sought to characterize the distribution of these adaptive mutations among resident microbial strains. We found that the majority of mutations distributed in Bacteroides spp. and Bi dobacterium spp., which were mostly identi ed as ecological competitors of Lp082 in the resident gut microbiota of either humans or mice. We further demonstrated that the putative ecological competitors (such as Bacteroides spp. and Bi dobacterium spp.) proactively responded to dietary probiotics by accumulating signi cantly more adaptive mutations than non-competitors (i.e. other bystanders or facilitators) of this probiotic (p<0.05, Wilcoxon rank-sum test, Fig. 4C, S4A). Remarkably, in the mouse model, mutations occurred in probiotic competitors were one to two orders of magnitude greater than those in the related species in the placebo group or the probiotic group over 28 days. Therefore, we provided the clue about the potentially causal links between evolutionary and population dynamics in the resident gut microbiota.
We next characterized how evolutionary dynamics in resident gut microbiota differed between humans and mice. Firstly, the resident microbes in the mouse gut accumulated much more adaptive mutations than those in the human gut (Fig. 4C, D). For example, Bacteroides xylanisolvens is one of the overlapped members in the gut microbiota of both humans and mice, being negatively correlated with Lp082 in both hosts. After 14 days of probiotic ingestion, we identi ed on average 596 within-host SNPs in the mouse gut microbiotas, whereas only 15 SNPs identi ed on average for a human individual in the gut microbiota (Fig. 4D). These may suggest probiotics intake imposes more intensive selection pressure to resident gut microbiota in the mouse than does in the human. Next, we observed a substantial functional difference in the adaptive-mutated genes of the resident gut microbiota between human and mouse models (Table  S9). The human-associated mutations primarily occurred in the genes related to cAMP-binding proteins and conjugative transposon protein, whereas mouse-associated mutations involved in more diverse functions, such as the degradation and utilization of the complex carbohydrates (Beta-galactosidase, Phosphoglycolate phosphatase, 3-oxoacyl-[acyl-carrier protein] reductase, and Galactoside Oacetyltransferase), Tetracycline resistance element mobilization regulatory protein, Asparaginyl-tRNA synthetase, DNA primase, conjugative transposon protein, and Transposase. There was only one gene (traG) that was consistently mutated in both the human and the mouse models, which related to the conjugative transfer of plasmid RP4 (TraG-like proteins are essential components of type IV secretion systems 25 ). Collectively, even though the engrafted probiotic applied universal strategy in the distinct hosts, the indigenous microbiome response to its colonization was highly divergent, underscoring the importance of careful interpretation of ecological and evolutionary patterns identi ed in the mouse model for human-related translational studies.

Discussion
Probiotics are live microorganisms that, in su cient dose, confer a bene t to the host, and are widely used for improving digestive health. However, their genomes and functional traits can vary during administration due to the myriad selection pressures in the host luminal environment 26 . The gastrointestinal tract harbors a large collection of microorganisms. The resident gut microbes are typically well adapted to their host, including the immune system and uctuations in environmental resources over a relatively long period, resulting in higher competitive tness over any potential invaders 17 . Changes in environmental selection pressures, such as competitions for limited carbon sources, force new microbes to dramatically reshape their genome and functional characteristics in a short-time scale.
In the present study, we demonstrated how the probiotic Lactobacillus plantarum HNU082 (Lp082) can evolve within the human and mouse gut by acquiring heritable changes to its genome. We found that Lp082 accumulated single-nucleotide mutations that modulate carbon utilization (i.e. rhamnose) and acid tolerance in two mammalian host models. The signi cant improvement in rhamnose utilization can be attributed to a nonsynonymous mutation in the gene related to bacterial alpha-L-rhamnosidase 6 hairpin glycosidase. Furthermore, the ndings also indicated that rhamnose can serve as a strain-speci c prebiotic 27  Knowledge of normal genomic variations within and between host species can further allow us to develop new probiotics or modify existing strains that can engraft and colonize a human host, adding desired functions, safely, effectively, and reliably 6,7 . For example, it is translational promising if probiotics can be genetically engineered by means of experimental evolution in the animal hosts where in vivo tness can be improved under the human gut selection pressure 3 . While many animal models (e.g. mice, ies, and zebra sh) have been employed to study adaptive evolution of probiotics 11,12 , very few studies have applied animal-gut-adapted probiotics to humans. There are several reasons. Firstly, it remains challenging to conduct in-deep study of in vivo adaptive evolutionary behaviors of individual probiotic strains. The probiotic characteristics are highly strain-speci c, the in vivo evolutionary behaviors can be more divergent. Unfortunately, the Next, the in vivo adaptation can be shaped by multiple intrinsic and extrinsic factors and thus highly unpredictable. A systematic investigation is required in the diverse schemes of probiotic strains, microbiome diversity levels, diet, and host factors. Crook et al. compared the genome adaptation of E. coli Nissle (EcN) in the host gut with different microbiome complexities 3 . As a gut commensal, this strain accumulated fewer mutations in the "high diversity" gut microbiome, while its genome varied more in the "low diversity" microbiome. By contrast, we found that Lp082 maintained a much higher genetic stability than EcN regardless of the radical differences in the microbiome complexity between human and mouse gut. Other than that, its adaptive evolutionary changes were not differentiated by the host factors from humans and mice. We argue that the universal adaptive strategy of Lp082 against selection pressures from multiple hosts could suggest the transferability of the tness advantages of this probiotic from mice to humans, which de nitely merits the further investigation.
Probiotics adapt and evolve in vivo in order to survive. They also in uence selection pressures on resident microbial strains. Here we systematically investigated the ecological and evolutionary impacts of a probiotic Lactobacillus plantarum on the resident gut microbiome in humans and mice. First, we observed a minimal ecological change (either compositional or functional changes) in the human gut due to probiotic ingestion which is consistent with previous studies 15, 28 . By contrast, the resident population of gut bacteria can rapidly and extensively evolve within 3-7 days after probiotic administration in both humans and mice. Most within-genome changes persisted for 28 days in our experiment. In the placebo group, very few evolutionary changes in resident microbes were found over time. Notably, strains competing with probiotic Lp082 account for the largest proportion (Human:74.36%; Mouse:77.68%) of evolutionary events in the whole resident gut microbial community. As described by the "Red Queen hypothesis" 6 , microbial competitions within our gut ecosystems can force the acceleration of microbial evolution possibly without any observable ecological changes, resulting in apparent stability. It also suggested that the common or even daily supplement of probiotics could rapidly become a strong driving force of evolution and ecology in the resident gut microbiome within a relatively short timescale, which has been largely overlooked. Given the substantial number of adaptive mutations accumulated, it is still challenging to identify any worrying ones in the abundant resident gut microbes, which have potential to alter the gut ecological features in the long-term run. Therefore, the ecological and evolutionary impact of dietary probiotics on gut microbiota should be seriously considered prior to administration 29 . Accordingly, future work should carefully re ect the limitations in the present evaluation system of probiotic e cacy and safety in order to better design a framework for systematic assessment of interactions between probiotics, gut microbiome, and host over time.

Methods And Materials
The experimental design In this study, we used Lp082 as a model probiotic strain to explore the effects of host-derived selection pressure (the mouse and human models were applied, the zebra sh model was used for further veri cation of adaptive mutations) on the genetic stability of the consumed probiotics and the impact of probiotic mutations on the indigenous gut microbiome of different hosts (the mouse and human model). First, we sequenced the complete genome of this model probiotic strain, including one chromosome and four plasmids. Secondly, we isolated the probiotics from the feces of hosts at different time points to identify genetic mutations using whole-genome resequencing. Simultaneously, the original strain was continuously inoculated in vitro and sequenced to assess genetic mutation in the absence of host selective pressure. Next, we employed the metagenomic sequencing method to characterize temporal dynamics of the abundance of probiotics strain and other gut-microbiota members and genetic variations in the community level after probiotic ingestion and con rming the impacts by comparing the results with that in the placebo group. Last but not least, the difference in genetic variations of Lp082 and its impact on the indigenous gut microbiota among hosts were highlighted in this study.
In the mice (C57BL/6, 5 weeks age) experiments, each animal was housed in a single cage. Room temperature was 26℃ and the padding was replaced once every day. All mice were divided into two groups, the probiotic group (n=6), and the placebo group (n=5). For the mice in the probiotic group, Lp082 (about 4 x 10 8 CFU/g with the fodder) was infused daily for 7 days. Fresh feces were collected from each cage for bacterial isolation on days 3, 7, 14, 21, and 28 after stopping probiotic feeding. Three pieces of fresh excrement were placed in a 7 mL sample tube, 4.5 ml saline (0.85%) was added after sterilization and then homogenized with homogenizer. Then the fecal samples were diluted and coated for Lp082 isolation as well as for shotgun metagenomic sequencing. For the mice in the placebo group, the same feeding method was performed only without the probiotic infused. The fresh feces were collected from each cage on days 0, 14 and 28 for shotgun metagenomic sequencing.
For human participants, each individual was informed of the experimental guidelines and details and consent obtained;12 volunteers agreed to participate in the experiment. They were randomly divided into 2 groups including the probiotic group (n=7, 4 females) and the placebo group (n=5, 2 females). During the experiment, the subjects were asked to avoid ingesting any probiotic product or antibiotic and to maintain their regular diet. All healthy participants nished the whole experiment including 6 females and 6 males aged from [18][19][20].49), they did not have in ammatory bowel disease or diabetes and had not used antibiotics for at least 3 months prior to sampling. For the seven volunteers in the probiotic group, they were asked to consume vacuum freeze-drying Lp082 powder (including 7 x 10 9 CFU live strains) 2g every day for 7 days 30 , and fecal samples were collected on days 3, 7, 14, 21 and 28 after stopping probiotic consumption. Fresh feces were placed in a sterile fecal sampling tube, then the diluted feces samples spread to the MRS agar plate for Lp082 isolation as well as for shotgun metagenomic sequencing. In the placebo group, the ve volunteers were requested to maintain a regular diet during the whole experiment, and their fresh feces were collected on days 0, 14, and 28 for shotgun metagenomic sequencing. The study was reviewed and approved by the Ethics Committee of Hainan University (HNU-2018037, Haikou, China), and informed consent was obtained from all volunteers before they enrolled in the study. The participants provided written informed consent to participate in the study. Sampling and all described subsequent steps were conducted in accordance with the approved guidelines.
Additionally, a zebra sh model was used to validate our ndings related to the impacts of host intestinal selective pressure on the genetic stability of the ingested probiotic. A total of 20 sh (15 weeks age) tanks (15-18 zebra sh in each tank) were used in this experiment. The water was changed and the sh were fed at 9: 00 am daily, the capacity of feed was calculated according to 3% of the body weight of each sh per day. After the adaptation period, Lp082 of 10 8 CFU/g was fed for 7 days with the fodder 23 . The intestines of 5 sh were taken and homogenized with 1 mL saline (0.85%) in 1.5 mL EP tube and in 3 days, 7 days, 14 days, 21 days and 28 days after stopping probiotic feeding. Then the intestine samples were diluted and coated for Lp082 isolation.

Isolation and con rmation of the ingested probiotic strain (Lp082) in the feces
The reference strain Lp082 used in the present study was rst isolated in traditional fermented food in Hainan province of China 21 and further certi ed as a probiotic due to the common characteristics of probiotic Lactobacillus spp. and speci c functions such as hyperlipidemia prevention 23 and neurotransmitter secretion disorder regulation. We obtained the complete genome of the strain with functional annotations in our previous research (PRJCA000348, PRJNA637783). By comparing with other genome sequenced Lactobacillus plantarum strains (Table S4), we identi ed strain-speci c primers of Lp082, which was used for genetic con rmation of isolates in the next step. The detailed experimental steps for Lp082 isolation from feces and con rmation were as follows: (1) Feces/intestinal contents were taken from the sample tube, and sterilized saline (0.85%) was added to the sample tube. After stirring with a homogenizer, the samples were diluted and coated. The mixture of 0.5 mL and 4.5 mL of 0.85% saline (NaCl) was recorded as 10-1, 10-2, 10-3, 10-4, 10-5, 10-6 were obtained.
(2) Each diluted gradient was selected for coating. About 100 μL mixed liquid was taken on solid de Man, Rogosa and Sharpe (MRS) medium and 50 μL 1280 μg/mL vancomycin solution and 50 μL 1280 ug/mL nor oxacin solution was added evenly. Date, number, weeks and concentration were marked on a petri dish, which was placed in an incubator at 37 °C for 48h.
(3) The single colony cultured in the medium was selected and again cultured in the test tube of 5 mL MRS broth medium. After being shaken uniformly, the bacteria were cultured in an incubator at 37 °C for 48h.
(4) Colony PCR was used to validate the bacterial solution by strain-speci c primers. Samples that could be ampli ed were preserved and whole-genome sequenced for further validation.
Isolate whole genome sequencing and data quality control.
Bacterial genomic DNA was extracted from the isolates for whole-genome sequencing. Using highthroughput sequencing, paired-end reads (2 x 150 bp) library of each single bacteria sample prepared on the Illumina Hiseq 2500 platform in Shanghai Personal Biotechnology company (Shanghai, China).
Quality control was carried out with FastQC, AdapterRemoval (v2.1.7) was used to remove the joint contamination 31 and SOAPec (v2.0) software was used to carry out quality correction based on Kmer frequency 32 .

SNP-Calling and construction of the phylogenetic tree of mutants
To estimate the evolutionary distance between isolates across the 3 models and phylogenetic tree construction, we aligned all short reads to the reference genome of Lactobacillus plantarum HNU082 for SNPs identi cation. Reads were aligned using Bowtie2 33 (Alignment parameters: bowtie2 -p -x --no-mixed --very-sensitive --n-ceil 0,0.01 -1 -2 | SAMtools sort -O bam -@ 24 -o -> *.bam). Candidate SNPs were identi ed and ltered with SAMtools 34 . In particular, genomic positions were considered to be candidate SNP positions if at least one pair of isolates were discordant on the called base and both members of the pair had: FQ scores (produced by Bftools) less than 60, at least 7 reads that aligned to each of the forward strands and reverse strand and a major allele frequency of at least 90%. If the median coverage across samples at a candidate position was less than 10 reads or if 33% or more of the isolates failed to meet the lters described above, this position was discarded 8 . We generated a neighbor-joining tree from the concatenated list of variable positions from conserved genomic regions present in all isolates from all samples by MEGA-X and iTOL software 35 . When computing the distance between each pair of isolates, we only used variable positions that had unambiguous nucleotide calls from both isolates.
Calculating the distance to the most recent common ancestor (dMRCA) and the mutation type To calculate dMRCA for each model at each time point, we counted the number of positions at which the called allele was different from the ancestral allele for each isolate, assessing only SNP positions that were polymorphic among isolates from the particular time point, and averaged the results. SNPs were categorized into 6 types, based on the chemical nature of the single nucleotide changes. we computed the mutation spectrum for each model and then computed the mean and standard deviation of each of the 6 types. The frequency of the mutation type of G-C to A-T was the most abundant in all three models, but the frequency of the mutation type of A-T to G-C was signi cantly higher in the zebra sh model than that in the other two models.
Identi cation of the mobile elements of Lp082 during ingestion To identify mobile elements of the strain Lp082, we performed the pan-genome analysis of all strains (N=109) isolated in this study, including the original reference strain and all isolates from different hosts and time points. In general, the whole genomic sequencing reads of each isolate was assembled by SPAdes with the parameters as below: spades.py -t 24 -k 21,33,55,77,89 --careful --only-assembler 36 . To delineate the core genome of all strains, we mapped the assembled contigs (including only those with the length of >500 bp) of each strain against the reference strain of Lp082 using BLASTn 37 with identity levels ≥ 90% and an e-value < 1e-5. Then, we obtained a set of strain-speci c contigs by excluding the contigs that mapped to the core-genome. Further, to remove the potential redundancy in the strainspeci c sequences, we pairwisely compared these sequences using BLAT 38 with an identity level ≥ 90% and a match length ≥ 85%, which resulted in a set of non-redundant contigs speci c to strains that can be also termed as the accessory genome. We reason that all the contigs in the accessory genome of this strain can be inferred as the "mobile elements" inserted into the reference genome (Lp082) during microbial colonization.
To investigate the presence/absence of each potential mobile element in a given isolate, we calculated the read coverage and depth of each isolate against the accessory genome: coverage >80% and depth >60x were considered to be present, while coverage <20% or depth <60x were regarded as absent 39 .
Accordingly, we calculated the coverage and depth of each mobile element in each isolate (Table S5).
Identi cation of the genes with parallel evolution among the models We counted a gene as under parallel evolution if, in at least one host, the gene had multiple independent SNPs and more than 1 SNP per 2,000 bp (to account for the fact that long genes are more likely to be mutated multiple times by chance). Cases in which two SNPs in the same gene that always occurred together in the same isolates were not included as parallel evolution 8 . based on these principles, a total of 5 genes related to carbohydrate utilization and transposase were identi ed under parallel evolution.
Evolutionary dynamics of the strain Lp082 in the gut While too few probiotic isolates (n=7) in the zebra sh model over the 28-days experiment, we only characterized the evolutionary dynamics of the probiotic strain within the human and mouse model. To increase the reproducibility of SNPs studied, we only focused on SNPs that should be identi ed from at least 4 isolates. For each of the 22 SNPs that met this criterion, we calculated the frequency of reads at each SNP position that agreed with the mutation (derived) allele. To ll the time points where no strain was isolated, we generated a continuous relative abundance of SNPs over time by continuous bezier interpolation 43 . To con rm the parent and child relationship of these SNPs, the phylogenetic trees were constructed for all strains isolated in the human and mouse model. To visualize parent and child lineages separately, we subtracted the relative abundance of a parent sublineage by the sum of relative abundances of its child sublineages. When the combined relative abundance of child sublineages exceeded that of their parent sublineage, we set the frequency of the parent sublineage to 0. At last, the muller plot was employed for dynamic visualization 40 . We did not build the model of the evolutionary dynamics of this probiotic in the zebra sh model, as too few isolates (only seven) were found.
Metagenomic DNA extraction and shotgun metagenomic sequencing and data quality control The QIAamp® DNA Stool Mini Kit (Qiagen, Hilden, Germany) was used for DNA extraction from the fecal samples. The quality of the extracted DNA was assessed by 0.8% agarose gel electrophoresis, and the OD 260/280 was measured by spectrophotometry. All of the DNA samples were subjected to shotgun metagenomic sequencing by using an Illumina HiSeq 2500 instrument in the Novogene Company (Beijing, China). Libraries were prepared with the paired-end reads (2 x 150 bp). The raw reads were trimmed using Sickle(https://github.com/najoshi/sickle) and subsequently aligned to the human genome to remove the host DNA fragments.
Identi cation of metagenomic species, microbial functional genes, and metabolic pathway annotation Bracken 41 was applied for metagenomic species identi cation and abundance estimation (Table S6). For metagenomic functional features and metabolic pathway annotation (Table S7), HUMAnN2 was performed by using the UniRef90 database 42 . More information was listed in ''code availability". Accordingly, we got the relative abundances of intestinal microbial taxonomic compositions, gene families and metabolic pathways, respectively. To perform the differential abundance analysis, we introduced a compositionality-aware statistical method as Morton et al demonstrated 43 . In the microbiota analysis, we performed a widely-used centered log ratio of microbial relative abundance to verify our analysis results. The R package "composition" was applied to implement the centered log-ratio (CLR) transformation for raw relative abundance pro les.
Co-evolution analysis based on shotgun metagenomic data of gut microbiota Based on the correlation of the strain Lp082 with the indigenous microbes, we employed the MIDAS (Metagenomic Intra-Species Diversity Analysis System) to perform intestinal microbiota mutations annotation 44 . Brie y, a reference genome database including 33 species with the abundance more than 0.1% and the species closely related to Lp082 was constructed. Then we mapped the shotgun metagenomic sequencing reads from each microbiome to the reference bacterial genomic database and quanti ed nucleotide variation along the entire genome (i.e. the read depth and observed alleles at each position). The data from multiple samples was used for SNP calling of individual intestinal species. The samples at baseline in each host were set as the reference for the calculation of bacterial mutations occurred within hosts at other time points. The SNPs pro les for these intestinal microbes among the human and mouse models both in probiotic and placebo groups (Table S8).
Experimental veri cation of SNPs and the phenotypic changes associated with SNPs For SNPs veri cation, we retrieved the upstream and downstream 100-bp sequences of each SNP in the nal set and designed primers (Table S3) for PCR ampli cation. The Sanger sequencing results of PCR products were used to verify the nucleotide status of SNPs acquired from the in silico analysis. After veri cation, 22 SNPs out of 71 putative SNPs were nally con rmed.
Further, we examined the potential phenotypic improvements of this strain caused by the selective singlenucleotide mutations in vitro, including the ability of rhamnose utilization, trehalose utilization (Fig. S2C) and acid tolerance. For the polysaccharides utilization experiments, the rhamnose/trehalose was set as the only source of carbohydrate in broth and the bacterial growth conditions including growth rate and the number of colonies were calculated. For the acid tolerance experiment, the mutational strains and the ancient Lp082 strain were inoculated into the broth with pH values of 2, 2.5, 3.0, and 7.0 respectively, and the strain survival rates were calculated.

Statistical analysis
All statistical analyses were performed using R software (version 3.5.1). PCoA analysis was performed in R using the ade4 package. CLR transformation was performed by the "zcomposition" package. The heatmap was constructed using the "pheatmap" package, and the evolutionary dynamics were built using the "ggmuller" package. The differential abundances of various pro les were tested with the Wilcoxon rank-sum test and were considered signi cantly different at p<0.05. For boxplot construction, the package "ggpubr" was used. The edges of the network were calculated by the "SpiecEasi" cooccurrence correlations 24 and visualized in Cytoscape (Version 3.7.1). For more details, please refer to the document "Coevolution project analysis code and gure reproducibility".

Data and code availability
The sequence data reported in this paper have been deposited in the NCBI database (genome resequencing and metagenomic sequencing data: PRJNA594992, PRJNA597371, PRJNA590026, PRJNA597969 and PRJNA588621). All the project analysis code had been deposited in Github:https://github.com/zhjch321123/Host-dependent-co-evolution-of-supplemented-probiotic Note: "Aa" represented the original single nucleotide; "Alt" represented the mutated single nucleotide; "MT" represented Mutation Type, Synonymous mutations or Non-synonymous mutations; "AAC" represented amino acid changes;