Probiotic Supplementation Impacts Gut Microbiota Development in Infants at Risk of Atopy: A Follow-Up Analysis of a Randomised, Double-Blind, Placebo-Controlled Study

1. Introduction

The community of microorganisms that reside in the human gastrointestinal tract, known collectively as the human gut microbiota, is increasingly recognised for its critical role in overall health, influencing metabolic, immunological, and neurological functions throughout life [1,2,3]. The establishment of the gut microbiota plays a pivotal role in shaping how the immature immune system develops, with implications for health and quality of life during childhood and beyond [4,5,6]. Recent studies have proposed a “typical” pattern of microbiota progression within the first year of life that supports healthy immune system development [7,8,9,10]. This progression is characterised by the dominance of bifidobacterium species, which utilise human milk oligosaccharides (HMOs) in breastmilk [5,11], help establish a low abundance of pathobionts [8], and achieve a balanced rate of diversification that is neither delayed nor accelerated [9,12]. Practices such as caesarean delivery, infant formula feeding, and antibiotic exposure may disrupt this “typical” pattern [7,13,14] and are associated with higher incidences of immunological disorders, including allergies, atopy and autoimmune diseases [9,12,15]. Given the potential lifelong implications of establishing a healthy gut microbiota during infancy, there is significant interest in exploring the use of interventions to support microbiota development [2,16,17].

Probiotics offer one option as an intervention for the modulation of the developing infant microbiota. Probiotics are defined by the World Health Organisation as ‘live microorganisms which when administered in adequate amounts confer a health benefit on the host’ [18] and several studies report beneficial effects in lowering the incidence of atopy [16,19] and eczema [20], especially when given both prenatally to mothers and postnatally to infants [21]. However, studies have shown mixed results in terms of whether probiotic supplementation can alter the composition of the gut microbiota of infants born after 37 completed weeks of gestation [22].

The ‘Probiotics in the Prevention of Atopy in Infants and Children’ (PROBAT) study (ISRCTN26287422), was a large, randomised, double-blind, placebo-controlled trial which, assessed the ability of the Lab4B probiotic to prevent the development of atopic conditions in a cohort of 454 new-borns [23]. The study cohort predominantly comprised infants with an anticipated increased risk of developing atopy (a first-degree relative with clinically diagnosed asthma, eczema, or allergic rhinitis). The probiotic reduced the incidence of atopic eczema and protected against sensitisation to cow’s milk and egg proteins. Details of study procedures are described elsewhere [24].

In the current study, faecal samples from a subset of the first 300 infants recruited to the PROBAT trial were selected for microbial analysis using traditional plate culture and DNA sequencing techniques. The objective was to identify differences in the composition of the developing gut microbiome of infants in the Lab4B probiotic and placebo groups.

2. Materials and Methods

2.1. Study Design and Sample Numbers

During the PROBAT trial, expectant mothers in the final month of pregnancy were randomly assigned (1:1) to receive a daily supplement of the Lab4B multi-strain probiotic or a matching placebo until delivery, with their babies receiving the same daily intervention until the age of 6-months. Randomisation was performed by external commercial partners using a computer-generated random sequence without blocks [23]. Participants and researchers were blind to study group allocation until data collection was completed and databases were locked [23]. Ethical approval was granted in February 2004 by the South West Wales Research Ethics Committee on behalf of NHS Wales (Project Ref. HE09 COL 1002) and recruitment took place between May 2005 and October 2007 [25]. Although the inclusion criteria required mothers carrying a foetus with a first-degree relative clinically diagnosed with atopy, this criterion was not always met.

The probiotic daily dose of 1 x 10¹⁰ colony-forming units (CFU) comprised 6.25 x 10⁹ CFU Ligilactobacillus salivarius CUL61 (National Collection of Industrial, Food and Marine Bacteria (NCIMB 30211) and 1.25 x 10⁹ CFU each of Lacticaseibacillus paracasei CUL08 (NCIMB 30154), Bifidobacterium animalis subsp. lactis CUL34 (NCIMB 30172) and Bifidobacterium bifidum CUL20 (NCIMB 30153) on a maltodextrin base. The placebo comprised maltodextrin only and was identical in appearance to the probiotic. The intervention was taken orally by the mothers and for the infants, the powder was either sprinkled into the infant’s open mouth or mixed with infant formula/expressed breast milk.

Carers of the infants were contacted by a research assistant at home and/or over the phone at 6, 12, 18, and 24 weeks to complete questionnaires which recorded information on common symptoms, compliance with the intervention (number of unused capsules), feeding method, visits to general practitioners, hospital admissions, medicines administered (including any antimicrobials), any adverse events, and the infant’s general health [24,25].

Carers were requested, but not obliged, to submit nappies containing their infants' faecal samples at birth, 2 weeks, 6 weeks, 12 weeks, 18 weeks, and 6 months - but provision was sporadic. Nappies were collected by the carers, transferred to Anaerogen bags (Sigma-Aldrich, UK), anaerobically sealed and stored in a refrigerator prior to transfer to the laboratory (typically within 3 days) for storage at -80°C pending further analysis. For this follow-up analysis, antibiotic-treated infants were excluded, and faecal samples from which high-quality DNA was extracted were analysed; 105 samples from 46 infants in the placebo group and 113 samples from 54 infants in the probiotic group.

2.2. Faecal Viable Numbers

Faecal samples from infants who had not received any antibiotics were enumerated within three months of their arrival at the laboratory using traditional culture techniques on selective agars to quantify viable Lactobacilli, Bifidobacteria, Enterobacter, Enterococci, Bacteroides, Staphylococci, Streptococci, Clostridia, as well as yeast, total aerobes, and total anaerobes, following the method of Madden et al. [26]. The detection limit was 5 log₁₀ CFU/g faecal dry weight. Any faecal material not used for this testing was stored at -80°C.

2.3. Genomic DNA Extraction and Quantification

In 2022, the faecal samples were removed from storage and genomic DNA was extracted. The samples were mechanically lysed using Matrix Lysing B tubes in conjunction with a FastPrep-24^TM bead beater (MP Biomedicals, USA) and three cycles of 30 seconds at 5 metres/second (with a 5-minute interval between each cycle). Genomic DNA was extracted from the lysed samples using the QIAamp Fast DNA stool kit (Qiagen, UK) and DNA concentrations were measured using a Qubit^® 2.0 Fluorometer (Invitrogen, USA), both in accordance with the manufacturer’s instructions. Prior to sequencing, the quality of extracted genomic DNA was assessed by running 20 ng of total isolated gDNA on a 0.7% agarose gel. Intact, non-sheared gDNA was evidenced by discrete, high molecular weight bands on the gel. Isolated gDNA was diluted to 4 ng/µl in nuclease-free water.

2.4. Genomic Analysis

Samples which yielded high-quality DNA were categorised by the age of the infant at the time of sampling; 0 to 6 weeks (T1), 7 to 12 weeks (T2), 13 to 18 weeks (T3) and 19 to 24 weeks (T4). A subset of these samples was analysed in more depth by metagenomic shotgun sequencing: 16 samples (8 placebo and 8 probiotic) from 16 infants up to 2 weeks old (Starting Point [SP]), and 18 samples (5 placebo and 13 probiotic) from 18 infants aged between 19 and 24 weeks (Ending Point [EP]).

2.5. Analysis of the Faecal Microbiota by 16S rDNA

2.5.1. Sequencing and Initial Processing

16S rDNA sequencing from the genomic DNA extracts was performed in accordance with a previously described method [27]. Briefly, sample libraries were prepared using a modified version of Illumina’s 16S Metagenomic Sequencing Library Preparation Protocol targeting the V1-V2 hypervariable region of the 16S rRNA gene [28]. Sample libraries were quantified using the NEBNext Library Quant Kit for Illumina (New England Biolabs, Hitchin, UK) and sequencing was performed on an Illumina MiSeq platform with the MiSeq Reagent Kit v3 (Illumina Inc., Saffron Walden, UK) using paired-end 2 × 300bp chemistry. Raw sequencing data was processed following the DADA2 [29] pipeline in R. Paired-end reads were filtered and trimmed to ensure all nucleotides had a Q score ≥ 30. Forward and reverse reads were merged and chimeric sequences removed. The resulting amplicon sequence variants (ASVs) were assigned taxonomy via alignment to the SILVA database v138 [30].

2.5.2. Bacterial Taxonomic Analysis

The tax_glom function of the R package Phyloseq [31] was used to collapse amplicon sequence variants (ASVs) to the level of family, genus or species, and counts were converted to relative abundance. The 10 most abundant bacterial families across all samples and all time points were identified by summing the relative abundance of all samples, and the remaining families (77) were combined into the ‘Other’ category. Relative abundance data was centre-log-ratio (CLR) transformed using the Microbiome package [32] before statistical analysis.

2.5.3. Differential Abundance

At each time point, differential genera were identified with MaAsLin2 [33] using relative abundance data and controlling for the random effect of the participant and with all covariates as fixed effects (genera present in 10% of samples were compared with the Compound Poisson Linear Model (CPLM) analysis method, default transformation and no normalisation). False discovery rate (FDR) was controlled with the Benjamini-Hochberg method with a significance threshold of Q ≤ 0.05.

2.5.4. Diversity Measures

Alpha and beta diversities were calculated with Phyloseq. Alpha diversity was quantified from untransformed species counts with Shannon and Simpson diversity indices. Beta diversity was analysed by generating Bray-Curtis dissimilarity matrices from relative abundance data at each time point and visualised with Non-Metric Multidimensional Scaling (NMDS) plots.

2.5.5. Neonatal Community State Type Analysis

Partitioning Around Medoids (PAM) analysis (Cluster R package; [34]) was applied to a Jensen-Shannon divergence (JSD) matrix derived from genus-level relative abundance data for the T1 samples. The optimal number of clusters (Community State Types, CSTs) was determined by maximising the median sample silhouette width (Si) and optimising the ratio of positive to negative Si scores across clusters, with k ranging from 2 to the square root of the number of samples. Cluster medoids were visualised with the JSD matrix on a Principal Coordinates Analysis (PCoA) plot using the "cmdscale" function (k = 2). Bacterial genera (relative abundance) and independent variables (study group and covariates) were assessed for their explanatory power on the ordination using the “envfit” function of the vegan R package [35]. The relationship between independent variables and each CST was assessed using CST membership as a dependent variable and following the approach outlined in the ‘statistical analysis’ section.

2.5.6. Microbial Networks

Data was split per study group at each time point and group sample numbers were balanced (using the ‘sample’ function in R) to avoid sample size disparity bias when comparing networks [36]. The NetCoMi R package [37] with SPRING method [38] was used to generate networks from CLR transformed relative abundance of genera present in at least 10% of samples. Each node (genus) in a network had a score calculated for the following metrics of centrality: degree (number of direct connections to other nodes), betweenness (frequency of appearing on the shortest path between nodes), closeness (proximity to other nodes) and eigenvector (influence on highly connected nodes). For each measure of centrality, the genera which scored in the top quartile of the network were considered ‘keystones’ and were compared between study groups at each time point. Network information was imported into Gephi [39] for visual styling and structured using the OpenOrd force-directed algorithm. Nodes were coloured by bacterial family and their sizes were scaled by ‘hub score’, which is the sum of normalised values for degree, betweenness, closeness and eigenvector centrality scores. Edge size represents the strength of an association between genera; green edges represent positive associations (mutualism), red edges represent negative associations (antagonism).

2.6. Analysis of the Faecal Microbiota by Metagenomics

2.6.1. Sequencing and Initial Processing

Shotgun metagenomic sequencing of the genomic DNA extracts was performed using the Illumina NovaSeq 6000 platform (Novogene, China) with 150-nucleotide long paired-end reads. Adaptor removal and quality control of raw reads was conducted by Novogene. Quality was confirmed using FASTQC [40], ensuring all base calls had a Q score of at least 30. Host sequence decontamination was performed using Bowtie2 v2.5.4 [41] (default parameters) against the GRCh38 genome.

2.6.2. Microbial Profiling and Gene Prediction

The relative abundance of microbial taxa was determined by mapping reads to the CHOCOPhlAn database with MetaPhlAn4 v4.1.1 [42] using default parameters. Reads were assembled into contigs using MEGAHIT [43] and open reading frames (ORFs) were predicted with Prodigal in metagenomic-mode [44]. Redundant ORFs were clustered and consolidated into representative sequences using CD-HIT v4.8.1 [45] with criteria of > 95% identity and > 90% coverage, and used to create a custom database. Metagenomic reads were aligned to the custom database using Bowtie2 with the ‘--very-sensitive-local’ parameters. Mapped and unmapped reads were separated, sorted, and indexed using SAMtools v1.2 [46]. The number of reads mapped to each gene was calculated, and the counts were normalised to reads per kilobase per million mapped reads (RPKM).

2.6.3. Antibiotic Resistance Gene, Mobile Genetic Element and Metabolic Pathway Annotation

Identification of Antibiotic Resistant Genes (ARGs) was performed using the Resistance Gene Identifier (RGI) provided by the Comprehensive Antibiotic Resistance Database (CARD; [47]). The abundance of mobile genetic elements (MGEs) was evaluated with ShortBred v.0.9.4 [48] using a reference database of MGEs [49]. Metabolic pathway profiling was performed using HUMAnN v3.0 with the MetaCyc database [50] as a reference. Differential abundance analysis of pathways at each time point was performed in DESeq2 [51] after filtering for pathways present in at least 10% of samples. P-values were adjusted using the Benjamini-Hochberg correction.

2.7. Statistical Analysis

All statistical analyses of the 16S rDNA data, except for beta diversity and microbial network comparisons, were conducted via Generalised Linear Mixed Models (GLMMs) with the glmmTMB R package [52]. Each model included the participant as a random effect, in addition to study group, time point, and their interaction as fixed effects. Any of the following covariates were included as fixed effects if they significantly reduced the Akaike Information Criterion (AIC): delivery mode (vaginal or caesarean), number of weeks of breastfeeding, infant’s intervention compliance as an interaction with group, sex of infant, Townsend index score, presence of a sibling, Good’s coverage (for 16S rDNA data) and infant’s age in days (if separate models were required for each time point). Model details can be found in Supplementary Tables S3, S5, S9 and S11. To assess the significance of fixed effects, the emmeans package [53] was used to calculate marginal means from model coefficients and perform contrasts with a Tukey’s adjustment to control the false discovery rate.

To analyse beta diversity, we first assessed the community dispersion (within-group variability) of each group at each time point using the betadisper and permutest functions [35]. If significant differences in community dispersion were detected between groups, we balanced group sizes using the sample function in R to minimise the risk of false positives when comparing community compositions using PERMANOVA [35,54]. PERMANOVA tests included covariates as fixed effects if they accounted for a significant amount of variance. Within-group community dispersion between time points was compared using pairwise permutests and community composition was assessed using pairwise PERMANOVA [55]. All PERMANOVAs were run with 999 permutations.

Significant differences in keystone taxa between group microbial networks were assessed using Jaccard’s dissimilarity index (NetCoMi package).

Between- and within-group comparisons of the abundance of ARGs and MGEs were assessed with the two tailed Mann-Whitney U test.

For all statistical analyses, * represents differences between groups at a time point and # the difference within groups at different time points. Significance was considered as */# p ≤ 0.05, **/## p ≤ 0.01, ***/### p ≤ 0.001.

3. Results

3.1. Study Population Characteristics

Table 1 provides a summary of the demographic, environmental, behavioural, and socio-economic characteristics of the population, stratified by group. Rates of health-related events during the 6 month study period for infants in this sub-study are presented in Supplementary Table S1.

3.2. Viable Microbial Numbers

The viable counts for bacteria and yeast enumerated on selective media are presented in Figure 1 (detailed data in Supplementary Table S2). Between-group differences were observed at T2 when the probiotic group had a higher abundance of both Bifidobacterium spp. and total bacterial numbers (7.87 vs 6.34 log₁₀ CFU/g, p = 0.008 and 9.36 vs 9.14 log₁₀ CFU/g, p = 0.040, respectively). At T4, Bacteroides spp. numbers were higher in the probiotic group than the placebo group (1.89 log₁₀ CFU/g and 0.89 log₁₀ CFU/g, p = 0.025). At T4, numbers of yeasts, although very low, were significantly higher in the placebo than probiotic (0.31 vs 0.2 log₁₀ CFU/g, p = 0.014, respectively).

Within-group differences in viable numbers indicated an increase of Bifidobacterium spp. from T2 to T3 in the placebo group (6.34 vs 8.38 log₁₀ CFU/g, p = 0.002) and a decrease of yeasts from T2 to T4 in the probiotic group (0.31 vs 0.2 log₁₀ CFU/g, p = 0.023).

3.3. 16S Analysis of Faecal Microbiota

3.3.1. Relative and Differential Abundance of Bacterial Taxa

The abundances of the 10 most prevalent bacterial families differed between the probiotic and placebo groups (Figure 2A, Supplementary Figure S1 and Supplementary Table S4). The probiotic group was significantly enriched in the Lactobacillaceae family throughout the study and had a higher prevalence of the Bifidobacteriaceae at T1. The placebo group had a higher prevalence of the Clostridiaceae and the Staphylococcaceae at T1 and T3, respectively.

At T1, Bifidobacteriaceae represented >70% of the bacterial population in the probiotic group, compared to 51.2% in the placebo group, which significantly increased in this group to 67.17% by T4 (p = 0.006). T1 abundances of the Staphylococcaceae were 5.96% in the placebo group and 1.24% in the probiotic but were not detectable in either group at the end of the intervention period (p = 0.002 and 0.006, respectively). Both groups had an increase in Lachnospiraceae abundance from T1 to T4 (placebo: 0.09% to 4.98%; p = 0.024, probiotic: 0% to 2.07%; p = 0.005).

Between-group differential abundance analysis of all bacterial genera (Figure 2B), indicated significantly higher proportions of Lacticaseibacillus at T1 and T2 and Ligilactobacillus at T1, T2 and T3 in the probiotic group. Statistical modelling focused on the abundance of the supplemented probiotic organisms (Table 2) confirmed the differential abundance analysis results, showing significant differences between groups for Ligilactobacillus and Lacticaseibacillus at T1, T2 and T3. Furthermore, the probiotic group was enriched in B. bifidum at all time points and B. animalis at T2 and T3.

3.3.2. Diversity Measures

There were no significant differences in alpha diversity between groups (Figure 3A and Figure 3B), although at T2 there was a trend (p = 0.059) of a higher Shannon diversity index for the probiotic cohort (median = 1.95 vs 1.81). For the probiotic group both Shannon and Simpson diversity indices significantly increased over time, whereas in the placebo cohort, Shannon diversity remained unchanged but Simpson diversity increased between T1 and T2.

Visualisation of beta diversity revealed that the probiotic and placebo group centroids occupied distinct positions in ordination space at each time point (Figure 4A and Supplementary Table S6), indicative of differences in microbial community composition. There was lower variability in the community composition (community dispersion) of samples within the probiotic group than the placebo, calculated from the distance of samples to the group’s centroid (p = 0.006 at T1, p ≤ 0.001 at T2 and at T3, p ≤ 0.001; Figure 4A and B and Supplementary Table S7). Within each group, the community composition of T1 was distinct from the other time points (Figure 4C and Supplementary Table S8). Additionally, community dispersion varied within groups between time points (Figure 4B and C and Supplementary Table S7). For both groups, community dispersion was highest at T1 and declined over time. For the placebo group, this reduction was only significant by T4 (average sample distance to centroid at T4 was 0.617 vs 0.659 at T1; p ≤ 0.001). Whereas, for the probiotic group, the community dispersion at T2 was significantly lower than T1 (average sample distance to centroid at T1 was 0.620 but 0.557 at T2; p = 0.003) and remained lower until T4, when community dispersion was comparable to T1.

3.3.3. Neonatal Community State Type Analysis

Analysis of the microbiotas of all the faecal samples at T1 identified three distinct community state types (CSTs), each characterised by predominance of a particular genus (Figure 5A to C). CST 1 was characterised by a high relative abundance of Bifidobacterium (median = 45.4%), CST 2 by Staphylococcus (82.9%), and CST 3 by Escherichia-Shigella (47.1%). The abundance of each of these genera explained a significant (after FDR adjustment) proportion of the variance in the ordination (Figure 5B); Escherichia-Shigella R² = 0.91 (p = 0.046), Staphylococcus R² = 0.74 (p = 0.046), Bifidobacterium R² = 0.65 (p = 0.046). Some variance in microbial community composition was explained by the following independent variables (Figure 5A): infant’s age (R² = 27.4%; p < 0.001), Good’s Coverage scores (R² = 13.7%; p = 0.013), infant’s compliance with the intervention (R² = 10.2%; p = 0.03) and whether the infant was in the probiotic or placebo group (R² = 5.8%; p = 0.021). CST 1 included 44 samples (15 placebo and 29 probiotic; 63.77% of all T1 samples) whereas CST 2 and 3 included 12 (7 placebo and 5 probiotic; 17.4% of T1 samples) and 13 (9 placebo and 4 probiotic; 18.84% of T1 samples), samples respectively.

Associations between independent variables and membership of each CST were evaluated (Supplementary Table S10). CST 1 was strongly associated with older infants (average age: 19 days, range: 1 to 42 days; Adjusted Odds Ratio (AOR): 1.27 [1.06, 1.52] per day increase; p = 0.009) and trended toward association with the probiotic group (AOR: 5.69 [0.95, 33.42]; p = 0.057). CST 2 showed no significant association with independent variables, but CST 3 was associated with: (i) younger infants (average age: 5 days, range: 0 to 16 days; AOR: 0.86 [0.75, 0.98] per one day increase in age; p = 0.028) (ii) higher rates of breastfeeding (AOR: 3.23 [1.27, 8.19] for each 20% increase in weeks breastfed; p = 0.041) (iii) lower compliance with the intervention (AOR: 0.43 [0.22, 0.83] for each 20% increase in compliance; p = 0.041).

3.3.4. Microbial Networks and Keystone Taxa

Bacterial networks for the groups constructed at each time point (Figure 6 and Supplementary Figure S2) indicated keystone taxa differed between groups at T1 (related to betweenness, closeness, and eigenvector centrality), at T2 (eigenvector centrality), at T3 (betweenness centrality) and all measures at T4 (Table 3 and Supplementary Table S12). Descriptive and comparative metrics of the whole and largest connected component of each network are provided in the supplementary material (Supplementary Table S13).

3.4. Metagenomic Analysis

3.4.1. Antibiotic Resistance Gene (ARG) and Mobile Genetic Element (MGE) Abundance

No differences in total ARG abundance were observed between the groups (Figure 7A, Table 4 and Supplementary Table S14), but in the placebo group, the genes associated with resistance to beta-lactams and cephalosporins were more abundant at ending point (EP) (p = 0.035 and 0.007, respectively).

Within the probiotic group, between SP and EP there was a decrease in the abundance of genes associated with resistance to beta-lactams (p = 0.010), elfacmycin (p = 0.003), macrolides (p = 0.013), multidrug (p = 0.008) and phosphonic acid (p = 0.037). In the placebo group, there were significant decreases in the prevalence of genes associated with resistance to disinfecting/antiseptic agents (p = 0.011) and penams (p = 0.045).

Analysis of the mobilome (Figure 7C/D, and Table 5) revealed no significant differences between groups at both time points.

3.4.2. Differentially Abundant Metabolic Pathways

Differential metabolic pathway abundance analysis (Figure 7E and Supplementary Table S15) revealed that, at SP the probiotic group was significantly enriched in the ‘glycogen degradation I’ and ‘dTDP-L-rhamnose biosynthesis’ pathways derived from Bifidobacterium bifidum, the ‘D-galactose degradation I’, ‘stachyose degradation’ and ‘5-aminoimidazole ribonucleotide biosynthesis I’ pathways from Bifidobacterium longum and a ‘L-ascorbate degradation I’ pathway derived from Lacticaseibacillus casei. In contrast, the placebo group had a higher SP abundance of pathways derived from Escherichia coli, namely ‘anaerobic sucrose degradation’, ‘heme b biosynthesis from glycine’, ‘mannosylglycerate biosynthesis I’, and four menaquinol biosynthesis pathways. At EP the probiotic group showed enrichment in ‘pyruvate fermentation to isobutanol’ from Bifidobacterium bifidum and several Ligilactobacillus salivarius pathways: ‘dTDP-L-rhamnose biosynthesis’, ‘stachyose degradation’, ‘UDP-N-acetyl-D-glucosamine biosynthesis I’, and three ‘5-aminoimidazole ribonucleotide biosynthesis’ pathways. The ‘dTDP-L-rhamnose biosynthesis’, ‘stachyose degradation’ and ‘5-aminoimidazole ribonucleotide biosynthesis I’ pathways were enriched at both SP and EP, though their microbial sources differed between time points.

4. Discussion

In this subset of infants from the PROBAT study [23], we observed differences in gut microbiota development between those receiving the Lab4B probiotic and those receiving the placebo.

The microbiota composition and maturation from birth for the infants in the placebo group is comparable with that observed in other full-term infants [7,8,56,57] with Bifidobacteriaceae as the most abundant family (51%) at T1 followed by Enterobacteriaceae (19%) and Streptococcaceae (13%). In the probiotic group, Bifidobacteriaceae constituted 73% of the T1 microbiota followed by Lactobacillaceae (12.5%; largely Lacticaseibacillus and Ligilactobacillus), and Streptococcaceae (7%). For both groups, Staphylococcaceae was highest at T1 and declined over time, and this has been seen in other studies [7,56]. Over time for the infants in both groups there was an increase in the presence of bacterial families known to be associated with a more mature microbiota, such as Lachnospiraceae, Bacteroidaceae and Ruminococcaceae [8,13,14,58].

The probiotic group exhibited an age-related microbial diversification that has been seen in healthy infants [12], with temporal increases in the richness and evenness of species (Shannon index) alongside a reduction in the domination by a few groups (Simpson index). In contrast, the placebo group displayed no sustained increase in alpha diversity, potentially reflecting the high atopy-risk population recruited for this study. Delays in microbiota diversification have been associated with the development of paediatric allergies [12] and it has been proposed that a 'microbiota age' which lags behind the infant's chronological age could serve as a biomarker for predicting a predisposition to food allergies [59]. The infants in the placebo group of the PROBAT trial were found to have a higher incidence of food allergy [23] and this may be related to the observed composition of their microbiota.

Beta diversity analysis showed that the gut microbiota composition of infants at T1 was disparate from later time points, and early differences in composition have been previously reported [7,56]. Probiotic supplementation appeared to support the development of a distinct enterotype with between-group differences in community composition at all time points and lower variance between microbiotas (community dispersion) for the probiotic group at T1, T2 and T3.

We observed that the gut microbiotas of infants at T1 were most effectively partitioned into three distinct community state types (CSTs) each defined by the dominance of a specific genus as reported by Shao et al. [15]. At T1, the infant’s age at the time of sampling was the primary driver of microbial differences, with older infants more likely to have a gut microbiota dominated by Bifidobacteria (CST 1) and younger infants enriched in Escherichia-Shigella (CST 3). It has been found that the neonatal gut microbiota is initially dominated by Enterobacteriaceae, with a transition toward a Bifidobacterium-rich community as the neonatal gut matures and becomes increasingly anaerobic [7,13]. There were indications that the probiotic supported such a transition, with CST 1 showing a near-significant association with the probiotic group. The placebo group appeared to align more closely with CST 3.

The bifidobacteria species present in the Lab4B probiotic (B. animalis and B. bifidum) were more abundant in the infants in the probiotic group. Infants with a high abundance of Bifidobacteria in their gut microbiota have been found to exhibit better health outcomes compared to those with lower levels; the benefits observed include a reduced risk of childhood obesity and metabolic disorders [5,11], enhanced resistance to pathogen colonisation [15], and a more tolerant immune system [60] which may contribute to the lower incidence of atopic diseases and cow’s milk allergy [12]. Bifidobacteria also play an important role in the maturation of the butyrogenic microbiota by cross-feeding acetate to support the establishment of butyrate-producing bacteria such as Faecalibacterium prausnitzii, Eubacterium rectale, and species of Anaerostipes, Eubacterium, and Roseburia [61]. Microbially derived butyrate is known to support gut epithelial integrity and reduce the risk of immune-mediated diseases, including food allergies [62,63].

In the probiotic group, Lacticaseibacillus spp. and Ligilactobacillus spp. (components of the probiotic) were detectable throughout the study. Members of the Lactobacillaceae family are typically enriched in vaginally delivered and breastfed infants compared with caesarean delivered or formula fed infants [5]. These bacteria are believed to contribute to positive health outcomes via production of metabolites which lower gastrointestinal pH and may inhibit enteric pathogens [5,64].

There were no between-group structural differences in the genus level microbial networks but keystone genera differed between the two networks at all time points. Keystone (hub) organisms have a large influence on the structure and function of ecological networks but are not necessarily the most abundant organisms [65]. They exert their influence through (i) a high number of direct connections to other organisms (degree centrality), (ii) frequently appearing on the shortest paths between pairs of organisms (betweenness centrality), (iii) being in close proximity to many other organisms (closeness centrality), or (iv) exerting influence over other highly connected organisms within the network (eigenvector centrality). None of the Lab4B probiotic organisms was identified specifically as a keystone organism.

Analysis of the gut resistome revealed no between-group differences in the overall abundance of antibiotic resistance genes (ARGs) but that there were indications of a age/time related decrease in abundance within both groups, as seen in previous studies [66,67,68]. Genes related to multidrug resistance were the most abundant in both groups, as found by Pärnänen et al. [49]; Casaburi et al. [69], and their abundance within the probiotic group decreased significantly between SP and EP. By EP, probiotic supplementation reduced the abundance of ARGs related to beta-lactams and cephalosporin in particular – antibiotics that are used as first-line treatments for serious infections in infants [70]. The abundance of mobile genetic elements was similar between the two groups throughout, suggesting that probiotic supplementation did not impact horizontal ARG transfer. Transposase genes constituted the largest fraction of mobile genetic elements in the gut metagenomes of our cohort, as has been reported for other infant populations [49].

Differential metabolic pathway analysis revealed that at SP, the gut microbiomes of the probiotic group exhibited an enrichment of metabolic pathways related to carbohydrate metabolism, including galactose, glycogen, and stachyose degradation. The enrichment of these pathways may be a reflection of the extensive repertoire of carbohydrate-metabolising genes possessed by the bifidobacteria [71] whose abundance was higher in the probiotic group at T1. Enhanced galactose degradation pathways may support an increased capacity for the gut microbiotas of the probiotic group to utilise galactose from breastmilk. The ‘stachyose degradation’ pathway, enriched in the probiotic group and attributed to B. longum at SP and L. salivarius at EP, enables the metabolism of stachyose (a prebiotic oligosaccharide) into acetate or propionate, which reduces intestinal pH and thereby can inhibit pathogen colonisation [72]. Metabolism of stachyose by the gut microbiota has been found to relieve inflammation [73] and reduce constipation in murine models [74]. Interestingly, Allen et al. [25] reported a reduced incidence of constipation in the infants in the Lab4B group. Peptidoglycan synthesis pathways were also more abundant in the probiotic group, including ‘dTDP-L-rhamnose biosynthesis’ at both time points and ‘UDP-N-acetyl-D-glucosamine biosynthesis I’ at EP. Peptidoglycan, a major component of bacterial cell walls, plays a critical role in shaping the developing innate immune system and influencing immune homeostasis [75]. Additionally, pathways for 5-aminoimidazole ribonucleotide biosynthesis were enriched at both SP and EP, potentially enhancing the synthesis of ribonucleotides—essential nutrients for infant gut and immune development, primarily derived from breastmilk but also synthesised by the gut microbiota [76,77,78]. The placebo group showed enrichment of metabolic pathways derived from E. coli at SP, including several menaquinol (vitamin K2) biosynthesis pathways. Vitamin K2 is primarily produced by gut bacteria and is vital for bone, cardiovascular, neural, and joint health, although its increased presence in faeces does not necessarily correlate with higher plasma levels [79]. Interestingly, the enrichment of ‘dTDP-L-rhamnose biosynthesis’ and D-galactose degradation I’ pathways, along with depletion of the ‘heme b biosynthesis from glycine’ pathway in the probiotic group aligned with findings from a meta-analysis studying the probiotic supplementation of preterm infants [80].

A strength of this study was that both viable microbial numbers alongside molecular analysis have been used to assess the progression of the intestinal microbiota during early infancy with and without probiotic supplementation. Our analysis was limited to bacteria and did not consider the potential impact of probiotic supplementation on the gut mycobiome and virome. The inconsistency in the frequency of faecal sample provision represented a weakness in this work that was addressed as far as possible through the statistical analysis methodology applied. Additionally, the limited size of the cohort available for faecal analysis prevented any correlations with the primary/secondary trial outcomes. Another observation is that our population subset had higher intervention compliance than the wider cohort, suggesting that carers who provided faecal samples were more likely to follow the recommended dosing protocol.

5. Conclusions

In summary, supplementation with the Lab4B probiotic influenced the development of the infant gut microbiota. The microbiome of the infants in the probiotic group had an increased abundance of Bifidobacteriaceae and Lactobacillaceae, an accelerated diversification, lower levels of some antimicrobial resistance genes, and predicted additional functional capabilities when compared with the control infants. Future longitudinal studies throughout infancy and beyond should be considered to gain a greater understanding of the longer term impacts of probiotic supplementation during neonatal development.