2.1 Experimental design
The impact of 2 dietary emulsifiers – TWEEN80 (P4780 – suitable for cell culture) and rhamnolipids (R90 – 90% pure, Sigma Aldrich) - on the human gut microbiota was investigated using the mucosal simulator of the human intestinal tract (M-SHIME®), a dynamic semi-continuous reactor-based model that mimics the different stages of the human gastro-intestinal tract (ProDigest-Ghent University). For this purpose, two independent M-SHIME experiments were executed, each inoculated with the faecal microbiota from a different human donor. In each experiment, two concentrations (0,05 m% and 0,5m%) of TWEEN80 and rhamnolipids were tested alongside a control condition, receiving water, in parallel sets of one stomach/small intestinal vessels and one colon vessel (Additional files 1 Figure 1). The mucosal environment was simulated through addition of mucin coated carriers in the colon compartments (Van den Abbeele et al., 2012). During the experiment, the double-walled vessels were continuously stirred by magnetic stirrers (Prosense) at 200 rpm and kept at 37°C through connection to a hot water bath (Julabo). The pH of the colon vessels was continuously monitored using pH electrodes (Consort) and controlled at 6,15 – 6,3 by semi-automatic addition of 0,5 M NaOH or 0,5 M HCl through pH controllers (Consort) and pumps (ProMinent).
2.2 Operation of the SHIME experiments
On the day of start-up, a faecal sample was harvested from a previously selected donor. The donors for these experiments were chosen from a pool of 10 individuals that were consulted for previous experiments (Miclotte et al., 2020). Based on the response of the gut microbiota from these donors to the dietary emulsifiers in terms of short chain fatty acid (SCFA) production, one high and one low responder was chosen. Both donors were female, one 24 and the other 29 years old. Both had not received any antibiotic treatment in the three months prior to the faecal donation. The faecal samples were harvested in plastic lidded containers, which were rendered anaerobic by addition of AnaerogenTM sachets (Oxoid Ltd.), and stored at 4°C for a maximum of 1 hour. A faecal slurry was then prepared as described in De Boever et al., (2000).
Before inoculation, the SHIME-model was operated on water for 24 hours. On the day of inoculation, all tubing and vessels were emptied, the stomach/small intestinal vessels were connected to pancreatic solution, containing 12,5 g/L NaHCO3 (Sigma Aldrich), 0,9 g/L porcine pancreatin (Sigma Aldrich) and 6 g/L Oxgall (BD), and nutritional SHIME feed (Prodigest PD – NM002A)(composition noted in Additional file 1 Table 1), acidified to a pH of 2,0. The colon vessels were filled with 360 mL of non-acidified nutritional SHIME medium, were then inoculated with 40 mL of faecal slurry, manually flushed with 100% N2 gas and left to acclimatize overnight.
The next day, the pumps were turned on and a first set of samples were taken, marking the official start of the experiment (0h). The feed – pumping scheme of the SHIME model is set out in Additional file 1 Table 2. Briefly, the stomach vessels were fed 140 mL of nutritional medium 3 times a day, as well as 60 mL of pancreatic juice, the feeding times 8 hours apart. From the stomach vessels, the solution was pumped to the colon vessels and after a residence time of 16h, discarded.
After at least 5 days of stable operating conditions, which was evaluated by means of pH-stability and visual inspection of the model, the microbial community was considered adapted to the SHIME-environment and the treatment was started. The pH was considered stable when it remained within a range of 6,15 -6,3 and visually, conditions were considered stable if the contents of the colon vessels remained turbid, indicating proper bacterial growth. Treatment solutions were prepared from rhamnolipids and TWEEN80 from Sigma Aldrich in distilled water and 10x more concentrated than the intended concentrations in the vessels. Three times a day, 40 mL of treatment solution was pumped into the respective stomach/small intestine vessels. After a treatment period of one week, the experiment was terminated.
To maintain anaerobic conditions during the experiments, all 14 SHIME-reactors were flushed for 5 min with N2 daily, between 9 and 10 am, except on days where mucin beads were replaced, in which case flushing occurred after bead replacement, around noon. Luminal samples were taken on days 2 (0h), 3, 5, 7, 8, 9, 10, 11 (twice: before and 1h after the start of the emulsifier treatment), 12, 13, 14, 16 and 18 (432h) for SHIME experiment 1 and days 2 (0h), 3, 5, 7, 8, 9, (twice: before and 1h after the start of the emulsifier treatment), 10, 11, 12, 14 and 16 (384h) for SHIME experiment 2. The full sampling scheme can be found in Additional file 1 Table 3 and 4. The procedure for taking luminal samples entailed attaching a 10 mL syringe to the sampling outlet on the lids of the colon reactors of the model, rinsing the syringe and the sampling tube of the colon vessel with the content of the vessel by filling and emptying the syringe multiple times, finally taking a 10 mL aliquot and dividing it in volumes of 1 mL over 1,5 mL Eppendorf tubes and one 10 mL tube for SCFA analysis.
Half of the 60 mucin beads (4 clusters of 15 beads) in each of the colon vessels were replaced every 2 days (except on weekends). First, mucin beads were coated in agar-mucin as described in Van den Abbeele et al., (2012). The beads were then sown together in clusters of 15 and stored in a sterilized container at 4°C until use. In the vessels, the beads were hanged inside a polyamide net, attached to the lid with a nylon wire to facilitate exchange. Beads were always replaced in the morning between 10 and 12 am in the following way. First, the lids of the vessels were connected with the N2 flush system and the vessels would be flushed to ascertain anaerobic conditions during the whole procedure. Next, the vessels were opened one by one, 2 of the 4 clusters in the vessels were replaced by 2 fresh ones, the vessels were closed and flushed with N2-gas for a further 5 min. Meanwhile, the beads were centrifuged for 3 minutes at 500g to collect the mucin.
2.3 SCFA analysis
The SCFA-concentrations were determined by means of diethylether extraction and capillary gas chromatography coupled to a flame ionization detector (Anderson et al., 2017; K. De Paepe et al., 2017). Briefly, 1 mL aliquots were were diluted 2x with 1 mL mili-Q water and SCFA were extracted by adding approximately 400 mg NaCl, 0,5 mL concentrated H2SO4, 400 µL of 2-methyl hexanoid acid internal standard and 2 mL of diethyl ether before mixing for 2 min in a rotator and centrifuging at 3000 g for 3 minutes. Upper layers were collected and analysed using a GC-2014 capillary gas chromatograph (Shimadzu), equipped with a capillary fatty acid-free EC-1000 Econo-Cap column (Alltech), 25 m × 0,53 mm; film thickness 1,2 μm, and coupled to a flame ionization detector and split injector. The resulting SCFA-levels were visualized in R (version 4.1.0) by use of lineplots created using ggplot2 (v3.3.3).
2.4 Total cell counts
The impact of the emulsifiers on total cell concentrations was assessed by first performing SYBR® green staining, on luminal SHIME-samples after which cells were counted on an Accuri C6+ flow cytometer from BDbiosciences Europe. Samples were analyzed every few days in batches from frozen at -20°C. Dilutions up to 10 000x were prepared in 96-well plates using 0,22 μm filtered sterile 0,01 M phosphate buffered saline (PBS) (HPO42- / H2PO4-, 0,0027 M KCl and 0,137 M NaCl, pH 7,4, at 25 °C) and these were subsequently stained with SYBR® green (SG) (100x concentrate SYBR® Green I, Invitrogen, in 0,22 μm-filtered dimethyl sulfoxide) (Props et al., 2016a; Van Nevel et al., 2013). After incubation for 25 minutes, the total populations were measured immediately with the flow cytometer, which was equipped with four fluorescence detectors (530/30 nm, 585/40 nm, >670 nm and 675/25 nm), two scatter detectors and a 20 mW 488 nm laser. The flow cytometer was operated with Milli-Q (MerckMillipore) as sheath fluid. The blue laser (488 nm) was used for the excitation of the stains and a minimum of 10 000 cells per sample were measured for accurate quantification. Settings used were an FLH-1 limit of 1000, a measurement volume of 25 µL and the measurement speed was set to ‘fast’. Cell counts were obtained by gating the total cell populations in R (version 4.1.0) according to the Phenoflow-package (v1.1.6)(Props et al., 2016b). Gates were verified using data from cell-free control samples (0,22 µm filter sterilized 0,01M PBS) (Additional file 2 Figure 1). Lineplots of total cell counts were created using ggplot2 (v3.3.3).
2.5 Amplicon sequencing
Luminal samples from day 2, 11, 12, 14, 18 from SHIME 1 and day 2, 9, 10, 12 and 16 from SHIME 2 and mucosal samples from day 4, 11, 14, 16 and 18 from SHIME 1 and day 4, 9, 11, 14 and 16 from SHIME 2 were selected for Illumina 16S rRNA gene amplicon sequencing. Extraction and quality verification of DNA, library preparation and 16S rRNA gene amplicon sequencing were performed as described in Miclotte 2020.
Recovery of operational taxonomic units (OTU’s) from the amplicon data was carried out using mothur software version 1.40.5 and guidelines (Kozich et al., 2013). First, contigs were assembled, resulting in 9 815 109 sequences, and ambiguous base calls were removed. Sequences with a length of 291 or 292 nucleotides were then aligned to the silva_seed nr.123 database, trimmed between positions 2 and 13 423 (Quast et al., 2013). After removing the sequences containing homopolymers longer than nine base pairs, 2 222 311 (99%) of the unique sequences were retained. A pre-clustering step was then performed, allowing only three differences between sequences clustered together and chimera.vsearch was used to remove chimeras, retaining 73% of the sequences. The sequences were then classified using a naïve Bayesian classifier against the Ribosomal Database Project (RDP) 16S rRNA gene training set version 16, with a cut-off of 85% for the pseudobootstrap confidence score. Sequences that were classified as Archaea, Eukaryota, Chloroplasts, unknown, or Mitochondria at the kingdom level were removed. Finally, sequences were split at the order level into taxonomic groups using the OptiClust method with a cut-off of 0.03. The data were classified at a 3% dissimilarity level into OTUs, resulting in a .shared (count table) and a .tax file (taxonomic classification).
For the entire dataset of 140 samples, 137 458 OTUs were detected in 199 genera. An OTU was here defined as a collection of sequences with a length between 291 and 292 nucleotides and with 97% or more similarity to each other in the V4 region of their 16S rRNA gene after applying hierarchical clustering.
Bioinformatics in R
The shared and taxonomy files resulting from the mothur pipeline were loaded into R for further processing. Absolute singletons (OTUs with only one read over all samples) were removed, resulting in 28 128 OTUs being retained (McMurdie & Holmes, 2014). Rarefaction curves were created to evaluate the sequencing depth (Additional file 2 Figure 2) (Oksanen et al., 2019). Relative and absolute abundances of the OTUs and genera were calculated from the read counts and were explored via bar plots using ggplot2 (v3.3.3) and via principle coordinate analysis (PCoA) on the abundance based Jaccard distance matrix using the prcomp-function in the stats (v4.1.0) package. Relative abundances were calculated as percentages of the total read counts per sample. Absolute abundances (cfr. quantitative microbial profiling) were calculated by multiplying the total cell counts obtained via flow cytometry with the relative abundances of the OTUs (similar to Vandeputte et al., (2017)).
The Chao1, Shannon and Simpson diversity indices were calculated for the microbial community based on the OTU-table using the SPECIES (v1.0) package and the diversity function in the vegan (v2.5-7) package. Indices were plotted using ggplot2.
To investigate the effects of the individual constraints on the microbial community, a series of distance based redundancy analyses (dbRDAs) was then performed on the scores obtained in the principle coordinate analysis (PcoA) on the Jaccard distance matrix using the capscale function in the vegan (v2.5-7) package. Permutation tests were used to evaluate the significance of the models and of the explanatory variables (De Paepe et al., 2018). DbRDAs were created for both the relative and absolute abundance dataset.The global model included the factors Treatment (a concatenation of the factors Emulsifier and Emulsifier concentration) Timepoint, and Donor as explanatory variables and the absolute abundances of the genera as explanatory variables. In three subsequent dbRDAs, the constrained variance attributed to each of the factors was investigated by each time conditioning out all but one factor. The results of the dbRDAs were plotted as Type II scaling correlation triplots showing the two first constrained canonical axes (labeled as dbRDA Dim 1/2) and the proportional constrained eigenvalues representing their contribution to the total (both constrained and unconstrained) variance. The sites were calculated as weighed sums of the scores of the response variables.
Finally, to evaluate which genera were likely differentially abundant between emulsifier treatments and control, the DESeq2 package (v 1.32.0) was applied on the read count-table at genus level. To this end, only the read counts from the samples in the treatment period were used.
To streamline the DESeq-process, pre-filtering according to McMurdie & Holmes, (2014) was first applied on the read count-table, after which a genus-level table was created using the aggregate function (stats package v4.1.0). In the generalized linear model, the factor Timepoint, Donor, and Treatment were included. A likelihood ratio test was employed within the DESeq function on the reduced model, containing only the factors Donor and Timepoint, to test for the significance of the model. Low count genera were subjected to an empirical Bayesian correction (Love et al., 2014). For pairwise comparison of TWEEN80 and rhamnolipid treatments versus controls, Wald tests were used after shrinkage of the Log2FoldChange (L2FC) values by means of the lfcShrink function with the adaptive shrinkage estimator from the ashr-package (Stephens, 2017). P-values were adjusted by means of a Benjamini-Hochberg procedure (Love et al., 2014). Results were visualized in volcanoplots, displaying the −log(adjusted p-value) versus the Log2FoldChange of each genus. Additionally, box plots were created showing the log-transformed pseudocounts extracted by the plotCounts function for each genus that showed significant differential abundance.
2.6 Metabolomics
Both targeted and untargeted polar metabolomics analyses of luminal SHIME-samples from days 2, 11, 12 ,14 and 18 from SHIME 1 and days 2, 9, 10, 12, 16 from SHIME 2 were performed in collaboration with the Lynn Vanhaecke lab according to previously described protocols (E. De Paepe et al., 2018; De Spiegeleer et al., 2020; Vanden Bussche et al., 2015).
Polar metabolomics analysis
Polar metabolites were first extracted from the frozen SHIME-samples by means of ultrapure water (0,055 mS cm -1) obtained via a purified water system (VWR International, Merck). For this purpose, 1,5 mL SHIME-sample was first spiked with 75 µL internal standard solution (100 ng/µL L-alanine-d3 and D-valine-d8) and then centrifuged for 5 min at 13 300g. The supernatant was then filtered over a 0,22 µm polyvinylidene difluoride filter (Millipore), diluted 1:5 using ultrapure water and transferred into an HPLC-vial (Vanden Bussche et al., 2015).
Polar metabolites were then separated and detected using an ultra-high performance liquid chromatograph system coupled to a Q-ExactiveTM mass spectrometer with heated electronspray ionisation (UHPLC-HRMS Thermo Fisher Scientific) according to validated methods as described in De Paepe et al 2018 and De Spiegeleer et al 2020. Separation was achieved using a Dionex Ultimate XRS 3000 system (Dionex) equipped with a Acquity UPLC HSS T3 column (1,8 mm particle size, 2,1 mm internal diameter, 150 mm length, Waters) kept at a constant 45°C during analysis. Samples were injected in randomized order at a volume of 10 µL and eluted using a polar to apolar gradient formed using ultrapure water (A) and LC/MS grade acetonitrile (Thermo Fisher Scientific) (B), both acidified with 0,1% LC/MS grade formic acid (Thermo Fisher Scientific). Elution occurred at a constant flow rate of 0,4 mL/min. The gradient profile comprised the following proportions of solvent A: program: 0 – 1,5 min at 98%, 1,5 – 7,0 min from 98% to 75%, 7,0 – 8,0 min from 75% to 40%, 8,0 – 12,0 min from 40% to 5%; 12,0 – 14,0 min at 5%, 14,0 – 14,1 min from 5 to 98%, followed by 4,0 min of re-equilibration (E. De Paepe et al., 2018; Vanden Bussche et al., 2015).
Subsequent detection of the metabolites occurred using a Q-exactive™ Orbitrap mass spectrometer (Thermo Fisher Scientific), with a heated electrospray ionization source (HESI – II) operated in polarity switching mode (0/B/1 position). The system was operated under the following settings: m/z range of 53,4 - 800 Da, a maximum injection time of 70 ms, a mass resolution of 140 000 full width at half maximum and an automatic gain control of 106 ions. The sheath, auxiliary and sweep gas (N2) flow rate for the HESI-II source were set at 50, 25 and 3 arbitrary units. The heater and capillary temperatures were set at 350 °C and 250 °C respectively, the S-lens RF-level was 50 V and the spray voltage was 4 kV (De Spiegeleer et al., 2020).
Before analysis of the samples, both ionization modes of the Q-Exactive system were calibrated using ready-to-use calibration mixtures (Thermo Fisher Scientific) for optimal mass accuracy. Sample analysis was also preceded by injection of a standard mixture followed by injection of a series of four quality control samples – obtained by pooling equal aliquots from all samples together – to allow system stabilization. Duplicate quality control samples were also added intermittently between series of 10 samples, to allow for correction of instrument instability and signal drift.
Data Processing
Targeted metabolomics
After data acquisition, peak identification and quantification was performed within Xcalibur 3.0 software (Thermo Fisher Scientific). Metabolite peaks were manually identified using an in-house retention time library, the peak location in the standard mixture and the C13/C12 ratio according to CD 2002/657/EC guidelines. A minimal peak intensity threshold of 500 000 a.u. was employed. Peak areas were gathered as one table in Excel 2016 and normalized through division by peak area values for the internal standard (D-valine-d3). This base table containing area ratio’s was then plugged into R (version 4.1.0) for further principle coordinate analysis or split up per treatment versus the control condition for analysis in MetaboAnalyst 5.0 (2021).
Disease signatures present in the metabolite dataset were explored using the quantitative Enrichment Analysis application in MetaboAnalyst 5.0 (Xia Lab, McGill University). Prior to these analyses, log transformation and Pareto scaling were performed via the website tool. The faecal library was chosen as reference database. The Statistical Analysis application on MetaboAnalyst was also employed to identify differential compounds. Compounds were considered differential when they were indicated as significantly different by the Wilcoxon Rank-sum test, the Variable Importance of Prediction (VIP) value was over 1 in partial least squares discriminant analysis (PLS-DA) and if p(1) > 0,1 and p(corr)[1] > 0,4 in orthogonal partial least squares discriminant analysis (OPLS-DA) for treatment with both concentrations with TWEEN80 and rhamnolipids.
Untargeted metabolomics
Compound Discoverer 3.0 (Thermo Fisher Scientific) was used first for automated peak extraction, peak alignment, deconvolution and noise removal with settings as given in Table 1 (E. De Paepe et al., 2018; De Spiegeleer et al., 2020). Suspect compounds were identified based on the Chemspider database (Royal Society of Chemistry).
Table 1: Settings utilized in Compound Discoverer for peak extraction and compound detection.
Settings in Compound Discoverer 3.0
|
Value
|
Peak intensity threshold
|
500 000 a.u.
|
Specific adducts
|
[M+H], [M+K],[M+Na], [M-H], [M-K], [M-Na]
|
Retention Time
|
0,5 - 16
|
S/N threshold
|
15
|
Max peak shift
|
0,2 min
|
Mass Tolerance
|
<5 ppm
|
Max peak width
|
0,55
|
Peak area data obtained from Compound Discoverer was plugged into SIMCA 17 (Umetrics) for multivariate statistical analysis. To this end, first QC-normalisation and QC-filtering based on a CV < 30% were performed to eliminate instrumental signal drift and metabolites for which detection was unstable. QC-normalisation was performed per compound, by dividing the peak area in a sample by the means of the interspersed duplicate QC-samples that follow the batch of 10 samples that sample belongs to.
In SIMCA, a principle component analysis (PCA) was first generated as an exploration of overall data patterns. Next, Orthogonal Partial Least Squares Discriminant Analysis (OPLS-DA) was used to generate predictive models comparing emulsifier treatment conditions to the emulsifier free control condition. For each of these comparisons a model with a Q2 > 0,5, R2Y ≈ 1 and R2X ≈ 1 were selected, of which the significance was verified by ANOVA and Permutation tests. Model creation involved compound-filtering based on the Variable Importance of Projection criterion: VIP >1. Compounds characteristic for the particular emulsifier treatment were then identified based on their VIP (>1), a positive Jackknife confidence interval and the S-plots (|p(1))| > 0,05 and relevance |p(corr(1))| > 0,4).