Gut microbiome contributes to impairment of immunity

Environmental Microbiology (2017) 00(00), 00–00

doi:10.1111/1462-2920.14015

Gut microbiome contributes to impairment of immunity in pulmonary tuberculosis patients by alteration of butyrate and propionate producers

Abhijit Maji,1,2† Richa Misra,1,3,4† Darshan B. Dhakan,2† Vipin Gupta,3 Nitish K. Mahato,3 Rituja Saxena,2 Parul Mittal,2 Nitin Thukral,5 Eshan Sharma,6 Anoop Singh,3 Richa Virmani,1,6 Mohita Gaur,3 Harshvardhan Singh,7 Yasha Hasija,5 Gunjan Arora,1 Anurag Agrawal,1 Anil Chaudhry,8 Jitendra P. Khurana,6 Vineet K. Sharma,2* Rup Lal3* and Yogendra Singh1,3* 1 Department of Microbial Pathogenesis, CSIR-Institute of Genomics & Integrative Biology (IGIB), Mall Road, Delhi, India. 2 Department of Biological Sciences, Indian Institute of Science Education and Research (IISER), Bhopal, India. 3 Department of Zoology, University of Delhi, Delhi, India. 4 Department of Zoology, Sri Venkateswara College, University of Delhi, India. 5 Department of Biotechnology, Delhi Technological University, Delhi, India. 6 Department of Plant Molecular Biology, University of Delhi, Delhi, India. 7 Department of Biochemistry, Hindu Rao Hospital, Malka Ganj, Delhi, India. 8 Department of TB and Chest, Rajan Babu Institute of Pulmonary Medicine and Tuberculosis (RBIPMT), Kingsway Camp, Delhi, India. Summary Tuberculosis (TB) is primarily associated with decline in immune health status. As gut microbiome (GM) is implicated in the regulation of host immunity and metabolism, here we investigate GM alteration in TB patients by 16S rRNA gene and whole-genome shotgun sequencing. The study group constituted of patients with pulmonary TB and their healthy household contacts as controls (HCs). Significant alteration Received 19 April, 2017; accepted 30 November, 2017. *For correspondence. E-mail: [email protected]; Tel: 191-9871095673; Email: [email protected]; Tel: 191-7552691401; Email: ruplal@ gmail.com; Tel: 191-11-27666254. †These authors contributed equally to the manuscript. C 2017 Society for Applied Microbiology and John Wiley & Sons Ltd V

of microbial taxonomic and functional capacity was observed in patients with active TB as compared to the HCs. We observed that Prevotella and Bifidobacterium abundance were associated with HCs, whereas butyrate and propionate-producing bacteria like Faecalibacterium, Roseburia, Eubacterium and Phascolarctobacterium were significantly enriched in TB patients. Functional analysis showed reduced biosynthesis of vitamins and amino acids in favour of enriched metabolism of butyrate and propionate in TB subjects. The TB subjects were also investigated during the course of treatment, to analyse the variation of GM. Although perturbation in microbial composition was still evident after a month’s administration of antiTB drugs, significant changes were observed in metagenome gene pool that pointed towards recovery in functional capacity. Therefore, the findings from this pilot study suggest that microbial dysbiosis may contribute to pathophysiology of TB by enhancing the anti-inflammatory milieu in the host.

Introduction Mycobacterium tuberculosis, the causative agent of Tuberculosis (TB), latently infects one-third of the world’s population, with predilection mainly for lungs (World Health Organization report, 2014). Extensive research over the years has increased our knowledge immensely about M. tuberculosis and its interaction with the host (Sachdeva et al., 2010; O’Garra et al., 2013; Stanley and Cox, 2013; Maji et al., 2015; Sajid et al., 2015), yet we are unable to fully understand the complexities of the disease process such as latency and reactivation. Epidemiological studies have revealed that socioeconomic conditions play an important role in the progression of the disease (Dye et al., 2009). Various factors, which are effectively associated with immune surveillance breakdown, such as malnutrition, human immunodeficiency virus (HIV) infection or administration of immunosuppressive medication, contribute to the manifestation of the active form of the disease. To better comprehend the TB disease process, it is imperative to

2 A. Maji et al. elucidate the critical factors governing the modulation of immune response. Over the recent past, with the help of next-generation sequencing technologies, we have begun to understand the contribution of gut microbiome (GM) in development and homeostasis of the immune system, host nutrition and metabolic disorders; however, its role in infectious diseases is still not well understood (Gill et al., 2006; Qin et al., 2010; Qin et al., 2012; Karlsson et al., 2013). Dysbiosis of GM composition and function has been implicated in disease development in lungs by the proposed ‘gut-lung axis’ (Samuelson et al., 2015). We are just beginning to understand this cross-talk between the gut and respiratory mucosal surfaces, which allows bidirectional passage of microbial products, metabolites and cytokines. Although evidence of GM influence has been highlighted in lung disorders such as asthma and cystic fibrosis, we know very little about respiratory infections (Marsland et al., 2015; Budden et al., 2017). The primary objective of this study was to determine whether individuals with active M. tuberculosis infection display alterations in their GM and as a secondary objective examine the short-term effects of the anti-TB drugs on the microbial composition. Here, we report for the first time an in-depth functional characterization of GM in Indian subjects, an ethnic cohort which features in the high burden list for TB in the World Health Organization (WHO) report, 2014. We compare the taxonomic and functional diversity of GM from faecal samples of six treatment-na€ıve TB patients and follow-up the same patients, receiving directly observed treatment, short-course (DOTS) regimen, one week and one month after the beginning of treatment. The six healthy household contacts (HCs) were chosen as the controls to rule out biases of genetic background and influence of the environmental conditions. Despite small cohort size, statistically significant differences in microbial composition were obtained between the groups. The results obtained from this pilot study provide insights into the role played by the GM in a pulmonary infectious disease. Results Clinical characterization reveals low BMI and cholesterol levels in TB patients Patients displaying symptoms of persistent coughing, fever and anorexia were prospectively recruited for this study at the Rajan Babu Institute of Pulmonary Medicine and Tuberculosis (RBIPMT), Delhi. These patients were subjected to sputum smear microscopy for detection of M. tuberculosis bacilli (acid-fast bacilli staining) and chest X-ray examination for the presence of any lesions. Smeartest positive cases were classified as pulmonary TB patients. X-ray results also confirmed the presence of cavitary lung lesions in all patients. The study exclusion criteria

for patients were any previously diagnosed serious medical condition such as any form of TB, HIV seropositive, cancer and recent (< 1 months prior) use of any antibiotics, to avoid the effect of confounding factors. Based on these criteria, the patient group consisted of six pulmonary TB patients, with a median age of 28.5 years (range 14–47 years). Faecal samples were collected prior to the start of DOTS treatment from these patients in sterile containers and termed as TBZ (TB group at day Zero). The patients were put on the standard DOTS regimen and samples were collected during the treatment at two time-points, a week and a month after the commencement of treatment and termed as TBW (TB group after one week of treatment) and TBM (TB group after one month of treatment) respectively. These early treatment time-points were chosen to understand the immediate effect of perturbations in GM caused by usage of anti-TB drugs. The six HCs, a blood-relative from each patient’s family sharing the same household, constituted the control group, to minimize the impact of genotype, dietary and environmental factors on the composition of GM and matched as closely as possible for age (range 19–35 years; median 25.5 years). The HCs showed no signs of any symptomatic illness (no coughing/weight loss/fever), as confirmed by the physician and were negative for the exclusion criteria set for the patient group. The X-ray screening and tuberculin skin test could not be performed on HCs since the ethics committee denied approval for it in the absence of clinical symptoms for control subjects. Among the clinical variables tested in the study, a significant decrease in body mass index (BMI) was observed in TB patients (p 5 0.0152). This result is in coherence with earlier observations, where a strong association between decrease in BMI and increased risk of active TB has been established (Lonnroth et al., 2010). The serum levels of total cholesterol (TC) and low-density lipoprotein (LDL) were also found to be significantly lower (p 5 0.002) in TB patients than in HCs. This is also consistent with an earlier report where low levels of TC, high-density lipoprotein (HDL) and LDL were observed in TB patients (Deniz et al., 2007). In a report by Akpovi et al., recovery from TB after completion of six months of treatment was accompanied with normalization of cholesterol levels (Akpovi et al., 2013); however, improvement in cholesterol levels was not observed in this study, perhaps due to the early treatment time-point chosen for the analysis. The clinical data and details for these participants are summarized in Table 1 and Supporting Information Additional File S1: Table S1, respectively. Dysbiosis of gut microbiota in TB patients To explore the association of GM with pulmonary TB, faecal microbiome profiling of all subjects was carried out

C 2017 Society for Applied Microbiology and John Wiley & Sons Ltd, Environmental Microbiology, 00, 00–00 V

Gut microbiome of TB patients 3 Table 1. Summary of clinical variables in study groups. Mann–Whitney U Test (p values) Subject details

HC (n 5 6) median (range)

TBZ (n 5 6) median (range)

Age (years) 25.5 (19 to 35) 28.5 (14 to 47) BMI (z-score) 0.72 (0.13 to 1.26) 21.13 (–1.39 to TC mmol/l (z score) 1.48 (1.20 to 1.68) 20.32 (–1.06 to HDL mmol/l (z score) 0.87 (0.04 to 1.57) 20.98 (–1.48 to LDL mmol/l (z score) 1.59 (1.35 to 1.90) 20.37 (–0.89 to

0.78) 0.35) 1.57) 0.23)

TBW (n 5 6) median (range)

TBM (n 5 6) median (range)

HC vs. HC vs. TBZ vs. TBZ TBM TBM

Same as TBZ Not available 20.28 (–1.73 to 0.05) 20.53 (–1.10 to 1.44) 20.61 (–0.97 to 20.09)

Same as TBZ Not available 20.67 (–1.39 to 0.16) 20.21 (–1.23 to 1.31) 20.85 (–0.93 to 20.01)

0.8528 0.0152 0.002 0.1234 0.002

– – 0.002 0.0195 0.002

– – 0.3312 0.4632 0.2511

HC, household contacts; TBZ, TB patients at day zero; TBW, TB patients after week treatment; TBM, TB patients after month treatment; BMI, body mass index; TC, total cholesterol; HDL, high-density lipoprotein; LDL, low-density lipoprotein. p 0.05 indicated by bold font.

using 16S ribosomal RNA (16S rRNA) gene sequencing and whole-genome shotgun sequencing. In total, we generated 16,168,634 high-quality paired-end reads (n 5 673,681 6 146,659 (mean 6 SD) reads per sample) for 16S rRNA V3 hypervariable region (see Supporting Information Additional File S2) and an average of 41,462,376 paired-end reads (6 4,564,434 SEM) totalling of 199.019 gigabases of high-quality whole genome DNA sequence data (see Supporting Information Additional File S3). Sequencing data were filtered for human and lowquality reads and further analysed as described in the Experimental Procedures section. The 16S rRNA gene reads were clustered into 3022 Operational Taxonomic Units (OTUs) at 97% identity. Since our cohort included only one vegetarian family and diet has been identified to influence microbial diversity, we tested top five principal components (PCs) for correlation with covariates such as diet, age, BMI and disease status (see Supporting Information Additional File S1). We found diet, BMI and disease status to significantly influence the variability in microbial composition based on UniFrac distance matrix profile (Supporting Information Additional File S1: Fig. S1a,b; Supporting Information Additional File S4: Table S4a), so the vegetarian subjects were excluded from 16S rRNA-based diversity and taxonomic analysis to remove confounding effect of diet on the GM. However, all samples were included for metagenomic data analysis since only TB disease status was a significant factor to influence the composition of Hellinger distance-based microbial gene diversity (Supporting Information Additional File S1: Fig. S1c,d; Supporting Information Additional File S4: Table S4b). An additional PERMANOVA (permutational multivariate analysis of variance) (McArdle and Anderson, 2001) analysis revealed that among the clinical variables, TC, LDL and TB disease status are the significant factors associated with the variation in UniFrac distances (Table 2). Next, the OTU table was rarefied at equal depths starting from 100 sequences to 270,000 sequences per

sample for maximum coverage of phylotypes. Principal component analysis (PCA) of unweighted UniFrac distances from the rarefied OTU table (depth 270,000 sequences/sample) showed that TBZ samples clustered separately from HC samples (Fig. 1A, see Supporting Information Additional File S5). We quantified the alpha diversity of all samples by employing Shannon index for biodiversity, OTU species count and phylogenetic diversity. The observed species and Shannon index indicated greater diversity and richness of species among the TBassociated microbiota (Table 3 and Fig. 1B). UniFrac analysis confirmed the dysbiosis of GM in TB patients, with significantly greater interpersonal variation seen among TB patients than in HCs (Fig. 1C). Hellinger distance measurements calculated from Kyoto Encyclopedia of Genes and Genomes (KEGG) Orthologues (KO) abundance showed that TB patient group (TBZ, TBW and TBM) are significantly different from one another than are HCs in terms of their functional repertoire (Fig. 1D); revealing dysbiosis at both phylogenetic and functional level of the GM in TB subjects.

TB-associated microbiota changes at the taxonomic level To determine the differential taxonomic abundance between HC and TBZ, Wilcoxon Rank-sum test was applied to the taxonomic profiles generated from 16S rRNA data set (see Supporting Information Additional File S6). At the phylum level, the majority of the sequences were distributed among three major phyla, Bacteroidetes, Table 2. PERMANOVA analyses of clinical variables with UniFrac distances. BMI

HDL

LDL

Total cholesterol Treatment Disease status

0.411 0.365 0.002 0.004

0.066

0.003

BMI, body mass index; HDL, high-density lipoprotein; LDL, lowdensity lipoprotein. p 0.05 indicated by bold font.


4 A. Maji et al. Table 3. Bacterial a–diversity indices from 16S amplicon sequencing data. Statistical comparison using Student’s t test (p values)

Mean 6 standard deviation Alpha diversity metrics

HC (n 5 5)

TBZ (n 5 5)

TBW (n 5 5)

TBM (n 5 5)

HC vs. TBZ

TBZ vs. TBW

TBZ vs. TBM

Observed species Shannon index Phylogenetic diversity

697.6 6 78.68 3.786 6 0.152 18.64 6 1.65

1038 6 94.45 4.66 6 0.33 20.54 6 1.33

1268 6 125.3 5.26 6 0.248 24.23 6 1.79

1074 6 109 4.93 6 0.322 20.89 6 1.88

0.0317 0.0317 0.22

0.22 0.158 0.22

0.801 0.66 0.309

HC, household contacts; TBZ, TB patients at day zero; TBW, TB patients after week treatment; TBM, TB patients after month treatment; p 0.05 indicated by bold font.

Firmicutes and Proteobacteria, which contributed more than 97% of total population in the Indian cohort. Although analysis showed no significant differences at the higher hierarchy levels (Fig. 2A inset), we observed significant increase in abundance of many short-chain fatty acid (SCFA) producers such as Faecalibacterium, Coprococcus, Phascolarctobacterium and Pseudobutyrivibrio in the TBZ group (TBZ > HC, false discovery rate (FDR) adjusted p 0.05, Fig. 2B, see Supporting Information Additional File S6). Also, of interest was the predominance of Prevotella in all the subjects, as recently observed in another Indian cohort (Bhute et al., 2016). However, Prevotella was significantly depleted in TB patients in this study. Moreover, we observed increased abundance of Bacteroides associated with TB patient group. Studies in the past have established an anti-correlation between Bacteroides and Prevotella in the GM (Ley, 2016), and it is suggested that Prevotella spp. are better equipped at consuming plantderived carbohydrates and absence of dietary fiber enable Bacteroides colonize the gut on decline of Prevotella populations (Koropatkin et al., 2012). Our results also show a similar anti-correlation trend and in view of the stark loss of appetite seen in TB patients, support the theory that reduction in overall food intake may allow selection of Bacteroides in TB patients capable of utilizing host mucus glycans more efficiently. To study the GM of subjects at a finer taxonomic resolution, the filtered Illumina-sequencing derived metagenomic reads were aligned to 4097 reference microbial genomes (see Supporting Information Additional File S7). To circumvent the inaccuracies that might arise due to non-taxa specific hits, we applied stringency criteria, wherein taxa were only included in further analysis if at least 10 reads hit the same taxon. In addition, we also utilized a conserved gene-based method MetaPhlAn2 (Truong et al., 2015), for taxonomic annotation of reads. Supporting Information Additional File S8 and Supporting Information Additional File S1: Figs. S2–S4 show the mean abundance of metagenomic families, genera and species obtained by MetaPhlAn2. Linear modelling analysis also revealed significant correlations of TBZ/HC ratio of metagenomic species obtained by both the methods (see Supporting

Information Additional File S1, Fig. S5). The taxa common between MetaPhlAn2 and read-mapping method were further used for analysis. The species composition of HC and TBZ subjects obtained by read-mapping method were compared by the Wilcoxon Rank-sum test. Supporting Information Additional File S9 lists the species-level bacterial taxa whose representation changes with disease and Fig. 2C represents most of the differentially abundant species in TB patients (FDR adjusted p 0.05). Similar to the 16S rRNA analysis, we observed a significant decrease in Prevotella and Bifidobacterium and increase in propionateproducer Phascolarctobacterium succinatutens in TBZ (Fig. 2C). Furthermore, a trend of increased mean abundance of butyrate producers Roseburia inulinivorans and Faecalibacterium prausnitzii was seen in TBZ samples. This confirmed the relevance of Prevotella and SCFA producers in discriminating the TB patient group from the HCs. Interestingly, we noted a significant abundance of M. tuberculosis in TBZ samples (Fig. 2C). Since the GI tract gets exposed to any pathogen or antigen that invades the respiratory system so obtaining mycobacterial species in faecal matter of TB patients was not very surprising. However, the presence of mycobacterial species in HC samples was startling. An increased presence of some other pathobionts such as Shigella sonnei, Escherichia coli, Streptococcus pneumonia and Streptococcus vestibularis was also seen in TBZ samples (see Supporting Information Additional File S9), signifying the increased fitness of these pathogens in dysbiotic GM. Metagenome-wide association study confirms enrichment of butyrate-producing bacteria in TB patients To overcome the limitations of taxonomic approaches and identify microbial species and functions independent of reference genomes, de novo assembly of metagenomic reads was performed (see Experimental Procedures). After filtering, a total of 638,069 non-redundant genes were obtained from the Indian gut samples. These data were incorporated in the Integrated gene catalogue (Li et al., 2014) to generate a Combined gene catalogue, which contributed 297,171 unique genes to the existing catalogue updating it to 10,199,488, as described in the


Gut microbiome of TB patients 5


6 A. Maji et al. Fig. 1. Gut microbiota of TB patients show altered bacterial diversity. A. PCA of unweighted UniFrac distances from 16S rRNA gene data set showed that TB patients (n 5 5, green colour) clustered separately from household contacts (n 5 5, red colour), p value 5 0.015. PC1 and PC2 represent the top two principal coordinates that captured most of the variance, with the fraction of variance captured by that coordinate shown in percent. B. Influence of 16S rRNA gene sequencing depth on rarefaction analysis. Rarefaction curves of different samples are shown, each with at least 27,000 V3–16S rRNA sequences per sample. Boxes represent the interquartile range between the first (25th percentile) and third quartiles (75th percentile) and the line within represents the median at each rarefaction depth for HC and TBZ (with n 5 5). C. Unweighted UniFrac distances between samples within a group calculated from the V3–16S rRNA data set are shown for HC, TBZ, TBW and TBM group samples. Mean values 6 SEM are plotted. D. Interpersonal variations were calculated for the whole genome shotgun sequence data using Hellinger distances derived from KO profile across the samples. Mean values 6 SEM are plotted. The statistical analyses across groups (in c and d) were performed using Student’s t test with multiple comparisons (Tukey’s method). FDR adjusted p values: *p 0.05, **p 0.01,***p 0.001.

experimental procedures section (see Supporting Information Additional File S1: Fig. S6). These genes were then grouped into metagenomic clusters (MGCs), based on metagenome-wide association study (MWAS), following the method employed by Qin et al. (2012). We only considered genes that were shared by at least six subjects and included 146,234 genes in our analysis. The abundance of the largest 195 MGCs with at least 50 genes in an MGC was calculated and the differentially abundant clusters (log odds ratio < 22, log odds ratio > 2) are listed in Supporting Information Additional File S10. We found a distinctive association of Prevotella-dominated MGCs with HCs; while, MGCs enriched in TB patients mostly included butyrate-producing bacteria, such as Eubacterium rectale, R. inulinivorans, F. prausnitzii apart from the enrichment of Bacteroides (Fig. 3). Metagenomic evidence for variation in microbial metabolic functions in TB patients To characterize the microbial functions, we annotated all the genes in our catalogue to the KO database. Differential abundance of KOs was calculated using an Odds Ratio test (Greenblum et al., 2012) as described in the experimental procedures section, to examine the association of metabolic functions in HC and TBZ (Odds Ratio > 1| Odds Ratio < 1, Fisher’s test p value < 0.05). Functional annotation of KO into metabolic pathways revealed significant enrichment of propionate and butyrate metabolism in TBZ, highlighting the fact that difference in species composition is reflected in the shift in microbial functions. The most notable pathways higher in TBZ samples were flagellar assembly, sulfur relay system, bacterial chemotaxis and ABC transporters, seen upregulated in other diseased conditions and commonly associated with virulence and pathogenicity (Deng et al., 2017). Metabolic functions involved in biosynthesis of amino acids and cofactors and vitamin metabolism were depleted in TB patients (Fig. 4A, see Supporting Information Additional File S11). Reporter feature algorithm (Oliveira et al., 2008) was also used in combination with KEGG metabolic network to identify pathways associated with TBZ and HC groups, which yielded

qualitatively similar results (see Supporting Information Additional File S12). Spearman’s rank correlation analysis between significant metagenomic species and enriched pathways showed that M. tuberculosis correlated positively with propionate metabolism while negatively correlated with vitamin B6 metabolism. Conversely, Prevotella correlated positively with vitamin metabolism and amino acid synthesis (see Supporting Information Additional File S1: Fig. S7 and Supporting Information Additional File S13). Since the central theme from our microbiome analysis appeared to be an overrepresentation of butyrate and propionate metabolism in TB patients, we show the enrichment of individual genes encoding enzymes in these KEGG pathways (Supporting Information Additional File S1: Fig. S8). The dominant source of butyrate and propionate production by gut bacteria is the acetyl coenzyme A pathway (Vital et al., 2014) and succinate pathway (Reichardt et al., 2014) respectively. Fig. 4B, C show the enrichment of enzymes involved in these pathways by odds ratio method, displaying clear association of enzymes involved in butyrate and propionate metabolism with TBZ samples. Effect of DOTS therapy on the gut microbiota of TB patients To explore the alteration in GM associated with the DOTS multi-drug regimen, comparisons were made using Wilcoxon Rank-sum test between TBZ and TBW/TBM samples (see Supporting Information Additional File S14). We did not observe a significant recovery in relative abundances of species altered due to disease state, apart from the significant decrease (FDR adjusted p 0.05) in abundance of mycobacterial pathogens after one month of treatment (Fig. 5A). As evident from analysis of mycobacterial species, the abundance of M. tuberculosis is comparable in TBZ and TBW samples even after a week of commencement of DOTS (Fig. 5B), and perhaps this is the reason for maximum perturbation observed in this sample group. These treatment time-points were inadequate to observe any long-term effects of DOTS therapy on the GM and longer follow-up studies are required to gain better insights into antibiotic-induced perturbations.



Fig. 2. Compositional differences in microbial population between TB patients and healthy controls. A. Taxonomic comparison of the abundant (mean proportion 1%) phyla (shown as inset) along with the differentially abundant families (FDR adjusted p value 0.05; Wilcoxon rank sum test) between HC and TBZ are shown as boxplots. B. Differential abundance of bacterial genera that differ between patients with TB (TBZ) and household contacts (HC), Wilcoxon Rank-sum test, FDR adjusted p 0.05. The boxes represent the IQR between the first (25th percentile) and third quartiles (75th percentile) and the line inside represents the median. Whiskers denote lowest and highest values within 1.5 times IQR from the first and third quartiles respectively. Circles denote the outliers. C. Scatter plot of log transformed median values of species abundance in HC (n 5 6) and TBZ (n 5 6) subjects. Red-blue colour gradient points represent the differentially abundant (FDR adjusted p 0.05; Wilcoxon Rank sum test) species between HC and TBZ with lower p values from red to blue, whereas grey points represent species not differentially abundant (FDR adjusted p value > 0.05, Wilcoxon rank-sum test). Phylum, family and genus level abundances were measured from taxonomic profiles generated from 16S rRNA data set and species abundance was measured using reference-genome-based mapping.


8 A. Maji et al.

Fig. 3. Identification of differentially abundant MGCs in TBZ and HCs. The enrichment of MGCs is shown with the log odds ratio of 195 MGCs plotted against –log10 (p values) on y axis. As shown in the figure, the MGCs with highly significant enrichment in either HC or TBZ are seen on the extremes of the x axis and y axis. The MGCs (in red) with log odds ratio > 2.0 are enriched in TBZ, whereas the MGCs (in green) with log odds ratio < 22 are enriched in HC.

Notably, the comparison of gut microbial functions in TBZ and TBM samples showed recovery in some microbial pathways in TBM such as significant reduction of flagellar assembly, bacterial chemotaxis and propanoate metabolism functions and enrichment of amino acid metabolism and vitamin metabolism (Fig. 5C and Supporting Information Additional File S15). Discussion This study describes previously unexplored human gut microbial composition of pulmonary TB patients. The study group consists of subjects from the Indian subcontinent, which is a relatively less explored metagenome pool and has not been functionally characterized yet (Shetty et al., 2013). We established an updated gene catalogue by adding data from Indian ethnicity that would provide additional information to existing reference set. Although attempts to study GM changes in TB have been made earlier, they have offered us limited insights, probably due to the design of the study (Dubourg et al., 2013; Winglee et al., 2014).

We observed microbial dysbiosis associated with active TB infection in patients. The bacterial diversity results indicated an increase in diversity in TB patients, even after one month of antibiotic administration, as compared to the HCs. This was not very surprising considering out of the four drugs in the DOTS regimen, three are specific to M. tuberculosis and thus may not strongly impact the other gut bacteria. Interestingly, previous studies have reported increased microbial diversity in sputum samples of TB patients than healthy participants (Cui et al., 2012; Krishna et al., 2016). These results suggest a dynamic response of the GM in face of any disease or antibiotic-induced dysbiosis. We observed the predominance of Prevotella in the study cohort disposing its categorization to Enterotype 2, one of the three stratifications of the GM (Arumugam et al., 2011), which corroborates with earlier reports in Indian population (Bhute et al., 2016). Although its prevalence was seen in all individuals, Prevotella was seen as a discriminatory taxon associated mainly with HCs in our cohort and possibly, the decrease in Prevotella abundance allows the increase in microbial diversity observed in TB group.



Fig. 4. Functional alteration associated with TB patients. A. A comparison of enriched KEGG pathways in TBZ and HC derived from metagenomic data using Odds ratio for enrichment and further Fisher’s exact test applied on proportion of enriched KOs belonging to each pathway. The log odds ratio (TBZ/HC) was calculated for significantly enriched pathways (with mean proportion 0.05) as identified using Fisher’s test and represented as bar plots, TBZ (red; OR 1.5) and HC (green; OR 0.75) samples. Asterisks depict significance of proportions of enriched KOs (Fisher’s test) observed for each pathway (p value 0.05, Fisher’s exact test; *p 0.05, **p 0.01, ***p 0.001). B. Genes encoding enzymes for butyrate metabolism that carries out metabolism of acetyl-CoA to butyrate are found enriched (Odds Ratio > 1.0) in gut metagenome of TB patients. Odds ratio of genes is shown on left and their corresponding pathway location is shown as flowchart on the right, with colours of the arrows showing enzyme catalysis reaction. C. Genes encoding enzymes for propionate metabolism that carries out metabolism of succinyl-CoA to propionate are found enriched (Odds Ratio > 1.0) in gut metagenome of TB patients. Odds ratio of genes is shown on left and corresponding pathway location is shown as a flowchart on the right, with colours of the arrows showing enzyme catalysis reaction.

The results display increased ratio of Firmicutes to Bacteroidetes in TB subjects, which is known to affect the concentration of SCFAs, GM-derived bacterial metabolites, viz. acetate, butyrate and propionate, recognized to regulate energy metabolism, cholesterol metabolism and

systemic inflammatory response in humans (den Besten et al., 2013). A prominent theme that emerged from the GM data of TB patients was the significant increase in butyrate and propionate-producing bacteria (F. prausnitzii, Roseburia, E. rectale, Butyrivibrio, Phascolarctobacterium)


10 A. Maji et al.

Fig. 5. Effect of DOTS therapy on metagenome profile of TB patients. A. Scatter plot of median species abundance in TBZ (n 5 6) and TBM (n 5 6) subjects. Colour gradient points (red-blue) represent differentially abundant (FDR adjusted p 0.05; Wilcoxon Rank sum test) species between groups TBZ and TBM with lower p values from red to blue, whereas grey points represent species which are not differentially abundant (FDR adjusted p > 0.05, Wilcoxon Rank-sum test). B. The bar plot shows the relative abundance (mean 6 SEM) of abundant mycobacterial species (mean relative abundance > 0.001%) in all samples (HC, TBZ, TBW and TBM). C. A comparison of enriched KEGG pathways in TBZ and TBM using Odds ratio for enrichment and further Fisher’s exact test applied on proportion of enriched KOs belonging to each pathway. The log odds ratio (TBZ/TBM) was calculated for significantly enriched pathways (with mean proportion 0.05) through Fisher’s test and represented as bar plots, TBZ (red; OR 1.5) and TBM (green; OR 0.75). Asterisks depict significance of proportions of enriched KOs (Fisher’s test) observed for each pathway (p value 0.05, Fisher’s exact test; *p 0.05, **p 0.01, ***p 0.001).

and their related metabolic functions (butyrate and propionate metabolism). This was intriguing since F. prausnitzii is usually associated with a healthy colon and its reduction is

often linked to diseased condition (Miquel et al., 2013). To our knowledge, only one study has emphatically challenged the protective role for F. prausnitzii in all disease


Gut microbiome of TB patients conditions and observed its increase in Crohn’s disease patients (Hansen et al., 2012). In view of our results, we also believe that upsurge of SCFA producers may not be beneficial in all health conditions and more elaborate studies are required to better understand their dynamic role in diverse dysbiotic conditions such as obesity and TB. Butyrate impacts human health via its role as an antiinflammatory agent, primarily by inhibition of nuclear factor-jB activation and interferon gamma (IFN-g) signalling (Segain et al., 2000; Canani et al., 2011). The nuclear factor-jB signalling pathway is central to body’s immune response to many pathogens and is involved in transcriptional regulation of many cytokine genes, including tumor necrosis factor-alpha (TNF-a). It is known that IFN-g and TNF-a play a critical role in the control of M. tuberculosis infection and is required for formation and maintenance of granuloma. Notably, subjects with latent TB on receiving anti-TNF therapy for treatment of rheumatoid arthritis often show increased rate of reactivation of active TB. Also, TNF neutralization during a study on macaques resulted in disseminated disease (O’Garra et al., 2013). By an independent mechanism, butyrate is also known to enhance regulatory T-cells (Treg) population in the gut (Smith et al., 2013). Although Treg limits collateral tissue damage caused by vigorous antimicrobial immune responses, it suppresses the protective proinflammatory T-cell response among TB patients via interleukin-10 (IL-10) and facilitates chronic infection. Interestingly, butyrate has been implicated as a mediator of altered immune response to M. tuberculosis in diabetes mellitus patients, which may influence susceptibility to TB (Lachmandas et al., 2016). It has also been observed that Treg increases in blood or at the site of infection in patients with active TB and significant decline in the frequency of various Treg subsets, along with an elevation of IFNg production, has been seen on completion of anti-TB therapy (Guyot-Revol et al., 2006; Chen et al., 2007; Singh et al., 2012). It is suggested that various M. tuberculosis strains may induce increased production of IL-10 as a way to evade immune evasion (Redford et al., 2011). Thus, a sharp increase in butyrate levels may have deleterious effects on host immune response to M. tuberculosis infection. This highlights the fact that while reduction in butyrate-producing bacteria might be unfavourable in gastrointestinal inflammatory disorders like inflammatory bowel disease, an upsurge can also prove detrimental to the host in response to bacterial infectious diseases. We postulate another way by which GM impacts host immunity in TB is through its regulation on the bioavailability of nutrients such as vitamins and amino acids to host (Gill et al., 2006; Morowitz et al., 2011; Brestoff and Artis, 2013). The microbial function analysis revealed a significant decline in processes such as biosynthesis of amino

11

acids and vitamin metabolism in TB patients. This depletion may impact circulating pool of essential amino acids and vitamins and since vitamin derivatives play an important role in the regulation of immune system via activation of mucosa-associated invariant T-cells (Brestoff and Artis, 2013), it reaffirms involvement of GM in TB disease process. Depletion of amino acids may also impact its contribution to host amino acid metabolism. In the absence of dietary supplementation, depletion of amino acids can manifest as protein-energy malnutrition, an important risk factor for TB, which has been shown in animal models for TB to confer devastating effects on the ability of the host to restrict M. tuberculosis proliferation (Cegielski and McMurray, 2004). Significant improvement in recovery of patients with kwashiorkor, a protein-energy malnutrition disease, is reported with the addition of oral antibiotics to the therapeutic food regimens (Trehan et al., 2013), implicating function of GM in the manifestation of proteinenergy malnutrition (Brestoff and Artis, 2013). Thus, it seems logical to provide increased dietary protein to TB patients to help in recovery process. Alteration of appetite-regulatory gut hormones has been implicated in wasting associated with TB (Chang et al., 2013). Butyrate and propionate are known to regulate the release of these gut hormones to protect against dietbased obesity (Lin et al., 2012). In an epidemiological study, low BMI was positively associated with risk of death from respiratory diseases in Indians/Asians (Zheng et al., 2011). In light of these observations from previous studies, we speculate that higher susceptibility for TB in low BMI subjects may be in part related to excess butyrate and propionate production that affects patient appetite and weight and suppresses the immune system. However, at the moment we cannot distinguish between the two possible conditions: whether TB changes the GM or the GM changes lead to increased risk of TB. Furthermore, another characteristic feature in TB patients, as also revealed in our results, is the low blood cholesterol levels, which may get regulated by butyrate and propionate by inhibiting intestinal cholesterol biosynthesis pathway (Arora et al., 2011; Canani et al., 2011). Clinical evidence suggests that cholesterol-rich diet accelerates recovery in TB patients (Perez-Guzman et al., 2005). In light of these observations, we must consider this hypothesis that while decrease in butyrate-producing bacteria might be harmful in metabolic disorders that accompany hypercholesterolemia (Canani et al., 2011), in the case of infectious diseases like TB, an upsurge of butyrate and propionate-producing bacteria might have grave consequences on host lipid homeostatic mechanisms. Here, an intriguing aspect that deserves attention is the apparent paradox whereby increase in serum levels of many acute-phase proteins such as C-reactive protein is observed in TB patients (Chegou et al., 2016); but, our


12 A. Maji et al. findings reveal GM-influenced anti-inflammatory conditions, which allow increased bacterial diversity in TB patients, including many opportunistic pathogens. Notably, retroviruses are known to employ the GM to elicit IL-10 for immune subversion and utilize it for its transmission (Kane et al., 2011). We speculate that M. tuberculosis might also utilize GM for its persistence, but further studies need to be made to explore this possibility. Although the current study reports GM changes associated with TB in a small cohort size, our results emphasize the importance of gut-lung axis and provide many promising leads that deserve further exploration. We must mention that analysis of the microbiome data revealed reads of mycobacterial species, including M. tuberculosis, in HCs, so although they were clinically ‘healthy’ individuals but exposed to M. tuberculosis. Nonetheless, we followed up (3 years after recruitment) on the health status of the participants and confirmed the absence of symptomatic features of TB in the control group since the time of study. Moreover, the follow-up also revealed that all the patients fully recovered after TB treatment and have not shown any relapse till now. We speculate that the HCs could be a reservoir for latent infection; and the microbial dysbiosis might be an accessory factor contributing towards the switch to the diseased state. This hypothesis gains support by a recent study wherein broad-spectrum antibiotic-driven microbial dysbiosis provokes susceptibility to TB in experimental TB mice model set-up, through an elevation in the number of Tregs (Khan et al., 2016). It has also been suggested in an earlier study that frequent oral exposure to environmental non-TB mycobacteria induces populations of mycobacteria-specific Treg cells in the gut, which potently restrict protective immune responses during TB (Shafiani et al., 2010). In light of these results, the significance of high numbers of M. smegmatis population in HC subjects merits further clinical investigation. The results from this study provide a novel framework for understanding the gut microbial changes accompanying M. tuberculosis exposed Indian population cohort. Detection of M. tuberculosis from faecal samples and its sharp decline after a month’s treatment in patients reaffirms the potential of using faecal samples for non-invasive diagnostic and monitoring purposes as suggested earlier (El Khechine et al., 2009). We hypothesize that the dysbiosis of the GM associated with TB, especially enrichment of propionate and butyrate-producing bacteria, contributes to impairment of immune response. This may also, in part, be responsible for appetite and weight loss as well as reduced cholesterol levels, accompanied with the disease condition. Our findings indicated recovery of some of the vital functional pathways in TBM group, which is of interest, given clinical improvement usually seen in most patients within 1–2 months of drug administration, including subjects in this study. Larger prospective follow-up studies would be

required to validate the insights obtained from the present findings. This knowledge will help us better appreciate the gut-lung axis and can have potential implications for development of better management strategies such as vaccines, nutrient supplements for controlling TB. Conclusions This study characterized for the first time the GM of pulmonary TB patients. The study group consisted of Indian subjects, which has not yet been functionally characterized. Significant taxonomic and functional alterations of GM were observed between HCs and TB patients. Apart from drastic reduction of M. tuberculosis population, no other major taxonomic recovery was observed till one month of DOTS treatment. Data analyses revealed enrichment of butyrate and propionate metabolism in TB patients, which given their role in host immunity and nutrient metabolism, strongly implicates GM in the pathophysiology of TB. Experimental procedures Subject enrollment and sample collection All the subjects were recruited for the study at the Rajan Babu Institute of Pulmonary Medicine and Tuberculosis (RBIPMT), Delhi, India between March and September 2013. Written informed consent was obtained from all the subjects. The institutional ethics committees of CSIR-Institute of Genomics and Integrative Biology (IGIB), Delhi and RBIPMT approved the recruitment criteria and study protocol. All the patients were prospectively recruited on the basis of TB symptomatic features. Following confirmation by sputum-smear microscopy, the patients were put on the DOTS regimen, as per the World Health Organization recommended strategy implemented by Revised National Tuberculosis Control Programme in India, where the patients receive antibiotics isoniazid, rifampicin, ethambutol and pyrazinamide. After the diagnosis was confirmed, faecal samples from TB subjects (n 5 6) were collected in sterile containers at three time points: before the start of treatment (TBZ), a week after treatment (TBW) and a month after treatment (TBM). One healthy blood relative of each patient was recruited for the control group (n 5 6). The faecal samples of the control group were only collected once (HC), corresponding to day zero of TB subjects. The collected faecal samples were frozen on procurement (within 5 h of sample collection by subjects). The weight and height measurements were done for all subjects at the start of the study, day zero. Also, blood samples were collected from all subjects on day zero, and from patients on day 7 and day 30 in vacutainer tubes (K2-EDTA). Serum was separated and stored at 2808C until further analysis. We had a brief telephonic contact with the recruited subjects in May 2016 to discuss recovery and presence of any symptomatic features in any of the subjects.

Analysis of serum lipid levels Total cholesterol, LDL and HDL levels were measured with commercially available kits from Wako Diagnostics and Cobas


Gut microbiome of TB patients

13

Integra automated laboratory workstation (Roche Diagnostics), following manufacturer’s instructions.

measures: observed species, Shannon index (Shannon, 1997) and phylogenetic diversity (Faith, 1992).

DNA extraction and sequencing

Whole genome shotgun sequencing

Total DNA from all faecal samples (250 mg) were extracted using QIAamp DNA Stool Mini Kit (Qiagen, CA) following manufacturer’s protocol with the inclusion of a mechanical disruption step with 0.1 mm zirconium beads using a beadbeater (BioSpec Products). DNA concentration and quality were estimated using a Nanodrop instrument (Thermo Scientific) and agarose gel electrophoresis. Illumina NextSeq500 sequencer (Illumina, CA) was used at IISER, Bhopal for 16S rRNA gene sequencing and the faecal whole genome shotgun sequencing was performed at Beijing Genomics Institute, China using Illumina HiSeq 2000 instrument. A detailed description of sequencing procedure is mentioned in Supporting Information (Supporting Information Additional File S1).

Quality filtering of reads. The reads with ambiguous bases were filtered out and the resulting reads were further quality filtered (reads having 100% bases above Q30) using NGSQC Toolkit (Patel and Jain, 2012).

Data analysis 16S rRNA analysis Variable region 3 (V3) of 16S rRNA amplified from 24 faecal samples resulted into a total of 16,168,634 high-quality paired-end reads {n 5 673,681 6 146,659 (mean 6 SD) reads per sample}. The paired-end reads were assembled into single end reads using FLASH version 1.2.11 (Magoc and Salzberg, 2011). The reads were further quality filtered using NGS QC Toolkit (Patel and Jain, 2012) with a stringent criteria of 90% bases above Q30 to select the reads for further analysis. The 16S rRNA sequences were clustered into species-level phylotypes using Closed reference OTU picking protocol (Edgar, 2010) against Greengenes database version 13_5 (DeSantis et al., 2006) at 97% identity using Quantitative Insights into Microbial Ecology suite of software tools version 1.8 (Caporaso et al., 2010). Of the total reads, 98.7% reads (15,961,140) matched a reference sequence in Greengenes database at 97% nucleotide sequence identity. The most abundant read from each OTU was selected as representative read for that OTU. The taxonomy associated with Greengenes database to which OTUs matched was assigned to OTUs. The OTU table was filtered to remove OTUs containing < 10 counts in all samples. A total of 3022 OTUs were obtained using closed reference OTU picking protocol. Biodiversity analysis The intrasample diversity (alpha) and inter-sample diversity (beta) using unweighted UniFrac distances were performed, as described previously (Lozupone and Knight, 2005; Qin et al., 2012; Yatsunenko et al., 2012), for rarefied OTU tables at a depth of 270,000 sequences/sample. The alpha diversity was calculated using three different

Phylogenetic assignment of reads. A total of 4097 reference microbial genomes were obtained from Human Microbiome Project (HMP) and National Centre for Biotechnology Information (NCBI) RefSeq databases on 5 December 2015. The databases were independently indexed into two Bowtie indexes using Bowtie-2 (Langmead and Salzberg, 2012). The metagenomic reads were aligned to reference microbial genomes using Bowtie-2. The mapped reads from both indexes were merged by selecting the alignment having higher percent identity (with at least 90% identity). The percent identity was calculated using the formula: % identity 5 100*(matches/total aligned length). The reference microbial genome abundance was calculated by summing the total number of reads aligned to the genome, normalized by genome length and the total number of reads in each dataset. For reads having hits to both indexes with equal identity, each genome was assigned 0.5 read count. The relative abundance of each genome was calculated by adding the normalized abundance of each genome divided by total abundance. A taxon was considered to be present in the metagenomic samples only if at least 10 reads mapped to the same taxon. In an independent approach, MetaPhlAn2 (Truong et al., 2015) was used using default parameters to perform taxonomic assignment of reads, which is a conserved gene marker based assignment method. Only those taxa that were predicted by MetaPhlAn2 were considered for the calculation in the abundance tables generated from the read-mapping method and further included in analysis. By this combined approach, 151 species were predicted, which had at least 10 read hits mapped to them and were common with the MetaPhlAn2 analysis. Updated metagenomic gene catalogue The de novo assembly of high-quality paired-end reads from each of 24 samples was performed using MetaVelvet (Namiki et al., 2012) with a kmer size of 77 bp. The genes from assembled contigs were predicted using MetaGeneMark (v2.8) (Zhu et al., 2010). The genes with read lengths < 100 bp were filtered out resulting in a total of 1,785,638 genes predicted from 24 samples. All predicted genes were aligned pairwise using BLAT and the genes


14 A. Maji et al. aligned at 95% identity and 90% alignment coverage were clustered together resulting into a non-redundant gene catalogue with a total of 638,069 non-redundant genes. The longest read in the gene cluster was selected as a representative in unique gene catalogue. To increase the coverage of the total gene cohort found in the human GM, the gene catalogue derived from 24 faecal samples of Indian subjects were combined with the Integrated gene catalogue (Li et al., 2014) consisting of 9,879,896 genes from 1267 gut metagenomes (1070 individuals) representing the European, Chinese and US populations to generate a Combined gene catalogue. Our gene catalogue of 24 Indian samples was combined with integrated gene catalogue by pairwise alignment using BLAT in a similar manner as mentioned above resulting into an updated gene catalogue of 10,199,488 genes.

sequences for KO IDs from KEGG server. Each amino acid sequence was assigned to eggNOG orthologous group and KO orthologous group based on the highest scoring hit containing at least one HSP above 60 bits and E value 1026. The relative abundance of each KO was calculated by adding up the abundance of genes mapping to the same KO ID, which was then normalized to calculate the relative abundance for each sample. A similar approach was used to calculate the relative abundance of eggNOGs in each sample. A total of 6567 KO orthologous groups and 50,294 eggNOG orthologous groups were identified from the repertoire of combined genes. The Hellinger distances were calculated between the samples using normalized values of eggNOG relative abundance and KO relative abundance in each sample. Effect of covariates

Quantification of gene content For quantification of gene content, we adopted the strategy used by Qin et al (Qin et al., 2012). The method was as follows: The high-quality reads were aligned against the generated combined gene catalogue using SOAP2 (Li et al., 2009) with a criterion of 90% identity match in SOAP2 aligner. In our sequence based profiling, the two types of alignments were considered: (i) The entire of paired-end read mapped to the gene. (ii) One end of paired-end read mapped to a gene and other end remained unmapped. In both cases the mapped read was counted as one copy. Further, the read count was normalized based on length of the gene as: bi 5 xi/Li The relative abundance of a gene within the sample was calculated as: x

i bi L ai 5 P 5P i xj j bj j L

j

ai: relative abundance of gene i in sample S; xi: The times which gene i can be detected in sample S (the number of mapped reads); Li: length of gene i; bi: copy number of gene i in sequenced data from sample S. Rjbj 5 Sum total of all gene abundances from sample S.

Functional annotation of genes The genes from combined gene catalogue were translated into amino acids sequences and were aligned against eggNOG v4.0 (Powell et al., 2014) and KEGG v60 (Kanehisa et al., 2012) databases using BLASTP. The KEGG v60 (year 2011) version was updated by retrieving new

To assess the effect of different covariates such as age, BMI, disease status and diet on 16S rRNA gene and whole genome shotgun sequence profiles, Pearson correlation coefficients of these covariates were calculated with factor loadings of each factor. The UniFrac Distance matrix was used for 16S rRNA taxonomic profiles, and Hellinger distance matrix was used to determine the effect of covariates on metagenomic profiles. PERMANOVA (McArdle and Anderson, 2001) was also performed to assess the association of the clinical variables and unweighted UniFrac distance matrices based on 16S rRNA data. The analysis was performed using the method implemented in R package-‘vegan’ (Dixon, 2003) with 10,000 permutations. Construction of MGCs for MWAS Using the gene profiles generated from the metagenomic read-mapping, we identified highly correlated genes by calculating pairwise Pearson’s correlation coefficient between each of the gene pairs. The genes present at 0.0001% relative abundance in at least 6 or more samples were used for the pairwise correlations to reduce noise from low abundant genes. The correlation coefficients of genes across subjects were calculated pairwise with FDR adjusted p values for multiple comparisons to control for FDRs. Clustering of genes was performed using Markov Clustering Algorithm implemented in MCL software (Enright et al., 2002) which uses flow simulation in the graph to identify robust clusters. Edges were created between the genes that had correlation coefficients 0.85. The clusters were generated using inflation parameters of 1.2, 2.0, 3.0, 4.0, 5.0 and 6.0. The inflation factor of 2.0 was finally chosen for further analysis as the number of clusters did not show much variation with higher inflation factor. The relative abundance of MGC in each sample


Gut microbiome of TB patients was estimated using relative abundance values of all genes from that MGC. The genes with top 1% and lowest 1% relative abundance were removed and Poisson distribution was fitted to the relative abundance values of the rest of genes. The mean estimated from Poisson distribution was assigned as the relative abundance of that MGC. The profile of MGCs was generated and used for further analysis. Taxonomic assignment of MGCs To determine the taxonomic origin of each MGC, all the genes were aligned against 4097 reference microbial genomes from HMP and NCBI using BLASTN. The alignment hits were filtered using E value 1026 and alignment coverage 80% of the gene length. In the case of multiple best hits with equal identity and scores, Lowest Common Ancestor method was used for taxonomic assignment. The genes were finally assigned taxa based on comprehensive parameters of sequence similarity across phylogenetic ranks as described in METAHIT (Qin et al., 2010). An identity threshold of 95% was used for assignment up to species level, 85% identity threshold for genus level and 65% identity for phylum level. The taxonomic assignment of MGCs was performed using the criteria that 50% of the total genes in each MGC should map to the same phylogenetic group. Statistical analyses All statistical analyses were carried out using the R packages (Team, 2014). Differential abundance of species was tested by Wilcoxon Rank sum test. When multiple hypotheses were considered simultaneously, p values were adjusted to control for the FDR with the method described previously (Benjamini et al., 2001). Enzyme enrichment To identify the enzymes associated with a host state, we used Odds Ratio as a measure of the enrichment of an enzyme in a host. The Odds Ratio were calculated as: OR P P P P P (k) 5 [ s5TBZ Ask/ s5TBZ( i6¼k Asi)]/[ s5HC Ask/ s5HC P ( i6¼k Asi)], where Ask denotes abundance of enzyme k in sample S. The odds ratio value > 1 were considered higher in TBZ whereas those < 1 were considered higher in HC. The percentage of KOs found to be enriched in TBZ or HC in each pathway were calculated and the Fisher’s exact test was used to determine the significance of the pathway as differentially abundant. In this analysis, only those pathways were considered which have at least 70% KO (out of all KOs) mapped to that pathway. Apart from that, Reporter features algorithm was used for gene-set analysis of significant pathways associated with different groups of

15

samples. The algorithm takes the adjusted p values and fold changes (in this case log odds ratio) as input for each KO. The gene statistic is calculated based on the significant association of KO and its direction of change through which the pathway is scored by calculating the global p value (Oliveira et al., 2008). Enrichment of SCFA producing enzymes All the enzymes in pathways of propionate and butyrate metabolism were mapped from substrate to final product and Odds ratio of all the enzymes in the pathway were calculated (as mentioned in the last section). Acknowledgements The whole genome shotgun sequencing was outsourced at BGI, China. The 16S rRNA sequencing work was carried out at the NGS facility, IISER Bhopal. The bioinformatic computations were performed on resources provided by the computing facilities at IISER, Bhopal and CSIR-IGIB, Delhi. We gratefully acknowledge all the subjects who participated in this study. This work was supported by the research grant from University of Delhi, J. C. Bose research grant (SERB), CSIR Task Force project BSC0123 and intramural funding received from IISER Bhopal. None of the funding organizations had a role in the study design, sample collection, data analysis or write up.

Availability of data and material DNA sequences have been deposited to SRA database (https://www.ncbi.nlm.nih.gov/sra) under SRA accession: SRP079687 and will be released at the time of publication. Conflict of interests All authors declare that they have no conflict of interests. Authors’ Contributions YS and RM conceived the project. RL, JPK, VKS and YS managed the project. AC helped in clinical sample procurement. HS generated data for clinical parameters. AM performed sample collection, extracted DNA and helped in 16S rRNA gene sequencing run. RS and PM carried out 16S rRNA sequencing data generation and preliminary analysis. VG, NKM, NT, ES, RV, MG and YH helped in data analysis in initial stages. DBD performed computational analysis and helped in manuscript writing. RM and DBD interpreted the data and RM drafted the manuscript. AS, VKS, GA, AA helped in discussion and improving the manuscript. All authors read and approved the final manuscript. AM, RM and DBD contributed equally to the study. References Akpovi, C.D., Gbaguidi, S.L.H., Anago, E., Affolabi, D., Dougnon, V.T., Faihun, F., and Anagonou, S. (2013)


16 A. Maji et al. Tuberculosis treatment raises total cholesterol level and restores high density lipoprotein cholesterol (HDL-C) in patients with pulmonary tuberculosis. Afr J Biotechnol 12: 6019–6024. Arora, T., Sharma, R., and Frost, G. (2011) Propionate. Anti-obesity and satiety enhancing factor? Appetite 56: 511–515. Arumugam, M., Raes, J., Pelletier, E., Le Paslier, D., Yamada, T., and Mende, D.R. (2011) Enterotypes of the human gut microbiome. Nature 473: 174–180. Benjamini, Y., Drai, D., Elmer, G., Kafkafi, N., and Golani, I. (2001) Controlling the false discovery rate in behavior genetics research. Behav Brain Res 125: 279–284. Bhute, S., Pande, P., Shetty, S.A., Shelar, R., Mane, S., Kumbhare, S.V., et al. (2016) Molecular characterization and meta-analysis of gut microbial communities illustrate enrichment of Prevotella and Megasphaera in Indian subjects. Front Microbiol 7: 660. Brestoff, J.R., and Artis, D. (2013) Commensal bacteria at the interface of host metabolism and the immune system. Nat Immunol 14: 676–684. Budden, K.F., Gellatly, S.L., Wood, D.L., Cooper, M.A., Morrison, M., Hugenholtz, P., and Hansbro, P.M. (2017) Emerging pathogenic links between microbiota and the gut-lung axis. Nat Rev Microbiol 15: 55–63. Canani, R.B., Costanzo, M.D., Leone, L., Pedata, M., Meli, R., and Calignano, A. (2011) Potential beneficial effects of butyrate in intestinal and extraintestinal diseases. World J Gastroenterol 17: 1519–1528. Caporaso, J.G., Kuczynski, J., Stombaugh, J., Bittinger, K., Bushman, F.D., Costello, E.K., et al. (2010) QIIME allows analysis of high-throughput community sequencing data. Nat Methods 7: 335–336. Cegielski, J.P., and McMurray, D.N. (2004) The relationship between malnutrition and tuberculosis: evidence from studies in humans and experimental animals. Int J Tuberc Lung Dis 8: 286–298. Chang, S.W., Pan, W.S., Lozano Beltran, D., Oleyda Baldelomar, L., Solano, M.A., Tuero, I., et al. (2013) Gut hormones, appetite suppression and cachexia in patients with pulmonary TB. PLoS One 8: e54564. Chegou, N.N., Sutherland, J.S., Malherbe, S., Crampin, A.C., Corstjens, P.L., Geluk, A., et al. (2016) Diagnostic performance of a seven-marker serum protein biosignature for the diagnosis of active TB disease in African primary healthcare clinic attendees with signs and symptoms suggestive of TB. Thorax 71: 785–794. Chen, X., Zhou, B., Li, M., Deng, Q., Wu, X., Le, X., et al. (2007) CD4(1)CD25(1)FoxP3(1) regulatory T cells suppress Mycobacterium tuberculosis immunity in patients with active disease. Clin Immunol 123: 50–59. Cui, Z., Zhou, Y., Li, H., Zhang, Y., Zhang, S., Tang, S., and Guo, X. (2012) Complex sputum microbial composition in patients with pulmonary tuberculosis. BMC Microbiol 12: 276. den Besten, G., van Eunen, K., Groen, A.K., Venema, K., Reijngoud, D.J., and Bakker, B.M. (2013) The role of short-chain fatty acids in the interplay between diet, gut microbiota, and host energy metabolism. J Lipid Res 54: 2325–2340.

ski, S.P., Jarek, M., Bhuju, S., and Deng, Z.-L., Szafran € bler, I. (2017) Dysbiosis in chronic periodontitis: Wagner-Do key microbial players and interactions with the human host. Sci Rep 7: 3703. Deniz, O., Gumus, S., Yaman, H., Ciftci, F., Ors, F., Cakir, E., et al. (2007) Serum total cholesterol, HDL-C and LDL-C concentrations significantly correlate with the radiological extent of disease and the degree of smear positivity in patients with pulmonary tuberculosis. Clin Biochem 40: 162–166. DeSantis, T.Z., Hugenholtz, P., Larsen, N., Rojas, M., Brodie, E.L., Keller, K., et al. (2006) Greengenes, a chimerachecked 16S rRNA gene database and workbench compatible with ARB. Appl Environ Microbiol 72: 5069–5072. Dixon, P. (2003) VEGAN, a package of R functions for community ecology. J Veg Sci 14: 927–930. Dubourg, G., Lagier, J.C., Armougom, F., Robert, C., Hamad, I., Brouqui, P., and Raoult, D. (2013) The gut microbiota of a patient with resistant tuberculosis is more comprehensively studied by culturomics than by metagenomics. Eur J Clin Microbiol Infect Dis 32: 637–645. Dye, C., Lonnroth, K., Jaramillo, E., Williams, B.G., and Raviglione, M. (2009) Trends in tuberculosis incidence and their determinants in 134 countries. Bull World Health Organ 87: 683–691. Edgar, R.C. (2010) Search and clustering orders of magnitude faster than BLAST. Bioinformatics 26: 2460–2461. El Khechine, A., Henry, M., Raoult, D., and Drancourt, M. (2009) Detection of Mycobacterium tuberculosis complex organisms in the stools of patients with pulmonary tuberculosis. Microbiology 155: 2384–2389. Enright, A.J., Van Dongen, S., and Ouzounis, C.A. (2002) An efficient algorithm for large-scale detection of protein families. Nucleic Acids Res 30: 1575–1584. Faith, D.P. (1992) Conservation evaluation and phylogenetic diversity. Biological Conservation 61: 1–10. Gill, S.R., Pop, M., Deboy, R.T., Eckburg, P.B., Turnbaugh, P.J., Samuel, B.S., et al. (2006) Metagenomic analysis of the human distal gut microbiome. Science 312: 1355– 1359. Greenblum, S., Turnbaugh, P.J., and Borenstein, E. (2012) Metagenomic systems biology of the human gut microbiome reveals topological shifts associated with obesity and inflammatory bowel disease. Proc Natl Acad Sci U S A 109: 594–599. Guyot-Revol, V., Innes, J.A., Hackforth, S., Hinks, T., and Lalvani, A. (2006) Regulatory T cells are expanded in blood and disease sites in patients with tuberculosis. Am J Respir Crit Care Med 173: 803–810. Hansen, R., Russell, R.K., Reiff, C., Louis, P., McIntosh, F., Berry, S.H., et al. (2012) Microbiota of de-novo pediatric IBD: increased Faecalibacterium prausnitzii and reduced bacterial diversity in Crohn’s but not in ulcerative colitis. Am J Gastroenterol 107: 1913–1922. Kane, M., Case, L.K., Kopaskie, K., Kozlova, A., MacDearmid, C., Chervonsky, A.V., and Golovkina, T.V. (2011) Successful transmission of a retrovirus depends on the commensal microbiota. Science 334: 245–249. Kanehisa, M., Goto, S., Sato, Y., Furumichi, M., and Tanabe, M. (2012) KEGG for integration and interpretation of large-


Gut microbiome of TB patients scale molecular data sets. Nucleic Acids Res 40: D109– D114. Karlsson, F.H., Tremaroli, V., Nookaew, I., Bergstrom, G., Behre, C.J., Fagerberg, B., et al. (2013) Gut metagenome in European women with normal, impaired and diabetic glucose control. Nature 498: 99–103. Khan, N., Vidyarthi, A., Nadeem, S., Negi, S., Nair, G., and Agrewala, J.N. (2016) Alteration in the gut microbiota provokes susceptibility to tuberculosis. Front Immunol 7: 529. Koropatkin, N.M., Cameron, E.A., and Martens, E.C. (2012) How glycan metabolism shapes the human gut microbiota. Nat Rev Microbiol 10: 323–335. Krishna, P., Jain, A., and Bisen, P.S. (2016) Microbiome diversity in the sputum of patients with pulmonary tuberculosis. Eur J Clin Microbiol Infect Dis 35: 1205–1210. Lachmandas, E., van den Heuvel, C.N., Damen, M.S., Cleophas, M.C., Netea, M.G., and van Crevel, R. (2016) Diabetes mellitus and increased tuberculosis susceptibility: the role of short-chain fatty acids. J Diabetes Res 2016: 6014631. Langmead, B., and Salzberg, S.L. (2012) Fast gapped-read alignment with Bowtie 2. Nat Methods 9: 357–359. Ley, R.E. (2016) Gut microbiota in 2015: Prevotella in the gut: choose carefully. Nat Rev Gastroenterol Hepatol 13: 69–70. Li, J., Jia, H., Cai, X., Zhong, H., Feng, Q., Sunagawa, S., et al. (2014) An integrated catalog of reference genes in the human gut microbiome. Nat Biotechnol 32: 834–841. Li, R., Yu, C., Li, Y., Lam, T.W., Yiu, S.M., Kristiansen, K., and Wang, J. (2009) SOAP2: an improved ultrafast tool for short read alignment. Bioinformatics 25: 1966–1967. Lin, H.V., Frassetto, A., Kowalik, E.J., Jr., Nawrocki, A.R., Lu, M.M., Kosinski, J.R., et al. (2012) Butyrate and propionate protect against diet-induced obesity and regulate gut hormones via free fatty acid receptor 3-independent mechanisms. PLoS One 7: e35240. Lonnroth, K., Williams, B.G., Cegielski, P., and Dye, C. (2010) A consistent log-linear relationship between tuberculosis incidence and body mass index. Int J Epidemiol 39: 149–155. Lozupone, C., and Knight, R. (2005) UniFrac: a new phylogenetic method for comparing microbial communities. Appl Environ Microbiol 71: 8228–8235. Magoc, T., and Salzberg, S.L. (2011) FLASH: fast length adjustment of short reads to improve genome assemblies. Bioinformatics 27: 2957–2963. Maji, A., Misra, R., Kumar Mondal, A., Kumar, D., Bajaj, D., Singhal, A., et al. (2015) Expression profiling of lymph nodes in tuberculosis patients reveal inflammatory milieu at site of infection. Sci Rep 5: 15214. Marsland, B.J., Trompette, A., and Gollwitzer, E.S. (2015) The Gut-Lung Axis in Respiratory Disease. Ann Am Thorac Soc 12(Suppl 2): S150–S156. McArdle, B.H., and Anderson, M.J. (2001) Fitting multivariate models to community data: a comment on distance-based redundancy analysis. Ecology 82: 290–297. Miquel, S., Martin, R., Rossi, O., Bermudez-Humaran, L.G., Chatel, J.M., Sokol, H., et al. (2013) Faecalibacterium prausnitzii and human intestinal health. Curr Opin Microbiol 16: 255–261. Morowitz, M.J., Carlisle, E.M., and Alverdy, J.C. (2011) Contributions of intestinal bacteria to nutrition and metabolism in the critically ill. Surg Clin North Am 91: 771–785.

17

Namiki, T., Hachiya, T., Tanaka, H., and Sakakibara, Y. (2012) MetaVelvet: an extension of Velvet assembler to de novo metagenome assembly from short sequence reads. Nucleic Acids Res 40: e155. O’Garra, A., Redford, P.S., McNab, F.W., Bloom, C.I., Wilkinson, R.J., and Berry, M.P. (2013) The immune response in tuberculosis. Annu Rev Immunol 31: 475–527. Oliveira, A.P., Patil, K.R., and Nielsen, J. (2008) Architecture of transcriptional regulatory circuits is knitted over the topology of bio-molecular interaction networks. BMC Syst Biol 2: 17. Patel, R.K., and Jain, M. (2012) NGS QC Toolkit: a toolkit for quality control of next generation sequencing data. PLoS One 7: e30619. Perez-Guzman, C., Vargas, M.H., Quinonez, F., Bazavilvazo, N., and Aguilar, A. (2005) A cholesterol-rich diet accelerates bacteriologic sterilization in pulmonary tuberculosis. Chest 127: 643–651. Powell, S., Forslund, K., Szklarczyk, D., Trachana, K., Roth, A., Huerta-Cepas, J., et al. (2014) eggNOG v4.0: nested orthology inference across 3686 organisms. Nucleic Acids Res 42: D231–D239. Qin, J., Li, R., Raes, J., Arumugam, M., Burgdorf, K.S., Manichanh, C., et al. (2010) A human gut microbial gene catalogue established by metagenomic sequencing. Nature 464: 59–65. Qin, J., Li, Y., Cai, Z., Li, S., Zhu, J., Zhang, F., et al. (2012) A metagenome-wide association study of gut microbiota in type 2 diabetes. Nature 490: 55–60. Redford, P.S., Murray, P.J., and O’garra, A. (2011) The role of IL-10 in immune regulation during M. tuberculosis infection. Mucosal Immunol 4: 261–270. Reichardt, N., Duncan, S.H., Young, P., Belenguer, A., McWilliam Leitch, C., Scott, K.P., et al. (2014) Phylogenetic distribution of three pathways for propionate production within the human gut microbiota. ISME J 8: 1323–1335. Sachdeva, P., Misra, R., Tyagi, A.K., and Singh, Y. (2010) The sigma factors of Mycobacterium tuberculosis: regulation of the regulators. FEBS J 277: 605–626. Sajid, A., Arora, G., Singhal, A., Kalia, V.C., and Singh, Y. (2015) Protein phosphatases of pathogenic bacteria: role in physiology and virulence. Annu Rev Microbiol 69: 527–547. Samuelson, D.R., Welsh, D.A., and Shellito, J.E. (2015) Regulation of lung immunity and host defense by the intestinal microbiota. Front Microbiol 6: 1085. Segain, J.P., Raingeard de la Bletiere, D., Bourreille, A., Leray, V., Gervois, N., Rosales, C., et al. (2000) Butyrate inhibits inflammatory responses through NFkappaB inhibition: implications for Crohn’s disease. Gut 47: 397–403. Shafiani, S., Tucker-Heard, G., Kariyone, A., Takatsu, K., and Urdahl, K.B. (2010) Pathogen-specific regulatory T cells delay the arrival of effector T cells in the lung during early tuberculosis. J Exp Med 207: 1409–1420. Shannon, C.E. (1997) The mathematical theory of communication. 1963. MD Comput 14: 306–317. Shetty, S.A., Marathe, N.P., and Shouche, Y.S. (2013) Opportunities and challenges for gut microbiome studies in the Indian population. Microbiome 1: 24. Singh, A., Dey, A.B., Mohan, A., Sharma, P.K., Mitra, D.K., and Unutmaz, D. (2012) Foxp31 regulatory T cells among tuberculosis patients: impact on prognosis and restoration


18 A. Maji et al. of antigen specific IFN-gamma producing T cells. PLoS One 7: e44728. Smith, P.M., Howitt, M.R., Panikov, N., Michaud, M., Gallini, C.A., Bohlooly-Y, M., Glickman, J.N., and Garrett, W.S. (2013) The microbial metabolites, short-chain fatty acids, regulate colonic Treg cell homeostasis. Science 341: 569–573. Stanley, S.A., and Cox, J.S. (2013) Host-pathogen interactions during Mycobacterium tuberculosis infections. Curr Top Microbiol Immunol 374: 211–241. Team, R.C. (2014) R: A Language and Environment for Statistical Computing. Vienna, Austria: R Foundation for Statistical Computing. Trehan, I., Goldbach, H.S., LaGrone, L.N., Meuli, G.J., Wang, R.J., Maleta, K.M., and Manary, M.J. (2013) Antibiotics as part of the management of severe acute malnutrition. N Engl J Med 368: 425–435. Truong, D.T., Franzosa, E.A., Tickle, T.L., Scholz, M., Weingart, G., Pasolli, E., et al. (2015) MetaPhlAn2 for enhanced metagenomic taxonomic profiling. Nat Methods 12: 902–903. Vital, M., Howe, A.C., and Tiedje, J.M. (2014) Revealing the bacterial butyrate synthesis pathways by analyzing (meta)genomic data. MBio 5: e00889. Winglee, K., Eloe-Fadrosh, E., Gupta, S., Guo, H., Fraser, C., Bishai, W., and Cardona, P.-J. (2014) Aerosol Mycobacterium tuberculosis infection causes rapid loss of diversity in gut microbiota. PLoS One 9: e97048. World Health Organization. (2014) Global tuberculosis report. Geneva: WHO. Yatsunenko, T., Rey, F.E., Manary, M.J., Trehan, I., Dominguez-Bello, M.G., Contreras, M., et al. (2012) Human gut microbiome viewed across age and geography. Nature 486: 222–227. Zheng, W., McLerran, D.F., Rolland, B., Zhang, X., Inoue, M., Matsuo, K., et al. (2011) Association between body-mass index and risk of death in more than 1 million Asians. N Engl J Med 364: 719–729. Zhu, W., Lomsadze, A., and Borodovsky, M. (2010) Ab initio gene identification in metagenomic sequences. Nucleic Acids Res 38: e132.

Supporting information Additional Supporting Information may be found in the online version of this article at the publisher’s web-site: Additional File S1. Word document file containing additional information; Supporting Information Table S1: Details of subjects chosen for the study; Supporting Information Fig. S1: Principal Component Analysis of unweighted UniFrac distances from 16S rRNA gene data set and Hellinger distances calculated from gene proportions; Supporting Information Fig. S2: Comparison of mean proportions of bacterial families in HC and TBZ obtained from MetaPhlAn2 analysis; Supporting Information Fig. S3: Comparison of

mean proportions of bacterial genera in HC and TBZ obtained from MetaPhlAn2 analysis; Supporting Information Fig. S4: Comparison of mean proportions of bacterial species in HC and TBZ obtained from MetaPhlAn2 analysis; Supporting Information Fig. S5: Correlation analysis of TBZ/ HC ratios obtained from whole-genome data set by readmapping and MetaPhlAn2 method; Supporting Information Fig. S6: The coverage of sequencing reads in Integrated gene catalogue, and Combined (Integrated 1 Indian) gene catalogue. Supporting Information Fig. S7: Spearman’s rank correlation of significant metagenomics species and enriched pathways; Supporting Information Fig. S8: KEGG pathway maps of (a) butyrate and (b) propionate metabolism. Additional File S2. Supporting Information Table S2, Summary of 16S rRNA gene sequencing reads Additional File S3. Supporting Information Table S3, Summary of whole genome shotgun sequencing reads Additional File S4. Supporting Information Table S4, PCA covariate correlation p values for (a) UniFrac distances and (b) Hellinger distances Additional File S5. Supporting Information Table S5, Unweighted UniFrac distances from the rarefied OTU table Additional File S6. Supporting Information Table S6, Taxa comparison between HC and TBZ based on 16S rRNA gene sequence data with p values calculated by Wilcoxon Rank sum test Additional File S7. Supporting Information Table S7, List of reference genomes Additional File S8. Supporting Information Table S8, Read-mapping comparison of reference genome-based and MetaPhlAn2 methods, Summary of bacterial abundance obtained by MetaPhlAn2 analysis Additional File S9. Supporting Information Table S9, List of metagenomic species obtained in Wilcoxon Rank sum test on comparison of HC and TBZ Additional File S10. Supporting Information Table S10, List of 195 MGCs obtained in MWAS analysis Additional File S11. Supporting Information Table S11, List of significantly enriched pathways in HC and TBZ as determined by Fisher’s exact test Additional File S12. Supporting Information Table S12, List of pathways associated with HC and TBZ as determined by Reporter features algorithm. Additional File S13. Supporting Information Table S13, Table showing Spearman’s rank correlation coefficients between metagenomic species and enriched pathways. Additional File S14. Supporting Information Table S14, List of metagenomic species obtained in Wilcoxon Rank sum test on comparison of TBZ and TBM/TBW. Additional File S15. Supporting Information Table S15, List of pathways in TBZ and TBM as determined by Fisher’s exact test.
