Gastrointestinal microbiota composition predicts peripheral inflammatory state during treatment of human tuberculosis

Abstract

The composition of the gastrointestinal microbiota influences systemic immune responses, but how this affects infectious disease pathogenesis and antibiotic therapy outcome is poorly understood. This question is rarely examined in humans due to the difficulty in dissociating the immunologic effects of antibiotic-induced pathogen clearance and microbiome alteration. Here, we analyze data from two longitudinal studies of tuberculosis (TB) therapy (35 and 20 individuals) and a cross sectional study from 55 healthy controls, in which we collected fecal samples (for microbiome analysis), sputum (for determination of Mycobacterium tuberculosis (Mtb) bacterial load), and peripheral blood (for transcriptomic analysis). We decouple microbiome effects from pathogen sterilization by comparing standard TB therapy with an experimental TB treatment that did not reduce Mtb bacterial load. Random forest regression to the microbiome-transcriptome-sputum data from the two longitudinal datasets reveals that renormalization of the TB inflammatory state is associated with Mtb pathogen clearance, increased abundance of Clusters IV and XIVa Clostridia, and decreased abundance of Bacilli and Proteobacteria. We find similar associations when applying machine learning to peripheral gene expression and microbiota profiling in the independent cohort of healthy individuals. Our findings indicate that antibiotic-induced reduction in pathogen burden and changes in the microbiome are independently associated with treatment-induced changes of the inflammatory response of active TB, and the response to antibiotic therapy may be a combined effect of pathogen killing and microbiome driven immunomodulation.

Trial registration: ClinicalTrials.gov NCT02684240.

Conflict of interest statement

M.F.W. is currently an employee and shareholder of Regeneron Pharmaceuticals, Inc. M.S.G. reports consulting fees and equity in Vedanta Biosciences, Inc., consulting fees from Takeda, and is on the SAB of PRL-NYC. V.B. is supported by a Sponsored Research Agreement from Vedanta Biosciences, Inc. The remaining authors declare no competing interests.

Figures

Fig. 1. Overview of cohorts, subjects, timepoints, samples, and hypotheses in this study.

A This study investigates microbiome-transcriptome relationships in three separate cohorts of individuals in Haiti. Cohort 1 (2-week longitudinal and interventional clinical trial) consists of secondary analysis of a randomized clinical trial of study volunteers, where we collected disease severity measurements (Mtb bacterial load, TTP), microbiome profiling, and peripheral transcriptomics in active TB patients at baseline, before randomization to either HRZE (arm 1) (standard of care TB treatment), or Nitazoxanide (NTZ) (arm 2). Cohort 2 (6 month longitudinal and observational study) consists of study volunteers who were followed throughout the course of 6 months of TB treatment, where we collected TTP, microbiome, and transcriptomics data. Finally, Cohort 3 (cross sectional and observational) consists of healthy volunteers. These healthy volunteers were enrolled separately. Around half are healthy and TB-negative household or family contacts (FC) of active TB-patients, and the other half are community controls (CC), with no know TB exposure. We performed microbiome profiling and peripheral transcriptomics on these individuals as well. B Numbers of individuals in this study. C Diagram showing the major questions investigated in this study. Supplementary Data 1 provides a table with dates of first and last enrollment for every Cohort.

Fig. 2. Both HRZE and NTZ perturb the gut microbiota after two weeks of therapy, but only HRZE reduces M. tuberculosis bacterial load.

A Schematic of the clinical trial comparing bactericidal effect of HRZE and NTZ. B Paired M. tuberculosis sputum time to positivity (TTP) at day 0 and day 14 for the NTZ treatment arm and HRZE treatment arm. Data are displayed as range (minimum and maximum) of two-three technical replicates; n = 16 biologically independent individuals for the HRZE arm and n = 19 biologically independent individuals for the NTZ arm. Linear mixed effect modeling was used to determine significance of difference in post/pre treatment in each arm as TTP~1+Sex+Age+Treatment+Time+Treatment:Time+1∣ID, where Treatment indicates the arm (NTZ or HRZE), Time indicates Pre or Post antibiotic administration, and : indicates the interaction term. For each individual we use as TTP measurement the average of the relative technical measurements. NTZ treatment is associated with no difference in TTP between day 0 and day 14 (p > 0.05 for the coefficient of variable Time, see Supplementary Data 2), whereas HRZE significantly reduces bacterial load (p < 0.05 for the coefficient of variable Treatment : Time see Supplementary Data 2). Data for TTP were obtained from Walsh et al.. C Principal Coordinate analysis (PCoA) with Bray–Curtis distance showing differences in microbiome community structure between individuals before and after 14 days of either HRZE or NTZ treatment. The gray line connects baseline and day 14 treatment paired samples. PCoA1 clearly discriminates samples post NTZ treatment (pink triangles) from those at baseline or after HRZE treatment. PERMANOVA analysis was used to reject (p = 0.001, see Supplementary Data 4) the null hypothesis that the centroids and dispersion of the groups (pretreatment, after NTZ and after HRZE) are equivalent for all groups (see Supplementary Data 4). D Microbiota alpha diversity plotted using the Inverse Simpson Diversity index; n = 16 biologically independent individuals for the HRZE arm and n = 19 biologically independent individuals for the NTZ arm. Linear mixed effect modeling was used to determine the significance of difference of treatment on diversity. We fitted the model Diversity~1+Sex+Age+Treatment:Time+1∣ID. The symbol: indicates the interaction term. HRZE was used as the reference level. No significant difference between the two treatment at baseline was observed. Both groups (p < 0.05 for the coefficient of variable Time corresponding to HRZE treatment and p < 0.05 for the coefficient of variable Treatment:Time corresponding to NTZ treatment in this model) display significantly reduced Inverse Simpson diversity after 14 days of treatment (see Supplementary Data 3). Source data are provided as a Source Data file.

Fig. 3. Overlapping and distinct microbiome perturbation induced by NTZ and HRZE.

A Volcano plots indicating the post (day 14) vs pretreatment (baseline) differences at the ASV level for HRZE and NTZ. The color of each ASV is according to the phylogenetic Order. A single linear mixed effect model for each ASV of the form ASVicounts~Sex+Batch+Group+1∣ID was fitted to determine differences due to treatment while accounting for sequencing batch and sex. ASVs significantly affected by the treatment were those determined to have a Benjamini-Hochberg false discovery rate (FDR) adjusted p-value less than 0.05 for the variable Group in the limma/voom model (see “Methods” section). The horizontal dotted lines indicate FDR < 0.05 and vertical dotted lines indicate |log2FC| > 1.5. B Within-arm unsupervised hierarchical clustering of the abundances of 404 ASVs found to be significantly affected by HRZE or NTZ treatment (FDR < 0.05, see Supplementary 5). The heatmap columns are split by arm membership (including baseline randomization group), and the heatmap rows are split by ASV phylogenetic Phylum, and within the Phylum, the Order is colored as in A. The right annotations (HRZE and NTZ) indicate whether each ASV was significantly perturbed by treatment. P value in y axis is adjusted according to Benjamini–Hochberg (FDR). Source data are provided as a Source Data file.

Fig. 4. Hallmark pathway gene set enrichment analysis and gene expression comparison in HRZE and NTZ treated arms.

A Hallmark gene pathway changes associated with 2 weeks of HRZE (A) or NTZ (B). Positive are pathways overrepresented at 2 weeks of therapy (up), and negative are pathways underrepresented at 2 weeks (down), both compared to baseline. All pathways are significant (FDR < 0.05, see Supplementary Data 6) with the size of the arrow indicating level of significance. Only pathways from MiSigDB Hallmark pathway set found to be significantly altered in this analysis are shown. Here NES stands for Normalized Enrichment Score. B, C TB-associated peripheral blood transcripts from ref. , highlighting post treatment vs baseline changes in gene expression for HRZE (B) or NTZ (C). Notably, HRZE renormalizes (i.e., towards a healthy control state) the expression of 144 validated TB inflammatory transcripts and exacerbates only 13 (see Supplementary Data 7). NTZ is found only to exacerbate four (see Supplementary Data 8). TB-associated peripheral blood transcripts significantly affected by each treatment were those determined to have a Benjamini–Hochberg false discovery rate (FDR) adjusted p-value less than 0.05 for the variable Group in the limma/voom model (see “Methods” section). D, E Effect of HRZE and NTZ on blood gene expression for a set of IBD-associated genes from Palmer, et al.. Both HRZE (D) and NTZ (E) cause different genes to either renormalize (HRZE 66, NTZ 34) (see Supplementary Data 9) or exacerbate (HRZE 55, NTZ 21) (see Supplementary Data 10). IBD-associated peripheral blood transcripts significantly affected by each treatment were those determined to have a Benjamini–Hochberg false discovery rate (FDR) adjusted p-value less than 0.05 for the variable Group in the limma/voom model (see “Methods” section) P value in y axis is adjusted according to Benjamini–Hochberg (FDR). Source data are provided as a Source Data file.

Fig. 5. Longitudinal profiling of HRZE treatment induced changes of microbiome composition and peripheral gene inflammatory expression.

A. Schematic diagram of Cohort. B. Time to positivity was measured at baseline (n = 19 biologically independent samples), day 14 (n = 14 biologically independent samples), one month (n = 5 biologically independent samples), and two months (n = 9 biologically independent samples). To determine statistical significance of differences in TTP at different time points we fit the linear mixed-effect model TTP~Sex+Age+Time+1∣ID. We inspected the p-value associated by running contrasts for the variable Time to determine significant (p-value for the contrast <0.05) differences in TTP. C Microbiome diversity was computed for each study volunteer at baseline (n = 20 biologically independent samples), day 7 (n = 7 biologically independent samples), day 14 (n = 19 biologically independent samples), one month (n = 13 biologically independent samples), two months (n = 13 biologically independent samples), and 6 months (n = 13 biologically independent samples). Microbiome α diversity was measured using the Inverse Simpson index. To determine statistical significance of differences in α diversity at different time points we fit the linear mixed-effect model Diversity~Sex+Age+Time+1∣ID. We inspected the p-value associated by running contrasts for the variable Time to determine significant (p-value for the contrast <0.05) differences in microbiota diversity D Volcano plots showing significance of differences in microbiome composition vs. fold change from baseline at Day 7, Day 14, Day 30, Day 56, and Day 180. E Normalized enrichment scores calculated for the Hallmark Pathway list for Day 14 vs. Baseline, Day 56 vs Day 14, and Day 56 vs. Baseline. F–H Volcano plot showing TB transcripts from ref. at Day 14 vs. Baseline at Day 56 vs. Day 14, and at Day 14 vs. Baseline. I–K Volcano plot showing transcripts from Palmer et al. of IBD cases vs. controls detected in this study for Day 14 vs. Baseline, Day 56 vs. Day 14, and Day 56 vs. Baseline. Source data are provided as a Source Data file.

Fig. 6. Use of Random Forest Regression Modeling to search for associations between immune-related peripheral blood gene signatures and changes in gastrointestinal microbiota and TTP.

The heatmap displays the sign of the derivative of the ALE curve (See Text). Blue/orange entries indicate features found to significantly associate with changes in a specific inflammatory pathway. Blue indicates a negative relationship, while orange a positive. For each immune pathway, a pathway-ASV association was determined significant if the Benjamini–Hochberg false discovery rate (FDR) adjusted p-value from the permutated importance analysis was found to be less than 0.05. Black symbols are used to identify the modeling-identified top important predictor (i.e., the predictor that if missing would lead to highest increase in mean squared error between model predictions and observations) for each specific host pathway. This analysis shows that reduction in TB burden and increased abundance of health-associated Cluster IV and XIVa Clostridia predicts inflammatory dampening. In contrast, increased abundance of oxygen-tolerant pathobionts including Enterococcus, Streptococcus, and E. coli predicts inflammatory exacerbation. A table reporting variable importance values, slope and intercepts from the ALE plot calculation and the related p values is provided in Supplementary Data 18. Source data are provided as a Source Data file.

Fig. 7. Analysis of microbiome and blood peripheral gene expression in an independent healthy control human cohort validates association between specific microbiome members and host peripheral gene expression.

A NES scores of 50 Hallmark pathways from the MiSigDB on a per-sample basis for all cohorts in this study. NES score was calculated using the variance stabilized transformed counts from DSEeq, calculated with the GSVA package in R, and plotted after scaling (Z score) across all samples. Columns are split based on arm or group membership and rows are split based on Hallmark pathway categorization. B Random forest regression results associating specific microbial taxa with Hallmark pathways. Only pathways identified in the RFR model are shown. The ‘Relation’, calculated by taking the first derivative of the ALE plot for each relationship, is positive if the pathway positively associates with a particular ASV, or negative if the pathway negatively associates with a particular ASV. Source data are provided as a Source Data file.