Inferring Mycobacterium Tuberculosis Drug Resistance and Transmission using Whole-genome Sequencing in a High TB-burden Setting in China

Tuberculosis (TB) remains a global threat to human health disproportionately affecting the poor and vulnerable populations. China ranked second in the global TB burden in 2020^[1,2]. The Xinjiang Autonomous Region has one of the most serious TB epidemics in China^[3], especially in the southern region of Xinjiang^[4]. Among the provinces wherein TB prevalence surveys were conducted after 2010, its prevalence in Xinjiang (430/100,000) was the highest^[5]. Of the reported cases in Xinjiang, 70% were from the southern region^[6] with Kashgar being one of the most densely distributed regions^[7].

In high-burden areas, transmission severely exacerbates the TB epidemic^[8]. In addition, high burden regions are plagued by the transmission of drug-resistant (DR) strains^[9]. In the past, it was thought that DR-TB was caused by the development of resistance during antituberculosis treatment (acquired resistance); however, recent studies have shown that the direct transmission of DR strains between individuals (transmitted resistance) is also significant^[10,11]. Transmission dynamics are heterogeneous in different regions^[12], and understanding local transmission patterns can help to develop effective prevention and control strategies. For example, in San Francisco, molecular epidemiological studies began after 1990 to identify clustered strains, high-risk groups for recent transmission (HIV-infected people and African Americans), and transmission sites (hospitals and bars). Targeted control measures were implemented, thereby resulting in a significant decrease in the incidence of TB (from 51.2/100,000 in 1992 to 29.8/100,000 in 1997)^[13-16].

The complex socioeconomic and ethnic composition of Xinjiang highlights the need to study TB transmission and drug resistance in this region. Therefore, we analyzed the whole-genome sequences of 1,668 clinical strains from the region in 2020 to reveal drug resistance characteristics, identify risk factors for recent transmission, and infer the proportion of transmitted drug resistance. This will provide necessary information on the TB situation in the region to assist in control activities to suppress the disease.

Kashgar Prefecture is located southwest Xinjiang, bordering Tajikistan, Afghanistan, and Pakistan, and is close to North India. The prefecture is divided into 11 counties and one county-level city. TB is diagnosed and treated at designated medical institutions (People’s Hospitals in each county and the TB dispensary in the county-level city). According to the Chinese National guidelines^[17], each patient who visited the hospital, underwent smear microscopy and sputum culture. Culture-positive strains were transported uniformly to the TB dispensary of Kashgar Prefecture for strain preservation, and the Chinese Center for Disease Control and Prevention collected patient strains from the dispensary in 2020. Demographic and clinical information of patients with TB were submitted to the National TB Reporting System. Unfavorable outcomes included treatment failure, death, and lack of follow-up.

Genomic DNA was extracted using the cetyltrimethylammonium bromide (CTAB) method^[18]. A 150 bp fragment length library was constructed for each sample, and sequencing was performed on an Illumina Novaseq6000 platform (Illumina Inc. San Diego, CA, USA), using a paired-end strategy. We used a suite of validated software for SNP calling^[19,20]. Each SNP was supported by at least five reads, with at least one read in each sequencing direction. Strains with a sequencing depth below 20× or a mapping coverage of less than 95% were also excluded.

Whole-genome sequencing has been proposed as the ultimate molecular diagnostic test for detecting all molecular determinants of resistance^[21]; therefore, we used whole-genome sequencing (WGS) to detect drug resistance in strains. The drug resistance profiles were predicted based on known drug resistance-conferring mutations^[22] for 14 anti-TB drugs. Pairwise SNP distances were calculated based on fixed mutations (frequency ≥ 75%)^[23,24]. Strains were grouped into genomic clusters based on the criterion that the maximum pairwise distance between at least two strains should be 12 or fewer single nucleotide polymorphisms (SNPs). Single SNPs related to drug resistance mechanisms were excluded when calculating pairwise distances.

Lineages of each strain were identified as described by Coll et al.^[25]. Phylogenetic trees were constructed using MEGA X^[26] with the maximum likelihood method. The trees were visualized using an Interactive Tree of Life (https://itol.embl.de/).

Transmitted resistance included both resistant strains from recent transmissions and resistant strains isolated from new cases. Resistant strains of recent transmission were screened from drug-resistant strains in genomic clusters using the location of DR mutations on the phylogenetic tree: first, mutations that were inherited from an inferred common ancestor indicated recent transmission, whereas mutations on terminal tips were independently acquired^[27]. Second, for clustered DR strains, as long as one resistance mutation was caused by recent transmission, the strain was considered to be resistant because of recent transmission (if all resistant mutations were acquired independently, the strain was not considered to be resistant because of recent transmission). This allowed the exclusion of strains that had independently acquired resistance mutations after recent transmission, as opposed to simply determining whether the DR strains were clustered^[28]. For resistant strains that were not clustered, or clustered but not resistant because of recent transmission, transmitted resistance was considered if the strain was from a new case^[23,28]. When calculating the proportion of transmitted resistance, strains that could not be identified owing to a lack of information on treatment history were excluded.

Categorical data were described as proportions (%). For categorical data, differences between groups were tested using the chi-square test or corrected chi-square test. Fisher's exact probability test was used when conditions of the chi-square test were not met. The Kruskal–Wallis test was used to compare the three groups of continuous non-normal data. Bonferroni-adjusted alpha levels were used for two-by-two comparisons. Risk factors were analyzed using univariate and multivariate logistic regression. Factors with a P-value less than 0.2 in the univariate regression were entered into the multivariate analysis. Patients with missing demographic and clinical information were excluded from the logistic analyses. All tests were two-sided, and statistical significance at 95% confidence level was defined as P < 0.05. All statistical analyses were performed using SPSS V25.0 (IBM Corp., Armonk, NY, USA). The geographical distribution of culture-positive TB cases was mapped using ArcGIS (version 10.8.1). Spatial aggregation was analyzed using kernel density estimation methods with Gaussian smoothing.

A total of 11,155 patients with pulmonary TB were reported from January 1 to December 31, 2020, of which 1,869 were culture-positive. Of the 1,869 patients, 201 were excluded because of infection with non-tuberculosis mycobacteria (NTM), contaminated strains, or ineligible whole-genome sequences, leaving 1,668 patients (Figure 1). For the 1,630 patients with demographic and clinical information (such information was lacking for the remaining 38 patients), the median age was 62 years (interquartile range, 42–71 years). Half of the patients (50.0%, 815/1,630) were male and the majority (81.7%, 1,331/1,630) were farmers (Table 1).

Figure 1.  Patient enrollment and study flowchart. NTM, non-tuberculosis mycobacteria; PTB, pulmonary TB.

The proportion of previously treated cases in Kashgar was 20.9% (340/1,630). Compared with new cases, previously treated cases tended to be older than 65 years (38.8% of new cases vs. 56.2% of previously treated cases, P < 0.001). In addition, they were associated with a total delay (period from symptom onset to treatment initiation) of > 12 weeks (9.8% vs. 16.5%, P = 0.003) and unfavorable treatment outcomes (4.7% vs. 9.7%, P < 0.001) (Table 1).

The whole genomes of 1,668 Mycobacterium tuberculosis (MTB) strains isolated from patients were sequenced with an average depth of 215.6 × (23.5× – 476×) and a coverage of 98.6% (97.8%–99.2%). Based on the phylogenetic analysis, these isolates were grouped into three lineages: 46.0% (768/1,668) belonged to lineage 2 (All were Beijing strains), 27.5% (458/1,668) belonged to lineage 3, and 26.5% (442/1,668) belonged to lineage 4 (Figure 2).

Figure 2.  The phylogenetic tree of 1,668 MTB strains isolated in Kashgar. The different colors on the branches indicate different lineages. The inner circle (marked with 1) indicates genomic clustered strains. The outside circle (marked with 2) indicates the drug resistance pattern.

Analysis of the genome sequences for drug resistance-conferring mutations showed that 13.4% (223/1,668) of the strains were resistant to at least one of the 13 drugs (Table 2). Thirty resistance profiles were identified, with the highest proportion of cases being resistant to only one drug (8.2%, 136/1,668). (Supplementary Table S1, available in www.besjournal.com). The drug with the highest single-drug resistance rate was isoniazid (7.4%, 124/1,668). Streptomycin ranked second (6.0%, 100/1,668). The third was rifampicin (3.3%, 55/1,668). The rate of rifampicin resistance was 1.8% (23/1,290) in the new cases and 9.4% (32/340) in the previously treated cases. The multidrug resistance rate was 0.8% (10/1,290) in newly diagnosed cases and 6.2% (21/340) in previously treated cases.

DR mutations were more frequently detected in lineage 2 (18.6%, 143/768) than in lineages 3 (8.7%, 40/458) and 4 (9.0%, 40/442; χ² = 33.899, P < 0.001). Resistance rates were not statistically different between lineages 3 and 4 (Supplementary Table S2, available in www.besjournal.com).

The clustering rate is used in epidemiological studies on TB to indicate recent transmission levels. All the strains formed 164 genomic clusters (each cluster contains 2 to 13 strains). The clustering rate was 25.9% (432/1,668). To screen for factors associated with clustered cases, differences in demographic, clinical, and bacteriological factors were compared between the clustered and non-clustered cases (Table 3). The proportion of previously treated cases among the clustered cases was not significantly different from that among the non-clustered cases, suggesting that previously treated cases were not at an increased risk of TB transmission. The clustering rate of lineage 2 (33.4%, 257/768) was significantly higher than those of lineages 3 (16.4%, 75/458) and 4 (22.6%, 100/442; χ² = 47.013, P < 0.001). No statistical differences were observed in the clustering rates between lineages 3 and 4. Notably, the clustering rate of the DR mutant strains was significantly higher (32.7%, 73/223) than that of the susceptible strains (24.8%, 346/1,412; χ² = 6.268, P = 0.012).

The screened factors were then substituted into the logistic models. Multivariate logistic regression analysis determined the following risk factors for recent transmission (Table 4): the risk of clustering was 1.760 times higher in men than in women [95% confidence interval (CI): 1.370–2.260]. It was more common for patients younger than 35 years to be in clusters than for those older than 36 years [adjusted odds ratio (aOR) = 2.230; 95% CI: 1.677–2.966]. Among the various occupations, the risk of clustering was higher for those engaged in housekeeping or the unemployed (aOR = 1.565, 95% CI: 1.076–2.275) and was highest for students (aOR = 1.887, 95% CI: 1.025–3.474). Beijing strains had a 1.681 (95% CI: 1.323–2.136) increase in the risk of forming genomic clusters compared to the other lineages. Drug resistant strains had a higher risk of forming genomic clusters than susceptible strains (aOR = 1.414, 95% CI: 1.023–1.954).

Kernel density maps were plotted to understand the spatial aggregation of the clustered cases. There were nine spatially concentrated areas (hotspots) of clustered patients, seven of which were also hotspots for nonclustered patients. The remaining two hotspots for clustered patients, one in Bachu County and one at the border between Payzawat and Shule counties, were not hotspots for non-clustered patients.

Resistant strains have a higher risk of clustering; therefore, they may have primarily been caused by transmission. To test this hypothesis, we quantified transmitted resistance based on the recent transmission of resistance mutations and the treatment history of the patients. Of the 73 clustered resistant strains, 66 had at least one resistance mutation from recent transmission. Mutations in the remaining seven strains were acquired independently. These seven strains, as well as 150 non-clustered strains, were uniformly judged based on the patient’s treatment history. We identified 105 strains isolated from new cases, 47 strains isolated from previously treated cases, and 5 strains isolated from cases with missing information on treatment history (these five strains were excluded from the calculation of proportions, thus the denominator was 218). Resistant strains from previously treated cases may have acquired resistance during the course of treatment, whereas resistant strains from new cases usually indicate direct transmission. Therefore, we inferred that the total proportion of transmitted resistance was 78.4% (171/218) (Figure 3). The percentage of resistance transmitted to all anti-TB drugs was higher than 50% (Table 5).

Figure 3.  Flowchart for calculating transmitted resistance. Unknown cases mean the cases with missing information on treatment history.

To the best of our knowledge, this is the first systematic population-based study conducted in Kashgar. We found that the proportion of previously treated cases in Kashgar (20.9%; 340/1,630) was almost twice as high as that reported in the latest national TB epidemiological survey (11.8%)^[29]. The rates of rifampicin resistance in new and previously treated cases were 1.8% and 9.4%, respectively, and the rates of multidrug resistance were 0.8% and 6.2% in the two categories, respectively. The clustering rate was 25.9%, indicating that the prefecture should be aware of the recent transmission of TB. Bacteriological factors (drug-resistant and Beijing strains) and demographic factors (sex, age, and occupation) were associated with recent transmission. The spatial analysis identified seven areas as hotspots shared by both clustered and non-clustered patients, suggesting a denser population distribution in these areas; two areas were hotspots for clustered patients only, thereby suggesting that these areas may be at a high risk for recent transmission. Transmitted drug resistance was 78.4%, suggesting that transmitted resistance is a more significant cause of drug resistance than acquired resistance in Kashgar.

Based on the information on the included cases, we found a high proportion of previously treated cases, which is a major public health concern. A survey of the treatment outcomes of new smear-positive cases in Kashgar between 2016 and 2018 showed a cure rate of 77%^[30], which is below the WHO target of at least an 85% cure rate set in 2005^[31]. This low cure rate may explain the high proportion of previously treated patients. Patients with TB in the prefecture are mainly located in rural areas and have low awareness of TB control information and hazards. Therefore, they are prone to discontinue treatment on their own after the initial symptoms have subsided^[30]. This finding suggests that Kashgar should strengthen the management of patient treatment processes to improve outcomes.

According to the WGS-based drug susceptibility results, we found that the rates of rifampicin resistance in Kashgar (1.8% in new cases and 9.4% in previously treated cases) were lower than those reported in China and globally. The results of the 2018 National Drug Resistance Surveillance published in November 2022 showed that the rates of rifampicin resistance were 5.08% in new cases and 23.31% in previously treated cases, which were more than twice as high as in Kashgar^[32]. WHO published the global estimates of rifampicin-resistant TB as 3%–4% in new cases and 18%–21% in previously treated cases, which were also higher than those in Kashgar^[1]. In addition, the resistance rates of the strains in Kashgar to isoniazid, ethambutol, streptomycin, and fluoroquinolones were all lower than those reported in the 2018 National Drug Resistance Surveillance. This suggests that the Kashgar Prefecture in China has low rates of drug resistance.

Based on the genetic distance between the strains, we found that the clustering rate obtained from a one-year sample was 25.9%. As the majority of patients with recent transmission develop active TB within two years of infection^[33,34], two years is an ideal time to collect samples. The short sample collection time in this study missed some clustered patients; therefore, the actual clustering rate in Kashgar was > 25.9%. To understand the relative levels of recent transmission in Kashgar, we selected recently published population-based WGS studies for comparison. The clustering rate in Kashgar may be higher than that in low-burden countries; for example, 15.9% from a six-year study in Oxfordshire, UK^[35], and lower than that in some high-burden countries in Africa; for example, 66% in Karonga District, Malawi, calculated using strains from three consecutive years^[36]. The clustering rate in Kashgar was relatively high when compared to that in Han-dominated sites within China. For example, a one-year study in Guizhou (southwest China) calculated a clustering rate of 13.1%^[37], which is lower than that in Kashgar. A 10-year study in Shanghai (eastern China) obtained a clustering rate of 27.25%, which was only 1.35% higher than that reported by Kashgar^[38]. A three-year sample from Shenzhen (southern China) yielded a clustering rate of 12.2%^[39]. Therefore, we consider that Kashgar is among the regions with a more severe recent transmission within China and that recent transmission should be a local public concern. One possible reason for this is that Kashgar is predominantly Uyghur, whereas studies in North America have found that ethnic minorities are at a higher risk of recent transmission than the dominant ethnic group^[40,41]. Another possible reason is that the lack of transportation in Kashgar leads to a longer total delay (time interval between the onset of symptoms and the start of treatment) in patients, which exacerbates the transmission of TB; for example, the proportion of patients with a total delay of more than 8 weeks was 10.8% in Shenzhen when compared to 31.9% in Kashgar^[24]. In addition, since 2013, the prefecture has intensified its efforts to alleviate poverty and develop labor-intensive production workshops in rural areas to increase farmers’ incomes, providing more opportunities for people to congregate, thus facilitating the transmission of TB.

Based on the analysis of risk factors for clustering, we found that the DR strains were more likely to be transmitted than the susceptible strains in Kashgar. Since the transmission of strains is influenced not only by bacteriological factors, but also by host and social environmental factors, there is no uniform conclusion as to whether the risk of transmission of DR strains is higher than that of the susceptible strains. If only bacteriological factors are considered, DR mutations can result in fitness costs in the form of slower growth rates^[42]. Thus, the prevailing view was that DR strains were less transmissible; however, this could not explain why DR strains are widely prevalent. Later studies found that the contacts of patients with DR have a higher risk of infection^[43], which may have been caused by the higher number of MTB cells in the aerosol produced when patients with DR cough^[44,45]. Animal models show significantly higher bacterial loads in the airways after infection with MDR strains, suggesting that DR strains may have a transmission advantage over susceptible strains^[46]. These results suggest that bacterial-host interactions also influence the transmission of strains. In addition, the level of healthcare and capacity of the public health system are also drivers of strain transmission^[47]. Therefore, the role of socio-environmental factors cannot be disregarded. Perhaps, because of this, DR strains were less likely to spread than susceptible strains in areas where TB is effectively controlled^[48,49]. In high-burden countries such as China, DR strains were found to be at a higher risk of transmission than susceptible strains^[23,50]. Recent studies have suggested that among the various influencing factors, bacteriological factors play a relatively minor role in transmission, whereas socio-environmental factors play a greater role^[20]. Our results provide information for an in-depth analysis of the transmission patterns of DR strains and suggest that the social environment in Kashgar should be improved to suppress the spread of DR strains.

Although it has been suggested that smear-positive patients are more likely to transmit TB, a positive smear result was not a high-risk factor for recent transmission in our study. This may be due to the low sensitivity of smears, resulting in the high proportion of “false” smear-negative patients and thus masking the transmission advantage of smear-positive patients^[50]. The smear-positive rate for culture-positive patients in Kashgar was only 30.8%, when compared to 76% in a similar study in Malawi^[36].

Consistent with a previous study in China, Beijing strains favored transmission^[50]. A study in South Africa observed an expansion of Beijing strains from almost none to 33% in 40 years^[51]. These results indicate that the Beijing strains have a transmission advantage. Previous studies have shown that Beijing strains are more virulent^[52,53] and are better adapted to the host immune system, such as by evading BCG-induced immune protection^[54] or favoring a secondary Th1 response that is less protective of the host^[55]. These features may accelerate the progression from infection to disease^[50], allowing the strain to exhibit a transmission advantage. Although the proportion of lineage 3 in this region is relatively high in China, the transmission advantage of lineage 3 is not evident. The same conclusion was obtained in another area, with a ratio of lineage 3 close to that of Kashgar^[56]. We speculate that the high proportion of lineage 3 in Kashgar may be due to the presence of the ancient Silk Road^[19]. Moreover, the prefecture is adjacent to Pakistan and North India, where the frequency of lineage 3 strains is high^[57].

We found that students accounted for only 3.3% of all patients, but had a higher risk of cluster formation than other occupations. A previous study suggested that this was due to diagnostic delays^[24]. Students worry about interrupting their studies because of illness or a sense of shame. Another possible reason is that this population receives extra attention from the disease prevention and control departments. Following the diagnosis of TB in a student, their classmates are also screened^[23], leading to a more timely detection of students who develop TB via recent transmission.

This study has several limitations. First, the short duration of the study meant that some genomic clusters may have been missed, resulting in an underestimation of recent transmission levels, and the proportion of transmitted drug resistance. The novel coronavirus pandemic has added a heavier burden to local primary care institutions and created a major obstacle to our continuous collection of strains; therefore, we have only been able to collect a one-year sample. Secondly, the culture positivity rate in this region needs to be improved. The high laboratory equipment and technical requirements for MTB culture have led to its implementation being less than optimal in developing regions with a high TB burden like Kashgar^[58]. We cannot deny that culture-negative patients are also a part of the transmission chain. Finally, the lack of epidemiological survey data in this study, such as household contact tracing, makes it difficult to analyze the populations and sites at high risk of transmission in detail.

In conclusion, control of TB transmission should be strengthened in the prefecture to reduce the prevalence of drug resistance. Continuation of WGS-based TB surveillance and epidemiological surveys in the region will help provide greater insight into the causes of the high TB burden in Kashgar.