This study analyzed the V3-V4 region of 16S rRNA in stool samples from the prospective, general-population cohort of the All Babies in Southeast Sweden (ABIS) study, which consisted of 17,055 children born between 1 October 1997 and 1 October 1999 in southeast Sweden (Ludvigsson et al., 2001). Biological specimens were longitudinally collected from birth to 13 years of age, creating a biobank of stool, urine, blood, and hair. Metadata was acquired in the form of regular questionnaires, first-year diaries, and human leukocyte antigen (HLA) genotyping. Parental report questionnaires and diaries included information about delivery mode, infection history, antibiotic use, duration of breastfeeding, introduction to and frequency of certain food consumption, living conditions, and additional environmental factors. Stool samples were analyzed from 1,602 children between the ages of 0.3 months and 37.2 months, with an average age of collection of 11.94 months ( Supplementary Table 1 ). This analysis included one stool sample from each child, with a total of 1,602 samples analyzed in a between-subjects design. The average read count for the 1,602 samples was 72,113, the minimum read count was 13,183, and the maximum read count was 623,767. To normalize the data, the read counts for the 1,602 samples analyzed were rarefied to the minimum read count of 13,183. Total abundance as reads per gram was calculated by multiplying the relative abundance and copies of 16S rRNA per gram of stool found from qPCR (Jian et al., 2020).

Upon analysis of the total population cohort, 1,602 children aged 0.3 months to 37.2 months at stool collection, it was determined that Lactobacillus did not colonize the majority of the pediatric population ( Table 1 ). An average total abundance of 224,378 (± 4,007,910; min: 0; max: 156,941,755) Lactobacillus reads per gram of stool were detected, with an average relative abundance of 0.29% (± 2.3%; min: 0%; max: 74.5%) in the total population. Zero Lactobacillus reads were detected in 68.2% of the total population. For comparison, the average total abundance of Bifidobacterium was 139,896,70 (± 95,585,095; min: 0; max: 2.4e+09) reads per gram of stool, with an average relative abundance of 16.6% (± 18.8%; min: 0%; max: 98.9%) and 98.5% prevalence. Lactobacillus ranked 50th in prevalence out of the 205 genera found and populated 31.8% of the entire ABIS cohort studied here ( Figure 1A ).

Within the genus, 25 Lactobacillus amplicon sequence variants (ASVs), each with at least one nucleotide difference, were found within the total population. Every Lactobacillus ASV was detected in a small percentage of the cohort. Thirteen ASVs were in fewer than 1% of the population and only three were in more than 5% of the population, at 9.3%, 9.1%, and 8.2% prevalence, respectively ( Figure 1B ). The top three ASVs in more than 5% of the population were identified using the SILVA_v138 database as Lactobacillus NA 5972, Lactobacillus delbrueckii 5913, and Lactobacillus NA 5971. The 16S sequences from each ASV were individually searched with NCBI BLAST and verified or identified with 100% query coverage and ≥ 99.7% identity. The top three ASVs were identified as Lacticaseibacillus paracasei, Lactobacillus delbrueckii, and Lacticaseibacillus rhamnosus ( Supplementary Table 2 ). The frequency of possessing multiple Lactobacillus ASVs was low, with only 11.8% of the population colonized by more than one Lactobacillus ASV—321 children (20.0%) had one ASV, 96 (5.9%) had two ASVs, 44 (3.7%) had three ASVs, 25 (1.6%) had four ASVs, 15 (0.9%) had five ASVs, six (0.4%) had six ASVs, and three (0.2%) had seven ASVs. No individual was colonized by more than seven Lactobacillus ASVs.

In the total population, Lactobacillus relative abundance and total abundance were negatively correlated with increasing age calculated with Spearman’s rank correlation (RA: rho = −0.096, p = 1.6e-4; TA: rho = −0.10, p = 8.9e-5). In addition, whether a child was actively breastfeeding when their stool sample was provided was the only variable significantly associated with both Lactobacillus relative abundance and/or total abundance when adjusted for false discovery rates using the Benjamini–Hochberg method (Kruskal–Wallis: p_adj = 2.4e-17 and p_adj = 4.2e-18, respectively) ( Figure 2 ). An individual was categorized as actively breastfeeding if their diet included breastmilk, either exclusively or supplementally, at the time of stool sample collection. An individual was categorized as not breastfeeding if their diet did not include breastmilk at the time of stool sample collection.

Further investigation into the defining features of those who were colonized by Lactobacillus, referred to here as the Lactobacillus presence cohort (LPC), was performed. The average total abundance of Lactobacillus in the LPC was 705,820 (± 7,089,251; min: 4; max: 156,941,755) reads per gram of stool, with 0.92% (± 3.9%; min: 0.0075%; max: 74.5%) average relative abundance. As with the total population in the LPC, Lactobacillus relative abundance decreased significantly with age ( Figure 3 ). In the 0 months to 6 months LPC, the mean relative abundance of Lactobacillus was 2.3% (± 10.7%). In the 6 months to 12 months LPC, the mean relative abundance of Lactobacillus was 0.79% (± 2.0%). In the 12 months to 18 months LPC, the mean relative abundance of Lactobacillus was 0.51% (± 1.4%). In the over 18 months LPC, the mean relative abundance of Lactobacillus was 0.02% (± 0.01%). The relative abundance of Bifidobacterium was shown over time for comparison: it significantly declines after 6 months but remains steady. Lactobacillus continued to significantly decrease as age increased. Through a chi-squared test of independence, Lactobacillus presence versus absence was most dependent on active breastfeeding (p_adj = 0.028) ( Supplementary Table 3 ).

The prevalence of Lactobacillus within the total population increased during the period from birth to 6 months, a time period when the majority of the population was actively breastfeeding ( Figure 4 ). Despite the majority of the population continuing to receive breastmilk up to 9 months, the prevalence of Lactobacillus colonization in the total population began to decline after 6 months and dropped starkly after 7 months. The 5 months to 6 months age group had the highest prevalence of Lactobacillus colonization at 61.5%. After 6 months, the prevalence of Lactobacillus colonization dropped from 53.8% at age 6–7 months to 33.3% at age 7–8 months, and 29.1% by 10–11 months. By 17–18 months, Lactobacillus colonization was at 11.1% prevalence.

While actively breastfeeding was a major determinant of Lactobacillus presence, relative abundance, and total abundance, only 55.4% of the total actively breastfeeding portion of the cohort (ABC) was colonized by Lactobacillus ( Table 1 ). Since 45.6% of the ABC remain uncolonized by Lactobacillus, statistical analysis was performed to investigate which variables impact whether an actively breastfeeding child is colonized or not. It was determined that microbial factors did not influence Lactobacillus colonization in the ABC. No other genera relative abundances or total abundances were significantly associated with Lactobacillus presence or absence, when Wilcoxon rank-sum tests were performed and p-values were adjusted for false discovery rate (data not shown). There were also no strong correlations between Lactobacillus relative abundance and other genera with Spearman’s rank correlations. Subsequently, it was determined that additional environmental factors did not affect the colonization of Lactobacillus in the ABC. Chi-squared tests and a Kruskal–Wallis one-way analysis of variance of the environmental factors in parental report questionnaires and diaries showed no significant association with Lactobacillus presence versus absence or relative abundance (data not shown).

Human leukocyte antigen genetics was the only factor tested that showed a significant relationship with whether an actively breastfeeding infant was colonized by Lactobacillus. Lactobacillus presence was determined to be dependent on HLA haplotypes DR5-DQ7 (p = 0.0045), DR15-DQ6.2 (p = 0.048), and DR8-DQ4 (p = 0.078) through a chi-squared test of independence. An odds ratio analysis indicated that individuals with HLA DR15-DQ6.2 were 3.4 times more likely to be colonized by Lactobacillus than those without that haplotype. Individuals with HLA DR5-DQ7 were more likely to not be colonized by Lactobacillus despite actively breastfeeding ( Figure 5 ).

The All Babies in Southeast Sweden (ABIS) study included a prospective general population cohort of 17,055 children born in southeast Sweden between 1 October 1997 and 1 October 1999. Parental consent followed oral, written, and video information about the study. Following parental consent, participating parents completed questionnaires and diaries from birth through the first year of life (Ludvigsson et al., 2001). Collected information included but was not limited to infection history, antibiotic use, duration of breastfeeding, introduction to or frequency of certain food consumption, living conditions, and additional environmental factors. Approval was obtained by the Research Ethics Committees of the Faculty of Health Science at Linköping University, Sweden, Ref. 1997/96,287 and 2003/03-092, and the Medical Faculty of Lund University, Sweden (Dnr 99227, Dnr 99321) as described previously (Russell et al., 2019). The microbial analysis performed at the University of Florida was approved by the University of Florida’s Institutional Review Board as an exempt study IRB201800903.

Pediatric stool samples were collected from 1,756 individual children aged 0.3 months to 37.2 months—the average age of collection was 11.94 months—and processed as previously reported (Russell et al., 2019). Following filtering, this analysis included one stool sample from each child, with a total of 1,602 samples analyzed in a between-subject effect design.

All available HLA genotypes determined from blood spots were included in the analysis. HLA DR-DQ haplotypes of interest here were defined as follows: DR5-DQ7 (DRB1*11-DQA1*05-DQB1*03); DR8-DQ4 (DRB1*08-DQA1*03:03-DQB1*04); and DR15-DQ6.2 (DRB1*15-DQA1*01:02-B1*06:02). Typing was performed with sequence-specific lanthanide-labeled oligonucleotide hybridization (Ilonen et al., 2016).

DNA extraction, 16S rRNA barcoded PCR, and V3-V4 16S rRNA Illumina sequencing were performed on stool samples collected as previously described (Russell et al., 2019). Universal 16S rRNA primers were used to perform bacterial quantification through quantitative polymerase chain reaction (qPCR) as previously described (Russell et al., 2021). Paired-end reads were merged, primers were removed, and ASV processing was performed as previously described (Milletich et al., 2022). Taxonomy was assigned via the assignTaxonomy() function using the SILVA_v138 database within the R package DADA2 (version 1.26) (Pruesse et al., 2007; Quast et al., 2013; Callahan et al., 2016). Species assignments for Lactobacillus ASVs that were below default minimum boostrapping confidence were determined via NCBI BLAST with a 100% query coverage and ≥ 99.7% identity ( Supplementary Table 2 ) (Altschul et al., 1990; Boratyn et al., 2012).

Filtering the original 1,756 children to exclude samples with fewer than 1,000 total reads resulted in the analysis of stool samples from 1,602 children. Filtering to remove ASVs with fewer than five reads using the filter_taxa function from phyloseq (version 1.42.0) resulted in 2,102 ASVs (McMurdie and Holmes, 2013). The average read count for the 1,602 samples analyzed was 72,113, the minimum read count was 13,183, and the maximum read count was 623,767. To normalize the data, the read counts for the 1,602 samples analyzed were rarefied to the minimum read count of 13,183 using rarefy_even_depth() from the phyloseq package. The relative abundance of samples was determined with transform_sample_counts, and the reads/g was calculated by multiplying the relative abundance and copies of 16S rRNA per gram of stool found from qPCR (Jian et al., 2020). ASVs were conglomerated into genera using the tax_glom() function from phyloseq.

All packages were run on RStudio (version 2023.03.0 + 386) (Posit team, 2023). Variables impacting binomial beta diversity were tested for using the permutational multivariate analysis of variance (PERMANOVA) test through the adonis() function in the package vegan (version 2.4-6) (Oksanen et al., 2020). Alpha diversity was calculated with the R function plot_richness() in the phyloseq package. Confounding factors of Lactobacillus abundance and relative abundance were determined via a non-parametric Kruskal–Wallis rank sum test using the R function kruskal.test(). Genera impacted by the presence or absence of Lactobacillus were determined via the non-parametric Wilcoxon rank sum test using the R function wilcox.test(), and p-values were corrected for false discovery rates (FDRs) using the Benjamini–Hochberg method using the R function p.adjust().

The R package PIME (version 0.1.0) determined core microbiomes of each cohort (Roesch et al., 2020). The plot_ordination() function constructed the associated principal coordinate analysis (PCoA) plot.

Relative abundances were calculated using the transform_sample_counts function. To calculate total abundance, the relative abundance values were multiplied by the total number of copies of 16s rRNA per gram of stool, as determined through qPCR. Prevalence was defined as the percentage of individuals in the cohort with non-zero abundance of the respective ASV or taxa. Prevalence was calculated using the getPrevalence() function from R package mia (version 1.4.0) (Ernst et al., 2023).

The odds ratio was calculated with the logistic regression model R function glm() and odds.ratio() function from the questionr package (version 0.7.8). The Spearman’s correlation was calculated with the R function cor(). The correlation plot was created using the R package corrplot (version 0.92) (Wei and Simko, 2023). Odds ratios were calculated using the oddsratio() function from the R package epitools (version 10.1) and MedCalc Software Ltd (Aragon et al., 2020; Schoonjans, n.d). The odds ratio plot was created using the or_plot() function from the finalfit package (version 1.0.6) (Harrison et al., 2023). The heatmaps were designed using the R package pheatmap (version 1.0.12) (Kolde, 2019). Boxplots were designed using the R package ggplot2 (version 3.4.0) (Wickham, 2016). The function stat_compare_means() was used to determine the Kruskal–Wallis and Wilcoxon p-values. The function ggarrange() was used to arrange each graph into a figure.