We included TB patients from an ongoing prospective cohort that studies the clinical and molecular epidemiology of TB in the Temeke district of Dar es Salaam, Tanzania (TB-DAR). Tanzania is among the top 20 countries with the highest TB incidence; 142,000 new TB cases were estimated in Tanzania during 2018, 28% of them were co-infected with HIV (WHO, 2019). The city of Dar es Salaam has the highest TB incidence in the country (20 % of TB cases notified during 2018) (NTLP, 2018). The Temeke district is a densely populated urban setting accounting for one third of the TB cases from Dar es Salaam (Said et al., 2017).

We conducted a cross-sectional study nested within the ongoing TB-DAR cohort. The study population and the procedures of this cohort have been previously described in detail (Mhimbira et al., 2016, 2017, 2019; Steiner et al., 2016; Hiza et al., 2017; Said et al., 2017; Hella et al., 2018; Sikalengo et al., 2018). Briefly, at the Temeke district hospital and since November 2013, TB-DAR have been recruiting sputum smear positive or Xpert MTB/RIF positive adult TB patients (≥18 years of age). TB was further confirmed by sputum culture in Lwenstein-Jensen (LJ) solid media; clinical isolates were further analyzed by lineage-specific allele probes in real-time PCR for singleplex SNP-typing according to standard protocols (Applied Biosystems, Carlsbad, USA) and as previously described (Stucki et al., 2018). In this study, we only included newly diagnosed TB patients, without previous history or diagnosis of TB, recruited between November 2013 and November 2015. Included TB cases were randomly chosen to obtain a representative subset of the entire cohort during the recruitment period.

At the time of TB diagnosis, participants were interviewed to collect data on socio-demographic characteristics, lifestyle, symptoms, previous use of medications, and health-seeking behavior (Said et al., 2017). Patients underwent physical and chest X-ray examination (CXR). Before initiation of TB treatment, biological specimens (sputum, naso-pharyngeal swabs, blood, urine, and stool) were collected for further investigation of anemia and co-infections, as previously described (Mhimbira et al., 2017, 2019; Hella et al., 2018). Investigated co-infections included HIV (Mhimbira et al., 2017; Hella et al., 2018), helminths (Mhimbira et al., 2017), and respiratory pathogens (Mhimbira et al., 2019). For HIV-positive patients, CD4+ T cell counts were also obtained (Mhimbira et al., 2017).

To ensure comparability and quality of sputum samples, health-care workers provided video-guided instructions to all participants and asked them to only submit specimens collected during early-morning. The use of the sputum submission instructional video, for improvement of sputum quality, was previously validated by Mhalu and colleagues (Mhalu et al., 2015). At the time of specimen reception, laboratory technicians visually assessed quality and volume; second samples were requested if salivary-like (transparent and watery specimen with bubbles) specimens were submitted. Accepted sputum samples were transported from Temeke district hospital at 4°C to the Bagamoyo Research and Training Center for processing. Total DNA was extracted from sputum specimens using the QIAamp DNA Kit (Qiagen, Germany) according to the supplier's instructions (Spin Protocol for DNA purification from Blood or Body Fluids).

We followed Illumina's 16S amplicon sequencing library preparation protocol (Illumina, 2013). The protocol included: (i) PCR amplification of the V3-V4 region of the 16S rRNA gene with primers PCR1_Forward (50 bp, 5′-TCGTCGGCAGCGTCAGATGTGTATAAGAGACAGCCTACGGGNGGCWGCAG-3′), and PCR1_Reverse (55 bp, 5′-GTCTCGTGGGCTCGGAGATGTGTATAAGAGACAGGACTACHVGGGTATCTAATCC-3′); (ii) library preparation with the Nextera XT Index Kit, and (iii) paired-end sequencing (2 × 300 bp, v3 chemistry) on the MiSeq platform (Illumina, San Diego, CA), samples were distributed across three sequencing runs. To account for potential sequencing batch effects, we included the sequencing run number as a co-factor in multivariate models.

DNA samples were first selected based on their quality (limited degradation and high concentration) as assessed by DNA gel electrophoresis (1:10 DNA dilutions run on 1% agarose gels at 90 volts for 45 min). DNA concentration was determined by PicoGreen quantification with Qubit and normalized to 0.2 ng/ul. Dual-indexed paired-end libraries were prepared using the Nextera XT DNA Library Prep kit; the standard protocol was followed at the sequencing facility of FISABIO (Foundation for the Promotion of Health and Biomedical Research in the Valencian Community). Paired-end sequencing (2 × 150 bp) on the HiSeq 2,500 platform was done at the Department of Biosystems Science and Engineering (D-BSSE) of ETH Zeich. To minimize batch effects, we used three sequencing runs to repeatedly sequence all samples and sample pools per lane combined samples picked in a random order. Background controls (water) and a mock community (20 strain even mix genomic material, ATCC® MSA-1002™) were included during the entire workflow.

Newly diagnosed TB cases were defined as TB patients that never received TB treatment or received TB treatment for <1 month. We further stratified those patients by multiple aspects of disease manifestations: mycobacterial load in their sputum, severity of clinical findings (signs and symptoms), and presence of radiologic signs (lung parenchymal infiltrates, cavities, lymphadenopathy, micronodules, pleural efussion, etc.) in chest X-rays (CXRs).

To define sputum mycobacterial load as low or high, we used sputum acid-fast bacilli (AFB) smear results. Following the World Health Organization/International Union Against Tuberculosis and Lung Disease (WHO/IUATLD) guidelines, sputum AFB smear results were expressed as quantitative categories (“Scanty,” “1+,” “2+,” and “3+”) that grade the number of AFB per number of microscopic fields (Mhimbira et al., 2017). As previously described (Mhimbira et al., 2019), we defined sputum mycobacterial burden as high if AFB smear results were “2+” (1–10 AFB per field, 50 fields) or “3+” (> 10 AFB per field, 20 fields); and as low if AFB smear results were “Scanty” (1–9 AFB in 100 fields) or “1+” (10–99 AFB in 100 fields).

To define clinical findings for pulmonary tuberculosis as mild or severe, we adapted the TB score (Mhimbira et al., 2017) defined by Wejse and colleagues (Wejse et al., 2008). The TB score was expressed as cumulative points (from 0 to 12) that quantify the number of clinical findings observed in a patient during physical examination. Twelve clinical findings were considered, each one counted as one if present: (i) cough, (ii) haemoptysis, (iii) dyspnoea, (iv) chest pain, (v) night sweating, (vi) anemic conjuctivae, (vii) positive finding at lung auscultation (crepitation, rhonci, subdued, or complete absence of respiratory sounds), (viii) fever (temperature > 37°C), (ix) mid upper arm circumference (MUAC) <220 mm, (x) MUAC <200 mm, (xi) body mass index (BMI) <18 kg m-2, and (xii) BMI <16 kg m-2 (see Figure 1). As in Mhimbira and colleagues (Mhimbira et al., 2017), we defined clinical findings for pulmonary tuberculosis as severe if the TB score was ≥6; and as mild, otherwise. To screen for lung abnormalities, double readings of CXRs were performed by board-certified radiologists, and discrepancies were resolved by an independent reader.

To account for differences in disease manifestations due to delay in diagnosis, duration of diagnostic delay was included based on the longest reported TB-related symptom and categorized into: ≤ 3 weeks and >3 weeks, as previously described by Said and colleagues (Said et al., 2017). Differences in clinical manifestations might be associated to differences in the genetic background of MTB, thus we included the SNP-typing results of the MTB clinical isolates.

Socio-demographic variables included sex and age. Physical health status parameters included BMI, hemoglobin (Hb) levels (kg dL-1), smoking, and alcohol abuse. To assess adult nutritional status, we followed WHO's BMI-based definitions (Bailey and Ferro-Luzzi, 1995), as follows: (i) “Underweight” (BMI <18.5), (ii) “Normal weight” (BMI, 18.5–24.9), and iii) “Obesity” (BMI ≥ 25). Similarly, we followed WHO's Hb cut-offs to define “Non-Anemia” (Hb ≥ 12 for women, and Hb ≥ 13 for men), “Mild” (11–11.9 for women, and 11–12.9 for men), “Moderate” (8.0–10.9 for women and men), and “Severe” (<8.0 for women and men) anemia (WHO, 2011).

16S-rRNA-gene amplicon sequencing (hereafter 16S-AS) reads were processed with QIIME2 plugins v2019.7 (Bolyen et al., 2019). We first carried out a quality inspection of paired-end reads with QIIME2's “demux” plugin which revealed reverse reads had poor quality at the 3'end and had to be trimmed off before merging read pairs with DADA2 v1.10 (Callahan et al., 2016); this resulted in a 90% loss of the data (analysis is available at https://git.scicore.unibas.ch/TBRU/tbdarbiome_cases/-/blob/master/notebooks/00_MOT_qiime2-qc-denoising-features_pairedreads.ipynb). Therefore, we decided to only use the forward reads of every pair, which were processed as follows. First, we run the denoise-single method of QIIME2's dada2 plugin, followed default parameter settings but adjusted the number of bases trimmed at the 5′- and 3′- ends. To create a phylogenetic tree of the resulting amplicon sequence variants (ASVs), we run the q2-fragment-insertion plugin which uses the SAT-enabled phylogenetic placement (SEPP, Mirarab et al., 2012) algorithm to insert our ASVs into the reference phylogeny of 16S rRNA gene sequences (Greengenes v13.8, reference sequences clustered at 99% sequence similarity) (Janssen et al., 2018); unplaced ASVs were removed. Finally, we used the q2-feature-classifier plugin (Bokulich et al., 2018) to assign taxonomic classifications to our ASVs; we first retrained the classifier with 16S V3-V4 fragments extracted from the Greengenes v13.8 reference sequences. The workflow is available as a Jupyter notebook at https://git.scicore.unibas.ch/TBRU/tbdarbiome_cases/-/blob/master/notebooks/01_MOT_qiime2-qc-dada2_denoising_features_read1.ipynb.

Whole-metagenome shotgun (WMS) reads were processed with an in-house analysis pipeline that performs quality pre-processing, taxonomic profiling, and functional profiling. Quality pre-processing included adapter/quality trimming, and filtering with fastp (Chen et al., 2018); removal of reads derived from human DNA with BBTools' (Bushnell, 2014) bbsplit; and removal of duplicated reads with BBTools' clumpify. Taxonomic profiling was performed with MetaPhlAn2 (Truong et al., 2015). The pipeline is available at https://git.scicore.unibas.ch/TBRU/MetagenomicSnake and details of its execution for this publication are available in a Jupyter notebook at https://git.scicore.unibas.ch/TBRU/tbdarbiome_cases/-/blob/master/notebooks/02_MOT_metasnk-wms-data-processing.ipynb. For each species detected, we obtained the reported mode of metabolism (i.e., oxygen requirement) by manually querying https://bacdive.dsmz.de, https://www.lgcstandards-atcc.org/, https://microbewiki.kenyon.edu/, http://www.homd.org, and Pubmed Central® (PMC); the resulting table is available at https://git.scicore.unibas.ch/TBRU/tbdarbiome_cases/-/blob/master/data/raw/metadata/species_oxygen_tolerance.csv.

We used R package decontam v1.6 (Davis et al., 2018) to identify contaminant DNA sequences. For the 16S-AS dataset, decontam used the DNA concentration of sputum samples to identify contaminant ASVs whose frequency varied inversely with total DNA concentration (frequency-based identification). For the WMS dataset, frequency-based identification was complemented with presence of contaminants in background controls (prevalence-based identification) and in a mock community (20 strain even mix genomic material, ATCC® MSA-1002™). The analysis is available at https://git.scicore.unibas.ch/TBRU/tbdarbiome_cases/-/blob/master/notebooks/03_MOT_decontam.ipynb.

To quantify diversity within each sputum sample (alpha diversity), we used the 16S-AS taxonomic profiles to compute Faith's phylogenetic diversity (PD) metric (Faith, 1992), as implemented in QIIME2 v2019.7 (Bolyen et al., 2019). We favored the PD metric as it considers phylogenetic differences between taxa which makes PD to account the limited contribution of closely related taxa to diversity (Matsen, 2015). To minimize the effect of differences in sequencing depth, we computed rarefaction curves to set a threshold for even sampling (i.e., equal number of sequences) across all samples before computing the PD metric. Analysis available at https://git.scicore.unibas.ch/TBRU/tbdarbiome_cases/-/blob/master/notebooks/06_MOT_alpha-rarefaction.ipynb.

To quantify differences in community composition among samples (beta diversity) and identify variability patterns, we computed Aitchison distances which acknowledge the compositional nature of the taxonomic profiles (Gloor et al., 2017; Quinn et al., 2018b). Aitchison distances correspond to Euclidean distances among taxonomic profiles (i.e., taxon abundances) that were centered log-ratio (CLR) transformed. The CLR transformation takes compositional vectors (i.e., relative abundances) from a constrained space (i.e., the unit simplex) to an unconstrained space of logratio vectors (logarithm of the ratio between the abundance of a taxon in a sample and the geometric mean of all taxon abundances in the same sample). Resulting CLR-transformed vectors are scale-invariant and subcompositionally coherent which implies that rarefying the taxonomic profiles to a constant value across samples will not have significant effects on the relationships between samples and between taxa; therefore we computed Aitchison distances on non-rarefied taxonomic profiles. To circumvent the limitation of undefined values when computing the logarithm of zero values, zero abundances were replaced following Martin-Fernndez et al.'s Bayesian-multiplicative replacement strategy (Martín-Fernández et al., 2014); implemented in cmultRepl of the R package zCompositions v1.3.4. To visualize community composition profiles across samples, we created heatmaps with sputum samples and taxa sorted according to agglomerative hierarchical clustering based on Aitchison distances; analysis available at https://git.scicore.unibas.ch/TBRU/tbdarbiome_cases/-/blob/master/notebooks/08_MOT_taxonomic-summaries.ipynb. To visualize relationships among taxa across sputum samples, we created relative variation biplots following Aitchison and colleagues rules for creation and interpretation of such biplots (Aitchison and Greenacre, 2002). Analysis is available at https://git.scicore.unibas.ch/TBRU/tbdarbiome_cases/-/blob/master/notebooks/10_MOT_inference-of-genus-genus-interactions.ipynb.

To infer associations among taxa, we created genus-level and species-level interaction networks with SParse InveresE Covariance Estimation for Ecological ASociation Inference (SPIEC-EASI); available as an R package (SpiecEasi v1.0.7). SPIEC-EASI uses CLR-transformed abundances to infer a graph model where nodes represent taxa and edges represent associations between taxa that cannot be explained by alternative paths in the graph (Kurtz et al., 2015). We selected the neighborhood selection method of SPIEC-EASI. Inferred interactions were also confirmed by Spearman's rank-order correlation. Additionally, we also used the SparCC method (Friedman and Alm, 2012), as implemented in the R package SpiecEasi.

To infer species-level interaction networks, we used the WMS-S taxonomic profiles and considered only those species present in at least 3 samples. To infer genus-level interaction networks, we used the 16S-AS taxonomic profiles and considered only those genera present in at least 5% of the samples included in the dataset; genera are more likely to be shared across samples, thus we increased the presence threshold. Analyses are available at https://git.scicore.unibas.ch/TBRU/tbdarbiome_cases/-/blob/master/notebooks/10_MOT_inference-of-genus-genus-interactions.ipynb and at https://git.scicore.unibas.ch/TBRU/tbdarbiome_cases/-/blob/master/notebooks/11_MOT_inference-of-spp-spp-interactions.ipynb.

We performed statistical analysis with the R software environment v3.6. Associations of categorical demographic and clinical characteristics with TB disease manifestations were assessed by chi-squared test or Fisher's exact test if expected frequencies were below five. Associations with continuous variables were assessed by student t-tests or Wilcoxon rank-sum test when normality could not be assumed; normality was evaluated with Shapiro-Wilk test and Q-Q plots. Single test significance level was Bonferroni adjusted for multiple comparisons. Analyses are available at https://git.scicore.unibas.ch/TBRU/tbdarbiome_cases/-/blob/master/notebooks/07_MOT_characteristics-of-cohort.ipynb.

To test associations of Faith's PD with TB disease manifestations (mycobacterial load, severity of clinical findings, and peresence of radiologic signs), we used the R-package car v3.0 to create two multi-way ANCOVA models with log transformed Faith's PD as response variable. With the first ANCOVA model, we simultaneously assessed the marginal effects of abnormal CXR-findings and high mycobacterial load of sputum; a decision based on previous analysis indicating that these clinical manifestations were not associated with each other. The model was adjusted for age, sex, physical health status parameters (underweight, anemia, smoking, and alcohol abuse), co-infections (HIV, helminths, viral, and bacterial pathogens), season, delay in diagnosis, non-TB medications, sequencing depth and sequencing batch; two-way interactions of CXR-findings or mycobacterial sputum burden with age, sex, physical health status parameters, co-infections and season were also evaluated. With the second model, we tested the marginal effect of clinical findings severity (as assessed by a TB score) while adjusting for delay in diagnosis, and the covariates included in the first model; except for BMI and anemia which are parameters used to compute the TB score. Interaction terms in both models were selected using stepwise-selection based on the Akaike information criterion (AIC), as implemented in the function step of the R package stats. On selected models, we tested normality of residuals, homoscedascity, non-multicolinearity, independence of errors and effect of influential observations. To assess the effect sizes of factors within selected models, we computed the partial eta squared statistic (hp2), which is the proportion of the Sum of Squares (SS) of the effect and the error that is attributable to the effect (Maher et al., 2013). Post-hoc tests included t-tests on adjusted means with R-package emmeans v1.4. The analysis is available at https://git.scicore.unibas.ch/TBRU/tbdarbiome_cases/-/blob/master/notebooks/12_MOT_alpha-diversity-and-clinical-features.ipynb.

Similarly, to investigate if compositional and structural changes of the microbial communities in the sputum of TB patients were associated with disease manifestations of pulmonary TB, we performed transformation-based redundant analysis (tb-RDA) with stepwise selection of interaction terms, and Permutational Analysis of Variance (PERMANOVA) as implemented in the R-package vegan v2.5. Response variables were CLR-transformed genus-level abundances. Continuous explanatory variables were always scaled. When comparing means among groups, p-values were adjusted for multiple comparisons following Holm-Bonferroni method. Analyses are available at https://git.scicore.unibas.ch/TBRU/tbdarbiome_cases/-/blob/master/notebooks/13_MOT_beta-diversity-and-clinical-features.ipynb.

From November 2013 to November 2015, 663 new TB cases were enrolled at the Temeke district hospital in Dar es Salaam, Tanzania. Although it was originally planned to include all the sputa stored during these 2 years of patient enrollment, we did not have enough DNA extraction kits and ended up choosing a subset. We selected at random 4 9 × 9 boxes of sputa stored at –20°C; this resulted in 324 sputa. We still had reagents for another 10 DNA extractions, therefore we selected 10 more sputa from a randomly-picked 5th stored box. Thus, we included 334 (53%) sputa to investigate the potential relationship between sputum microbial composition and three aspects of TB-disease manifestations, at time of diagnosis, which include: (i) high mycobacterial load in sputum, (ii) severe clinical findings (signs and symptoms), and (iii) chest X-ray findings (see Figure 1). Sputum volumes were consistent across categories of TB-disease manifestations (Supplementary Figure S1).

In our study population and across categories of TB-disease manifestations, the distributions of demographics, physical health status, co-infections, non-TB medication, delay in diagnosis, season, and MTB genetic background are summarized in Table 1. Consistent with previous studies on this cohort (Mhimbira et al., 2017, 2019; Hella et al., 2018), most TB cases were males (72%), had anemia (74%), and were underweight (53%). A large proportion of TB cases were recruited during the Dry season (42%). Also, significant proportions were co-infected with HIV (24%), helminths (34%), viral (21%) and bacterial (34%) respiratory pathogens; 16% were smokers and 19% were alcohol abusers. Additionally, 69% reported a delay in TB diagnosis of more than 3 weeks, and 94% reported to have received non-TB medications before diagnosis.

Pulmonary TB manifests in CXRs as parenchymal infiltrates/consolidations, cavities, pleural effusion, lymphadenopathy, and micronodules. As summarized in Table 2, in this cohort, parenchymal infiltrates were the most common (63%, n = 157), followed by cavities (36%, n = 89), pleural effusion (12%, n = 29), lymphadenopathy (10%, n = 25), and micronodules (5%, n = 13).

We analyzed the total DNA of all sputum samples (334 patients) by 16S rRNA gene amplicon sequencing (16S-AS). We obtained about 24.6 million reads, median reads per sample was 55,672 (IQR: 36,666–98,143). After filtering, denoising and chimera removal with DADA2, about 12.4 million sequences remained; 28,231 (IQR: 19,227–45,916) median reads per sample. Samples had a median of 171 (IQR:111–252) amplicon sequence variants (ASVs). Additional filtering removed 8 ASVs flagged as potential contaminants by Decontam, as described in section Identification of Potential Contaminants (Supplementary Figure S2), and 27 observations which had <12,000 reads (threshold at which most samples approached plateau in richness rarefaction curves; Supplementary Figure S3), which resulted in 307 16S-AS abundance profiles.

Based on DNA concentration and quality (limited fragmentation, see methods), 125 samples were selected for Whole-Metagenome Shotgun sequencing (WMS-S). We obtained about 3.96 × 10⁹ reads, median reads per sample was 28,948,058 (IQR: 22,116,534–37,234,300), and the percentage of reads annotated as human ranged from 2.6 to 92.1. After filtering out human and low-quality reads, about 1.11 × 10⁹ reads remained; median reads per sample was 4,836,788 (IQR: 3,366,538–7,818,888). Additional filtering removed five species identified as contaminants (Supplementary Figure S4) by Decontam, three species not expected as part of a mock community (Supplementary Figure S5), 19 samples with more than 20% reads assigned to contaminant species, two samples with 100% of unclassified reads, and 15 samples with <50% of reads assigned to human DNA; thus, 89 WMS-S abundance profiles were included in further analyses.

16S-AS and WMS-S provided complementary snapshots of the microbial composition in the sputum of TB patients (see Figure 2). Although 16S-AS was able to profile a larger number of samples, species-level resolution was limited; most (86%) of the ASVs detected by DADA2 did not have a species classification. WMS-S on the other hand, provided a cross-domain and species-level resolution, evidenced by the detection of viruses and fungi, and by the accurate taxonomic profiling of a mock community composed of even abundances of 20 species (ATCC® MSA-1002^TM; Supplementary Figure S5A). However, WMS-S resolution came at the cost of losing detection of low abundant taxa; WMS-S detected 78 bacterial genera while 16S-AS detected 319 genera, and genera not-detected by WMS-S had mean relative abundances generally below 0.01% (Supplementary Figure S6).

Based on these observations, we decided to use the 16S-AS dataset for genus-level analysis as this dataset included more samples and captured 2.5 times more genera than the WMS-S dataset. Considering the limited species-level accuracy in our 16S-AS and the accurate taxonomic profiling of a mock community by our WMS-S processing pipeline, we used the WMS-S dataset for species-level analysis.

Taxonomic assignments derived from 89 WMS abundance profiles showed that microbial communities in the sputum of TB patients were primarily composed of bacteria (90.6, 83.7–97%; mean, 95% CI) and viruses (9.4, 5.6–19.6%; mean, 95% CI) (Figure 2A). The presence of members of the Eukaryota or Archaea domains was rather rare, or below detection limit: only one sample, among those screened by WMS-S, had reads assigned to Candida albicans and four samples, among those screened by 16S-AS, had Metanobrevibacter ASVs.

Taxonomic assignments derived from 307 16S-AS abundance profiles showed that the bacterial content of the sputum in TB patients was dominated by 10 phyla (Figure 2B): Firmicutes, Proteobacteria, Bacteroidetes, Fusobacteria, Actinobacteria, TM7 (Saccharibacteria), Spirochaetes, GN02 (Gracilibacteria), SR1 (Absconditabacteria), Tenericutes; together they accounted for 99.3% of the total bacterial relative abundance measured by 16S-AS. At the genus level (Figure 2C), 16 genera which had mean relative abundances of at least 1% covered approximately 82% of the total bacterial relative abundance. In agreement with previous studies (Krishna et al., 2016; Hong et al., 2018; Sala et al., 2020), the top 10 most abundant genera included: Streptococcus, Neisseria, Veillonella, Prevotella, Haemophilus, Porphyromonas, Fusobacterium, Campylobacter, Rothia, and Lautropia.

DNA Viruses were detected in 37% of the samples screened by WMS-S; the most common viruses were Human Herpesvirus 4 (HHV-4, also known as Epstein-Barr virus) and Torque Teno Virus (TTN) (Figure 2D). A few sputum samples (n = 7) had excessive viral relative abundance (> 51 %, equivalent to mean+2SD), dominated by one or few species (Streptococcus phages, HHV-4, and unclassified Roseolovirus). Although, both HHV-4 and TTN viruses are generally carried asymptomatically, these viruses can expand or cause (respiratory) infections, specially in immuno-compromised individuals (Thom and Petrik, 2007; Reid et al., 2016). Thus, we investigated if the detection of viruses in the sputum of TB patients, which in this cohort is dominated by both HHV-4 and TTN viruses, was linked to HIV. The proportion of HIV in TB patients with detectable viral DNA in their sputa was 43%, as compared to 21% of TB patients with no viral DNA detected in their sputa. However, a chi-squared test of independence showed a borderline significant association between the two variables (p = 0.06).

To look for associations between the identified viruses and bacterial species, we included virus abundances when reconstructing the network of microbial associations with SpiecEasi (Supplementary Figure S21A). SpiecEasi inferred direct positive associations between TTVs, Moraxella catarrhalis and Streptococcus vestibularis (Supplementary Figure S21A). However, we have to be skeptical about these associations because they were derived from the detection of both M. catarrhalis and S. vestibularis in few sputum samples (n = 3), see Supplementary Figure S21B.

Although our observations of the viral fraction of the microbial communities in the sputum of TB patients are exploratory, they highlight the overlooked importance of the virome in the airways of TB patients co-infected with HIV.

Our findings should be considered in the light of some limitations. First, a large proportion (94%) of TB patients in the cohort reported the previous intake of non-TB medications; patients mostly took penicillin-derivatives (78%). Similar to the previously reported lack of effect of TB treatment on sputum microbial composition (Sala et al., 2020), the distribution of non-TB medications neither differed among the levels of TB-disease manifestations nor were they associated with compositional variation of the sputum microbial communities. Therefore, we are confident that the intake of non-TB medications did not bias the observed associations. Nonetheless, given the well-known impact of medications on human microbiomes, and the serious burden self-medication or miss-treatment poses in low and middle-income countries, we think previous medication-intake should be considered in study designs; excluding participants that received previous medications might neglect the role of medication-intake as a potential mediator or modifying factor in the microbiome-disease interplay in particular study settings.

Second, we followed sputum collection practices to ensure that sputum specimens were produced and collected consistently across participants. These practices included the use of instructional videos, collection of early morning sputum, and qualitative assessment of the specimen's color and viscosity at the time of collection by experienced laboratory technicians. Although these practices were proven to reduce the collection of salivary-like specimens and improve the quality of the sputum for detection of MTB (Mhalu et al., 2015), there are additional quantitative measurements that can improve the discrimination of sputum from salivary-like specimens. For instance, the number of squamous epithelium cells or leukocytes, or their ratio, quantified by smear microscopy (Wong et al., 1982). We encourage future studies to incorporate such quantitative measurements to improve the assessment of sputum quality for microbiome studies.

Third, our study observed no correlation between Mycobacterial reads abundance and AFB smear results; likely as a result of our DNA extraction methodology which did not specifically targeted for MTB DNA or for removal of host's DNA. We made this decision because both MTB DNA enrichment and human DNA depletion would result in distorted compositional profiles of the microbial communities that we are interested in characterizing. For instance, MTB enrichment would inevitably affect the estimation of the MTB abundance relative to other microorganisms. On the other hand, due to enzymatic/chemical/mechanical steps to differentially filter/lyse host cells and degrade exposed DNA, human DNA depletion steps can cause an overall loss of microbial DNA or bias toward microorganisms that are less susceptible to such steps. Nonetheless, we must acknowledge that mycobacterial relative abundances estimated from both 16S-AS and WMS-S reads might be an under-representation or miss-representation of the actual MTB load in the sputum.

Fourth, even though sputum originates in the lower airways, its passage through the upper airways during expectoration inevitably results in a sample mixed with microbes from the upper respiratory tract and the oral cavity. Therefore, future studies with respiratory samples better rendering the lung lesion environment must verify the reported associations with abnormal CXRs (lung parenchymal infiltrates in particular). Fifth, we report microbial interactions based on the relationships of taxon abundances across a large set of sputum microbial assemblages. However, further experimentation is needed to imply ecological processes such as antagonism or co-existence; for instance, correlations based on absolute abundances measured by quantitative PCRs would confirm the microbial associations reported in this study.

And sixth, our analysis of the WMS-S dataset was limited by the high host-to-microbial DNA ratio. On average 90% of the reads were mapped to a human reference genome and filtered out; remaining reads were not enough for a good coverage of the microbial community and therefore biased against low-abundant members. This bias affects the sensitivity of functional profiling methods to estimate the contribution of microbial species with low-covered genomes to a functional pathway or to a gene family. For this reason we decided to limit our analysis of the WMS-S dataset to perform only a species-level taxonomic profiling that could complement the 16S-AS profiles.

In summary, this study sheds light on the relationship between pulmonary TB and the host's airway microbiome. We have identified specific microbial interactions responsible for structuring the microbial communities in the sputum of TB patients. More importantly, our results suggest underweight nutritional status (BMI ≤ 18.5) as a determinant factor for sputum microbial associations with hypoxic lesions in pulmonary TB. Thus, our hypothesis that variations in TB disease manifestations are associated with microbial composition variation in the airways of TB patients, partly holds; the association is restricted to pulmonary damage and determined by underweight status. Finally, but importantly, our study points out a non explored nexus between TB and HIV, the airways virome. These new insights should guide further research to unravel the underlying mechanisms behind the observations presented here.