Intestinal microbiome composition influences the peripheral 1 inflammatory state during treatment of human tuberculosis 2

Although the composition of the intestinal microbiota influences systemic immune responses, the contribution 41 of this relationship to infectious disease pathogenesis and to the resolution of infectious diseases by antibiotic 42 therapy is poorly understood. This question is rarely examined in humans due to the difficulty in dissociating 43 the immunologic effects of antibiotic-induced pathogen clearance and salutary microbiome alteration. To 44 address these questions in the context of Tuberculosis, we analyzed three independent human datasets from 45 Haiti (two longitudinal treatment and one cross-sectional control), a prototypical infection remarkable for 46 chronic inflammation, in which we had measured sputum TB bacterial load, gut microbiota composition, and 47 peripheral blood transcriptomics. Using data from the longitudinal datasets combined with inflammatory 48 pathway enrichment analysis, we determined that antibiotic treatment of TB, despite significantly perturbing 49 the gut microbiota, dampens the proinflammatory signature characteristic of active TB. Contrarily, an 50 investigational TB treatment that failed to clear TB, but that caused similar microbiota perturbations, 51 exacerbated peripheral inflammation. To decouple the effects of antibiotic induced changes in the microbiota 52 from Mtb sterilization as predictors of normalization of TB associated inflammation, we applied random forest 53 regression to the microbiome-transcriptome-sputum data from the two longitudinal datasets. We found 54 inflammatory renormalization is positively affected by both pathogen sterilization and by the abundance of 55 health-associated Cluster IV and Cluster XIa Clostridia. Oppositely, increases in the abundance of commonly 56 known pathobionts such as Bacilli and Proteobacteria clusters predict inflammatory exacerbation. We 57 independently investigated and validated these microbiota-peripheral inflammatory signature associations by 58 applying machine learning to the peripheral gene expression and microbiota profiling in an independent human 59 cohort of 52 healthy control individuals. Together, our findings indicate that antibiotic-induced reduction in 60 pathogen burden and changes in the microbiome are independently associated with treatment-induced 61 changes of the inflammatory response of active TB, and more broadly indicate that response to antibiotic 62 therapy may be a combined effect of pathogen killing and microbiome driven immunomodulation. Our Causal relationships between microbiome composition and host peripheral gene expression have been studied and validated extensively in mice but have not been demonstrated in humans. In this work, we investigate microbiome-periphery gene expression relationships in three unique human cohorts from Haiti, in the context 72 of a multicenter TBRU collaboration. The objectives of this work were to quantitate the relationships in two longitudinal and one cross sectional cohorts of Tuberculosis treatment. To do this we applied machine learning methods for high-dimensional data to learn associations in humans between these commonly profiled tissues. Collectively, our results provide support to the hypothesis that there exists clear links between microbiome composition and host peripheral gene expression, that can be biologically elucidated using common and well validated molecular pathway analysis. hypothesis that there exists clear regulatory relationships between gut microbiome composition and peripheral gene expression, at both the immune 5 , and gene-pathway regulatory levels in humans.


Introduction 80
There is mounting evidence that the gut microbiome has an important role in the modulation of host physiology, 81 with a wealth of studies having associated microbiome composition and functions with differential 82 inflammatory, neurological, and even behavioral activity 1 . Gastrointestinal colonization by specific taxa with 83 particular metabolic capacities has been shown to differentially modulate host biology2. For example, 84 colonization by a subset of Clostridia enhanced anti-inflammatory phenotypes in mice 3 , and enrichment in 85 specific members of the Bacteroides and Parabacteroides genera induced CD8+ T cell responses and anticancer 86 activity in mice and marmosets 4 , as well as correlating with the abundance of these immune effectors in 87 for systemic changes in inflammatory responses in humans. This knowledge gap is due in part to the difficulty 93 of isolating the microbiome dependent effects from other aspects of human physiology and in discerning the 94 direction of causality in human studies. As microbial communities in the gut promote the development and 95 maintenance of innate and adaptive immune responses 8 , including microbiota-educated immune cells and 96 many small molecules that circulate throughout the periphery 9 , we would expect to observe both localized and 97 systemic host effects upon major microbiome alterations such as treatment with antibiotics, which we can 98 measure using technologies such as shotgun RNA sequencing (RNAseq) 10 . 99 100 Infection by Mycobacterium tuberculosis (Mtb) (the 9 th leading cause of death on Earth 11 ) has been shown to 101 have a markedly different systemic gene expression profile compared to people with latent disease, other 102 respiratory diseases, or no known infection [11][12][13] . Specifically, infection with Mtb leads to heightened expression 103 of inflammatory pathways, most notably the Type I and Type II interferon pathways 14-17 , with this pattern 104 resolving with antibiotic therapy 14,17,18 . A recent meta-analysis combining microarray and RNAseq data from 105 studies aimed at identifying active TB transcriptional signatures, confirmed the findings about a specific set of 106 peripheral blood transcripts that are biomarkers of active TB disease, relative to healthy individuals or those 107 with latent TB infection (LTBI) 19 . Antibiotic treatment for active TB involves combination therapy with narrow 108 spectrum and prodrug agents with mostly Mycobacterial-specific targets, (HRZE) is given for two months and is 109 then followed by an HR-only administration for an additional four months, in order to achieve over 95% 110 likelihood of Mtb clearance 20 . The disruptive effect of HRZE therapy on the intestinal microbiome was 111 5 demonstrated in a longitudinal study in mice 21 and cross-sectional study in humans 22 , which indicated that the 112 major phyla perturbed are from the class Clostridia, a group of obligate anaerobes in the gut with well described 113 immunomodulatory effects on the host 2,3,23,24 . 114 115 Given that HRZE treatment causes GI microbiota shifts that include the depletion of many Clostridia species, 116 and given the role that these species play in modulation of host biology in mice and humans, we reasoned that 117 there could be a connection between the microbiome alterations observed during HRZE therapy and the 118 resolution of systemic inflammatory responses to TB. However, because HRZE therapy rapidly reduces the 119 burden of Mtb in the early phase of treatment, it is difficult to uncouple the immunologic effects of pathogen 120 killing from microbiota perturbation without a control group that has either pathogen killing or microbiome 121 perturbation, but not both. consists of a randomized clinical trial of study volunteers, where we collected disease severity measurements 127 (TTP), microbiome profiling, and peripheral transcriptomics at baseline, before active TB patients were 128 randomized to either HRZE (standard of care TB treatment), or NTZ. Cohort 2 (longitudinal) consists of study 129 volunteers who were followed throughout the course of 6 months of TB treatment, where we collected TTP, 130 microbiome, and transcriptomics data. Finally, Cohort 3 (cross sectional) consists of healthy volunteers. These 131 6 healthy volunteers were enrolled separately through the Tri-I TBRU. Around half are healthy and TB-negative 132 household contacts of active TB-patients, and the other half are community controls, with no know TB exposure. 133 We performed microbiome profiling and peripheral transcriptomics on these individuals as well. B. Numbers of 134 individuals in this study. C. Causal inference diagram demonstrating the major hypothesis this study tests.

136
To address this problem, we combined three independent clinical datasets for which we had gathered 137 microbiome profiling via 16S rRNA sequencing, peripheral blood transcriptomics, and where relevant Mtb 138 abundance in the sputum. The first dataset ( Figure 1A,B) consists of the secondary endpoint data from a clinical 139 trial (NCT02684240) that compared the early bactericidal effect (EBA) of standard tuberculosis (TB) therapy 140 isoniazid (H), rifampin (R), pyrazinamide (Z), and ethambutol (E) (HRZE) to the antiparasitic drug nitazoxanide 141 (NTZ), shown to possess antimycobacterial activity in vitro 25,26 . We thought this dataset to be ideal for testing 142 our hypothesis because (1) the clinical trial demonstrated that as opposed to HRZE, NTZ did not show any 143 antimycobacterial activity in humans 26 ( Figure 2B), and (2) preliminary microbiome analysis by our group 144 showed that NTZ caused microbiome dysbiosis similar to that caused by HRZE, thus providing us with the desired 145 dysbiosis in the absence of any pathogen killing or aggravation of disease. We then found that HRZE and NTZ 146 had distinct effects on host peripheral gene expression, with HRZE causing a resolution of inflammatory gene 147 signatures, and NTZ causing an exacerbation of them in the absence of disease progression. Because we had 148 used secondary endpoint observations from a clinical trial dataset that was powered for its primary endpoint 149 (reduction of Mtb in the sputum) and that only spanned a two-week treatment time frame, we additionally 150 collected data from an independent 6-month longitudinal HRZE treatment cohort ( Figure 1A,B). This cohort not 151 only shows similar resolution of disease with treatment compared to the two-week HRZE arm, but also allowed 152 us to characterize for the first time the long-term effects HRZE treatment on the intestinal microbiota and 153 peripheral gene expression. We took advantage of these datasets to train Random Forest Regression models to 154 assess the changes in expression of peripheral inflammatory pathways as a function of changes in microbiota 155 species abundances and simultaneous changes in Mtb in the sputum. Our modeling predicted a number of 156 bacteria to differentially affect proinflammatory responses, which were consistent with previously published 157 experimental reports. Specifically, we found that the increase in abundance of SCFA-producing Clostridia species 158 predicts inflammatory renormalization; oppositely, we found that an increase in abundance of oxygen-tolerant 159 pathobionts such as E. coli, Klebsiella, Enterococcus, and Streptococcus species associate with exacerbation of 160 proinflammatory pathways ( Figure 1C). Finally, we investigated these microbiota-peripheral gene expression 161 relationships using a cohort of healthy Haitian community controls and healthy household contacts of TB 162 patients, previously described 5 ( Figure 1A,B) and confirmed the existence of these microbiome-peripheral 163 immune gene expression signatures. Overall, we believe that these results provide support to the often stated 164 hypothesis that there exists clear regulatory relationships between gut microbiome composition and peripheral 165 gene expression, at both the immune 5 , and gene-pathway regulatory levels in humans. 166

167
Gut microbiome diversity is depleted after two weeks of HRZE or NTZ treatment 168 As detailed elsewhere, the GHESKIO centers in Port au Prince, Haiti conducted a prospective, randomized, early 169 bactericidal activity (EBA) study in treatment-naive, drug-susceptible adult patients with uncomplicated 170 pulmonary tuberculosis (TB) (ClinicalTrials.gov Identifier: NCT02684240) 26 . Thirty-four participants were 171 randomized to receive either NTZ, 1000 mg po (oral) twice daily, or standard oral therapy with isoniazid 300 mg 172 daily, rifampin 600 mg daily, pyrazinamide 25 mg/kg daily, and ethambutol 15 mg/kg daily (referred to as HRZE) 173 for 14 days ( Figure 1A). The primary endpoint of the trial was sputum bacterial load (measured by time to culture 174 positivity, TTP) in a BACTEC liquid culture system, a sensitive method to detect disease progression and severity 175 in We have previously reported 22 that HRZE therapy depletes members of the order Clostridiales, but the cross-188 sectional design of that study did not allow for conclusions about the rapidity of this effect, and most 189 importantly, did not include pretreatment samples to allow for the assessment of baseline microbiome 190 composition. To investigate microbiome changes induced by HRZE or NTZ, we extracted and amplified bacterial 191 and archaeal DNA using V4 -V5 16S rDNA sequencing (see Methods). Stool samples were collected before the 192 start of treatment (baseline) and on day 14 of therapy ( Figure 1A). Using principal components analysis (PCoA) 193 with Bray-Curtis distances, we found that that the component accounting for the greatest variation in the 194 microbiome data qualitatively represented changes in microbiome community structure that occur after two 195 weeks of NTZ treatment ( Figure 2C). Inspecting Axis 2 of Figure 2C, we found that the observed separation 196 correlates with sequencing batch. Therefore, for any subsequent statistical modeling analysis (i.e., differential 197 9 microbiota and gene expression modeling), sequencing batch information has been controlled for by including 198 it as fixed effect in the modeling.

200
To compare the effect of the two treatments to microbiome alpha diversity we calculated the Inverse Simpson

201
Diversity Index for every microbiome sample 28 . We then regressed the Shannon Diversity Index via linear mixed-202 effect modeling as Methods). We found that there were no differences in alpha diversity at baseline between the two arms, while 204 both treatments significantly reduced microbial diversity (p<0.01, see Supplementary Table S1), with NTZ 205 treatment not having a significantly different effect compared to HRZE (p>0.05, for the interaction term, See 206 Supplementary Table S1) ( Figure 1D).   sequencing batch-dependent effects in addition to establishing effects that are due to treatment group (pre-233 treatment, HRZE, NTZ). We used the 1| random effect to control for baseline differences among individuals. 234 ASVs significantly affected by the treatments were determined using a Benjamini-Hochberg false discovery rate 235 (FDR) adjusted p-value of 0.05 (see Methods).

248
Taken together these results demonstrate that NTZ, despite having no effect on TTP or disease severity, causes 249 a perturbation of the intestinal microbiota which is more pronounced compared to that caused by HRZE. This  250 includes the depletion of a large number of Clostridia, and the selection for known disease-associated, oxygen-251 tolerant pathobionts.  As the above results highlighted a differential effect on TB disease severity and intestinal microbiome 267 composition of the two antimicrobials tested, we recognized that this presented a unique opportunity to infer 268 possible relationship between microbiome composition, Mtb bacterial load, and peripheral gene expression. As 269 the patients were randomized at baseline before being assigned to the two treatments, we used linear mixed-270 effect modeling with Limma/Voom to model the abundance of each host transcript as which are consistent with the broad immunologic effects of antibiotic induced reduction in the levels of a 280 bacterial pathogen 14,17,36 ( Figure 4A). In contrast, and to our surprise, NTZ treatment showed the opposite effect. 281 Inflammatory signaling pathways reduced by HRZE, including TNF signaling, IFN signaling, and type 1 282 interferon signaling were all significantly enriched by NTZ treatment at day 14 ( Figure 4A). Several other 283 pathways such as hypoxia, apoptosis, and reactive oxygen species (ROS) that are considered hallmarks of 284 immune dysregulation 37 , were also enriched by NTZ treatment. As NTZ was found to perturb the microbiome 285 while keeping Mtb bacterial load in the sputum (TTP) substantially unchanged, we hypothesized that the NTZ 286 effect was likely not completely due to prognostic disease progression, and rather at least a partial function of 287 microbiome alteration (Supplementary Table S6). 288 289  subjects. We defined three classes of changes to these transcripts with two weeks of HRZE or NTZ treatment: 1) 310 renormalization (transcripts whose pre-post HRZE/NTZ fold change in expression displays the same sign (or 311 direction) of the previously-reported fold-change between active TB and control/LTBI; 2) unchanged (transcripts 312 with no change in expression between pre-post HRZE/NTZ administration); and 3) exacerbation (genes whose 313 pre-post HRZE/NTZ fold-change sign is opposite to the previously-reported fold-change between active TB and 314 control/LTBI). 157 of these active TB signature genes were significantly affected by HRZE ( < 0.05) ( Figure  315 4B, Supplementary Table S3). Of these 157 affected genes, 144 (92%) were found to renormalize with the 316 treatment (i.e., displaying the same direction of the fold change reported for active TB vs. control individuals), 317 while 13 (8%) were found to exacerbate (i.e., opposite direction of the fold change). On the other hand, only 318 four of these TB-related inflammatory genes were found to be affected by NTZ, and all of them (100%) displayed 319 an exacerbation trend ( Figure 4C). 320 321 Because of the disruption in microbiome composition in the absence of a change in disease severity observed 322 in NTZ-treated people, we hypothesized that there may be other subsets of host genes that are linked to 323 microbiome-dependent immunity that could be responsive to the observed antimicrobial perturbations. To test 324 this, considering the already well-established link between microbiome dysbiosis and autoinflammatory 325 14 conditions such as Inflammatory Bowel Disease (IBD) 39 , we selected a recently published panel of 880 genes 326 differentially expressed in colon biopsies from IBD patients and asymptomatic controls 38 . Of these 880 genes, 327 364 were detected in our dataset (expected given that we are profiling whole blood transcriptomics, rather than 328 colon biopsy tissue). Despite the more limited effect on the microbiome compared to NTZ, HRZE administration 329 was found to be responsible for a change in 117 of these genes ( < 0.05) ( Figure 4D, Supplementary Table  330 S3), while NTZ was found to be responsible for a change in 55 of them (FDR< 0.05) ( Figure 4E, Supplementary 331 S3). Similarly, we defined as renormalization those genes that have a post/pre fold change due to antimicrobial 332 treatment of the same sign as control/active IBD. Interestingly, while NTZ affected a smaller subset of these IBD 333 immune genes, it caused greater exacerbation compared to HRZE (77.8% vs 42.0% genes exacerbating).

335
Longitudinal profiling in an early bactericidal activity HRZE treatment cohort

336
To validate the observations obtained after two weeks of HRZE and to identify effects that may occur with 337 increased treatment length we enrolled a longitudinal treatment cohort (20 individuals) to measure EBA over 338 the course of the 6 months of HRZE treatment, with periodic sampling of sputum for TTP, stool for microbiome 339 profiling, and whole blood for peripheral RNA sequencing. All participants received standard of care HRZE 340 therapy (isoniazid 300 mg daily, rifampin 600 mg daily, pyrazinamide 25 mg/kg daily, and ethambutol 15 mg/kg 341 daily). Sputum mycobacterial load (TTP) was collected at baseline, day 7, day 14, one month, two months, and 342 six months at the completion of treatment. Stool from microbiome profiling was collected at each of these 343 timepoints as well. Whole blood was collected at baseline, day 14, and two months. Stool and peripheral blood 344 samples were processed for microbiome and transcriptomic profiling as for the above described clinical trial 345 (See Methods).

347
Similarly, to what we observed in the HRZE arm of the clinical trial, TTP is significantly higher compared to 348 baseline after two weeks of treatment (p <0.05). TTP also significantly increases after two months compared to 349 the two weeks timepoint confirming the presence of an additional sterilization effect of the antibiotic on Mtb 350 (p <0.05, linear mixed-effect modeling, see Methods) ( Figure 5A). Interestingly, microbiome diversity drops after 351 just 7 days of treatment, increases after 1 month of treatment but remains significantly lower compared to 352 baseline at the 6-month follow-up timepoint (p < 0.05, linear mixed-effect modeling, see Methods) ( Figure 5B).

354
We performed Limma/voom differential analysis to determine effects of treatment and time on the microbiome 355 and peripheral host transcriptomics using the same approach as above. As observed in the clinical trial and in 356 our previous work 22 , HRZE was found to depress Clostridiales after one week of treatment, with most of these 357 Clostridiales remaining significantly depressed compared to baseline, even at the 6-month follow-up time point 358 ( Figure 5C). At day 7 compared to baseline, 19 ASVs were depleted and 6 were increased in abundance, at day 359 14 compared to baseline, 74 ASVs were depleted and 2 were increased in abundance, at one month compared 360 to baseline, 97 ASVs were depleted and none were increased in abundance. Thus, during the first month of 361 treatment, microbiome depletion relative to individual baseline samples was evident, however, later in the 362 course of treatment, we observed a different trajectory. Relative to baseline, at two months of HRZE, we 363 observed 13 ASVs depleted in abundance, while 95 ASVs increased in abundance. At the 6-month mark of HRZE 364 treatment, we observed 3 ASVs depleted in abundance while 96 ASVs increased in abundance ( Figure 5C). The 365 increased abundance of specific ASVs at the two month and 6-month mark appeared to be relatively 366 heterogeneous between individuals, and overall, this longitudinal analysis suggests that after one month of 367 HRZE therapy, most ASVs that would be affected by therapy are depleted (Supplementary Table S4).

369
With respect to the peripheral host inflammatory profiling, we observed distinct changes in gene signatures at 370 two weeks (day 14) and two months of treatment, compared to baseline (Supplementary Table S6). We 371 observed a similar decrease in common inflammatory pathways in the Hallmark pathway dataset described in 372 the clinical trial ( Figure 5D). Interestingly, comparing Day 56 to baseline or to Day 14, we see additional reduction 373 in inflammatory gene signatures, potentially explained in part by the further reduction in TTP at Day 56 ( Figure  374 5E-I). Collectively, this EBA dataset confirms many of the findings observed in the HRZE arm of the trial, 375 comparing Day 14 to baseline in both cases, and adding additional treatment information at these later 376 timepoints. Importantly, it can be used in addition to the clinical trial data for training the host-microbiome-Mtb 377 model described subsequently. 378  . We then regressed this quantity against the 400 corresponding fold change in abundance of for every ASV v in the same interval +1 ] =1… and against the 401 corresponding fold change in TTP, +1 . Using change from baseline values accounts for the random effect of 402 each subject without having to incorporate this into the model statement. We solved this regression problem 403 using Random Forest regression as in 41 . To train the models we used all the observations from the clinical trial 404 and from the longitudinal EBA cohort for a total of 34 paired samples (Figure 1). We fit a model for each 405 inflammatory using all the data from the three patients' group (HRZE clinical trial, NTZ clinical trial and HRZE 406 EBA) because we wanted to find patterns that are general across multiple datasets. Each model was trained 407 using 5000 trees and a with train-validation partitioning of 80-20% of the data. We reasoned that this approach 408 was appropriately suited for this type of "large p, small n" multi-omics dataset common in clinical research 42 . 409 Importantly RFR modeling has significant advantages compared to traditional multi-linear regression 410 techniques, because it is agnostic to model structure (e.g. non-parametric regression), it does not need to meet 411 common assumptions underlying classical regression techniques and is able to intrinsically perform ranked 412 feature selection. Importantly, while the interpretation of RFR is apparently less immediate compared to 413 traditional regression (e.g. there are per-se no regression coefficients or betas), downstream analysis, which 414 includes Permutated Importance 43 and Accumulated Local Effects (ALE) calculations 44 (see Methods) allows for 415 the estimation of the significance of predictors (e.g. TTP, microbiome constituents, etc.) and of their effects on 416 the dependent variable (e.g. host peripheral inflammatory markers). 417

418
When plotting the average slope of the ALE curves for predictors with significant Permutated Importance value 419 (p < 0.05), our analysis identifies the increase in TTP (and therefore a decrease of Mtb in the sputum) and the 420 increase in abundance of ASVs from Clostridia, especially members from the Cluster IV and XVIa groups 45 , which 421 has been shown to induce anti-inflammatory responses (e.g. Treg-induction) 2,3 and also comprehends SCFA-422 producing species 46  in abundance of oxygen-tolerant pathobionts including Enterococcus, Streptococcus, and E. coli is found to 449 predict inflammatory exacerbation. 450

Relationship of the microbiome and peripheral gene expression in a healthy control validation cohort 451
The results from our machine learning modeling on the data from both longitudinal treatment cohorts provide 452 support to the hypothesis that specific intestinal microbiota members are associated with immune-related 453 peripheral blood gene signatures in humans. Specifically, our modeling predicts that higher abundance of 454 Clostridia is negatively associated with inflammation (e. have broad qualitative differences in distinct biological pathways ( Figure 7A). While these pathways are not by 473 any means completely representative of everything happening biologically in these individual, we do feel that 474 they highlight clear trends across the course of the TB spectrum of disease, compared to a control population 475 of healthy volunteers.

477
We performed RFR (see Methods) for each pathway against the microbiome space, and summarized the findings 478 in Figure 7B. Surprisingly, we found the abundance of a large number of Firmicutes and particularly Clostridia to 479 be associated with a number of the characterized Hallmark molecular pathways (Supplementary Table S7). Even 480 more intriguingly we found that higher abundance of ASVs that are mapped to the health-associated F. if the pathway positively associates with a particular ASV, or negative if the pathway negatively associates with 498 a particular ASV. 499 500

501
Since the advent of high-throughput microbiome characterization, it has become clear that antibiotics are one 502 of the most common and severe perturbing influences on human microbiome composition, with both acute and 503 longer lasting effects 53,54 . It also has become evident that the specific microbiome constituents have specific 504 effects on host immunity, including the abundance and function of immune cell subsets 24 . Significant prior data 505 have documented the effects of antibiotics on microbiome composition and function and the consequent 506 influence of these microbiome factors on immune cell populations 55 , with the majority of these findings derived 507 using in vivo mouse models. While there is no doubt that microbiota dynamics affect host immunity 6 , it remains 508 unknown to what degree antibiotic induced perturbation of the microbiome may modify the outcome of 509 treatment of infection, or what relationships exist in humans between gut microbiota composition and 510 peripheral gene expression. It is conceivable that antibiotics work to clear infection both due to direct pathogen 511 killing and by immune modulation through the microbiome. It is also possible that the pathogen killing effect of 512 antibiotics may be partially counteracted by detrimental immune effects induced by microbiome perturbation. 513 Such dynamics may be particularly relevant to the treatment of chronic infections such at tuberculosis, in which 514 antibiotic therapy is prolonged and the disease manifestations reflect a mixture of pathogen burden and the 515 balance of inflammatory mediators that cause tissue destruction and chronic symptomatology 56,57 . 516 517 Antibiotic sensitive tuberculosis is treated with six months of antibiotics with predominantly mycobacterial 518 specific agents. In this study we report the early and late microbiome effects of HRZE therapy in subjects with 519 active TB and demonstrate that the same changes observed in a human cross-sectional study of TB treatment 22 520 (comparing vs cured and LTBI individuals) were present after just two weeks of treatment. As previously 521 shown 22 , we conclude that HRZE treatment has a rapid and narrow effect on the intestinal Class of Clostridia, a 522 finding that was also demonstrated in mice 21,58 . We note that given the mycobacterial-specific nature of TB 523 drugs, and the combinatorial nature in which small molecules interact to affect the microbiome, it was difficult 524 to predict that primarily Clostridia, in the Phylum Firmicutes, would be targeted by HRZE therapy, whereas 525 Actinobacteria, the phylum to which Mtb belongs, are relatively unaffected. Experiments in mice have 526 demonstrated that this anti-Clostridia effect is primarily driven by rifampicin 21 . Clostridia are immunologically 527 active components of the microbiota through their production of metabolites such as short chain fatty acids 528 and other compounds 2,5,6,59,60 . 529 530 The clinical trial allowed us to dissect the relative contributions of pathogen killing and microbiome perturbation 531 to disease resolution because one treatment arm, standard therapy, both reduced Mtb bacterial burden and 532 25 perturbed the microbiome, whereas NTZ had no effect on average Mtb burden, but did perturb the microbiome 533 in a fashion that overlapped with HRZE. We were able to extend these findings in the EBA dataset with both an 534 investigation into microbiome changes after day 7, day 14, one month, two months, and 6 months post 535 treatment, compared to Baseline, as well as peripheral gene expression and pathway changes after day 14 and 536 two months of HRZE therapy. The findings from these analyses may provide support to the hypothesis that 537 antibiotic perturbation of the microbiome has systemic effects on peripheral gene expression. Further, we 538 speculate that the large heterogeneity in the rapidity of treatment response in TB may be a partial function of 539 heterogeneity in the effects of antibiotics on the microbiome. 540 541 To validate the inferred microbiome-host inflammatory relationship, we mined microbiome and blood 542 transcriptomic profiling from an independent human cohort of healthy Haitian individuals. Remarkably, despite 543 the reduced peripheral levels of inflammatory pathways compared to subjects with active TB, we again observe 544 that members of Clostridium IV and XIVa predict a reduction in the expression in pro-inflammatory pathways. 545 This validation data strongly supports our conclusion that microbiome composition sets the tone of systemic 546 inflammation, both in disease states and in homeostatic conditions. Further, it is consistent with the prior 547 findings in both humans and animals that Clostridia have been associated with induction of anti-inflammatory 548 states 61,62 . 549 550 Finally, given the challenge of explaining the relationship between microbiome composition and peripheral gene 551 expression with paired samples, randomized to drug combinations with vastly different effects on both body 552 systems, we strove to use appropriate mathematical approaches for this type of analysis. For single-omic 553 microbiome and RNAseq data, we chose to use limma/voom to model changes in these data given their ability 554 to model many effects while using subject's baseline values as random effects. For multi-omics integration, we 555 used Random Forest Regression. While there are a variety of statistical and machine learning techniques able 556 to investigate relationships between complex multiparametric datasets ("large p": microbiome composition, 557 RNAseq data, clinical metadata, randomization cohort, paired-sample baseline normalization, etc., and a "small 558 n" of individuals in early phase clinical trials), Random Forests are adequate for microbiome purposes, as they 559 have been shown to outperform Support Vector Machines in some instances, especially for continuous variable 560 data, and need initialization of a smaller set of parameters compared to other deep-learning methods. We 561 believe that our results highlight the utility of these models to: 1. Provide evidence for or against a particular 562 hypothesis about clinically significant relationships between many potentially related parameters, and 2. To 563 26 provide hypothesis generating relationships between the multi-omic constituents (i.e., features) of these 564 models, which can be further tested in mice, validation cohorts, or other model systems.

566
Our data indicates that within the first 14 days of treatment of tuberculosis, resolution of the active 567 inflammatory response of TB (as measured by peripheral blood transcriptomics) may be strongly affected both 568 by reducing Mtb burden as well as through antimicrobial-induced microbiome perturbations that may act 569 directly on systemic immune function. Among the pathways tightly correlated with both factors are the 570 signature activated pathways of active TB disease: IFN, type I interferon, and TNF 14 . There is growing evidence 571 that the outcome of active TB reflects a mixture of pathogen burden and cytokine networks that include IL-1 572 and IFN, with the latter acting to exacerbate disease 56 . Our findings indicate that the microbiome perturbation 573 that accompanies TB treatment is a predictor of the normalization of these same pathways during early 574 treatment, suggesting that microbiome perturbation could modify or predict the rapidity of disease resolution. 575 In the first two weeks of treatment, pathogen killing is the dominant factor, but microbiome dependent 576 modulation of inflammatory responses during treatment may assume an important role during the later phases 577 of treatment when pathogen killing slows. The validation of the relationships between microbiome composition 578 and peripheral gene expression in a healthy control cohort, especially for the collective expression of these same 579 pro-and anti-inflammatory pathways, suggests that these relationships may extend into other populations. 580 Whether these relationships are causal, or biomarkers of another state will remain at the forefront of future 581 study design. Future work will be directed to applying the analytical tools and study design presented here to 582 later time points in the TB treatment course to examine whether microbiome perturbation during treatment 583 associates with clinically relevant treatment outcomes, and whether the abundance of Clostridia correlates with 584 rapidity of Mtb sterilization or the resolution of the inflammatory response that accompanies active TB. Such 585 data might help support trials to test microbiome monitoring as a predictor of TB treatment outcome or help 586 understand interindividual heterogeneity in treatment outcomes. 587 588

Materials and Methods 589
Ethical statement and study approval 590 All volunteers provided written informed consent to participate in this study. • Batch represents the sequencing batch information 669 Alpha diversity indices were computed using phyloseq package in R and the implementation of linear mixed 670 effect models were carried out using nlme package in R.

671
For the EBA longitudinal cohort, to identify the significance of the influence of sex, age, and treatment time on 672 time to positivity (TTP), we implemented linear mixed effect model as Similarly, to identify the significance of the influence of sex, age, and treatment time on microbiota diversity on 674 microbiota diversity (Inv Simpson), we performed linear mixed effect modeling Differential analysis for microbial ASVs and host genes: Both raw 16Sr rRNA microbiota ASVs and peripheral 678 blood RNAseq gene-expression counts were modeled using the limma/voom pipeline. 29 This allowed us to use 679 linear mixed-effect modeling of gene/ASV counts as of ~ + ℎ + + 1| . This model 680 statement enables quantifying sex and sequencing batch-dependent effects in addition to establishing effects 681 that are due to treatment group (pre-treatment, HRZE, NTZ). For EBA cohort, we used similar differential 682 analysis approach as clinical trial by modeling gene/ASV counts as Time represents the Day 0, Day 14 and Day 56 time points. The advantage of limma is that we could use 684

1|
as a random effect to control for baseline differences among individuals, important in this clinical 685 setting. Significance of ASVs affected by the treatment were determined using a Benjamini-Hochberg false 686 discovery rate (FDR) adjusted p-value of 0.05 from the modeling-produced contrast lists (e.g. HRZE vs. Pre, NTZ 687 vs Pre, Day 0 -EBA vs. Day 14 -EBA) 29 .

689
To determine how the anti-TB treatment affects both microbiome and peripheral gene expression profiles we 690 performed differential analysis on the counts data obtained by microbiome DNA and peripheral blood RNA 691 sequencing. As the primary endpoint of the clinical trial was powered to determine differences in Mtb load 692 (TTP), we determined the statistical power available to identify significant differences in the abundance of both 693 microbiota ASVs, and in the expression of peripheral genes. We ran power calculations to determine that with 694 16 pre and 16 post treatment microbiome samples and 8 pre and 8 post treatment RNAseq samples for the 695 HRZE cohort, with 80% power at a significance level () of 0.05, we could detect a fold change of 1.4 for 696 microbiome difference and a fold change of 1.8 for mRNA transcripts. In the NTZ cohort, with 18 pre and 18 697 post treatment microbiome samples and 14 pre and 14 post treatment RNAseq samples, with 80% power at 698 <0.05, we can detect a fold change of 1.4 for microbiome differences and a fold change of 1.6 for mRNA 699 transcripts. In the EBA cohort with 20 baseline, 10 day 7, 18 day 14, 13 one month, 13 two month, and 11 six 700 month follow up microbiome samples, with 80% power at a significance level () of 0.05, we can detect a fold 701 change of 1.4 (day 7), 1.36 (day 14), 1.4 (one, two, and 6 months). In the EBA cohort with 19 baseline, 19 day 702 14, and 13 two-month RNAseq samples, with 80% power at a significance level () of 0.05, we can detect a fold 703 change of 1.4 for mRNA transcripts at day 14, and 1.5 at two months. Power calculations were performed with 704 the RNAseqPower package in R. For microbiome data we calculated a biological coefficient of variation of 0.3, 705 and for RNAseq, we used a coefficient of variation of 0.4. We estimated the expected minimum fold change that 706 we could observe for each group based on the sample size, sequencing depth, and an  of 0.05. To visualize 707 trends of transcript fold changes, we used scatter plots and calculated post-vs-pre fold changes for all transcripts 708 throughout.

711
We assessed the relative contribution of the intestinal microbiota and Mtb dynamics towards peripheral gene 712 expressions using Random forest regression (RFR). Instead of modeling each gene/transcript profile as a 713 function of microbiome and TTP, we mapped our gene expression data to a set of 50 Hallmark Pathways via 714 GSVA. To avoid having correlated samples from same individual in a model, we instead modeled the changes in 715 normalized enrichment score of hallmark pathways at two time points (Day 0 and Day 14) as a function of change 716 Variance stabilized transformed (vst) counts derived from DESeq2 were used as input into the GSVA function in 749 the GSVA R package with default parameters (kcdf="Gaussian") and scaled Normalized Enrichment Scores (NES) 750 were plotted as heatmaps. Importantly, unlike classical GSEA, this analysis is agnostic to sample phenotype. 751 Identification of differential pathways post antibiotic treatment in both the clinical trial and EBA cohort was 752 performed using linear mixed effect model where we modeled the normalized enrichment score (NES) of each 753 pathways as