Drug Resistance of Mycobacterium tuberculosis Based on Whole-Genome Sequencing in the Yi Ethnic Group, Sichuan Province, China

This study investigated drug-resistant tuberculosis (DR-TB) in the Yi ethnic group. The study was designed to identify risk factors for DR-TB and its relationship with HIV/AIDS. To establish the resistance to antituberculosis drugs, whole-genome sequencing (WGS) was performed using culture-positive Mycobacterium tuberculosis samples collected from people of the Yi ethnic group from March 2019 to March 2021. Baseline characteristics were obtained from China's tuberculosis surveillance system. A total of 116 M. tuberculosis strains were included in the final analysis. Lineage 2.2 (75.86%) was the dominant sublineage, followed by lineage 4.5 (18.97%) and lineage 4.4 (5.17%). The rates of rifampicin-resistant (RR-TB), multidrug-resistant (MDR-TB), and preextensively drug-resistant TB (pre-XDR-TB) were 18.97%, 10.34%, and 6.03%, respectively. Drug-resistant strains were not found in the elderly (age ≥ 65 years). The proportions of RR/MDR-TB and pre-XDR-TB cases among re-treatment patients were higher than those among new patients (χ2 = 12.155, P = 0.003; χ2 = 22.495, P = 0.001, respectively). The pre-XDR-TB case proportions were higher among female patients than among males and higher among referred patients (χ2 = 5.456, P = 0.032; χ2 = 15.134, P = 0.002, respectively). The rates of RR/MDR-TB and pre-XDR-TB did not differ appreciably among groups with different HIV infection statuses nor lineage populations. DR-TB poses a serious challenge to the Yi ethnic group. Re-treatment patients, women, and referred patients were at high risk of MDR/RR-TB or pre-XDR-TB while HIV and lineage 2 had negligible association with drug resistance. Whole-genome sequencing should be used to guide the design of treatment regimens and to tailor public interventions.


Introduction
Tuberculosis (TB) is a chronic communicable disease and remains a serious public threat worldwide. Although the global burden is stable, drug-resistant TB (DR-TB) is a major threat to the WHO End TB Strategy due to the limited treatment options, longer treatment duration, lower cure rate, and the need for more second-line anti-TB drugs with greater side effects and higher prices [1]. According to the WHO, China is estimated to have the second-highest inci-dence of TB; multidrug-resistant TB or rifampicin-resistant TB (RR/MDR-TB) is only lower than India [2].
Liangshan Yi Autonomous Prefecture (Liangshan) of Sichuan province, located in southwest China, is a mountainous area populated by the Yi ethnic group and has an underdeveloped economy. Our previous study found that the TB epidemic is particularly serious in the Yi ethnic group, especially in regions with high HIV/AIDS prevalence [3,4]. Butuo County (Butuo) is one of the counties in Liangshan most plagued by HIV/AIDS; it reportedly has the highest incidence of the disease in China, according to a county-level based survey. However, no study has investigated the situation of DR-TB and its relationship with HIV/AIDS in the Yi ethnic group.
The culture-based phenotypic drug sensitivity test, which detects the growth of Mycobacterium tuberculosis in the presence of anti-TB drugs, is limited by the slow mycobacterial duplication time and lack of a uniform minimum inhibitory concentration analysis or evaluation. Genotypic laboratory tests can quickly detect the presence of DNA mutations conferring resistance but can only identify a limited number of resistance genes. Recently, whole-genome sequencing (WGS) has been widely used in testing for drug resistance because it offers more accurate prediction and better determination of resistance to most drugs. WGS also has high predictive power to infer resistance profiles in areas with high TB burdens [5]. Furthermore, WGS can be used to understand evolution, lineages, and genomic variations [6]. Compared with analysis of the restriction fragment length polymorphisms of the IS6110 gene (IS6110-RFLP), spacer oligonucleotide typing (spoligotyping), and mycobacterial interspersed repetitive unit-variable number of DNA tandem repeats, WGS provides high-resolution data that can be used to determine genetic relationships, especially between strains that are very close at the genetic level [7,8].
In this study, we chose Butuo County as the research site to investigate the genetic diversity and drug-resistant profile of local circulating strains in the local Yi ethnic group, using WGS for the first time, to provide a scientific basis for TB prevention and control. Patients with pulmonary tuberculosis (PTB) were diagnosed following the "Criteria for PTB Diagnosis" of China (WS288-2017). All cases of PTB among the Yi ethnic group were included in this research. Patients from other ethnic groups were excluded. Sputum samples were collected from all study participants. Each sample was liquefied using 4% NaOH and inoculated into tubes with acidified Löwenstein-Jensen medium for further culture. The identification of M. tuberculosis included p-nitrobenzoic acid testing. Subsequently, M. tuberculosis strains were tested for WGS.

Patient
Information. Demographic information (gender, age, occupation, and census register) and clinical characteristics of patients (HIV infection status, history of previous TB treatment, visiting hospital delay, health system delay, and patient source) were obtained from the Tuberculosis Information Management System of China. Visiting hospital delay was defined as the time interval between the date of symptom onset and the patient's first consultation with a TB-designated hospital, and health system delay was defined as the time interval between the patient's first consultation with TB-designated hospitals and the date of diagnosis. Pulmonary TB patients were detected via clinic visits due to symptoms, referral, tracing, or health examination.

Whole-Genome
Sequencing. Extraction and purification of genomic DNA were carried out following Bacterial DNA Extraction Kit (Gene-Optimal, 60300 K-50 T) protocols. Libraries were constructed on the Illumina platform using an FS DNA Lib Prep Kit V6 (RK20259). The samples were sequenced using an Illumina NovaSeq 6000 sequencer. Base-calling was performed using PE 150 software. All whole-genome sequencing procedures were performed by Shanghai Gene-Optimal Science & Technology Co. Ltd. (Shanghai, China). Cutadapt (v1.15) was used to trim adapter sequences at the tail of sequencing reads. Then, sequencing reads were aligned to reference genome (H37Rv, NC000962.3) using BWA (v0.7.15). Duplicated reads were marked by Picard (v2.4.1). Qualimap (v2.2.1) was used to calculate base quality metrics, genome mapping rate, and the coverage of targeted regions. Mycobacterium tuberculosis samples were selected for variant detection and annotation. The selection criteria were as follows: detected by 100% coverage of reads mapping to M. tuberculosis specific sequences (H37Rv region 315947-316534), with an average depth > 5. SNP and InDel calling was performed following the SAMtools (v1.6) and VarScan (v2.3.9). Large deletions were detected using Delly (v0.8.7). Variant filtering and annotating were done using a finely tuned inhouse script. The WGS data obtained in our study were compared with available databases of M. tuberculosis drug resistance-associated gene mutations [9,10]. Strains with a sequence depth less than 20× or a genome coverage less than 95% were excluded from the analysis. The phylogenetic tree was constructed using the maximum likelihood method and reconstructed using RAxML-NG12 (v.1.0.2) with "-model GTR + G + ASC_LEWIS". All strains were phylogenetically classified according to the SNPs barcode nomenclature proposed by Napier G [11]. The phylogenetic tree was visualized and modified using iTOL (https://itol.embl.de/).

Statistical Analysis.
Chi-squared or Fisher's exact tests were used for categorical data. Correlation analysis was performed to explore the risk factors for high-drug resistance. All statistical analyses were performed in SPSS version 20 software (SPSS Inc., Chicago, Illinois.). P values less than 0.05 were considered statistically significant.

Definitions.
Rifampicin-resistant (RR-TB) or RR/MDR-TB was classified as M. tuberculosis resistance to at least rifampicin. MDR-TB was classified as M. tuberculosis resistance to at least isoniazid and rifampicin. Pre-XDR-TB was defined as MDR-TB with additional resistance to any fluoroquinolones (moxifloxacin or ofloxacin), based on a new definition released by WHO in January 2021 [12].
2.6. Ethics Approval. The study protocol was approved by the Ethical Review Committee at the Biomedical Ethics Committee, West China Hospital, Sichuan University (No. 2019-151). All participants provided written informed consent after reviewing the description of the study.

Prediction of Drug Resistance and Mutation
Characteristics. As predicted, one or more drug resistance mutations were detected in 33 strains, which means 28.45% of the strains were drug-resistant; the rest were considered as susceptible strains. In 116 strains, rifampicinresistant (RR-TB), multidrug-resistant (MDR-TB), and preextensively drug-resistant (pre-XDR-TB) strains were found in 22, 12, and 7 patients, accounting for 18.97%, 10.34%, and 6.03%, respectively.
The resistances of 116 strains to different drugs are shown in Table 2. In this study, 12 strains resistant to antituberculosis drugs were detected, and the drug resistance rates for rifampicin (RIF: 22 The Sm resistance gene mutations occurred in rpsL, rrs, and gid, and a total of five resistance mutation types were generated, mainly rpsL_K88R (41.67%, 5/12). A total of ten mutation types were generated in the resistance genes of fluoroquinolones, mainly gyr_A90V (40.00, 4/10). Ethambutol resistance was found in eight strains and a total of five mutation types were produced; the main mutation types Note: nine cases were missing information for census register, health system delay, and visiting hospital delay. Seven cases were missing information for patient source.
3 Journal of Immunology Research were embB_M306V (37.50%, 3/8). One (12.50%, 1/8) of the EMB-resistant strains has combined mutations at two loci. The resistance gene mutations to other drugs are shown in Table 2.

Discussion
To the best of our knowledge, this is the first study to use whole-genome sequencing technology to predict drug resistance and analyse genetic diversity of TB among members of the Yi ethnic group with high HIV/AIDS prevalence. Our findings help understand the drug resistance spectrum and toxicity of local epidemic strains and provide a scientific basis for TB prevention and control.
The epidemic situation of DR-TB in the Yi ethnic group is serious, with a 10.34% incidence rate of MDR-TB, higher than the average level in Sichuan province (8.14%) [13] and higher than that in other regions of China [14][15][16][17]. Furthermore, the MDR-TB rates among new and re-treated patients of the Yi ethnic group (7.55% and 40.00%, respectively) are  Journal of Immunology Research much higher than the baselines reported in the national survey in 2007 (5.7% and 25.6%, respectively) [18] and higher than the global average in the last decade (3-4% and 18-21%, respectively) [2]. Although the Yi ethnic group already has a serious epidemic of drug-resistant TB, the situation has the potential to deteriorate further. The overall RR/MDR-TB rate was 18.97%, while the rates in new and re-treatment patients were 15.09% and 60.00%, respectively, which are higher than those reported in a previous Indian study [19]. In addition to MDR, 7.55% (8/106) of new patients and 20.00% (2/10) of re-treated patients were resistant to rifampicin. If these patients are not effectively treated, they can easily develop MDR-TB. In addition, 31.82% (7/22) of all RR-TB patients suffered from pre-XDR-TB, a proportion higher than the global rate (20%) [20] and higher than the Chinese rate (25%) [18]. These patients are easily susceptible to developing XDR-TB, which would further worsen the local TB control.
The most common drug-resistant mutations (RIF-, fluoroquinolone-, and EMB-resistant strains) in the Yi ethnic group were consistent with previous reports, but differed in the frequency of resistance. In terms of the MDR-TB strain in all China, the majority mutation types of INH and Sm resistance-associated genes were katG315T and rpsLK43R, respectively, which differed from the results of our study [21]. These differences may be related to the study sample size, strain genotype, medication habits, and patient compliance.
The results showed strong associations between pre-XDR-TB and treatment history, gender, and source of patients; meanwhile, MDR/RR-TB was associated with treatment history. In previous studies, treatment history has been reported to be the most influential risk factor for developing drug tolerance [22][23][24]. Patients with a history of antituberculosis treatments might have problems with irregular treatment or interruption of treatment, leading to increased selection pressure for drug resistance mutations and ultimately drug tolerance [24].

Journal of Immunology Research
In addition, the proportion of pre-XDR-TB among female patients was higher than that in males, which was possibly related to lower incomes among females, resulting in poor treatment conditions, and the fact that females were more prone to irregular medication [25]. The pre-XDR-TB rate among referred patients was higher than that among patients from other sources. Patients referred for treatment in nondesignated medical institutions or outpatient clinics could not get the correct diagnosis nor standardized medication in a timely manner; in particular, abuse of fluoroquinolones may have occurred, resulting in the generation of pre-XDR-TB.
The global epidemiology of drug-resistant TB in HIVinfected persons is not known owing to a lack of information. The global project of antituberculosis drug resistance explored the interaction between HIV infection and drugresistant TB in seven settings, none of which had a high prevalence of HIV infection. The results showed significant associations only in two settings, in Latvia and Ukraine [26]. Eldholm et al. and Masenga et al. reported that HIV coinfection was not significantly associated with drug-resistant TB [27,28]. Our result for this study area with high prevalence of HIV infection echoed their finding. In our study, the drug-resistant strain was not found in the elderly (aged at least 65 years), which was a similar result to that of the Indian study and may be explained by reduced survival among those with TB infections [19].
Twenty-three percent of the global TB burden is ascribed to lineage 2 strains, which are mainly distributed across East Asia and have been associated with drug resistance, hypervirulence, and high transmissibility [29]. The prevalent TB bacteria in the Yi ethnic group were lineage 2 (all Beijing genotype (lineage 2.2)), accounting for 75.21% of cases in this study, which is close to the average level of China (70%) [30]. However, lineage 2 was not the major contributor for drug resistance in our study, which was consistent with the conclusion of another study conducted in southern China [31]. More research is required to clarify the underlying reasons for the lack of association of lineage 2 with the prevalence of drug resistance TB in the Yi ethnic group.

Conclusions
In summary, the prevalent strains in the Yi ethnic group belonged to lineage 2, and our efforts to mitigate the impacts of TB should be focused on this genotype. The Yi ethnic group are facing a serious situation of DR-TB. Retreatment patients were at high risk of MDR/RR-TB and pre-XDR-TB. Measures should be taken to reduce acquired drug resistance during treatment through standardizing therapy and management of patients. Women and referred patients were more susceptible to pre-XDR-TB. Those key populations should receive more attention in the management of TB. The drug-resistant strain was not found in the elderly (≥65 years of age). HIV coinfection and lineage 2 were not found to be significantly associated with drugresistant TB. Whole-genome sequencing, as the most promising tool for predicting drug resistance and analysing genetic diversity, should be used to guide the design of treatment regimens and tailor public interventions.

Data Availability
The raw sequence data reported in this paper have been deposited in the Genome Sequence Archive (Genomics, Proteomics & Bioinformatics 2021) in the National Genomics Data Center (Nucleic Acids Res 2022), China National Center for Bioinformation/Beijing Institute of Genomics, Chinese Academy of Sciences (GSA: CRA008002), and are publicly accessible at https://ngdc.cncb.ac.cn/gsa.

Conflicts of Interest
The authors declare no competing interests.