Whole Genomic Analysis Revealed High Genetic Diversity and Drug-Resistant Characteristics of Mycobacterium tuberculosis in Guangxi, China

Introduction

Tuberculosis (TB), caused by etiological Mycobacterium tuberculosis (Mtb), remains a public health worldwide and the leading cause of death from a single infectious agent until the COVID-19 pandemic.¹ Previous studies revealed that M. tuberculosis originated in Africa and spread with human migration worldwide, it developed to a contemporary phylogeny consisting of nine recognized lineages (L1–L9) with vary in their distribution between countries and continents.² In China, Lineage 2 (L2) and Lineage 4 (L4) are the most widely distributed throughout the country, whereas Lineage 1 (L1) and Lineage 3 (L3) are most prevalent in Southeast and Northwest, respectively.³

To date, the emergence of drug-resistant tuberculosis poses a great challenge for TB control and prevention. The WHO Global Tuberculosis Report of 2022 estimated that the burden of drug-resistant TB has increased between 2020 and 2021, with 450000 (95% UI: 399 000–501 000) new cases of rifampicin resistant in 2021.⁴ The treatment success rate for rifampicin-resistant TB is typically in the range of 50–75%.¹ It is urgent to learn about the drug-resistant profile of the M. tuberculosis isolates in order to reduce the risk of infection and transmission. Whole-genome sequencing (WGS) is a useful tool for predicting the susceptibility to anti-TB drugs in clinical specimen.⁵ In some countries, WGS of M. tuberculosis is increasingly used in routine care settings for species identification, drug resistance profile determination and transmission investigation.^6–8 Traditionally, a great number of single nucleotide polymorphisms (SNPs) mutation in some genes have been identified as associated with drug resistance.⁹

Although TB control has made noticeable progress over the past years, China remains a high burden country worldwide.¹ Meanwhile, the prevalent genotype and drug resistance of the Mtb isolates show a great difference among diverse regions. Guangxi, located on the border of China and borders Vietnam, is one of the high TB-burden regions in South China. Previous studies have reported the epidemiology of Mtb in hot and cold spots of Guangxi,^10,11 but the detailed lineage distribution evolution and drug-resistant profile of Mtb are limited. In this study, we analyzed the genome of M. tuberculosis isolates and the clinical system from pulmonary TB patients in Guangxi, which can provide insights into the genetic diversity and drug-resistant profile of M. tuberculosis isolates.

Methods

Study Setting and Sample Enrollment

In this retrospective cross-sectional investigation, five locations (Guigang, Fangchenggang, Baise, Chongzuo, and Guilin) in Guangxi were selected as study sites. Geographically, these locations are dispersed over the east, south, west, and north of Guangxi. To be specific, Baise, Chongzuo, and Fangchenggang bordering Vietnam represent the border monitoring points, and Guilin and Guigang belong to non-border monitoring points (Figure 1).

Figure 1 Distribution of the study site in Guangxi, China. The five cities located in border and non-border of Guangxi, which belongs to the southwest region of China.

During February 2016 and August 2017, we collected 1638 Mtb strains isolated from pulmonary tuberculosis patients with sputum smears positive, who attended local designated hospitals. Then, 600 isolates were randomly selected for drug sensitivity testing and whole genome sequencing. Due to issues like unsuccessful re-culture, contamination, unsuccessful genome extraction, and failure to match the epidemiologic information, a total of 462 clinical isolates were ultimately included in the in-depth study and further analysis (Figure S1 and Table S1). The epidemiologic investigations were collected from the patient after written informed consent. The study received ethical approval, which complies with the Declaration of Helsinki, from the Ethics Committee of the Guangxi Center for Disease Control and Prevention (IRB00001594).

Drug-Susceptibility Testing

Drug-susceptibility testing (DST) was performed by conventional proportion methods as the described previously,¹² with the following drug concentration: isoniazid (0.2 μg/mL), rifampin (40 μg/mL), ethambutol (2 μg/mL), streptomycin (4 μg/mL), kanamycin (30 μg/mL), ofloxacin (2 μg/mL). The isolate was regarded as resistant if the proportion of growth on the medium containing drug was 1% of that observed in drug-free medium. The reference isolate H37Rv (ATCC 27294) as a quality control.

DNA Extraction and Sequencing

All the M. tuberculosis isolates were cultured on Lowenstein Jensen media and genome DNA was extracted following the protocol of the cetyltrimethylammonium bromide (CTAB) method.¹³ The whole genome sequence was performed on purified DNA using Illumina HiSeq 2500 platform. All the whole genome sequencing operations were completed by Annoroad Gene Technology company (Beijing, China).

Bioinformatic Analysis

The WGS data analysis was carried out using procedures that were similar to those reported previously. The paired-end reads were examined using FastQC (v0.11.9) and trimmed using Trimmomatic (v0.39).¹⁴ The qualified reads were mapped against the H37Rv genome (NC_000962.3) by BWA-MEM algorithm (v0.7.16). Samtools (v1.13) and bcftools (v1.13) were used to variant calling from the sorted mapped sequences. The variant filtration was performed based on the following criteria: mapping quality greater than 30 and minimum read depth of 10. SnpEff (v4.3) was used to annotate the filter variant SNPs. And the isolate mutation conferring resistance to six anti-tuberculosis drugs was identified in TBProfiler.¹⁵

Phylogenetic Analysis

Phylogenetic analysis was performed by snippy pipeline (v4.3.6) for alignment of the core SNPs (https://github.com/tseemann/snippy). SNPs’ positions within PE/PPE genes or other repetitive regions with low mapping scores were removed. Then, the recombinant region of the alignment was filtered with Gubbins. The maximum-likelihood phylogenetic tree based on the core SNPs was construct by RAxML, with 1000 bootstrap iterations and the GTR+G model. The lineage of the M. tuberculosis isolates was identified using fast-lineage-caller v1.0 (https://github.com/farhat-lab/fast-lineage-caller). And we used Coll et al and Shitikov et al SNP schemes to define the sub-lineages and genotype of lineage 2, respectively.^16,17 The genomic cluster was defined as strains with a genetic distance of 12 SNPs or less, suggesting they were the result of recent transmission.¹⁸ The genomic pairwise distance was calculated using the program snp-dists (https://github.com/tseemann/snp-dists).

Statistical Analysis

All statistical analysis was used with SPSS (v18.0). The Chi-square test or Fisher’s exact test was performed for this study data. P < 0.05 was considered statistically significant.

Results

Demographic Characteristics of the Patients

Among the enrolled patients, 78.6% (363/462) were male, and the age of 41.8% of patients belong to 45–65 years. The majority of patients are living in rural (83.5%) and being Han (68.0%). And 34.9% and 40.5% of the patients have primary education and junior education, respectively. In addition, most patients do not have diabetes (77.7%) or hepatitis (82.7%). We also investigated the clinical symptoms of the patients. Most patients have the symptoms of cough (98.5%) and expectoration (92.4%), and various percentages of patients also have symptoms like hemoptysis (11.7%), chest pain (8%), fever (9.9%), weak (17.5%), night sweat (10.2%), and other conditions (6.7%). Other detailed information is shown in Table 1.

Table 1 The Demographic Information and Clinical Characteristics of the Mtb Infection Patients

Population Structure of M. tuberculosis

The population structure of Mtb isolates in Guangxi was investigated (Figure 2). The maximum likelihood phylogenetic tree and the barcode-SNP–based genotyping show that majority of isolates (58.8%, 272/462) belong to lineage 2.2.1 (L2.2.1), and 2.2% (10/462) isolates to lineage 2.2.2 (L2.2.2), which are the subgroup of Beijing lineage (L2). Unexpectedly, 6.5% (30/462) belong to the proto-Beijing sub-lineage (L2.1), which is not found in most areas of China and other regions of the world.^3,19 Within lineage 4 (L4), 48.2% (67/139) isolates belong to lineage 4.4, followed by lineage 4.5 (28.1%, 39/139) and lineage 4.2 (23.0%, 32/139), only one isolate is lineage 4.1. In addition, 2.2% (10/462) isolates belong to lineage 1 (L1) and one isolate belongs to lineage 3 (L3) in this study.

Figure 2 Phylogenetic tree of the 462 Mtb isolates in this study. The maximum likelihood phylogenetic tree was constructed by Raxml, with 1000 bootstraps replicates. The sub-lineages, drug resistant type and drug resistant profile of the isolates are shown on the tree (from inner to outer circles), according to the color legend shown on the right. The red stars indicate the cluster isolates.

Meanwhile, we further classified the isolates of lineage 2 according to Shitikov’s schemes,¹⁷ which divided the L2 isolates into 11 genotypes based on genetic markers. The results showed that 69.2% (216/312) of L2 isolates were divided into seven genotypes in our study (Figure 2). Specifically, 21.2% (66/312) of L2 is Asia Ancestral 3, followed by Asian African 2 (13.1%, 41/312), Asian African 3 (12.2%, 38/312), Pacific RD150 (9.9%, 31/312) and Proto-Beijing (9.3%, 29/312). And small amounts of isolates belong to Asia Ancestral 2 (2.9%, 9/312) and Asian African 1 (0.6%, 2/312). However, there are some L2 isolates that do not belong to these genotypes and more detailed classification may be needed to identify them. Overall, our results indicate that Mtb isolates in Guangxi with high genetic diversity, and majority of TB cases were caused by L2 and L4 isolates.

The Transmission Dynamic Analysis

The cluster analysis employing a threshold of a maximum distance of 12 SNPs grouped 58 (12.6%) of the 462 isolates into 28 clusters, ranging in size from two to three isolates (Figure 3). The genetic differences between isolates in each cluster ranged from 1 to 11 SNPs; 16 pairs differed by less than 8 SNPs and 12 pairs by 10 to 11 SNPs. The isolates from different cities are interspersed in the phylogenetic tree, but the clustering rate in border areas is higher than that in non-border areas (P < 0.01). Among all the cluster isolates, 25 pairs (89.3%) isolate from the same city and only 3 pairs (10.7%) isolate from separate cities (Figure 3). We also noted that there was a statistical difference in the clustering rate between the Beijing genotype and non-Beijing genotype (P < 0.05).

Figure 3 Phylogenetic tree of the 58 Mtb isolates with 28 clusters in this study. The lineages, cities and genotypic drug resistance of isolates are shown on the tree according to the color legend.

Phenotypic Drug-Resistant Profile

Among all the isolates in our study, 21.1% (98/462) were resistant to isoniazid, 13.2 (61/462) to rifampicin, 6.7% (31/462), 12.1% (56/462), 6.7% (31/462) and 1.9% (9/462) isolates were resistant to ethambutol, streptomycin, ofloxacin and kanamycin, respectively (Figure 4A). According to the latest definition of drug-resistant TB,²⁰ 1.3% (6/462) isolates are RR-TB (rifampicin-resistant tuberculosis), 7.1% (33/462) isolates are MDR-TB (multidrug-resistant tuberculosis) and 3.0% (14/462) are Pre-XDR (pre-extensively drug-resistant), whereas 7.1% (33/462) and 11.0% (51/462) isolates belong to MR-TB (mono-resistant tuberculosis) and Hr-TB (isoniazid-resistant tuberculosis), respectively (Figure 4B). These results showed that high phenotypic resistance proportions of the Mtb isolates in Guangxi, especially for the isoniazid resistance.

Figure 4 The drug-resistant profiles of the isolates in this study. (A and B) The phenotypic drug resistance and different phenotypic drug type of the isolates in this study. (C and D) Comparison of drug-resistant profiles between Beijing phenotype and non-Beijing phenotype. (E and F) Comparison of drug-resistant profiles between Border monitoring points and non-border monitoring points.

Unlike the previous report showed that high prevalence of XDR-TB in proto-Beijing (L2.1),¹⁹ we found that only 13.3% (4/30) isolates belong to Hr-TB and one isolate is MR-TB in L2.1 of this study (Figure 2). By analyzing, the prevalence of drug resistance of non-Beijing genotype was higher than that of Beijing genotype strains, especially for isoniazid resistance (Figure 4C), but there was no significant difference between these two genotypes, including the drug types (Figure 4D). In addition, our data showed that the isolates of border monitoring points are more likely to be resistant to rifampicin, isoniazid, streptomycin and ofloxacin resistant than those of no-border monitoring points (P < 0.01) (Figure 4E). RR-TB and MDR-TB isolates are more likely to be present in border monitoring points (p < 0.05) (Figure 4F).

Drug-Resistance Related Gene Mutations

The drug resistance profiles based on the whole genome sequence were assessed. Of the 61 phenotypic rifampicin resistant isolates, the most prevalent mutation was Ser450 Leu (37.7%, 23/61), followed by His445Asn (9.8%, 6/61) and His445Tyr (6.5%, 4/61). Combination mutations were detected in 7 (11.5%, 7/61) isolates, and two of these isolates also carried mutations in the rpoC gene (Table S2). Eighty-five (18.3%, 85/462) isolates detected mutations related to isoniazid resistance genes including katG, ahpC, fabG1and inhA (Table S3). Of these mutations, KatG Ser315Thr (62.2%, 61/98) was the most common mutation, followed by fabG1 −15C>T (10.2%, 10/98). And eight isolates (8.2%, 8/98) carried combined mutations. Thirty-one isolates predicted ethambutol resistance were mainly related to the mutation embB Met306Val (45.2%, 14/31) and Met306 Ile (22.6%, 7/31) (Table S4). Mostly, phenotypically streptomycin resistant isolates were detected with Lys43Arg (53.57%, 30/56) mutation in rpsL gene (Table S5). In addition, 54.8% (17/31) phenotypically fluoroquinolone-resistant isolates were identified with mutations in the gyrA gene like Asp94Gly (29.0%, 9/31) and Ala90Val (22.6%, 7/31), but no isolates have mutation in gyrB gene (Table S6). Among kanamycin-resistant isolates, 66.7% of isolates (6/9) have 1401a>g mutation in rrs gene (Table S7).

Discussion

Tuberculosis remains a major public health problem in China, and the population characteristics and drug resistance profiles in different regions showed a great difference.^3,21,22 In this study, we provide the population characteristics, transmission dynamics and drug resistance characteristics of Mtb in Guangxi, China. Our results show an initial picture of the Mtb population structure in Guangxi, which has high genetic diversity including lineage 1 to lineage 4. Our data show that 67.5% and 30.1% of Mtb isolates belonged to lineage 2 and Lineage 4, respectively, which is in line with the previous report that lineage 2 and lineage 4 are most prevalent and distributed throughout in China.^3,21,23,24 Although most isolates belonged to Beijing lineage (L2), 9.6% of them are proto-Beijing genotype (L2.1), which is not common in other regions of China and is a significant character of the Mtb isolates in Guangxi.³ Previous studies have shown that L2.1 isolates are also present in parts of Vietnam and Thailand,^19,25 especially in Thailand where L2.1 has a great propensity to transmit. Overall, the sub-lineage L2.1 is geographically restricted and need more widespread surveillance in other areas.

Drug resistance was formerly assumed to be related to the Beijing genotype,²⁶ but a recent study from China showed that the Beijing genotype is less associated with drug resistance.²⁷ Conversely, the non-Beijing genotype was more likely to be resistant to isoniazid, rifampicin and ofloxacin in this study, which showed the findings have been inconsistent because of different geographic settings and isolates populations. This question requires more research, particularly from the standpoint of molecular evolution. Our data showed that the isolates belonging to L2.1 exhibit low drug resistance, which is inconsistent with the recent report that proto-Beijing genotype (L2.1) isolates have a great propensity to be XDR-TB in Thailand.¹⁹ Here, we also found that the drug-resistance in border regions was higher than in non-border regions, which suggests that we should pay more attention to generating efficient tuberculosis prevention and control strategies in border regions.

Our data show a high rate of isoniazid resistance among the Mtb isolates in Guangxi, and the prevalence of Hr-TB is significantly higher than the prevalence of rifampicin resistance, which includes MDR-TB. Importantly, the rate of Hr-TB is higher than the previous report in China (5.7%)²⁸ and the global prevalence of Hr-TB (7.4% in new TB patients and 11.4% in retreated TB patients).²⁹ According to a meta-analysis, patients with INH-resistant TB who received first-line HRZE therapy had high rates of unfavorable outcomes, such as failure, relapse, and acquired multidrug resistance (11%, 10%, and 8%, respectively).²⁹ To date, isoniazid is one of the most powerful drugs for the treatment of TB, but the procedures are not standardized and effective enough for the diagnosis and promptly modified treatment for Hr-TB in some areas, which emphasizes the urgent necessity to prioritize its rapid detection and effective therapy to curb the future epidemics of MDR-TB.

Concerning the molecular mechanism of resistance to anti-TB drugs, the associated gene mutations were screened by the whole genome sequence. In concordance with previous studies, the rifampicin resistance-determining region (RRDR) of rpoB and katG 315 codon mutations clearly predominated in the assessment of rifampicin- and isoniazid-resistant mutations.^30–32 And we also found that there are some combinations of mutations like secondary mutations in rpoC and non-RRDR of rpoB, which may have a compensatory function of the isolates.³³ The study showed that EMB resistance was associated with embCAB mutations and embB was the most frequently mutated gene.³⁴ Previous studies show that SM resistance has been shown to be associated with rpsL, rrs, and gid genes, of which rpsL mutations were predominant,³⁵ and this study showed 64.3% of isolates have mutation in rpsL. Mutation in gyrA and gyrB have been known to be related to ofloxacin, the current study has found 77.4% of the isolates have mutation in gyrA, but no isolates in gyrB. This is consistent with the report that gyrA mutations are more common than gyrB mutations.³² More importantly, some Mtb strains with phenotypical resistance against anti-TB drugs (eg, rifampicin, isoniazid, ofloxacin, etc.) had no detectable gene mutations in our study, suggesting that there may be other genomic variants, such as insertions and deletions,³⁶ which need further clarification. In addition, alternative mechanisms, such as drug efflux pump and decreased cell wall permeability to drugs, may also be related to drug resistance in MTBC.^37,38

Our study also found that 12.6% (58/462) of the pulmonary TB cases had genomic-clustered isolates in Guangxi during a period of fewer than two years. The clustering rate was similar to Luodian (13.1%) in a one-year study,²¹ lower than Shenzhen (15.2%) in five years study.²⁸ The clustering rate in the current study may have been underestimated because molecular epidemiological studies of tuberculosis typically utilize a sampling period of two years or longer. Thus, further long-term investigation is required to comprehensively characterize TB transmission. Of note, most isolates within a cluster were collected from the same county in our study, indicating that the recent transmission predominantly occurred in the local area. Interestingly, some isolates of clusters showed different drug-resistant types, suggesting the progression and accumulation of mutations linked to drug resistance.³⁹ But some of them were isolated in the same year, this result could be partially explained by the fact that not all cases from the possible cluster were enrolled in this study, and some index cases with primary resistance may have been overlooked in the chosen population.³⁹

There are several limitations in the current study. First, the proportion of recent transmission events might be underestimated since our collection of isolates is a period of less than two years. Second, the isolates were only selected from five sites, which may not be representative of the overall situation of the Guangxi area. Thus, the generalizability of the results obtained in this study may be limited. Third, the clinical information collected in this study was relatively limited, with only basic clinical symptoms, thus the correlation between clinical symptoms and genomic variants was not analyzed. In addition, due to the retrospective nature of this study, treatment outcome was not collected, which prevented us from exploring the correlation between drug-resistance profile and clinical outcomes.

Conclusion

In summary, we obtained the diversity of Mtb isolates and the trend of drug resistance, and get a better understanding of the epidemic characteristic of TB in Guangxi, China. We found that the genomic of Mtb isolates is high diversity including four major lineages (L1-L4), and it has rare proto-Beijing isolates (L2.1). The genomic clusters contained 58 isolates, with the clustering rate of 12.6%. In addition, our data show the high rate of isoniazid resistance among the Mtb isolates in Guangxi, and the prevalence of Hr-TB is significantly higher than the prevalence of rifampicin resistance, including MDR-TB, which emphasizes the urgent necessity to prioritize its rapid detection and effective therapy for Hr-TB and to curb the future epidemics of MDR-TB. The isolates in border monitoring points are more likely to be drug resistant than no-border monitoring points, which needs to be noted. Overall, this work provided a deep analysis of the genomic population characteristics and drug resistance of Mtb in Guangxi, which could contribute to developing effective TB prevention and control strategies.