The study was conducted in accordance with the Declaration of Helsinki, and approved by the Ethics Committee of the University Hospital Ghent (reference number B670201836585).

The effects of EpiCor postbiotic on the gut fermentation characteristics, and gut microbiome and metabolome were investigated making use of the SHIME^® model (Van de Wiele et al., 2015), a dynamic semi-continuous reactor-based model that mimics the different stages of the human gastrointestinal tract (ProDigest-Ghent University, Ghent, Belgium). EpiCor was provided by Cargill (USA) and tested at an in vitro dose of 500 mg per day, an equivalent dose in an adult human.

The reactor configuration for this experiment was adapted from the setup described by Molly et al. (1993), and consisted of a succession of three reactors simulating the different regions of the gastrointestinal tract, i.e., upper gastrointestinal tract including subsequent simulation of stomach and small intestine, proximal colon (PC; pH = 5.75–5.95; 500 mL) and distal colon (DC; pH = 6.6–6.9; 800 mL), respectively. At the start of the experiment, the colonic reactor vessels were inoculated with the fecal material from one healthy human donor, who had not received any antibiotic treatment in the 4 months prior to donation. Inoculum preparation, feeding regime, temperature settings and reactor feed composition were previously described by Possemiers et al. (2004) and adopted during the current study. To maintain anoxic conditions, all reactors of the SHIME^® setup were flushed once every day for 5 min with N₂-gas.

The timeline of the SHIME^® -experiment comprised a 2-week stabilization period, followed by a 2-week control period and a 3-week treatment period (Supplementary Figure 1). During the stabilization phase, the gut microbiota was allowed to differentiate in the different reactors depending on the local environmental conditions, while during the control period the baseline microbial activity, and during the 3-week treatment period, the effect of repeated daily EpiCor administration on the microbiome was investigated. To assess fermentation characteristics, microbiome composition and function, and metabolomic profile, samples were collected from each colonic reactor vessel (PC and DC) three times per week during the 2-week control period (n = 3 samples/week × 2 weeks = 6/reactor) and 3-week treatment period (n = 3 samples/week × 3 weeks = 9/reactor).

Determination of SCFA concentrations, including acetate, propionate, butyrate and branched chain fatty acids (bCFA; sum of isobutyrate, isovalerate, and isocaproate) was executed as previously described by De Weirdt et al. (2010). Lactate levels were determined using a commercially available enzymatic assay kit (R-Biopharm, Darmstadt, Germany) according to manufacturer’s instructions and measured on an iCubio iMagic-M9 (Shenzhen iCubio Biomedical Technology Co. Ltd., Shenzhen, China). Ammonium analysis was run on an AQ300 Discrete Analyzer (Seal-Analytical, Rijen, Netherlands) using the indophenol blue spectrophotometric method according to manufacturer’s instructions.

DNA was isolated as previously described by Van den Abbeele et al. (2018b), starting from pelleted cells originating from 1 mL SHIME^® samples, with minor modifications. A bead-beating step was included using the NGS PréMax Homogenizer device (BioSPX, Abcoude, Netherlands), thereby performing the homogenization procedure twice for 40 s at 6 m/s with the sample being allowed to rest for 5 min between shakings. Following extraction, DNA concentrations were determined using the Quant-iT kit (Thermo Fisher Scientific, USA) following manufacturer’s instructions. Next, libraries were constructed with the SQK-RPB004 Rapid PCR Barcoding kit (Oxford Nanopore Technology (ONT), Oxford, UK). Blanks were included during library preparation for quality control. Shotgun metagenomic sequencing was performed using R9.4.1 FLO-MIN 106 flow cells on the GridION platform (ONT, Oxford, UK), multiplexing 12 samples in each flow cell. Each sequencing run lasted 70 h. MinKNOW ONT software (v 4.2.5) with Guppy basecaller (v 4.3.4) was used for sequencing using the high-accuracy basecalling setting, followed by de-multiplexing, adapter trimming, and quality control using default settings. The sequence data can be found at: https://www.ncbi.nlm.nih.gov/, Submission ID: SUB13991422 under BioProject ID: PRJNA1045629.

Fastq files obtained from the MinKNOW ONT workflow were used for microbial taxonomic classification. First, host DNA was removed by mapping fastq files to the human genome using Minimap2 (Li, 2018) followed by the removal of any reads matching the host genome using Samtools (Li et al., 2009). Remaining reads were assumed to be from microbial DNA and used for taxonomic assignation, which was performed using Kraken2 (Wood et al., 2019) with a custom database containing high quality genomes from RefSeq database (Tatusova et al., 2014), and published metagenome-assembled genomes (Almeida and De Martinis, 2021; Hiseni et al., 2021), classified according to the Genome Database Taxonomy (Chaumeil et al., 2019). Estimation of abundance at different taxonomic levels (i.e., species, genus, family, and phylum) was done with Bracken (Bayesian Reestimation of Abundance with KrakEN) (Lu et al., 2017). First, a Bracken file was created from our customized Kraken database using a read length flag of 500. Bracken was then run to estimate abundance at multiple taxonomic levels from the Kraken report, using default settings.

To assess microbiome functional potential, we identified gene orthologs annotated by KEGG (Kanehisa et al., 2004), and carbohydrate-active enzymes (CAZy) (Drula et al., 2022). First, genomes from microbial species identified with Kraken2 with relative abundance higher than 10^–4 were annotated using PROKKA (Seemann, 2014), followed by additional assessment of gene function using EggNOG-mapper v2 (Cantalapiedra et al., 2021). After the annotation process was completed, KEGG Orthology genes (KOs) and CAZys were compiled for each genome, weighted based on the presence of each species abundance within the community, generating the accumulated KEGG KO and CAZy potential for all the microbes identified in each sample.

First, laser-assisted rapid evaporative ionization mass spectrometry (LA-REIMS) metabolic fingerprinting was performed as a first-line rapid discriminative metabotyping to inform decision making for which samples to select for additional in-depth LC-MS-based metabolomics. Briefly, 200 μL sample was transferred into a well of a 96-well plate and subjected to LA-REIMS analysis as described by Van Meulebroek et al. (2020).

Next, in-depth LC-MS-based metabolomics, including ultra-high pressure liquid chromatography high-resolution mass spectrometry (UHPLC-HRMS) analysis of samples for polar metabolomics and lipidomics was performed on a subset of the collected samples, including three samples collected during the final control week and five samples collected during the final treatment weeks from the DC reactor. For polar metabolomics analysis, samples were subjected to a liquid extraction protocol and subsequent instrumental analysis based on the work of Vanden Bussche et al. (2015) and De Paepe et al. (2018). For lipidomics analysis, samples were subjected to a liquid-liquid extraction protocol and subsequent instrumental analysis as described by Van Meulebroek et al. (2017). Prior to analysis, the MS system was calibrated according to manufacturer’s guidelines (Thermo Fisher Scientific, USA) to warrant accurate mass measurements (<5 ppm mass deviation) in both positive and negative ionization mode. In addition, the chromatographic and MS performances were evaluated by injecting standard mixtures that contained the target metabolites. Furthermore, an internal quality control (iQC) sample was assembled by combining equal shares from the extracts that were obtained from the collected samples and employed to condition the LC-MS instrument as well as to perform continuous monitoring of the instrument performance by repeated analysis of the iQC-sample throughout the sequence. A metabolite was positively identified in the LC-MS raw data taking into account the (i) accurate m/z-value of the molecular ion (allowed mass deviation of 5 ppm), (ii) 13C isotope pattern, and (iii) relative retention time (taking into account the retention time of the nearest eluting internal standard with an allowed time deviation of 2.5%). For identification and quantification of metabolites, Xcalibur software (v 4.0.27.21; Thermo Fisher Scientific, USA) was used. In general, a metabolite was considered below the limit of detection/quantification if the metabolite’s observed peak area was below 250,000 arbitrary units (AU). Following peak integration, the area ratio was determined for each metabolite by calculating the ratio between the area of the metabolite and that of the most suited internal standard. Besides the sample-specific normalization that was based on the internal standard, an intra-batch normalization strategy was also applied based on the iQC samples.

To measure the potential gut barrier strengthening and immunomodulatory effects of EpiCor, an intestinal co-culture model consisting of Caco-2 and THP1-cells was exposed to samples collected at the end of the control period and at the end of the treatment period from each colonic reactor vessel of the SHIME^® experiment. Transepithelial electrical resistance (TEER) was used as a measure for gut barrier integrity, while measurements of NF-kB and cytokines [IL-6, IL-10, IL-1b, IL-8, tumor necroses factor (TNF)-alpha, Cys-X-Cys motif chemokine ligand (CXCL) 10, and monocyte chemoattractant protein (MCP)-1] activity were performed to determine immunomodulatory effects. The co-culture was set up and the different measurements performed as previously described (Ghyselinck et al., 2021). As reported by Ghyselinck et al. (2021) experimental control were included in the setup, using a medium control when investigating TEER and a control in which cells were stimulated at the basolateral side with Caco-2 complete medium containing ultrapure LPS (LPS + control) for the other parameters.

Comparison of data obtained on fermentation characteristics for control and treatment periods was performed by means of Wilcoxon Rank Sum tests in R. Differences were considered significant if p < 0.05 and were considered tendencies if 0.05 < p < 0.1.

Before analyzing the microbiome composition data, rare genera with relative abundance less than 10^–5 were removed, decreasing the dataset from 193 to 133 genera. Alpha diversity metrics (i.e., Shannon indices) were calculated in R using the Phyloseq package (McMurdie and Holmes, 2013) using species count tables from Bracken as input normalized using rarefaction. Statistical analysis for alpha diversity metrics was performed with the lme4 package from R, using the fit linear mixed-effect model function lmer (De Boeck et al., 2011) with the formula: ∼ Treatment × Colon region. This model was used to test for differences in the microbiome between control and treatment samples for the PC and DC. Beta diversity was performed with Principal Component Analysis (PCA) of centered log-ratio (CLR) transformed genus counts using the Phyloseq package. Beta diversity was tested with Permutational Multivariate Analysis of Variance (PERMANOVA) using the Adonis function in the vegan package. The model included Aitchison distances (Aitchison et al., 2000) calculated based on CLR transformed values as the dependent variable (Anderson, 2001), assessing associations between treatment, colon region, and their interaction with the composition of the microbiome. Differential abundance analysis was performed with the R package LinDA (Zhou et al., 2022) which was modified to save each individual model for every taxon, CAZy gene family, and KEGG KO, rather than only the coefficient estimates to allow for the use of more complex post hoc contrasts using the R package emmeans (Lenth, 2022). The same models used in the alpha diversity analyses were used in the differential abundance analysis. The emmeans package was used to conduct one post hoc pairwise comparisons after the models were fit, comparing treatments in a pairwise method at each colon region separately. P-values were corrected for multiple hypothesis testing using the Benjamini-Hochberg method (Benjamini and Hochberg, 2018). Overall, differences were considered significant if p < 0.01. To assess ecological dynamics, microbial network analyses were performed with the R package NetCoMi (Peschel et al., 2021) using phylum level microbiome data. First, rare species with relative abundance less than 10^–4 were removed from the analysis. Network construction was normalized using CLR and limited the taxa that was present in a minimum of 75% of samples, using an absolute Pearson correlation ≥0.5. This was done separately for PC and DC samples. Each comparison was done using the “netCompare” function.

The LA-REIMS data were first processed using MassLynx^® V4.1 Progenesis^® Bridge (Waters Corporation, UK) whereby extracted ion chromatograms were created as well as background subtraction was performed. Following ion chromatogram extraction, bridged data were processed by Progenesis^® QI V2.3 (Waters Corporation, UK) in order to list the detected m/z features and their relative abundances. To this end, the automatic sensitivity mode for peak picking (default value) was applied. The abundances of the detected m/z features were normalized based on the total ion abundance. Normalized data were subjected to multivariate statistical analysis using SIMCA^® version 17 (Sartorius, Germany). Here, data were pre-processed, including log-transformation to induce normal distributions and unit variance scaling (1/standard deviation) to standardize the range of signal intensities. Unsupervised principal component analysis (PCA-X) was executed to assess the natural patterning of samples and reveal potential outliers (based on the Hoteling’s T² criterion). Orthogonal partial least squares discriminant analysis (OPLS-DA) was used to differentiate samples according to experimental conditions in a supervised fashion. Validity of the OPLS-DA models was verified by permutation testing (n = 100), cross-validated analysis of variance (p < 0.05), and the quality parameter Q²(Y) (≥0.5).

Processing of the metabolite data obtained using the in-depth LC-MS-based metabolomics approach was founded on a sequential approach of multivariate and univariate statistical analysis. As part of the applied multivariate statistical analysis, PCA-X and OPLS-DA modeling were performed using SIMCA^® version 17. Validity of the OPLS-DA models was assessed based on the same three parameters as used above. If all three parameters were compliant with the set criteria, the model was considered valid and key metabolites that were responsible for the separation between the investigated groups were selected based on the variable importance in projection score (VIP-value), Jack-knifed confidence interval, and S-plot coordinates. Selected metabolites were then subjected to univariate statistical analysis. As a preliminary check for univariate statistics, normality of relative abundances was evaluated using a Shapiro-Wilk test (p < 0.05) and a Quantile-Quantile plot for every metabolite of interest. Since many of the detected metabolites showed deviations from normality, non-parametric tests were used during further data analysis, i.e., Wilcoxon rank-sum tests (p < 0.05).

Finally, to evaluate differences in TEER and immune markers between the control and treatment, a two-way ANOVA with Sidak’s multiple comparisons test was performed using GraphPad Prism version 9.1.2 for Windows (GraphPad Software, San Diego, CA, USA), with differences considered significant if p < 0.05.

Average across weeks within control and treatment periods, acetate levels slightly decreased during treatment compared to the control in both colonic regions (−1.9 and −4.5% in PC and DC, respectively), though only reaching significance in DC (p = 0.012) (Figure 1 and Supplementary Table 1). Butyrate levels, on the other hand, were increased in both colonic regions upon EpiCor supplementation compared to the control period (+ 36.7% and + 5.3% in PC and DC, respectively), with strongest effects observed in PC, reaching significance in that colonic region (p = 0.0004). Administration of Epicor tended to change the propionate production in both colonic regions (p = 0.088 and p = 0.050 in PC and DC, respectively), with increased propionate levels in PC (+ 6.7%), but reduced levels in DC (−10.4%). The bCFA and ammonium concentrations were significantly increased following repeated administration with EpiCor compared to the control period (p = 0.0004 for bCFA in PC and DC; p = 0.0008 and p = 0.0004 for ammonium in PC and DC, respectively).

Average across weeks within the control period, Shannon alpha diversity indices of 2.31 ± 0.09 and 3.12 ± 0.10 were observed in PC and DC regions, respectively (Figure 2A). Following treatment with EpiCor, alpha diversity tended to increase in PC (2.43 ± 0.15; p = 0.057) but was not affected in DC (3.09 ± 0.11; p = 0.0607).

Beta diversity revealed distinct clustering patterns of samples based on PC and DC regions, and between control period and EpiCor treatment (Figure 2B). The first principal components explained 59.2% of the variability, separating samples by region, while the second principal component, which explained 10.6% of the variability, separated samples by treatment. The PERMANOVA analysis confirmed above observations (Region, p = 0.001; Treatment, p = 0.001, Region × Treatment, p = 0.021), demonstrating EpiCor differentially affected microbiome composition in PC vs. DC. Differential abundance analysis (Figure 3) identified significantly (p < 0.01) stimulated bacterial genera following EpiCor supplementation. Significant genera that were affected by a CLR difference > 1 included Bifidobacterium, Phytobacter, Escherichia, Anaerobutyricum, Mediterraneibacter, Eubacterium_I, and Pseudomonas in the PC; Intestinimonas, Harryflintia, Parasutterella, Roseburia, Anaerotruncus, Raoultibacter, NK3B98, QAKS01, 51–20, AF33-28, and CAG-81 in the DC; and Blautia_A, Enterococcus_D, and Agathobaculum in both the PC and DC. On the other hand, the Planococcus and Chryseomicrobium genera were consistently and significantly more abundant during the control compared to the treatment period in both colonic regions, while Sutterella was statistically reduced following repeated EpiCor administration in the PC, with similar effects observed for Alistipes_A, Barnesiella, Butyricicoccus, Eggerthella, Faecalibacterium, Odoribacter, Veillonella, CAG-451, CAG-245, UBA2821, and OF09-33XD in the DC. It was observed that most treatment effects were already established during the first week following EpiCor exposure (Supplementary Figures 2–5). However, while for Roseburia and Anaerotruncus an initial stimulation was observed in the DC following EpiCor supplementation, their abundance reduced during the final weeks of treatment (Supplementary Figures 2–5), indicating gradual establishment of novel cross-feeding mechanisms following repeated administration of the EpiCor postbiotic. Furthermore, when considering the time effects during the SHIME^® experiment (Supplementary Figures 2–5), the observed treatment effects on Pseudomonas and Blautia_A in the PC, and Harryflintia and CAG-245 in the DC were less evident, as a linear increase over time was observed throughout the complete experimental run.

To further explore the effect EpiCor on the structure of microbial communities, microbiome correlation networks were constructed to determine whether treatment with EpiCor was affecting cross-feeding interactions between bacterial groups in the community (Figure 4). This network analysis confirmed that the PC microbiome was significantly more interconnected following treatment with EpiCor as compared to the control microbiome in this colonic area (p = 0.021). A specific trend concerned the enhanced microbial interaction of the Actinobacteria phylum with the microbial community in both colonic regions, especially showing strong interactions with the Firmicutes phylum.

Upon identifying the CAZy potential of the microbiome (Figure 5), it was observed that the profile of enzymes that cleave complex carbohydrates shifted following EpiCor supplementation, with the strongest significant effects observed in PC, while only limited statistical differences were found between the control and treatment period in DC. Stimulated CAZy in PC following EpiCor treatment included glycoside hydrolase families (GH) 4, GH8, GH39, GH59, GH103, glycosyl transferase families (GT) 36 [re-classified in the CAZy database (Drula et al., 2022) to GH94], GT87 and carbohydrate-binding module families (CBM) 73, with strongest effects (indicated by the highest CLR difference) observed for GH4, GT36, and GT87. On the other hand, GT10 was significantly reduced in PC with the same observations for GT14 and GT89 in DC following repeated EpiCor administration.

At the gene level, we found 395 KEGG KOs with significant differences in PC (p < 0.01), while no significant differences were found between the control and treatment period in DC. Gene functions with increase abundance in the PC during the treatment period included 24 ABC transporters, and genes involved in starch and sucrose metabolism, propionate, butyrate, and purine metabolism. Conversely, during the control period in the PC, we observed increase abundance of genes involved in biosynthesis of amino acids, and glycine, serine, and threonine metabolism (Supplementary Table 2).

The PCA-X modeling was performed to assess the natural patterning of the SHIME^® samples and define the analytical performance of the LA-REIMS analysis (Figure 6). The PCA-X model was constructed based on the SHIME^® samples as well as the iQC samples and was composed out of 12 principal components, whereby the first two principal components explained together 64.0% of the present variance (41.2% was enclosed by the first principal component) (Figure 6A). Three sample clusters were noted, of which one was related to the iQC samples and two to the SHIME^® samples, which related to the simulated PC and DC regions. When investigating the colon regions separately, it was observed that PCA-X modeling for the PC samples (6 principal components, with the first two explaining 51.7% of the variance) did not show a clear metabolic impact of EpiCor treatment (Figure 6B). Additional investigation of the treatment effect was achieved through supervised OPLS-DA modeling whereby no valid model could be generated for the comparison of EpiCor treatment with its control period. The PCA-X modeling for the DC samples (6 principal components, with the first two explaining 72.6% of the variance) revealed a clear differentiation between the control and treatment periods (Figure 6C). This was confirmed with OPLS-DA modeling, where a significant metabolic impact of EpiCor treatment was observed (R²X 0.776, R²Y 0.950, Q²Y 0.764, p = 0.0308).

The metabolic impact of EpiCor in the DC was further assessed using UHPLC-HRMS polar and lipidomics profiling (Supplementary Figure 6). A PCA-X model was constructed, whereby metabolic shifts were observed in both the polar (Supplementary Figure 6A) and lipidomics (Supplementary Figure 6B) profiling following treatment with EpiCor compared to the respective control period. The OPLS-DA modeling confirmed these findings (R²X 0.704, R²Y 1, Q²Y 0.974, p = 0.0104 for the polar metabolites; R²X 0.490, R²Y 0.888, Q²Y 0.821, p = 0.0135 for the lipids). Upon assessing the metabolites that were responsible for the observed significant metabolic segregations, univariate Wilcoxon-Mann-Whitney tests showed that the levels of 46 metabolites significantly altered following EpiCor supplementation compared to the control period (Table 1). Of these metabolites, 9 were decreased and 37 were increased following repeated EpiCor administration. Several metabolites from the tryptophan metabolism pathway were significantly altered, with 3-hydroxyanthranilic acid, 3-hydroxykynurenine, kynurenine, nicotinamide, tryptamine, and xanthurenic acid being stimulated following EpiCor supplementation, while the combined measurement of nicotinic acid and picolinic acid was reduced. Metabolic shifts were also observed with respect to carbohydrate homeostasis, with EpiCor treatment resulting in increased levels of 2-aminoadipic acid, lactose, maltose, sucrose and trehalose, and reduced levels of xylose. Furthermore, bile acid metabolism was affected following EpiCor supplementation, as shown by significant reductions of cholic acid (CA) and 7-ketodeoxycholic acid, while ursodeoxycholic acid (UDCA) levels were enhanced. Finally, other metabolites that were significantly affected and showed at least a twofold change following treatment with EpiCor, included 1,3-propanediol, 3,4-dihydroxyphenylacetic acid, 4-hydroxyproline, 4-methyl-2-oxovaleric acid, acetylcarnitine, dehydroisoandrosterone, glutaric acid, hippuric acid, malonic acid, malonylcarnitine, alanine, and sphinganine, with the latter two compounds being elevated in the control period, while the former metabolites were enhanced during the treatment period.

When assessing effects of EpiCor on intestinal barrier integrity in vitro (Figure 7), it was observed that all colonic SHIME^® samples increased TEER compared to the experimental control, keeping its levels at the initial values, with no additional protective effects being observed following repeated EpiCor administration.

Furthermore, most immune markers (except TNF-alpha, MCP-1 and IL-8) were stimulated following administration of the colonic SHIME^® samples (from both the control and treatment period) to the Caco-2/THP1-Blue™ co-culture compared to the experimental control (Figure 7). Treatment with EpiCor only numerically increased NF-kB activity compared to the control period in each colon region. Furthermore, secretion of IL-10 (p = 0.043) and IL-1b (p = 0.004) was significantly stimulated in the PC and DC, respectively, following repeated EpiCor administration.