The study was a retrospective study that included all the culture-positive patients diagnosed with RR-TB at local TB dispensaries in Ningbo, China from Jan 1, 2021, to Dec 31, 2023. Patients aged above 18 years with sputum culture-positive, pulmonary, RR-TB were assessed for eligibility, and those with consent to standardized RR-TB regimen were included. Patients were excluded if they were pregnant or infected with HIV, hepatitis B or C virus, or refused to participate. Records related to demographics, clinical and microbiology were retrieved from the national TB information management system.

The response to treatment in this study was evaluated by: 2-month sputum culture conversion as a marker of early treatment response, 6-month culture conversion (previously reported to be predictive of treatment outcome). Sputum culture conversion was defined as two consecutive negative cultures of samples taken at least 30 days apart (Mirzayev et al., 2021). This definition was strictly maintained for all timepoint evaluations (2/6 months).

Drug susceptibility tests of four first-line anti-TB drugs and five second-line drugs were conducted on solid media (Lowenstein-Jensen) based on WHO recommendations (Bagcchi, 2023). The drug concentrations are isoniazid (INH) 0.2 µg/ml, rifampicin (RIF) 40 µg/ml, ethambutol (EMB) 2 µg/ml, streptomycin (SM) 4 µg/ml, levofloxacin (LVX) 2 µg/ml, amikacin (AM) 30 µg/ml, capreomycin (CM) 40 µg/ml, prothionamide (PTO) 40 µg/ml, and para-aminosalicylic acid (PAS) 1 µg/ml. H37RV strains were used as a reference for quality control.

MTB culture products were inactivated, and genomic DNA was isolated using a bacterial DNA extraction kit (QIAGEN Inc., Dusseldorf, Germany), according to the manufacturer’s instructions. The isolated and purified DNA products were transported via a cold chain to a sequencing facility. The purified genomic DNA was quantified using a TBS-380 fluorometer (Turner BioSystems Inc., Sunnyvale, CA, USA) to ensure that the DNA met the quality requirements for library preparation, sequencing, and detection. At least 1 μg of genomic DNA per sample was used as the input material for DNA sample preparation. The DNA samples were treated and fragmented to a size of ~400 bp. Sequencing libraries were generated using the NEXTflex™ Rapid DNA-Seq Kit. The prepared library was multiplexed and loaded on Illumina NovaSeq 6000 PE150 system (San Diego, CA92122, USA). Sequencing was carried out using a 2×150 paired-end (PE) configuration. Raw sequencing data were processed using fastp (v0.20.1) to remove adapter sequences and filter out low-quality bases (Chen et al., 2018). High-quality sequence data were then input into Kraken v2 for species identification, and samples identified as other species or with an MTB proportion below 90% were rejected as contaminated samples (Wood et al., 2019). Finally, the sequencing data from the remaining samples were aligned to the H37Rv reference genome (NC_000962.3) using BWA (v0.7.17) (Li, 2013). Samples with an average sequencing depth > 20 × and average genome coverage >95% were selected for subsequent data analysis. The SAMtools/BCFtools suite was used for calling fixed (frequency ≥90%) SNPs at loci where the alternate alleles were supported by at least five reads (including both forward and reverse reads) (Danecek et al., 2021).

Clean sequencing data were input into the local version of TB-Profiler (v6.5.0) with its strict implementation of the V2 catalog (database name: who) to identify the genotype of resistance-associated mutations and detect the resistance profile of 15 anti-TB drugs, including Amikacin (AM), Bedaquiline (BDQ), Capreomycin (CM), Clofazimine (CFZ), Delamanid (DLM), Ethambutol (EMB), Ethionamide (ETO), Isoniazid (INH), Kanamycin, Levofloxacin (LVX), Linezolid (LZD), Moxifloxacin (MFX), Pyrazinamide (PZA), Rifampicin (RIF), Streptomycin (SM). Mutations with a frequency of less than 10% were excluded. WGS-based drug susceptibility testing (DST) results were determined by assessing the presence or absence of resistance-associated mutations in a WHO-recommended database, which classified mutations into Tier 1 (those most likely to confer resistance) and Tier 2 (genes with a reasonable pretest probability of resistance) (Walker et al., 2022). Hetero-resistance was defined based on the frequency of resistant alleles in the sequence reads <99% in this study.

The fixed SNPs, excluding those in SNPs in PE/PPE areas, insertion elements, repetitive regions, and drug resistance-associated genes, were combined into a concatenated alignment (Luo et al., 2014). Maximum-likelihood (ML) phylogenetic trees were inferred from the concatenated alignment using IQ-Tree v2 (Nguyen et al., 2015). The best-scoring ML tree were rooted using M. canettii (RefSeq: NC_015848.1) as the outgroup and visualized with the Interactive Tree of Life (iTOL) (Letunic and Bork, 2021).

Statistical analyses were performed in R package gtsummary (Sjoberg et al., 2021). To compare multiple categorical variables, binary logistic regression was employed where appropriate, with results reported as odds ratios (ORs) and 95% confidence intervals (CIs). Variables showing a significant association (P < 0.05) in univariate logistic regression were included in the multivariate analysis. Forward stepwise logistic regression was then performed to identify whether these statistically significant covariates were independently associated with SCC.

A Random Forest (RF) model was implemented in R v4.0 with a fixed seed of 27 to identify genetic mutations and clinical factors associated with delayed SCC at 2 and 6 months in RR-TB. RF was chosen for its ability to handle high-dimensional, nonlinear data, robustness to overfitting, and capacity to rank feature importance, making it well-suited for integrating clinical and genomic predictors compared to other models like logistic regression or support vector machines. Mutations in DR related genes and clinical metadata from 150 isolates were merged. Clinical and phenotypic DST features were encoded, and mutations imputed with 0. Mutations in ≥5% of samples were selected. Class imbalance was addressed using SMOTE (DMwR package) (Kumari et al., 2020). The RF model was trained using the caret package’s train function (randomForest method, 100 trees) with an 80/20 train-test split and 5-fold cross-validation (ROC-AUC). Hyperparameter tuning was performed with tuneLength = 3, optimizing the number of variables sampled at each split (mtry). Feature importance in the RF models was assessed using mean decrease in accuracy (MDA), with standard deviations estimated from 100 bootstrap iterations to quantify variability. Test-set AUC and its standard deviation were estimated using 100 bootstrap iterations. Feature importance and ROC curves were visualized using ggplot2 (Villanueva and Chen, 2019).

Between January 2021 and September 2023, a total of 2,704 patients were diagnosed with sputum culture-positive tuberculosis, among whom 225 were identified as having RR-TB. From 177 initially enrolled patients, 48 were excluded: 11 lost to follow-up pre-treatment, 22 declined treatment, 2 died pre-enrollment, and 13 transferred. During WGS processing, 27 additional cases were excluded: 17 from culture failure, six due to contamination (<90% MTBC sequences), and four confirmed as NTM infections (three M. intracellulare, one M. abscessus) by Kraken2 analysis. Ultimately, 150 patients with successful WGS results were included in the final analysis ( Figure 1 ). Among the 150 RR-TB patients, 97 patients (64.7%) achieved SCC after 2 months of treatment. By 6 months of treatment, the number of patients who successfully achieved SCC increased to 123 (82.0%).

As shown in Table 1 , the mean age of patients was 45 years, and 117 (78.0%) were male. A total of 72 patients (48.0%) were smear-positive. Retreatment cases accounted for 44.0% ( n = 66 ). Underlying health conditions were reported in 62 patients (42.7%), with diabetes (53.2%, 33/62) and hypertension (35.5%, 22/62) being the most common. Additionally, 24 patients were smokers, and 21 reported alcohol consumption. Migrant individuals comprised 46.0% ( n = 60 ) of the cohort. Regarding occupation, the largest groups were workers (32.7%, n = 49 ), unemployed individuals (29.3%, n = 44 ), and farmers (20.0%, n = 30 ).

Excluding RIF, phenotypic DST revealed that resistance rates for the remaining eight anti-tuberculosis drugs ranged from highest to lowest as follows ( Table 1 ): isoniazid (INH), 74.0% (111/150); streptomycin (SM), 62.0% (93/150); levofloxacin (LVX), 36.0% (54/150); ethambutol (EMB), 33.3% (50/150); para-aminosalicylic acid (PAS), 16.7% (25/150); prothionamide (PTO), 15.3% (23/150); amikacin (AM), 8.7% (13/150); and capreomycin (CM), 6.7% (10/150).

A total of 133 mutations were identified by WGS across 19 antibiotic resistance genes related to RIF, INH, EMB, PZA, ETO, SM, PAS, aminoglycosides, and fluoroquinolones resistance ( Supplementary Table S1 ). Genomic analysis revealed RIF-resistance-associated mutations in 98.7% (148/150) of isolates, predominantly in rpoB (146 strains), with rare rpoC mutations (2 strains). The most frequent rpoB mutation was Ser450Leu, accounting for 58.2% (85/146) of cases, followed by Leu452Pro (12.3%, 18/146), with other low-frequency mutations (e.g., Leu430Pro, Asp435Ala). Additionally, two strains carried rpoC mutations (Ile435Thr, Phe452Ser). Beyond RIF, mutations conferring resistance to other drugs were detected, including embB_p.Met306Ile (EMB), katG_p.Ser315Thr (INH), gyrA_p.Asp94Gly (fluoroquinolones), and rrs_n.1401A>G (aminoglycosides), reflecting the multidrug-resistant nature of the cohort. Based on WGS results and WHO’s updated definitions, the 150 RR-TB strains were further categorized into 65 MDR-TB, 51 pre-XDR-TB, and 34 RR-TB cases. These findings highlight the genetic diversity of resistance in RR-TB and underscore the utility of WGS in identifying critical mutations driving treatment challenges.

To elucidate the evolutionary relationships of the 150 RR-TB strains, we constructed a ML phylogenetic tree based on concatenated sequences from non-redundant SNP loci ( Figure 2 ). Genotyping analysis revealed that 90.0% (135/150) of the strains were classified into lineage 2 (L2), with the majority falling under the L2.2.1 ( n = 125 ) and L2.2.2 ( n = 9 ) sublineages (Beijing family), accounting for a total of 134 strains. The remaining 10.0% (15/150) belonged to lineage 4 (L4), comprising sublineages L4.4 ( n = 7 ), L4.5 ( n = 6 ), and L4.2 ( n = 2 ), respectively.

Univariate analysis of 150 RR-TB patients identified factors associated with the SCC failure at 2- and 6-month post-treatment initiation ( Table 2 ). At 2-month, SCC failure showed significant associations ( O R = 3.15 ,   95 %   C I :   1.58 − 6.46 ,   P = 0.001 ) with smear-positive status, phenotypic resistance to (INH O R = 3.24 ,   95 %   C I :   1.38 − 8.55 ,   P = 0.011 ), AM ( O R = 7.29 ,   95 %   C I :   2.11 − 33.7 ,   P = 0.004 ), CM ( O R = 8.44 ,   95 %   C I :   2.02 − 57.5 ,   P = 0.009 ), LVX ( O R = 4.53 ,   95 %   C I :   2.23 − 9.45 ,   P < 0.001 ), and genotypic pre-XDR-TB ( O R = 4.70 ,   95 %   C I :   1.80 − 13.5 ,   P = 0.002 ). However, by 6-month, only smear-positive status and LVX resistance remained significantly associated with delayed conversion ( O R = 3.25 ,   95 %   C I :   1.39 − 7.86 ,   P = 0.007 ), indicating that other resistance profiles became less influential over time.

To identify independent predictors of delayed SCC, we constructed multivariable logistic regression models. For 2-month conversion failure, the model was adjusted for smear-positive status, phenotypic resistance to INH, AM, CM and LVX, as well as genotypic pre-XDR-TB. For 6-month failure, adjustments included smear-positive status and LVX resistance. After adjustment, smear positivity ( a O R = 3.32 ,   95 %   C I : 1.54 − 7.44 ,   P = 0.003 ) and LVX resistance ( a O R = 6.09 ,   95 %   C I : 1.31 − 35.0 ,   P = 0.026 ) were significantly associated with 2-month SCC failure. Other factors in the 2-month model, including INH, AM, CM resistance, and genotypic pre-XDR-TB were not statistically significant. For 6-month SCC failure, smear positivity ( a O R = 2.90 ,   95 %   C I : 1.19 − 7.61 ,   P = 0.023 ) and LVX resistance ( a O R = 3.01 ,   95 %   C I : 1.27 − 7.40 ,   P = 0.013 ) remained significant predictors. The strong association of smear positivity and LVX resistance with SCC failure suggested that higher bacterial burden and fluoroquinolone resistance significantly impede early and sustained treatment response.

The RF model demonstrated robust performance in predicting SCC failure at both 2 and 6 months among patients with RR-TB, yielding ROC-AUC values of 0.86 ± 0.06 and 0.763 ± 0.103, respectively, based on 5-fold cross-validation ( Figures 3A, B ). Feature importance analysis measured as MDA, highlighted key determinants of delayed SCC with variability estimated from 100 bootstrap iterations ( Figures 3C, D ). For the 2-month endpoint, the most influential predictors included LVX resistance ( M D A   = 6.37   ±   1.05 ), AM resistance ( M D A = 5.95   ±   0.46 ), and the embB_p.Met306Ile mutation (EMB resistance, M D A = 5.94   ±   0.47 ), followed by rrs_n.1401A>G (AM resistance, M D A = 5.83   ±   0.30 ) and INH resistance ( M D A = 5.61   ±   0.58 ). At 6 months, the top contributors were LVX resistance ( M D A = 5.46   ±   1.33 ), baseline sputum smear positivity ( M D A = 5.00   ±   1.50 ), and the katG_p.Ser315Thr mutation (INH resistance, M D A = 4.85   ±   0.92 ), with genotypic DR type ( M D A = 4.16   ±   0.96 ) and gyrA_p.Asp94Gly (fluoroquinolone resistance, M D A = 3.43   ±   0.62 ) also playing significant roles. Notably, katG_p.Ser315Thr and gyrA_p.Asp94Gly emerged as stronger predictors at the 6-month mark, underscoring their association with prolonged treatment failure due to INH and fluoroquinolone resistance, respectively.