The CiPP study recruited 78 children aged 2–11 years with ongoing CD autoimmunity, i.e., defined as screening-identified persistent positivity for tTGA in two consecutive samples, enrolled between March 2012 and August 2015 (15). The enrolled children were identified among carriers of HLA-genotypes associated with CD (DR3-DQ2/DR3-DQ2, DR3-DQ2/DR4-DQ8, DR4-DQ8/DR4-DQ8, and/or DR4-DQ8/DR8-DQ4).

Participating children were invited to a randomisation and baseline visit (visit 0) and scheduled for follow-up visits ~3 (visit 1) and 6 (visit 2) months later. Participants were randomised at a 1:1 ratio to either placebo or intervention group. Among the 78 enrolled children (placebo, n = 38; intervention, n = 40) (characteristics listed in Table 1), 63 (81%) provided faeces samples for all three visits.

The study product was an equal mixture of Lactiplantibacillus plantarum HEAL9, formerly classified as Lactobacillus plantarum, and Lacticaseibacillus paracasei 8700:2, formerly classified as Lactobacillus paracasei, in a total bacterial dose of 1 × 10¹⁰ CFU/sachet in maltodextrin (1.0 g). The placebo product consisted of maltodextrin only. The combination of the two Lactobacillus strains was chosen due to their different physiological effects, i.e., Lactiplantibacillus plantarum HEAL9 targets the permeability of the mucosa, and Lactocaseibacillus paracasei 8700:2 targets the immune system (26, 27). All enrolled children followed a regular gluten-containing diet during the study.

Stool sample collection was carried out at home by the study participant's caregiver. Samples were stored at −20°C until they were transported to the lab, where they were kept at −80°C until the analysis. The faecal microbiome was previously analysed using 16S rDNA sequencing (16). The ensuing sequencing data were reprocessed for this study, using the DADA2 pipeline (28) with the non-redundant Silva database version 138 (29).

The faecal aliquots were prepared accordingly to Jaimes et al. (30). All chemicals and reagents used were of analytical grade and were purchased from Sigma-Aldrich (Merck, Darmstadt, DE). The ¹H NMR spectra were recorded on a Bruker Avance III HD spectrometer equipped with a broadband fluorine observation SmartProbe^TM with z-axis gradients (Bruker BioSpin GmbH, Rheinstetten, Germany) operating at the proton frequency of 500.18 MHz. All samples were acquired using a 1D NOESY pulse sequence with presaturation and calibrated to the internal standard (3-(trimethylsilyl) propionic-2,2,3,3-d₄ acid sodium salt) at 0.0 ppm, manually phased in TopSpin 3.6.4 (Bruker Biospin GmbH, Rheinstetten, Germany). The spectra were pre-processed with an in-house script under MATLAB^® R2020a (MathWorks, Natick, MA, USA) consisting of multipoint baseline correction in user-defined segments, ensuring the same pre-processing for all the spectra. Spectra between δ 0.5 and 9.0 ppm (excluding the residual water region, δ 5.1–4.6 ppm) were reduced into defined buckets; each bin representing a spin system or a part of a spin system that was ideally pure, distinct, and quantitative—in most cases, one bin for each metabolite. Ranges for the bins were chosen after annotation of a subset of spectra in the software Chenomx ver. 8.6, using the built-in spectral library, our in-house database, and published annotated stool spectra (30–33). A detailed description of the workflow is shown in Supplementary Figure 1. For statistical processing, buckets were normalised using probabilistic quotient normalisation (34).

The metabolome was evaluated using the principal components analysis (PCA) algorithm implemented in PLS-Toolbox 8.9 (Eigenvector Research, Wenatchee, WA, USA) under MATLAB^® R2020a environment including buckets of annotated and unknown peaks. PCA was run with concentrations of all annotated and unknown peaks identified in the spectrum, while all other tests were run using only data of annotated features. Linear mixed-effects models were built to characterise the changes in individual metabolite abundance during the three visits; only the annotated buckets were included. The models included individuals as random effect varying intercept only, and the fixed effects were clinic visits and intervention arm. The interaction term between clinic visits and interventions was considered of interest with the aim of showing the impact of dietary intervention in time. The analysis was carried out using the package 'lmerTest' v. 3.1–3 (35) in R v. 4.2.1 (36). Additionally, the Wilcoxon sum rank test and the Wilcoxon signed-rank sum test were applied for pairwise comparison between interventions at each visit as well as for comparison within visit for each intervention, for the subjects who completed all three visits. Correlations between faecal metabolome and faecal microbiome were tested using Spearman's rank correlation test considering ρ > |0.5| using family taxonomic levels that were detected at least in 40% of the samples; the correlations were visualised as a heatmap using Euclidean distance for clustering. Associations between metabolome and 16S rRNA microbiome profiles were analysed on paired metabolome-microbiome datasets of all samples regardless of the group and visit. Only samples from children that provided faecal samples at all clinic visits (visit 0, visit 1, and visit 2) were used for the Wilcoxon tests. After the exclusion, 189 samples (n = 63) were compared (intervention group, n=32; placebo group, n = 31). For all other analyses, all faecal samples were included (intervention group at visit 0, n = 38; intervention group at visit 1, n = 37; intervention group at visit 2, n = 35; placebo group at visit 0, n = 36; placebo group at visit 1, n = 34; placebo at visit 2, n = 36). Furthermore, to uncover physiological patterns, a pathway analysis between the groups at visit 2 was conducted using MetaboAnalyst 5.0 and the KEGG metabolic library using the genus E. coli as a model organism and a proxy for the whole gut microbiota (37, 38). The analysis was performed considering the metabolite overrepresentation within a pathway using the global test, and the influence of changed metabolites on the pathway's function through relative betweenness centrality.

Forty-six metabolites were identified using a semi-targeted approach covering primarily molecules of short-chain fatty acids, branched-chain fatty acids, amino acids, and sugars. The most abundant metabolites were propylene glycol, acetate, butyrate, leucine, and alanine. Additionally, 145 unknown spectral features were uncovered (Supplementary Table 1). An exploratory data approach based on PCA was used as a first step to assess the progress of individuals through the study. Neither common tendency nor clear distinction between the subjects based on the placebo or intervention group was observed at the baseline, visit 1, and visit 2 (Figure 1; Supplementary Figures 2A–C). Therefore, a PCA model showing subjects at visit 1 and visit 2 after having subtracted the registered concentrations at visit 0 on a subject-by-subject basis was considered to highlight the trends possibly hidden by the differences naturally occurring among subjects. No common tendency was observed (Supplementary Figure 2D). On the other hand, the multilevel modelling showed a significant association (p < 0.05) of the interaction between time (visits) and the intervention and placebo with eight faecal metabolites (Figure 2). These included a decrease mean slope in the intervention group of threonine (p = 6.9 × 10⁻⁵, Figure 2A), methionine (p = 0.015, Figure 2B), leucine (p = 0.022, Figure 2C), valine (p = 0.027, Figure 2D), isoleucine (p = 0.043, Figure 2E), phenylalanine (p = 0.046, Figure 2F), and marginally fumarate (p = 0.049, Figure 2G), and a positive slope for the placebo group (Figure 2). Oppositely, 4-hydroxyphenylacetate (p = 0.001, Figure 2H) increased within the intervention group.

In addition to the linear mixed-effects models, univariate statistics was performed with pairwise comparison for each visit separately. This approach has uncovered differences at the baseline between fumarate production, higher in the intervention group (p = 0.021) caused by an outlier in the treatment group. Thus, the two groups were considered equal in terms of metabolomic profile. The stool composition of donors under probiotic supplementation differed in lower ethanol (p = 0.017) and glycerol (p = 0.046) at 3 months when compared to the placebo group; nevertheless, a similar difference did not occur at 6 months visit. At 6 months, probiotic supplementation caused a significant increase of 4-hydroxyphenylacetate (p = 0.019, Figure 3A), while a decrease was reported in aspartate (p = 0.037, Figure 3B), lactate (p = 0.027, Figure 3C), and threonine (p = 0.001, Figure 3D) when compared to the baseline.

The pathway analysis showed depletion in six metabolic pathways; glycine, serine, and threonine metabolism (p = 0.013); cyanoamino acid metabolism (p = 0.024); methane metabolism (p = 0.025), pyruvate metabolism (p = 0.027); cysteine and methionine metabolism (p = 0.033); and nicotinate and nicotinamide metabolism (p = 0.049). The perturbed metabolic pathways in the faecal samples are shown in Supplementary Table 2 and Supplementary Figure 3.

The relative abundance of identified families and the concentration of the annotated metabolites showed that families Pasteurellaceae, Monoglobaceae, Ruminococcaceae, Lachnospiraceae, Carnobacteriaceae, and Aerococcaceae were positively associated with saccharolytic metabolites such as glucose and revealed negative correlation with proteolytic metabolites (Figure 4). Of particular interest was the inverse correlation between the family Pasteurellaceae and glucose (ρ = 0.58, p < 2.2 × 10⁻¹⁶; Figure 5A). The reverse was observed for families Rikenellaceae, Oscillospiraceae, Christensenellaceae, Oscillospirales family UCG-010, Oscillospirales family UCG-011, Marinifilaceae, Barnesiellaceae, Akkermansiaceae, Eubacterium coprostanoligenes group, Anaerovoracaceae, Defluviitaleaceae, Peptococcaceae, Tannerellaceae, and Desulfovibrionaceae, which were positively associated with an environment rich in proteolytic metabolites such as phenylacetate and negatively correlated with saccharolytic metabolites (Figure 4). An inverse correlation was observed between the family Rikenellaceae and glucose (ρ = −0.58, p < 2.2 × 10⁻¹⁶; Figure 5B). The abundance of Rickenellaceae diminished exponentially with increased glucose. On the other hand, the family Rikenellaceae was positively correlated with metabolite isovalerate (ρ = 0.52, p < 2.2 × 10⁻¹⁶; Figure 5C). The family Lactobacillaceae did not manifest any significant correlation.