This study sequenced the genomes of 89 L. reuteri isolates obtained from animal, food, and human intestinal sources, provided by the Lactic Acid Bacteria Collection Center at Inner Mongolia Agricultural University. L. reuteri isolates were activated, cultured, and genomes sequenced. The specific strain information is shown in Supplementary Table S1. The data that support the findings of this study have been deposited into National Genomics Data Center (Ma et al., 2025), Beijing Institute of Genomics (Members C.-N. and Partners, 2025), Chinese Academy of Sciences/China National Center for Bioinformation, under accession number PRJCA042243 that is publicly accessible at https://ngdc.cncb.ac.cn/gwh. Additionally, genome data for 92 L. reuteri were retrieved from the National Center for Biotechnology Information (NCBI) database (as of October 2023). Combining these datasets resulted in a total of 181 genomes (Supplementary Table S2), comprising 55 of animal origin, 85 of food origin, 36 of human intestinal origin, and 5 of unknown origin. More specifically, 65 isolates were from dairy products, 35 from human intestinal tract, 21 from rodents, 17 from fermented foods, 15 from mammals, 11 from ruminants, 8 from birds, and 9 from other sources (including 3 from fruits, 1 from fish, and 4 from unknown sources).

A total of 181 publicly available or newly sequenced L. reuteri genomes were initially retrieved. We used CheckM (version 1.1.2) to assess the integrity and contamination of each genome, with a completeness level greater than 95% and a contamination level below 5%, with all parameters set to default values. Average nucleotide identity was calculated against the L. reuteri type strain, and genomes with ANI < 95% were excluded, consistent with commonly applied prokaryotic species delineation thresholds. After applying all quality and taxonomic filters, 176 high-quality genomes remained and were used for all subsequent comparative analyses.

The isolates, stored in ampoules, were initially inoculated into 5 mL of de Man Rogosa and Sharpe broth (Oxoid, United Kingdom) and cultured at 37 °C for 24 h. After two generations of subculturing, the bacterial cultures were expanded in 50 mL of the same culture medium for genomic DNA extraction.

Bacterial cells were washed twice with phosphate-buffered saline, and genomic DNA was extracted using a bacterial genomic DNA extraction kit (TIANGEN, Beijing, China) following the manufacturer’s instructions.

A sequencing library was constructed for each isolate, and paired-end sequencing (PE150) was performed using the Illumina NextSeq 2000 platform (Novogene, Beijng, China). Raw reads underwent quality control to remove low-quality data and adapter sequences. High-quality clean-data were assembled into contigs using SOAPdenovo (v2.04; Xie et al., 2014), followed by single-base correction. Gaps in the assembled genomes were filled using GapCloser.¹

Average nucleotide identity was calculated using fastANI software (Jain et al., 2018) to assess genetic diversity and species delineation. A fragment length of 1,000 bp was used for comparisons, and ANI values below 95% were considered indicative of different species. The 95% ANI threshold is widely used for bacterial species definition, consistent with the 70% threshold for DNA–DNA hybridization (DDH), and has been widely adopted in taxonomic studies as a conserved threshold for species definition. Applying this threshold ensures that only genuine L. reuteri genes are retained (Olm et al., 2020). The resulting ANI matrices were clustered using the unweighted pair group method with arithmetic mean.

Gene prediction was performed using Prokka software (V1.14.6; Seemann, 2014) and potential probiotic-related genes were identified. Functional annotation of the core and pan-genome sets was conducted using the Clusters of Orthologous Groups (COG) database (Huerta-Cepas et al., 2017).² Differences in COG functional profiles among genetic lineages were analyzed to explore functional genomic variation.

The pan-genome of L. reuteri was analyzed using Roary software (v3.6.1; Page et al., 2015), which constructed core and accessory gene sets at the subspecies level. A minimum BLASTP protein identity of 95% (−i 95) was used to cluster genes into orthologous groups. The core genome was stringently defined as genes present in 100% of the analyzed genomes (−cd 100). Roary was run using the following parameters: -e --mafft -p 4 -r -t 11 -cd 100, the remaining parameters are the default values.

Phylogenetic trees were constructed using the Neighbor-Joining method based on the core gene set sequences with or without using three strains of Limosilactobacillus fermentum as the outgroup (strain ID is EFEL6800, ATCC 14931, DSM20052). The tree was visualized using the Interactive Tree of Life (iTOL) platform,³ enabling comparison of genetic backgrounds and evolutionary relationships among the L. reuteri genomes.

The HMMER package (Minoru et al., 2016) was used to annotate CAZyme in the genome sequences, with reference to the CAZyme database.⁴ Detailed CAZyme family information was retrieved from the CAZy database.⁵

Drug resistance genes were identified using the ResFinder database⁶ with thresholds of > 90% sequence similarity and > 60% coverage. Virulence factors were annotated using the VirulenceFinder database⁷ with the same thresholds (Marchin et al., 2005). Additionally, virulence factors of 176 L. reuteri genomes were detected by the Virulence Factor Database.⁸

Secondary metabolite synthesis gene clusters, including those encoding bacteriocins, were predicted using the Prediction Informatics for Secondary Metabolomes (PRISM) online tool (Krug and Muller, 2014) and the BAGEL4 online platform (van Heel et al., 2018).

The CRISPR-Cas online (Grissa et al., 2007) was used to search for CRISPR arrays and Cas genes in 176 L. reuteri genomes. The number of CRISPR gene clusters, Cas gene types, and spacer sequences was recorded. Variations in CRISPR-Cas systems among these genomes were analyzed to explore their diversity and functional implications.

We performed the Shapiro–Wilk normality test and the Levene homogeneity of variance test. ANOVA was used for comparisons between groups, with a post-hoc Tukey HSD test; non-parametric tests (Kruskal-Wallis test, with a post-hoc Dunn test) were used where necessary. Benjamini-Hochberg correction was used for multiple comparisons. Effect sizes (Hedges’ g) and their 95% confidence intervals are given. All analyses were performed in R (v4.5.0).

The L. reuteri genomes sequenced and assembled in this study exhibited high quality, these results confirm the reliability of the sequencing and assembly data for downstream bioinformatics analysis. The 176 L. reuteri genomes analyzed in this study exhibited an average genome size of 2.08 ± 1.23 Mb, a GC content of 38.75% ± 0.20%, and an average scaffold number of 143.6 ± 129.3. The reference strain L. reuteri DSM 20016^T, isolated from the intestine of an adult human, served as the type strain for comparison (Sriramulu et al., 2008). Its genome size and GC content were consistent with data available in the NCBI database, validating the accuracy of our sequencing and assembly pipeline. Comparative analysis of genomic characteristics across isolates from different regions and sources revealed significant differences in coding sequence counts and genome size (p < 0.05; Figure 1). Foodborne strains exhibited smaller genomes and higher GC content, likely due to gene loss associated with host adaptation (Sheppard et al., 2018). In contrast, animal-derived isolates had significantly larger genomes (p < 0.001), while human intestinal-derived isolates showed intermediate genome sizes but fewer coding sequences. Genome size is a fundamental biological characteristic that exhibits variation among species and even within strains of the same species (Westoby et al., 2021). This variation is influenced by environmental factors and evolutionary processes (Ngugi et al., 2023). Generally speaking, the more complex the environment in which the strain is located, the larger the genome size, and the less the gene decay experienced (Kraaijeveld, 2010). No significant differences in GC content were observed among isolates from different sources.

We then compared the genomic characteristics of isolates obtained from eight ecological sources, including dairy products, fermented foods, rodents, mammals, ruminants, birds, the human intestinal tract, and miscellaneous environments were compared (Figures 1d–f). A substantial degree of variation was observed in genome size, GC content, and protein-coding sequence (CDS) numbers among the various ecological groups. The genome size ranged from approximately 1.8 to 2.2 Mbp. Strains isolated from dairy products and fermented foods exhibited significantly smaller genomes compared to those from rodents, mammals, ruminants, and birds (p < 0.05 to p < 0.001). Strains from host-associated environments, ecially rodents and mammals tended to have larger genomes, whereas human intestinal tract isolates displayed intermediate genome sizes. Similarly, the GC content showed significant ecological stratification, ranging from approximately 38.5 to 40.5%. Isolates from dairy products and fermented foods had the lowest GC content, while rodent and mammalian isolates showed the highest GC content (p < 0.001). Intermediate GC content values were observed in isolates from birds and ruminants. The number of predicted CDSs mirrored genome size trends. Dairy-derived isolates possessed the lowest CDS counts, whereas host-adapted groups rodent and mammal isolates had significantly higher CDS counts (p < 0.01). Only minor variation in CDS numbers was observed among isolates from birds and the human intestinal tract, which formed a middle range between host-adapted and food-origin groups.

These findings suggest that L. reuteri has undergone distinct evolutionary processes to thrive in its respective environments. The larger genome size of animal-derived isolates (p < 0.05) may reflect their need for additional genetic flexibility to adapt to diverse host conditions. Overall, the observed genomic variations among L. reuteri strains are likely influenced by a combination of environmental factors, natural selection, genetic variation, and host-specific lifestyles. Because these source-related genome size and CDS differences were apparent at the whole-genome level, we next examined whether the same ecological partitioning could also be detected at the nucleotide identity level using ANI.

The results of the ANI analysis conducted on the genomes of L. reuteri isolates, including a comparison with the type strain DSM20016^T, are shown in Figure 2. The ANI values across all genomes averaged 96.33% ± 1.59%. Following the species delineation threshold proposed by Richter and Rossello-Mora (2009), an ANI range greater than 95% indicates the same species. Based on this criterion, 176 genomes analyzed in this study were confirmed to belong to L. reuteri. The five genomes with ANI values below 95% were excluded from further analysis.

The ANI results also revealed a clear clustering pattern among genomes from the same source or geographic location. Genomes from the same isolation source or region exhibited high nucleotide consistency, suggesting a close genetic relationship. This clustering trend highlights the influence of isolation source and geographic origin on the intraspecific genetic diversity. Annotation bars show that isolation source explains most clustering, whereas geographic labels are broadly interleaved, indicating weak phylogeographic structure relative to ecology (Pasolli et al., 2019).

The ANI heat map and hierarchical clustering of 176 L. reuteri genomes from different ecological sources showed (Figure 2) that rodents constitute a large cohesive block with uniformly high within-group ANI; narrow dark bands connecting to mammal and human blocks indicate close genomic affinity and host-proximal diversification (Velez et al., 2023). Mammals are distributed in two to three sub-blocks adjacent to rodents with gradual transitions toward rodents and human blocks, compatible with shared ancestry and partial mixing (Li et al., 2023). Ruminants are embedded within a subset block of the mammal and rodents region, exhibiting close similarity and moderate connectivity, suggesting their status as a diet-related sub-structure nested within the mammalian backbone (Yu et al., 2018). Birds constitute a distinct, smaller block with clearer boundaries from the mammal and food block, pointing to a more divergent, host-restricted cohort and limited recent exchange (Zhang et al., 2020). Dairy products form a compact, high-ANI community nested within the mammalian block, suggesting recent differentiation following domestication in stable, nutrient-rich environments (Somerville et al., 2022a). Fermented foods comprise one or two compact communities adjacent to the dairy block; these exhibit high ANI values but slightly increased heterogeneity, consistent with multiple introductions or broader fermentation ecologies (Park and Mannaa, 2025). Human intestinal tract is widely distributed across multiple microbial blocks intersecting with rodent/mammalian blocks, indicating the presence of multiple lineages adapted to humans rather than a single human-specific population (Li et al., 2023). Isolated strains scattered at the periphery of other mixed microbial blocks reflect heterogeneous ecological origins, lacking a single cohesive lineage.

A core gene is a gene that is ubiquitous in all or most strains and is typically linked to fundamental cellular processes such as metabolism, replication, and survival. These genes are highly conserved during species evolution and form the basis of pan-genome analysis.

The core genome of a species comprises conserved homologous genes shared by all strains, which are essential for the biological functions and phenotypic stability of the species (Caputo et al., 2019). Using the Roary result, we analyzed the pan-genome and core-genome of 176 L. reuteri isolates (Supplementary Figure S1). The pan-genome consisted of 16,814 genes, while the core-genome contained 553 genes. Supplementary Figure S1 shows that the number of core genes in L. reuteri tends to remain stable with an increase in the number of strains, while the number of pan-genes increases with the rise in the number of strains, indicating that the genetic material of the species is open (Rajput et al., 2023).

To elucidate the functional roles of the 553 core genes, we performed functional annotation using the COG database. The COG database, maintained by NCBI, classifies gene products based on homology and facilitates the identification of orthologous genes through extensive protein sequence comparisons. Functional annotation revealed 21 COG categories (Figure 3a). Among these, the top four abundant categories were: function unknown [S] (12.58%), general function prediction only [R] (12.14%), replication, recombination and repair [L] (11.23%), and amino acid transport and metabolism [E] (9.7%). Notably, the distribution of these COG categories varied significantly among isolates from different sources (Figures 3b–e). Generally, genomes from animal and human intestinal sources exhibited a richer diversity of functional genes compared to those from food sources, likely due to differences in ecological niches. COG functional categories that showed significant differences among genomes from different isolate sources included information storage and processing, cellular processes and signaling, metabolism, and poorly characterized across isolation sources (p < 0.05; Figures 3f–i). To account for testing several major categories at once, we also ran a multivariate analysis (ANOVA) across these four dominant COG categories, which confirmed a significant overall effect of isolation source (p < 0.001).

To explore functional genomic adaptations of L. reuteri to different ecological environments, genes were assigned to the COG database and compared across eight isolation sources: rodents, mammals, ruminants, birds, dairy products, fermented foods, the human intestinal tract, and other environments (Figures 3j–m). The different COG functional patterns between isolation sources may further reflect their ecological niche (Xiao et al., 2021). Strains isolated from rodents and birds exhibited the highest number of genes in the L category, significantly exceeding those isolated from dairy products and fermented foods (p < 0.001). In contrast, dairy-associated isolates showed the fewest L category genes, suggesting reduced genomic capacity for DNA repair and recombination. In the E category, strains from rodents and mammals contained significantly more genes involved in amino acid metabolism (p < 0.01), whereas dairy and fermented food isolates showed lower gene richness in this category. For R and S categories, which reflect accessory and hypothetical gene diversity, host-associated strains (rodents, mammals, ruminants) had significantly higher gene counts compared to dairy and fermented food isolates (p < 0.001). These patterns suggest reduced functional diversity in food-origin strains compared to host-associated strains. Figures 3n–q indicates that the annotated genes are distributed across multiple functional COG categories, showing clear variation in abundance. Several categories exhibit prominently higher representation than others. Specifically, categories corresponding to “General function prediction only” and “Function unknown” appear to have the largest bar heights, suggesting that a substantial proportion of annotated sequences are either poorly characterized or involved in broad, multifunctional processes. Moderately represented categories include those related to “Amino acid transport and metabolism,” “Carbohydrate transport and metabolism,” and “Energy production and conversion.” These bars are noticeably taller than those representing specialized functions, indicating their key roles in the metabolic landscape of the organism or sample studied. Conversely, categories linked to cellular structure, defense mechanisms, and signal transduction display relatively lower bar heights, suggesting either lower gene representation or less activity within these functional pathways under the examined conditions.

This highlights the evolutionary plasticity of L. reuteri and underscores the importance of ecological context in shaping bacterial functional genomes. These results further underscore the adaptive changes in gene functions that enable L. reuteri to thrive in diverse environments. To determine how these source-enriched gene pools correspond to evolutionary lineages, we constructed a core-gene phylogeny of the same 176 genomes.

Phylogenetic trees were constructed using 553 core genes from the 176 genomes analyzed with and without using Limosilactobacillus fermentum as the outgroup (Supplementary Figure S2). The tree revealed that animal-derived isolates are positioned closer to the root, suggesting they may represent ancestral lineages of L. reuteri. Additionally, genomes from similar isolation sources tended to cluster together, indicating a closer genetic relationship. This clustering suggests that the habitats of these isolates have influenced the formation of distinct phylogenetic branches.

Based on the evolutionary relationships depicted in the phylogenetic tree, the 176 L. reuteri isolates were divided into two major branches and seven distinct clades (Clade I to VII; Supplementary Figure S2). Clade I consists of 10 isolates, with six derived from food sources and four from animals. Clade II also contains 10 isolates, including seven animal-derived, two human intestinal -derived, and one from an unknown origin. Clade III comprises 19 isolates, which include six food-derived, four animal-derived, six human intestinal-derived, and three from unknown sources. Clade IV consists of 13 isolates, with six food-derived, five animal-derived, one human intestinal-derived, and one from an unknown source. Clade V contains 40 isolates: 14 animal-derived, six food-derived, and 20 human-derived. Clade VI includes 20 isolates, made up of 13 animal-derived, six human-derived, and one food-derived. Finally, Clade VII is the largest, featuring 62 isolates, all of which are food-derived, exhibiting a strong clustering trend. The circular phylogeny resolves several well-defined clades with partial enrichment by isolation source.

Rodent isolates are abundant and form one major clade plus several satellites. Some subclades intermix with human intestinal strains, consistent with host-proximal diversification and occasional host shifts (Duar et al., 2017a; Li et al., 2023). Mammal isolates occur in two to three subclades adjacent to rodent lineages; several cluster tightly with rodent genomes, compatible with shared ancestry or recent cross-host movement (Zhang et al., 2024). Ruminants form small, cohesive sublineages embedded within the broader mammal/rodent backbone, indicating close evolutionary affinity with mammalian gut lineages (Li et al., 2023). Bird-derived genomes form a small, discrete clade with relatively long internal branches, suggesting greater divergence and/or undersampling and limited recent exchange with mammalian or food lineages (Greppi et al., 2020). Dairy products strains concentrate in a compact, shallow subclade nested within mammalian lineages, consistent with recent diversification following domestication (Somerville et al., 2024). Fermented foods are two main clusters are observed: one adjacent to the dairy clade and another embedded among mammalian lineages. Diversity is slightly higher than dairy but still restricted relative to host-associated groups (Somerville et al., 2022b). Human intestinal tract isolates are widely distributed across the tree, including a source-enriched clade and several mixed subclades with rodent and mammal genomes, indicating multiple human-adapted lineages rather than a single human-specific cluster. Other isolates are sparsely scattered at the edges of mixed clades, reflecting heterogeneous origins without a single ecological lineage.

The phylogenetic tree also included three L. fermentum strains as an outgroup, which showed no clustering with L. reuteri genomes, confirming their distant genetic relationship (Figure 4). The type strain DSM20016^T clustered with food-derived isolates, while most food-derived strains formed a distinct branch, with a few distributed across other clades. Isolates from similar geographic locations, particularly those derived from animal- and human intestinal-derived sources, exhibited a noticeable aggregation trend, suggesting a link between phylogeny and geographic origin. Notably, sourdough-derived isolates, despite their diverse phylogenetic origins, clustered together, a pattern also observed among food-derived strains in this study. This finding aligns with a previous study reporting significant differences in genome sizes among L. reuteri lineages (Zheng J. et al., 2015b). The results demonstrate that genomes from the same isolation source tend to cluster together, likely due to the influence of geographic location and niche-specific adaptations. This source-enriched phylogenetic pattern mirrors the ANI blocks described, confirming that ecological structuring is detectable at both sequence-similarity and tree-based evolutionary scales.

The clustering trend of food isolates in the phylogenetic tree may indicate that these strains have undergone convergent evolution or adaptive evolution in a specific environment, emphasizing the complex interaction between phylogeny, source, and environmental adaptation.

Carbohydrate utilization is an important factor in determining the colonization and abundance of L. reuteri in the host. Previous sections showed clear ecological structuring of genomes and lineages, we next asked whether CAZymes repertoires follow the same source-related pattern. A total of 72 CAZyme families were predicted across the 176 L. reuteri genomes, including 13 carbohydrate-binding modules (CBMs), eight carbohydrate esterases (CEs), 31 glycoside hydrolases (GHs), and 17 glycosyl transferases (GTs). Among these, GHs and CBMs were the most abundant, followed by GTs. The GH family plays a key role in carbohydrate utilization, enabling L. reuteri to degrade complex carbohydrates. All 176 L. reuteri isolates encode GT2 and GT4, which are mainly synthase enzymes involved in the synthesis of disaccharides such as sucrose, lipopolysaccharides, and cellulose. This suggests that despite differences in isolation sources, the isolates share conserved CAZyme genes, possibly due to common ancestry or horizontal gene transfer.

Notably, the GH73, GH25, and GH13 families were particularly abundant. These families are involved in the degradation of indigestible carbohydrates found in staple foods and dairy products, and GH13 enzymes are involved in starch hydrolysis to hydrolyze digestible polysaccharides (Wei et al., 2023). The CBM family enhances the affinity of GHs for their substrates through a tandem arrangement with GH catalytic domains, facilitating targeted binding to polysaccharides and increasing enzymatic efficiency. The high prevalence of CBM50, which recognizes explicitly lactose, galactose, and their derivatives (Ficko-Blean and Boraston, 2012), underscores the capability of L. reuteri to effectively utilize host-derived carbohydrates. This includes breast milk oligosaccharides and mucin, which are typically challenging for the human body to digest.

Further analysis revealed significant differences in CAZyme profiles among genomes from different isolation sources (p < 0.05; Figure 5). One-way ANOVA indicated that GH13, GH43, GH3, GH36, GT2, GT4, CE10, CE1, and CE2 were all significantly affected by isolation source (p < 0.05). Tukey’s HSD multiple-comparison test, as declared in the figure caption, showed that rodent, mammal, and human-intestinal strains formed the high-abundance group, whereas dairy and fermented-food isolates formed the low-abundance group (Tukey-adjusted p = 0.000–0.047). Because multiple CAZy families were examined simultaneously, we report the adjusted p-value range rather than a single value. Several CAZyme families remained uniquely or predominantly enriched in animal-derived isolates, further supporting host-driven metabolic adaptation. Key differences included variations in GT2, GT4, GT8, GT19, GT28, GH25, GH13, GH68, GH36, GH8, CBM50, CE10, CE1, CE2, and CE7, among others. Additionally, several CAZyme families were uniquely present in specific isolation sources. For example, GT39, GT44, GT11, GT45, GT7, GH30, GH5, GH28, GH115, GH39, and GH67 were exclusively present among animal-derived isolates. GT27 was exclusively identified in food-derived isolates. This enzyme plays a crucial role in the initial step of O-glycan biosynthesis by catalyzing the transfer of N-acetyl-D-galactosamine to serine or threonine residues on protein receptors. Moreover, GT27 demonstrates activity toward the EA2 peptide substrate, although its activity is relatively weak against Muc2 or Muc1b substrates (Wojczyk et al., 2003). AA4, a member of the vanillyl-alcohol oxidase family, is exclusively found in human-derived isolates. It catalyzes the conversion of phenolic compounds that possess side chains at the para-position of the aromatic ring (Gygli et al., 2018).

Comparative CAZy functional classification of L. reuteri genomes from eight different isolation sources (Rodents, Mammals, Ruminants, Birds, Dairy products, Fermented foods, Human intestinal tract, and Others) revealed both a conserved core of carbohydrate-active enzymes and marked source-specific functional differences (Figures 5A–O). Across all groups, the dominant CAZy categories included Glycoside Hydrolases (GH) and Glycosyl Transferases (GT), followed by Carbohydrate Esterases (CE), whereas Polysaccharide Lyases (PL), Auxiliary Activities (AA), and Carbohydrate-Binding Modules (CBM) appeared at lower frequencies. GH families consistently represented the largest fraction of annotated CAZy enzymes, highlighting the central role of carbohydrate degradation in L. reuteri metabolism.

Rodent-derived strains possessed the most extensive CAZyme repertoire, particularly in GH and GT families. GH13, GH43, GH3, and GH36 were predominant, indicating a strong capacity for degrading diverse plant-derived polysaccharides and host mucins. Several GH families (e.g., GH31, GH32) and CE families were uniquely enriched, suggesting metabolic versatility in the complex gut environment of rodents. Human intestinal and fermented food isolates also showed a rich GH and GT repertoire, reflecting a strong adaptation to carbohydrate-rich environments. Fermented food strains shared similar CAZyme reductions with dairy isolates but maintained slightly higher counts of GT2, GT4, and GH73, indicating preservation of cell-wall modification and biofilm-forming capabilities. Their limited GH diversity suggests long-term adaptation to stable, nutrient-sufficient fermentation environments. Human intestinal isolates displayed intermediate CAZyme abundance, dominated by GH13, GH3, and GT2 families, reflecting a generalist carbohydrate-utilization strategy. The presence of several host-glycan-degrading enzymes (GH29, GH95) suggests adaptation to mucin utilization in the human gut. Dairy isolates exhibited marked genome contraction in CAZyme content, particularly within GH and CE families. GH13 and GH25 were the most conserved families, primarily associated with starch and peptidoglycan hydrolysis, respectively. The absence or reduction of hemicellulase-related GH43 and GH51families supports a pattern of functional specialization toward milk carbohydrate metabolism.

By contrast, bird-derived isolates and the “Others” group displayed comparatively lower representation across most CAZy categories. Isolates from ruminants displayed intermediate enzyme profiles, with moderately elevated GT and CE abundances, whereas dairy product isolates showed a balanced GH/GT distribution but low PL and AA content. CBM modules remained uniformly scarce across all groups.

These findings highlight the diversity and host-specificity of CAZymes in L. reuteri and emphasize that CAZy functional composition in L. reuteri is shaped by the ecological niche and host environment, reflecting a pattern of host-driven evolutionary adaptation. This plasticity in carbohydrate-active enzyme content likely contributes to its successful colonization of diverse vertebrate hosts and food matrices.

The L. reuteri genomes from different sources were annotated for probiotic-related genes (Figure 6e). Annotation of functional genes revealed widespread presence of ribF, ribX, ribU, and ribZ (riboflavin biosynthesis), luxS (bioactive peptide synthesis), atpC and atpH (adhesion), pyK and ldh (lactic acid production), rfbX (immunomodulation), and pduC (reuterin production). While all fractions contained these genes, their relative abundance varied, with some sources enriched in adhesion and metabolism-related functions. These findings align with known probiotic mechanisms in L. reuteri (Spinler et al., 2008; Frese et al., 2024; Yu et al., 2018). Rodents and Mammals host-associated groups showed enrichment for luxS and pduC, together with detectable riboflavin genes (ribU, ribZ) in subsets. Ruminants profiles emphasized riboflavin-pathway genes with moderate luxS and pduC, consistent with versatile nutrient metabolism. Birds isolates carried narrower panels, with comparatively fewer luxS and pduC positives. Dairy products genomes were consistently high for glycolysis markers (pyk, ldh) while showing reduced luxS/pduC and patchy riboflavin genes. Fermented foods similar to dairy but with slightly broader riboflavin gene presence in a subset, consistent with multiple introductions into food chains.

Among human, animal, and food-derived, there were no significant differences (p > 0.05) among riboflavin-generating genes (ribF, ribZ, ribU), bioactive peptide genes (LuxS), adhesion genes (atpC, atpH), and lactic acid-producing genes (pyk). However, there were significant differences in the lactic acid-producing gene ldh, the immunomodulatory gene rfbX, and the encodes reuterin gene pduC (p < 0.05). One-way ANOVA followed by Tukey’s HSD showed that human- and animal-derived (host-associated) isolates had higher detection frequencies of pduC and rfbX than dairy and fermented-food isolates (Tukey-adjusted p = 0.004–0.039), while ldh remained broadly conserved but reached higher counts in dairy genomes because of glycolysis specialization. The ldh gene encodes lactate dehydrogenase, which catalyzes the conversion of pyruvic acid to lactic acid and is the core gene for acid production in lactic acid bacteria (Nikkhah et al., 2024). Isolates may regulate acid production capacity through differences in ldh expression levels due to different pH requirements of the environment they are in, such as the human intestine or fermented foods, thereby inhibiting competing microorganisms (Hul et al., 2024). rfbX is involved in the synthesis of O antigen in lipopolysaccharide (LPS) of bacterial cell walls and affects the diversity of surface antigens. The differences in the recognition of O antigen by the immune systems of different hosts (humans, animals, and foods) may lead to other natural selection pressures of the rfbX gene (Klee et al., 2020).

The pduC gene, responsible for the production of the antimicrobial compound reuterin, showed significant variation among isolation sources. The pduC gene is significantly more prevalent in human- and animal-derived isolates compared to those from food sources. The reason for this might be because the pduC gene may be involved in specific metabolic pathways (such as 1, 2-propanediol utilization), which play a key role in the colonization or competition of bacteria in the intestinal environment of human or animal hosts. For example, specific metabolic substrates may exist in the host’s intestine, promoting the selective retention of strains carrying the pduC gene. However, food source isolates may not retain the gene due to the lack of relevant substrates in the environment or different metabolic requirements (Li et al., 2024). This aligns with the ability of many human intestinal and animal-derived L. reuteri to produce reuterin, a well-known antibacterial compound (Schaefer et al., 2010). It exhibits broad-spectrum antimicrobial activity, particularly against L. reuteri.

Comparative genomic analysis was used to assess the safety of L. reuteri at the genomic level by identifying genes related to drug resistance, pathogenicity, and environmental adaptation. This approach provides valuable insights into the biosafety of L. reuteri as a probiotic (Plavec and Berlec, 2020). Strain-specific safety assessments are necessary prior to clinical or food applications.

Analysis of 176 L. reuteri strains revealed the presence of drug resistance genes in different derived isolates. A total of 17 species were annotated in this study, with animal-derived L. reuteri being the most abundant and food-derived being the least abundant (Supplementary Table S3). At the source level (ARG-positive genomes only), Mammal contributed the most hits, followed by Ruminants, birds, human intestinal tract, dairy products, rodents, and fermented foods. Functionally, tetracycline resistance genes tet(W/N/W), tet(W), tet(O/W), tet(O/W/O) dominated the dataset across sources, primarily erm(B) and lincosamide genes (lnuA) also present. A smaller number of chloramphenicol acetyltransferases and streptogramin A acetyltransferases (vatE) were detected.

Notably, the tet and vatE were identified in human intestinal-derived isolates SD2112 and SD-LREI-IT. ErmB was mainly detected in the food-derived HC5-4. Additionally, the lincomycin resistance gene lnuA was found in genomes from all isolation sources (animal, human, and food). Lincomycin, a lincosamide antibiotic, is widely used in livestock to treat respiratory and gastrointestinal infections. Its extensive use has contributed to the emergence and spread of resistance genes, with Campylobacter being a major global cause of bacterial foodborne diseases. An interesting finding is the presence of the fexA gene, which confers resistance to florfenicol, in the animal-derived isolate AN417. This gene has been previously identified in chicken feces and Campylobacter jejuni isolates from chickens. It was first reported on Campylobacter plasmids, providing evidence of horizontal gene transfer between Gram-positive and Gram-negative bacteria (Trieu-Cuot et al., 1985). Further studies indicate that fexA is located on a resistance gene island tet(L)-fexA-catA-tet(O) and can be horizontally transferred, enhancing the adaptability of Campylobacter jejuni under florfenicol selection pressure (Zheng B. et al., 2015a).

Bacteriocins are antimicrobial peptides produced by bacteria that inhibit the growth of other microorganisms (Breukink et al., 1999). Unlike conventional antibiotics, bacteriocins are proteinaceous and typically exhibit narrow-spectrum activity, targeting strains of the same or closely related species. They are also associated with minimal resistance development, as spontaneous mutations conferring resistance are rare (Soltani et al., 2021). When resistance does occur, it is often due to modifications in the cell envelope, such as alterations in charge or thickness (Sakayori et al., 2003; McBride and Sonenshein, 2011).

Bacteriocin prediction analysis across L. reuteri isolates from different sources revealed that animal-derived genomes had the highest number of predicted bacteriocin gene clusters, while food-derived strains had the fewest (p < 0.001; Figure 6a). This may reflect a stronger antimicrobial and probiotic activity of L. reuteri isolates from animal sources (Hernandez-Gonzalez et al., 2021).

Host-associated groups (human, rodents, other mammals) collectively carry the majority of predicted bacteriocin loci, whereas food-origin genomes (dairy, fermented foods) harbor fewer, with dairy showing the lowest count (Figure 6c). Two main gene clusters encoding bacteriocins were identified in L. reuteri genomes (Supplementary Table S4), including enterolysin A: a class III bacteriocin originally purified from Enterococcus fecalis LMG Some studies have also shown that the lack of sufficient 2,333. It exhibits bacteriolytic activity against selected strains of Enterococcus, Micrococcus, Lactococcus, and Lactobacillus, with a dose-dependent killing effect (Nilsen et al., 2003). Sactipeptides: a class of ribosomally synthesized and post-translationally modified peptides characterized by sulfur-to-alpha carbon thioether cross-links. Sactipeptides display a broad spectrum of biological activities, including antibacterial, spermicidal, and hemolytic properties (Chen et al., 2021).

Enterolysin_A like loci are broadly distributed but are most frequent in human, rodent, and mammal isolates, mirroring the overall source ranking. Sactipeptides occur as minor, scattered hits across mammalian sources; they are rare or absent in dairy. Across sources, genomes carried 1–2 BGCs per strain with modest between-group dispersion. Human intestinal tract and fermented-food isolates showed slightly higher average BGC counts, whereas dairy and bird isolates tended to be lower. Rodents and mammal groups were intermediate with wider variance, consistent with mixed ecological pressures on specialized-metabolite capacity. These findings suggest that L. reuteri has the potential to produce bacteriocins, which could contribute to its probiotic and antimicrobial properties. The higher bacteriocin-cluster loads in host-associated sources are in line with their richer probiotic-gene repertoires and with the greater ARG diversity observed in some animal-derived isolates, pointing to stronger competitive and defensive pressures in animal and human intestinal niches.

The CRISPR-Cas system is an adaptive immune system in prokaryotes that protects against the invasion of exogenous genetic elements, such as phages and plasmids (Xu and Li, 2020). This system consists of CRISPR arrays, which store fragments of foreign DNA, and Cas (CRISPR-associated) genes, which encode proteins responsible for targeting and cleaving invading DNA. Together, they form a highly conserved defense mechanism in bacteria.

In this study, the CRISPR-Cas systems of 176 L. reuteri isolates were annotated using online tools (Supplementary Table S5). While differences were observed in CRISPR arrays, Cas gene cluster types, and spacer numbers among genomes from animal, human intestinal tract, food isolation sources, these differences were not statistically significant (Figure 6b). However, after carefully classifying the sources of separation, counts of CRISPR arrays, spacers, and Cas-type diversity varied by the specific source (Figure 6d). Human intestinal tract isolates displayed the highest Cas-type richness and spacer loads, with significant pairwise differences from dairy and fermented-food cohorts. Rodents and mammal genomes showed intermediate CRISPR diversity, while food-origin (dairy food and fermented food) and bird groups were generally lower and more variable. CRISPR array length and spacer number varied across sources, indicating differences in adaptive immune capacity. Fractions with longer CRISPR arrays and higher spacer counts may reflect historical exposure to a broader spectrum of mobile genetic elements such as bacteriophages and plasmids (McCarty et al., 2020). Of the 176 strains, only 41 (23.3%) contained CRISPR-Cas systems, with the following distribution: 66% from animals, 20% from food sources, 12% from humans, and 2% from unknown origins. This suggests that animal-derived L. reuteri isolates are more likely to possess this defense mechanism, potentially reflecting stronger selective pressure from phages in animal hosts. The low prevalence of CRISPR-Cas systems in L. reuteri aligns with previous reports. While CRISPR-Cas systems are present in 85.2% of archaea and 42.3% of bacteria (Liu et al., 2023), only 17% of L. reuteri strains have been reported to contain these systems (Plagens et al., 2015). The relatively low occurrence of CRISPR-Cas systems in L. reuteri is consistent with its niche-specific evolutionary adaptations.

Notably, the L. reuteri genome contains a rich diversity of Type II-A CRISPR-Cas systems, which are known for their strong anti-phage activity. This further supports the hypothesis that animal-derived L. reuteri strains have evolved robust defense mechanisms against phage predation (Ksiezarek et al., 2022). The results of this study are consistent with findings in other lactic acid bacteria. For example, Crawley et al. reported that some bacterial species contain multiple CRISPR-Cas systems, while others, such as Lactobacillus acidophilus, encode CRISPR arrays but lack associated cas genes (Crawley et al., 2018). This parallels the observations in L. reuteri, where CRISPR-Cas systems are present in only a subset of isolates, highlighting the variability in CRISPR-Cas distribution across bacterial species. Together with the source-enriched bacteriocin clusters and the uneven ARG profiles, the CRISPR-Cas patterns support a unified view that ecological origin not only shapes core genome architecture but also modulates accessory, defense, and probiotic traits in L. reuteri.