Peripheral inflammatory response in human tuberculosis treatment is predicted by a combination of pathogen sterilization and microbiome dysbiosis

Gut microbiome diversity is depleted after two weeks of HRZE or NTZ treatment

As detailed elsewhere, the GHESKIO centers in Port au Prince, Haiti conducted a prospective, randomized, early bactericidal activity (EBA) study in treatment-naive, drug-susceptible adult patients with uncomplicated pulmonary tuberculosis (TB) (ClinicalTrials.gov Identifier: NCT02684240)²³. Participants were randomized to receive either NTZ, 1000 mg po (oral) twice daily, or standard oral therapy with isoniazid 300 mg daily, rifampin 600 mg daily, pyrazinamide 25 mg/kg daily, and ethambutol 15 mg/kg daily (referred to as HRZE) for 14 days (Figure 1A). The primary endpoint of the trial was sputum bacterial load (measured by time to culture positivity, TTP) in a BACTEC liquid culture system. Sputum was collected from 6pm to 9am every other day to quantify mycobactericidal activity of each treatment regimen. As reported²³, HRZE resulted in a predictable increase in the TTP (corresponding to reduced bacterial load) over the first two weeks of therapy compared to pretreatment TTP, consistent with its known potent bactericidal activity. However, NTZ, despite potent in vitro activity²⁴, did not have any significant effect on TTP after 14 days (Figure 1B)²³. This lack of NTZ efficacy was traced to a failure of the drug to penetrate the sputum²³. All patients were subsequently switched to HRZE standard treatment.

• Download figure
• Open in new tab

Figure 1: Both HRZE and NTZ perturb the gut microbiome after two weeks of therapy.

A. Schematic showing active TB patients randomized to either HRZE (standard of care) or NTZ. B. Paired time to positivity (TTP) at day 0 and day 14 for the NTZ treatment cohort and HRZE treatment cohort. C. Principal components 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 grey line connects pre and day 14 treatment paired samples. D. Microbiota alpha diversity plotted using the Inverse Simpson index. There was no significant difference between the pretreatment groups, and both groups had significantly (p < 0.01, Wilcoxon signed-rank test) reduced alpha diversity after 14 days of treatment.

We have reported¹⁹ that HRZE therapy depletes members of the order Clostridiales, but the cross-sectional design of that study did not allow for conclusions about the rapidity of this effect and most importantly, did not include pretreatment samples that will allow for assessment of baseline microbiome composition. To investigate microbiome changes induced by NTZ or HRZE, we extracted and amplified bacterial and archaeal DNA using V4 – V5 16S rDNA sequencing. Stool samples were collected before the start of treatment and on day 14 of therapy (Figure 1A). Using principal components analysis (PCoA) with Bray Curtis distances, we determined that antibiotic administration induced rapid changes in microbiome community structure after two weeks in both the NTZ and HRZE groups (PERMANOVA: p<0.001 and p<0.02, respectively), compared with pretreatment (Figure 1C). HRZE samples clustered closer to pretreatment samples than did NTZ to pretreatment. Both treatments were characterized by a significant drop in alpha diversity (the Inverse Simpson index) when compared to pretreatment samples, and despite expected interindividual heterogeneity, NTZ treatment samples had a significantly lower diversity compared to HRZE (Figure 1C). There was no significant difference between the day 0 diversity metrics between the pre-randomized active TB groups (Figure 1D).

Overlapping taxonomic alterations in microbiome composition induced by NTZ and HRZE

Taxonomic profiling of differential amplicon sequence variants (ASVs) obtained from the 16S rDNA profiling between the pretreatment and antibiotic groups revealed that, whereas all trial participants began with typical heterogeneous community structure, after two weeks of treatment with HRZE, members of the class Clostridia were depleted, leaving most other clades unaltered (Figure 2A,B). In contrast, NTZ had a much more pronounced and broad effect compared to HRZE (Figure 2A,C). Specifically, NTZ not only reduced the relative abundance of a greater number of Clostridia compared to HRZE, but also led to the increase of aerobic and facultative anaerobic pathobiont organisms such as Escherichia and Klebsiella (Proteobacteria). Additionally, several participants became dominated (with domination defined as >30% relative abundance according to Taur, et al²⁵) by single clades of Actinobacteria or Bacilli in the NTZ arm of the trial (Figure 2A,C and Supplementary Figure S1 and S2).

• Download figure
• Open in new tab

Figure 2: Overlapping and distinct microbiome perturbation induced by NTZ and HRZE.

A. Relative abundance showing microbiota composition in each individual patient as stacked paired samples corresponding to pre and post treatment. Only paired samples were used for these analyses. B. Unsupervised hierarchical clustering of the abundances of ASVs identified to be significantly affected by HRZE treatment (p<0.01); VST indicates variance stabilized transformed counts from DESeq2. C. Unsupervised hierarchical clustering of the abundances of ASVs identified to be significantly affected by NTZ treatment (p<0.01). D. Number of Clostridia genera depleted (i.e. reduced in relative abundance) by HRZE and NTZ treatment.

Comparison of the overlap of the two antibiotic perturbations with respect to number of Clostridia genera affected (i.e., testing for differences pre-post HRZE/NTZ administration on bacterial abundances grouped at the genus level) revealed that 83% (10/12) of the HRZE-depleted Clostridia were also depleted by NTZ (Figure 2D). HRZE uniquely affected Clostridia genera Acutalibacter and Hespellia, whereas NTZ affected 22 unique Clostridia genera, including Blautia, Dorea, Eubacterium, Faecalibacterium, Oscillibacter, Ruminococcus, and Lachnospiraceae. Both treatments affected the genera Clostridium, Fusicatenibacter, Hungatella, Intestinibacter, Intestinimonas, Kineothrix, Roseburia, Sporobacter, Terrisporobacter, and an uncultured member of the Family Ruminococcaceae. Taken together, these data demonstrate that NTZ had a more severe disruptive effect on the intestinal microbiota than HRZE and that most of the HRZE effects on the microbiota (e.g. loss of Class Clostridia) were also evident in NTZ treated subjects.

HRZE and NTZ uniquely affect host peripheral gene expression

To determine how inflammatory transcriptional signatures in peripheral blood are affected by TB treatment, we applied unsupervised (PCA) and supervised learning (Support Vector Machine and Random Forest Classification) methods to the DESeq-normalized RNA transcript abundances from RNA derived from peripheral blood (Supplementary Figure S3). This analysis showed that, overall, HRZE treatment caused a more substantial change (in terms of the number of genes affected) in expression profile than NTZ when comparing pretreatment to 2 week treatment transcriptome samples, as measured by random forest/support vector machine ROC curves or principal component analysis (Supplementary Figure S3A-B).

To identify specific RNA transcripts affected by each treatment, we used DESeq analysis on the paired-sample transcript abundance of RNA and compared pretreatment samples to samples obtained at day 14 after HRZE (n=8) or NTZ treatment (n=14) independently. The paired nature of samples was considered in the analysis to account for person-specific baseline normalization. We found that 2 weeks of HRZE treatment was associated with changes in 1374 transcripts at a p-value cutoff of p<0.05 and 503 at a p-value cutoff of p<0.001 (throughout, all p-values from DESeq are adjusted using the Benjamini-Hochberg method for multiple comparisons). Repeating the same analysis for NTZ-treated individuals, we identified 811 differentially expressed genes at a p-value cutoff of p<0.05 and only 15 at p<0.001.

We determined the functional pattern of treatment-induced changes in overall transcript abundance, by performing gene set enrichment analysis (GSEA)²⁶ (See Methods) using the ranked DESeq gene expression data for the differentially expressed genes in each arm separately. In the HRZE arm we observed reduced expression of the pathways of inflammatory response, IFNγ response, IFNα response, TNFα signaling via NFκB, and IL6 JAK STAT3, all of which are consistent with the immunologic effects of antibiotic induced reduction in the levels of the bacterial pathogen, given the demonstrated relevance of these signaling pathways to pathogenesis^12,15,27 (Figure 3A). In contrast, NTZ treatment, which perturbed the microbiome without significantly affecting Mtb bacterial load in the sputum, had the opposite effect. Inflammatory signaling pathways reduced by HRZE, including TNFα signaling, IFNγ signaling, and type 1 interferon signaling, were all enriched by NTZ treatment (Figure 3A). Several other pathways such as hypoxia, apoptosis, and reactive oxygen species (ROS), which are considered hallmarks of immune dysregulation²⁸, were also enriched by NTZ treatment (Figure 3B). When comparing the list of transcripts significantly affected by either HRZE (1374) or NTZ (811), after taking into account person-specific gene normalization, we found that only 86 genes were affected by both treatments (p < 0.05). Eighty of these genes had the same pre/post treatment pattern in HRZE and NTZ-treated individuals, thus leaving 1294 (94.2%) and 731 (90.1%) to be HRZE and NTZ-specific changes to peripheral gene expression (Figure 3C). Of the 86 genes affected by both HRZE and NTZ, only 10 were reported to be peripheral blood transcriptomic markers of active tuberculosis ^12,13,15. Interestingly, six active-TB associated genes belong to the group of 80 genes that similarly respond to HRZE and NTZ (7.5%). The remaining four TB-related genes (P2RY14, ADM, CARD16, DHRS9) belong to the group of six genes (P2RY14, NRN1, ADM, JAG1, CARD16, DHRS9) that were reduced in expression and renormalized with HRZE but increased and exacerbated TB signature with NTZ (67%) (Figure 3C).

• Download figure
• Open in new tab

Figure 3: Hallmark pathway gene set enrichment analysis and gene expression comparison in HRZE and NTZ treated cohorts.

A,B. Hallmark gene pathway changes associated with 2 weeks of HRZE (A) or NTZ (B). Positive are pathways overrepresented at 2 weeks of therapy, and negative are pathways underrepresented at 2 weeks, both compared to pretreatment. All pathways are significant (p<0.05) with the associated Normalized Enrichment Score (NES) shown on the × axis, which considers pathway size. C. Overlap in genes that are differentially altered in both HRZE and NTZ treatments compared to pretreatment. Transcripts annotated in the left column in red (Up) rise with HRZE and NTZ. Transcripts annotated in blue (Down) are suppressed by both HRZE and NTZ. Genes annotated in green are down in HRZE and up in NTZ. The second column of vertical metadata scores each gene as either present or absent in prior active TB transcriptional signatures.

To focus our analysis on validated transcriptomic markers of active TB from prior studies, we examined the recently reported list of 373 transcripts that have been associated and validated in multiple human cohorts on multiple sequencing platforms (microarray and RNAseq) to be differentially abundant between active and healthy control and/or active and LTBI individuals¹⁷ in a comprehensive meta-analysis. In our study, we detected 361 of these 373 transcripts in pretreatment active TB subjects. We defined three classes of changes to these transcripts with two weeks of HRZE or NTZ treatment: 1) renormalization (transcripts whose pre-post HRZE/NTZ foldchange in expression displays the same sign (or direction) of the previously-reported fold-change between active TB and control/LTBI; 2) unchanged (transcripts with no change in expression between pre-post HRZE/NTZ administration); and 3) exacerbation (genes whose pre-post HRZE/NTZ fold-change sign is opposite to the previously-reported fold-change between active TB and control/LTBI). Of the 361 transcripts, 173 (48%) were found to be affected by HRZE. Of these 173, 151 (87%) renormalize with HRZE treatment, whereas 22 (13%) exacerbate (Supplementary Table S1, Figure 4). In contrast, NTZ was found not only to have a smaller overall effect on the active TB signature, but also to only contribute to exacerbation. Specifically, only 28 genes were affected by NTZ (8%), of which 26 (96%) were in the exacerbation category (Supplementary Table S1, Figure 4).

• Download figure
• Open in new tab

Figure 4: Effect of HRZE and NTZ on active TB transcriptional signatures.

A. Summary of the number of active TB genes that renormalize (blue), do not change (grey), or exacerbate (red) as a result of each drug administration. B. Comparison of log2 Fold Change in gene expression from¹⁷ where the 373 active TB signature was first introduced with the log2 Fold Change of the same genes pre-post HRZE (left) or NTZ (right) from this study. A significant positive correlation (Pearson’s p<0.01) is observed for the 151 genes that renormalize with HRZE treatment.

Taken together, our data suggest that both HRZE and NTZ lead to changes in the peripheral transcriptomic pattern of active TB, with HRZE having a more dramatic effect, at least as measured by the number of genes differentially expressed before and after two weeks of therapy compared to NTZ, both overall as well as when considering a previously-reported signature of active TB. Additionally, our analysis suggests that two weeks of HRZE treatment normalizes almost 50% of the peripheral gene expression profile characteristic of active tuberculosis, consistent with prior studies, whereas NTZ had the opposite effect: exacerbating the transcriptional pattern for a subset of transcripts that reflect disease activity.

Relationships between gene expression and changes in the microbiome in the longitudinal treatment cohort

Resolution of disease manifestations by antibiotic therapy of a chronic, pro-inflammatory infection such as tuberculosis is likely to reflect a complex interplay of pathogen killing and resolution of inflammatory responses, which may occur with different temporal profiles²⁹. Although a primary effect of antibiotics is killing of the pathogen, perturbation of the microbiome during therapy might also affect disease resolution due to an independent effect on peripheral inflammation, but these two effects are difficult to uncouple and therefore remain hypothetical. Because we observed changes in peripheral transcriptomic profiles by both HRZE and NTZ in active TB subjects, in some cases divergent between HRZE and NTZ, despite the fact that NTZ did not significantly reduce Mtb bacterial load in the sputum, this dataset provided a unique opportunity to investigate potential effects of antibiotic induced microbiome changes on disease driven inflammation. We hypothesized that TB therapy-induced perturbations of the gut microbiota that might be predictive of the observed changes in host gene expression.

To decouple the effect of microbiome perturbation from Mtb-killing in peripheral inflammatory transcriptomics, we used a machine learning approach. Specifically, we built random forest regression (RFR) models³⁰ to predict the DESeq-normalized pre- and post-treatment expression profiles of each significantly differentially abundant transcript as a function of both Mtb bacterial load (TTP) (at day 0 and day 14) and of DESeq-normalized abundance of ASVs (microbiome components) (also at day 0 and day 14) found to be affected by each TB treatment, separately. We modeled data from each cohort independently because the treatment effects at the single gene and ASV level were distinct between each cohort and because at baseline there is no significant difference in microbiome and gene expression due to arm membership (PERMANOVA p > 0.05). We used the Boruta algorithm ³¹ for feature selection in predicting the expression profile of each gene (see Methods). We reasoned that this approach was appropriately suited for this type of “large p, small n” multi-omics dataset common in clinical research³². Several advantages of RFR modeling include: being agnostic to model structure (e.g. non-parametric regression), not having to meet common assumptions underlying classical regression techniques, and being able to intrinsically perform ranked feature selection. Importantly, while the interpretation of RFR is apparently less immediate compared to traditional regression (e.g. there are per-se no regression coefficients or betas), downstream analysis, which includes Permutated Importance ³³ and accumulated local effects calculations³⁴ (see Methods) allows for the estimation of the significance of predictors (e.g. TTP, microbiome constituents, etc.) as well as their effects on the dependent variable (e.g. host transcriptomic markers).

For both sets of HRZE- and NTZ-affected transcripts, we estimated the frequency of each feature to be identified as important by the RFR model. We performed this analysis considering all differentially expressed genes in both datasets (Supplementary Figure S4) as well as focusing only the 361 genes belonging to the list of 373 transcripts previously reported as an active TB signature. Remarkably, this analysis identified bacterial load and the abundance of members of Clostridia as the most important predictors of gene expression (Supplementary Figure S6A and 5B). For the HRZE arm (Supplementary Figure S6A), change in bacterial load (TTP) was the top predictor of the change in transcript abundance with treatment, which is found to be important in predicting the change in 75% of the 173 active TB transcripts and HRZE affected transcripts (Supplementary Figure S6). This is consistent with the idea that HRZE normalizes the active TB signature by reducing the pathogen burden in the lung (Supplementary Figure S6C). Surprisingly, even though sputum bacterial load did not change in aggregate in the NTZ cohort, TTP is still found to be the second most important feature in predicting changes in active TB gene expression, as it is found to contribute to the abundance of 48% of the 27 active TB and NTZ affected transcripts (Figure 3). We confirmed that the relationship between bacterial load and transcript abundance in the NTZ cohort (obtained via RFR) was a biological signal by performing linear-mixed effect modeling (lme in R) to predict the expression of each of these NTZ associated transcripts as a function of TTP and using patient ID as a random effect (Supplementary Figure S4) (model: lme(∼ TTP, random = ∼1|Patient.ID)). We found that despite the lack of clinical efficacy of the NTZ drug as measured by average TTP in the treated group, transcripts that are known to respond to TB bacterial load are altered in a direction that correlates with the TTP signal. This reflects the fact that while the magnitude of TTP change was not clinically significant across the entire NTZ cohort, person-specific changes in TTP were present in the NTZ cohort for TB-associated transcripts.

With respect to microbiome effects, Clostridia are predicted by both models to be the most important microbiome component predicting the change in the inflammatory transcriptomic pattern that reflects active TB and normalizes with treatment. Specifically, in the HRZE arm (Supplementary Figure S6A) E. hallii (ASV 336) and S. termiditis (ASV 71) are the second and third-best predictors of the abundance of over 50% of the HRZE-affected transcripts associated with active TB. In the NTZ arm (Supplementary Figure S6B) G. formicilis (ASV 27) is identified as the most important predictor with over 50% of NTZ-affected active TB transcripts predicted to co-vary with it. Of the other important microbiome features selected by the models, we observed the presence of a number of well characterized short chain fatty acid-producing Clostridia (including F. prausnitizii, Clostrdium spp., B. lutii in the HRZE group, or D. longicatena and B. lutii in the NTZ group), which are all depleted by the anti-TB treatment (Supplementary Figure S6C-D). Interestingly both models identified the same ASV, Clostridium XIVa F. saccharivorans (ASV 57) to be an important predictor of transcriptomic markers that are depressed by both treatments. Peculiar to the NTZ-data based model is the identification as important variables of S. alactolyticus (ASV 40), E. faecium (ASV 222) and K. pneumoniae (ASV 59) which are observed to significantly increased in abundance and dominate several of the NTZ-treated individuals (Figure 2).

Variable importance estimation techniques such as Boruta allow determination of the presence of significant associations among variables, but not their directionality (positive/negative effect). Therefore, we determined the direction of the contribution of each important feature to the abundance of each of the modeled active TB signature transcripts from the RFR modeling, by adapting the concept of accumulated local effects (ALE) calculations (see Methods)³⁴ and using the sign of the average of the first-order derivative of the ALE plot as a binary indicator for direction. As most of these plots are monotonic and continuous curves, the first order derivative is a good approximation of the quantitative contribution of a hypothesized biologically relevant predictor to transcript abundance, in this case microbiota constituents or change in Mtb bacterial load (TTP). We displayed the sign of the features x genes matrices containing the first order derivative values (or slopes) as clustered heatmaps (Figure 5A-B). For both the HRZE and NTZ dataset-derived models, we observed that sputum bacterial load (TTP) clustered independently of the model-selected microbiome features (Figure 5A-B). The effect of reducing sputum bacterial load was predominantly to normalize gene expression in the relevant inflammatory pathways such as IFNγ and IFNα (negative contribution in heatmap for TTP means that rising TTP (i.e., lower bacterial load) is associated with normalizing transcripts) (Figure 5A).

• Download figure
• Open in new tab

Figure 5: Hallmark Pathway Analysis of peripheral gene transcriptional data and its relationship to gastrointestinal microbiota and M. tuberculosis bacterial load (TTP).

We used Random Forest regression to model the change in gene expression as a function of gastrointestinal microbiota and reduction in bacterial load (TTP). Top predictors are shown for both HRZE (A) and NTZ (B) study arms. The contribution of each predictor to a gene’s expression profile is estimated by calculating the derivative (or slope) of the accumulated local effects (ALE) plots between each predictor and each modeled gene (see Methods). The sign of the derivative is displayed in the heatmaps. On the right, the presence or absence of a particular gene in a Hallmark gene set is indicated. If that gene is in the pathway, it is marked as red. Some genes are in multiple pathways, and only the significantly enriched pathways from the differential analysis (of the full transcriptome) are shown. The Normalized Enrichment Score (NES) is how significantly overrepresented a particular pathway is in the treated case (comparing pretreatment vs treatment in both instances) and takes into account pathway size. A majority of the organisms are obligate anaerobic Clostridia, from the following clusters: Clostridium disporicum (Cluster I), Faecalibacterium prausnitzii (Cluster IV), Eubacterium hallii (Cluster XIVa), Sporobacter termitidis (Cluster XI), Fusicatenibacter saccharivorans (Cluster XIVa), Clostridium fimetarium (obligate anaerobe), Clostridium innocuum (obligate anaerobe), Gemmiger formicilis (obligate anaerobe), Blautia luti (Cluster XIVa), Dorea longicatena (obligate anaerobe), Hungatella hathewayi (obligate anaerobe).

In the model trained on the NTZ treatment data, we find two opposing effects of NTZ induced microbiome perturbation. Depletion of Clostridia by NTZ is predicted to contribute to the exacerbation of inflammatory pathways of active TB observed in the transcriptomic data (Figure 5B: negative contribution in heat map=lower Clostridia>>higher inflammatory transcripts) consistent with the anti-inflammatory properties of these microbiota members. In contrast, the model predicts that pathobionts, whose abundance is enhanced by NTZ (E. faecium, K. pneumoniae and S. alactolyticus, Figure 2C), have the opposite effect and exacerbate inflammatory transcripts (Figure 5B: positive contribution in heat map=higher pathobionts>>higher inflammatory transcripts). Pathobionts are potentially pathogenic symbionts of the microbiota, and include E. coli, K. pneumoniae and S. alactolyticus).

Taken together, our data and related computational analyses show that the changes in inflammatory gene expression that accompany treatment of TB may be mediated both by the anti-microbial activity of the drugs that lead to pathogen clearance and by antibiotic induced changes in microbiome composition. Our model identifies two modules of microbiome-inflammatory effects. The first is the exacerbation of TB associated inflammation by depletion of Clostridia, which is evident in both the HRZE and NTZ models. Additionally, the enhancement of pathobionts such as Klebsilla, which only occurs with NTZ, also exacerbates inflammatory pathways. Based on this modelling, we predict that successful disease resolution will be associated with preservation of Clostridia, whereas their depletion and consequent enhancement of Proteobacteria and Bacilli pathobionts might slow resolution or even support inflammatory exacerbation.

Relationship of the microbiome and peripheral gene expression in a healthy control validation cohort

The results from our primary analysis in the longitudinal treatment cohort demonstrate that specific gastrointestinal bacteria are associated with proinflammatory peripheral blood gene signatures in humans. Specifically, we predict that higher abundance of Clostridia is negatively associated with inflammation (e.g. INFα, INFγ, IL6/JAK/STAT3, Inflammatory Response gene signatures) while high abundance of commonly known pathobionts including E. coli, Klebsiella, Citrobacter, Streptococcus and Enteroccocus promotes exacerbation of these signatures. To determine the generality of these results we analyzed a control set of human data from two healthy cohorts. A subset of these data were previously reported in our work⁵, and come from a cross-sectional study of healthy household contacts of active pulmonary TB patients (termed Family Contacts, FC) and healthy unexposed donors from the same community in Haiti (termed Community Controls, CC) (see Methods). For these two new cohorts we have a total of 52 healthy control individuals (18 FC and 36 CC) for which we gathered both microbiome 16S rRNA sequencing data as well as peripheral blood transcriptomics for the same individuals.

We analyzed the blood transcriptomics and performed correspondence analysis on the DSEeq2 normalized abundance of the transcripts in each sample (Figure 6A). We observed pretreatment and treated HRZE and NTZ sample clustering away from FC and CC samples, highlighting the fact that individuals with active TB, even while on treatment, represent the transcriptional profile of active TB patients has not returned to a profile of healthy controls (Supplementary Figure S7). More importantly we see that FC and CC samples do not separate in this ordination, suggesting that both cohorts can be used as healthy controls in the subsequent analysis.

• Download figure
• Open in new tab

Figure 6: 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. t-Distributed Stochastic Neighbor Embedding (t-SNE) ordination all RNAseq samples used in this analysis. This highlights the transcriptional differences between active TB patients, even those on HRZE treatment, and the healthy control cohorts. B. Within-sample GSEA (ssGSEA) 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, and plotted after scaling across all samples. C. Four representative pathways that correlated with microbiome composition in all cohorts analyzed in this study. D. Pearson’s correlations of ASVs with Hallmark pathways in an unbiased fashion highlighting negative correlation of many Clostridia and the proinflammatory pathways highlighted in 6C.

To link transcript abundance to immune pathway enrichment we utilized gene set variation analysis (GSVA) (see Methods)³⁵. GSVA enables computing an enrichment score for any defined list of genes (taken as a surrogate for a biological pathway) in each sample. For each sample, each pathway generates a normalized enrichment score (NES) that can be used for downstream analysis to compare pathway profiling across samples. Importantly, this method is completely agnostic to the phenotype label of the sample.

We performed unsupervised clustering on the samples-by-pathway NES scores, and found that individuals from FC and CC cluster separately from individuals in the longitudinal treatment cohort (Figure 6B). We observed that the NES score signature comparisons between the FC and CC individuals, and individuals before treatment and after two-weeks of HRZE or NTZ have qualitatively reduced enrichment in INFα, INFγ, IL6/JAK/STAT3, and Inflammatory Response (Figure 6B). As discussed, these are the pathways that renormalize with HRZE treatment and exacerbate with NTZ in the longitudinal dataset (Figure 3). These are also the pathways negatively associated with Clostridia abundance and positively with Proteobacteria and Bacilli abundance in the random forest modeling analysis (Figure 5). Interestingly, post HRZE treatment individuals have an enrichment distribution for this pathway which is in between FC and CC individuals and NTZ/Pretreatment samples thus reinforcing the result that HRZE treatment promotes renormalization (Figure 6C).

Even though FC and CC individuals show an overall lower inflammatory enrichment compared to active TB samples, and samples post HRZE and NTZ treatment (Figure 6C), we hypothesize that the FC and CC intragroup inflammatory pathway enrichment distributions would correlate with the abundance of microbes that we found to predict the expression of genes associated with these pathways in the Random Forest Modeling analysis performed on the HRZE and NTZ microbiome/transcriptomic data (Figure 5 and Figure S6). We performed Pearson’s correlation analysis between these NES values, and the abundance of ASVs from bacterial orders having representative species found to be predictive of the abundance of these inflammatory pathways in the longitudinal treatment cohorts (Figure 6D). In these two human cohorts, we again observe a negative association between Clostridia ASVs and INFα, INFγ, IL6/JAK/STAT3, and Inflammatory Response pathways, as well as a positive association between these pathways and the abundance of pathobionts including Enterococcus, Streptococcus, E. coli and Klebsiella. We interpret these correlations as independent validation of the observation that specific microbiome members are strongly correlated with peripheral inflammatory response pathways in both inflammatory disease states and in homeostatic conditions.

Ethical statement and study approval

All volunteers provided written informed consent to participate in this study. All human studies were reviewed and approved by the IRBs of both Weill Cornell Medicine and Groupe Haitien d’etude du Sarcome de Kaposi et des Infections Opportunistes (GHESKIO) Centers (Port-au-Prince, Haiti). Participants provided informed consent prior to peripheral blood draw for whole blood collection and stool collection for 16S rDNA sequencing. All methods and procedures were performed in accordance with the relevant institutional guidelines and regulations.

Donor recruitment and protection of human subjects

Microbial DNA extraction from stool

DNA extraction from stool was performed as described.¹⁹ Stool specimens were collected and stored for less than 24 hours at 4°C, aliquoted (∼2 ml each), frozen at –80°C, and shipped to WCM/MSKCC. About 200 – 500 mg of stool from frozen samples was suspended in 500 μl of extraction buffer (200 mM Tris-HCl [Thermo Fisher Scientific], pH 8.0; 200 mM NaCl [Thermo Fisher Scientific]; 20 mM EDTA [MilliporeSigma]), 210 μl of 20% SDS, 500 μl of phenol/chloroform/isoamyl alcohol (25:24:1; MilliporeSigma), and 500 μl of 0.1-mm–diameter zirconia/silica beads (Biospec Products). Samples were lysed via mechanical disruption with a bead beater (Biospec Products for 2 minutes, followed by 2 extractions with phenol/chloroform/isoamyl alcohol [25:24:1]). DNA was precipitated with ethanol and sodium acetate at –80°C for at least 1 hour, resuspended in 200 μl of nuclease-free water, and further purified with QIAamp DNA Mini Kit (Qiagen) according to the manufacturer’s protocols. DNA was eluted in 200 μl of nuclease-free water and stored at –20°C.

16S rDNA sequencing

Primers used to amplify rDNA were: 563F (59-nnnnnnnn-NNNNNNNNNNNN-AYTGGGYDTAAAGN G-39) and 926R (59-nnnnnnnn-NNNNNNNNNNNN-CCGTCAATTYHTTTR AGT-39). Each reaction contained 50 ng of purified DNA, 0.2 mM dNTPs, 1.5 μM MgCl2, 1.25 U Platinum TaqDNA polymerase, 2.5 μl of 10× PCR buffer, and 0.2 μM of each primer. A unique 12-base Golay barcode (Ns) preceded the primers for sample identification after pooling amplicons. One to 8 additional nucleotides were added before the barcode to offset the sequencing of the primers. Cycling conditions were the following: 94°C for 3 minutes, followed by 27 cycles of 94°C for 50 seconds, 51°C for 30 seconds, and 72°C for 1 minute, where the final elongation step was performed at 72°C for 5 minutes. Replicate PCRs were combined and were subsequently purified using the Qiaquick PCR Purification Kit (Qiagen) and Qiagen MinElute PCR Purification Kit. PCR products were quantified and pooled at equimolar amounts before Illumina barcodes and adaptors were ligated on using the Illumina TruSeq Sample Preparation procedure. The completed library was sequenced on an Illumina Miseq platform per the Illumina recommended protocol.

16S rDNA bioinformatics analysis

Forward and reverse 16S MiSeq-generated amplicon sequencing reads were dereplicated and sequences were inferred using dada2.⁴⁶ Potentially chimeric sequences were removed using consensus-based methods. Taxonomic assignments were made using BLASTN against the NCBI refseq rna database. These files were imported into R and merged with a metadata file into a single Phyloseq object.

Deposition of data. 16S rDNA sequencing data is deposited with the SRA under accession no. PRJNA445968 (https://www.ncbi.nlm.nih.gov/bioproject/PRJNA445968). Code used for 16S analysis is available at https://wipperman.github.io/TBRU/.

Computational Analysis

Raw ASV count data was normalized using DESeq2⁴⁷. DESeq2 was used to make differential abundance and expression comparisons between treatment cohorts and between individual pre and post-treatment samples for both microbiota and host gene count data.

To determine how the anti-TB treatment affects both microbiome and peripheral gene expression profiles we performed differential analysis on the counts data obtained by microbiome DNA and peripheral blood RNA sequencing. As the primary endpoint of the clinical trial was powered to determine differences in Mtb load (TTP), we had to determine the statistical power available to identify significant differences in the abundance of both microbiota members and in the expression of peripheral genes. We ran power calculations (see Methods) to determine that with 16 microbiome samples and 8 RNAseq samples for the HRZE cohort, with 80% power at α<0.05, we could detect a fold change of 1.4 for microbiome difference and a fold change of 1.8 for mRNA transcripts. In the NTZ cohort, with 18 microbiome samples and 14 RNAseq samples, with 80% power at α<0.05, we can detect a fold change of 1.4 for microbiome differences and a fold change of 1.6 for mRNA transcripts. Additionally, we utilized DESeq in our analyses using baseline normalization within an individual as the gold standard. Power calculations were performed with the RNAseqPower package in R. For microbiome data we calculated a biological coefficient of variation of 0.3, and for RNAseq, we used a coefficient of variation of 0.4. We estimated the expected minimum fold change that we could observe for each group based on the sample size, sequencing depth, and an α < 0.05.

To test whether there were differences between groups, we employed PERMANOVA using the Adonis function in the Vegan R package (https://cran.r-project.org/web/packages/vegan/vegan.pdf), which partitions a distance matrix of ASV count data and runs one-way ANOVA between groups of samples.

To predict the post-pre fold change in abundance of HRZE and NTZ-affected host genes as a function of the corresponding fold change in abundance of microbiota ASVs and TTP using Random Forest Regression.³⁰ To perform feature selection for each gene and rank features based on their prediction importance we used Boruta⁴⁸. Boruta is a RF classification and regression wrapper for feature selection that allows identification of variables important for the prediction task while also removing redundant ones. Boruta creates a copy of each independent variable and shuffles them to remove correlation with the original variables (shadow variables). Using this augmented set Boruta builds a RF model and performs a Variable Importance Estimation of all the independent variables (both original and shadow). For every variable it then computes a normalized accuracy score. A true variable is important if its normalized accuracy score is significantly greater than that of shadow variables.

Accumulated Local Effect (ALE) plots were implemented to compute how the change in the expression of each modeled host gene is affected by the level of each predictor identified as important by the forest model. Code to analyze the data and to reproduce all the figures and results is available on Github at https://wipperman.github.io/TBRU/.

Within sample GSEA analysis

The ssGSEA (single sample gene set enrichment analysis) method⁴⁹ was used to profile within-sample differences between pathways from the MiSigDB Hallmark pathways list⁵⁰ with the GSVA package in R³⁵. The MiSigDB Hallmark pathways list is a well validated set of general curated biological pathways that can give insight into specific biological and cellular processes. Additionally, we obtained a list of well validated active TB signatures from the TBSignatureProfilier R package (David Jenkins, Yue Zhao, W. Evan Johnson and Aubrey Odom (2020). TBSignatureProfiler: Profile RNA-Seq Data Using TB Pathway Signatures. R package version 0.0.0.9005. https://github.com/compbiomed/TBSignatureProfiler. Variance stabilized transformed (vst) counts derived from DESeq2 were used as input into the GSVA function in the GSVA R package with default parameters and scaled Normalized Enrichment Scores (NES) were plotted as heatmaps. Importantly, unlike classical GSEA, this analysis is agnostic to sample phenotype.