Study design
Three clinical studies were conducted to assess the effects of fructooligosaccharides (FOS) and polydextrose (PDX) on the gut microbiome of healthy male and female human participants age 18–45 years. FOS used in the studies is a brand of FOS called Orafti® P95 dry powder manufactured by BENEO. PDX used in the studies is a brand called Litesse® Ultra™ manufactured by Danisco. The first two studies (K001 and K002) had similar designs except that K001 used FOS as the study product, and K002 used PDX as the study product (Additional file 1: Figure S1). Both were randomized, open-label studies that included three cohorts, which received different feeding amounts of the study products. Following a 2-week baseline/run-in period, participants were randomized via computer-generated codes on study day 0 to one of the three cohorts. Participants consumed the study products during two feeding periods of 2 weeks each (feeding1 and feeding2). There was a 4-week washout period following the feeding periods during which participants were monitored but no study product was consumed. In K001, participants in cohort 1 were instructed to consume 2.5 g BID of FOS per day during feeding1 and feeding2. In cohort 2, 2.5 g BID during feeding1 and 5 g BID during feeding2. In cohort 3, 5 g BID during feeding1 and 10 g BID during feeding2. In K002, participants in cohort 1 were instructed to consume 7.5 g BID of PDX per day during feeding1 and feeding2. In cohort 2, 10 g BID during feeding1 and 20 g BID during feeding2. In cohort 3, 20 g BID during feeding1 and 30 g BID during feeding2. The study products were consumed twice daily, dissolved in water. The study product amounts were comparatively greater in K002 as compared to K001 due to reported higher tolerability for PDX than FOS [16, 26]. A third randomized, double-blinded crossover study (K003) was conducted to assess whether effects on the microbiome were reproducible and/or dependent on the order of dosing. Following a 2-week run-in period, healthy participants were randomized into one of two cohorts on study day 0. Participants in cohort 1 received 5 g BID of FOS during the first 2-week feeding period (feeding1), discontinued consumption of FOS during the subsequent 4-week period (washout1), received 25 g BID of PDX during the second 2-week feeding period (feeding2), and were followed for a second 4-week period (washout2) after discontinuation of the consumption of PDX. Cohort 2 followed an identical plan except that participants in the cohort received the study products in the reverse order across the two feeding periods. Study products were consumed as in the two earlier studies. The differences between the intended and actual consumption amount are summarized in Additional file 1: Figure S2. Participants were instructed to maintain their normal diet throughout the study. For each study, we randomized 10 participants per cohort. The demographic information of participants is listed in Additional file 1: Table S1. Participants with a major protocol deviation or completely missing data from one of the study periods were excluded from downstream analyses.
Sample collection and sequencing
Fecal samples were collected in order to characterize the gut microbiome of the participants. Participants were provided with fecal sample collection instructions. Briefly, participants were instructed to place a stool collection bowl and holder under toilet seat in the center rear of the toilet and lower seat, to ensure that the collection unit was securely in place. Participants then completed bowel movement on the toilet, into the sample collection unit. To generate the swab sample, participants opened the red capped swab collection tube (BD BBL CultureSwab EZ II single head foam swab), stuck swab into the center of the stool to heavily coat swab head, returned swab back into the collection tube, and twisted cap back onto collection tube to seal completely. Swab samples were then placed immediately in the freezer or were temporarily placed in the provided insulated cooler bag with four frozen freezer packs until participants could access their freezer. To bring collected samples back to the clinical site, participants placed fecal swabs in the provided insulated cooler containing the freezer packs.
Fecal microbiomes were characterized using a shotgun DNA sequencing-based approach similar to one previously described [27]. Briefly, fecal swabs were removed from the collection device quickly and cut into Qiagen’s DNeasy PowerSoil extraction plates. Extraction plates remained on dry ice until all fecal swabs were transferred. Plates were stored at – 80 °C until analysis. Extraction plates were shipped on dry ice to CoreBiome, Inc. for downstream extraction and library preparation. Extracts were quantified using the Quant-iT PicoGreen dsDNA assay (Thermo Fisher). Libraries were prepared using the NexteraXT kit and a HiSeq 1 × 150-cycle v3 kit (Illumina) was used to sequence samples. The taxa count tables from raw sequencing data were generated using the SHOGUN pipeline [27] and have been made publicly available (see the “Availability of data and materials” section). The database used here was generated by selecting up to the first 20 strains per species in RefSeq v87 by first choosing genomes with assembly level annotated as “Complete Genome,” then “Chromosome,” then “Scaffold,” and then “Contig.” All calculations and visualizations were conducted on an operational taxonomic unit (OTU) table that excluded OTUs that could not be resolved to the species level. This table was rarefied to 10,000 without replacement and used for downstream analyses.
Microbiome data analysis
We assessed the effects of FOS or PDX on microbiome diversity by comparing median alpha (Shannon) diversity values between each of the feeding periods and the baseline/run-in period and testing the significance of differences via paired Wilcoxon signed-rank tests. P values were false discovery rate (FDR) corrected for multiple comparisons across periods.
To assess the overall effect of FOS or PDX on microbiome community composition, we first calculated Bray-Curtis dissimilarities between every sample pair within each participant based on the square-root transformed relative abundances of the bacterial species. We determined whether the compositions during the feeding periods differed from the baseline/run-in period using permutational analysis of variance (PERMANOVA) as implemented by the “adonis” function in the vegan package [28]. Participant identity was included for the “strata” parameter to limit permutations to within participants. We next assessed FOS or PDX’s effect on community dissimilarity. We tested the differences between the dissimilarities between one subject’s every sample from baseline/run-in and every sample from another period and the dissimilarities within samples from baseline/run-in using Kruskal-Wallis test followed by Dunn’s post hoc test (as implemented in the dunn.test R package [29] with Benjamini-Hochberg correction).
To evaluate whether the high-dose cohort induced higher community shift than the lower dose cohorts, we again looked at the pair-wise dissimilarity between every baseline/run-in sample and every feeding-period sample across three cohorts and tested their differences between cohorts using a Kruskal-Wallis test followed by Dunn’s post hoc test with Benjamini-Hochberg correction.
We assessed whether individual taxa significantly differed between feeding periods and the run-in period using a linear mixed effect model as in [30]. Specifically, for each taxon, we rank transformed its relative abundance data and fit a linear mixed effect model (using the lme function from nlme R package [31]) with period as a fixed effect and participant as a random effect. The test results for all taxa were corrected using FDR. Only taxa with relative abundances greater than or equal to 0.1% in either the run-in or feeding periods were examined. All aforementioned analyses were performed using R version 3.5.3 [32].
Differential abundance analysis of genes encoding CAZymes
We obtained functional information for each species using a custom annotation pipeline derived from [33]. For each species under analysis, all genomes in NCBI corresponding to that species were downloaded. We used the dbCAN database for CAZymes [34, 35]. The annotation pipeline was applied to the genomes, resulting in a map from strains to CAZymes. A map from species to functional units was created by merging the results for all strains within each species. For each species, we thus obtained a set of CAZymes. For each CAZyme, the abundances of species that passed filtering criteria (see the “Initial filtering of taxa” section below) in each participant and annotated as having the gene were aggregated together and log2 fold changes (baseline versus feeding period) were calculated. Aggregated fold changes were similarly calculated for taxa not annotated to have the gene. These fold changes were then compared using the two-sided Wilcoxon rank-sum test (with p values adjusted for multiple hypothesis testing using the Benjamini-Hochberg method and FDR < 0.01.) To be precise: Consider a CAZyme c. For each subject s, let As,i, i = 1,...Ns denote the series of Ns timepoints of aggregated relative abundance of all species in subject s that are annotated as having CAZyme c, and let Bs,i similarly denote the time series of aggregated relative abundance of all species in s not annotated to have CAZyme c. Fold changes for taxa annotated as having the gene are calculated according to the following:
FCA(s) = median(As,i for i in the feeding periods)/median(As,i for i in the baseline periods).
FCB(s) is calculated analogously. The Wilcoxon test is then performed on the set of log2 FCA(s) versus the set of log2 FCB(s).
MC-TIMME2 model and software
Software
MC-TIMME2 was implemented in Python 3 using the Scipy ecosystem and Numba [36,37,38,39]. The software and source code to reproduce analyses in the manuscript is publicly available under the GNU General Public License [40]. Configuration settings for the model and inference are specified in a text-based configuration file, and the software is run from the command line. Data is loaded from three files: an OTU table containing the counts for each taxon in each sample, a Sample Info file containing the time and participant of each sample, and a Treatment Info file containing the time and quantity of each administered dose of the compound. Additionally, a phylogenetic tree for the taxa being analyzed must be provided. In order to learn the parameters of the sequencing measurement noise model, replicate data must be provided. This consists of another OTU table containing the count data for a set of technical replicates. Here we used 10 technical replicates of shotgun sequencing from a single human stool sample. DNA was extracted using Qiagen’s PowerFecal kit and DNA quantification, library preparation, and sequencing were done at CoreBiome, Inc. as described earlier. Inference for the parameters of the measurement noise model was performed upfront and separately from the rest of the model.
As the software performs posterior inference, it produces an output file containing the posterior samples for all parameters in HDF5 format. Once inference is complete, the software includes functionality to generate output files from this file, including visualizations for each time-series, reconstructed perturbation groups, and text-based tables showing the carrying capacities, pharmacokinetic parameters, onset times, and cluster memberships for each time-series.
Initial filtering of taxa
Before running the model, we apply two initial filters to remove very low abundance taxa trajectories. The first filter removes any species which does not achieve a mean abundance (across the trajectory) of at least 0.0005 in at least 25% of the participants. The next filter is applied to each individual trajectory and removes any trajectory that does not achieve at least 15 counts in any three consecutive timepoints.
Model details
Double Dirichlet process clustering
Two levels of Dirichlet process mixture model clustering are employed. Individual species-participant time-series are probabilistically assigned in the model to microbe groups, with group assignments denoted by zso for species o in participant s. Each microbe group k is itself assigned either to a null perturbation effect or to a perturbation group in the top level of the double Dirichlet process. Let ek denote the assignment to a perturbation group for microbe group k. Phylogenetic information is incorporated into clustering using a potential function ψ:
$$ \psi \left(o,{\lambda}_k\right)=\exp \left(-{\zeta}_0d\left(o,{\lambda}_k\right)+{\zeta}_1\right) $$
(1)
Each microbe group probabilistically selects a representative member λk. Then, the likelihood for each time-series being assigned to the microbe group k incorporates the phylogenetic distance between that species and the representative. Letting Ωm indicate the perturbation parameters for perturbation cluster m, we write the probability of assigning species o in participant s as follows:
$$ P\left({z}_{so}=k|{x}_{so},{e}_k,{\lambda}_k\right)\propto P\left({z}_{so}=k\right)P\left({x}_{so}|{\Omega}_{e_k}\right)\ \psi \left(o,{\lambda}_k\right) $$
(2)
Pharmacokinetic model
The inputs to the model include the participant-specific doses dsi and their times of administration asi. The level of the compound over time, cs, is modeled using first-order pharmacokinetics, with a participant-specific elimination rate κs.
$$ {c}_s(t)={\sum}_i{d}_{si}\exp \left(-{\kappa}_s\left(t-{a}_{si}\right)\right)I\left(t>{a}_{si}\right) $$
(3)
Dynamical model
We use a stochastic logistic growth model for baseline microbial dynamics. In the absence of a perturbation, each species is assumed to follow a stochastic logistic growth model with a growth rate α, a self-limiting term β, and a process variance σso2 specific to each time-series so.
$$ {\mu}_{so}\left({t}_{s,i}\right)={x}_{so}\left({t}_{s,i-1}\right)+{\alpha}_{so}{x}_{so}\left({t}_{s,i-1}\right){\Delta}_{s,i}+{\beta}_{so}{x}_{so}{\left({t}_{s,i-1}\right)}^2{\Delta}_{s,i} $$
(4)
$$ {x}_{so}\left({t}_{s,i}\right)\sim N\left({\mu}_{so}\left({t}_{s,i}\right),{\Delta}_{s,i}{\sigma}_{so}^2\left({t}_{s,i}\right)\right) $$
(5)
The process variance takes different values on and off the perturbation period, allowing the model to capture changes in variance which may accompany the administration of the compound.
$$ {\sigma}_{so}^2(t)=\left\{\begin{array}{c}{c}_{so}\kern1.25em \mathrm{t}\ \mathrm{off}\ \mathrm{perturbation}\\ {}{\hat{c}}_{so}\kern1.5em \mathrm{t}\ \mathrm{on}\ \mathrm{perturbation}\end{array}\right. $$
(6)
Auxiliary trajectory
Inference efficiency depends on the use of a normal distribution (not truncated) in the stochastic dynamical Eq. (5). However, it is desirable to avoid negative values for the trajectory x, which are physically unrealistic and unstable. This is achieved using a relaxation method which we previously developed [41]. This method introduces a strictly positive auxiliary trajectory which is closely coupled to the dynamical trajectory. The effect of this coupling is to tend to keep x away from negative values without forcing a computationally intractable hard constraint.
Model of dose-dependent perturbation
The baseline growth rate αso in (4) is modulated by a perturbation as follows:
$$ {\alpha}_{so}^{\prime}\left({t}_{s,i}\right)={\alpha}_{so}\ \left(1+{\gamma}_{1m}{w}_{sm}\left({t}_{s,i}\right)h\left({c}_s\left({t}_{s,i}\right);{r}_m\right)+{b}_{so}{\gamma}_{2 so}{p}_{so}\left({t}_{s,i}\right)\right) $$
(7)
In (7), w is a step function specifying when the perturbation is active; w is defined by an on-time in days and a duration relative to the dose administration period. Here, h is a sigmoidal transfer function which reshapes the participant-specific pharmacokinetic trajectories cs into a perturbation magnitude over time.
$$ h\left({c}_s(t);{r}_m\right)=\frac{2}{1+\exp \left(-{r}_m{c}_s(t)\right)}-1 $$
(8)
Each perturbation cluster also has the parameter γ1m controlling the strength of the perturbation.
The parameters bso, γ2so, and pso allow for a persistent effect in which the growth rate continues to be modulated even after the compound is excreted. Occurrences of persistent effects are probabilistically selected to be on or off, and when they are selected to be on, the parameters are specific to an individual time-series. Here, bso is the indicator variable that probabilistically selects whether this effect is on, with the prior on bso favoring no persistence. Here, γ2so is the magnitude of the persistent effect and pso is a step function that turns on at the last dose administration time and continues for tsopersist days.
Measurement model
The observed data are sequencing counts y. The counts are assumed to be drawn from a negative binomial distribution, which has been used extensively to model microbiome sequencing counts data [42]. In order to model the noise properties of the shallow shotgun sequencing data used in these studies, we used a mixture model for the replicates data (separate from the main model of microbiome responses). The replicates data consists of counts yio for species o in replicate i, with each species assumed to have one latent mean relative abundance μo across all replicates. Thus, the sample-specific mean ηio equals νi multiplied by μo. For the overdispersion, we use a Dirichlet process mixture model of negative binomial components, where each component has a single overdispersion parameter. Therefore, conditional on indicator variables zo,
$$ {y}_{io}\mid {z}_o\sim NB\left({\eta}_{io},{\epsilon}_{z_o}\right) $$
(9)
Inference is performed for this model, resulting in a consensus set of negative binomial components, each with an overdispersion value and a mixture weight. To use this result in the main model, assignments to the consensus clusters are marginalized out:
$$ {y}_{so}\left({t}_s\right)\sim {\sum}_j{\pi}_j^{\mathrm{noise}} NB\left({\eta}_{so}\left({t}_s\right),{\epsilon}_j\right) $$
(10)
Phylogenetic tree
MC-TIMME2 takes as input a phylogenetic tree between all species in the model. For each species in the database used for taxonomic calling described earlier, the aligned 16S rRNA sequences for all strains corresponding to that species were selected from RDP [43]. Several species in the database were not present in RDP. For those species, unaligned 16S rRNA sequences were obtained from NCBI RefSeq [44] and aligned to the RDP alignment using QIIME/PyNAST [45]. Using this alignment of 16S rRNA sequences, FastTree was used to generate a phylogenetic tree [46].
Comparison of onset times
For comparison of onset times, we restricted analyses to time-series with strong evidence of response (Bayes factor > 100). Additionally, for clarity, we tested only pairs of bacterial families if each family had at least 10 responding time-series. We used the Mann-Whitney test to compare the two sets of posterior median onset times. P values are corrected for multiple comparisons using the Benjamini-Hochberg procedure.
Consensus clustering
To obtain a single consensus clustering from the MCMC samples, we used an agglomerative clustering method similar to that used in [47]. Reconstruction was applied to those time-series for which there was strong evidence of a response (Bayes factor > 100). First, the microbe co-clustering frequencies were used to form an affinity matrix between each time-series. Agglomerative clustering up to the mode of the number of microbe groups was performed to create a set of consensus microbe groups. To form consensus perturbation groups, an affinity matrix between the consensus microbe groups is formed, with the affinity between each pair of consensus microbe groups given by the average of all pairwise perturbation co-clustering frequencies among the members of the consensus microbe groups.
Visualization and filtering of perturbation groups
For visualization, we applied post hoc filtering to the consensus perturbation groups. First, we filtered species, keeping only those species which appeared in any perturbation group at least 6 times (corresponding to a response in 25% (FOS) or 23% (PDX) of participants). Next, we filtered out perturbation groups which did not have among their members at least 15% of participants represented. This filtered set of perturbation groups was used for visualization (Fig. 4). For functional enrichment analysis (Additional file 1: Figure S6), all species that were originally in the group were used, including those species that were filtered out in the first step for visualization. In Fig. 4, we applied an ordering to the phylogenetic tree that tends to place the more frequently responding species near the top. Each species in the tree was associated with an importance score given by the number of times a time-series of that species appeared in any perturbation group. Each slot in the grid was numbered in reverse order (the bottom slot has value 1, the second-to-bottom slot has value 2, and so on). We then maximized the dot product of the vector of species importance scores and the vector of slot values over valid configurations of the phylogenetic tree using Monte-Carlo optimization.
Model inference, full mathematical description, and sensitivity analyses
MC-TIMME2 performs approximate posterior inference using a custom Markov chain Monte Carlo (MCMC) method. Full mathematical details of the inference algorithm and model, as well as descriptions of sensitivity analyses of hyperparameters, are given in Additional file 2: Supplemental Methods.