In this retrospective, observational study, we used the BioBank of Tuberculosis System, which collected standardized data of all MTB isolates from culture-positive cases, who sought healthcare in the Beijing Chest Hospital, a national TB-designated hospital in China. Eligible patients met all of the following criteria: (i) completion of anti-TB treatment following the guidelines of National Tuberculosis Control Program; (ii) successful outcomes indicated by sputum culture conversion plus clinical findings; (iii) two or more episodes of TB which occurred in the period between January 2008 and December 2019; and (iv) a minimum time interval of 12months between two episodes (Zong et al., 2018). For patients experiencing multiple recurrences, only the first recurrent episode was included in our analysis. Data on TB cases in our cohort came from the Electronic Medical Record System in our hospital. The sociodemographic and clinical characteristics were analyzed in this study: age, sex, ethnicity, local authority of residence, TB recurrence, previous TB treatment, and site of TB disease. The study protocol was approved by the Ethics Committee of Beijing Chest Hospital, Capital Medical University and written informed consent was obtained from each participant prior to the study.

All MTB isolates were stored in 7H9 medium supplemented with 25% glycerin at −80°C refrigerator. Prior to DNA extraction and in vitro drug susceptibility testing (DST), the frozen bacteria were subcultured into Löwenstein Jensen (L-J) medium. After 4weeks of incubation at 37°C, the fresh colonies were scraped from the surface of L-J medium and inactivated in a water bath at 85°C for 30min. The bacteria in the pellet were subjected to DNA extraction by using the MasterPure™ Complete DNA and RNA Purification Kit (Epicenter, Madison, WI, United States) following the manufacturer’s instructions. The quality of the extracted DNA was assessed by NanoDrop 2000c spectrophotometer (Thermo Fisher Scientific, United States). A single absorbance peak at 260nm plus a 260/280 absorbance ratio of 1.8–2.0 indicated high-quality intact DNA. In addition, the concentration was determined by Qubit 2.0 fluorometer (Invitrogen, Thermo Fisher Scientific, United States). One microliter of purified DNA was added to reaction mixture for species identification (Zeesan Biotech, Xiamen, China). A decontamination at 50°C for 2min was followed by initial denaturation at 95°C for 10min, 55cycles of denaturation at 95°C for 15s, annealing at 57°C for 20s, and plus a final extension at 78°C for 20s. The amplicons were analyzed by multicolor melting curve analysis as previously reported (Xu et al., 2019; Tan et al., 2021).

Because the previous DST results were produced by the absolution concentration method which was not endorsed by the World Health Organization (Organization, W.H, 2018), the fresh MTB cultures were collected for phenotypic DST with the BACTEC MGIT system (BD, Franklin Lakes, NJ). Briefly, the suspension was prepared to a turbidity equivalent to that of 0.5 McFarland standard; 0.5ml of a 1:5 dilution was used for inoculation into the MGIT containing drugs, and 0.5ml of a 1:500 dilution was used for inoculation into the drug-free control MGIT tube. After 6–14days incubation in the MGIT 960, the results were reported automatically by comparing the number of growth units (GUs) in the control tube (GU=400) and the number of GUs in the tubes containing drugs. The concentrations of drugs were as follows: isoniazid 0.1μg/ml, rifampicin 1.0μg/ml, ethambutol 5.0μg/ml, streptomycin 1.0μg/ml, levofloxacin 1.0μg/ml, moxifloxacin 0.25μg/ml, amikacin 1.0μg/ml, and capreomycin 2.5μg/ml. Isolates resistant to isoniazid and rifampicin were defined as multidrug-resistant tuberculosis (MDR-TB).

Whole-genome sequencing was performed using the TIANSeq direct rapid DNA library construction kit (Tiangen, China) and the NovaSeq Illumina sequencer (Illumina, United States). After quantification and quality inspection, the extracted genomic DNA was broken into nucleic acid fragments with a length of 200–300bp by an ultrasonic interrupter (Covaris, United States) with Peak Incident Power: 175; duty factor: 10%; cycle/burst: 200; temperature: <8°C; and time: 120s. The 5′ and 3′ protruding ends of the broken DNA fragments were flattened and connected to adaptors by T4 DNA ligase. Fragment selection and purification are performed by magnetic bead purification. An amplification method was used for enrichment and magnetic bead purification was used for purification to obtain a DNA library for sequencing. Qubit Flurometer 3.0 fluorometer (Thermo Fisher Scientific, United States) and Agilent 2100 (Agilent, United States) were used to measure the concentration and length distribution of the library, respectively. The paired-end sequencing program (PE150) for the qualified library was performed on the NovaSeq Illumina sequencing platform. FastQC quality control software was used for quality control and removal of low-quality data in Raw Data, Quality filter: q, qualified_quality_phred, set the base quality value to be no less than a qualified base, the default base quality value is 15, that is, the default base quality ≥15 is qualified base, <15 is unqualified Base. u, unqualified_percent_limit, set the allowable percentage of unqualified bases, remove this read, the default is 40, that is, when the default percentage of unqualified bases >40%, remove the read; Q, disable_quality_filtering, set this parameter to disable the default quality filtering parameters (q, u). Length filter: l, length required, set the minimum length of read, the default is 15, that is, reads with length <15 are removed; L, length limit, sets the maximum length of read, the default is 0, there is no maximum length limit; and Low complexity filtering: Y, complexity threshold, set the complexity filter threshold of read, the default is 30. The MTB whole-genome sequence analysis platform (SAM-TB) software system was used to compare the sequencing data with MTB H37Rv whole-genome sequences and evaluate the sequencing depth and coverage. FastQC was used for quality control of SAM-TB, via adopting Cutadapt (v1.15) software to remove linkers and low-quality sequences. The retained high-quality sequence was aligned to the reference genome H37Rv (NC_000962.2) using BWA (v0.7.15). After removing duplicate reads with Picard (v2.0.1), we used SAMtools (v1.6)/VarScan (v2.3.6) to identify SNPs and short indels, and used Delly (v0.8.7) software to identify larger deletions (≥50bp). By comparison with the reference genome H37Rv, a series of information, such as single-nucleotide polymorphism (SNP) and insertion–deletion (InDel) of the target genome, were obtained and a mutation matrix of all the mutation sites of each strain was generated. SNP distance matrix, cluster strains, and genetic distance were analyzed according to the mutation matrix. SNP differences were determined between the pairs of isolates using WGS. In addition, PE/PPEs/PGRs genes, insert sequences, and other highly repetitively related positions are filtered. Reads of these regions are usually mapped to the wrong genomic coordinates, which will lead to misrecognition of SNPs. Therefore, these regions should not be considered in SNP calling. and recurrence was classified as being due to relapse if there was a difference of ≤6 SNPs (Cock et al., 2010; Jiang et al., 2011; Chen et al., 2020; Shanmugam et al., 2021). A minimum spanning tree was constructed based on pairwise whole-genome SNP differences using BioNumerics (Version5. 0, Applied Maths, Sint-Martens-Latem, Belgium) software (Mazuet et al., 2017; Rouard et al., 2019). Genotypic drug resistance prediction was predicted using the TB-Profiler (Faksri et al., 2019). The average variant depth for each isolate was ≥10 and the variant calls for drug resistance were based on the algorithms of TB-Profiler. Resistance was predicted for a total of 14 anti-TB drugs, including rifampicin, isoniazide, ethambutol, pyrazinamide, and streptomycin, and the second-line drugs Linezolid, clofazimine, bedaquiline, fluoroquinolone, capreomycin, amikacin, kanamycin, PAS, and ethionamide.

For the patients’ information, Epidata 3.0 software with dual-entry method was used to establish the database, and SPSS 22.0 (SPSS Inc.)¹ statistical software was used to perform statistical analysis on the data. For continuous variables that were normally distributed, the mean±standard deviation (SD) was used; for data that were not normally distributed, the median [interquartile range; M (IQR)] was used. Counting data were described in terms of rate or composition ratio (%). Univariate and multivariate unconditional logistic regression models were used to determine the influencing factors of patients’ recurrence. Wilcoxon rank sum test was used to evaluate gene mutation rate. The Kaplan–Meier curve was generated to describe and compare the time to relapse between the different groups. p<0.05 was considered statistically significant.

Between January 2008 and December 2019, a total of 89 patients meeting the inclusion criteria were included in our analysis. The majority of patients was adult men with a median age of 51years (range: 16–86years). Of the 89 subjects in the study, diabetes (n=29, 32.6%) was the most frequently observed comorbidity, followed by chronic disease (n=19, 21.4%), malnutrition (n=16, 18.0%), and anemia (n=14, 15.7%; Table 1). Figure 1 illustrates a flowchart of the enrolled samples. Among 89 recurrent TB patients, 12 (13.5%) were excluded due to subculture failure of any paired isolate. Isolates from the remaining 77 cases were analyzed by the WGS to distinguish relapse and reinfection. Of 77 cases, nine were excluded from analysis due to identification of nontuberculous mycobacteria (NTM) by WGS and multicolor melting curve analysis, leaving 68 patients with paired patients in final analysis.

Lineage 2 was predominant among recurrent TB patients in the present study, accounting for 97.0% (132/136) of MTB strains, followed by Lineage 4 (3.0%, 4/136, Supplementary Figure S1). A comparison with SNP differences is shown in Figure 2. Two main groups were identified: 62 pairs of isolates had six or fewer SNP differences, and six pairs of samples had a dramatically higher number of SNP differences (range: 15–360). Ultimately, of the confirmed recurrence, 62 (91.2%) were due to relapse (Supplementary Table 1), whereas the remaining six (8.8%) were due to reinfection. Of 62 relapses, 35 (56.5%) occurred within 12months of treatment completion and 54 (87.1%) within 24months (Figure 3). We further assessed the time intervals to TB relapse stratified to drug susceptibility profiles and comorbidities. As shown in Figure 4, there was a strong association between earlier relapse and underlying chronic diseases (p=0.033). In addition, no significant differences were noted between MDR- and non-MDR-TB groups.

Interestingly, there was one isolate that mapped closely to other patients isolates on the tree, and this merited closer attention to see if there were nosocomial infection or community transmission. As demonstrated in Figure 2, the isolate of second episode of Patient 41 showed high sequence similarity to that of Patient 42, and only one SNP was found between two isolates from different individuals. However, the two isolates of patient 41 differed by 75 SNPs, suggesting recurrence due to reinfection. We firstly reviewed their hospitalization records. Nosocomial infection was excluded since no overlapping hospitalization period occurred. Of note, both two patients lived in the neighboring blocks. Thus, we conducted interview to determine whether the community transmission could contribute to this cluster. Epidemiological data revealed that these two retired elderly individuals routinely completed morning exercises in the same park. This combined with their genotyping data supported the potential community transmission between Patient 41 and Patient 42 (Figure 5).

We further analyzed the rate of microevolution of MTB from the first and second episodes of the longitudinally sampled patients (Figure 6). A rate of change in DNA sequences of 0.40 SNPs per genome per year [95% confidential interval (CI) 0.24–0.56] was noted in our study. We recorded a significant higher mutation rate in non-diabetic patients (0.53 SNPs per genome per year) compared with that in diabetic patients (0.14 SNPs per genome per year, p =0.012). By contrast, no difference was observed in mutation rate between MDR-TB and non-MDR-TB patients (0.39 SNPs per genome per year for MDR-TB vs. 0.42 SNPs per genome per year for non-MDR-TB, p =0.919).

We further analyzed the differences in drug susceptibility profiles between first and second TB episode. As shown in Figure 7, levofloxacin resistance was the most frequently observed for 8 (12.9%) of 62 relapse cases, followed by rifampicin (9.7%, n =6), streptomycin (9.7%, n =6), and amikacin (9.7%, n =6). For moxifloxacin, another fluoroquinolone (FQ), 2 (3.2%) paired isolates exhibited changes from sensitivity to resistance. Statistical analysis revealed that this proportion was significantly lower than that for levofloxacin (p =0.048), indicating that MTB isolates were more prone to develop levofloxacin resistance than moxifloxacin resistance after FQ exposure. A deep analysis was conducted to compare the drug susceptibility profiles between two episodes by WGS. As shown in Supplementary Table 2, FQ resistance was recorded in 22.6% (14/62) of relapse cases, followed by kanamycin (11.3%, 7/62), and pyrazinamide (8.1%, 5/62). In contrast, no bedaquiline and linezolid resistance were noted in these relapse cases.