This study was approved by the Ethics Review Committee of Affiliated Dongyang Hospital of Wenzhou Medical University (2022-YX-268), and written informed consent was obtained from all participants.

A total of 69 pulmonary tuberculosis patients and 63 healthy controls were initially enrolled from the Affiliated Dongyang Hospital of Wenzhou Medical University, Zhejiang, China, from May 2021 to June 2022. All patients were recruited on the basis of TB symptomatic features (WS288-2008). Healthy controls (HCs) were recruited in hospital health examination center, who did not have close contact with TB patients at least one year. All participants had not taken any broad-spectrum antibiotics in the previous six months and without other basic diseases or major illness.

Blood samples were collected for biochemical analysis using the automated equipment and standard methods. The level of Albumin (Alb), Serum fasting blood glucose (Glu), Total Cholesterol (TC), Triglycerides (TG), High-Density Lipoprotein cholesterol (HDL), Low-Density Lipoprotein cholesterol (LDL), and Creatinine (Cr) were examined using a Hitachi 7600 automatic biochemical analyzer (HITACHI, Japan). The Hemoglobin (Hb) and White Blood Cell (WBC) were examined by Sysmex pocH-100i automatic blood Analyzer (SYSMEX, Japan).

Fresh stool samples (morning samples, 200 mg/aliquot) were collected using GUHE Flora Storage buffer (GUHE Laboratories, Hangzhou, China), and stored at room temperature until DNA extraction. DNA was extracted with nucleic acid extraction or purification reagent (GUHE Laboratories, GHFDE100, Hangzhou, China), and the quantity and quality of extracted DNAs were measured using the NanoDrop ND-1000 spectrophotometer (Thermo Fisher Scientific, Waltham, MA, USA) and agarose gel electrophoresis, respectively.

The V4 hypervariable regions of the bacterial 16S rRNA gene were amplified from each sample using the 515F(5’-GTGCCAGCMGCCGCGGTAA-3’) and 806R(5’-GGACTACHVGGGTWTCTAAT-3’) primers as previously described prior to sequencing. Phusion High-Fidelity PCR Master Mix with HF Buffer (PHUSION high fidelity DNA polymerase) was used for PCR. The following primary PCR cycling conditions were used: (1) Pre-denaturatio at 98℃ for 30 s; (2) 25 cycles of PCR at 98℃ for 15 s, 58℃ for 15 s, and 72℃ for 15 s; (3) a final extension at 72℃ for 1 min. PCR products were purified with AMPure XP Beads (Beckman Coulter, Indianapolis, IN), and the library was quantified with Qubit. After quantification, the Illumina NovaSeq6000 platform at GUHE Info Technology (GUHE, Hangzhou, China) was used for sequencing. The sequence data of this study have been uploaded in the GenBank Sequence Read Archive (SRA) of NCBI under the accession code BioProject PRJNA901399.

The original data of each sample were split according to the Barcode sequence and the primer sequence. The low-quality sequences were filtered through following criteria (Gill et al., 2006; Chen and Jiang, 2014): sequences that had a length of <150 bp, sequences that had average Phred scores of <20, sequences that contained ambiguous bases, and sequences that contained mononucleotide repeats of >8 bp. Paired-end reads were assembled using Vsearch V2.4.4 (–fastq_mergepairs –fastq_minovlen 5), than Operational taxonomic unit (OTU) picking using Vsearch v2.15.0, included dereplication(–derep_fulllength), cluster(-cluster_unoise), detectection of chimeras(-uchime3_denovo) (Rognes et al., 2016). The final effective tags were obtained by removing the chimeric sequences.

The filtered clean reads were clustered to get ASV statistical results. The representative sequences of OTUs were selected using default parameters. The representative sequences were annotated through QIIME2 (v2020.6) based on the SILVA138 database to further generate the OTU list. At each classification level, phylum, family, and genus were used to count the community composition of each sample.

Sequence reads for 16S rRNA gene amplicons were analyzed using QIIME2 software and R package (v3.2.0). QIIME2 was used to calculate alpha diversity index in OTU level, including Chao1, ACE, Shannon, Simpson index. Beta diversity analysis was calculated through Qimme software based on UniFrac distance metric (Lozupone and Knight, 2005; Lozupone et al., 2007). “Vegan” from R package was used to evaluate the markers of microbial community structure differentiation between groups through PERMANOVA (Permanent multivariate analysis of variance) (McArdle and Anderson, 2001). The classification and abundance of bacteria were visualized by MEGAN software and GraPhlAn. And the R package “VennDiagram” was used to generate a Venn diagram (Zaura et al., 2009). Taxa abundances at the phylum, class, order, family, genus and species levels were statistically compared among samples or groups by Kruskal test from R stats package (metagenomeSeq packages). LEfSe (Linear discriminant analysis effect size) was performed to detect differentially abundant taxa across groups using the default parameters (Segata et al., 2011). Microbial functions were predicted by PICRUSt2 (Phylogenetic investigation of communities by reconstruction of unobserved states, https://github.com/picrust/picrust2/), based on high-quality sequences. For each sample, functions were profiling by using the Omixer-RPM version 1.0 (https://github.com/raeslab/omixer-rpm) with KO redundants predicted from PICRUSt2. Spearman’s rank correlation between the intestinal flora and biochemical indexes was calculated with our own python scripts (Ramsey, 1989). Random forest analysis was applied to discriminating the samples from different groups using the R package “random Forest” with 1,000 trees and all default settings. Values were expressed as the mean ± SEM. All other comparisons were calculated by two-tailed Student’s t-test or one-way ANOVA followed by Tukey’s test using GraphPad Prism 6 software.

Table 1 showed the detailed clinical characteristics of the study participants. 69 PTB patients and 63 age- and sex-matched healthy controls (HCs) were enrolled in this study. Among the clinical biochemical indexes tested in the study, a significant increase in serum fasting blood glucose (Glu) was observed in PTB patients (P<0.05). The levels of Albumin (Alb), Creatinine (Cr), Hemoglobin (Hb), Triglycerides (TG), Total cholesterol (TC), High-density lipoprotein cholesterol (HDL) and Low-density lipoprotein cholesterol (LDL), were all found to be significantly lower in PTB patients than in HCs (P<0.001). There was no statistical difference in White blood cells (WBC) between the two groups.

In total, we generated 16,259,048 high-quality reads, and the reads were clustered into 7296 operational taxonomic units (OTUs) based on a similarity of 97%. A Venn diagram showed that 1141 OTUs and 732 OTUs were unique in the healthy controls and PTB patients, respectively ( Figure 1A ). The rarefaction curves showed that the curve entered in plateau, indicating that the sequencing data was large enough to reflect the bacterial diversity of the intestinal flora in HCs was higher than that in the PTB patients ( Figure 1B ). The alpha diversity indexes, including Shannon, Simpson, Chao1 and ACE reflected that the diversity of fecal microbiota in PTB patients was significantly lower than in HCs (P<0.01, P<0.05, P<0.001, P<0.001, respectively) ( Figure 1C ). The principal coordinate analysis (PCoA), based on Bray-Curtis, unweighted UniFrac and weighted UniFrac distances were performed to identify the difference in fecal microbiota composition between PTB patients and HCs. As shown in Figure 1D , there was a tendency of separation and had a partial overlap between two groups, which demonstrated that the total community of gut microbiota in the two groups was different.

The compositions of the fecal microbiota in the PTB patients and the HCs were assessed at different taxonomic levels, and the top 20 relative abundance of bacteria were shown in Figure 2A . At the phylum level, we observed an enrichment of Bacteroidota(33.10% vs. 27.08%), Proteobacteria(12.99% vs. 9.83%) and Fusobacteriota(0.58% vs. 0.56%, P<0.05), and a remarkable decrease of Firmicutes(48.07% vs. 55.89%, P<0.05) and Actinobacteriota(2.44% vs. 5.83%, P<0.01) in PTB patients compared with HCs ( Figure 2B ). At the family level, the relative abundance of Bacteroidaceae (24.59% vs. 16.14%, P<0.05), Tannerellaceae (2.88% vs. 0.90%, P<0.05) and Oscillospiraceae (2.19% vs. 1.18%, P<0.05) increased significantly in PTB patients, while Bifidobacteriaceae (1.7% vs. 5.03%, P<0.01), Butyricicoccaceae (0.45% vs. 0.81%, P<0.05) and Ruminococcaceae (10.74% vs. 22.20%, P<0.001) decreased drastically ( Figure 2B ). At the genus level, compared with the HCs, the PTB patients had a significantly higher relative abundance of Bacteroides (25.59% vs. 16.14%, P<0.05), Parabacteroides (2.88% vs. 0.90%, P<0.05) and Veillonella (4.24% vs. 0.84%, P<0.05), and lower relative abundance of Faecalibacterium (5.31% vs. 14.80%, P<0.05), Bifidobacterium (1.67% vs. 5.02%, P<0.05), Agathobacter (0.86% vs. 3.32%, P<0.05) and CAG-352 (0.60% vs. 2.21%, P<0.05) ( Figure 2B ). Our data indicated the altered gut microbiota composition in PTB patients.

LEfSe analysis demonstrated that many key taxa were clearly different between PTB patients and HCs. According to the LDA score, the significantly enriched taxa in the gut microbiome of the PTB patients were Lactobacillus, Blautia, Erysipelatoclostridium and Fusobacterium. While the abundance of Faecalibacterium, Bifidobacterium, Butyricicoccus, Roseburia and Lachnospiraceae_NK4A136_group were higher in HCs ( Figure 2C ). The structure of the intestinal flora was determined by dynamic interactions between these community members. Based on the relative abundance of OTUs, PTB patients showed a more complex network of interactions than the HCs, with more positive and negative correlations among the microbiota ( Figure S1 ).

A heat map was used to examine the associations between intestinal bacteria and the biochemical index ( Figure 3 ). TC, HDL, LDL and Hb were positively correlated with 5 genera that were enriched in HCs, namely Lachnospiraceae_ND3007_group, Haemophilus, Faecalibacterium, Roseburia, Lachnospira (P<0.05). The levels of TG, Hb, LDL and TC were negatively associated with Lactobacillus, Erysipelatoclostridium and Megasphaera (P<0.05). Particularly, Alb was positively correlated with SCFA producers which were enriched in HCs, such as Lachnospiraceae_NK4A136_group, Faecalibacterium, Dorea, Butyricicoccus and Bifidobacterium (P<0.001), and negatively associated with Lactobacillus, Erysipelatoclostridium, Megasphaera and Sellimonas, most of which were enriched in PTB patients. Moreover, Cr was positively correlated with some bacteria enriched in HCs, and negatively correlated with bacteria enriched in PTB group (P<0.05). Besides, WBC had no direct relation with the genera in our study. It is also unusual in that Glu only had a negative correlation with Incertae_Sedis, Faecalibacterium, [Eubactrium]_eligens_group, Collinsella, Ruminococcus, CAG-352, Bilophila, Akkermanisia and Clostridia_UCG-014 (P<0.05). Taken together, there was an intimate relationship between the deferential genera of PTB patients and the biochemical index, indicating that the altered gut microbiota may be the key pathophysiology of PTB.

PICRUSt2 was used to identify the metabolic and functional changes in the fecal microbiota between the PTB patients and the HCs. The results of PICRUSt2 based on the KEGG classification showed that the predominant predicted bacterial functions were related to cellular processes, genetic system, metabolism and organismal system ( Figure 4A ). Top 20 pathways were shown in Figure 4B , among which 18 were overrepresented in the PTB patients, while only 2 were enriched in healthy controls. Compared with HCs, 7 biosynthesis pathways related to ubiquinone, polyketide, siderophore, folate, lipopolysaccharide, N-glycan and steroid hormone were highly represented in PTB patients. Interestingly, 5 metabolism pathways related to lipoic acid, ascorbate, taurine, fructose and biotin were strikingly increased in PTB patients. Moreover, the degradation capacities for different types of substrates, such as amino acid and glycan were increased in patients, while the dioxin degradation pathway was increased in HCs. Notably, PTB patients had more pathways involved in toxin degradation and energy circulation, specifically reflected in the rise of peroxisome and TCA cycle pathways. However, HCs obviously had more proteasome, indicating a superior protein utilization capacity ( Figure 4B ). Thus, these results further confirmed that Mycobacterium tuberculosis infection would disturb the functional of the fecal microbiota and might participate in the pathogenesis and development of PTB.

To further explore the diagnostic value of the gut microbiota for PTB, we constructed a random forest model to accurately distinguish PTB patients from HCs. 30 genus markers were selected from the model. The training of the random forest classifier involved the calculation of the feature importance of these 30 genera and their relative abundances ( Table S1 ). The vital genera identified based on the model were described based on a mean decrease in accuracy ( Figure 5A ). ROC analysis was used to verify the diagnostic ability of these biomarkers. The area under the curve (AUC) obtained using the training datasets was 0.91, showing that patients with PTB could be successfully distinguished from the HCs. The AUC values for six genera, Lactobacillus (AUC=0.758), Faecalibacterium (AUC=0.803), Roseburia (AUC=0.829), Dorea (AUC=0.851), Monnoglobus (AUC=0.867) and [Eubacterium]_ventriosum_group (AUC=0.88), had an adequate diagnostic efficacy, respectively. Importantly, the model including all six genera had a relatively better diagnostic ability (AUC=0.90) ( Figure 5B ), indicating the combination of six genera may be a set of new diagnostic biomarkers for PTB.