Prevalence and Molecular Characteristics Based on Whole Genome Sequencing of Mycobacterium tuberculosis Resistant to Four Anti-Tuberculosis Drugs from Southern Xinjiang, China.

OBJECTIVE: Drug-resistant tuberculosis is a major public health problem, especially in the southern region of Xinjiang, China; however, there is little information regarding drug resistance profiles and mechanism of Mycobacterium tuberculosis in this area. The aim of this study was to determine the prevalence and molecular characteristics of M. tuberculosis resistant to four anti-tuberculosis drugs from this area.
METHODS: Three hundred and forty-six isolates from the southern region of Xinjiang, China were included and used to perform phenotypic drug susceptibility testing and whole genome sequencing (WGS). Mutations in seven loci associated with drug resistance, including rpoB for rifampicin (RMP), katG, inhA promoter and oxyR-ahpC for isoniazid (INH), rrs 530 and 912 loops and rpsL for streptomycin (STR), and embB for ethambutol (EMB), were characterized.
RESULTS: Among 346 isolates, 106, 60, 70 and 29 were resistant to INH, RMP, STR and EMB, respectively; 132 were resistant to at least one of the four anti-tuberculosis drugs and 51 were multi-drug resistant (MDR). Beijing genotype and retreated patients showed a significantly increased risk for developing MDR tuberculosis. Compared with the phenotypic data, the sensitivity and specificity for WGS to predict resistance were 96.7% and 98.6% for RMP, 75.5% and 97.1% for INH, 68.6% and 99.6% for STR, 93.1% and 93.7% for EMB, respectively. The most common mutations conferring RMP, INH, STR and EMB resistance were Ser450Leu (51.7%) in rpoB, Ser315Thr (44.3%) in katG, Lys43Arg (35.7%) in rpsL and Met306Val (24.1%) in embB.
CONCLUSION: This study provides the first information on the prevalence and molecular characters of drug resistant M. tuberculosis in the southern region of Xinjiang, China, which will be helpful for choosing early detection methods for drug resistance (ig, molecular methods) and subsequently initiation of proper therapy of tuberculosis in this area.

Introduction

Tuberculosis, caused by Mycobacterium tuberculosis, remains one of the top 10 causes of death worldwide and the leading cause of death from a single infectious agent (ranking above HIV/AIDS).1 In 2020, World Health Organization (WHO) reported that 7.1 million people with tuberculosis were newly diagnosed and notified in 2019, up from 7.0 million in 2018 and a large increase from 6.4 million in 2017.2 China contributes 8.4% of the global incident tuberculosis cases, ranked third behind Indonesia (8.5%). A latest national tuberculosis epidemiology survey3 (performed in 2010) reported that the prevalence of tuberculosis in western China is significantly higher than that in central or eastern China, and Xinjiang province in the northwestern China is one of the most serious tuberculosis endemic areas of China. The evolution and spread of rifampicin-/multidrug-resistant tuberculosis (RR-/MDR-TB) poses a major obstacle for controlling tuberculosis disease in the world and China. It is estimated that close to half a million people developed RR-TB, of which 78% had MDR-TB globally.2 In a system review on studies from China, the prevalence of MDR-TB was 4.8% for new cases, 26.3% for retreatment cases in recent years,4 meanwhile isoniazid (INH) and rifampicin (RMP) resistance in retreated cases was found to be the most common with prevalence of 40.0% and 33.3%, respectively. Molecular drug susceptibility testing (DST) methods have advantages on closing the gap between detection and treatment of drug-resistant tuberculosis. Understanding the drug resistance mechanism of M. tuberculosis in a certain area is helpful to choose an appropriate DST. However, little is known on the mechanism of drug-resistant M. tuberculosis in southern Xinjiang, China, stressing the urgency to investigate the prevalence and molecular mechanism of drug resistance in this area.

Previous studies indicated that the most common molecular mechanisms in M. tuberculosis have been associated with mutations in the rpoB, katG, inhA promoter, rpsL, rrs and embB genes.5–9 Mutations in RMP resistant determined region (RRDR) of rpoB is the major cause for RMP resistance, with codons 450, 445 and 435 being the predominant sites.10 INH resistance is mainly related to the mutations in katG, followed by that in the inhA promoter and oxyR-ahpC intergenic region.11 Mutations in rpoB and katG were found attributed to 91.1–97.7% RMP resistance5,12 and 65–88.5% INH resistance,13,14 respectively. Resistance to streptomycin (STR) is mainly explained by nucleotide changes in rpsL and rrs genes particularly at codons rpsL43, rpsL 88 and in rrs 530 and 912 loops, accounting for 50–95% STR resistance.6,15,16 Mutations in embB were identified to be responsible for 38–73% ethambutol (EMB) resistant strains.15–17 Nevertheless, there is no recent data on the molecular nature of drug-resistant M. tuberculosis from southern region of Xinjiang, China.

Whole-genome sequencing (WGS) was recognized as a rapid and reliable method for determining drug resistance of M. tuberculosis and used as a routine investigation in a few high-income, low-tuberculosis burden countries such as England.18 WGS allows simultaneous identification of all known resistance-associated loci, resulting in the decrease of missed diagnosis of drug resistance of more drugs than other molecular methods, such as line probe assay and GeneXpert, which both include limited numbers of loci and are incapable to differentiate silent mutations from those affecting drug efficacy, leading to false negative and false-positive results, respectively.19,20

In this study, we analyze the resistant characterization and the sequence polymorphisms in seven chosen genes or regions associated with RMP, INH, STR and EMB resistance based on WGS data of 346 M. tuberculosis isolates from southern region of Xinjiang, China and evaluate the ability of WGS to predict resistance and susceptibility compared with phenotypic DST. The sequences or regions were chosen on the basis of their demonstrated association with drug resistance according to previous literatures.21–24 The results in the present study will be helpful to choose early detection methods of drug resistance (ig, molecular methods) and subsequent initiation of proper therapy of tuberculosis in southern region of Xinjiang, China.

Materials and Methods

Patients and Mycobacterial Isolates

This study was carried out between Sep 2017 and Dec 2018 in Xinjiang Uygur Autonomous Region Chest Hospital, Kashi Lung Hospital, Kuche County Infectious Disease Hospital and Wushi County Peoples Hospital, which all serve as the designated tuberculosis hospitals in southern region of Xinjiang. All of the pulmonary tuberculosis patients aged ≥16 years with positive cultures identified as M. tuberculosis complex (MTBC) and lived in southern Xinjiang, China for more than 10 years were interviewed and enrolled during the study period. Totally, 346 isolates were collected. Only one isolate per patient was collected and tested.

Drug Susceptibility Testing

The DST for all strains was performed in Xinjiang Uygur Autonomous Region Chest Hospital. Before DST, the strains were recovered on Löwenstein Jensen (L-J) medium for 4 weeks at 37°C. All positive cultures were tested the growth status in the 7H9 broth contained specific drug concentrations of 9 drugs, which were prepared in 96-well plates from Encode Medical Engineering Co., Ltd, Zhuhai, China. The critical concentrations for the four studied drugs indicate resistance are defined by the laboratory of Xinjiang Uygur Autonomous Region Chest Hospital and shown as following: INH, 0.4 µg/mL; RMP, 4.0 µg/mL; STR, 4.0 µg/mL; EMB, 5.0 µg/mL.

For the drug-susceptible isolates diagnosed by Xinjiang Uygur Autonomous Region Chest Hospital carrying gene mutations, the susceptibility of these isolates was repeated by proportional method using L-J slants in the National Institute for Communicable Disease Control and Prevention, Chinese Center for Disease Control and Prevention. The critical concentrations for each drug indicating resistance are shown as following: INH, 0.2 µg/mL; RMP, 40.0 µg/mL; STR, 4.0 µg/mL; EMB, 2.0 µg/mL.25

DNA Isolation

Genomic DNA was extracted from M. tuberculosis colonies on L-J medium using the cetyltrimethylammonium bromide (CTAB) method.26 Briefly, after the strains were killed at 80°C for 30 min, 100 μL of 10% sodium dodecyl sulfate (SDS) and 5 μL proteinase K (10 mg/mL) was added, and samples were vortexed for a few seconds and incubated for 10 min at 65°C. Then, the samples were added with 100 μL solution of CTAB-NaCl (4.1% NaCl and 10% CTAB) and then mixed with an equal volume of chloroform-isoamyl alcohol (24:1 [vol/vol]; 700 μL) and centrifuged for 15 min at 13,000 × g in an Eppendorf centrifuge. The aqueous phase (650 μL) was then separated and mixed with an equal volume of isopropanol. The samples were left at −20°C for 30 min and then centrifuged for 15 min at 13,000 × g. The DNA pellet was washed once with 70% ethanol, air dried, and resuspended in a final volume of 100 μL Tris-EDTA (TE, pH8.0).

Genome Sequencing

DNA libraries were prepared with genomic DNA using kits as instructed by the manufacturer. DNA libraries were then selected to perform cluster growth and 150 bp paired-end sequencing on DNB SEQ-2000 instrument (Beijing Genomics Institute, China). The raw FASTQ sequence reads were filtered by removing the adapter sequences, duplicate reads, low-quality reads that had a quality score below 20 in more than 30% of the bases. The clean reads were mapped to the genome of H37Rv (GenBank accession number, NC_000962.2) using in-house softwares Bowtie2 (Version 2.3.4.1) and samtools (Version 1.7). VarScan (Version 2.4.4) was used for single-nucleotide polymorphisms (SNPs) finding. All the genome wide SNPs were identified by the VarScan software by parsing the mapping genome sequence data, and then the SNPs related to phylogeny or located in PE/PPE gene family regions were filtered out. An average of 15.4 million sequence reads were acquired per genome at a depth of 500× and with coverage of 98.0%.

Mycobacterium Species Identification

A fast K-mer algorithm method (KmerFinder) was used to predict bacterial species in this study. After the Software KmerFinder 3.227 was downloaded, the local Strain Type library was used to identify the species. The Query coverage ≥85% was considered as MTBC species. All isolates in the present study were identified as MTBC.

Identify Mutations in Drug Resistance–Associated Genes or Regions

Identification of resistance-causing SNPs from genome-wide sequence is challenging. We chose 7 known resistance genes and regions on the basis of their demonstrated association with drug resistance and according to previous literatures21–24 (Table 1). All mutations in these genes and regions were compared with the pan-susceptible reference genome (H37Rv, accession number: NC_000962.2) at the level of SNPs in promoter regions or intergenic regions, amino acids in genes, or insertions and deletions. The phenotypic and genotypic results were compared to determine the specificity and sensitivity for each gene with WGS to predict resistance.

Table 1

Mutations Within Seven Drug Resistance-Associated Loci in Drug Resistant Mycobacterium tuberculosis Isolates

Drug (n*)	Gene/Region	Codon Change(s)	Amino Acid/Nucleotide Changes	No. of Mutated Isolates
RMP (60)	rpoBa	450(TCG-TTG)	Ser450Leu	27
		450(TCG-TTG)+150(GAC-GGC)	Ser450Leu+Asp150Gly	1
		450(TCG-TTG)+534(GTG-ATG)	Ser450Leu+Val534Met	1
		450(TCG-TTG)+545(GAC-GAG)	Ser450Leu+Asp545Glu	1
		450(TCG-TTG)+571(GAC-TAC)	Ser450Leu+Asp571Tyr	1
		450(TCG-ATG)	Ser450Met	4
		450(TCG-TGG)	Ser450Trp	1
		450(TCG-TGG)+432(CAA-GAA)	Ser450Trp+Gln432Glu	1
		430(CTG-CCG)	Leu430Pro	6
		430(CTG-CCG)+445(CAC-CAG)	Leu430Pro+His445Gln	1
		430(CTG-CCG)+491(ATC-CTC)	Leu430Pro+Ile491Leu	1
		445(CAC-TAC)	His445Tyr	3
		445(CAC-GAC)	His445Asp	1
		445(CAC-GAC)+435(GAC-GGC)	His445Asp+Asp435Gly	1
		445(CAC-GAC)+512(AAG-GAG)	His445Asp+Lys512Glu	1
		445(CAC-AAC)	His445Asn	2
		445(CAC-AAC)+435(GAC-GGC)	His445Asn+Asp435Gly	1
		445(CAC-GGC)	His445Gly	1
		435(GAC-TAC)	Asp435Tyr	1
		435(GAC-TTC)	Asp435Phe	1
		441(TCG-TTG)	Ser441Leu	1
INH (106)	katG	315(AGC-ACC)	Ser315Thr	38
		315(AGC-ACC)	Ser315Thr	2b
		315(AGC-ACC)+590(AAG-GAG)	Ser315Thr+Lys590Glu	4
		315(AGC-AAC)	Ser315Asn	2
		315(AGC-CGC)	Ser315Arg	1
		98(TAC-TCC)	Tyr98Ser$	1
		105(ATG-AAG)	Met105Lys$	1
		138(AAC-AGC)	Asn138Ser	1
		138(AAC-CAC)	Asn138His	1
		139(GCC-CCC)	Ala139Pro	1
		139(GCC-CCC)+140(AGC-AAC)	Ala139Pro+Ser140Asn	1
		145(CGC-CCC)	Arg145Pro$	1
		151(GTC-TTC)	Val151Phe$	1
		183(TTC-CTC)	Phe183Leu$	1
		189(GAC-GGC)	Asp189Gly	1
		191(TGG-GGG)	Trp191Gly	2
		234(GGG-GAG)	Gly234Glu$	1
		249(CGC-CAC)	Arg249His$	1c
		298(TTG-TCG)	Leu298Ser$	1
		382(CTC-CGC)	Leu382Arg$	1
		394(ACG-GCG)	Thr394Ala	1b
		481(TCG-TTG)	Ser481Leu$	1b
		498(CGC-CAC)	Arg498His$	1b
		619(CTC-CGC)	Leu619Arg$	1
		630(GGC-GTC)	Gly630Val$	1b
		632(CGC-CAC)	Arg632His$	1b
		711(TAT-GAT)	Tyr711Asp$	2
		728(TGG-CGG)	Trp728Arg$	1b
	inhA promoter	–	C-15T	6d
		–	T-8A	1e
		–	C-93T	1
	oxyR-ahpc intergenic region	–	C-52T	1e
		–	C-54T	4f
		–	C-57T	2e
		–	C-72T	1e
		–	C-73A	1
		–	T-77G	1e
STR (70)	rrs 530 and 912 loops	–	A514C	8
		–	C517T	1
		–	A908C	4
		–	T950C	1
	rpsL	43(AAG-AGG)	Lys43Arg	25
		88(AAG-AGG)	Lys88Arg	9
EMB (29)	embB	306(ATG-CTG)+983(CCG-CGG)	Met306Leu+Pro983Arg	1
		306(ATG-ATA)	Met306Ile	4
		306(ATG-GTG)	Met306Val	7
		306(ATG-ATC)	Met306Ile	3
		306(ATG-ATT)+655(CCG-GCG)+1024(GAC-AAC)	Met306Ile+Pro655Ala+Asp1024Asn	1
		406(GGC-GCC)	Gly406Ala	2
		406(GGC-GAC)	Gly406Asp	1
		406(GGC-AGC)	Gly406Ser	1
		246(GGC-CGC)	Gly246Arg	1
		328(GAT-GGT)	Asp328Gly	1
		328(GAT-GGT)+354(GAC-GCC)	Asp328Gly+Asp354Ala	1
		328(GAT-GGT)+1024(GAC-AAC)	Asp328Gly+Asp1024Asn	1
		354(GAC-GCC)	Asp354Ala	1
		495(GCC-ACC)	Ala495Thr	1
		575(ATG-ATA)	Met575Ile	1

Notes: aThe amino acid number is based on M. tuberculosis H37Rv codon number; bMeans that combined mutations in the inhA promoter or oxyR-ahpC intergenic region; cMeans that combined mutations in inhA promoter and oxyR-ahpC intergenic region; dMeans that 1 out of 6 isolates carried mutations in katG and 1 out of 6 isolates carried mutations in katG and oxyR-ahpC intergenic region; eCombined with mutations in katG; fMeans 1 out of 4 isolates carried mutations in katG and 1 out of 4 isolates carried mutations in katG and inhA promoter; *n means the No. of corresponding drug resistant isolates; $Novel mutations in katG.

Abbreviations: RMP, rifampicin; INH, isoniazid; STR, streptomycin; EMB, ethambutol.

Spoligotyping

All isolates were analyzed by spoligotyping, performed according to a standard protocol.28 Simply, the DNA was amplified with primers DRa (5ʹ- Biotin - GGTTTTGGGTCTGACGAC −3ʹ) and DRb (5ʹ- GGTTTTGGGTCTGACGAC −3ʹ). Then, the amplified products were hybridized with PALL Biodyne membrane prepared with 43 spacer oligonucleotide probes. After washing with 2×SSPE solution (360 mM NaCl, 20 mM NaH2PO4, 2 mM EDTA [pH 7.2]) supplemented with 0.5% SDS, the membrane was hybridized with streptavidin-peroxidase conjugate. The final image was detected with a chemiluminescence system, including the ECL detection liquid (Amersham, Buckinghamshire, United Kingdom) and ECL-Hyperfilm (Kodak, Rochester, NY).

Statistical Analysis

The spoligotyping results were entered into an Excel spreadsheet in binary format and compared with the spoligotyping database SpolDB4 (), the results were also analyzed using BioNumerics software (Version 5. 0, Applied Maths, Sint-Martens-Latem, Belgium). Cluster analysis was performed and a dendrogram was generated in Bionumerics using the Dice similarity coefficient and UPGMA coefficient. In addition, SPSS 16.0 (SPSS Inc., Chicago, IL, United States) was used to perform chi-square test and logistic regression analysis. A P value less than 0.05 was defined as significant. The extent of association was shown as an odds ratio (OR) and 95% confidence interval (95% CI).

Results

Demographic Information

A total of 346 patients diagnosed with pulmonary tuberculosis were enrolled with a mean age of 50.4 (±19.3) years. Of these, 48.3% (167/346) were males and 51.7% (179/346) were females; 53.5% (185/346) were new cases and 46.5% (161/346) were retreated cases.

Drug Resistance Profiles

According to the DST results acquired from Xinjiang Uygur Autonomous Region Chest Hospital, a total of 78 drug-susceptible isolates were found carried mutations in the 7 genes or regions which maybe attributed to the higher critical drug concentrations defined by Xinjiang Uygur Autonomous Region Chest Hospital, so DSTs on these isolates were repeated in National Institute for Communicable Disease Control and Prevention, Chinese Center for Disease Control and Prevention and then the susceptibility patterns of these isolates was adopted in the present study.

Phenotypic DST results against the four anti-tuberculosis drugs showed that among 346 isolates, 214 were fully drug susceptible, 132 (38.2%; 95% CI, 33.0–43.3%) were resistant to at least one drug; the isolates resistant to INH, RMP, STR and EMB were 106 (30.6%; 95% CI, 25.8–35.5%), 60 (17.3%; 95% CI, 13.3–21.4%), 70 (20.2%; 95% CI, 16.0–24.5%) and 29 (8.4%; 95% CI, 5.4–11.3%), respectively. In total, 51 isolates (14.7%; 95% CI, 11.0–18.5%) were identified as MDR. Thirty-eight out of 51 (74.5%) MDR isolates were from retreated patients. A full susceptibility profile for all strains is shown in Table 2.

Table 2

Drug Susceptibility Patterns of 346 Clinical Mycobacterium tuberculosis Isolates

Susceptibility or Resistance	Number of Strains
Fully susceptible*	214
H-monoresistant	37
R-monoresistant	4
S-monoresistant	16
E-monoresistant	1
Over all poly-resistant	23
HS	17
RS	3
HSE	1
RSE	2
MDR	51
HR	15
HRS	11
HRE	5
HRSE	20

Note: *Means the isolates were simultaneously susceptible to isoniazid, rifampicin, streptomycin, and ethambutol.

Abbreviations: H, isoniazid; R, rifampicin; S, streptomycin; E, ethambutol; MDR, multi-drug resistant.

Genotype Distribution of the M. tuberculosis Isolates

Among the M. tuberculosis isolates for WGS, 202 (58.4%) belonged to the Beijing genotype, while 144 (41.6%) were non-Beijing family, which included the CAS1-Delhi family (47), Ural-2 family (31), T1 family (12), EAI family (11), T2 family (5), T3 family (3), T family (2), H3 family (2), LAM9 family (2), H1 family (1), Ural-1 family (1), and un-defined genotypes (27).

Among the 346 isolates, a total of 53 spoligotypes were identified. Of these, 33 spoligotypes were previously represented as Shared International Types (SITs) according to SpolDB4.0, while the other 34 were reported for the first time (Table 3). After clustering with BioNumerics software, 312 (90.2%) isolates were classified into 19 clusters containing 2 or more strains. Additionally, 34 (9.8%) strains did not form clusters (Figure 1, Table 3). Of the 19 clusters, about 76.3% (264/346) of all isolates were contained in 6 predominant clusters, including SIT 1 (Beijing, 181 isolates), SIT 127 (Ural-2, 30 isolates), SIT 25 (CAS1-Delhi, 16 isolates), SIT 26 (CAS1-Delhi, 13 isolates), SIT 357 (CAS1-Delhi, 13 isolates) and SIT 27 (EAI, 11 isolates).

Table 3

Numbers and Frequencies of Mycobacterium tuberculosis Isolates Clustered by Spoligotyping

Parameter	Value
No. of isolates studied	346
No. of clusters	19
No. of new found spoligotypes	30
Mean no. of isolates per cluster	16.2
No.(%) of clustered isolates	312 (90.2%)
No.(%) of unclustered isolates	34 (9.8%)

Figure 1

Minimum spanning tree generated with spoligotypes of Mycobacterium tuberculosis.

Factors Linked to Drug-Resistant Tuberculosis

Risk factors for MDR and drug-resistant (but not MDR) tuberculosis were also analyzed in the present study. As shown in Table 4, the Beijing genotype infected with and retreated patients presented significantly higher risks for developing MDR-TB than pan-susceptible (simultaneously susceptible to INH, RMP, STR and EMB) tuberculosis with ORs of 3.05 (95% CI, 1.4 to 6.6; P< 0.01) and 6.46 (95% CI, 3.1 to 13.6; P< 0.01), respectively.

Table 4

Factors Associated with Drug-Resistant Tuberculosis

	No.(%) of Isolates			Drug-Resistant but Not MDR TB vs Pan-Susceptible& TB		MDR TB vs Pan-Susceptible& TB
Factor	Susceptible TB	Drug-Resistant but Not MDR TB	MDR- TB	Odds Ratio (95% CI)	P value	Odds Ratio (95% CI)	P value
Sex
Male	102	39	26	Reference		Reference
Female	112	42	25	1.05(0.63–1.77)	0.85	0.85(0.43–1.68)	0.64
Age group
<30 yrs	47	14	15	Reference		Reference
30–59 yrs	77	30	21	0.62(0.28–1.35)	0.23	2.44(0.92–6.41)	0.07
≧60 yrs	90	37	15	0.94(0.53–1.67)	0.83	1.95(0.87–4.39)	0.11
Occupation
Others	75	48	22	Reference		Reference
Farmer	139	33	29	0.79(0.41–1.54)	0.49	0.75(0.34–1.71)	0.49
Residence area
Rural	162	58	41	Reference		Reference
Urban	52	23	10	1.11(0.56–2.20)	0.76	0.48(0.19–1.22)	0.12
Treatment
New cases	137	35	13	Reference		Reference
Retreated cases	77	46	38	1.23(0.72–2.10)	0.45	6.46(3.07–13.60)	<0.01*
Genotypes
Non-Beijing	100	33	11	Reference		Reference
Beijing	114	48	40	1.32(0.77–2.24)	0.31	3.05(1.40–6.64)	<0.01*

Notes: *P<0.05 (significant); &Pan-susceptible means that the Mycobacterium tuberculosis isolates were simultaneously susceptible to isoniazid, rifampicin, streptomycin and ethambutol.

Abbreviations: CI, confidence interval; TB, tuberculosis; MDR, multi-drug resistant.

RMP and rpoB

RMP resistance-conferring mutations in rpoB in M. tuberculosis isolates were identified by WGS. A total of 21 genotype patterns in rpoB were identified. A 96.7% (58/60) RMP-resistant M. tuberculosis isolates carried mutations in the RRDR of rpoB. The most frequently mutated codons were 450, 445, 430 and 435 with mutation frequencies of 61.7% (37/60 isolates), 18.3% (11/60 isolates), 13.3% (8/60 isolates), and 6.7% (4/60 isolates) (Table 1). Two novel mutations Asp150Gly and Asp571Tyr were only found in RMP-resistant isolates and both combined with mutations in the rpoB RRDR region (Table 1). Four out of 286 (1.4%) RMP-susceptible isolates carried nonsynonymous mutations within the whole sequence of the rpoB gene (), of which 1 isolate carried Gln432Glu, which is in the RRDR of rpoB. The detection of mutations in rpoB by WGS resulted in 96.7% sensitivity and 98.6% specificity compared with phenotypic DST (Table 5).

Table 5

The Ability of Whole-Genome Sequencing Analysis for Drug Resistance Prediction in Comparison with the Phenotypic Drug Susceptibility Testing

Drugs	Genes	No. of Resistant Isolates (%)		No. of Susceptible Isolates (%)		χ2	P	Sensitivity (%)	Specificity (%)
		With Mutations	Without Mutations	With Mutations	Without Mutations
RMP	rpoB	58	2	4	282	306.0	<0.01	96.7	98.6
INH	katG	72	34	5	235	184.2	<0.01	67.9	97.9
	inhA promoter	8	98	1	239	12.1	<0.01	7.5	99.6
	oxyR-ahpC intergenic region	10	96	1	239	16.6	<0.01	9.4	99.6
	katG and inhA promoter	77	29	6	234	198.4	<0.01	72.6	97.5
	katG, inhA promoter and oxyR-ahpC intergenic region	80	26	7	233	205.6	<0.01	75.5	97.1
STR	rrs 530 and 912 loops	14	56	1	275	47.3	<0.01	20.0	99.6
	rpsL	34	36	0	276	148.7	<0.01	48.6	100.0
	rrs 530 and 912 loops or rpsL	48	22	1	275	208.4	<0.01	68.6	99.6
EMB	embB	27	2	19	298	166.4	<0.01	93.1	94.0

Abbreviations: RMP, rifampicin; INH, isoniazid; STR, streptomycin; EMB, ethambutol.

INH and katG, inhA Promoter and oxyR-Ahpc Intergenic Region

According to the WGS results, among 106 INH-resistant and 240 INH-susceptible isolates, 72 and 5 carried mutations in katG, respectively, whilst 8 and 1 in the inhA promoter and 10 and 1 in the oxyR-ahpC intergenic region, respectively (Table 1 and ). The katG Ser315Thr was the most dominant mutation found in 44.3% (47/106) phenotypic INH resistant strains; among the isolates carried this mutation, only two combined mutations in the inhA promoter and/or the oxyR-ahpC intergenic region. Among 25 INH-resistant isolates, which carried mutations in katG non-315, seven combined mutations in the inhA promoter and/or the oxyR-ahpC intergenic region. Among 34 INH-resistant isolates, which carried the wild-type of katG, eight carried mutations in the inhA promoter or the oxyR-ahpC intergenic region. Five isolates only carried mutation in the inhA promoter, whilst three only carried mutations in the oxyR-ahpC intergenic region (Table 1 and ). There were still 26 INH-resistant strains with the wild-type of three sequenced gene and regions.

Of 24 additional mutations except katG315, 23 were found only in phenotypic INH-resistant isolates (Table 1, ). Eight of the 23 combined mutations in katG, inhA promoter and/or oxyR-ahpC intergenic region () and 16 of the 23 were novel mutations in katG (Table 1). A total of 7 INH susceptible isolates were also found mutations in katG, inhA promoter and/or oxyR-ahpC intergenic region ().

Compared with the phenotypic DST, WGS predicts INH resistance based on mutations in katG combined with inhA promoter and oxyR-ahpC intergenic region with higher sensitivity of 75.5% (VS 67.9%, 9.4% and 7.5%) than that based on the mutations from a single gene, however with a lower specificity of 97.1% (VS 99.6% and 97.9%) (Table 5).

STR and rrs 530 and 912 Loops and rpsL

The WGS analysis on mutations in rrs 530 and 912 loops and rpsL associated with STR resistance showed that, among 70 STR-resistant isolates, 14 carried mutations in rrs 530 and 912 loops, 34 in rpsL. As shown in Table 1, for rrs, the most frequent mutations were A514C and A908C, which occurred in 11.4% (8/70) and 5.7% (4/70) phenotypic STR-resistant isolates. For rpsL, the most frequent codons were 43 and 88, which were observed in 35.7% (25/70) and 12.9% (9/70) phenotypic STR-resistant isolates. Mutations in rrs and rpsL genes accounted for 68.6% of the STR-resistant isolates. We also found that, among 276 STR-susceptible isolates, only one carried mutation C517T in rrs (). The sensitivity and specificity of WGS to predict STR resistance according to mutations in rpsL combined with rrs were 68.6% and 99.6% (Table 5).

EMB and embB

The whole embB sequence of 346 M. tuberculosis isolates was analyzed in this study. Among the 29 isolates resistant to EMB, 27 isolates (93.1%) carried mutations in embB. The most prevalent codon was Met306 (16, 55.2%), where the codon ATG (Met) was replaced with GTG (Val, 7, 24.1%), ATA (Ile, 4, 13.8%), ATC (Ile, 3, 10.3%), ATT (Ile, 1, 3.4%) and CTG (Leu, 1, 3.4%), respectively. Four of EMB-resistant isolates carried Gly406 codon mutations, where the codon GGC (Gly) was replaced with GCC (Ala, 2, 6.9%), GAC (Asp, 1, 3.4%) and AGC (Ser, 1, 3.4%) (Table 1). Other mutations at codon 328 (n =3), codon 354 (n = 2, one was combined with condon 328), codon 246 (n = 1), codon 495 (n = 1) and codon 575 (n = 1) were found in EMB resistant isolates, as listed in Table 1. Among the EMB-susceptible isolates, mutations in Met306 and Gly406 were also the most predominant in embB gene, furthermore, 12 mutations except Met306 and Gly406 were found (shown in ). WGS analysis on embB had a sensitivity and specificity of 93.1% and 94.0%, respectively (Table 5).

Discussion

The prevalence of RR-TB and MDR-TB patients was 17.3% (60/346) and 14.7% (51/346) in the southern Xinjiang, China, in the present study, both comparable to a previous study from Xinjiang (19.4% and 13.2%),29 however, both higher than the data from a national survey of drug-resistant tuberculosis in 2007 in China (5.1% and 10.2%),30 indicating a serious epidemic of drug-resistant tuberculosis in the southern Xinjiang, China. Yet the epidemic of drug-resistant tuberculosis in southern Xinjiang, China maybe overestimated due to that the isolates of this study were collected from tuberculosis designated hospitals rather than random survey.

The Beijing genotype family is the dominant lineage in southern Xinjiang, accounting for 58.4% of strains in the present study, which is in line with Yuan et al’s report from Xinjiang, China31 and the data from the south of China (53.2%),32 but lower than that from the north of China (76.5%).32 According to subtype analysis with SITs, the most predominant clusters were SIT 1 (Beijing), SIT 127 (Ural-2), SIT 25 (CAS1-Delhi), SIT 26 (CAS1-Delhi), SIT 357 (CAS1-Delhi) and SIT 27 (EAI). Three subtypes SIT 357, 25 and 26 of CAS-Delhi family were evenly distributed in this area (all about 3%, Figure 1). A previous study from Xinjiang province has demonstrated that Beijing genotype strains might be correlated with INH and EMB resistance.33 Another epidemiological study reported that the Beijing genotype showed greater correlation with RMP and ofloxacin resistance and MDR phenotypes in China.34 We also observed that the Beijing genotype exhibited a significantly higher risk for developing MDR-TB compared to non-Beijing genotype, suggesting that Beijing genotype be responsible for the spread and emergence of MDR-TB in this region. As reported in previous studies,35,36 retreated patients in this study were also found with high risk for developing MDR-TB, stressing the importance of early diagnosis, and timely, long-term and standardized medication, especially among patients treated within the hospital system in this region.

Our results present the first WGS-based molecular characterization of M. tuberculosis from southern region of Xinjiang in China. We focused on genes known to confer resistance to four anti-tuberculosis drugs in M. tuberculosis. WGS could offer a rapid and comprehensive SNP identifications, which help to understand the drug resistance mechanisms and choose suitable molecular DSTs suitable for this region leading to quicker and more appropriate treatment.

Previous reports show that mutations in rpoB are the main cause of RMP resistance.14,37 Our results showed a sensitivity of 96.7% based on rpoB mutations for predicting RMP resistance by WGS, consistent with the data from Jiangxi province, China14 and France.37 The most common mutations are in the rpoB RRDR region, particularly in codons 450, 445, 430 and 435,14,22,37 which are in line with the results that 86.7% of RMP resistant isolates had at least one of these four codon mutations in the present study. Two novel mutations Asp150Gly and Asp571Tyr combined with mutations in the rpoB RRDR and only existed in RMP-resistant isolates. It is unclear if these mutations play a direct role in RMP resistance; therefore, their possible functions should be verified in further studies.

Mutations in the katG gene are often considered to be the major contributors for INH resistance, followed by that in the inhA promoter and oxyR-ahpC intergenic region.12 Accordingly, our study indicated that mutations in katG accounted for 67.9% INH resistance, lower than results from Hunan province, China,36 America38 and Kyrgyz Republic.39 The differences can be attributed to geographical variations. In the present study, 23 additional mutations except katG315 were found in INH resistant isolates, of which 16 were novel mutations, and only eight novel mutations combined with katG315 or inhA promoter or oxyR-ahpC intergenic region mutations. All of these novel mutations were found only in phenotypic INH-resistant isolates, suggesting that these mutations were resistance-associated but needed to be further verified by in vitro mutagenesis experiments. In addition, 8 out of 34 INH-resistant isolates possessed wild-type katG carried inhA C (−15) T, T (−8) A and C-93T and/or mutations in the region of oxyR-ahpC from −77 to −52. So, the inhA (−15), T (−8) A and oxyR-ahpC −77 to −52 combined with katG315 can make a preferable set for INH-resistance diagnoses. The mutation of inhA C (−93) T is first reported by us, its role on INH resistance needs to be further confirmed. The combination of mutations in katG, inhA, and the oxyR-ahpC intergenic region were found in 75.5% of INH-resistant isolates, which were far lower than 92.7% reported in another study from China.6 There may be alternative mechanisms involved in the INH resistance of these isolates, and possibly the mutations occur in other structural genes or other loci. A review about mechanisms of isoniazid resistance suggested that mutations in kasA and ndh genes were the tertiary cause for INH resistance, followed by iniABC, fadE and furA.11 Besides, several recent studies pointed out the contribution of efflux pump to INH resistance, and this could explain why the remaining 20–30% of phenotypically resistant isolates that do not contain any genotypic mutation.40,41 Seven INH susceptible isolates carried mutations in katG, inhA promoter or oxyR-ahpC intergenic region, which may be conferred with low-level INH resistance and easily lead to misdiagnosis or not be associated with INH resistance.

Resistance to STR is due to alterations within rpsL, rrs (530 loop and the 912 loops).42 Two mutations rpsL Lys43Arg and Lys88Arg were the most common in STR resistant isolates in the present study, which were similar to the data from other areas in China17 and Pakistan.23 The most frequent mutation in rrs in STR resistant isolates was A514C, followed by C517T and A908C, in line with other reports.6,42 There was no isolate carried rrs mutation combined with rpsL mutation in this study, suggesting that mutations in rrs and rpsL genes were mutually exclusive. Mutations in rrs and rpsL genes accounted for 68.6% STR resistance, higher than the data from the tertiary care tuberculosis hospital in China35 and Mexico,43 but lower than that from other regions,44,45 indicating that there were regional differences in the mutations associated with STR resistance and additional resistance mechanisms correlated with STR resistance, such as gidB which was reported to contribute 2.4% to 37.5% STR resistance.42,44,46–48

EMB resistance was mostly associated with mutations within the embB gene (codon 306 as the dominant).50 In the present study, embB mutations were found in 93.1% EMB resistant isolates, showed a high mutation prevalence compared to previous reports, ranged from 20% to 90.9%.21,49,50 The most prevalent mutation at locus 306 was detected, which was also the most predominant codon in EMB-susceptible isolates. The phenomenon was also reported in previous studies.8,51 Overall 6.0% of EMB-susceptible isolates also carried mutations in embB, which is close to the data from a previous study (6.5%),52 which maybe attributed to the narrow range of EMB critical concentrations for differentiating resistant and susceptible strains, and the possible presence of microcolonies that are difficult to detect visually.53 A previous study shows that, using the EMB concentration with 1.6 µg/mL instead of 2.0 µg/mL in L-J slants by the proportional method, more than 90% EMB-susceptible isolates that carried embB306 mutations could be successfully recognized as EMB-resistant isolates, while the EMB susceptible isolates with wild-type embB were not changed.54

Studies on EMB resistance showed that mutations outside embB codon 306 and 406, do occur but are quite rare.51,55 However, substitutions at other codons were identified in 12 EMB-susceptible and eight EMB-resistant isolates. Two novel mutations embB Ala495Thr and Met575Ile were found only in EMB resistant isolates. The role of these mutant types remained unknown, therefore required further exploration.

Two out of 29 EMB resistant isolates were not found mutations in embB. Mutations in embC and embA, which were shown to be involved in EMB resistance development,15,49 were not included in the present study. However, Wan et al6 found that mutations in embC or embA always combined with mutations in embB, and both EMB resistant and susceptible isolates had the same low frequencies of mutations in embC or embA.

A systematic review about WGS in M. tuberculosis for detection of drug resistance pointed out that the sensitivities of WGS compared with phenotypic DST for the following anti-tuberculosis drugs were high but also varied: RMP (89.2% to 100%), INH (90% to 100%), EMB (71.4% to 95.8%), STR (57.1% to 89.9%).56 In the present study, the sensitivities between the WGS and phenotypic DST for the four anti-tuberculosis drugs were similar to previous studies, ranging from 68.6% to 96.7%.

Conclusions

Comprehensively identifying mutations by WGS could provide a comprehensive understanding on the drug resistance mechanism caused by chromosome mutations. The data in the present study raises our understanding of the prevalence and molecular characteristics of M. tuberculosis resistant to INH, RMP, EMB and STR in southern region of Xinjiang China, which would be benefited for choosing early detection methods of drug resistance and subsequently initiation of proper therapy of tuberculosis in southern region of Xinjiang, China.